Astronomical Model Comparison Scripts

This collection of Python functions is designed for comparing different astronomical models by generating comparison charts. It focuses on evolutionary tracks, solar calibration, luminosity-radius tracks, and opacity comparisons of stellar models, using matplotlib for visualization and tqdm for progress indication.

from pysep.dsep.utils import interpolate_model_to_age
from pysep.dsep import load_model
from pysep.dm.filetypes import lowtempopac_file
import os
 
import matplotlib.pyplot as plt
 
from tqdm import tqdm
 
def evolutionary_track_comparision():
    fergModel = load_model("./fergOutput/model.dsep")
    aesoModel1 = load_model("./aesopusOutput/model.dsep")
    aesoModel2 = load_model("./aesopusExtendOutput/model.dsep")
 
    fig, ax = plt.subplots(1,1,figsize=(10,7))
    ax.plot(fergModel['iso']['Log_T'], fergModel['iso']['Log_g'], color='black', label='ferg GS98')
    ax.plot(aesoModel1['iso']['Log_T'], aesoModel1['iso']['Log_g'], color='green', label='aesopus GS98 (Glue)')
    ax.plot(aesoModel2['iso']['Log_T'], aesoModel2['iso']['Log_g'], color='red', linestyle='dashed', label='aesopus GS98 (Extend)')
    ax.invert_xaxis()
    ax.invert_yaxis()
    ax.legend(loc='best')
 
    ax.set_xlabel("Teff")
    ax.set_ylabel("Log g")
 
    plt.savefig("ModelTrackComparison.pdf", bbox_inches="tight")
 
def solar_calib_evolutionary_track_comparision():
    fergModel = load_model("./fergCalibOutput/model.dsep")
    aesoModel = load_model("./aesopusCalibOutput/model.dsep")
 
    fergModelAtAge = interpolate_model_to_age(fergModel, 4.57)
    aesoModelAtAge = interpolate_model_to_age(aesoModel, 4.57)
 
    fig, ax = plt.subplots(1,1,figsize=(10,7))
    ax.plot(fergModel['track']['log_Teff'], fergModel['track']['log_g'], color='black', label='ferg GS98')
    ax.plot(aesoModel['track']['log_Teff'], aesoModel['track']['log_g'], color='green', label='aesopus GS98')
    ax.plot(fergModelAtAge['log_Teff'], fergModelAtAge['log_g'], marker='*', markersize=4, color='orange')
    ax.plot(aesoModelAtAge['log_Teff'], aesoModelAtAge['log_g'], marker='*', markersize=4, color='blue')
    ax.invert_xaxis()
    ax.invert_yaxis()
    ax.legend(loc='best')
 
    ax.set_xlabel("Teff")
    ax.set_ylabel("Log g")
 
    plt.savefig("SolarCalibratedModelTrackComparison.pdf", bbox_inches="tight")
 
def solar_calib_LR_track_comparision():
    fergModel = load_model("./fergCalibOutput/model.dsep")
    aesoModel = load_model("./aesopusCalibOutput/model.dsep")
 
    fergModelAtAge = interpolate_model_to_age(fergModel, 4.57)
    aesoModelAtAge = interpolate_model_to_age(aesoModel, 4.57)
 
    fig, ax = plt.subplots(1,1,figsize=(10,7))
    ax.plot(fergModel['track']['log_L'], fergModel['track']['log_R'], color='black', label='ferg GS98')
    ax.plot(aesoModel['track']['log_L'], aesoModel['track']['log_R'], color='green', label='aesopus GS98')
    ax.plot(fergModelAtAge['log_L'], fergModelAtAge['log_R'], marker='*', markersize=4, color='orange')
    ax.plot(aesoModelAtAge['log_L'], aesoModelAtAge['log_R'], marker='*', markersize=4, color='blue')
    ax.invert_xaxis()
    ax.invert_yaxis()
    ax.legend(loc='best')
 
    ax.set_xlabel("log(L)")
    ax.set_ylabel("Log(R)")
 
    plt.savefig("SolarCalibratedLRComparison.pdf", bbox_inches="tight")
 
def opacity_comparison():
    FergLowTempRoot = "./inputs/opac/low/ferg04"
    AesoLowTempRoot1 = "./inputs/opac/low/aesopus"
    AesoLowTempRoot2 = "./inputs/opac/low/aesopusExtend"
 
    FergFiles = sorted(os.listdir(FergLowTempRoot), key=lambda x: float(f"0.{x.split('.')[1]}"))
    AesoFiles1 = sorted(os.listdir(AesoLowTempRoot1), key=lambda x: float(f"0.{x.split('.')[1]}"))
    AesoFiles2 = sorted(os.listdir(AesoLowTempRoot2), key=lambda x: float(f"0.{x.split('.')[1]}"))
 
    FergPaths = [os.path.abspath(os.path.join(FergLowTempRoot, x)) for x in FergFiles]
    AesoPaths1 = [os.path.abspath(os.path.join(AesoLowTempRoot1, x)) for x in AesoFiles1]
    AesoPaths2 = [os.path.abspath(os.path.join(AesoLowTempRoot2, x)) for x in AesoFiles2]
 
    ferg = [lowtempopac_file(x) for x in FergPaths]
    aeso1 = [lowtempopac_file(x) for x in AesoPaths1]
    aeso2 = [lowtempopac_file(x) for x in AesoPaths2]
 
    LogR = ferg[0].LogR
    LogT = ferg[0].LogT
 
    opacFigDir = "./opacFigs"
    if not os.path.exists(opacFigDir):
        os.mkdir(opacFigDir)
 
    for f, a in tqdm(zip(ferg, aeso1), total=len(ferg)):
        assert f.X == a.X
        for ft, at in tqdm(zip(f, a), total=15, leave=False):
            for ridx, R in tqdm(enumerate(LogR), total=len(LogR), leave=False):
                fname = os.path.join(opacFigDir, f"opacCompare_Glue_{R}_{f.X}_{ft['Z']}.pdf")
 
                fig, ax = plt.subplots(1, 1, figsize=(10,7))
                ax.plot(LogT, at['Kappa'][:, ridx], color='red', label="Aesopus")
                ax.plot(LogT, ft['Kappa'][:, ridx], color='black', label="Ferg")
 
                ax.set_xlabel("Log T")
                ax.set_ylabel("Kappa")
                ax.legend(loc='best')
 
                plt.savefig(fname, bbox_inches='tight')
                plt.close(fig)
    for f, a in tqdm(zip(ferg, aeso2), total=len(ferg)):
        assert f.X == a.X
        for ft, at in tqdm(zip(f, a), total=15, leave=False):
            for ridx, R in tqdm(enumerate(LogR), total=len(LogR), leave=False):
                fname = os.path.join(opacFigDir, f"opacCompare_Extend_{R}_{f.X}_{ft['Z']}.pdf")
 
                fig, ax = plt.subplots(1, 1, figsize=(10,7))
                ax.plot(LogT, at['Kappa'][:, ridx], color='red', label="Aesopus")
                ax.plot(LogT, ft['Kappa'][:, ridx], color='black', label="Ferg")
 
                ax.set_xlabel("Log T")
                ax.set_ylabel("Kappa")
                ax.legend(loc='best')
 
                plt.savefig(fname, bbox_inches='tight')
                plt.close(fig)
 
if __name__ == "__main__":
   evolutionary_track_comparision()
   solar_calib_evolutionary_track_comparision()
   solar_calib_LR_track_comparision()
   # opacity_comparison()