thermal_emission crashes with IndexError on energy grids with limited bins - breaks Fitter for low-energy spectra

[cypherdavy]

Describe the bug

The thermal_emission function raises an IndexError when called with certain valid energy grids that contain relatively few bins. This also affects the Fitter workflow because the thermal model cannot be evaluated successfully.

Error

IndexError: index 84 is out of bounds for axis 0 with size 84

Minimal Reproduction

from sunkit_spex.models import thermal
import numpy as np
import astropy.units as u

engs = np.arange(1.01, 5.0, 0.01) << u.keV
temperature = 10e6 << u.K
emission_measure = 1e46 << u.cm**(-3)

thermal.thermal_emission(
    engs,
    temperature,
    emission_measure,
)

Expected Behavior

thermal_emission() should either successfully evaluate the spectrum or return a clear validation error if the supplied energy grid is unsupported.

Actual Behavior

The function raises:

IndexError: index 84 is out of bounds for axis 0 with size 84

which indicates an out-of-bounds array access during evaluation.

Additional Information

This appears related to the previously reported Issue #99, which notes that running the thermal model with certain energy grids can trigger an index error.

poc_bug.py

To Reproduce

To Reproduce

Reproduction 1: thermal_emission() raises IndexError

#!/usr/bin/env python3

import numpy as np
import astropy.units as u

from sunkit_spex.models import thermal

# Energy grid
engs = np.arange(1.01, 5.0, 0.01) << u.keV

temperature = 10e6 << u.K
emission_measure = 1e46 << u.cm**(-3)

# Expected: thermal emission spectrum
# Actual: IndexError
thermal_model = thermal.thermal_emission(
    engs,
    temperature,
    emission_measure,
)

print(thermal_model)

Expected Behavior

The function should return a thermal emission spectrum (or return a clear validation error if the input is unsupported).

Actual Behavior

IndexError: index 84 is out of bounds for axis 0 with size 84

---

Reproduction 2: Through the Fitter interface

import numpy as np

from sunkit_spex.legacy.fitting.fitter import Fitter

energy_bins = np.linspace(0.5, 10, 20)
counts = np.exp(-energy_bins[:-1] / 2.0) * 1000
counts_err = np.sqrt(counts)
srm = np.eye(len(energy_bins) - 1)

spec_data = {
    "count_channel_bins": energy_bins.reshape(1, -1),
    "counts": counts.reshape(1, -1),
    "count_error": counts_err.reshape(1, -1),
    "effective_exposure": 10.0,
    "srm": srm,
}

spec = Fitter(spec_data)
spec.loglikelihood = "poisson"
spec.model = "C*f_vth"

# Restrict fitting to energies above the model minimum
spec.energy_fitting_range = [1.83, 15.0]

spec.fit()

Actual Behavior

ValueError:
Lower bound of the input energy must be within the range
1.0002920302956426--200.15819869050395 keV.

Notes

• The first example reproduces the IndexError in thermal_emission().
• The second example demonstrates a validation error encountered through the Fitter interface when fitting a spectrum containing lower-energy bins.
• These may be related, but I have not confirmed whether they share the same root cause.

Screenshots

System Details

• sunpy version: 6.1.9
• astropy version: 8.0.1
• relevant package version(s):
  • sunkit-spex: 0.5.1.dev9+gdce9c42c4 (development version from GitHub)
  • numpy: 2.5.1
  • scipy: 1.18.0
  • matplotlib: 3.11.1
  • emcee: 3.1.6
  • corner: 2.3.0
  • numdifftools: installed
• Python version: 3.13.2
• OS information: Kali Linux 2025.2 (Kernel 6.11.2-amd64)

Installation method

No response

[settwi]

10 eV bins are quite fine in the first example. I agree errors relating to indexing should be more descriptive. I think this issue might be due to how some of the array slicing is happening in the thermal model interpolations. @DanRyanIrish might be the one to fix that.

The second example is the expected behavior and an unrelated issue. That is because the CHIANTI grid is only defined above 1.0003 keV, as the error message says. @ianan has worked on lowering the defined grid.

[cypherdavy]

Thanks for taking a look. That makes sense.

I agree that the IndexError appears to be separate from the CHIANTI energy-range validation issue. I'll treat the Fitter behavior as expected and keep this issue focused on the out-of-bounds indexing in thermal_emission().

If it helps, I'm happy to test any proposed patch or provide additional debugging information (full traceback, intermediate array shapes, or try different energy grids) to help verify the fix.