Updating non-LTE, LTE, and AGB stellar libraries

[karllark]

The Tlusty non-LTE models have been updated in 2025 to have an enhanced line list and better coverage in the IR. This stellar library supersedes the existing tlusty library.

The LTE model grid has been updated in the BOSZ 2024 work. This updated library supercedes the existing Kuruz library.

The AGB model grid has been updated to have the same updated spectral resolution and wavelength range as the other two grids. This uses the Aringer 2016 grid. Much of the code for this was adapted from Martha Boyer's code that she shared with me. These models do not cover the full wavelength range, so a black body with the stellar Teff was used to extend the coverage. This was done at the shorter wavelengths below 3600 A and longer wavelengths beyond 25 micron.

A common resolution and wavelength range is required if used of multiple libraries at the same time is desired. While there is code to handle non-similar wavelength grids, it has failures at the edges at this point.

Code is included that generates the libraries as well as plotting the coverage and boundaries.

Spectra extend from 500 A to 32 micron. The spectral resolution has been set to 2000.

As part of this work, it was discovered that the original tlusty grid had significantly different logT coverage for different metallicities. Specifically, at Z = 2e-5, 2e-4, 4e-4, and 6.6e-4 the lowest logT was 4.44. This causes an issue with the interpolation routines that assume the coverage in logT, logg is the same at all metallicities. There is no clear fix for the old tlusty grid and it should be no longer be used, the new tlusty2025 grid is better in a number of ways.

The new Tlusty2025 grid shows related issue in that at Z=0.000544 the highest logT values are 4.8 and at Z = 0.000561, the lowest logT values are 4.44. These two metallicities are close together and code has been added to the Tlusty2025 stellib to set the Z values for these to values to be the average of the two (i.e. 0.000553) so that full coverage is available at this Z value and the interpolation routines will work correctly.

Closes #180.
Closes #669.
Closes #858.

[codecov]

Codecov Report

❌ Patch coverage is 40.57971% with 41 lines in your changes missing coverage. Please review.
✅ Project coverage is 39.68%. Comparing base (bffcdfc) to head (eb624a9).
⚠️ Report is 23 commits behind head on master.

Files with missing lines	Patch %	Lines
beast/physicsmodel/stars/stellib.py	40.57%	41 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #860      +/-   ##
==========================================
+ Coverage   39.44%   39.68%   +0.24%     
==========================================
  Files         110      109       -1     
  Lines       10755    10742      -13     
==========================================
+ Hits         4242     4263      +21     
+ Misses       6513     6479      -34     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[karllark]

Here is an example of the incomplete coverage in the old tlusty grid for some metallicities.

[karllark]

Here is the full coverage - like it should be for all metallicities.

[karllark]

doc_test failure is unrelated to this PR.

[karllark]

Coverage of the three grids.

[karllark]

Coverage of the original Tlusty and Kurucz grids.

[galaxyumi]

Thanks for this updating! I left some specific comments for the changed files, but I do have a general question. Do we account for alpha-abundance mapping at all?

[galaxyumi]

With new models, do we still stick to this old Zsun value? For example, the PARSEC models use Zsun=0.0152.

[karllark]

The solar_z value for each model is different. I've gone back to the papers to figure out the correct value for each grid. This is given with references in each script that generates the grids. We want the absolute metallicity, not relative to solar as each grid uses a different solar value - as the solar value changes with time. Exciting.

I will update this script to compare models from any two stellar libraries using the grids themselves. This is needed to ensure that the scaling is correct for all the models. And provides a great way to compare models from two grids for the same stellar parameters. Will not resolve this comment till the expanded script is ready.

[galaxyumi]

This line contains a hardcoded absolute path pointing to your specific directory. To ensure portability across different environments and for other users, please replace this with a relative path or leverage the standard beast data directory configuration paths. Also the corresponding documentation should be provided.

[karllark]

See comment above.

[galaxyumi]

Same comment as above.

[karllark]

Ditto. ;-)

[galaxyumi]

Do the boundaries of new lower-resolution bins align with the pixel edges of the high-resolution input spectrum? I guess not. If an input spans the boundary (e.g., 40% of its width is in bin A and 60% is in bin B), this code assigs the entire flux of that pixel to one bin and zero to the other based on where its center point of edge falls. It might not be a significant fraction of flux, but something to think about.

[galaxyumi]

I tried to come up with a block of code that could replace the rest of this function if this is doing what you intended, which has the following benefit:

1. By sampling the cumulative distribution function at exactly twave, pixels that are split across a boundary have their integrated energy split perfectly proportionally down to decimals. No flux is lost or duplicated.
2. By explicitly calculating and carrying in_widths (Δλ_input), it will correctly downgrade models even if the source library switches sampling strategies between different regimes (e.g., if a library samples line cores more densely than the continuum).

# Define edge boundaries for the input high-res pixels
# If the input grid contains center coordinates, estimate its bin edges
in_edges = np.zeros(len(wave) + 1)
in_edges[1:-1] = 0.5 * (wave[0:-1] + wave[1:])
in_edges[0] = wave[0] - (wave[1] - wave[0]) * 0.5
in_edges[-1] = wave[-1] + (wave[-1] - wave[-2]) * 0.5

# Calculate total energy (Flux * delta_lambda) per input bin
in_widths = np.diff(in_edges)
energy = flux * in_widths

# Compute the cumulative integrated energy distribution
cum_energy = np.zeros(len(in_edges))
cum_energy[1:] = np.cumsum(energy)

# Interpolate cumulative energy onto the new output grid edges
# Use linear interpolation because cum_energy is stepwise linear (if flux constant in bins)
# or piecewise quadratic (if trapezoidal), but linear on step integrals keeps exact step conservation.
energy_interp = interp1d(
    in_edges, cum_energy, kind="linear", bounds_error=False, fill_value=0.0
)

new_cum_energy = energy_interp(twave)

# Total energy integrated inside each new bin window
new_bin_energy = np.diff(new_cum_energy)

# Convert energy back to flux density (divide by new bin widths)
new_widths = np.diff(twave)
full_flux = new_bin_energy / new_widths

# Track contributing pixels & catch out-of-bounds
full_npts = np.zeros(npts, dtype=int)
for k in range(npts):
    # Kept for compatibility with BEAST's diagnostic tracking downstream
    full_npts[k] = np.sum((wave >= full_wave_min[k]) & (wave < full_wave_max[k]))

# If the output grid expands outside the range of the input library,
# those bins will have accumulated 0 energy and have zero widths or false zeros.
# We flag them as non-finite or zero out-of-bounds where the input didn't reach.
out_of_bounds = (twave[:-1] < in_edges[0]) | (
    twave[1:] > in_edges[-1]
)
full_flux[out_of_bounds] = 0.0
full_npts[out_of_bounds] = 0

return (full_wave, full_flux, full_npts)

[galaxyumi]

Clean up unused example lists.