Source code for as4012_sstr.reader
# Python module to handle MESA-Web output files
# Original Code by
# (c) 2021-2023 Rich Townsend (rhtownsend@me.com)
#
# This module is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, version 3.
#
# This module is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# 2024
# Updated for St Andrews AS4012 by Ian Czekala (ic95@st-andrews.ac.uk).
# Removed plain Python dict option for return type
# Reformatted docstrings to Google/NumPy style.,
import os.path as op
import re
import numpy as np
from astropy.table import Table
[docs]
def read_history(filename):
"""Read data from a MESA-Web history file
Parameters
----------
filename
string giving name of history file
Returns
-------
:class:`astropy.Table` object containing header and history data.
Notes
-----
See `Astropy.table <https://docs.astropy.org/en/stable/table/>`_ for more
information on working with the returned object.
>>> table = read_history("trimmed_history.data")
>>> log_L = table["log_L"]
**Header Data**: The following keys/value pairs in the returned Table.meta dict
contain header data -- i.e., scalars describing time-independent properties
of the star. Where applicable, units are given in square brackets [].
* version_number -- version number of MESA
* compiler -- name of compiler used to build MESA
* build -- version of compiler used to build MESA
* MESA_SDK_version -- version of Software Development Kit used to build MESA
* date -- date on which MESA-Web calculation began
* burn_min1 -- 1st limit for reported burning [erg/g/s]
* burn_min2 -- 2nd limit for reported burning [erg/g/s]
**History Data**: The following keys/value pairs in the returned
Table contain history data -- i.e., arrays describing global properties of the
star over a sequence of time-steps. Where applicable, units are given in
square brackets [].
* model_number -- model number
* star_age -- stellar age [years]
* star_mass -- stellar mass [Msun]
* log_L -- log10(stellar luminosity [Lsun])
* log_R -- log10(stellar radius [Rsun])
* log_Teff -- log10(effective temperature [K])
* log_center_T -- log10(center temperature [K])
* log_center_Rho -- log10(center density [g/cm^3])
* log_center_P -- log10(center pressure [dyn/cm^2])
* center_h1 -- center 1H mass fraction
* center_he3 -- center 3He mass fraction
* center_he4 -- center 4He mass fraction
* center_c12 -- center 12C mass fraction
* center_n14 -- center 14N mass fraction
* center_o16 -- center 16O mass fraction
* center_ne20 -- center 20Ne mass fraction
* center_mg24 -- center 24Mg mass fraction
* center_si28 -- center 28Si mass fraction
* center_s32 -- center 32S mass fraction
* center_ar36 -- center 36Ar mass fraction
* center_ca40 -- center 40Ca mass fraction
* center_ti44 -- center 44Ti mass fraction
* center_cr48 -- center 48Cr mass fraction
* center_fe52 -- center 52Fe mass fraction
* center_fe54 -- center 54Fe mass fraction
* center_fe56 -- center 56Fe mass fraction
* center_ni56 -- center 56Ni mass fraction
* center_degeneracy -- center electron degeneracy parameter [kB*T]
* center_ye -- center average charge per baryon [e]
* center_entropy -- center entropy [kB]
* compactness_parameter -- (m/Msun)/(R(m)/1000km) for m = 2.5 Msun
* dynamic_timescale -- dynamical timescale [s]
* kh_timescale -- Kelvin-Helmholtz timescale [s]
* nuc_timescale -- nuclear timescale [s]
* pp -- log10(total pp luminosity [Lsun])
* cno -- log10(total CNO luminosity [Lsun])
* tri_alfa -- log10(total triple-alpha luminosity [Lsun])
* log_LH -- log10(total H-burning luminosity, excluding
neutrinos [Lsun])
* log_LHe -- log10(total He-burning luminosity, excluding
neutrinos [Lsun])
* log_LZ -- log10(total metal-burning luminosity, excluding
neutrinos [Lsun])
* log_Lneu -- log10(total neutrino luminosity [Lsun])
* he_core_mass -- mass of helium core [Msun]
* c_core_mass -- mass of carbon core [Msun]
* o_core_mass -- mass of oxygen core [Msun]
* si_core_mass -- mass of silicon core [Msun]
* fe_core_mass -- mass of iron core [Msun]
* he_core_radius -- radius of helium core [Rsun]
* c_core_radius -- radius of carbon core [Rsun]
* o_core_radius -- radius of oxygen core [Rsun]
* si_core_radius -- radius of silicon core [Rsun]
* fe_core_radius -- radius of iron core [Rsun]
* max_abs_v_velocity -- maximum absolute velocity
* surf_avg_omega_div_omega_crit -- surface average rotation angular frequency
[Omega_crit]
* log_total_angular_momentum -- log10(total angular momentum [cm^2 g/s]
* surf_avg_omega -- surface average rotation angular
frequency [rad/s]
* surf_avg_v_rot -- surface average rotation velocity [km/s]
* star_mdot -- mass-loss rate [Msun/year]
"""
return __read_data(filename)
[docs]
def read_profile(filename):
"""Read data from a MESA-Web profile file
Parameters
----------
filename: string
name of profile file
Returns
-------
:class:`astropy.Table` containing header and profile data
>>> table = read_profile("profile123.data")
>>> pressure = table["pressure"]
Notes
-----
**Header data**: The following keys/value pairs in the returned Table.meta dict
contain header data -- i.e., scalars describing position-independent
properties of the star. Where applicable, units are given in square
brackets [].
* star_mdot -- mass-loss rate [Msun/year]
* model_number -- model number
* num_zones -- number of zones
* initial_mass -- initial mass [Msun]
* initial_z -- initial metal mass fraction
* star_age -- stellar age [years]
* time_step -- current time-step [s]
* Teff -- effective temperature [K]
* photosphere_L -- photospheric luminosity [Lsun]
* photosphere_r -- photospheric radius [Rsun]
* center_eta -- center electron chemical potential [kB*T]
* center_h1 -- center 1H mass fraction
* center_he3 -- center 3He mass fraction
* center_he4 -- center 4He mass fraction
* center_c12 -- center 12C mass fraction
* center_n14 -- center 14N mass fraction
* center_o16 -- center 16O mass fraction
* center_ne20 -- center 20Ne mass fraction
* star_age -- stellar age [years]
* star_mass -- stellar mass [Msun]
* star_mdot -- mass-loss rate [Msun/year]
* star_mass_h1 -- total mass in 1H [Msun]
* star_mass_he3 -- total mass in 3He [Msun]
* star_mass_he4 -- total mass in 4He [Msun]
* star_mass_c12 -- total mass in 12C [Msun]
* star_mass_n14 -- total mass in 14N [Msun]
* star_mass_o16 -- total mass in 16O [Msun]
* star_mass_ne20 -- total mass in 20Ne [Msun]
* he_core_mass -- mass of helium core [Msun]
* c_core_mass -- mass of carbon core [Msun]
* o_core_mass -- mass of oxygen core [Msun]
* si_core_mass -- mass of silicon core [Msun]
* fe_core_mass -- mass of iron core [Msun]
* neutron_rich_core_mass -- mass of neutron-rich core [Msun]
* tau10_mass -- mass coordinate of optical depth 10 [Msun]
* tau10_radius -- radius coordinate of optical depth 10 [Rsun]
* tau100_mass -- mass coordinate of optical depth 100 [Msun]
* tau100_radius -- radius coordinate of optical depth 100 [Rsun]
* dynamic_timescale -- dynamical timescale [s]
* kh_timescale -- Kelvin-Helmholtz timescale [s]
* nuc_timescale -- nuclear timescale [s]
* log_LH -- log10(total H-burning luminosity, excluding
neutrinos [Lsun])
* log_LHe -- log10(total He-burning luminosity, excluding
neutrinos [Lsun])
* log_LZ -- log10(total metal-burning luminosity, excluding
neutrinos [Lsun])
* power_nuc_burn -- total nuclear burning luminosity, excluding
photodisintegrations [Lsun]
* power_h_burn -- total H-burning luminosity, excluding neutrinos [Lsun]
* power_he_burn -- total He-burning luminosity, excluding neutrinos [Lsun]
* power_neu -- total neutrino luminosity [Lsun]
* burn_min1 -- 1st limit for reported burning [erg/g/s]
* burn_min2 -- 2nd limit for reported burning [erg/g/s]
**Profile data**: The following keys/value pairs in the returned Table contain
profile data -- i.e., describing local properties of the star over a
sequence of spatial zones. Where applicable, units are given in
square brackets [].
* mass -- mass coordinate [Msun]
* radius -- radius coordinate [Rsun]
* luminosity -- luminosity [Lsun]
* pressure -- pressure [dyn/cm^2]
* logRho -- log10(density [g/cm^3])
* logT -- log10(temperature [K])
* energy -- log10(specific internal energy [erg/g])
* entropy -- log10(specific entropy [kB*N_avo])
* cp -- specific heat at constant pressure [erg/K/g]
* gamma1 -- first adiabatic exponent
* grada -- adiabatic temperature gradient
* mu -- mean molecular weight
* free_e -- free electrons per nucleon
* ye -- average charge per baryon [e]
* pgas -- gas pressure [dyn/cm^2]
* prad -- radiation pressure [dyn/cm^2]
* gradr -- radiative temperature gradient
* gradT -- physical temperature gradient
* velocity -- velocity [km/s]
* conv_vel -- convective velocity [km/s]
* opacity -- opacity [cm^2/g]
* eps_nuc -- nuclear energy release rate, excluding neutrinos [erg/s/g]
* pp -- pp energy release rate [erg/s/g]
* cno -- CNO energy release rate [erg/s/g]
* tri_alfa -- triple-alpha energy release rate [erg/s/g]
* eps_nuc_neu_total -- energy loss rate as nuclear neutrinos [erg/s/g]
* non_nuc_neu -- energy loss rate as non-nuclear neutrinos [erg/s/g]
* eps_grav -- thermal energy release rate [erg/s/g]
* h1 -- 1H mass fraction
* he3 -- 3He mass fraction
* he4 -- 4He mass fraction
* c12 -- 12C mass fraction
* n14 -- 14N mass fraction
* o16 -- 16O mass fraction
* ne20 -- 20Ne mass fraction
* mg24 -- 24Mg mass fraction
* si28 -- 28Si mass fraction
* s32 -- 32S mass fraction
* ar36 -- 36Ar mass fraction
* ca40 -- 40Ca mass fraction
* ti44 -- 44Ti mass fraction
* cr48 -- 48Cr mass fraction
* fe52 -- 52Fe mass fraction
* fe54 -- 54Fe mass fraction
* fe56 -- 56Fe mass fraction
* ni56 -- 56Ni mass fraction
* eta -- electron chemical potential [kB*T]
* log_omega -- log10(rotation angular velocity [rad/s])
* v_rot -- rotation velocity [km/s]
* j_rot -- specific angular momentum [cm^2/s]
* dynamo_log_B_r -- log10(dynamo-generated radial field strength [gauss])
* dynamo_log_B_phi -- log10(dynamo-generated azimuthal field strength [gauss])
* log_D_conv -- log10(convective diffusivity [cm^2/s])
* log_D_semi -- log10(semiconvective diffusivity [cm^2/s])
* log_D_ovr -- log10(overshoot diffusivity [cm^2/s])
* log_D_thrm -- log10(thermohaline diffusivity [cm^2/s])
"""
return __read_data(filename, rev=True)
# Find and read a MESA-Web profile file
[docs]
def find_read_profile(filename, model_number, nearest=False):
"""Find the MESA-Web profile file corresponding to a given model number, and
read data from it
Parameters
----------
filename: string
name of profile index file (usually, 'profiles.index')
model_number: int
model number
nearest: bool
flag indicating whether the profile with the nearest
model number should be read, IF the exact model can't be
found (default False)
Returns
-------
:class:`astropy.Table` object containing header and profile data (see
read_profile() for details)
"""
# Open the index file and determine the profile file mapping
mod_num, prof_num = np.loadtxt(
filename, skiprows=1, unpack=True, usecols=(0, 2), dtype=int
)
# Find the closest model number
i = np.argmin(np.abs(mod_num - model_number))
if not nearest:
if mod_num[i] != model_number:
raise Exception(
"A model with the desired model number could not be found."
+ "Try using the 'nearest' flag"
)
return read_profile(op.join(op.dirname(filename), f"profile{prof_num[i]}.data"))
def __read_data(filename, rev=False):
with open(filename) as file:
# Read header data
def __num(s):
try:
return int(s)
except ValueError:
sr = re.sub(r"D", "E", s)
sr = re.sub(r"(\d)([+-])", r"\1E\2", sr)
try:
return float(sr)
except ValueError:
return s
file.readline()
header_names = file.readline().split()
header_values = [__num(value) for value in file.readline().split()]
# Read array data
file.readline()
file.readline()
column_names = file.readline().split()
lines_values = []
while True:
line = file.readline()
if not line:
break
lines_values.append([__num(value) for value in line.split()])
# Create data structure
data = Table(meta=dict(zip(header_names, header_values)))
# Populate the structure
if rev:
for i in range(0, len(column_names)):
data[column_names[i]] = np.array(
[line_values[i] for line_values in reversed(lines_values)]
)
else:
for i in range(0, len(column_names)):
data[column_names[i]] = np.array(
[line_values[i] for line_values in lines_values]
)
return data