diff --git a/beast/config.py b/beast/config.py
index 2223818ae..a0c683433 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/kurucz/__init__.py b/beast/physicsmodel/stars/kurucz/__init__.py
deleted file mode 100644
index e69de29bb..000000000
diff --git a/beast/physicsmodel/stars/kurucz/generate.py b/beast/physicsmodel/stars/kurucz/generate.py
deleted file mode 100644
index e11048cbb..000000000
--- 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)
diff --git a/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py b/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
new file mode 100644
index 000000000..b1f71cd41
--- /dev/null
+++ b/beast/physicsmodel/stars/stellar_libraries/check_slib_coverage.py
@@ -0,0 +1,62 @@
+import argparse
+import numpy as np
+import matplotlib.pyplot as plt
+
+from beast.physicsmodel.stars import stellib
+from helpers import get_stellib_boundaries
+
+if __name__ == "__main__":  # pragma: no cover
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        "--grid",
+        default="BOSZ2024",
+        choices=["Tlusty2025", "BOSZ2024", "Aringer2016", "Kurucz", "Tlusty"],
+        help="Grid to show",
+    )
+    parser.add_argument("--png", action="store_true", help="save figures to png files")
+    args = parser.parse_args()
+
+    if args.grid == "BOSZ2024":
+        slib = stellib.BOSZ2024()
+    elif args.grid == "Tlusty2025":
+        slib = stellib.Tlusty2025()
+    elif args.grid == "Aringer2016":
+        slib = stellib.Aringer2026()
+    elif args.grid == "Kurucz":
+        slib = stellib.Kurucz()
+    elif args.grid == "Tlusty":
+        slib = stellib.Tlusty()
+
+    # useful when new grids are generarted
+    # results are used to define the bounding box in a model definition
+    bound = get_stellib_boundaries(slib, dlogT=0.0, dlogg=0.0)
+    print("Algorithmically generated bounds")
+    print(bound)
+
+    # setup the plots
+    fontsize = 12
+    font = {"size": fontsize}
+    plt.rc("font", **font)
+    plt.rc("lines", linewidth=2)
+    plt.rc("axes", linewidth=2)
+    plt.rc("xtick.major", width=2)
+    plt.rc("ytick.major", width=2)
+
+    for cz in np.unique(slib.Z.data):
+        gvals = slib.Z.data == cz
+        plt.plot(slib.logT[gvals], slib.logg[gvals], "ko")
+        plt.plot(bound[:, 0], bound[:, 1], "g-", label="algorithmly generated")
+        slib.plot_boundary(
+            dlogT=0.0, dlogg=0.0, label="defined in class def", alpha=0.5, color="b"
+        )
+        plt.title(f"{args.grid}; z = {cz:.2e}")
+        plt.xlabel("log(Teff)")
+        plt.ylabel("log(g)")
+        plt.legend()
+
+        if args.png:
+            plt.savefig(f"slib_coverage_{args.grid}_z_{cz:.2e}.png")
+            plt.close()
+        else:
+            plt.show()
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 000000000..08c6cfb04
--- /dev/null
+++ b/beast/physicsmodel/stars/stellar_libraries/comp_tlusty25_bosz24.py
@@ -0,0 +1,68 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import astropy.units as u
+from astropy.io import ascii
+
+from helpers import rebin_spectrum
+
+if __name__ == "__main__":  # pragma: no cover
+
+    solar_z = 0.02
+
+    # read the kurucz spectrum
+    # wavelength first
+    wave_filename = "/home/kgordon/Python/extstar_data/Models/BOSZ2024/r2000/bosz2024_wave_r2000.txt"
+    mspec_lte_wave = ascii.read(
+        wave_filename,
+        format="no_header",
+        # fast_reader={"exponent_style": "D"},
+        names=["Wave"],
+    )
+
+    model_filename = "/home/kgordon/Python/extstar_data/Models/BOSZ2024/r2000/m+0.00/bosz2024_ap_t15000_g+4.0_m+0.00_a+0.00_c+0.00_v2_r2000_resam.txt.gz"
+    mspec_lte = ascii.read(
+        model_filename,
+        format="no_header",
+        # fast_reader={"exponent_style": "D"},
+        names=["SFlux", "SCont"],
+    )
+    mspec_lte["Wave"] = mspec_lte_wave["Wave"] * u.angstrom
+
+    # read the tlusty spectrum
+    model_filename = (
+        "/home/kgordon/Python/extstar_data/Models/Tlusty_2023/z100t15000g400v2.spec.gz"
+    )
+    # read in the model spectrum
+    mspec = ascii.read(
+        model_filename,
+        format="no_header",
+        fast_reader={"exponent_style": "D"},
+        names=["Wave", "SFlux"],
+    )
+
+    # convert the type to float
+    mspec["SFlux"] = mspec["SFlux"].astype(float)
+
+    # set the units
+    fluxunit = u.erg / (u.s * u.cm * u.cm * u.angstrom)
+    mspec["Wave"].unit = u.angstrom
+    mspec["SFlux"].unit = fluxunit
+
+    # now extract the wave and flux colums
+    mwave = mspec["Wave"]
+    mflux = mspec["SFlux"]
+
+    # rebin to R=4000 for speed
+    rbres = 4000.0
+    wave_rebin, flux_rebin, npts_rebin = rebin_spectrum(
+        mwave.value, mflux.value, rbres, [200.0, 310000.0]
+    )
+    wave_rebin *= u.angstrom
+    print(len(wave_rebin))
+
+    # plot
+    plt.plot(wave_rebin.to(u.micron), flux_rebin * 4 * np.pi)
+    plt.plot(mspec_lte["Wave"].to(u.micron), mspec_lte["SFlux"] * 4 * np.pi)
+    plt.yscale("log")
+    plt.xscale("log")
+    plt.show()
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 000000000..2db858dbb
--- /dev/null
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_aringer16.py
@@ -0,0 +1,203 @@
+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)
+    model_params["Z"] = (10.0**FeH) * 0.02  # see comments above
+
+    return model_params
+
+
+if __name__ == "__main__":  # pragma: no cover
+
+    solar_z = 0.02  # see comments in decode_params
+
+    # get all vtrub=2 models
+    # 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])
+
+    # 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
+
+    # first go though and get the parameters for all the models
+    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)
+
+    # 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]
+        # (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
+
+    # 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
new file mode 100644
index 000000000..aff94876e
--- /dev/null
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_bosz24.py
@@ -0,0 +1,131 @@
+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 bosz filenames for the model parameters.
+
+    Parameters
+    ----------
+    filename : str
+        name of file to decode
+
+    """
+    model_params = {}
+
+    slashpos = filename.rfind("/")
+
+    gpos = filename.find("_g", slashpos)
+    zpos = filename.find("_m", gpos)
+    tpos = filename.find("_t", slashpos)
+    vpos = filename.find("_v", slashpos)
+
+    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
+
+    # from Asplund et al. (2005)/Grevasse et al. (2007) as referenced in BOSZ2024 paper
+    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)
+
+    # fmt: off
+    outtab = QTable(
+        names=["Z", "Teff", "logg", "vturb", "logT", "logz"],
+        dtype=["f", "f", "f", "f", "f", "f"]
+    )
+    # fmt: on
+
+    kk = 0
+    for k, cfile in enumerate(files):
+        print(cfile)
+
+        # get model parameters
+        model_params = decode_params(cfile)
+
+        # 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")
+
+        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"]),
+                model_params["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=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]
+            )
+
+            if k == 0:
+                outspec = np.zeros((len(files) + 1, len(wave_rebin)))
+            outspec[kk, :] = flux_rebin * 4.0 * np.pi
+            kk += 1
+
+    # trim the grid to remove all the models not included
+    outspec = outspec[0 : kk + 1, :]
+
+    # add the wavelengths to the end
+    outspec[-1, :] = wave_rebin
+
+    # 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
new file mode 100644
index 000000000..546e72e13
--- /dev/null
+++ b/beast/physicsmodel/stars/stellar_libraries/gen_tlusty25.py
@@ -0,0 +1,109 @@
+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("/")
+    periodpos = filename.rfind(".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.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)
+
+    # 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)
+
+        # 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=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]
+        )
+
+        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"]),
+            model_params["Z"],
+        )
+        print(row)
+        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 000000000..1b8d46b5f
--- /dev/null
+++ b/beast/physicsmodel/stars/stellar_libraries/helpers.py
@@ -0,0 +1,109 @@
+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)
+
+
+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))
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 000000000..c0224e405
--- /dev/null
+++ b/beast/physicsmodel/stars/stellar_libraries/plot_multiple_slib_coverage.py
@@ -0,0 +1,51 @@
+import argparse
+import matplotlib.pyplot as plt
+
+from beast.physicsmodel.stars import stellib
+
+if __name__ == "__main__":  # pragma: no cover
+
+    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)
+
+    plt.xlim(3.3, 4.8)
+    plt.ylim(-1.5, 6.0)
+    plt.xlabel("log(Teff)")
+    plt.ylabel("log(g)")
+    plt.legend()
+
+    if args.old:
+        extstr = "_old"
+    else:
+        extstr = ""
+    if args.png:
+        plt.savefig(f"slib_mulicoverage{extstr}.png")
+    else:
+        plt.show()
diff --git a/beast/physicsmodel/stars/stellib.py b/beast/physicsmodel/stars/stellib.py
index fd2e57400..433c99fa5 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
@@ -16,7 +17,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
@@ -33,7 +33,9 @@
     "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",
 }
 
 __all__ = [
@@ -45,26 +47,28 @@
     "Munari",
     "Elodie",
     "BaSeL",
-    "Aringer",
+    "Aringer2016",
+    "BOSZ2024",
+    "Tlusty2025",
 ]
 
 
 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__():
@@ -102,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
 
@@ -195,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.
@@ -282,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):
@@ -297,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.
@@ -327,7 +331,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
 
@@ -357,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
 
@@ -454,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 ),
@@ -485,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
@@ -540,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
 
@@ -560,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
@@ -664,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
         # ===========================
@@ -794,7 +796,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)
@@ -816,10 +818,10 @@ 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.asarray([osl.wavelength[:] for osl in self._olist]))
+        lambs = np.unique(np.concatenate([osl.wavelength[:] for osl in self._olist]))
         return lambs
 
     @property
@@ -908,6 +910,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):
@@ -916,7 +919,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
@@ -954,7 +957,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.
@@ -1023,7 +1026,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
 
@@ -1105,7 +1108,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
 
@@ -1130,7 +1133,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
@@ -1228,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
@@ -1277,7 +1278,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
 
@@ -1300,7 +1301,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
         ----------
@@ -1382,7 +1383,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
@@ -1411,7 +1412,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
         ----------
@@ -1514,7 +1515,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
@@ -1522,7 +1525,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
         ----------
@@ -1620,7 +1623,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
         ----------
@@ -1708,7 +1711,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
         ----------
@@ -1797,7 +1800,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
         ----------
@@ -1847,7 +1850,7 @@ def logZ(self):
         return self.grid["logZ"]
 
 
-class Aringer(Stellib):
+class Aringer2016(Stellib):
     """Aringer C+M+K giants Library
 
     References
@@ -1871,7 +1874,7 @@ class Aringer(Stellib):
 
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
-        self.name = "Aringer"
+        self.name = "Aringer2016"
         self.source = config["aringer"]
         self._load_()
 
@@ -1883,7 +1886,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
         ----------
@@ -1903,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),
@@ -1943,3 +1944,181 @@ def Z(self):
     @property
     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 = "BOSZ2024"
+        if filename is None:
+            self.source = config["bosz2024"]
+        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.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)
+
+    @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"]
+
+
+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
+
+        # 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
+
+        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.760 + dlogg),
+            (4.176 - dlogT, 1.750 - 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.760 + dlogg),
+            (4.176 - dlogT, 4.760 + 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"]