This Python script serves as a wrapper around a Fortran program named 'newpoly.py', which is used to generate polytropic models of stars. The script includes functionalities to handle command line arguments, format input for the Fortran program, and manage the file environment for the program's execution. It provides features such as estimating effective temperature and luminosity based on input parameters like metal mass fractions and solar mass units, setting up the environment for the Fortran program to run by creating symlinks to specified input and output files, and invoking the Fortran program as a subprocess. The script is designed to facilitate the interaction with the 'newpoly' program by abstracting parameter file formatting and environmental setup tasks, making it easier for users to generate polytropic star models by providing simple command line arguments.

# runNewpoly.py
# Thomas M. Boudreaux
# Last Modified: May 10th 2021
# python wrapper around newpoly.py fortran file
#  newpoly generates a polytropic model of a star
#  this program takes in command line arguments and formates
#  the namelist file for that program so that it is easier to call
import generateNML as gNML
import numpy as np
import argparse
import os
 
 
def setup_env(inputFile, outputFile):
    """setup_env
    desc: setup the enviroment for the fortran program newpoly to run
    Inputs:
        - inputFile: (str) path to the input namelist file
        - outputFile: (str) path to file to save the polytropic model to
    PreState:
        N/A
    PostState:
        - If fort.11 and fort.12 exist in the local directory remove them
        - Symlink the inputfile to fort.11
        - Symlink the outputfile to fort.12
    Returns: bool
        - True if both symlinks were created sucsessfully
        - False otherwise
    """
    if os.path.islink("fort.11"):
        os.remove("fort.11")
    if os.path.islink("fort.12"):
        os.remove("fort.12")
 
    os.symlink(inputFile, "fort.11")
    os.symlink(outputFile, "fort.12")
 
    if os.path.islink("fort.11") and os.path.islink("fort.12"):
        return True
    else:
        return False
 
def estimate_teff_lum(zz, zsolar, mixlen, sumass):
    """estimate_teff_lum
    desc: empirically calibrated relations between metal mass fraction,
          mixing length, mass and Temperature + Luminosity
    Inputs:
        - zz: (float) Metal Mass Fraction
        - zsolar: (float) Metal mass fraction of the sun
        - mixlen: (float) Mixing length
        - sumass: (float) Mass in solar mass units
    PreState:
        N/A
    PostState:
        N/A
    Returns (teff: float, lum: float):
        - (float) teff - estimate log effective temperature in solar units
        - (float) lum - estimated log luminosity in solar units
 
    """
    # log(zz/zsolar) ~> [Fe/H]
    Tadjust = -0.010*np.log(zz/zsolar)/np.log(10) + 0.1*(mixlen-1.9)
    Ladjust = -0.015*np.log(zz/zsolar)/np.log(10)
 
    # Handel stars of mass less than 0.39 solar masses differently
    #  this is about the fully convective boundary
    if (sumass <= .390):
        teff = 0.430 * sumass + 3.410
        lum = 6.20 * sumass - 1.02
    else:
        teff = 3.64 + (4.0 - sumass) * 0.01
        lum = 2.0 + 0.4 * (sumass - 2.0)
    teff = teff + Tadjust
    lum = lum + Ladjust
    return teff, lum
 
 
def run_program(prog):
    """run_program
    desc: call a subprocess of some program
    Inputs:
        - prog: (str) some path to some program to instantiate as a subprocess
    PreState:
        N/A
    PostState:
        - Prog has been run
    Returns:
        N/A
    """
    os.popen(prog)
 
if __name__ == "__main__":
    parser = argparse.ArgumentParser(description="Generate name list files for "
                                    "the newpoly program incoperating correct "
                                    "isotropic abundances")
    parser.add_argument("prog", help="path to new poly program", type=str)
    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("--sumass", help="mass in solar mass units",
                        type=float, default=1.0)
    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("--polyOutput", help="path to save nml file too",
                        default='poly.nml')
    parser.add_argument("--compName",
                        help="Composition Name to add to output file",
                        default="comp", type=str)
 
    args = parser.parse_args()
 
    teff, lum = estimate_teff_lum(args.Z, 0.0134, args.CMIXL, args.sumass)
 
    eS, eES = gNML.format_elem(args)
    polynml = gNML.generate_nml(args, eS, eES, args.polyOutput, teff, lum)
 
    outputFile = f"m{int(args.sumass*100)}{args.compName}.prems"
    envOK = setup_env(args.polyOutput, outputFile)
    if envOK:
        run_program(args.prog)
    else:
        raise OSError("Enviroment not setup correctly! Unit files cannot be found!")