Too few energy bins breaks f_vth/thermal_emission?

[ianan]

Describe the bug

From some quick testing, I need more than 83 energy bins for the thermal model to work.

Should be able to supply a single energy bin and the model should still work.

The following code doesn't work. If change upper energy range to say mini, maxi, step = 1.6, 5.0, 0.04 and hence >83 bins it then works.

To Reproduce

import numpy as np
from sunxspex.sunxspex_fitting.photon_models_for_fitting import f_vth

mini, maxi, step = 1.6, 4.6, 0.04
chan_bins = np.stack((np.arange(mini,maxi, step), np.arange(mini+step,maxi+step, step)), axis=-1)
phmod=f_vth(10,1e3,energies=chan_bins)

What happened?

---

IndexError Traceback (most recent call last)
/var/folders/zp/m_v8v84d2ws3m21zq90wp2t40000gn/T/ipykernel_70315/4007065872.py in
6 chan_bins = np.stack((np.arange(mini,maxi, step), np.arange(mini+step,maxi+step, step)), axis=-1)
7 print(chan_bins.shape)
----> 8 phmod=f_vth(10,1e3,energies=chan_bins)

/opt/anaconda3/lib/python3.9/site-packages/sunxspex-0.1.dev131+gebc7f49-py3.9.egg/sunxspex/sunxspex_fitting/photon_models_for_fitting.py in f_vth(temperature, emission_measure46, energies)
49 temperature = temperature1e6 << u.K
50 emission_measure = emission_measure461e46 << u.cm**(-3)
---> 51 return thermal_emission(energies, temperature, emission_measure).value
52
53

/opt/anaconda3/lib/python3.9/site-packages/astropy/units/decorators.py in wrapper(*func_args, **func_kwargs)
302 # Call the original function with any equivalencies in force.
303 with add_enabled_equivalencies(self.equivalencies):
--> 304 return_ = wrapped_function(*func_args, **func_kwargs)
305
306 # Return

/opt/anaconda3/lib/python3.9/site-packages/sunxspex-0.1.dev131+gebc7f49-py3.9.egg/sunxspex/thermal.py in thermal_emission(energy_edges, temperature, emission_measure, abundance_type, relative_abundances, observer_distance)
219 # Calculate fluxes.
220 continuum_flux = _continuum_emission(energy_edges_keV, temperature_K, abundances)
--> 221 line_flux = _line_emission(energy_edges_keV, temperature_K, abundances)
222 flux = ((continuum_flux + line_flux) * emission_measure /
223 (4 * np.pi * observer_distance**2))

/opt/anaconda3/lib/python3.9/site-packages/sunxspex-0.1.dev131+gebc7f49-py3.9.egg/sunxspex/thermal.py in _line_emission(energy_edges_keV, temperature_K, abundances)
436 # the dimensionality of the line_intensities array.
437 line_peaks_keV = LINE_GRID["line peaks keV"][energy_roi_indices]
--> 438 split_line_intensities, line_spectrum_bins = _weight_emission_bins_to_line_centroid(
439 line_peaks_keV, energy_edges_keV, line_intensities)
440

/opt/anaconda3/lib/python3.9/site-packages/sunxspex-0.1.dev131+gebc7f49-py3.9.egg/sunxspex/thermal.py in _weight_emission_bins_to_line_centroid(line_peaks_keV, energy_edges_keV, line_intensities)
633
634 if len(pos_deviation_indices) > 0:
--> 635 pos_line_intensities, pos_neighbor_intensities, pos_neighbor_iline = _weight_emission_bins(
636 line_deviations_keV, pos_deviation_indices,
637 energy_center_diffs, line_intensities, iline, negative_deviations=False)

/opt/anaconda3/lib/python3.9/site-packages/sunxspex-0.1.dev131+gebc7f49-py3.9.egg/sunxspex/thermal.py in _weight_emission_bins(line_deviations_keV, deviation_indices, energy_center_diffs, line_intensities, iline, negative_deviations)
675 # fraction of the distance between the spectrum bin center and the
676 # center of the nearest neighboring bin.
--> 677 wghts = np.absolute(line_deviations_keV[deviation_indices]) / energy_center_diffs[iline[deviation_indices+a]]
678 # Tile/replicate wghts through the other dimension of line_intensities.
679 wghts = np.tile(wghts, tuple([line_intensities.shape[0]] + [1] * wghts.ndim))

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

Expected behavior

Return photon flux model for parameters supplied.

Screenshots

No response

System Details

General
#######
OS: Mac OS 10.16
Arch: 64bit, (i386)
sunpy: 4.0.0
Installation path: /opt/anaconda3/lib/python3.9/site-packages/sunpy-4.0.0.dist-info

Required Dependencies
#####################
astropy: 5.0.4
numpy: 1.20.3
packaging: 21.0
parfive: 1.5.1

Optional Dependencies
#####################
asdf: 2.8.3
asdf-astropy: 0.2.1
beautifulsoup4: 4.10.0
cdflib: 0.3.20
dask: 2021.10.0
drms: 0.6.2
glymur: 0.9.7.post1
h5netcdf: 0.13.0
h5py: 3.2.1
matplotlib: 3.4.3
mpl-animators: 1.0.0
pandas: 1.3.4
python-dateutil: 2.8.2
reproject: 0.8
scikit-image: 0.18.3
scipy: 1.7.1
sqlalchemy: 1.4.22
tqdm: 4.62.3
zeep: 4.1.0

Installation method

git checkout

[ianan]

There is an issue here, but f_vth should only be called by the fitting routines. When manually wanting to do emission models call them directly, i.e. from sunxspex.thermal import thermal_emission