Python Walkthrough
This section walks through using the MSG Python interface
to evaluate spectra and photometric colors for a model of Sirius. The code fragments are presented as a
sequence of Jupyter notebook cells, but pasting all of them
into a Python script should work also. Alternatively, the entire notebook (including this text) is available in the source distribution at $MSG_DIR/docs/source/user-guide/walkthroughs/python-walkthrough.ipynb
.
Preparation
Before starting Jupyter, make sure that you have installed the pymsg
package as described in the Quick Start and/or Installation chapters. Also, ensure that the NumPy, Matplotlib and
Astropy packages are installed on your
system. Finally, confirm that the MSG_DIR
environment
variable is set.
Initialization
Launch Jupyter, create a new Python3 notebook, and initialize the notebook with the following statements:
[1]:
# Import standard modules
import os
import numpy as np
import astropy.constants as con
import astropy.units as unt
import matplotlib.pyplot as plt
# Import pymsg
import pymsg
# Set plot parameters
%matplotlib inline
plt.rcParams.update({'font.size': 16})
Loading & Inspecting the Demo Grid
For demonstration purposes MSG includes a low-resolution specgrid (spectroscopic grid) file, $MSG_DIR/data/grids/sg-demo.h5
[1]. Load this demo grid into memory[2] by creating a
new pymsg.SpecGrid
object:
[2]:
# Load the SpecGrid
MSG_DIR = os.environ['MSG_DIR']
GRID_DIR = os.path.join(MSG_DIR, 'data', 'grids')
specgrid_file_name = os.path.join(GRID_DIR, 'sg-demo.h5')
specgrid = pymsg.SpecGrid(specgrid_file_name)
This object has a number of (read-only) properties that inform us about the parameter space spanned by the grid:
[3]:
# Inspect grid parameters
print('Grid parameters:')
for label in specgrid.axis_labels:
print(f' {label} [{specgrid.axis_x_min[label]} -> {specgrid.axis_x_max[label]}]')
print(f' lam [{specgrid.lam_min} -> {specgrid.lam_max}]')
Grid parameters:
Teff [3500.0 -> 50000.0]
log(g) [0.0 -> 5.0]
lam [2999.9999999999977 -> 9003.490078514782]
Here, Teff
and log(g)
correspond (respectively) to the
\(\Teff/\kelvin\) and \(\log_{10}(g/\cm\,\second^{-2})\) photospheric parameters, while
lam
is wavelength \(\lambda/\angstrom\).
Plotting the Flux
With the grid loaded, let’s evaluate and plot a flux spectrum for Sirius A. First, store the photospheric parameters for the star in a dict:
[4]:
# Set Sirius A's fundamental parameters
M = 2.02 *con.M_sun
R = 1.71 * con.R_sun
L = 25.4 * con.L_sun
p = 0.379
d = 1/p * con.pc
# Set photospheric parameters dict
Teff = (L/(4*np.pi*R**2*con.sigma_sb))**0.25/unt.K
logg = np.log10(con.G*M/R**2/(unt.cm/unt.s**2))
x = {'Teff': Teff, 'log(g)': logg}
print(Teff, logg)
9909.16258994498 4.277426775965449
(the mass, radius and luminosity data are taken from Liebert et al., 2005, while the parallax comes from Simbad). Then set up a uniformly-spaced wavelength abcissa for a spectrum spanning the visible range, \(3\,000\,\angstrom\) to \(7\,000\,\angstrom\).
[5]:
# Set up the wavelength abscissa
lam_min = 3000.
lam_max = 7000.
lam = np.linspace(lam_min, lam_max, 501)
lam_c = 0.5*(lam[1:] + lam[:-1])
Here, the array lam
defines the boundaries of 500 wavelength bins
\(\{\lambda_{i},\lambda_{i+1}\}\) (\(i=1,\ldots,500\)) and the
array lam_c
stores the central wavelength of each bin.
With all our parameters defined, evaluate the observed flux using a
call to the pymsg.SpecGrid.flux()
method, and then plot it:
[6]:
# Evaluate the observed flux
F_obs = (R/d)**2*specgrid.flux(x, lam)
# Plot
plt.figure(figsize=[8,8])
plt.plot(lam_c, F_obs)
plt.xlabel(r'$\lambda ({\AA})$')
plt.ylabel(r'$F^{\mathrm{o}}_{\lambda}\ ({\rm erg\,cm^{-2}\,s^{-1}}\,\AA^{-1})$')
[6]:
Text(0, 0.5, '$F^{\\mathrm{o}}_{\\lambda}\\ ({\\rm erg\\,cm^{-2}\\,s^{-1}}\\,\\AA^{-1})$')

(See equation (7) of the MSG Fundamentals chapter to understand the origin of the expression for F_obs
). The plot looks about right — we can clearly see the Balmer series, starting with H\(\alpha\) at \(6563\,\angstrom\).
Plotting the Intensity
Sometimes we want to know the specific intensity of the radiation
emerging from each element of a star’s surface; an example might be when we’re
modeling a rotationally distorted star, whose surface geometry and properties vary with latitude. For this, we can use the pymsg.SpecGrid.intensity()
method.
Here’s a demonstration of this method in action, plotting the specific intensity across the H\(\alpha\) line profile for ten different values of the emergence direction cosine \(\mu=0.1,0.2,\ldots,1.0\).
[7]:
# Set up the wavelength abscissa
lam_min = 6300.
lam_max = 6800.
lam = np.linspace(lam_min, lam_max, 100)
lam_c = 0.5*(lam[1:] + lam[:-1])
# Loop over mu
plt.figure(figsize=[8,8])
for mu in np.linspace(1.0, 0.1, 10):
# Evaluate the intensity
I = specgrid.intensity(x, mu, lam)
# Plot
if mu==0.1 or mu==1.0:
label=r'$\mu={:3.1f}$'.format(mu)
else:
label=None
plt.plot(lam_c, I, label=label)
plt.xlabel(r'$\lambda ({\AA})$')
plt.ylabel(r'$I_{\lambda}\ ({\rm erg\,cm^{-2}\,s^{-1}}\,\AA^{-1}\,sr^{-1})$')
plt.legend()
[7]:
<matplotlib.legend.Legend at 0x7ff643465670>

We can clearly see that limb-darkening in the line core is much weaker than in the continuum — exactly what we expect from such a strong line.
Evaluating Magnitudes & Colors
As a final step in this walkthrough, let’s evaluate the magnitude and
colors of Sirius A in the Johnson photometric system. We can do this by
creating a new pymsg.PhotGrid
object for each passband:
[8]:
# Load the PhotGrids
PASS_DIR = os.path.join(MSG_DIR, 'data', 'passbands')
filters = ['U', 'B', 'V']
photgrids = {}
for filter in filters:
passband_file_name = os.path.join(PASS_DIR, f'pb-Generic-Johnson.{filter}-Vega.h5')
photgrids[filter] = pymsg.PhotGrid(specgrid_file_name, passband_file_name)
(for convenience, we store the pymsg.PhotGrid
objects in a
dict, keyed by filter name). In the calls to the object constructor
pymsg.PhotGrid()
, the first argument is the name of a
specgrid file (here, the demo grid), and the second argument is the
name of a passband file; a limited set of these files is
provided in the $MSG_DIR/data/passbands
subdirectory [3].
The observed fluxes of Sirius A are then found using the pymsg.PhotGrid.flux()
method:
[9]:
# Evaluate the observed fluxes (each normalized to the passband
# zero-point flux)
F_obs = {}
for filter in filters:
F_obs[filter] = (R/d)**2*photgrids[filter].flux(x)
To convert these into apparent magnitudes, we apply Pogson’s logarithmic formula:
[10]:
# Evaluate apparent magnitudes and print out magnitude & colors
mags = {}
for filter in filters:
mags[filter] = -2.5*np.log10(F_obs[filter])
print(f"V={mags['V']}, U-B={mags['U']-mags['B']}, B-V={mags['B']-mags['V']}")
V=-1.4466779855804992, U-B=-0.05964681835037444, B-V=0.001121934559611848
Reassuringly, the resulting values are within 10 millimag of Sirius A’s apparent V-band magnitude and colors, as given by Simbad.
Footnotes