From 3aed43262add254b2aae77882265128dacee0248 Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Tue, 28 Apr 2026 11:54:37 -0400
Subject: [PATCH 01/32] 1st draft of code to generate the tlusty 2025 stellar
 library

---
 .../stellar_libraries/comp_tlusty25_bosz24.py | 137 ++++++++++++++++++
 .../stars/stellar_libraries/gen_tlusty25.py   |  95 ++++++++++++
 .../stars/stellar_libraries/helpers.py        |  71 +++++++++
 3 files changed, 303 insertions(+)
 create mode 100644 beast/physicsmodel/stars/stellar_libraries/comp_tlusty25_bosz24.py
 create mode 100644 beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
 create mode 100644 beast/physicsmodel/stars/stellar_libraries/helpers.py

diff --git a/beast/physicsmodel/stars/stellar_libraries/comp_tlusty25_bosz24.py b/beast/physicsmodel/stars/stellar_libraries/comp_tlusty25_bosz24.py
new file mode 100644
index 00000000..36c6b8c1
--- /dev/null
+++ b/beast/physicsmodel/stars/stellar_libraries/comp_tlusty25_bosz24.py
@@ -0,0 +1,137 @@
+from scipy.interpolate import interp1d
+import matplotlib.pyplot as plt
+import numpy as np
+from astropy.table import QTable
+from astropy.io import fits
+import astropy.units as u
+from astropy.io import ascii
+
+
+def rebin_spectrum(wave, flux, resolution, wave_range):
+    """
+    Rebin spectrum to input resolution over input wavelength range
+    High to lower resolution only
+
+    Parameters
+    ----------
+    wave : vector
+        wavelengths of spectrum
+
+    flux : vector
+        spectrum in flux units
+
+    resolution : float
+        resolution of output spectrum
+
+    wave_range : [float, float]
+        wavelength range of output spectrum
+
+    Outputs
+    ------
+    wave, flux, npts : tuple of vectors
+        the model wavelength, flux, and npts at the requested wavelength
+    """
+    npts = int(
+        np.log10(wave_range[1] / wave_range[0])
+        / np.log10((1.0 + 2.0 * resolution) / (2.0 * resolution - 1.0))
+    )
+
+    twave = np.logspace(
+        np.log10(wave_range[0]), np.log10(wave_range[1]), num=npts + 1, endpoint=True
+    )
+    full_wave_min = twave[0:-1]
+    full_wave_max = twave[1:]
+    full_wave = 0.5 * (full_wave_min + full_wave_max)
+
+    full_flux = np.zeros((npts))
+    full_npts = np.zeros((npts), dtype=int)
+
+    for k in range(npts):
+        (indxs,) = np.where((wave >= full_wave_min[k]) & (wave < full_wave_max[k]))
+        n_indxs = len(indxs)
+        if n_indxs > 0:
+            full_flux[k] = np.sum(flux[indxs])
+            full_npts[k] = n_indxs
+
+    # divide by the # of points to create the final rebinned spectrum
+    (indxs,) = np.where(full_npts > 0)
+    if len(indxs):
+        full_flux[indxs] = full_flux[indxs] / full_npts[indxs]
+
+    # interpolate to fill in missing points in the rebinned spectrum
+    #  e.g., the model spectrum is not computed at a high enough resolution
+    #        at all the needed wavelengths
+    (zindxs,) = np.where(full_npts <= 0)
+    if len(zindxs):
+        ifunc = interp1d(
+            full_wave[indxs], full_flux[indxs], kind="linear", bounds_error=False
+        )
+        full_flux = ifunc(full_wave)
+        full_npts[zindxs] = 1
+        (nanindxs,) = np.where(~np.isfinite(full_flux))
+        if len(nanindxs):
+            full_flux[nanindxs] = 0.0
+            full_npts[nanindxs] = 0
+
+    return (full_wave, full_flux, full_npts)
+
+
+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()
diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
new file mode 100644
index 00000000..f16d7d15
--- /dev/null
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
@@ -0,0 +1,95 @@
+import glob
+
+import numpy as np
+
+from astropy.table import QTable
+from astropy.io import fits
+import astropy.units as u
+from astropy.io import ascii
+
+from helpers import rebin_spectrum
+
+
+def decode_params(filename):
+    """
+    Decode the tlusty filenames for the model parameters
+    """
+    model_params = {}
+
+    slashpos = filename.rfind("/")
+    periodpos = filename.rfind(f".spec")
+
+    zpos = filename.find("z", slashpos)
+    tpos = filename.find("t", slashpos)
+    gpos = filename.find("g", slashpos)
+    vpos = filename.find("v", slashpos)
+
+    if tpos - zpos > 4:
+        model_params["Z"] = float(filename[zpos + 1 : tpos]) * 0.001
+    else:
+        model_params["Z"] = float(filename[zpos + 1 : tpos]) * 0.01
+    model_params["Teff"] = float(filename[tpos + 1 : gpos])
+    model_params["logg"] = float(filename[gpos + 1 : vpos]) * 0.01
+    model_params["vturb"] = float(filename[vpos + 1 : periodpos])
+
+    return model_params
+
+
+if __name__ == "__main__":  # pragma: no cover
+
+    solar_z = 0.02
+
+    # get all vtrub=2 models
+    path = "/home/kgordon/Python/extstar_data/Models/Tlusty_2023/"
+    files = glob.glob(f"{path}/*v2.spec.gz")
+    files = files[0:2]
+    n_files = len(files)
+
+    # fmt: off
+    outtab = QTable(
+        names=["Z", "Teff", "logg", "vturb"],
+        dtype=["f", "f", "f", "f"]
+    )
+    # fmt: on
+
+    for k, cfile in enumerate(files):
+
+        # read in the model spectrum
+        # wave units are angstrom
+        # flux units are H, F = 4*pi*H in ergs/(cm^2 s A)
+        mspec = ascii.read(
+            cfile,
+            format="no_header",
+            fast_reader={"exponent_style": "D"},
+            names=["Wave", "SFlux"],
+        )
+
+        # convert the type to float
+        mspec["SFlux"] = mspec["SFlux"].astype(float)
+
+        # now extract the wave and flux colums
+        mwave = mspec["Wave"]
+        mflux = mspec["SFlux"]
+
+        # rebin to R=4000 for speed, common res and wave range with BOSZ LTE models
+        rbres = 4000.0
+        wave_rebin, flux_rebin, npts_rebin = rebin_spectrum(
+            mwave.value, mflux.value, rbres, [500.0, 320000.0]
+        )
+
+        if k == 0:
+            outspec = np.zeros((len(wave_rebin), len(files) + 1))
+            outspec[:, -1] = wave_rebin
+        outspec[:, k] = flux_rebin
+
+        model_params = decode_params(cfile)
+        Z = model_params["Z"] * solar_z
+        row = (Z, model_params["Teff"], model_params["logg"], model_params["vturb"])
+        outtab.add_row(row)
+
+    # output the stellar library in the beast format
+    spec_hdu = fits.PrimaryHDU(data=outspec)
+    table_hdu = fits.BinTableHDU(data=outtab)
+    table_hdu.name = "TLUSTY"
+    hdul = fits.HDUList([spec_hdu, table_hdu])
+    hdul.writeto("tlusty2025.grid.fits", overwrite=True)
diff --git a/beast/physicsmodel/stars/stellar_libraries/helpers.py b/beast/physicsmodel/stars/stellar_libraries/helpers.py
new file mode 100644
index 00000000..4b469089
--- /dev/null
+++ b/beast/physicsmodel/stars/stellar_libraries/helpers.py
@@ -0,0 +1,71 @@
+from scipy.interpolate import interp1d
+import numpy as np
+
+
+def rebin_spectrum(wave, flux, resolution, wave_range):
+    """
+    Rebin spectrum to input resolution over input wavelength range
+    High to lower resolution only
+
+    Parameters
+    ----------
+    wave : vector
+        wavelengths of spectrum
+
+    flux : vector
+        spectrum in flux units
+
+    resolution : float
+        resolution of output spectrum
+
+    wave_range : [float, float]
+        wavelength range of output spectrum
+
+    Outputs
+    ------
+    wave, flux, npts : tuple of vectors
+        the model wavelength, flux, and npts at the requested wavelength
+    """
+    npts = int(
+        np.log10(wave_range[1] / wave_range[0])
+        / np.log10((1.0 + 2.0 * resolution) / (2.0 * resolution - 1.0))
+    )
+
+    twave = np.logspace(
+        np.log10(wave_range[0]), np.log10(wave_range[1]), num=npts + 1, endpoint=True
+    )
+    full_wave_min = twave[0:-1]
+    full_wave_max = twave[1:]
+    full_wave = 0.5 * (full_wave_min + full_wave_max)
+
+    full_flux = np.zeros((npts))
+    full_npts = np.zeros((npts), dtype=int)
+
+    for k in range(npts):
+        (indxs,) = np.where((wave >= full_wave_min[k]) & (wave < full_wave_max[k]))
+        n_indxs = len(indxs)
+        if n_indxs > 0:
+            full_flux[k] = np.sum(flux[indxs])
+            full_npts[k] = n_indxs
+
+    # divide by the # of points to create the final rebinned spectrum
+    (indxs,) = np.where(full_npts > 0)
+    if len(indxs):
+        full_flux[indxs] = full_flux[indxs] / full_npts[indxs]
+
+    # interpolate to fill in missing points in the rebinned spectrum
+    #  e.g., the model spectrum is not computed at a high enough resolution
+    #        at all the needed wavelengths
+    (zindxs,) = np.where(full_npts <= 0)
+    if len(zindxs):
+        ifunc = interp1d(
+            full_wave[indxs], full_flux[indxs], kind="linear", bounds_error=False
+        )
+        full_flux = ifunc(full_wave)
+        full_npts[zindxs] = 1
+        (nanindxs,) = np.where(~np.isfinite(full_flux))
+        if len(nanindxs):
+            full_flux[nanindxs] = 0.0
+            full_npts[nanindxs] = 0
+
+    return (full_wave, full_flux, full_npts)
\ No newline at end of file

From 787540888329c9ae43991177c39d21a1d68fa365 Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Tue, 28 Apr 2026 11:59:30 -0400
Subject: [PATCH 02/32] adding stellib entry

---
 beast/physicsmodel/stars/stellib.py | 86 +++++++++++++++++++++++++++++
 1 file changed, 86 insertions(+)

diff --git a/beast/physicsmodel/stars/stellib.py b/beast/physicsmodel/stars/stellib.py
index fd2e5740..9d1c1c3c 100644
--- a/beast/physicsmodel/stars/stellib.py
+++ b/beast/physicsmodel/stars/stellib.py
@@ -1943,3 +1943,89 @@ def Z(self):
     @property
     def logZ(self):
         return self.grid["logz"]
+
+
+class Tlusty2025(Stellib):
+    """
+    Tlusty 2025 O and B stellar atmospheres.  Updated in 2025 with better line lists
+    and extended IR wavelength coverage
+
+    * NLTE
+    * Plane Parallel
+    * line blanketing
+
+    References
+    ----------
+    Hubeny 1988 for initial reference
+    Lanz, T., & Hubeny, I. (2003) for NLTE developments
+    Hubeny, Bohlin, Gordon, & Fitzpatrick (2025) for updated models
+
+    files are available at: https://stsci.box.com/v/tlustyOB2025
+    """
+
+    def __init__(self, filename=None, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        self.name = "Tlusty2025"
+        if filename is None:
+            self.source = config["tlusty2025"]
+        else:
+            self.source = filename
+        self._load_()
+
+    def _load_(self):
+        g = SpectralGrid(self.source, backend="memory")
+        self.wavelength = g.lamb
+        self.grid = g.grid
+        self.header["NAME"] = "tlusty2025"
+        self.spectra = g.seds
+
+    def bbox(self, dlogT=0.05, dlogg=0.25):
+        """ Boundary of Tlusty library
+
+        Parameters
+        ----------
+        dlogT: float
+            log-temperature tolerance before extrapolation limit
+
+        dlogg: float
+            log-g tolerance before extrapolation limit
+
+        Returns
+        -------
+        bbox: ndarray
+            (logT, logg) edges of the bounding polygon
+        """
+        bbox = [
+            (4.176 - dlogT, 4.749 + dlogg),
+            (4.176 - dlogT, 1.750 - dlogg),
+            (4.176 + dlogT, 1.750 - dlogg),
+            (4.255 + dlogT, 2.000 - dlogg),
+            (4.447 + dlogT, 2.750 - dlogg),
+            (4.478 + dlogT, 3.000 - dlogg),
+            (4.544 + dlogT, 3.250 - dlogg),
+            (4.740 + dlogT, 4.000 - dlogg),
+            (4.740 + dlogT, 4.749 + dlogg),
+            (4.176 - dlogT, 4.749 + dlogg),
+        ]
+
+        return np.array(bbox)
+
+    @property
+    def logT(self):
+        return self.grid["logT"]
+
+    @property
+    def logg(self):
+        return self.grid["logg"]
+
+    @property
+    def Teff(self):
+        return self.grid["Teff"]
+
+    @property
+    def Z(self):
+        return self.grid["Z"]
+
+    @property
+    def logZ(self):
+        return np.log10(self.grid["Z"])

From 90d1069cde203294e1f27422263fac2de9a01788 Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Tue, 28 Apr 2026 12:36:39 -0400
Subject: [PATCH 03/32] code for checking the boundaries

---
 .../stars/stellar_libraries/check_tlusty25.py | 46 +++++++++++++++++++
 beast/physicsmodel/stars/stellib.py           |  6 +++
 2 files changed, 52 insertions(+)
 create mode 100644 beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py

diff --git a/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py b/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py
new file mode 100644
index 00000000..8d8e12d1
--- /dev/null
+++ b/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py
@@ -0,0 +1,46 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+from beast.physicsmodel.stars import stellib
+
+
+def get_stellib_boundaries(s, dlogT=0.1, dlogg=0.3, closed=True):
+    """Returns the closed boundary polygon around the stellar library with
+    given margins
+
+    INPUTS:
+        s   Stellib     Stellar library object
+
+    KEYWORDS:
+        dlogT   float       margin in logT
+        dlogg   float       margin in logg
+        closed  bool        if set, close the polygon
+
+    OUTPUTS:
+        b   ndarray[float, ndim=2]  (closed) boundary points: [logg, Teff]
+    """
+    #        use "points_inside_poly" to test wether a point is inside the limits
+    #        >>> data = np.array([iso.data['logg'], iso.data['logT']]).T
+    #        >>> aa = points_inside_poly(data, leftb)
+    #    """
+    leftb = [(k, np.max(s.logT[s.logg == k]) + dlogT) for k in np.unique(s.logg)]
+    leftb += [(leftb[-1][0] + dlogg, leftb[-1][1])]
+    leftb = [(leftb[0][0] - dlogg, leftb[0][1])] + leftb
+    rightb = [(k, np.min(s.logT[s.logg == k]) - dlogT) for k in np.unique(s.logg)[::-1]]
+    rightb += [(rightb[-1][0] - dlogg, rightb[-1][1])]
+    rightb = [(rightb[0][0] + dlogg, rightb[0][1])] + rightb
+    b = leftb + rightb
+    if closed:
+        b += [b[0]]
+    return np.array(b)
+
+
+slib = stellib.Tlusty2025()
+bound = get_stellib_boundaries(slib, dlogT=0.01, dlogg=0.01)
+bbox = np.array(slib.bbox())
+print(bbox[:, 0])
+
+plt.plot(slib.logT, slib.logg, "ko")
+plt.plot(bound[:, 1], bound[:, 0])
+plt.plot(bbox[:, 0], bbox[:, 1])
+plt.show()
diff --git a/beast/physicsmodel/stars/stellib.py b/beast/physicsmodel/stars/stellib.py
index 9d1c1c3c..7a81850e 100644
--- a/beast/physicsmodel/stars/stellib.py
+++ b/beast/physicsmodel/stars/stellib.py
@@ -34,6 +34,7 @@
     "btsettl_medres": __ROOT__ + "bt-settl.medres.grid.fits",
     "munari": __ROOT__ + "atlas9-munari.hires.grid.fits",
     "aringer": __ROOT__ + "Aringer.AGB.grid.fits",
+    "tlusty2025": __ROOT__ + "tlusty2025.grid.fits",
 }
 
 __all__ = [
@@ -46,6 +47,7 @@
     "Elodie",
     "BaSeL",
     "Aringer",
+    "Tlusty2025"
 ]
 
 
@@ -2022,6 +2024,10 @@ def logg(self):
     def Teff(self):
         return self.grid["Teff"]
 
+    @property
+    def logT(self):
+        return np.log10(self.grid["Teff"])
+
     @property
     def Z(self):
         return self.grid["Z"]

From 653b4e06459018fb2125d79a9a04fb1eec8fbb92 Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Tue, 28 Apr 2026 12:47:35 -0400
Subject: [PATCH 04/32] fixing various codestyle issues

---
 .../stars/stellar_libraries/check_tlusty25.py | 20 ++---
 .../stellar_libraries/comp_tlusty25_bosz24.py | 80 ++-----------------
 .../stars/stellar_libraries/gen_tlusty25.py   |  5 +-
 .../stars/stellar_libraries/helpers.py        |  2 +-
 beast/physicsmodel/stars/stellib.py           |  4 -
 5 files changed, 20 insertions(+), 91 deletions(-)

diff --git a/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py b/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py
index 8d8e12d1..6f15b904 100644
--- a/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py
+++ b/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py
@@ -35,12 +35,14 @@ def get_stellib_boundaries(s, dlogT=0.1, dlogg=0.3, closed=True):
     return np.array(b)
 
 
-slib = stellib.Tlusty2025()
-bound = get_stellib_boundaries(slib, dlogT=0.01, dlogg=0.01)
-bbox = np.array(slib.bbox())
-print(bbox[:, 0])
-
-plt.plot(slib.logT, slib.logg, "ko")
-plt.plot(bound[:, 1], bound[:, 0])
-plt.plot(bbox[:, 0], bbox[:, 1])
-plt.show()
+if __name__ == "__main__":  # pragma: no cover
+
+    slib = stellib.Tlusty2025()
+    bound = get_stellib_boundaries(slib, dlogT=0.01, dlogg=0.01)
+    bbox = np.array(slib.bbox())
+    print(bbox[:, 0])
+
+    plt.plot(slib.logT, slib.logg, "ko")
+    plt.plot(bound[:, 1], bound[:, 0])
+    plt.plot(bbox[:, 0], bbox[:, 1])
+    plt.show()
diff --git a/beast/physicsmodel/stars/stellar_libraries/comp_tlusty25_bosz24.py b/beast/physicsmodel/stars/stellar_libraries/comp_tlusty25_bosz24.py
index 36c6b8c1..c3796005 100644
--- a/beast/physicsmodel/stars/stellar_libraries/comp_tlusty25_bosz24.py
+++ b/beast/physicsmodel/stars/stellar_libraries/comp_tlusty25_bosz24.py
@@ -1,79 +1,9 @@
-from scipy.interpolate import interp1d
 import matplotlib.pyplot as plt
 import numpy as np
-from astropy.table import QTable
-from astropy.io import fits
 import astropy.units as u
 from astropy.io import ascii
 
-
-def rebin_spectrum(wave, flux, resolution, wave_range):
-    """
-    Rebin spectrum to input resolution over input wavelength range
-    High to lower resolution only
-
-    Parameters
-    ----------
-    wave : vector
-        wavelengths of spectrum
-
-    flux : vector
-        spectrum in flux units
-
-    resolution : float
-        resolution of output spectrum
-
-    wave_range : [float, float]
-        wavelength range of output spectrum
-
-    Outputs
-    ------
-    wave, flux, npts : tuple of vectors
-        the model wavelength, flux, and npts at the requested wavelength
-    """
-    npts = int(
-        np.log10(wave_range[1] / wave_range[0])
-        / np.log10((1.0 + 2.0 * resolution) / (2.0 * resolution - 1.0))
-    )
-
-    twave = np.logspace(
-        np.log10(wave_range[0]), np.log10(wave_range[1]), num=npts + 1, endpoint=True
-    )
-    full_wave_min = twave[0:-1]
-    full_wave_max = twave[1:]
-    full_wave = 0.5 * (full_wave_min + full_wave_max)
-
-    full_flux = np.zeros((npts))
-    full_npts = np.zeros((npts), dtype=int)
-
-    for k in range(npts):
-        (indxs,) = np.where((wave >= full_wave_min[k]) & (wave < full_wave_max[k]))
-        n_indxs = len(indxs)
-        if n_indxs > 0:
-            full_flux[k] = np.sum(flux[indxs])
-            full_npts[k] = n_indxs
-
-    # divide by the # of points to create the final rebinned spectrum
-    (indxs,) = np.where(full_npts > 0)
-    if len(indxs):
-        full_flux[indxs] = full_flux[indxs] / full_npts[indxs]
-
-    # interpolate to fill in missing points in the rebinned spectrum
-    #  e.g., the model spectrum is not computed at a high enough resolution
-    #        at all the needed wavelengths
-    (zindxs,) = np.where(full_npts <= 0)
-    if len(zindxs):
-        ifunc = interp1d(
-            full_wave[indxs], full_flux[indxs], kind="linear", bounds_error=False
-        )
-        full_flux = ifunc(full_wave)
-        full_npts[zindxs] = 1
-        (nanindxs,) = np.where(~np.isfinite(full_flux))
-        if len(nanindxs):
-            full_flux[nanindxs] = 0.0
-            full_npts[nanindxs] = 0
-
-    return (full_wave, full_flux, full_npts)
+from helpers import rebin_spectrum
 
 
 if __name__ == "__main__":  # pragma: no cover
@@ -86,7 +16,7 @@ def rebin_spectrum(wave, flux, resolution, wave_range):
     mspec_lte_wave = ascii.read(
         wave_filename,
         format="no_header",
-        #fast_reader={"exponent_style": "D"},
+        # fast_reader={"exponent_style": "D"},
         names=["Wave"],
     )
 
@@ -94,13 +24,15 @@ def rebin_spectrum(wave, flux, resolution, wave_range):
     mspec_lte = ascii.read(
         model_filename,
         format="no_header",
-        #fast_reader={"exponent_style": "D"},
+        # 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"
+    model_filename = (
+        "/home/kgordon/Python/extstar_data/Models/Tlusty_2023/z100t15000g400v2.spec.gz"
+    )
     # read in the model spectrum
     mspec = ascii.read(
         model_filename,
diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
index f16d7d15..609fee82 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
@@ -4,7 +4,6 @@
 
 from astropy.table import QTable
 from astropy.io import fits
-import astropy.units as u
 from astropy.io import ascii
 
 from helpers import rebin_spectrum
@@ -17,7 +16,7 @@ def decode_params(filename):
     model_params = {}
 
     slashpos = filename.rfind("/")
-    periodpos = filename.rfind(f".spec")
+    periodpos = filename.rfind(".spec")
 
     zpos = filename.find("z", slashpos)
     tpos = filename.find("t", slashpos)
@@ -80,7 +79,7 @@ def decode_params(filename):
         if k == 0:
             outspec = np.zeros((len(wave_rebin), len(files) + 1))
             outspec[:, -1] = wave_rebin
-        outspec[:, k] = flux_rebin
+        outspec[:, k] = flux_rebin * 4.0 * np.pi
 
         model_params = decode_params(cfile)
         Z = model_params["Z"] * solar_z
diff --git a/beast/physicsmodel/stars/stellar_libraries/helpers.py b/beast/physicsmodel/stars/stellar_libraries/helpers.py
index 4b469089..7386468c 100644
--- a/beast/physicsmodel/stars/stellar_libraries/helpers.py
+++ b/beast/physicsmodel/stars/stellar_libraries/helpers.py
@@ -68,4 +68,4 @@ def rebin_spectrum(wave, flux, resolution, wave_range):
             full_flux[nanindxs] = 0.0
             full_npts[nanindxs] = 0
 
-    return (full_wave, full_flux, full_npts)
\ No newline at end of file
+    return (full_wave, full_flux, full_npts)
diff --git a/beast/physicsmodel/stars/stellib.py b/beast/physicsmodel/stars/stellib.py
index 7a81850e..11a761bc 100644
--- a/beast/physicsmodel/stars/stellib.py
+++ b/beast/physicsmodel/stars/stellib.py
@@ -2012,10 +2012,6 @@ def bbox(self, dlogT=0.05, dlogg=0.25):
 
         return np.array(bbox)
 
-    @property
-    def logT(self):
-        return self.grid["logT"]
-
     @property
     def logg(self):
         return self.grid["logg"]

From 2a295905644628296e3533eeafb05e501ca4750d Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Tue, 28 Apr 2026 12:55:04 -0400
Subject: [PATCH 05/32] fixing docstrings

---
 .../stars/stellar_libraries/check_tlusty25.py | 32 +++++++++++--------
 .../stars/stellar_libraries/gen_tlusty25.py   |  7 +++-
 .../stars/stellar_libraries/helpers.py        |  2 +-
 3 files changed, 26 insertions(+), 15 deletions(-)

diff --git a/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py b/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py
index 6f15b904..12236f5e 100644
--- a/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py
+++ b/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py
@@ -5,19 +5,25 @@
 
 
 def get_stellib_boundaries(s, dlogT=0.1, dlogg=0.3, closed=True):
-    """Returns the closed boundary polygon around the stellar library with
-    given margins
-
-    INPUTS:
-        s   Stellib     Stellar library object
-
-    KEYWORDS:
-        dlogT   float       margin in logT
-        dlogg   float       margin in logg
-        closed  bool        if set, close the polygon
-
-    OUTPUTS:
-        b   ndarray[float, ndim=2]  (closed) boundary points: [logg, Teff]
+    """
+    Returns the closed boundary polygon around the stellar library with
+    given margins.
+
+    Parameters
+    ----------
+    s :  Stellib
+        Stellar library object
+    dlogT :  float
+        margin in logT
+    dlogg  : float
+        margin in logg
+    closed : bool
+        if set, close the polygon
+
+    Returns
+    ------
+    b : ndarray[float, ndim=2]
+        (closed) boundary points: [logg, Teff]
     """
     #        use "points_inside_poly" to test wether a point is inside the limits
     #        >>> data = np.array([iso.data['logg'], iso.data['logT']]).T
diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
index 609fee82..341cf327 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
@@ -11,7 +11,12 @@
 
 def decode_params(filename):
     """
-    Decode the tlusty filenames for the model parameters
+    Decode the tlusty filenames for the model parameters.
+
+    Parameters
+    ----------
+    filename : str
+        name of file to decode
     """
     model_params = {}
 
diff --git a/beast/physicsmodel/stars/stellar_libraries/helpers.py b/beast/physicsmodel/stars/stellar_libraries/helpers.py
index 7386468c..1e5dbf97 100644
--- a/beast/physicsmodel/stars/stellar_libraries/helpers.py
+++ b/beast/physicsmodel/stars/stellar_libraries/helpers.py
@@ -5,7 +5,7 @@
 def rebin_spectrum(wave, flux, resolution, wave_range):
     """
     Rebin spectrum to input resolution over input wavelength range
-    High to lower resolution only
+    High to lower resolution only.
 
     Parameters
     ----------

From 4ff2f81352661fa471368c886ace3ea49f167b71 Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Tue, 28 Apr 2026 13:00:35 -0400
Subject: [PATCH 06/32] updating path

---
 beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
index 341cf327..29a15749 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
@@ -44,9 +44,8 @@ def decode_params(filename):
     solar_z = 0.02
 
     # get all vtrub=2 models
-    path = "/home/kgordon/Python/extstar_data/Models/Tlusty_2023/"
+    path = "/home/kgordon/Python/extstar_data/Models/Tlusty_2025/"
     files = glob.glob(f"{path}/*v2.spec.gz")
-    files = files[0:2]
     n_files = len(files)
 
     # fmt: off
@@ -57,6 +56,7 @@ def decode_params(filename):
     # fmt: on
 
     for k, cfile in enumerate(files):
+        print(cfile)
 
         # read in the model spectrum
         # wave units are angstrom

From 693e1548f9e6586080918cea4ba8ac504cb768cb Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Fri, 1 May 2026 11:16:00 -0400
Subject: [PATCH 07/32] update tlusty file format and fix bug in stelllib

---
 .../stars/stellar_libraries/gen_tlusty25.py         | 13 +++++++------
 beast/physicsmodel/stars/stellib.py                 |  2 +-
 2 files changed, 8 insertions(+), 7 deletions(-)

diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
index 29a15749..03d14365 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
@@ -50,8 +50,8 @@ def decode_params(filename):
 
     # fmt: off
     outtab = QTable(
-        names=["Z", "Teff", "logg", "vturb"],
-        dtype=["f", "f", "f", "f"]
+        names=["Z", "Teff", "logg", "vturb", "logT", "logz"],
+        dtype=["f", "f", "f", "f", "f", "f"]
     )
     # fmt: on
 
@@ -82,13 +82,14 @@ def decode_params(filename):
         )
 
         if k == 0:
-            outspec = np.zeros((len(wave_rebin), len(files) + 1))
-            outspec[:, -1] = wave_rebin
-        outspec[:, k] = flux_rebin * 4.0 * np.pi
+            outspec = np.zeros((len(files) + 1, len(wave_rebin)))
+            outspec[-1, :] = wave_rebin
+        outspec[k, :] = flux_rebin * 4.0 * np.pi
 
         model_params = decode_params(cfile)
         Z = model_params["Z"] * solar_z
-        row = (Z, model_params["Teff"], model_params["logg"], model_params["vturb"])
+        row = (Z, model_params["Teff"], model_params["logg"], model_params["vturb"], np.log10(model_params["Teff"]),
+               np.log10(Z))
         outtab.add_row(row)
 
     # output the stellar library in the beast format
diff --git a/beast/physicsmodel/stars/stellib.py b/beast/physicsmodel/stars/stellib.py
index 11a761bc..ea6a974f 100644
--- a/beast/physicsmodel/stars/stellib.py
+++ b/beast/physicsmodel/stars/stellib.py
@@ -821,7 +821,7 @@ def wavelength(self):
         """ return a common wavelength sampling to all libraries. This can be
         used to reinterpolate any spectrum onto a common definition """
 
-        lambs = np.unique(np.asarray([osl.wavelength[:] for osl in self._olist]))
+        lambs = np.unique(np.concatenate([osl.wavelength[:] for osl in self._olist]))
         return lambs
 
     @property

From 4b47bc983b4bc9fd5e3d1942ebc08ee0fa3ea0fa Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Fri, 1 May 2026 11:30:14 -0400
Subject: [PATCH 08/32] fixing some codacy issues

---
 .../physicsmodel/stars/stellar_libraries/check_tlusty25.py | 7 ++-----
 beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py | 1 +
 beast/physicsmodel/stars/stellar_libraries/helpers.py      | 1 +
 3 files changed, 4 insertions(+), 5 deletions(-)

diff --git a/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py b/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py
index 12236f5e..233678ff 100644
--- a/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py
+++ b/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py
@@ -21,14 +21,11 @@ def get_stellib_boundaries(s, dlogT=0.1, dlogg=0.3, closed=True):
         if set, close the polygon
 
     Returns
-    ------
+    -------
     b : ndarray[float, ndim=2]
         (closed) boundary points: [logg, Teff]
+
     """
-    #        use "points_inside_poly" to test wether a point is inside the limits
-    #        >>> data = np.array([iso.data['logg'], iso.data['logT']]).T
-    #        >>> aa = points_inside_poly(data, leftb)
-    #    """
     leftb = [(k, np.max(s.logT[s.logg == k]) + dlogT) for k in np.unique(s.logg)]
     leftb += [(leftb[-1][0] + dlogg, leftb[-1][1])]
     leftb = [(leftb[0][0] - dlogg, leftb[0][1])] + leftb
diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
index 03d14365..bd39e607 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
@@ -17,6 +17,7 @@ def decode_params(filename):
     ----------
     filename : str
         name of file to decode
+
     """
     model_params = {}
 
diff --git a/beast/physicsmodel/stars/stellar_libraries/helpers.py b/beast/physicsmodel/stars/stellar_libraries/helpers.py
index 1e5dbf97..e2a36d91 100644
--- a/beast/physicsmodel/stars/stellar_libraries/helpers.py
+++ b/beast/physicsmodel/stars/stellar_libraries/helpers.py
@@ -25,6 +25,7 @@ def rebin_spectrum(wave, flux, resolution, wave_range):
     ------
     wave, flux, npts : tuple of vectors
         the model wavelength, flux, and npts at the requested wavelength
+
     """
     npts = int(
         np.log10(wave_range[1] / wave_range[0])

From a173edcf1bafeefa4352c016360629619f90be97 Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Fri, 1 May 2026 12:27:57 -0400
Subject: [PATCH 09/32] adding in bosz stellar library code

---
 beast/config.py                               |   2 +
 .../stars/stellar_libraries/check_tlusty25.py |  41 +------
 .../stars/stellar_libraries/gen_bosz24.py     | 113 ++++++++++++++++++
 .../stars/stellar_libraries/gen_tlusty25.py   |   2 +-
 .../stars/stellar_libraries/helpers.py        |  37 ++++++
 beast/physicsmodel/stars/stellib.py           |  96 ++++++++++++++-
 6 files changed, 251 insertions(+), 40 deletions(-)
 create mode 100644 beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py

diff --git a/beast/config.py b/beast/config.py
index 2223818a..e3c9073a 100644
--- a/beast/config.py
+++ b/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",
diff --git a/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py b/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py
index 233678ff..33f195bc 100644
--- a/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py
+++ b/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py
@@ -2,50 +2,17 @@
 import matplotlib.pyplot as plt
 
 from beast.physicsmodel.stars import stellib
-
-
-def get_stellib_boundaries(s, dlogT=0.1, dlogg=0.3, closed=True):
-    """
-    Returns the closed boundary polygon around the stellar library with
-    given margins.
-
-    Parameters
-    ----------
-    s :  Stellib
-        Stellar library object
-    dlogT :  float
-        margin in logT
-    dlogg  : float
-        margin in logg
-    closed : bool
-        if set, close the polygon
-
-    Returns
-    -------
-    b : ndarray[float, ndim=2]
-        (closed) boundary points: [logg, Teff]
-
-    """
-    leftb = [(k, np.max(s.logT[s.logg == k]) + dlogT) for k in np.unique(s.logg)]
-    leftb += [(leftb[-1][0] + dlogg, leftb[-1][1])]
-    leftb = [(leftb[0][0] - dlogg, leftb[0][1])] + leftb
-    rightb = [(k, np.min(s.logT[s.logg == k]) - dlogT) for k in np.unique(s.logg)[::-1]]
-    rightb += [(rightb[-1][0] - dlogg, rightb[-1][1])]
-    rightb = [(rightb[0][0] + dlogg, rightb[0][1])] + rightb
-    b = leftb + rightb
-    if closed:
-        b += [b[0]]
-    return np.array(b)
-
+from helpers import get_stellib_boundaries
 
 if __name__ == "__main__":  # pragma: no cover
 
     slib = stellib.Tlusty2025()
     bound = get_stellib_boundaries(slib, dlogT=0.01, dlogg=0.01)
     bbox = np.array(slib.bbox())
-    print(bbox[:, 0])
+    print(bbox)
+    print(bound)
 
     plt.plot(slib.logT, slib.logg, "ko")
-    plt.plot(bound[:, 1], bound[:, 0])
+    plt.plot(bound[:, 0], bound[:, 1])
     plt.plot(bbox[:, 0], bbox[:, 1])
     plt.show()
diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
new file mode 100644
index 00000000..4ef46a5a
--- /dev/null
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
@@ -0,0 +1,113 @@
+import glob
+
+import numpy as np
+
+from astropy.table import QTable
+from astropy.io import fits
+from astropy.io import ascii
+
+from helpers import rebin_spectrum
+
+
+def decode_params(filename):
+    """
+    Decode the tlusty filenames for the model parameters.
+
+    Parameters
+    ----------
+    filename : str
+        name of file to decode
+
+    """
+    model_params = {}
+
+    slashpos = filename.rfind("/")
+
+    zpos = filename.find("_m", slashpos)
+    tpos = filename.find("_t", slashpos)
+    gpos = filename.find("_g", slashpos)
+    vpos = filename.find("_v", slashpos)
+
+    # print(filename[zpos + 2 : zpos + 7])
+    # print(filename[tpos + 2 : gpos])
+    # print(filename[gpos + 2 : gpos + 6])
+    # print(filename[vpos + 2 : vpos + 3])
+    # exit()
+
+    model_params["Z"] = float(filename[zpos + 2 : zpos + 7])
+    model_params["Teff"] = float(filename[tpos + 2 : gpos])
+    model_params["logg"] = float(filename[gpos + 2 : gpos + 6])
+    model_params["vturb"] = float(filename[vpos + 2 : vpos + 3])
+
+    return model_params
+
+
+if __name__ == "__main__":  # pragma: no cover
+
+    solar_z = 0.02
+
+    # get all vtrub=2 models
+    path = "/home/kgordon/Python/extstar_data/Models/BOSZ2024/r5000/r5000/"
+    files = glob.glob(f"{path}/*/*a+0.00_c+0.00_v2_r5000_resam.txt.gz")
+    files = files[0:3]
+    n_files = len(files)
+
+    # fmt: off
+    outtab = QTable(
+        names=["Z", "Teff", "logg", "vturb", "logT", "logz"],
+        dtype=["f", "f", "f", "f", "f", "f"]
+    )
+    # fmt: on
+
+    for k, cfile in enumerate(files):
+        print(cfile)
+
+        # get model parameters
+        model_params = decode_params(cfile)
+        Z = (10 ** model_params["Z"]) * solar_z
+        row = (Z, model_params["Teff"], model_params["logg"], model_params["vturb"], np.log10(model_params["Teff"]),
+               np.log10(Z))
+        outtab.add_row(row)
+
+        # read in the model spectrum
+        # wave units are angstrom
+        # flux units are H, F = 4*pi*H in ergs/(cm^2 s A)
+        # wavelength first
+        wave_filename = f"{path}/bosz2024_wave_r5000.txt"
+        mspec_lte_wave = ascii.read(
+            wave_filename,
+            format="no_header",
+            names=["Wave"],
+        )
+
+        mspec = ascii.read(
+            cfile,
+            format="no_header",
+            names=["SFlux", "SCont"],
+        )
+        mspec["Wave"] = mspec_lte_wave["Wave"]
+
+        # convert the type to float
+        mspec["SFlux"] = mspec["SFlux"].astype(float)
+
+        # now extract the wave and flux columns
+        mwave = mspec["Wave"]
+        mflux = mspec["SFlux"]
+
+        # rebin to R=4000 for speed, common res and wave range with BOSZ LTE models
+        rbres = 4000.0
+        wave_rebin, flux_rebin, npts_rebin = rebin_spectrum(
+            mwave.value, mflux.value, rbres, [500.0, 320000.0]
+        )
+
+        if k == 0:
+            outspec = np.zeros((len(files) + 1, len(wave_rebin)))
+            outspec[-1, :] = wave_rebin
+        outspec[k, :] = flux_rebin * 4.0 * np.pi
+
+    # output the stellar library in the beast format
+    spec_hdu = fits.PrimaryHDU(data=outspec)
+    table_hdu = fits.BinTableHDU(data=outtab)
+    table_hdu.name = "BOSZ"
+    hdul = fits.HDUList([spec_hdu, table_hdu])
+    hdul.writeto("bosz2024.grid.fits", overwrite=True)
diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
index bd39e607..6b9ab3b9 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
@@ -72,7 +72,7 @@ def decode_params(filename):
         # convert the type to float
         mspec["SFlux"] = mspec["SFlux"].astype(float)
 
-        # now extract the wave and flux colums
+        # now extract the wave and flux columns
         mwave = mspec["Wave"]
         mflux = mspec["SFlux"]
 
diff --git a/beast/physicsmodel/stars/stellar_libraries/helpers.py b/beast/physicsmodel/stars/stellar_libraries/helpers.py
index e2a36d91..51a8c390 100644
--- a/beast/physicsmodel/stars/stellar_libraries/helpers.py
+++ b/beast/physicsmodel/stars/stellar_libraries/helpers.py
@@ -70,3 +70,40 @@ def rebin_spectrum(wave, flux, resolution, wave_range):
             full_npts[nanindxs] = 0
 
     return (full_wave, full_flux, full_npts)
+
+
+def get_stellib_boundaries(s, dlogT=0.1, dlogg=0.3, closed=True):
+    """
+    Returns the closed boundary polygon around the stellar library with
+    given margins.
+
+    Parameters
+    ----------
+    s :  Stellib
+        Stellar library object
+    dlogT :  float
+        margin in logT
+    dlogg  : float
+        margin in logg
+    closed : bool
+        if set, close the polygon
+
+    Returns
+    -------
+    b : ndarray[float, ndim=2]
+        (closed) boundary points: [logg, Teff]
+
+    """
+    leftb = [(k, np.max(s.logT[s.logg == k]) + dlogT) for k in np.unique(s.logg)]
+    leftb += [(leftb[-1][0] + dlogg, leftb[-1][1])]
+    leftb = [(leftb[0][0] - dlogg, leftb[0][1])] + leftb
+    rightb = [(k, np.min(s.logT[s.logg == k]) - dlogT) for k in np.unique(s.logg)[::-1]]
+    rightb += [(rightb[-1][0] - dlogg, rightb[-1][1])]
+    rightb = [(rightb[0][0] + dlogg, rightb[0][1])] + rightb
+    b = leftb + rightb
+    if closed:
+        b += [b[0]]
+    b = np.array(b)
+    outb = [b[:, 1], b[:, 0]]
+
+    return np.transpose(np.array(outb))
\ No newline at end of file
diff --git a/beast/physicsmodel/stars/stellib.py b/beast/physicsmodel/stars/stellib.py
index ea6a974f..b3e2fa40 100644
--- a/beast/physicsmodel/stars/stellib.py
+++ b/beast/physicsmodel/stars/stellib.py
@@ -16,7 +16,6 @@
 from astropy.table import Table
 
 from beast.physicsmodel.grid import SpectralGrid
-# from beast.external.eztables import Table
 from beast.config import __ROOT__, __NTHREADS__
 from beast.physicsmodel.stars.include import __interp__
 from beast.tools.helpers import nbytes
@@ -34,6 +33,7 @@
     "btsettl_medres": __ROOT__ + "bt-settl.medres.grid.fits",
     "munari": __ROOT__ + "atlas9-munari.hires.grid.fits",
     "aringer": __ROOT__ + "Aringer.AGB.grid.fits",
+    "bosz2024": __ROOT__ + "BOSZ2024.grid.fits",
     "tlusty2025": __ROOT__ + "tlusty2025.grid.fits",
 }
 
@@ -47,7 +47,8 @@
     "Elodie",
     "BaSeL",
     "Aringer",
-    "Tlusty2025"
+    "BOSZ2024",
+    "Tlusty2025",
 ]
 
 
@@ -1947,6 +1948,97 @@ def logZ(self):
         return self.grid["logz"]
 
 
+
+class BOSZ2024(Stellib):
+    """
+    BOSZ 2024 LTE stellar atmospheres for A stars and cooler.
+    Updated in 2024 with better line lists and extended IR wavelength coverage
+
+    * LTE
+    * Plane Parallel
+
+    References
+    ----------
+    Bohlin, Meszaros, et al. 2017, AJ, 153, 234
+    Meszaros, Bohlin et al. 2024, A&A, 688, 197
+
+    files are available at: https://archive.stsci.edu/hlsp/bosz
+    """
+
+    def __init__(self, filename=None, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        self.name = "Tlusty2025"
+        if filename is None:
+            self.source = config["tlusty2025"]
+        else:
+            self.source = filename
+        self._load_()
+
+    def _load_(self):
+        g = SpectralGrid(self.source, backend="memory")
+        self.wavelength = g.lamb
+        self.grid = g.grid
+        self.header["NAME"] = "tlusty2025"
+        self.spectra = g.seds
+
+    def bbox(self, dlogT=0.05, dlogg=0.25):
+        """ Boundary of Tlusty library
+
+        Parameters
+        ----------
+        dlogT: float
+            log-temperature tolerance before extrapolation limit
+
+        dlogg: float
+            log-g tolerance before extrapolation limit
+
+        Returns
+        -------
+        bbox: ndarray
+            (logT, logg) edges of the bounding polygon
+        """
+        bbox = [
+            (3.54406 - dlogT, 5.000 + dlogg),
+            (3.55403 - dlogT, 0.000 - dlogg),
+            (3.778, 0.000 - dlogg),
+            (3.778 + dlogT, 0.000),
+            (3.875 + dlogT, 0.500),
+            (3.929 + dlogT, 1.000),
+            (3.954 + dlogT, 1.500),
+            (4.146, 2.000 - dlogg),
+            (4.146 + dlogT, 2.000),
+            (4.279 + dlogT, 2.500),
+            (4.415 + dlogT, 3.000),
+            (4.491 + dlogT, 3.500),
+            (4.591 + dlogT, 4.000),
+            (4.689 + dlogT, 4.500),
+            (4.699 + dlogT, 5.000 + dlogg),
+            (3.544 - dlogT, 5.000 + dlogg),
+        ]
+
+        return np.array(bbox)
+
+    @property
+    def logg(self):
+        return self.grid["logg"]
+
+    @property
+    def Teff(self):
+        return self.grid["Teff"]
+
+    @property
+    def logT(self):
+        return np.log10(self.grid["Teff"])
+
+    @property
+    def Z(self):
+        return self.grid["Z"]
+
+    @property
+    def logZ(self):
+        return np.log10(self.grid["Z"])
+
+
 class Tlusty2025(Stellib):
     """
     Tlusty 2025 O and B stellar atmospheres.  Updated in 2025 with better line lists

From 2e5e2cc78714d182f7daa037b80a4e41287334da Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Fri, 1 May 2026 12:39:22 -0400
Subject: [PATCH 10/32] cleanup

---
 beast/config.py                               |  2 +-
 .../stars/stellar_libraries/helpers.py        |  2 +-
 beast/physicsmodel/stars/stellib.py           | 92 ++++++++++---------
 3 files changed, 49 insertions(+), 47 deletions(-)

diff --git a/beast/config.py b/beast/config.py
index e3c9073a..a0c68343 100644
--- a/beast/config.py
+++ b/beast/config.py
@@ -34,7 +34,7 @@
     vega="vega.hd5",
     filters="filters.hd5",
     kurucz04="kurucz2004.grid.fits",
-    bosz24="BOSZ2024.grid.fits",
+    bosz24="bosz2024.grid.fits",
     tlusty09="tlusty.lowres.grid.fits",
     tlusty25="tlusty2025.grid.fits",
     hstcovar="hst_whitedwarf_frac_covar.fits",
diff --git a/beast/physicsmodel/stars/stellar_libraries/helpers.py b/beast/physicsmodel/stars/stellar_libraries/helpers.py
index 51a8c390..1b8d46b5 100644
--- a/beast/physicsmodel/stars/stellar_libraries/helpers.py
+++ b/beast/physicsmodel/stars/stellar_libraries/helpers.py
@@ -106,4 +106,4 @@ def get_stellib_boundaries(s, dlogT=0.1, dlogg=0.3, closed=True):
     b = np.array(b)
     outb = [b[:, 1], b[:, 0]]
 
-    return np.transpose(np.array(outb))
\ No newline at end of file
+    return np.transpose(np.array(outb))
diff --git a/beast/physicsmodel/stars/stellib.py b/beast/physicsmodel/stars/stellib.py
index b3e2fa40..05f38dc2 100644
--- a/beast/physicsmodel/stars/stellib.py
+++ b/beast/physicsmodel/stars/stellib.py
@@ -7,6 +7,7 @@
 The interpolation is implemented from the pegase.2 fortran converted algorithm.
 (this may not be pythonic though)
 """
+
 import numpy as np
 from scipy.interpolate import interp1d
 from numpy.lib import recfunctions
@@ -53,21 +54,21 @@
 
 
 def isNestedInstance(obj, cl):
-    """ Test for sub-classes types
-        I could not find a universal test
+    """Test for sub-classes types
+    I could not find a universal test
 
-        keywords
-        --------
-        obj: object instance
-            object to test
+    keywords
+    --------
+    obj: object instance
+        object to test
 
-        cl: Class
-            top level class to test
+    cl: Class
+        top level class to test
 
-        returns
-        -------
-        r: bool
-            True if obj is indeed an instance or subclass instance of cl
+    returns
+    -------
+    r: bool
+        True if obj is indeed an instance or subclass instance of cl
     """
     tree = []
     for k in cl.__subclasses__():
@@ -105,7 +106,7 @@ def __interpMany__(
     pool=None,
     nthreads=__NTHREADS__,
 ):
-    """ run interp on a list of inputs and returns reduced results
+    """run interp on a list of inputs and returns reduced results
 
     Interpolation of the T,g grid at Z0 metallicity
 
@@ -198,7 +199,7 @@ def __interpMany__(
 
 
 def interp(T0, g0, Z0, L0, T, g, Z, dT_max=0.1, eps=1e-6, weight=1.0):
-    """ Interpolation of the T,g grid
+    """Interpolation of the T,g grid
 
     Interpolate on the grid and returns star indices and
     associated weights, and Z.
@@ -285,14 +286,14 @@ def interp(T0, g0, Z0, L0, T, g, Z, dT_max=0.1, eps=1e-6, weight=1.0):
                 weights[:8] *= 0.5
 
     ind = np.where(weights > 0)
-    return index[ind].astype(int), 10 ** L0 * weight * weights[ind]
+    return index[ind].astype(int), 10**L0 * weight * weights[ind]
 
 
 class Stellib(object):
-    """ Basic stellar library class """
+    """Basic stellar library class"""
 
     def __init__(self, *args, **kargs):
-        """ Contructor """
+        """Contructor"""
         self.header = {}
 
     def _load_(self):
@@ -300,11 +301,11 @@ def _load_(self):
 
     @property
     def nbytes(self):
-        """ return the number of bytes of the object """
+        """return the number of bytes of the object"""
         return nbytes(self)
 
     def interp(self, T0, g0, Z0, L0, dT_max=0.1, eps=1e-6):
-        """ Interpolation of the T,g grid
+        """Interpolation of the T,g grid
 
         Interpolate on the grid and returns star indices and
         associated weights, and Z.
@@ -360,7 +361,7 @@ def interpMany(
         pool=None,
         nthreads=__NTHREADS__,
     ):
-        """ run interp on a list of inputs and returns reduced results
+        """run interp on a list of inputs and returns reduced results
 
         Interpolation of the T,g grid at Z0 metallicity
 
@@ -457,7 +458,7 @@ def points_inside(self, xypoints, dlogT=0.1, dlogg=0.3):
         return p.contains_points(xypoints)
 
     def get_radius(self, logl, logt):
-        r""" Returns the radius of a star given its luminosity and temperature
+        r"""Returns the radius of a star given its luminosity and temperature
 
         Assuming a black body, it comes:
                 R ^ 2 = L / ( 4 \pi \sigma T ^ 4 ),
@@ -488,7 +489,7 @@ def get_radius(self, logl, logt):
         )
 
     def get_boundaries(self, dlogT=0.1, dlogg=0.3, **kwargs):
-        """ Returns the closed boundary polygon around the stellar library with
+        """Returns the closed boundary polygon around the stellar library with
         given margins
 
         Parameters
@@ -543,7 +544,7 @@ def get_boundaries(self, dlogT=0.1, dlogg=0.3, **kwargs):
         return self._bound[0]
 
     def genQ(self, qname, r, **kwargs):
-        """ Generate a composite value from a previously calculated
+        """Generate a composite value from a previously calculated
             interpolation
             Works on 1 desired star or a population of stars
 
@@ -563,7 +564,7 @@ def genQ(self, qname, r, **kwargs):
         return (self.grid[qname][r[:, 0].astype(int)] * r[:, 1]).sum()
 
     def genSpectrum(self, T0, g0=None, Z0=None, weights=None, **kwargs):
-        """ Generate a composite sprectrum
+        """Generate a composite sprectrum
         Does the interpolation or uses a previously calculated
         interpolation
         Works on 1 desired star or a population of stars
@@ -797,7 +798,7 @@ def __repr__(self):
 
 
 class CompositeStellib(Stellib):
-    """ Generates an object from the union of multiple individual libraries """
+    """Generates an object from the union of multiple individual libraries"""
 
     def __init__(self, osllist, *args, **kwargs):
         super().__init__(*args, **kwargs)
@@ -819,8 +820,8 @@ def __radd__(self, other):
 
     @property
     def wavelength(self):
-        """ return a common wavelength sampling to all libraries. This can be
-        used to reinterpolate any spectrum onto a common definition """
+        """return a common wavelength sampling to all libraries. This can be
+        used to reinterpolate any spectrum onto a common definition"""
 
         lambs = np.unique(np.concatenate([osl.wavelength[:] for osl in self._olist]))
         return lambs
@@ -919,7 +920,7 @@ def __repr__(self):
         )
 
     def get_boundaries(self, dlogT=0.1, dlogg=0.3, **kwargs):
-        """ Returns the closed boundary polygon around the stellar library with
+        """Returns the closed boundary polygon around the stellar library with
         given margins
 
         Parameters
@@ -957,7 +958,7 @@ def get_boundaries(self, dlogT=0.1, dlogg=0.3, **kwargs):
         return self._bound[0]
 
     def interp(self, T0, g0, Z0, L0, dT_max=0.1, eps=1e-6, bounds={}):
-        """ Interpolation of the T,g grid
+        """Interpolation of the T,g grid
 
         Interpolate on the grid and returns star indices and
         associated weights, and Z.
@@ -1026,7 +1027,7 @@ def interpMany(
         pool=None,
         nthreads=__NTHREADS__,
     ):
-        """ run interp on a list of inputs and returns reduced results
+        """run interp on a list of inputs and returns reduced results
 
         Interpolation of the T,g grid at Z0 metallicity
 
@@ -1108,7 +1109,7 @@ def interpMany(
         return g
 
     def genQ(self, qname, r, **kwargs):
-        """ Generate a composite value from a previously calculated
+        """Generate a composite value from a previously calculated
             interpolation
             Works on 1 desired star or a population of stars
 
@@ -1133,7 +1134,7 @@ def genQ(self, qname, r, **kwargs):
         return vals
 
     def genSpectrum(self, T0, g0=None, Z0=None, weights=None, **kwargs):
-        """ Generate a composite sprectrum
+        """Generate a composite sprectrum
             Does the interpolation or uses a previously calculated
             interpolation
             Works on 1 desired star or a population of stars
@@ -1280,7 +1281,7 @@ def gen_spectral_grid_from_given_points(
 
 
 class Elodie(Stellib):
-    """ Elodie 3.1 stellar library derived class
+    """Elodie 3.1 stellar library derived class
 
     This library matches BaSeL 2.2 grid definition
 
@@ -1303,7 +1304,7 @@ def _load_(self):
         self.spectra = g.seds
 
     def bbox(self, dlogT=0.05, dlogg=0.25):
-        """ Boundary of Elodie library
+        """Boundary of Elodie library
 
         Parameters
         ----------
@@ -1385,7 +1386,7 @@ def NHeII(self):
 
 
 class BaSeL(Stellib):
-    """ BaSeL 2.2 (This library is used in Pegase.2)
+    """BaSeL 2.2 (This library is used in Pegase.2)
         This library is used in Pegase.2
 
     The BaSeL stellar spectral energy distribution (SED) libraries are
@@ -1414,7 +1415,7 @@ def _load_(self):
         self.spectra = g.seds
 
     def bbox(self, dlogT=0.05, dlogg=0.25):
-        """ Boundary of Basel 2.2 library
+        """Boundary of Basel 2.2 library
 
         Parameters
         ----------
@@ -1517,7 +1518,9 @@ def __init__(self, filename=None, *args, **kwargs):
             self.source = filename
         self._load_()
 
-    def _load_(self,):
+    def _load_(
+        self,
+    ):
         g = SpectralGrid(self.source, backend="memory")
         self.wavelength = g.lamb
         self.grid = g.grid
@@ -1525,7 +1528,7 @@ def _load_(self,):
         self.spectra = g.seds
 
     def bbox(self, dlogT=0.05, dlogg=0.25):
-        """ Boundary of Kurucz 2004 library
+        """Boundary of Kurucz 2004 library
 
         Parameters
         ----------
@@ -1623,7 +1626,7 @@ def _load_(self):
         self.spectra = g.seds
 
     def bbox(self, dlogT=0.05, dlogg=0.25):
-        """ Boundary of Tlusty library
+        """Boundary of Tlusty library
 
         Parameters
         ----------
@@ -1711,7 +1714,7 @@ def _load_(self):
         self.spectra = g.seds
 
     def bbox(self, dlogT=0.05, dlogg=0.25):
-        """ Boundary of BT-Settl library
+        """Boundary of BT-Settl library
 
         Parameters
         ----------
@@ -1800,7 +1803,7 @@ def _load_(self):
         self.spectra = g.seds
 
     def bbox(self, dlogT=0.05, dlogg=0.25):
-        """ Boundary of Munari library
+        """Boundary of Munari library
 
         Parameters
         ----------
@@ -1886,7 +1889,7 @@ def _load_(self):
         self.spectra = g.seds
 
     def bbox(self, dlogT=0.05, dlogg=0.25):
-        """ Boundary of Aringer library for C+M+K giants
+        """Boundary of Aringer library for C+M+K giants
 
         Parameters
         ----------
@@ -1948,7 +1951,6 @@ def logZ(self):
         return self.grid["logz"]
 
 
-
 class BOSZ2024(Stellib):
     """
     BOSZ 2024 LTE stellar atmospheres for A stars and cooler.
@@ -1982,7 +1984,7 @@ def _load_(self):
         self.spectra = g.seds
 
     def bbox(self, dlogT=0.05, dlogg=0.25):
-        """ Boundary of Tlusty library
+        """Boundary of Tlusty library
 
         Parameters
         ----------
@@ -2074,7 +2076,7 @@ def _load_(self):
         self.spectra = g.seds
 
     def bbox(self, dlogT=0.05, dlogg=0.25):
-        """ Boundary of Tlusty library
+        """Boundary of Tlusty library
 
         Parameters
         ----------

From 49e42078593fab6fc333e9442944bdd99feba144 Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Fri, 1 May 2026 13:54:10 -0400
Subject: [PATCH 11/32] updating to handle ATLAS vs MARCS

---
 .../stars/stellar_libraries/gen_bosz24.py     | 106 ++++++++++--------
 1 file changed, 59 insertions(+), 47 deletions(-)

diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
index 4ef46a5a..d4127973 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
@@ -23,9 +23,9 @@ def decode_params(filename):
 
     slashpos = filename.rfind("/")
 
-    zpos = filename.find("_m", slashpos)
-    tpos = filename.find("_t", slashpos)
     gpos = filename.find("_g", slashpos)
+    zpos = filename.find("_m", gpos)
+    tpos = filename.find("_t", slashpos)
     vpos = filename.find("_v", slashpos)
 
     # print(filename[zpos + 2 : zpos + 7])
@@ -64,50 +64,62 @@ def decode_params(filename):
 
         # get model parameters
         model_params = decode_params(cfile)
-        Z = (10 ** model_params["Z"]) * solar_z
-        row = (Z, model_params["Teff"], model_params["logg"], model_params["vturb"], np.log10(model_params["Teff"]),
-               np.log10(Z))
-        outtab.add_row(row)
-
-        # read in the model spectrum
-        # wave units are angstrom
-        # flux units are H, F = 4*pi*H in ergs/(cm^2 s A)
-        # wavelength first
-        wave_filename = f"{path}/bosz2024_wave_r5000.txt"
-        mspec_lte_wave = ascii.read(
-            wave_filename,
-            format="no_header",
-            names=["Wave"],
-        )
-
-        mspec = ascii.read(
-            cfile,
-            format="no_header",
-            names=["SFlux", "SCont"],
-        )
-        mspec["Wave"] = mspec_lte_wave["Wave"]
-
-        # convert the type to float
-        mspec["SFlux"] = mspec["SFlux"].astype(float)
-
-        # now extract the wave and flux columns
-        mwave = mspec["Wave"]
-        mflux = mspec["SFlux"]
-
-        # rebin to R=4000 for speed, common res and wave range with BOSZ LTE models
-        rbres = 4000.0
-        wave_rebin, flux_rebin, npts_rebin = rebin_spectrum(
-            mwave.value, mflux.value, rbres, [500.0, 320000.0]
-        )
-
-        if k == 0:
-            outspec = np.zeros((len(files) + 1, len(wave_rebin)))
-            outspec[-1, :] = wave_rebin
-        outspec[k, :] = flux_rebin * 4.0 * np.pi
+
+        # skip if ATLAS ("ap") and 8000 K or lower
+        # MARCS models cover these temps
+        usemod = True
+        if ("ap" in cfile) and (model_params["Teff"] <= 8000.0):
+            usemod = False
+            print("skipped")
+            exit()
+
+        usemod = False
+        if usemod:
+            Z = (10 ** model_params["Z"]) * solar_z
+            row = (Z, model_params["Teff"], model_params["logg"], model_params["vturb"], np.log10(model_params["Teff"]),
+                np.log10(Z))
+            outtab.add_row(row)
+            print(row)
+
+            # read in the model spectrum
+            # wave units are angstrom
+            # flux units are H, F = 4*pi*H in ergs/(cm^2 s A)
+            # wavelength first
+            wave_filename = f"{path}/bosz2024_wave_r5000.txt"
+            mspec_lte_wave = ascii.read(
+                wave_filename,
+                format="no_header",
+                names=["Wave"],
+            )
+
+            mspec = ascii.read(
+                cfile,
+                format="no_header",
+                names=["SFlux", "SCont"],
+            )
+            mspec["Wave"] = mspec_lte_wave["Wave"]
+
+            # convert the type to float
+            mspec["SFlux"] = mspec["SFlux"].astype(float)
+
+            # now extract the wave and flux columns
+            mwave = mspec["Wave"]
+            mflux = mspec["SFlux"]
+
+            # rebin to R=4000 for speed, common res and wave range with BOSZ LTE models
+            rbres = 4000.0
+            wave_rebin, flux_rebin, npts_rebin = rebin_spectrum(
+                mwave.value, mflux.value, rbres, [500.0, 320000.0]
+            )
+
+            if k == 0:
+                outspec = np.zeros((len(files) + 1, len(wave_rebin)))
+                outspec[-1, :] = wave_rebin
+            outspec[k, :] = flux_rebin * 4.0 * np.pi
 
     # output the stellar library in the beast format
-    spec_hdu = fits.PrimaryHDU(data=outspec)
-    table_hdu = fits.BinTableHDU(data=outtab)
-    table_hdu.name = "BOSZ"
-    hdul = fits.HDUList([spec_hdu, table_hdu])
-    hdul.writeto("bosz2024.grid.fits", overwrite=True)
+    # spec_hdu = fits.PrimaryHDU(data=outspec)
+    # table_hdu = fits.BinTableHDU(data=outtab)
+    # table_hdu.name = "BOSZ"
+    # hdul = fits.HDUList([spec_hdu, table_hdu])
+    # hdul.writeto("bosz2024.grid.fits", overwrite=True)

From 17fd50ca7fc330ea68111f20f2b5b5f8366cdd9f Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Fri, 1 May 2026 13:56:10 -0400
Subject: [PATCH 12/32] all files!

---
 beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
index d4127973..412fbd23 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
@@ -49,7 +49,6 @@ def decode_params(filename):
     # get all vtrub=2 models
     path = "/home/kgordon/Python/extstar_data/Models/BOSZ2024/r5000/r5000/"
     files = glob.glob(f"{path}/*/*a+0.00_c+0.00_v2_r5000_resam.txt.gz")
-    files = files[0:3]
     n_files = len(files)
 
     # fmt: off

From f2e8e75471de46526c4db8719426cf57aa314ae1 Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Fri, 1 May 2026 14:07:37 -0400
Subject: [PATCH 13/32] going for the full grid

---
 beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
index 412fbd23..2bf54cd3 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
@@ -70,9 +70,7 @@ def decode_params(filename):
         if ("ap" in cfile) and (model_params["Teff"] <= 8000.0):
             usemod = False
             print("skipped")
-            exit()
 
-        usemod = False
         if usemod:
             Z = (10 ** model_params["Z"]) * solar_z
             row = (Z, model_params["Teff"], model_params["logg"], model_params["vturb"], np.log10(model_params["Teff"]),

From 6e46f0338afd35e97d6c528adc05c5ca4c4ff871 Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Fri, 1 May 2026 14:42:41 -0400
Subject: [PATCH 14/32] updating to get the right assumed solar abundances

---
 .../stars/stellar_libraries/gen_bosz24.py           | 13 ++++++++++---
 .../stars/stellar_libraries/gen_tlusty25.py         |  2 +-
 2 files changed, 11 insertions(+), 4 deletions(-)

diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
index 2bf54cd3..fee7da88 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
@@ -44,7 +44,8 @@ def decode_params(filename):
 
 if __name__ == "__main__":  # pragma: no cover
 
-    solar_z = 0.02
+    # from Asplund et al. (2005)/Grevasse et al. (2007) as referenced in BOSZ2024 paper
+    solar_z = 0.012
 
     # get all vtrub=2 models
     path = "/home/kgordon/Python/extstar_data/Models/BOSZ2024/r5000/r5000/"
@@ -73,8 +74,14 @@ def decode_params(filename):
 
         if usemod:
             Z = (10 ** model_params["Z"]) * solar_z
-            row = (Z, model_params["Teff"], model_params["logg"], model_params["vturb"], np.log10(model_params["Teff"]),
-                np.log10(Z))
+            row = (
+                Z,
+                model_params["Teff"],
+                model_params["logg"],
+                model_params["vturb"],
+                np.log10(model_params["Teff"]),
+                np.log10(Z),
+            )
             outtab.add_row(row)
             print(row)
 
diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
index 6b9ab3b9..6bc10bd2 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
@@ -42,7 +42,7 @@ def decode_params(filename):
 
 if __name__ == "__main__":  # pragma: no cover
 
-    solar_z = 0.02
+    solar_z = 0.017  # from Grevasse & Sauval (1998) - TBC
 
     # get all vtrub=2 models
     path = "/home/kgordon/Python/extstar_data/Models/Tlusty_2025/"

From 9aebf668f4905915302bcb49c8b38e27ade4e366 Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Sat, 2 May 2026 11:25:22 -0400
Subject: [PATCH 15/32] gotta actually output the resuls - lol

---
 .../physicsmodel/stars/stellar_libraries/gen_bosz24.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
index fee7da88..6ed8f705 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
@@ -122,8 +122,8 @@ def decode_params(filename):
             outspec[k, :] = flux_rebin * 4.0 * np.pi
 
     # output the stellar library in the beast format
-    # spec_hdu = fits.PrimaryHDU(data=outspec)
-    # table_hdu = fits.BinTableHDU(data=outtab)
-    # table_hdu.name = "BOSZ"
-    # hdul = fits.HDUList([spec_hdu, table_hdu])
-    # hdul.writeto("bosz2024.grid.fits", overwrite=True)
+    spec_hdu = fits.PrimaryHDU(data=outspec)
+    table_hdu = fits.BinTableHDU(data=outtab)
+    table_hdu.name = "BOSZ"
+    hdul = fits.HDUList([spec_hdu, table_hdu])
+    hdul.writeto("bosz2024.grid.fits", overwrite=True)

From 1d0659c46c058049b738586ea345226eb2feaa1a Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Sun, 24 May 2026 07:25:30 -0400
Subject: [PATCH 16/32] continuing debugging

---
 .../stars/stellar_libraries/check_tlusty25.py | 12 ++--
 .../stars/stellar_libraries/gen_bosz24.py     |  4 +-
 .../stars/stellar_libraries/gen_tlusty25.py   | 60 ++++++++--------
 beast/physicsmodel/stars/stellib.py           | 70 +++++++++----------
 4 files changed, 75 insertions(+), 71 deletions(-)

diff --git a/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py b/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py
index 33f195bc..4bd6363c 100644
--- a/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py
+++ b/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py
@@ -7,12 +7,14 @@
 if __name__ == "__main__":  # pragma: no cover
 
     slib = stellib.Tlusty2025()
-    bound = get_stellib_boundaries(slib, dlogT=0.01, dlogg=0.01)
-    bbox = np.array(slib.bbox())
-    print(bbox)
+    bound = get_stellib_boundaries(slib, dlogT=0.0, dlogg=0.0)
+    #bbox = slib.get_boundaries(dlogT=0.01, dlogg=0.01)
+    #print(bbox)
     print(bound)
 
     plt.plot(slib.logT, slib.logg, "ko")
-    plt.plot(bound[:, 0], bound[:, 1])
-    plt.plot(bbox[:, 0], bbox[:, 1])
+    plt.plot(bound[:, 0], bound[:, 1], "b-", label="new")
+    slib.plot_boundary(dlogT=0.0, dlogg=0.0)
+    #plt.plot(bbox[:, 0], bbox[:, 1], "g-", label="in lib")
+    plt.legend()
     plt.show()
diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
index 6ed8f705..aaae2161 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
@@ -80,7 +80,7 @@ def decode_params(filename):
                 model_params["logg"],
                 model_params["vturb"],
                 np.log10(model_params["Teff"]),
-                np.log10(Z),
+                model_params["Z"],
             )
             outtab.add_row(row)
             print(row)
@@ -111,7 +111,7 @@ def decode_params(filename):
             mflux = mspec["SFlux"]
 
             # rebin to R=4000 for speed, common res and wave range with BOSZ LTE models
-            rbres = 4000.0
+            rbres = 2000.0
             wave_rebin, flux_rebin, npts_rebin = rebin_spectrum(
                 mwave.value, mflux.value, rbres, [500.0, 320000.0]
             )
diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
index 6bc10bd2..a8154140 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
@@ -59,38 +59,40 @@ def decode_params(filename):
     for k, cfile in enumerate(files):
         print(cfile)
 
-        # read in the model spectrum
-        # wave units are angstrom
-        # flux units are H, F = 4*pi*H in ergs/(cm^2 s A)
-        mspec = ascii.read(
-            cfile,
-            format="no_header",
-            fast_reader={"exponent_style": "D"},
-            names=["Wave", "SFlux"],
-        )
-
-        # convert the type to float
-        mspec["SFlux"] = mspec["SFlux"].astype(float)
-
-        # now extract the wave and flux columns
-        mwave = mspec["Wave"]
-        mflux = mspec["SFlux"]
-
-        # rebin to R=4000 for speed, common res and wave range with BOSZ LTE models
-        rbres = 4000.0
-        wave_rebin, flux_rebin, npts_rebin = rebin_spectrum(
-            mwave.value, mflux.value, rbres, [500.0, 320000.0]
-        )
-
-        if k == 0:
-            outspec = np.zeros((len(files) + 1, len(wave_rebin)))
-            outspec[-1, :] = wave_rebin
-        outspec[k, :] = flux_rebin * 4.0 * np.pi
+        # # read in the model spectrum
+        # # wave units are angstrom
+        # # flux units are H, F = 4*pi*H in ergs/(cm^2 s A)
+        # mspec = ascii.read(
+        #     cfile,
+        #     format="no_header",
+        #     fast_reader={"exponent_style": "D"},
+        #     names=["Wave", "SFlux"],
+        # )
+
+        # # convert the type to float
+        # mspec["SFlux"] = mspec["SFlux"].astype(float)
+
+        # # now extract the wave and flux columns
+        # mwave = mspec["Wave"]
+        # mflux = mspec["SFlux"]
+
+        # # rebin to R=4000 for speed, common res and wave range with BOSZ LTE models
+        # rbres = 2000.0
+        # wave_rebin, flux_rebin, npts_rebin = rebin_spectrum(
+        #     mwave.value, mflux.value, rbres, [500.0, 320000.0]
+        # )
+
+        # if k == 0:
+        #     outspec = np.zeros((len(files) + 1, len(wave_rebin)))
+        #     outspec[-1, :] = wave_rebin
+        # outspec[k, :] = flux_rebin * 4.0 * np.pi
 
         model_params = decode_params(cfile)
         Z = model_params["Z"] * solar_z
-        row = (Z, model_params["Teff"], model_params["logg"], model_params["vturb"], np.log10(model_params["Teff"]),
-               np.log10(Z))
+        row = (Z, model_params["Teff"], model_params["logg"], model_params["vturb"],
+               np.log10(model_params["Teff"]),
+               model_params["Z"])
+        print(row)
         outtab.add_row(row)
 
     # output the stellar library in the beast format
diff --git a/beast/physicsmodel/stars/stellib.py b/beast/physicsmodel/stars/stellib.py
index 05f38dc2..fbda9004 100644
--- a/beast/physicsmodel/stars/stellib.py
+++ b/beast/physicsmodel/stars/stellib.py
@@ -34,7 +34,7 @@
     "btsettl_medres": __ROOT__ + "bt-settl.medres.grid.fits",
     "munari": __ROOT__ + "atlas9-munari.hires.grid.fits",
     "aringer": __ROOT__ + "Aringer.AGB.grid.fits",
-    "bosz2024": __ROOT__ + "BOSZ2024.grid.fits",
+    "bosz2024": __ROOT__ + "bosz2024.grid.fits",
     "tlusty2025": __ROOT__ + "tlusty2025.grid.fits",
 }
 
@@ -702,7 +702,15 @@ def gen_spectral_grid_from_given_points(
         ):
             if bound_cond[mk]:
                 s = np.array(self.interp(rT, rg, rZ, 0.0)).T
+                # print("logT, logg, Z")
+                # print(rT, rg, rZ)
+                # print(s)
+                # print("logT", self.grid["logT"][s[:, 0].astype(int)].data)
+                # print("logg", self.grid["logg"][s[:, 0].astype(int)].data)
+                # print("Z", self.grid["Z"][s[:, 0].astype(int)].data)
                 specs[mk, :] = self.genSpectrum(s) * weights[mk]
+                # print(specs[mk, :])
+                # exit()
 
         # Step 4: filter points without spectrum
         # ======================================
@@ -1969,9 +1977,9 @@ class BOSZ2024(Stellib):
 
     def __init__(self, filename=None, *args, **kwargs):
         super().__init__(*args, **kwargs)
-        self.name = "Tlusty2025"
+        self.name = "BOSZ2024"
         if filename is None:
-            self.source = config["tlusty2025"]
+            self.source = config["bosz2024"]
         else:
             self.source = filename
         self._load_()
@@ -1999,23 +2007,21 @@ def bbox(self, dlogT=0.05, dlogg=0.25):
         bbox: ndarray
             (logT, logg) edges of the bounding polygon
         """
+
         bbox = [
-            (3.54406 - dlogT, 5.000 + dlogg),
-            (3.55403 - dlogT, 0.000 - dlogg),
-            (3.778, 0.000 - dlogg),
-            (3.778 + dlogT, 0.000),
-            (3.875 + dlogT, 0.500),
-            (3.929 + dlogT, 1.000),
-            (3.954 + dlogT, 1.500),
-            (4.146, 2.000 - dlogg),
-            (4.146 + dlogT, 2.000),
-            (4.279 + dlogT, 2.500),
-            (4.415 + dlogT, 3.000),
-            (4.491 + dlogT, 3.500),
-            (4.591 + dlogT, 4.000),
-            (4.689 + dlogT, 4.500),
-            (4.699 + dlogT, 5.000 + dlogg),
-            (3.544 - dlogT, 5.000 + dlogg),
+            (3.447 - dlogT, 5.500 + dlogg),
+            (3.447 - dlogT, -0.51 - dlogg),
+            (3.677 + dlogT, -0.51 - dlogg),
+            (3.740 + dlogT, 0.000 - dlogg),
+            (3.760 + dlogT, 0.500 - dlogg),
+            (3.845 + dlogT, 1.000 - dlogg),
+            (3.845 + dlogT, 1.500 - dlogg),
+            (4.079 + dlogT, 2.000 - dlogg),
+            (4.079 + dlogT, 2.500 - dlogg),
+            (4.204 + dlogT, 3.000 - dlogg),
+            (4.204 + dlogT, 5.000 + dlogg),
+            (3.602 + dlogT, 5.500 + dlogg),
+            (3.447 - dlogT, 5.500 + dlogg),
         ]
 
         return np.array(bbox)
@@ -2036,10 +2042,6 @@ def logT(self):
     def Z(self):
         return self.grid["Z"]
 
-    @property
-    def logZ(self):
-        return np.log10(self.grid["Z"])
-
 
 class Tlusty2025(Stellib):
     """
@@ -2092,16 +2094,18 @@ def bbox(self, dlogT=0.05, dlogg=0.25):
             (logT, logg) edges of the bounding polygon
         """
         bbox = [
-            (4.176 - dlogT, 4.749 + dlogg),
+            (4.176 - dlogT, 4.760 + dlogg),
             (4.176 - dlogT, 1.750 - dlogg),
-            (4.176 + dlogT, 1.750 - dlogg),
-            (4.255 + dlogT, 2.000 - dlogg),
-            (4.447 + dlogT, 2.750 - dlogg),
-            (4.478 + dlogT, 3.000 - dlogg),
-            (4.544 + dlogT, 3.250 - dlogg),
+            (4.278 + dlogT, 2.000 - dlogg),
+            (4.342 + dlogT, 2.250 - dlogg),
+            (4.380 + dlogT, 2.500 - dlogg),
+            (4.544 + dlogT, 3.000 - dlogg),
+            (4.574 + dlogT, 3.250 - dlogg),
+            (4.677 + dlogT, 3.500 - dlogg),
+            (4.699 + dlogT, 3.750 - dlogg),
             (4.740 + dlogT, 4.000 - dlogg),
-            (4.740 + dlogT, 4.749 + dlogg),
-            (4.176 - dlogT, 4.749 + dlogg),
+            (4.740 + dlogT, 4.760 + dlogg),
+            (4.176 - dlogT, 4.760 + dlogg),
         ]
 
         return np.array(bbox)
@@ -2121,7 +2125,3 @@ def logT(self):
     @property
     def Z(self):
         return self.grid["Z"]
-
-    @property
-    def logZ(self):
-        return np.log10(self.grid["Z"])

From 439a3dc9d2f402d4704f7fd958e63343671defc8 Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Sun, 24 May 2026 07:54:23 -0400
Subject: [PATCH 17/32] found a bug!

---
 .../stars/stellar_libraries/gen_bosz24.py     | 12 ++++-
 .../stars/stellar_libraries/gen_tlusty25.py   | 54 +++++++++----------
 beast/physicsmodel/stars/stellib.py           | 28 ++++++----
 3 files changed, 56 insertions(+), 38 deletions(-)

diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
index aaae2161..e9cf4aa4 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
@@ -59,6 +59,8 @@ def decode_params(filename):
     )
     # fmt: on
 
+    kk = 0
+    files = files[0:2]
     for k, cfile in enumerate(files):
         print(cfile)
 
@@ -118,8 +120,14 @@ def decode_params(filename):
 
             if k == 0:
                 outspec = np.zeros((len(files) + 1, len(wave_rebin)))
-                outspec[-1, :] = wave_rebin
-            outspec[k, :] = flux_rebin * 4.0 * np.pi
+            outspec[kk, :] = flux_rebin * 4.0 * np.pi
+            kk += 1
+
+    # trim the grid to remove all the models not included
+    outspec = outspec[0:kk+2, :]
+    outspec[-1, :] = wave_rebin
+    print(outspec[kk+1, :])
+    exit()
 
     # output the stellar library in the beast format
     spec_hdu = fits.PrimaryHDU(data=outspec)
diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
index a8154140..b2147c05 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
@@ -59,33 +59,33 @@ def decode_params(filename):
     for k, cfile in enumerate(files):
         print(cfile)
 
-        # # read in the model spectrum
-        # # wave units are angstrom
-        # # flux units are H, F = 4*pi*H in ergs/(cm^2 s A)
-        # mspec = ascii.read(
-        #     cfile,
-        #     format="no_header",
-        #     fast_reader={"exponent_style": "D"},
-        #     names=["Wave", "SFlux"],
-        # )
-
-        # # convert the type to float
-        # mspec["SFlux"] = mspec["SFlux"].astype(float)
-
-        # # now extract the wave and flux columns
-        # mwave = mspec["Wave"]
-        # mflux = mspec["SFlux"]
-
-        # # rebin to R=4000 for speed, common res and wave range with BOSZ LTE models
-        # rbres = 2000.0
-        # wave_rebin, flux_rebin, npts_rebin = rebin_spectrum(
-        #     mwave.value, mflux.value, rbres, [500.0, 320000.0]
-        # )
-
-        # if k == 0:
-        #     outspec = np.zeros((len(files) + 1, len(wave_rebin)))
-        #     outspec[-1, :] = wave_rebin
-        # outspec[k, :] = flux_rebin * 4.0 * np.pi
+        # read in the model spectrum
+        # wave units are angstrom
+        # flux units are H, F = 4*pi*H in ergs/(cm^2 s A)
+        mspec = ascii.read(
+            cfile,
+            format="no_header",
+            fast_reader={"exponent_style": "D"},
+            names=["Wave", "SFlux"],
+        )
+
+        # convert the type to float
+        mspec["SFlux"] = mspec["SFlux"].astype(float)
+
+        # now extract the wave and flux columns
+        mwave = mspec["Wave"]
+        mflux = mspec["SFlux"]
+
+        # rebin to R=4000 for speed, common res and wave range with BOSZ LTE models
+        rbres = 2000.0
+        wave_rebin, flux_rebin, npts_rebin = rebin_spectrum(
+            mwave.value, mflux.value, rbres, [500.0, 320000.0]
+        )
+
+        if k == 0:
+            outspec = np.zeros((len(files) + 1, len(wave_rebin)))
+            outspec[-1, :] = wave_rebin
+        outspec[k, :] = flux_rebin * 4.0 * np.pi
 
         model_params = decode_params(cfile)
         Z = model_params["Z"] * solar_z
diff --git a/beast/physicsmodel/stars/stellib.py b/beast/physicsmodel/stars/stellib.py
index fbda9004..48b7de5a 100644
--- a/beast/physicsmodel/stars/stellib.py
+++ b/beast/physicsmodel/stars/stellib.py
@@ -8,6 +8,8 @@
 (this may not be pythonic though)
 """
 
+import matplotlib.pyplot as plt
+
 import numpy as np
 from scipy.interpolate import interp1d
 from numpy.lib import recfunctions
@@ -331,7 +333,7 @@ def interp(self, T0, g0, Z0, L0, dT_max=0.1, eps=1e-6):
             i2 (resp. i1) is not used.
             (see below for namings)
 
-        eps: foat
+        eps: float
             temperature sensitivity under which points are considered to
             have the same temperature
 
@@ -702,15 +704,23 @@ def gen_spectral_grid_from_given_points(
         ):
             if bound_cond[mk]:
                 s = np.array(self.interp(rT, rg, rZ, 0.0)).T
-                # print("logT, logg, Z")
-                # print(rT, rg, rZ)
-                # print(s)
-                # print("logT", self.grid["logT"][s[:, 0].astype(int)].data)
-                # print("logg", self.grid["logg"][s[:, 0].astype(int)].data)
-                # print("Z", self.grid["Z"][s[:, 0].astype(int)].data)
+                print("logT, logg, Z")
+                print(rT, rg, rZ)
+                print(s)
+                print("logT", self.grid["logT"][s[:, 0].astype(int)].data)
+                print("logg", self.grid["logg"][s[:, 0].astype(int)].data)
+                print("Z", self.grid["Z"][s[:, 0].astype(int)].data)
+
+                indxs = s[:, 0].astype(int)
+                for kkk in indxs:
+                    plt.plot(self.wavelength, self.spectra[kkk, :])
+                plt.xscale("log")
+                plt.yscale("log")
+                plt.show()
+
                 specs[mk, :] = self.genSpectrum(s) * weights[mk]
-                # print(specs[mk, :])
-                # exit()
+                print(specs[mk, :])
+                exit()
 
         # Step 4: filter points without spectrum
         # ======================================

From 2f8c6de4b92e8bd64774df0f725dd636e415e5cc Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Sun, 24 May 2026 08:00:33 -0400
Subject: [PATCH 18/32] fixed bosz 2024 grid generation

---
 beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
index e9cf4aa4..d2dafe79 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
@@ -124,10 +124,10 @@ def decode_params(filename):
             kk += 1
 
     # trim the grid to remove all the models not included
-    outspec = outspec[0:kk+2, :]
+    outspec = outspec[0:kk+1, :]
+
+    # add the wavelengths to the end
     outspec[-1, :] = wave_rebin
-    print(outspec[kk+1, :])
-    exit()
 
     # output the stellar library in the beast format
     spec_hdu = fits.PrimaryHDU(data=outspec)

From a024d691e9de358e0bab87be4c4dcc99a1148deb Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Sun, 24 May 2026 08:02:54 -0400
Subject: [PATCH 19/32] removing debug statement

---
 beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
index d2dafe79..f6953903 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
@@ -60,7 +60,6 @@ def decode_params(filename):
     # fmt: on
 
     kk = 0
-    files = files[0:2]
     for k, cfile in enumerate(files):
         print(cfile)
 

From 9f5aa9d70d561c3cfdbb53a8c59ab387939b67f4 Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Sun, 24 May 2026 20:16:15 -0400
Subject: [PATCH 20/32] works!  Tweak to the new tlusty2025 metallicities.

---
 .../stars/stellar_libraries/check_tlusty25.py | 27 +++++++-----
 beast/physicsmodel/stars/stellib.py           | 41 ++++++++++++-------
 2 files changed, 43 insertions(+), 25 deletions(-)

diff --git a/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py b/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py
index 4bd6363c..faed95f6 100644
--- a/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py
+++ b/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py
@@ -7,14 +7,21 @@
 if __name__ == "__main__":  # pragma: no cover
 
     slib = stellib.Tlusty2025()
-    bound = get_stellib_boundaries(slib, dlogT=0.0, dlogg=0.0)
-    #bbox = slib.get_boundaries(dlogT=0.01, dlogg=0.01)
-    #print(bbox)
-    print(bound)
+    # slib = stellib.Tlusty()
 
-    plt.plot(slib.logT, slib.logg, "ko")
-    plt.plot(bound[:, 0], bound[:, 1], "b-", label="new")
-    slib.plot_boundary(dlogT=0.0, dlogg=0.0)
-    #plt.plot(bbox[:, 0], bbox[:, 1], "g-", label="in lib")
-    plt.legend()
-    plt.show()
+    print(np.unique(slib.Z))
+
+    for cz in np.unique(slib.Z.data):
+        bound = get_stellib_boundaries(slib, dlogT=0.0, dlogg=0.0)
+        #bbox = slib.get_boundaries(dlogT=0.01, dlogg=0.01)
+        #print(bbox)
+        #print(bound)
+
+        gvals = slib.Z.data == cz
+        plt.plot(slib.logT[gvals], slib.logg[gvals], "ko")
+        plt.plot(bound[:, 0], bound[:, 1], "b-", label="new")
+        slib.plot_boundary(dlogT=0.0, dlogg=0.0)
+        #plt.plot(bbox[:, 0], bbox[:, 1], "g-", label="in lib")
+        plt.title(cz)
+        plt.legend()
+        plt.show()
diff --git a/beast/physicsmodel/stars/stellib.py b/beast/physicsmodel/stars/stellib.py
index 48b7de5a..1a304bfd 100644
--- a/beast/physicsmodel/stars/stellib.py
+++ b/beast/physicsmodel/stars/stellib.py
@@ -704,23 +704,23 @@ def gen_spectral_grid_from_given_points(
         ):
             if bound_cond[mk]:
                 s = np.array(self.interp(rT, rg, rZ, 0.0)).T
-                print("logT, logg, Z")
-                print(rT, rg, rZ)
-                print(s)
-                print("logT", self.grid["logT"][s[:, 0].astype(int)].data)
-                print("logg", self.grid["logg"][s[:, 0].astype(int)].data)
-                print("Z", self.grid["Z"][s[:, 0].astype(int)].data)
-
-                indxs = s[:, 0].astype(int)
-                for kkk in indxs:
-                    plt.plot(self.wavelength, self.spectra[kkk, :])
-                plt.xscale("log")
-                plt.yscale("log")
-                plt.show()
+                # print("logT, logg, Z")
+                # print(rT, rg, rZ)
+                # print(s)
+                # print("logT", self.grid["logT"][s[:, 0].astype(int)].data)
+                # print("logg", self.grid["logg"][s[:, 0].astype(int)].data)
+                # print("Z", self.grid["Z"][s[:, 0].astype(int)].data)
+
+                # indxs = s[:, 0].astype(int)
+                # for kkk in indxs:
+                #     plt.plot(self.wavelength, self.spectra[kkk, :])
+                # plt.xscale("log")
+                # plt.yscale("log")
+                # plt.show()
 
                 specs[mk, :] = self.genSpectrum(s) * weights[mk]
-                print(specs[mk, :])
-                exit()
+                # print(specs[mk, :])
+                # exit()
 
         # Step 4: filter points without spectrum
         # ======================================
@@ -930,6 +930,7 @@ def which_osl(self, xypoints, dlogT=0.0, dlogg=0.0):
                     ind = np.atleast_1d(np.squeeze(np.where(res == 0)))
                     r = ok.points_inside(xy[ind], dlogT=dlogT, dlogg=dlogg)
                     res[ind[r]] = ek + 1
+
         return res
 
     def __repr__(self):
@@ -2087,6 +2088,16 @@ def _load_(self):
         self.header["NAME"] = "tlusty2025"
         self.spectra = g.seds
 
+        # the logT covearge at Z=0.000544 stops at logT = 4.8
+        # the logT coverage at Z=0.000561 starts at logT = 4.44
+        # all the other Z values have the expected logT,logg coverage
+        # To avoid issues in the interpolation routine, the metallicities
+        # for both these values are set to the average of Z=0.000553.
+        gvals = self.grid["Z"] == 0.000544
+        self.grid["Z"][gvals] = 0.000553
+        gvals = self.grid["Z"] == 0.000561
+        self.grid["Z"][gvals] = 0.000553
+
     def bbox(self, dlogT=0.05, dlogg=0.25):
         """Boundary of Tlusty library
 

From 6c964ef99ddcf7170bc9d47ae740a14fec45d0d0 Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Mon, 25 May 2026 12:09:38 -0400
Subject: [PATCH 21/32] cleanup

---
 ...eck_tlusty25.py => check_slib_coverage.py} |  0
 .../stellar_libraries/comp_tlusty25_bosz24.py |  1 -
 .../stars/stellar_libraries/gen_bosz24.py     |  2 +-
 .../stars/stellar_libraries/gen_tlusty25.py   | 11 +++++--
 beast/physicsmodel/stars/stellib.py           | 30 +++----------------
 5 files changed, 13 insertions(+), 31 deletions(-)
 rename beast/physicsmodel/stars/stellar_libraries/{check_tlusty25.py => check_slib_coverage.py} (100%)

diff --git a/beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py b/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
similarity index 100%
rename from beast/physicsmodel/stars/stellar_libraries/check_tlusty25.py
rename to beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
diff --git a/beast/physicsmodel/stars/stellar_libraries/comp_tlusty25_bosz24.py b/beast/physicsmodel/stars/stellar_libraries/comp_tlusty25_bosz24.py
index c3796005..08c6cfb0 100644
--- a/beast/physicsmodel/stars/stellar_libraries/comp_tlusty25_bosz24.py
+++ b/beast/physicsmodel/stars/stellar_libraries/comp_tlusty25_bosz24.py
@@ -5,7 +5,6 @@
 
 from helpers import rebin_spectrum
 
-
 if __name__ == "__main__":  # pragma: no cover
 
     solar_z = 0.02
diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
index f6953903..c9bbe66a 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
@@ -123,7 +123,7 @@ def decode_params(filename):
             kk += 1
 
     # trim the grid to remove all the models not included
-    outspec = outspec[0:kk+1, :]
+    outspec = outspec[0 : kk + 1, :]
 
     # add the wavelengths to the end
     outspec[-1, :] = wave_rebin
diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
index b2147c05..535ca4a5 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
@@ -89,9 +89,14 @@ def decode_params(filename):
 
         model_params = decode_params(cfile)
         Z = model_params["Z"] * solar_z
-        row = (Z, model_params["Teff"], model_params["logg"], model_params["vturb"],
-               np.log10(model_params["Teff"]),
-               model_params["Z"])
+        row = (
+            Z,
+            model_params["Teff"],
+            model_params["logg"],
+            model_params["vturb"],
+            np.log10(model_params["Teff"]),
+            model_params["Z"],
+        )
         print(row)
         outtab.add_row(row)
 
diff --git a/beast/physicsmodel/stars/stellib.py b/beast/physicsmodel/stars/stellib.py
index 1a304bfd..7e97b1b4 100644
--- a/beast/physicsmodel/stars/stellib.py
+++ b/beast/physicsmodel/stars/stellib.py
@@ -8,8 +8,6 @@
 (this may not be pythonic though)
 """
 
-import matplotlib.pyplot as plt
-
 import numpy as np
 from scipy.interpolate import interp1d
 from numpy.lib import recfunctions
@@ -670,10 +668,8 @@ def gen_spectral_grid_from_given_points(
                 for key in pts.dtype.names:
                     _grid[key] = pts[key]
             else:
-                raise AttributeError(
-                    "pts is expected to have named items \
-                                     (keys or dtype.names)"
-                )
+                raise AttributeError("pts is expected to have named items \
+                                     (keys or dtype.names)")
 
         # Step 1: Avoid Extrapolation
         # ===========================
@@ -704,23 +700,7 @@ def gen_spectral_grid_from_given_points(
         ):
             if bound_cond[mk]:
                 s = np.array(self.interp(rT, rg, rZ, 0.0)).T
-                # print("logT, logg, Z")
-                # print(rT, rg, rZ)
-                # print(s)
-                # print("logT", self.grid["logT"][s[:, 0].astype(int)].data)
-                # print("logg", self.grid["logg"][s[:, 0].astype(int)].data)
-                # print("Z", self.grid["Z"][s[:, 0].astype(int)].data)
-
-                # indxs = s[:, 0].astype(int)
-                # for kkk in indxs:
-                #     plt.plot(self.wavelength, self.spectra[kkk, :])
-                # plt.xscale("log")
-                # plt.yscale("log")
-                # plt.show()
-
                 specs[mk, :] = self.genSpectrum(s) * weights[mk]
-                # print(specs[mk, :])
-                # exit()
 
         # Step 4: filter points without spectrum
         # ======================================
@@ -1251,10 +1231,8 @@ def gen_spectral_grid_from_given_points(
                 elif hasattr(pts, "dtype"):
                     keys = pts.dtypes.names
                 else:
-                    raise AttributeError(
-                        "Input pts is expected to have \
-                                         named fields"
-                    )
+                    raise AttributeError("Input pts is expected to have \
+                                         named fields")
                 for k in keys:
                     _pts[k] = pts[k][ind]
                 # keep track of the spectra library that is selected

From 7fd8448eda8ee2c0ec55cca809e7ffeceff658c5 Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Mon, 25 May 2026 12:14:51 -0400
Subject: [PATCH 22/32] last codestlye

---
 .../stars/stellar_libraries/check_slib_coverage.py        | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py b/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
index faed95f6..ceac249e 100644
--- a/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
+++ b/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
@@ -13,15 +13,15 @@
 
     for cz in np.unique(slib.Z.data):
         bound = get_stellib_boundaries(slib, dlogT=0.0, dlogg=0.0)
-        #bbox = slib.get_boundaries(dlogT=0.01, dlogg=0.01)
-        #print(bbox)
-        #print(bound)
+        # bbox = slib.get_boundaries(dlogT=0.01, dlogg=0.01)
+        # print(bbox)
+        # print(bound)
 
         gvals = slib.Z.data == cz
         plt.plot(slib.logT[gvals], slib.logg[gvals], "ko")
         plt.plot(bound[:, 0], bound[:, 1], "b-", label="new")
         slib.plot_boundary(dlogT=0.0, dlogg=0.0)
-        #plt.plot(bbox[:, 0], bbox[:, 1], "g-", label="in lib")
+        # plt.plot(bbox[:, 0], bbox[:, 1], "g-", label="in lib")
         plt.title(cz)
         plt.legend()
         plt.show()

From 5db671acedab1da9d93151133fa4479c2b7f5e52 Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Mon, 25 May 2026 16:04:46 -0400
Subject: [PATCH 23/32] minor doc update

---
 beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py | 8 +-------
 1 file changed, 1 insertion(+), 7 deletions(-)

diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
index c9bbe66a..33c4c15c 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
@@ -11,7 +11,7 @@
 
 def decode_params(filename):
     """
-    Decode the tlusty filenames for the model parameters.
+    Decode the bosz filenames for the model parameters.
 
     Parameters
     ----------
@@ -28,12 +28,6 @@ def decode_params(filename):
     tpos = filename.find("_t", slashpos)
     vpos = filename.find("_v", slashpos)
 
-    # print(filename[zpos + 2 : zpos + 7])
-    # print(filename[tpos + 2 : gpos])
-    # print(filename[gpos + 2 : gpos + 6])
-    # print(filename[vpos + 2 : vpos + 3])
-    # exit()
-
     model_params["Z"] = float(filename[zpos + 2 : zpos + 7])
     model_params["Teff"] = float(filename[tpos + 2 : gpos])
     model_params["logg"] = float(filename[gpos + 2 : gpos + 6])

From 10e72b616d3fbc1485764b91f0ae2149daab4998 Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Mon, 25 May 2026 17:35:00 -0400
Subject: [PATCH 24/32] adding aringer 2016

---
 .../stellar_libraries/check_slib_coverage.py  |   1 +
 .../stars/stellar_libraries/gen_aringer16.py  | 184 ++++++++++++++++++
 .../stars/stellar_libraries/gen_bosz24.py     |   2 +-
 .../stars/stellar_libraries/gen_tlusty25.py   |   2 +-
 beast/physicsmodel/stars/stellib.py           |   6 +-
 5 files changed, 190 insertions(+), 5 deletions(-)
 create mode 100644 beast/physicsmodel/stars/stellar_libraries/gen_aringer16.py

diff --git a/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py b/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
index ceac249e..8c4e730b 100644
--- a/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
+++ b/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
@@ -8,6 +8,7 @@
 
     slib = stellib.Tlusty2025()
     # slib = stellib.Tlusty()
+    slib = stellib.Aringer2016()
 
     print(np.unique(slib.Z))
 
diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_aringer16.py b/beast/physicsmodel/stars/stellar_libraries/gen_aringer16.py
new file mode 100644
index 00000000..11225254
--- /dev/null
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_aringer16.py
@@ -0,0 +1,184 @@
+import matplotlib.pyplot as plt
+
+import numpy as np
+
+from astropy.table import QTable
+from astropy.io import fits
+from astropy.io import ascii
+from astropy import units as u
+from astropy.modeling import models
+
+from helpers import rebin_spectrum
+
+
+def decode_params(filename):
+    """
+    Decode the tlusty filenames for the model parameters.
+
+    Parameters
+    ----------
+    filename : str
+        name of file to decode
+
+    """
+    model_params = {}
+
+    slashpos = filename.rfind("/")
+
+    tpos = filename.find("t", slashpos)
+    gpos = filename.find("g", slashpos)
+    masspos = filename.find("_m", slashpos)
+    carbonpos = filename.find("_c", slashpos)
+    oxygenpos = filename.find("_o", slashpos)
+    ironpos = filename.find("_fe", slashpos)
+    vpos = filename.find("_xi", slashpos)
+    cpos = filename.find("_c", vpos)
+
+    model_params["Teff"] = float(filename[tpos + 1 : gpos - 1])
+    model_params["logg"] = float(filename[gpos + 1 : masspos]) * 0.01
+    model_params["vturb"] = float(filename[vpos + 3 : cpos]) * 0.1
+    model_params["Mass"] = float(filename[masspos + 2 : oxygenpos]) * 0.01
+    model_params["logO"] = np.log10(float(f"0.{filename[oxygenpos + 2 : carbonpos]}"))
+    model_params["logC"] = np.log10(float(f"0.{filename[carbonpos + 2 : ironpos]}"))
+    model_params["C_to_O"] = (10 ** (model_params["logC"] - 12)) / (
+        10 ** (model_params["logO"] - 12)
+    )
+
+    # metallicity from Fe (comments from Martha's code)
+
+    # Aringer+2019 says "Since we use here the common definition, where log(eps(H)) is set
+    # constant to 12, and no COMARCS models with special deviations of the Fe content exist,
+    # the quantities [Z/H], [Fe/H], and log(Z/Zsun) are always equal."  So then all that
+    # matters is what solar value I use.  The latest one that the CMD interface uses
+    # is 0.0207... maybe I should use that?
+
+    # SEE UPDATE!  KEEPING THIS TEXT FOR POSTERITY. "For the solar composition we adopted the
+    # values from Anders & Grevesse (1989) except for C, N and O where we took the data from
+    # Grevesse & Sauval (1994). This is in agreement with our previous work
+    # (e.g. Aringer et al. 1999) and results in a Z⊙ of approximately
+    # 0.02." --> So I will use 0.01922
+
+    # UPDATE !!!!!! Asplund et al. 2021 say Zsun/Xsun = 0.0187.  Lodders et al. (2019)
+    # say 0.0200.  Leo's CMD page says 0.0207.  I will use 0.02 here b/c Bernhard mentioned
+    # 0.02 in his 2016 text AND because if you do the math on the Kurucz models, they use
+    # Zsun=0.02 exactly.
+
+    Fenum = float(f"0.{filename[ironpos + 3 : vpos]}") * 10.0
+    # solar value is listed as 7.56 in Aringer+2019 (MNRAS, 487,2133)
+    Fesun = 10.0 ** (7.56 - 12)
+    Fe = 10.0 ** (Fenum - 12)
+    FeH = np.round(np.log10(Fe) - np.log10(Fesun), 1)
+    Z = (10.0**FeH) * 0.02  # see comments below
+    model_params["Z"] = (10.0**FeH) * 0.02  # see comments below
+
+    return model_params
+
+
+if __name__ == "__main__":  # pragma: no cover
+
+    solar_z = 0.02  # see comments in decode_params
+
+    # get all vtrub=2 models
+    path = "/home/kgordon/Python/extstar_data/Models/Aringer16/"
+    mspec = ascii.read(f"{path}/mstar_spec_list.dat", data_start=0, format="no_header")
+    files = [f"{path}{cfile}" for cfile in mspec["col1"].data]
+
+    n_files = len(files)
+
+    # fmt: off
+    outtab = QTable(
+        names=["Z", "Teff", "logg", "vtrub", "logT", "logz", "Mass", "logO", "logC", "R_Rs"],
+        dtype=["f", "f", "f", "f", "f", "f", "f", "f", "f", "f"]
+    )
+    # fmt: on
+
+    for k, cfile in enumerate(files):
+        print(cfile)
+
+        model_params = decode_params(cfile)
+        Z = model_params["Z"] * solar_z
+
+        # will also need the radius assumed when creating the spectrum
+        # calculated from log g = log M − 2 log R + 4.437
+
+        # radius in solar units
+        R_Rs = 10.0 ** (
+            (model_params["logg"] - 4.437 - np.log10(model_params["Mass"])) / (-2.0)
+        )
+        # radius in cm
+        Reff = R_Rs * 6.957e10  # that's sun's radius in cm
+
+        row = (
+            Z,
+            model_params["Teff"],
+            model_params["logg"],
+            model_params["vturb"],
+            np.log10(model_params["Teff"]),
+            model_params["Z"],
+            model_params["Mass"],
+            model_params["logO"],
+            model_params["logC"],
+            R_Rs,
+        )
+        print(row)
+        outtab.add_row(row)
+
+        # read in the model spectrum
+        # wave units are angstrom
+        # flux is in nu*L_nu [erg/s]
+        # (radius has been used, unlike other stellar atmosphere results)
+        mspec = ascii.read(
+            cfile,
+            names=["Wave", "d1", "SFlux", "d2"],
+        )
+
+        # convert the type to float
+        mspec["SFlux"] = mspec["SFlux"].astype(float)
+
+        # now extract the wave and flux columns
+        mwave = mspec["Wave"]
+        mflux = mspec["SFlux"]
+
+        # spectra was in nu*L_nu [erg/s].  This conversion puts in erg/s/cm^2/AA
+        mflux /= (4 * np.pi * (Reff**2.0)) / mwave
+
+        # rebin to R=2000 for speed, common res and wave range with BOSZ LTE models
+        rbres = 2000.0
+        wave_rebin, flux_rebin, npts_rebin = rebin_spectrum(
+            mwave.value, mflux.value, rbres, [500.0, 320000.0]
+        )
+
+        # add a scaled blackbody for the shorter and longer wavelengths that were not calculated
+        bb = models.BlackBody(temperature = model_params["Teff"]*u.K,
+                              scale=1.0 * u.erg / (u.cm * u.cm * u.s * u.AA * u.sr)
+                              )
+        modbb = bb(wave_rebin)
+
+        # shorter wavelengths
+        gvals = np.absolute(wave_rebin - 3370.0) < 10.0
+        modbb *= np.nanmean(flux_rebin[gvals]) / np.nanmean(modbb[gvals])
+        gvals = (wave_rebin < 5000.0) & (flux_rebin == 0.0)
+        flux_rebin[gvals] = modbb[gvals]
+
+        # longer wavelengths
+        gvals = np.absolute(wave_rebin - (250000.0 - 150.0)) < 100.0
+        modbb *= np.nanmean(flux_rebin[gvals]) / np.nanmean(modbb[gvals])
+        gvals = (wave_rebin > 5000.0) & (flux_rebin == 0.0)
+        flux_rebin[gvals] = modbb[gvals]
+
+        if k == 0:
+            outspec = np.zeros((len(files) + 1, len(wave_rebin)))
+            outspec[-1, :] = wave_rebin
+        outspec[k, :] = flux_rebin
+
+        plt.plot(wave_rebin, flux_rebin)
+    plt.xscale("log")
+    plt.yscale("log")
+    plt.show()
+
+    # output the stellar library in the beast format
+    spec_hdu = fits.PrimaryHDU(data=outspec)
+    table_hdu = fits.BinTableHDU(data=outtab)
+    table_hdu.name = "Aringer16"
+    hdul = fits.HDUList([spec_hdu, table_hdu])
+    hdul.writeto("aringer2016.grid.fits", overwrite=True)
diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
index 33c4c15c..6e4221e3 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
@@ -105,7 +105,7 @@ def decode_params(filename):
             mwave = mspec["Wave"]
             mflux = mspec["SFlux"]
 
-            # rebin to R=4000 for speed, common res and wave range with BOSZ LTE models
+            # rebin to R=2000 for speed, common res and wave range with BOSZ LTE models
             rbres = 2000.0
             wave_rebin, flux_rebin, npts_rebin = rebin_spectrum(
                 mwave.value, mflux.value, rbres, [500.0, 320000.0]
diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
index 535ca4a5..d12e0b69 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
@@ -76,7 +76,7 @@ def decode_params(filename):
         mwave = mspec["Wave"]
         mflux = mspec["SFlux"]
 
-        # rebin to R=4000 for speed, common res and wave range with BOSZ LTE models
+        # rebin to R=2000 for speed, common res and wave range with BOSZ LTE models
         rbres = 2000.0
         wave_rebin, flux_rebin, npts_rebin = rebin_spectrum(
             mwave.value, mflux.value, rbres, [500.0, 320000.0]
diff --git a/beast/physicsmodel/stars/stellib.py b/beast/physicsmodel/stars/stellib.py
index 7e97b1b4..ea4bee5f 100644
--- a/beast/physicsmodel/stars/stellib.py
+++ b/beast/physicsmodel/stars/stellib.py
@@ -33,7 +33,7 @@
     "btsettl": __ROOT__ + "bt-settl.lowres.grid.fits",
     "btsettl_medres": __ROOT__ + "bt-settl.medres.grid.fits",
     "munari": __ROOT__ + "atlas9-munari.hires.grid.fits",
-    "aringer": __ROOT__ + "Aringer.AGB.grid.fits",
+    "aringer": __ROOT__ + "aringer2016.grid.fits",
     "bosz2024": __ROOT__ + "bosz2024.grid.fits",
     "tlusty2025": __ROOT__ + "tlusty2025.grid.fits",
 }
@@ -47,7 +47,7 @@
     "Munari",
     "Elodie",
     "BaSeL",
-    "Aringer",
+    "Aringer2016",
     "BOSZ2024",
     "Tlusty2025",
 ]
@@ -1850,7 +1850,7 @@ def logZ(self):
         return self.grid["logZ"]
 
 
-class Aringer(Stellib):
+class Aringer2016(Stellib):
     """Aringer C+M+K giants Library
 
     References

From 4d85c455fc0daa6b40f1fbd076deb6febe5d94fa Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Tue, 26 May 2026 07:24:34 -0400
Subject: [PATCH 25/32] cleanup

---
 .../stellar_libraries/check_slib_coverage.py   |  2 +-
 .../stars/stellar_libraries/gen_aringer16.py   | 18 +++++++++---------
 2 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py b/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
index 8c4e730b..0e91c420 100644
--- a/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
+++ b/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
@@ -20,7 +20,7 @@
 
         gvals = slib.Z.data == cz
         plt.plot(slib.logT[gvals], slib.logg[gvals], "ko")
-        plt.plot(bound[:, 0], bound[:, 1], "b-", label="new")
+        # plt.plot(bound[:, 0], bound[:, 1], "b-", label="new")
         slib.plot_boundary(dlogT=0.0, dlogg=0.0)
         # plt.plot(bbox[:, 0], bbox[:, 1], "g-", label="in lib")
         plt.title(cz)
diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_aringer16.py b/beast/physicsmodel/stars/stellar_libraries/gen_aringer16.py
index 11225254..96e5a822 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_aringer16.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_aringer16.py
@@ -68,8 +68,7 @@ def decode_params(filename):
     Fesun = 10.0 ** (7.56 - 12)
     Fe = 10.0 ** (Fenum - 12)
     FeH = np.round(np.log10(Fe) - np.log10(Fesun), 1)
-    Z = (10.0**FeH) * 0.02  # see comments below
-    model_params["Z"] = (10.0**FeH) * 0.02  # see comments below
+    model_params["Z"] = (10.0**FeH) * 0.02  # see comments above
 
     return model_params
 
@@ -149,9 +148,10 @@ def decode_params(filename):
         )
 
         # add a scaled blackbody for the shorter and longer wavelengths that were not calculated
-        bb = models.BlackBody(temperature = model_params["Teff"]*u.K,
-                              scale=1.0 * u.erg / (u.cm * u.cm * u.s * u.AA * u.sr)
-                              )
+        bb = models.BlackBody(
+            temperature=model_params["Teff"] * u.K,
+            scale=1.0 * u.erg / (u.cm * u.cm * u.s * u.AA * u.sr),
+        )
         modbb = bb(wave_rebin)
 
         # shorter wavelengths
@@ -171,10 +171,10 @@ def decode_params(filename):
             outspec[-1, :] = wave_rebin
         outspec[k, :] = flux_rebin
 
-        plt.plot(wave_rebin, flux_rebin)
-    plt.xscale("log")
-    plt.yscale("log")
-    plt.show()
+    #     plt.plot(wave_rebin, flux_rebin)
+    # plt.xscale("log")
+    # plt.yscale("log")
+    # plt.show()
 
     # output the stellar library in the beast format
     spec_hdu = fits.PrimaryHDU(data=outspec)

From ef27fe70e9eabb2d9a1cbc4fd50d38ae52aa2ea2 Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Tue, 26 May 2026 07:25:12 -0400
Subject: [PATCH 26/32] changing Aringer directory location

---
 beast/physicsmodel/stars/stellar_libraries/gen_aringer16.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_aringer16.py b/beast/physicsmodel/stars/stellar_libraries/gen_aringer16.py
index 96e5a822..522f4102 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_aringer16.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_aringer16.py
@@ -78,7 +78,8 @@ def decode_params(filename):
     solar_z = 0.02  # see comments in decode_params
 
     # get all vtrub=2 models
-    path = "/home/kgordon/Python/extstar_data/Models/Aringer16/"
+    # path = "/home/kgordon/Python/extstar_data/Models/Aringer16/"
+    path = "/astro/mboyer/Science/WINGS/Aringer16/"
     mspec = ascii.read(f"{path}/mstar_spec_list.dat", data_start=0, format="no_header")
     files = [f"{path}{cfile}" for cfile in mspec["col1"].data]
 

From f2c6c36243c93589bca7e973a11a65b7f8a2b9a0 Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Tue, 26 May 2026 09:14:30 -0400
Subject: [PATCH 27/32] improved plot

---
 .../stars/stellar_libraries/check_slib_coverage.py    | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py b/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
index 0e91c420..be8b0c62 100644
--- a/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
+++ b/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
@@ -6,9 +6,10 @@
 
 if __name__ == "__main__":  # pragma: no cover
 
-    slib = stellib.Tlusty2025()
-    # slib = stellib.Tlusty()
-    slib = stellib.Aringer2016()
+    # slib = stellib.Tlusty2025()
+    slib = stellib.BOSZ2024()
+    slib = stellib.Kurucz()
+    # slib = stellib.Aringer2016()
 
     print(np.unique(slib.Z))
 
@@ -24,5 +25,7 @@
         slib.plot_boundary(dlogT=0.0, dlogg=0.0)
         # plt.plot(bbox[:, 0], bbox[:, 1], "g-", label="in lib")
         plt.title(cz)
-        plt.legend()
+        plt.xlabel("log(Teff)")
+        plt.ylabel("log(g)")
+        # plt.legend()
         plt.show()

From 1d24f44510c27c11ba1e5c98f10e88f1f32bb3dc Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Tue, 26 May 2026 14:39:13 -0400
Subject: [PATCH 28/32] updating Aringer grid coverage to be more uniform

---
 .../stellar_libraries/check_slib_coverage.py  |  6 +--
 .../stars/stellar_libraries/gen_aringer16.py  | 38 ++++++++++++++-----
 beast/physicsmodel/stars/stellib.py           |  4 +-
 3 files changed, 32 insertions(+), 16 deletions(-)

diff --git a/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py b/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
index be8b0c62..008469b4 100644
--- a/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
+++ b/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
@@ -7,9 +7,9 @@
 if __name__ == "__main__":  # pragma: no cover
 
     # slib = stellib.Tlusty2025()
-    slib = stellib.BOSZ2024()
-    slib = stellib.Kurucz()
-    # slib = stellib.Aringer2016()
+    # slib = stellib.BOSZ2024()
+    # slib = stellib.Kurucz()
+    slib = stellib.Aringer2016()
 
     print(np.unique(slib.Z))
 
diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_aringer16.py b/beast/physicsmodel/stars/stellar_libraries/gen_aringer16.py
index 522f4102..79073087 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_aringer16.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_aringer16.py
@@ -1,5 +1,3 @@
-import matplotlib.pyplot as plt
-
 import numpy as np
 
 from astropy.table import QTable
@@ -81,9 +79,7 @@ def decode_params(filename):
     # path = "/home/kgordon/Python/extstar_data/Models/Aringer16/"
     path = "/astro/mboyer/Science/WINGS/Aringer16/"
     mspec = ascii.read(f"{path}/mstar_spec_list.dat", data_start=0, format="no_header")
-    files = [f"{path}{cfile}" for cfile in mspec["col1"].data]
-
-    n_files = len(files)
+    files = np.array([f"{path}{cfile}" for cfile in mspec["col1"].data])
 
     # fmt: off
     outtab = QTable(
@@ -92,6 +88,7 @@ def decode_params(filename):
     )
     # fmt: on
 
+    # first go though and get the parameters for all the models
     for k, cfile in enumerate(files):
         print(cfile)
 
@@ -123,6 +120,32 @@ def decode_params(filename):
         print(row)
         outtab.add_row(row)
 
+    # now remove the models that were for a special case for the paper
+    # this means the logg is not on a regular grid point
+    uvals, ucounts = np.unique(outtab["logg"], return_counts=True)
+
+    gvals_files = np.full(len(files), True)
+    gvals = ucounts < 5
+    for cval in uvals[gvals]:
+        bvals = cval == outtab["logg"]
+        gvals_files[bvals] = False
+
+    # now a cut on logT to remove two errant points
+    # and the slight extension for some loggs for only one metallicity
+    bvals = outtab["logT"] > 3.70
+    gvals_files[bvals] = False
+
+    # now cut the 7.981e-5 metallicity models
+    # only half the logT, logg region covered
+    bvals = np.absolute(np.log10(outtab["Z"]) - np.log10(7.981e-5)) < 0.1
+    gvals_files[bvals] = False
+
+    # update the output table and the file list
+    outtab = outtab[gvals_files]
+    files = files[gvals_files]
+
+    # now get the spectra
+    for k, cfile in enumerate(files):
         # read in the model spectrum
         # wave units are angstrom
         # flux is in nu*L_nu [erg/s]
@@ -172,11 +195,6 @@ def decode_params(filename):
             outspec[-1, :] = wave_rebin
         outspec[k, :] = flux_rebin
 
-    #     plt.plot(wave_rebin, flux_rebin)
-    # plt.xscale("log")
-    # plt.yscale("log")
-    # plt.show()
-
     # output the stellar library in the beast format
     spec_hdu = fits.PrimaryHDU(data=outspec)
     table_hdu = fits.BinTableHDU(data=outtab)
diff --git a/beast/physicsmodel/stars/stellib.py b/beast/physicsmodel/stars/stellib.py
index ea4bee5f..433c99fa 100644
--- a/beast/physicsmodel/stars/stellib.py
+++ b/beast/physicsmodel/stars/stellib.py
@@ -1874,7 +1874,7 @@ class Aringer2016(Stellib):
 
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
-        self.name = "Aringer"
+        self.name = "Aringer2016"
         self.source = config["aringer"]
         self._load_()
 
@@ -1906,8 +1906,6 @@ def bbox(self, dlogT=0.05, dlogg=0.25):
             (3.41497 - dlogT, 5.00 + dlogg),
             (3.69897 + dlogT, 5.00 + dlogg),
             (3.69897 + dlogT, 3.50 + dlogg),
-            (3.71600 + dlogT, 3.50 + dlogg),
-            (3.71600 + dlogT, 2.50 - dlogg),
             (3.69897 + dlogT, 2.50 - dlogg),
             (3.69897 + dlogT, 0.00 - dlogg),
             (3.57978 + dlogT, 0.00 - dlogg),

From 7b40d9a631cbcc568455dcbfa25c7954d139fb31 Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Tue, 26 May 2026 14:55:03 -0400
Subject: [PATCH 29/32] removing old kurucz generation code

---
 beast/physicsmodel/stars/kurucz/__init__.py |   0
 beast/physicsmodel/stars/kurucz/generate.py | 201 --------------------
 2 files changed, 201 deletions(-)
 delete mode 100644 beast/physicsmodel/stars/kurucz/__init__.py
 delete mode 100644 beast/physicsmodel/stars/kurucz/generate.py

diff --git a/beast/physicsmodel/stars/kurucz/__init__.py b/beast/physicsmodel/stars/kurucz/__init__.py
deleted file mode 100644
index e69de29b..00000000
diff --git a/beast/physicsmodel/stars/kurucz/generate.py b/beast/physicsmodel/stars/kurucz/generate.py
deleted file mode 100644
index e11048cb..00000000
--- a/beast/physicsmodel/stars/kurucz/generate.py
+++ /dev/null
@@ -1,201 +0,0 @@
-""" This module gives tools to generate Kurucz grid from original downloads """
-from beast.physicsmodel import grid
-from beast.physicsmodel.stars import stellib, isochrone
-from astropy.table import Table
-
-import astropy.io.fits as pyfits
-import numpy as np
-
-# import glob
-# from matplotlib.nxutils import points_inside_poly
-import matplotlib.path as mpltPath
-import sys
-
-
-def __treatSingleFile__(fname="ckp00_10000.fits"):
-    """ Grab the useful data from a single file
-
-    Used in Kurucz_to_Stellib
-
-    INPUTS:
-        fname   str file to process
-    OUTPUTS:
-        logg, teff, logz, lamb, data
-    """
-
-    with pyfits.open(fname) as f:
-        logz = f[0].header["log_Z"]
-        teff = f[0].header["TEFF"]
-        d = f[1].data
-        logg = [float(k[1:]) * 0.1 for k in d.dtype.names if k[0] == "g"]
-        d = np.array(d.tolist())
-        lamb = d[:, 0]
-        data = d[:, 1:].T
-        nspec = data.shape[0]
-        logz = np.asarray([logz] * nspec)
-        teff = np.asarray([teff] * nspec)
-    return logg, teff, logz, lamb, data
-
-
-def Kurucz_to_Stellib(lst):
-    """ Extract SED parameters and spectra from a single file and collapse
-    them into a MemoryGrid object
-
-    INPUTS:
-        lst list    list of files to process
-                e.g. lst = glob.glob('ck*/*fits')
-
-    OUTPUTS:
-        g   MemoryGrid  stellib grid in the original units
-    """
-    d = dict(logg=[], teff=[], logz=[])
-    specs = []
-
-    for kfile in lst:
-        logg, teff, logz, lamb, data = __treatSingleFile__(kfile)
-        d["logg"].append(logg)
-        d["teff"].append(teff)
-        d["logz"].append(logz)
-        specs.append(data)
-
-    # post process
-    for k in d:
-        d[k] = np.ravel(d[k])
-    specs = np.vstack(specs)
-    t = Table(d)
-    g = grid.MemoryGrid(lamb, specs, t)
-    return g
-
-
-def gen_spectral_grid_from_kurucz(outfile, osl, oiso, Z=0.02):
-    """ Reinterpolate a given stellar spectral library on to an Isochrone grid
-    INPUTS:
-        outfile     str         fits file to export to
-        osl     stellib.stellib     a stellar library
-        oiso        isochrone.Isochrone an isochrone library
-        Z       float           metallicity to use
-
-    OUTPUTS:
-        None
-
-        only write into outfile
-    """
-    if not (
-        grid.isNestedInstance(osl, stellib.Stellib)
-        and grid.isNestedInstance(oiso, isochrone.Isochrone)
-    ):
-        raise AssertionError()
-    specs = np.zeros((oiso.data.nrows + 1, len(osl.wavelength)), dtype=float)
-    specs[-1] = osl.wavelength[:]
-
-    def get_radius(logl, logt):
-        # get the radius of a star given its luminosity and temperature
-        # (assuming a black body)
-        lsun = 3.839e26  # W
-        sig = 5.67037321 * 1e-8  # W m**-2 K**-4
-        return np.sqrt((10 ** logl) * lsun / (4.0 * np.pi * sig * ((10 ** logt) ** 4)))
-
-    bounds = get_stellib_boundaries(osl, dlogT=0.1, dlogg=0.3, closed=True)
-
-    progress = 0
-    data = np.array([oiso.data["logg"], oiso.data["logT"]]).T
-    # converted to use matplotlib.path.Path.contains_points
-    #   (points_inside_poly deprecated and removed from numpy)
-    # bound_cond = points_inside_poly(data, bounds)
-    bounds_path = mpltPath.Path(bounds)
-    bound_cond = bounds_path.contains_points(data)
-
-    del data
-    radii = get_radius(oiso.data["logL"], oiso.data["logT"])
-    weights = (
-        4.0 * np.pi * (radii * 1e2) ** 2
-    )  # denorm models are in cm**-2 (4 * pi * rad)
-    for k in range(oiso.data.nrows):
-        p = int(100 * (k + 1) / oiso.data.nrows)
-        if progress < p:
-            progress = p
-            sys.stdout.write("progress... %d / 100\r" % progress)
-        if bound_cond[k] is True:
-            r = np.array(
-                osl.interp(oiso.data["logT"][k], oiso.data["logg"][k], Z, 0.0)
-            ).T
-            specs[k, :] = osl.genSpectrum(r) * weights[k]
-        else:
-            specs[k, :] = np.zeros(len(osl.wavelength), dtype=float)
-    sys.stdout.write("progress... %d / 100" % progress)
-
-    specs = specs[bound_cond, :]
-
-    pyfits.writeto(outfile, specs)
-
-    # copy pars
-    data = {}
-    for k in list(oiso.data.keys()):
-        data[k] = oiso.data[k][bound_cond]
-    data["radius"] = radii[bound_cond] / 6.955e8  # Rsun
-    pars = Table(data, name="Reinterpolated stellib grid")
-    pars.header["stellib"] = osl.source
-    pars.header["isoch"] = oiso.source
-
-    pars.write(outfile, append=True)
-
-
-def get_stellib_boundaries(s, dlogT=0.1, dlogg=0.3, closed=True):
-    """ Returns the closed boundary polygon around the stellar library with
-    given margins
-
-    INPUTS:
-        s   Stellib     Stellar library object
-
-    KEYWORDS:
-        dlogT   float       margin in logT
-        dlogg   float       margin in logg
-        closed  bool        if set, close the polygon
-
-    OUTPUTS:
-        b   ndarray[float, ndim=2]  (closed) boundary points: [logg, Teff]
-    """
-    #        use "points_inside_poly" to test wether a point is inside the limits
-    #        >>> data = np.array([iso.data['logg'], iso.data['logT']]).T
-    #        >>> aa = points_inside_poly(data, leftb)
-    #    """
-    leftb = [(k, np.max(s.logT[s.logg == k]) + dlogT) for k in np.unique(s.logg)]
-    leftb += [(leftb[-1][0] + dlogg, leftb[-1][1])]
-    leftb = [(leftb[0][0] - dlogg, leftb[0][1])] + leftb
-    rightb = [(k, np.min(s.logT[s.logg == k]) - dlogT) for k in np.unique(s.logg)[::-1]]
-    rightb += [(rightb[-1][0] - dlogg, rightb[-1][1])]
-    rightb = [(rightb[0][0] + dlogg, rightb[0][1])] + rightb
-    b = leftb + rightb
-    if closed:
-        b += [b[0]]
-    return np.array(b)
-
-
-# =========================== TEST UNITS ==================================
-# def test_stellib_boundaries():
-#    """ Test get_stellib_boundaries function """
-#    import pylab as plt
-#
-#    s = stellib.Kurucz()
-#    iso = isochrone.padova2010()
-#    leftb = get_stellib_boundaries(s, 0.1, 0.3, True)
-#
-#    plt.plot(s.Teff, s.logg, ".", color="k")
-#    plt.plot(leftb[:, 1], leftb[:, 0], "o-")
-#
-# leftb = [logg, teff]
-#    plt.plot(iso.data["logT"], iso.data["logg"], ",", color="r")
-#    data = np.array([iso.data["logg"], iso.data["logT"]]).T
-# Needs to be converted to use matplotlib.path.Path.contains_points
-# aa = points_inside_poly(data, leftb)
-#    plt.plot(iso.data["logT"][aa], iso.data["logg"][aa], ",", color="g")
-#    plt.xlabel("log(Teff)")
-#    plt.ylabel("log(g)")
-#    plt.xlim(plt.xlim()[::-1])
-#    plt.ylim(plt.ylim()[::-1])
-
-
-# def main_last_grid():
-#    s = stellib.Kurucz()
-#    iso = isochrone.padova2010()
-#    gen_spectral_grid_from_kurucz("tmp.fits", s, iso)

From 2681eb2614b63e504cc8d051b863138c4dac77ab Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Tue, 26 May 2026 15:17:30 -0400
Subject: [PATCH 30/32] adding coverage plot for multiple grids

---
 .../plot_multiple_slib_coverage.py            | 29 +++++++++++++++++++
 1 file changed, 29 insertions(+)
 create mode 100644 beast/physicsmodel/stars/stellar_libraries/plot_multiple_slib_coverage.py

diff --git a/beast/physicsmodel/stars/stellar_libraries/plot_multiple_slib_coverage.py b/beast/physicsmodel/stars/stellar_libraries/plot_multiple_slib_coverage.py
new file mode 100644
index 00000000..7995657f
--- /dev/null
+++ b/beast/physicsmodel/stars/stellar_libraries/plot_multiple_slib_coverage.py
@@ -0,0 +1,29 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+from beast.physicsmodel.stars import stellib
+
+if __name__ == "__main__":  # pragma: no cover
+
+    slibs = [stellib.Tlusty2025(),
+             stellib.BOSZ2024(),
+             stellib.Aringer2016()]
+    labels = ["Tlusty2025",
+              "BOSZ2025",
+              "Aringer2016"]
+    # slibs = [stellib.Tlusty(),
+    #          stellib.Kurucz()]
+    # labels = ["Tlusty",
+    #           "Kurucz"]
+    colors = ["r", "b", "g"]
+
+    plt.plot([1.0, 1.0], [1.1, 1.1])
+    for slib, ccol, clab in zip(slibs, colors, labels):
+        slib.plot_boundary(dlogT=0.0, dlogg=0.0, color=ccol, label=clab, alpha=0.5)
+
+    plt.xlim(3.3, 4.8)
+    plt.ylim(-1.5, 6.0)
+    plt.xlabel("log(Teff)")
+    plt.ylabel("log(g)")
+    plt.legend()
+    plt.show()
\ No newline at end of file

From eb624a92205ce9f4780f814ba5b96b0d4a15de4e Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Wed, 27 May 2026 10:42:41 -0400
Subject: [PATCH 31/32] better stellar lib coverage plots

---
 .../stellar_libraries/check_slib_coverage.py  | 61 ++++++++++++++-----
 1 file changed, 45 insertions(+), 16 deletions(-)

diff --git a/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py b/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
index 008469b4..b175fda2 100644
--- a/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
+++ b/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
@@ -1,3 +1,4 @@
+import argparse
 import numpy as np
 import matplotlib.pyplot as plt
 
@@ -6,26 +7,54 @@
 
 if __name__ == "__main__":  # pragma: no cover
 
-    # slib = stellib.Tlusty2025()
-    # slib = stellib.BOSZ2024()
-    # slib = stellib.Kurucz()
-    slib = stellib.Aringer2016()
+    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()
 
-    print(np.unique(slib.Z))
+    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()
 
-    for cz in np.unique(slib.Z.data):
-        bound = get_stellib_boundaries(slib, dlogT=0.0, dlogg=0.0)
-        # bbox = slib.get_boundaries(dlogT=0.01, dlogg=0.01)
-        # print(bbox)
-        # print(bound)
+    # 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], "b-", label="new")
-        slib.plot_boundary(dlogT=0.0, dlogg=0.0)
-        # plt.plot(bbox[:, 0], bbox[:, 1], "g-", label="in lib")
-        plt.title(cz)
+        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()
-        plt.show()
+        plt.legend()
+
+        if args.png:
+            plt.savefig(f"slib_coverage_{args.grid}_z_{cz:.2e}.png")
+            plt.close()
+        else:
+            plt.show()
\ No newline at end of file

From 3e8c96e0166a374046138ba9afd6857db80fb84a Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Wed, 27 May 2026 10:54:57 -0400
Subject: [PATCH 32/32] updating multi library coverage plot

---
 .../stellar_libraries/check_slib_coverage.py  |  6 ++-
 .../stars/stellar_libraries/gen_aringer16.py  |  2 +-
 .../stars/stellar_libraries/gen_bosz24.py     |  1 +
 .../stars/stellar_libraries/gen_tlusty25.py   |  1 +
 .../plot_multiple_slib_coverage.py            | 46 ++++++++++++++-----
 5 files changed, 41 insertions(+), 15 deletions(-)

diff --git a/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py b/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
index b175fda2..b1f71cd4 100644
--- a/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
+++ b/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
@@ -47,7 +47,9 @@
         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")
+        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)")
@@ -57,4 +59,4 @@
             plt.savefig(f"slib_coverage_{args.grid}_z_{cz:.2e}.png")
             plt.close()
         else:
-            plt.show()
\ No newline at end of file
+            plt.show()
diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_aringer16.py b/beast/physicsmodel/stars/stellar_libraries/gen_aringer16.py
index 79073087..2db858db 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_aringer16.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_aringer16.py
@@ -76,7 +76,7 @@ def decode_params(filename):
     solar_z = 0.02  # see comments in decode_params
 
     # get all vtrub=2 models
-    # path = "/home/kgordon/Python/extstar_data/Models/Aringer16/"
+    # change to the specific location the models have been downloaded to
     path = "/astro/mboyer/Science/WINGS/Aringer16/"
     mspec = ascii.read(f"{path}/mstar_spec_list.dat", data_start=0, format="no_header")
     files = np.array([f"{path}{cfile}" for cfile in mspec["col1"].data])
diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
index 6e4221e3..aff94876 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
@@ -42,6 +42,7 @@ def decode_params(filename):
     solar_z = 0.012
 
     # get all vtrub=2 models
+    # change to the specific location the models have been downloaded to
     path = "/home/kgordon/Python/extstar_data/Models/BOSZ2024/r5000/r5000/"
     files = glob.glob(f"{path}/*/*a+0.00_c+0.00_v2_r5000_resam.txt.gz")
     n_files = len(files)
diff --git a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
index d12e0b69..546e72e1 100644
--- a/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
@@ -45,6 +45,7 @@ def decode_params(filename):
     solar_z = 0.017  # from Grevasse & Sauval (1998) - TBC
 
     # get all vtrub=2 models
+    # change to the specific location the models have been downloaded to
     path = "/home/kgordon/Python/extstar_data/Models/Tlusty_2025/"
     files = glob.glob(f"{path}/*v2.spec.gz")
     n_files = len(files)
diff --git a/beast/physicsmodel/stars/stellar_libraries/plot_multiple_slib_coverage.py b/beast/physicsmodel/stars/stellar_libraries/plot_multiple_slib_coverage.py
index 7995657f..c0224e40 100644
--- a/beast/physicsmodel/stars/stellar_libraries/plot_multiple_slib_coverage.py
+++ b/beast/physicsmodel/stars/stellar_libraries/plot_multiple_slib_coverage.py
@@ -1,22 +1,36 @@
-import numpy as np
+import argparse
 import matplotlib.pyplot as plt
 
 from beast.physicsmodel.stars import stellib
 
 if __name__ == "__main__":  # pragma: no cover
 
-    slibs = [stellib.Tlusty2025(),
-             stellib.BOSZ2024(),
-             stellib.Aringer2016()]
-    labels = ["Tlusty2025",
-              "BOSZ2025",
-              "Aringer2016"]
-    # slibs = [stellib.Tlusty(),
-    #          stellib.Kurucz()]
-    # labels = ["Tlusty",
-    #           "Kurucz"]
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        "--old",
+        action="store_true",
+        help="Show the old Tlusty/Kurucz grids",
+    )
+    parser.add_argument("--png", action="store_true", help="save figure to a png file")
+    args = parser.parse_args()
+
+    if args.old:
+        slibs = [stellib.Tlusty(), stellib.Kurucz()]
+        labels = ["Tlusty", "Kurucz"]
+    else:
+        slibs = [stellib.Tlusty2025(), stellib.BOSZ2024(), stellib.Aringer2016()]
+        labels = ["Tlusty2025", "BOSZ2025", "Aringer2016"]
     colors = ["r", "b", "g"]
 
+    # 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)
+
     plt.plot([1.0, 1.0], [1.1, 1.1])
     for slib, ccol, clab in zip(slibs, colors, labels):
         slib.plot_boundary(dlogT=0.0, dlogg=0.0, color=ccol, label=clab, alpha=0.5)
@@ -26,4 +40,12 @@
     plt.xlabel("log(Teff)")
     plt.ylabel("log(g)")
     plt.legend()
-    plt.show()
\ No newline at end of file
+
+    if args.old:
+        extstr = "_old"
+    else:
+        extstr = ""
+    if args.png:
+        plt.savefig(f"slib_mulicoverage{extstr}.png")
+    else:
+        plt.show()