sunpy
Too few energy bins breaks f_vth/thermal_emission?
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
/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
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