In [1]:

    

# Import some numerical math and plotting packages
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

# Import the MesaProfile package
from MesaProfile import MesaProfile

# Import the OrderedDict package for working with MesaProfile outputs
from collections import OrderedDict


    

In [2]:

    

# Read the input MESA profile                                                                                                                                                                                      
mesa = MesaProfile('profile64.data')
mstar = mesa.getStar()


    

In [3]:

    

# List the data fields available in the profile
for k in mstar.keys():
    print(k)


    

    




zone
logT
logRho
logP
logR
luminosity
logL
velocity
entropy
conv_mixing_type
csound
v_div_csound
v_div_r
eta
mu
logdq
dq_ratio
q
radius
temperature
tau
logtau
pressure
pgas_div_ptotal
logPgas
grada
cp
logS
gam
free_e
chiRho
chiT
abar
zbar
z2bar
ye
log_opacity
eps_nuc
non_nuc_neu
eps_grav
mlt_mixing_length
log_D_mix
log_conv_vel
conv_vel_div_csound
log_mlt_D_mix
pressure_scale_height
gradT
gradr
dlnd_dt
dlnT_dt
mass
mmid
logxq
h1
he3
he4
c12
n14
o16
ne20
mg24
si28
pp
cno
tri_alfa
burn_c
burn_n
burn_o
burn_ne
burn_na
burn_mg
burn_si
burn_s
burn_ar
burn_ca
burn_ti
burn_cr
burn_fe
c12_c12
c12_o16
o16_o16
pnhe4
photo
other
extra_heat
gradr_sub_grada
brunt_N2
brunt_B
brunt_nonB
sign_brunt_N2
lamb_S2
lamb_S
log_brunt_nu
log_lamb_Sl2
log_lamb_Sl3
brunt_N_div_r_integral
k_r_integral
brunt_N2_sub_omega2
sl2_sub_omega2
logQ
radiuscm
model_number
num_zones
initial_mass
initial_z
star_age
time_step
Teff
photosphere_L
photosphere_r
center_eta
center_h1
center_he3
center_he4
center_c12
center_n14
center_o16
center_ne20
star_mass
star_mdot
star_mass_h1
star_mass_he3
star_mass_he4
star_mass_c12
star_mass_n14
star_mass_o16
star_mass_ne20
he_core_mass
c_core_mass
o_core_mass
si_core_mass
fe_core_mass
neutron_rich_core_mass
tau10_mass
tau10_radius
tau100_mass
tau100_radius
dynamic_time
kh_timescale
nuc_timescale
power_nuc_burn
power_h_burn
power_he_burn
power_neu
burn_min1
burn_min2

In [4]:

    

# Find out what isotopes the profile contains
isotope_list = mesa.getIsotopes()  
nisotopes = len(isotope_list)
print('The MESA profile contains {} isotopes.'.format(nisotopes))


    

    




The MESA profile contains 9 isotopes.

In [5]:

    

# Plot the abundances in the profile    
# Yellow is high mass number isotopes, purple is low mass number isotopes
fig, ax = plt.subplots()
colormap = plt.get_cmap('viridis')
isocolors = [colormap(i) for i in np.linspace(0, 1, nisotopes)]

legend_handles = []
for i, isotope in enumerate(isotope_list):
    handle = ax.plot(mstar['mass'], mstar[isotope], color=isocolors[i], label=isotope)
    legend_handles.append(handle)

plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
ax.set_xlabel('Mass ($M_\\odot$)')
ax.set_ylabel('Mass Fractions')


    

    
Out[5]:





<matplotlib.text.Text at 0x1119f2198>






    

In [6]:

    

# Plot the electron fraction (Ye) profile
#Number of electrons that you have for every baryon in your plasma. If you have extra neutrons, 
# the fraction will go below .5
fig, ax = plt.subplots()
ax.plot(mstar['mass'],mstar['ye'],'b')
ax.set_xlabel('Mass ($M_\\odot$)')
ax.set_ylabel('Electron Fraction ($\\mathrm{Y_e}$)')


    

    
Out[6]:





<matplotlib.text.Text at 0x114332a58>






    

In [40]:

    

# Plot the density profile
fig, ax = plt.subplots()
x = mstar['mass']
y = 10.0**mstar['logRho']
ax.plot(x,y,'b')
ax.set_xlabel('Mass ($M_\\odot$)')
ax.set_ylabel('Density ($g/cm^3$)')
print(x[0], y[0])


    

    




8.91139889744e-09 4625064900.09






    

In [38]:

    

# Plot the temperature profile
fig, ax = plt.subplots()
x = mstar['mass']
y = 10.0**mstar['logT']
ax.plot(x,y,'b')
ax.set_xlabel('Mass ($M_\\odot$)')
ax.set_ylabel('Temperature ($K$)')
print(x[0], y[0])


    

    




8.91139889744e-09 806887508.625






    

In [39]:

    

# Plot the mass-radius profile
fig, ax = plt.subplots()
x = mstar['radiuscm']
y = mstar['mass']
ax.plot(x,y,'b')
ax.set_xlabel('Radius (cm)')
ax.set_ylabel('Mass ($M_\\odot$)')
print(x[0], y[0])


    

    




97015.1752895 8.91139889744e-09






    

In [11]:

    

# Plot the luminosity profile
fig, ax = plt.subplots()
ax.plot(mstar['mass'],mstar['luminosity'],'b')
ax.set_xlabel('Mass ($M_\\odot$)')
ax.set_ylabel('Luminosity ($J/s$)')


    

    
Out[11]:





<matplotlib.text.Text at 0x1192acdd8>






    

In [12]:

    

# Plot the entropy profile
fig, ax = plt.subplots()
ax.plot(mstar['mass'],mstar['entropy'],'b')
ax.set_xlabel('Mass ($M_\\odot$)')
ax.set_ylabel('Entropy ($J/K$)')


    

    
Out[12]:





<matplotlib.text.Text at 0x1193aa400>






    

In [16]:

    

# Plot the dr profile
fig, ax = plt.subplots()
mstar["dr"] = np.diff(mstar["radiuscm"])
ax.plot(mstar['mass'][:-1],mstar['dr'],'b')
ax.set_xlabel('Mass ($M_\\odot$)')
ax.set_ylabel('dr ($cm$)')


    

    
Out[16]:





<matplotlib.text.Text at 0x1195d2f98>






    

In [ ]: