Script to Generate Namelist Files for Stellar Modelling

This Python script automates the generation of namelist files for stellar modelling programs. It processes isotropic abundances from composition files and incorporates these into the required namelist format. The script accepts various arguments related to the stellar composition and outputs a formatted .nml file, which is used by the newpoly program to model stellar properties such as mass, temperature, and luminosity.

# Thomas M. Boudreaux
# April 2021
# Generate namelist files for the newpoly file from 
#  composition files
import argparse
import os
import f90nml
 
def format_percents(args):
    # Taken from Table 3 of 
    # Asplund, M.; Grevesse, N.; Sauval, A.J.; Scott, P. The chemical composition of the Sun. Ann. Rev. Astron. Astrophys.2009, 47, 481–522.
    isoPercents = {
            'He3':(0, 1-(args.X+args.Z), 0.0166),
            'C12':(1, args.C, 98.8938),
            'C13':(2, args.C, 1.1062),
            'N14':(3, args.N, 99.771),
            'N15':(4, args.N, 0.229),
            'O16':(5, args.O, 99.7621),
            'O17':(6, args.O, 0.0379),
            'O18':(7, args.O, 0.2000)
            }
    isoExtendPercents = {
            'D':(0, args.X, 0.002)
            }
    if args.Li:
        isoExtendPercents['Li6'] = (1, args.Li, 7.59)
        isoExtendPercents['Li7'] = (2, args.Li, 92.41)
    if args.Be:
        isoExtendPercents['Be9'] = (3, args.Be, 100)
 
    return isoPercents, isoExtendPercents
 
 
def format_elem(args):
    # Works in version of python where dicts are not guarenteed to be ordered
    iP, iEP = format_percents(args)
 
    elem = [(iP[x][0], iP[x][1]*iP[x][2]/100) for x in iP]
    extenElem = [(iEP[x][0], iEP[x][1]*iEP[x][2]/100) for x in iEP]
 
    elemSorted = map(lambda x: x[1], sorted(elem, key=lambda x: x[0]))
    extenElemSorted = map(lambda x: x[1], sorted(extenElem, key=lambda x: x[0]))
 
    return list(elemSorted), list(extenElemSorted)
 
 
def generate_nml(args, elem, extenElem, output, teff, lum):
    polynml = {
            'DATA':{
                'SUMASS': args.sumass,
                'TEFFL1': teff,
                'SULUML': lum,
                "X": args.X,
                "Z": args.Z,
                "ELEM": elem,
                "EXTELEM": extenElem,
                "CMIXL": args.CMIXL,
                "BETA": args.BETA,
                "FMASS1": args.FMASS1,
                "FMASS2": args.FMASS2,
                "DDAGE": args.DDAGE,
                "PN": args.PN,
                "LEXCOM": not(args.NLEXCOM)
                }
            }
    if len(polynml['DATA']['EXTELEM']) != 4:
        polynml['DATA'].pop('EXTELEM', None)
    with open(output, 'w') as nml_file:
        f90nml.write(polynml, nml_file)
    return polynml
 
 
if __name__ == "__main__":
    parser = argparse.ArgumentParser(description="Generate name list files for "
                                    "the newpoly program incoperating correct "
                                    "isotropic abundances")
    parser.add_argument("X", help="Hydrogen Mass Fraction", type=float)
    parser.add_argument("Z", help="Metal Mass Fraction", type=float)
    parser.add_argument("C", help="Carbon Metal Mass Fraction", type=float)
    parser.add_argument("N", help="Nitrogen Metal Mass Fraction", type=float)
    parser.add_argument("O", help="Oxygen Metal Mass Fraction", type=float)
    parser.add_argument("-Li", help="Lithium Metal Mass Fraction", type=float)
    parser.add_argument("-Be", help="Berylium Metal Mass Fraction", type=float)
    parser.add_argument("--summass", help="mass in solar mass units",
                        type=float, default=1.0)
    parser.add_argument("--teffl1", help="effective temperature in T3 K",
                        type=float, default=3.60206)
    parser.add_argument("--suluml", help="Integrated luminosity", type=float,
                        default=1.5)
    parser.add_argument("--CMIXL", help="Mixing length", type=float,
                        default=1.9126)
    parser.add_argument("--BETA", help="Beta", type=float, default=1.0)
    parser.add_argument("--FMASS1", help="FMASS1", type=float, default=9.9e-6)
    parser.add_argument("--FMASS2", help="FMASS2", type=float, default=0.9999)
    parser.add_argument("--DDAGE", help="Age?", type=float, default=10000)
    parser.add_argument("--PN", help="PN", type=float, default=1.5)
    parser.add_argument("--NLEXCOM", help="! LEXCOM", action="store_false")
    parser.add_argument("-o", "--output", help="path to save nml file too",
                        default='poly.nml')
 
    args = parser.parse_args()
 
    eS, eES = format_elem(args)
    polynml = generate_nml(args, eS, eES, args.output)