Updating non-LTE, LTE, and AGB stellar libraries

beast/config.py

@@ -34,7 +34,9 @@
     vega="vega.hd5",
     filters="filters.hd5",
     kurucz04="kurucz2004.grid.fits",
+    bosz24="bosz2024.grid.fits",
     tlusty09="tlusty.lowres.grid.fits",
+    tlusty25="tlusty2025.grid.fits",
     hstcovar="hst_whitedwarf_frac_covar.fits",
     mist1="MIST_FeH-4.00_vvcrit0.4.fits",
     mist2="MIST_FeH-3.50_vvcrit0.4.fits",

beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py

@@ -0,0 +1,62 @@
+import argparse
+import numpy as np
+import matplotlib.pyplot as plt
+
+from beast.physicsmodel.stars import stellib
+from helpers import get_stellib_boundaries
+
+if __name__ == "__main__":  # pragma: no cover
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        "--grid",
+        default="BOSZ2024",
+        choices=["Tlusty2025", "BOSZ2024", "Aringer2016", "Kurucz", "Tlusty"],
+        help="Grid to show",
+    )
+    parser.add_argument("--png", action="store_true", help="save figures to png files")
+    args = parser.parse_args()
+
+    if args.grid == "BOSZ2024":
+        slib = stellib.BOSZ2024()
+    elif args.grid == "Tlusty2025":
+        slib = stellib.Tlusty2025()
+    elif args.grid == "Aringer2016":
+        slib = stellib.Aringer2026()
+    elif args.grid == "Kurucz":
+        slib = stellib.Kurucz()
+    elif args.grid == "Tlusty":
+        slib = stellib.Tlusty()
+
+    # useful when new grids are generarted
+    # results are used to define the bounding box in a model definition
+    bound = get_stellib_boundaries(slib, dlogT=0.0, dlogg=0.0)
+    print("Algorithmically generated bounds")
+    print(bound)
+
+    # setup the plots
+    fontsize = 12
+    font = {"size": fontsize}
+    plt.rc("font", **font)
+    plt.rc("lines", linewidth=2)
+    plt.rc("axes", linewidth=2)
+    plt.rc("xtick.major", width=2)
+    plt.rc("ytick.major", width=2)
+
+    for cz in np.unique(slib.Z.data):
+        gvals = slib.Z.data == cz
+        plt.plot(slib.logT[gvals], slib.logg[gvals], "ko")
+        plt.plot(bound[:, 0], bound[:, 1], "g-", label="algorithmly generated")
+        slib.plot_boundary(
+            dlogT=0.0, dlogg=0.0, label="defined in class def", alpha=0.5, color="b"
+        )
+        plt.title(f"{args.grid}; z = {cz:.2e}")
+        plt.xlabel("log(Teff)")
+        plt.ylabel("log(g)")
+        plt.legend()
+
+        if args.png:
+            plt.savefig(f"slib_coverage_{args.grid}_z_{cz:.2e}.png")
+            plt.close()
+        else:
+            plt.show()

beast/physicsmodel/stars/stellar_libraries/comp_tlusty25_bosz24.py

@@ -0,0 +1,68 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import astropy.units as u
+from astropy.io import ascii
+
+from helpers import rebin_spectrum
+
+if __name__ == "__main__":  # pragma: no cover
+
+    solar_z = 0.02
+
+    # read the kurucz spectrum
+    # wavelength first
+    wave_filename = "/home/kgordon/Python/extstar_data/Models/BOSZ2024/r2000/bosz2024_wave_r2000.txt"
+    mspec_lte_wave = ascii.read(
+        wave_filename,
+        format="no_header",
+        # fast_reader={"exponent_style": "D"},
+        names=["Wave"],
+    )
+
+    model_filename = "/home/kgordon/Python/extstar_data/Models/BOSZ2024/r2000/m+0.00/bosz2024_ap_t15000_g+4.0_m+0.00_a+0.00_c+0.00_v2_r2000_resam.txt.gz"
+    mspec_lte = ascii.read(
+        model_filename,
+        format="no_header",
+        # fast_reader={"exponent_style": "D"},
+        names=["SFlux", "SCont"],
+    )
+    mspec_lte["Wave"] = mspec_lte_wave["Wave"] * u.angstrom
+
+    # read the tlusty spectrum
+    model_filename = (
+        "/home/kgordon/Python/extstar_data/Models/Tlusty_2023/z100t15000g400v2.spec.gz"
+    )
+    # read in the model spectrum
+    mspec = ascii.read(
+        model_filename,
+        format="no_header",
+        fast_reader={"exponent_style": "D"},
+        names=["Wave", "SFlux"],
+    )
+
+    # convert the type to float
+    mspec["SFlux"] = mspec["SFlux"].astype(float)
+
+    # set the units
+    fluxunit = u.erg / (u.s * u.cm * u.cm * u.angstrom)
+    mspec["Wave"].unit = u.angstrom
+    mspec["SFlux"].unit = fluxunit
+
+    # now extract the wave and flux colums
+    mwave = mspec["Wave"]
+    mflux = mspec["SFlux"]
+
+    # rebin to R=4000 for speed
+    rbres = 4000.0
+    wave_rebin, flux_rebin, npts_rebin = rebin_spectrum(
+        mwave.value, mflux.value, rbres, [200.0, 310000.0]
+    )
+    wave_rebin *= u.angstrom
+    print(len(wave_rebin))
+
+    # plot
+    plt.plot(wave_rebin.to(u.micron), flux_rebin * 4 * np.pi)
+    plt.plot(mspec_lte["Wave"].to(u.micron), mspec_lte["SFlux"] * 4 * np.pi)
+    plt.yscale("log")
+    plt.xscale("log")
+    plt.show()