From cd3a615b229ce346c16fd2583d127e08d997b4a6 Mon Sep 17 00:00:00 2001 From: Kristopher Cooper Date: Mon, 27 Jul 2026 11:49:16 -0500 Subject: [PATCH 1/9] Trying to fix scaling of models --- sunkit_spex/models/physical/nonthermal.py | 200 ++++------------------ sunkit_spex/models/physical/thermal.py | 3 + sunkit_spex/models/scaling.py | 5 +- 3 files changed, 42 insertions(+), 166 deletions(-) diff --git a/sunkit_spex/models/physical/nonthermal.py b/sunkit_spex/models/physical/nonthermal.py index 1e95e29a..81c907bf 100644 --- a/sunkit_spex/models/physical/nonthermal.py +++ b/sunkit_spex/models/physical/nonthermal.py @@ -5,6 +5,7 @@ import astropy.units as u from astropy.modeling import FittableModel, Parameter +from sunkit_spex.models.scaling import norm_thick_target_eflux_units from sunkit_spex.legacy import constants as const from sunkit_spex.legacy.integrate import gauss_legendre @@ -88,7 +89,7 @@ class ThickTarget(FittableModel): ) total_eflux = Parameter( - name="total_eflux", default=1.5, unit=u.electron * u.s**-1, description="Total electron flux", fixed=True + name="total_eflux", default=1.5e35, unit=(u.electron * u.s**-1), description="Total electron flux", fixed=True ) _input_units_allow_dimensionless = True @@ -119,27 +120,23 @@ def __init__( def evaluate(self, energy_edges, p, break_energy, q, low_e_cutoff, high_e_cutoff, total_eflux): energy_centers = energy_edges[:-1] + 0.5 * np.diff(energy_edges) - if ( - hasattr(break_energy, "unit") - or hasattr(energy_centers, "unit") - or hasattr(low_e_cutoff, "unit") - or hasattr(high_e_cutoff, "unit") - or hasattr(total_eflux, "unit") - ): - flux = thick_fn( - energy_centers.value, - p, - break_energy.value, - q, - low_e_cutoff.value, - high_e_cutoff.value, - total_eflux.value, - self.integrator, - ) - else: - flux = thick_fn( - energy_centers, p, break_energy, q, low_e_cutoff, high_e_cutoff, total_eflux, self.integrator - ) + energy_centers <<= u.keV + break_energy <<= self.break_energy.unit + low_e_cutoff <<= self.low_e_cutoff.unit + high_e_cutoff <<= self.high_e_cutoff.unit + print(total_eflux) + total_eflux <<= norm_thick_target_eflux_units + print(total_eflux) + + flux = bremsstrahlung_thick_target( + energy_centers.value, + p, + break_energy.value, + q, + low_e_cutoff.value, + high_e_cutoff.value, + self.integrator, + ) * total_eflux.value return flux @@ -255,28 +252,22 @@ def __init__( def evaluate(self, energy_edges, p, break_energy, q, low_e_cutoff, high_e_cutoff, total_eflux): energy_centers = energy_edges[:-1] + 0.5 * np.diff(energy_edges) - - if ( - hasattr(break_energy, "unit") - or hasattr(energy_centers, "unit") - or hasattr(low_e_cutoff, "unit") - or hasattr(high_e_cutoff, "unit") - or hasattr(total_eflux, "unit") - ): - flux = thin_fn( - energy_centers.value, - p, - break_energy.value, - q, - low_e_cutoff.value, - high_e_cutoff.value, - total_eflux.value, - self.integrator, - ) - else: - flux = thin_fn( - energy_centers, p, break_energy, q, low_e_cutoff, high_e_cutoff, total_eflux, self.integrator - ) + + energy_centers <<= u.keV + break_energy <<= self.break_energy.unit + low_e_cutoff <<= self.low_e_cutoff.unit + high_e_cutoff <<= self.high_e_cutoff.unit + total_eflux <<= self.total_eflux.unit + + flux = bremsstrahlung_thin_target( + energy_centers.value, + p, + break_energy.value, + q, + low_e_cutoff.value, + high_e_cutoff.value, + self.integrator, + ) * total_eflux.value * 1e55 return flux @@ -299,127 +290,6 @@ def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): } -def thick_fn(energy_centers, p, break_energy, q, low_e_cutoff, high_e_cutoff, total_eflux, integrator): - """Calculates the thick-target bremsstrahlung radiation of a dual power-law electron distribution. - - [1] Brown, Solar Physics 18, 489 (1971) (https://link.springer.com/article/10.1007/BF00149070) - [2] https://hesperia.gsfc.nasa.gov/ssw/packages/xray/doc/brm_thick_doc.pdf - [3] https://hesperia.gsfc.nasa.gov/ssw/packages/xray/idl/brm2/brm2_thicktarget.pro - - Parameters - ---------- - - energy_edges : 1d array - Edges of energy bins in units of keV. - - total_eflux : int or float - Total integrated electron flux, in units of 10^35 e^- s^-1. - Need to take care here as the model returns units of cm-2 sec-1 as the scaling factor of 1e35 is hidden. - So actual units are 1.0d35 e^- s^-1. - - p : int or float - Power-law index of the electron distribution below the break. - - break_energy : int or float - Break energy of power law. - - q : int or float - Power-law index of the electron distribution above the break. - - low_e_cutoff : int or float - Low-energy cut-off of the electron distribution in units of keV. - - high_e_cutoff : int or float - High-energy cut-off of the electron distribution in units of keV. - - - Returns - ------- - A 1d array of thick-target bremsstrahlung radiation in units - of ph s^-1 keV^-1. - """ - - # hack = np.round([p, break_energy, q, low_e_cutoff, high_e_cutoff, total_eflux], 15) - # p, break_energy, q, low_e_cutoff, high_e_cutoff, total_eflux = hack[0], hack[1], hack[2], hack[3], hack[4], hack[5] - - # energies = np.mean(energies, axis=1) # since energy bins are given, use midpoints though - - # we want a single power law electron distribution, - # so set break_energy == high_e_cutoff at a high value. - # we don't care about q at E > break_energy. - # high_break = energies.max() * 10 - - output = bremsstrahlung_thick_target(energy_centers, p, break_energy, q, low_e_cutoff, high_e_cutoff, integrator) - - output[np.isnan(output)] = 0 - output[~np.isfinite(output)] = 0 - - # convert to 1e35 e-/s - return output * total_eflux * 1e35 - - -# def thin_fn(total_eflux, index, e_c, energies=None): -def thin_fn(energy_centers, p, break_energy, q, low_e_cutoff, high_e_cutoff, total_eflux, integrator): - """Calculates the thin-target bremsstrahlung radiation of a dual power-law electron distribution. - - [1] Brown, Solar Physics 18, 489 (1971) (https://link.springer.com/article/10.1007/BF00149070) - [2] https://hesperia.gsfc.nasa.gov/ssw/packages/xray/doc/brm_thick_doc.pdf - [3] https://hesperia.gsfc.nasa.gov/ssw/packages/xray/idl/brm2/brm2_thicktarget.pro - - Parameters - ---------- - energy_edges : 1d array - Edges of energy bins in units of keV. - - total_eflux : int or float - normalization factor in units of 1.0d55 cm-2 sec-1, - i.e. plasma density * volume of source * integrated nonthermal electron flux density - Need to take care here as the model returns units of cm-2 sec-1 as the scaling factor of 1e55 is hidden. - So actual units are 1.0d55 cm-2 sec-1. - - p : int or float - Power-law index of the electron distribution below the break. - - break_energy : int or float - Break energy of power law. - - q : int or float - Power-law index of the electron distribution above the break. - - low_e_cutoff : int or float - Low-energy cut-off of the electron distribution in units of keV. - - high_e_cutoff : int or float - High-energy cut-off of the electron distribution in units of keV. - - - - Returns - ------- - A 1d array of thin-target bremsstrahlung radiation in units - of ph s^-1 keV^-1. - """ - - # hack = np.round([total_eflux, index, e_c], 15) - # total_eflux, index, e_c = hack[0], hack[1], hack[2] - - # energies = np.mean(energies, axis=1) # since energy bins are given, use midpoints though - # energies = energy_centers - # we want a single power law electron distribution, - # so set break_energy == high_e_cutoff at a high value. - # we don't care about q at E > break_energy. - # high_break = energies.max() * 10 - output = bremsstrahlung_thin_target( - energy_centers, p, break_energy, q, low_e_cutoff, high_e_cutoff, total_eflux, integrator - ) - - output[np.isnan(output)] = 0 - output[~np.isfinite(output)] = 0 - - # convert to 1e35 e-/s - return output * total_eflux * 1e55 - - class BrokenPowerLawElectronDistribution: """ A broken or double power law electron flux distribution and integral. diff --git a/sunkit_spex/models/physical/thermal.py b/sunkit_spex/models/physical/thermal.py index 73d09897..f3c8073e 100644 --- a/sunkit_spex/models/physical/thermal.py +++ b/sunkit_spex/models/physical/thermal.py @@ -237,6 +237,9 @@ def evaluate( ca, fe, ): + energy_edges <<= u.keV + temperature <<= self.temperature.unit + emission_measure <<= self.emission_measure.unit line_flux = self.line.evaluate( energy_edges, temperature, diff --git a/sunkit_spex/models/scaling.py b/sunkit_spex/models/scaling.py index 549b92ff..29d1ff09 100644 --- a/sunkit_spex/models/scaling.py +++ b/sunkit_spex/models/scaling.py @@ -4,9 +4,12 @@ from astropy.modeling import FittableModel, Parameter from astropy.units import Quantity -__all__ = ["Constant", "InverseSquareFluxScaling"] +__all__ = ["Constant", "InverseSquareFluxScaling", "scaled_thick_target_eflux_units"] + +norm_thick_target_eflux_units = u.def_unit("scaled_eflux_units", 1e-35 * (u.electron * u.s**-1)) + class InverseSquareFluxScaling(FittableModel): """ InverseSqaureFluxScaling model converts luminosity output of physical models to a distance scaled flux. From a949b0febcfe13b68ff5febf84b806ef39ed9bd6 Mon Sep 17 00:00:00 2001 From: Kristopher Cooper Date: Mon, 27 Jul 2026 16:42:37 -0500 Subject: [PATCH 2/9] Units and norm scaling is fixed in the thermal, thick, thin models --- sunkit_spex/models/physical/nonthermal.py | 31 ++++++++++--------- .../models/physical/tests/test_nonthermal.py | 25 +++++++++++++++ .../models/physical/tests/test_thermal.py | 27 ++++++++++++---- sunkit_spex/models/physical/thermal.py | 29 +++++++++-------- sunkit_spex/models/scaling.py | 10 +++--- 5 files changed, 82 insertions(+), 40 deletions(-) diff --git a/sunkit_spex/models/physical/nonthermal.py b/sunkit_spex/models/physical/nonthermal.py index 81c907bf..fde58137 100644 --- a/sunkit_spex/models/physical/nonthermal.py +++ b/sunkit_spex/models/physical/nonthermal.py @@ -5,7 +5,7 @@ import astropy.units as u from astropy.modeling import FittableModel, Parameter -from sunkit_spex.models.scaling import norm_thick_target_eflux_units +from sunkit_spex.models.scaling import norm_thick_target_eflux_units, norm_thin_target_eflux_units from sunkit_spex.legacy import constants as const from sunkit_spex.legacy.integrate import gauss_legendre @@ -89,7 +89,7 @@ class ThickTarget(FittableModel): ) total_eflux = Parameter( - name="total_eflux", default=1.5e35, unit=(u.electron * u.s**-1), description="Total electron flux", fixed=True + name="total_eflux", default=1.5, unit=norm_thick_target_eflux_units, description="Total electron flux", fixed=False ) _input_units_allow_dimensionless = True @@ -107,6 +107,8 @@ def __init__( ): self.integrator = integrator + total_eflux <<= norm_thick_target_eflux_units + super().__init__( p=p, break_energy=break_energy, @@ -124,9 +126,7 @@ def evaluate(self, energy_edges, p, break_energy, q, low_e_cutoff, high_e_cutoff break_energy <<= self.break_energy.unit low_e_cutoff <<= self.low_e_cutoff.unit high_e_cutoff <<= self.high_e_cutoff.unit - print(total_eflux) total_eflux <<= norm_thick_target_eflux_units - print(total_eflux) flux = bremsstrahlung_thick_target( energy_centers.value, @@ -136,14 +136,14 @@ def evaluate(self, energy_edges, p, break_energy, q, low_e_cutoff, high_e_cutoff low_e_cutoff.value, high_e_cutoff.value, self.integrator, - ) * total_eflux.value + ) * total_eflux.decompose().value - return flux + return flux * self.return_units[self.outputs[0]] @property def input_units(self): # The units for the 'energy_edges' variable should be an energy (default keV) - return {self.inputs[0]: u.keV} + return {self.inputs[0]: u.keV} @property def return_units(self): @@ -154,10 +154,9 @@ def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): "break_energy": u.keV, "low_e_cutoff": u.keV, "high_e_cutoff": u.keV, - "total_eflux": u.electron * u.s**-1, + "total_eflux": norm_thick_target_eflux_units, } - class ThinTarget(FittableModel): r"""Calculates the thin-target bremsstrahlung radiation of a dual power-law electron distribution. @@ -222,7 +221,7 @@ class ThinTarget(FittableModel): ) total_eflux = Parameter( - name="total_eflux", default=1.5, unit=u.s**-1 * u.cm**-2, description="Total electron flux", fixed=True + name="total_eflux", default=1.5, unit=norm_thin_target_eflux_units, description="Total electron flux", fixed=True ) _input_units_allow_dimensionless = True @@ -240,6 +239,8 @@ def __init__( ): self.integrator = integrator + total_eflux <<= norm_thin_target_eflux_units + super().__init__( p=p, break_energy=break_energy, @@ -257,7 +258,7 @@ def evaluate(self, energy_edges, p, break_energy, q, low_e_cutoff, high_e_cutoff break_energy <<= self.break_energy.unit low_e_cutoff <<= self.low_e_cutoff.unit high_e_cutoff <<= self.high_e_cutoff.unit - total_eflux <<= self.total_eflux.unit + total_eflux <<= norm_thin_target_eflux_units flux = bremsstrahlung_thin_target( energy_centers.value, @@ -266,10 +267,10 @@ def evaluate(self, energy_edges, p, break_energy, q, low_e_cutoff, high_e_cutoff q, low_e_cutoff.value, high_e_cutoff.value, - self.integrator, - ) * total_eflux.value * 1e55 + integrator=self.integrator, + ) * total_eflux.to((u.electron * u.cm ** (-2) * u.s**-1)).value - return flux + return flux * self.return_units[self.outputs[0]] @property def input_units(self): @@ -285,7 +286,7 @@ def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): "break_energy": u.keV, "low_e_cutoff": u.keV, "high_e_cutoff": u.keV, - "total_eflux": u.s**-1 * u.cm**-2, + "total_eflux": norm_thin_target_eflux_units, # "total_eflux": u.electron * u.s**-1, } diff --git a/sunkit_spex/models/physical/tests/test_nonthermal.py b/sunkit_spex/models/physical/tests/test_nonthermal.py index ac5419f8..4f175bef 100644 --- a/sunkit_spex/models/physical/tests/test_nonthermal.py +++ b/sunkit_spex/models/physical/tests/test_nonthermal.py @@ -4,6 +4,7 @@ import astropy.units as u from sunkit_spex.models.physical import nonthermal +from sunkit_spex.models.scaling import norm_thick_target_eflux_units, norm_thin_target_eflux_units SSW_INTENSITY_UNIT = u.ph / u.cm**2 / u.s / u.keV @@ -144,3 +145,27 @@ def test_thin_target_against_ssw(ssw): output = model(energy_edges) expected_value = expected.to_value(output.unit) np.testing.assert_allclose(output.value, expected_value, rtol=0.035) + +def test_thick_target_flux_scaling(): + """Test thick target flux units being scaled.""" + energy_edges = np.arange(2, 15, 0.1) << u.keV + for _flux in np.arange(1, 10, 0.5): + eflux = (_flux * 1e35) << (u.electron * u.s**-1) + s_eflux = _flux << norm_thick_target_eflux_units + model = nonthermal.ThickTarget(total_eflux=eflux) + s_model = nonthermal.ThickTarget(total_eflux=s_eflux) + np.testing.assert_allclose(model(energy_edges).value, s_model(energy_edges).value) + np.testing.assert_allclose(model.evaluate(energy_edges, *model.parameters).value, + s_model.evaluate(energy_edges, *s_model.parameters).value) + +def test_thin_target_flux_scaling(): + """Test thin target flux units being scaled.""" + energy_edges = np.arange(2, 15, 0.1) << u.keV + for _flux in np.arange(1, 10, 0.5): + eflux = (_flux * 1e55) << (u.electron * u.cm**-2 * u.s**-1) + s_eflux = _flux << norm_thin_target_eflux_units + model = nonthermal.ThinTarget(total_eflux=eflux) + s_model = nonthermal.ThinTarget(total_eflux=s_eflux) + np.testing.assert_allclose(model(energy_edges).value, s_model(energy_edges).value) + np.testing.assert_allclose(model.evaluate(energy_edges, *model.parameters).value, + s_model.evaluate(energy_edges, *s_model.parameters).value) diff --git a/sunkit_spex/models/physical/tests/test_thermal.py b/sunkit_spex/models/physical/tests/test_thermal.py index 89fb2c53..a0d837ce 100644 --- a/sunkit_spex/models/physical/tests/test_thermal.py +++ b/sunkit_spex/models/physical/tests/test_thermal.py @@ -5,6 +5,7 @@ import astropy.units as u +from sunkit_spex.models.scaling import norm_thermal_emission_measure_units from sunkit_spex.models.physical import thermal # Manually load file that was used to compile expected flux values. @@ -50,7 +51,7 @@ def fvth_simple(): """ energy_edges = np.arange(3, 28.5, 0.5) * u.keV temperature = 6 * u.MK - emission_measure = 1e-5 / u.cm**3 + emission_measure = 1e44 / u.cm**3 abundance_type = DEFAULT_ABUNDANCE_TYPE observer_distance = (1 * u.AU).to(u.cm) # fmt: off @@ -111,7 +112,7 @@ def chianti_kev_cont_simple(): """ energy_edges = np.arange(3, 28.5, 0.5) * u.keV temperature = 6 * u.MK - emission_measure = 1e-5 / u.cm**3 + emission_measure = 1e44 / u.cm**3 abundance_type = DEFAULT_ABUNDANCE_TYPE observer_distance = (1 * u.AU).to(u.cm) # fmt: off @@ -172,7 +173,7 @@ def chianti_kev_lines_simple(): """ energy_edges = np.arange(3, 28.5, 0.5) * u.keV temperature = 6 * u.MK - emission_measure = 1e-5 / u.cm**3 + emission_measure = 1e44 / u.cm**3 abundance_type = DEFAULT_ABUNDANCE_TYPE observer_distance = (1 * u.AU).to(u.cm) # fmt: off @@ -234,7 +235,7 @@ def fvth_Fe2(): """ energy_edges = np.arange(3, 28.5, 0.5) * u.keV temperature = 6 * u.MK - emission_measure = 1e-5 / u.cm**3 + emission_measure = 1e44 / u.cm**3 abundance_type = DEFAULT_ABUNDANCE_TYPE observer_distance = (1 * u.AU).to(u.cm) # fmt: off @@ -297,7 +298,7 @@ def chianti_kev_cont_Fe2(): """ energy_edges = np.arange(3, 28.5, 0.5) * u.keV temperature = 6 * u.MK - emission_measure = 1e-5 / u.cm**3 + emission_measure = 1e44 / u.cm**3 abundance_type = DEFAULT_ABUNDANCE_TYPE observer_distance = (1 * u.AU).to(u.cm) # fmt: off @@ -360,7 +361,7 @@ def chianti_kev_lines_Fe2(): """ energy_edges = np.arange(3, 28.5, 0.5) * u.keV temperature = 6 * u.MK - emission_measure = 1e-5 / u.cm**3 + emission_measure = 1e44 / u.cm**3 abundance_type = DEFAULT_ABUNDANCE_TYPE observer_distance = (1 * u.AU).to(u.cm) # fmt: off @@ -513,3 +514,17 @@ def test_abundances_should_not_change(): after_models = thermal.DEFAULT_ABUNDANCES[thermal.DEFAULT_ABUNDANCE_TYPE].data assert np.allclose(after_models.data, orig.data) + +def test_thermal_emission_measure_scaling(): + """Test thermal emission measure units being scaled.""" + energy_edges = np.arange(2, 15, 0.1) << u.keV + for _em in np.arange(1, 10, 0.5): + em = (_em * 1e49) << (u.cm**-3) + s_em = _em << norm_thermal_emission_measure_units + model = thermal.ThermalEmission(emission_measure=em) + s_model = thermal.ThermalEmission(emission_measure=s_em) + np.testing.assert_allclose(model(energy_edges).value, s_model(energy_edges).value) + np.testing.assert_allclose(model.evaluate(energy_edges, *model.parameters).value, + s_model.evaluate(energy_edges, *s_model.parameters).value) + +test_thermal_emission_measure_scaling() \ No newline at end of file diff --git a/sunkit_spex/models/physical/thermal.py b/sunkit_spex/models/physical/thermal.py index f3c8073e..9047fd73 100644 --- a/sunkit_spex/models/physical/thermal.py +++ b/sunkit_spex/models/physical/thermal.py @@ -15,6 +15,7 @@ load_chianti_lines_lite, load_xray_abundances, ) +from sunkit_spex.models.scaling import norm_thermal_emission_measure_units # The default elemental abundance values correspond to coronal values DEFAULT_ABUNDANCE_TYPE = "sun_coronal_ext" @@ -144,7 +145,7 @@ class ThermalEmission(FittableModel): emission_measure = Parameter( name="emission_measure", default=1, - unit=(u.cm ** (-3)), + unit=norm_thermal_emission_measure_units, description="Emission measure of the observer", fixed=False, ) @@ -185,6 +186,8 @@ def __init__( if abundance_type != DEFAULT_ABUNDANCE_TYPE: mg, al, si, s, ar, ca, fe = _initialize_abundances(DEFAULT_ABUNDANCES[abundance_type]) + emission_measure <<= norm_thermal_emission_measure_units + self.line = LineEmission( temperature=temperature, emission_measure=emission_measure, @@ -278,7 +281,7 @@ def return_units(self): return {self.outputs[0]: u.ph / u.keV * u.s**-1} def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): - return {"temperature": u.MK, "emission_measure": (u.cm ** (-3))} + return {"temperature": u.MK, "emission_measure": norm_thermal_emission_measure_units} class ContinuumEmission(FittableModel): @@ -332,7 +335,7 @@ class ContinuumEmission(FittableModel): emission_measure = Parameter( name="emission_measure", default=1, - unit=(u.cm ** (-3)), + unit=norm_thermal_emission_measure_units, description="Emission measure of the observer", fixed=False, ) @@ -422,7 +425,7 @@ def return_units(self): return {self.outputs[0]: u.ph / u.keV * u.s**-1} def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): - return {"temperature": u.MK, "emission_measure": (u.cm ** (-3))} + return {"temperature": u.MK, "emission_measure": norm_thermal_emission_measure_units} class LineEmission(FittableModel): @@ -471,8 +474,8 @@ class LineEmission(FittableModel): emission_measure = Parameter( name="emission_measure", - default=1e50, - unit=(u.cm ** (-3)), + default=1, + unit=norm_thermal_emission_measure_units, description="Emission measure of the observer", fixed=False, ) @@ -511,7 +514,7 @@ def __init__( if abundance_type != DEFAULT_ABUNDANCE_TYPE: mg, al, si, s, ar, ca, fe = _initialize_abundances(DEFAULT_ABUNDANCES[abundance_type]) - + super().__init__( temperature=temperature, emission_measure=emission_measure, @@ -562,7 +565,7 @@ def return_units(self): return {self.outputs[0]: u.ph / u.keV * u.s**-1} def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): - return {"temperature": u.MK, "emission_measure": (u.cm ** (-3))} + return {"temperature": u.MK, "emission_measure": norm_thermal_emission_measure_units} def setup_continuum_parameters(filename=None): @@ -725,7 +728,7 @@ def continuum_emission( # Calculate flux. flux = _continuum_emission(energy_edges_keV, temperature_K, abundances) - flux *= emission_measure * 1e49 + flux *= emission_measure if temperature_K.isscalar and emission_measure.isscalar: flux = flux[0] @@ -769,8 +772,8 @@ def line_emission( abundances = _calculate_abundances(abundance_type, mg, al, si, s, ar, ca, fe) flux = _line_emission(energy_edges_keV, temperature_K, abundances) - - flux *= emission_measure * 1e49 + + flux *= emission_measure if temperature_K.isscalar and emission_measure.isscalar: flux = flux[0] @@ -1226,14 +1229,14 @@ def _sanitize_inputs(energy_edges, temperature, emission_measure): # If they were not already Quantities, the parameters get the default units. energy_edges <<= u.keV temperature <<= u.K - emission_measure <<= u.cm**-3 + emission_measure <<= norm_thermal_emission_measure_units energy_edges_keV = energy_edges.to(u.keV) temperature_K = temperature.to(u.K) if temperature.isscalar: temperature_K = np.array([temperature_K.value]) * temperature_K.unit - + emission_measure = emission_measure.to(u.cm**-3) if emission_measure.isscalar: emission_measure = np.array([emission_measure.value]) * emission_measure.unit diff --git a/sunkit_spex/models/scaling.py b/sunkit_spex/models/scaling.py index 29d1ff09..a66ae4c8 100644 --- a/sunkit_spex/models/scaling.py +++ b/sunkit_spex/models/scaling.py @@ -4,11 +4,13 @@ from astropy.modeling import FittableModel, Parameter from astropy.units import Quantity -__all__ = ["Constant", "InverseSquareFluxScaling", "scaled_thick_target_eflux_units"] +__all__ = ["Constant", "InverseSquareFluxScaling", "norm_thick_target_eflux_units", "norm_thin_target_eflux_units", "norm_thermal_emission_measure_units"] -norm_thick_target_eflux_units = u.def_unit("scaled_eflux_units", 1e-35 * (u.electron * u.s**-1)) +norm_thick_target_eflux_units = u.def_unit("scaled_thick_eflux_units", 1e35 * (u.electron * u.s**-1)) +norm_thin_target_eflux_units = u.def_unit("scaled_thin_eflux_units", 1e55 * (u.electron * u.cm ** (-2) * u.s**-1)) +norm_thermal_emission_measure_units = u.def_unit("scaled_em_units", 1e49 * (u.cm ** (-3))) class InverseSquareFluxScaling(FittableModel): """ @@ -21,9 +23,6 @@ class InverseSquareFluxScaling(FittableModel): observer_distance: Distance of the observer from the source. - - - Examples ======== .. plot:: @@ -104,7 +103,6 @@ class Constant(FittableModel): constant : A constant value which populates the output array - Examples ======== .. plot:: From e0972f93916c3d3c8441d3c3ec6cff0372fbe4c4 Mon Sep 17 00:00:00 2001 From: Kristopher Cooper Date: Tue, 28 Jul 2026 15:55:44 -0500 Subject: [PATCH 3/9] Found reason why the thick target model would return NaNs --- sunkit_spex/models/physical/nonthermal.py | 15 ++++++++------ .../models/physical/tests/test_nonthermal.py | 20 +++++++++++++++++++ sunkit_spex/models/physical/thermal.py | 8 ++++---- 3 files changed, 33 insertions(+), 10 deletions(-) diff --git a/sunkit_spex/models/physical/nonthermal.py b/sunkit_spex/models/physical/nonthermal.py index fde58137..99c83b55 100644 --- a/sunkit_spex/models/physical/nonthermal.py +++ b/sunkit_spex/models/physical/nonthermal.py @@ -4,6 +4,7 @@ import astropy.units as u from astropy.modeling import FittableModel, Parameter +from astropy.modeling.functional_models import FLOAT_EPSILON from sunkit_spex.models.scaling import norm_thick_target_eflux_units, norm_thin_target_eflux_units from sunkit_spex.legacy import constants as const @@ -30,6 +31,8 @@ __all__ = ["ThickTarget", "ThinTarget"] +FLOAT_EPSILON_FOR_POWER_LAW = 1+FLOAT_EPSILON*1e30 + class ThickTarget(FittableModel): r"""Calculates the thick-target bremsstrahlung radiation of a dual power-law electron distribution. @@ -74,22 +77,22 @@ class ThickTarget(FittableModel): n_inputs = 1 n_outputs = 1 - p = Parameter(name="p", default=2, description="Slope below break", fixed=False) + p = Parameter(name="p", default=2, description="Slope below break", fixed=False, bounds=(FLOAT_EPSILON_FOR_POWER_LAW, None)) - break_energy = Parameter(name="break_energy", default=100, unit=u.keV, description="Break Energy", fixed=False) + break_energy = Parameter(name="break_energy", default=100, unit=u.keV, description="Break Energy", fixed=True, bounds=(FLOAT_EPSILON, None)) - q = Parameter(name="q", default=5, min=0.01, description="Slope above break", fixed=True) + q = Parameter(name="q", default=5, description="Slope above break", fixed=True, bounds=(FLOAT_EPSILON_FOR_POWER_LAW, None)) low_e_cutoff = Parameter( - name="low_e_cutoff", default=7, unit=u.keV, description="Low energy electron cut off", fixed=False + name="low_e_cutoff", default=7, unit=u.keV, description="Low energy electron cut off", fixed=False, bounds=(FLOAT_EPSILON, None) ) high_e_cutoff = Parameter( - name="high_e_cutoff", default=1500, unit=u.keV, description="High energy electron cut off", fixed=True + name="high_e_cutoff", default=1500, unit=u.keV, description="High energy electron cut off", fixed=True, bounds=(FLOAT_EPSILON, None) ) total_eflux = Parameter( - name="total_eflux", default=1.5, unit=norm_thick_target_eflux_units, description="Total electron flux", fixed=False + name="total_eflux", default=1.5, unit=norm_thick_target_eflux_units, description="Total electron flux", fixed=False, bounds=(0, None) ) _input_units_allow_dimensionless = True diff --git a/sunkit_spex/models/physical/tests/test_nonthermal.py b/sunkit_spex/models/physical/tests/test_nonthermal.py index 4f175bef..28fc4413 100644 --- a/sunkit_spex/models/physical/tests/test_nonthermal.py +++ b/sunkit_spex/models/physical/tests/test_nonthermal.py @@ -169,3 +169,23 @@ def test_thin_target_flux_scaling(): np.testing.assert_allclose(model(energy_edges).value, s_model(energy_edges).value) np.testing.assert_allclose(model.evaluate(energy_edges, *model.parameters).value, s_model.evaluate(energy_edges, *s_model.parameters).value) + +def test_thick_target_parameter_check(): + """Test non-physical, erroneous values for parameters.""" + energy_edges = np.arange(2, 15, 0.1) << u.keV + model = nonthermal.ThickTarget() + + with pytest.warns(RuntimeWarning): + # produce a division by zero when calculating the internal n0 + p = 1 + model.evaluate(energy_edges, p, *model.parameters[1:]) + + with pytest.warns(RuntimeWarning): + # produce a division by zero when calculating the internal n0 + q = 1 + model.evaluate(energy_edges, *model.parameters[:2], q, *model.parameters[3:]) + + with pytest.warns(RuntimeWarning): + # produce a division by zero when calculating the internal n0 + ec = 0 + model.evaluate(energy_edges, *model.parameters[:3], ec, *model.parameters[4:]) diff --git a/sunkit_spex/models/physical/thermal.py b/sunkit_spex/models/physical/thermal.py index 9047fd73..7d98bd64 100644 --- a/sunkit_spex/models/physical/thermal.py +++ b/sunkit_spex/models/physical/thermal.py @@ -135,8 +135,6 @@ class ThermalEmission(FittableModel): temperature = Parameter( name="temperature", default=10, - min=1, - max=100, unit=u.MK, description="Temperature of the plasma", fixed=False, @@ -147,7 +145,8 @@ class ThermalEmission(FittableModel): default=1, unit=norm_thermal_emission_measure_units, description="Emission measure of the observer", - fixed=False, + fixed=False, + bounds=(0, None), ) mg = Parameter(name="Mg", default=8.15, min=6.15, max=10.15, description="Mg relative abundance", fixed=True) @@ -187,6 +186,8 @@ def __init__( mg, al, si, s, ar, ca, fe = _initialize_abundances(DEFAULT_ABUNDANCES[abundance_type]) emission_measure <<= norm_thermal_emission_measure_units + self.temperature.bounds = ((np.min([CONTINUUM_GRID["temperature range K"][0], LINE_GRID["temperature range K"][0]])< Date: Tue, 28 Jul 2026 16:25:56 -0500 Subject: [PATCH 4/9] pre-commit --- sunkit_spex/models/physical/nonthermal.py | 101 ++++++++++++------ .../models/physical/tests/test_nonthermal.py | 27 +++-- .../models/physical/tests/test_thermal.py | 14 ++- sunkit_spex/models/physical/thermal.py | 23 ++-- sunkit_spex/models/scaling.py | 10 +- 5 files changed, 118 insertions(+), 57 deletions(-) diff --git a/sunkit_spex/models/physical/nonthermal.py b/sunkit_spex/models/physical/nonthermal.py index 99c83b55..46a561e8 100644 --- a/sunkit_spex/models/physical/nonthermal.py +++ b/sunkit_spex/models/physical/nonthermal.py @@ -6,9 +6,9 @@ from astropy.modeling import FittableModel, Parameter from astropy.modeling.functional_models import FLOAT_EPSILON -from sunkit_spex.models.scaling import norm_thick_target_eflux_units, norm_thin_target_eflux_units from sunkit_spex.legacy import constants as const from sunkit_spex.legacy.integrate import gauss_legendre +from sunkit_spex.models.scaling import norm_thick_target_eflux_units, norm_thin_target_eflux_units const = const.Constants() @@ -31,7 +31,7 @@ __all__ = ["ThickTarget", "ThinTarget"] -FLOAT_EPSILON_FOR_POWER_LAW = 1+FLOAT_EPSILON*1e30 +FLOAT_EPSILON_FOR_POWER_LAW = 1 + FLOAT_EPSILON * 1e30 class ThickTarget(FittableModel): @@ -77,22 +77,48 @@ class ThickTarget(FittableModel): n_inputs = 1 n_outputs = 1 - p = Parameter(name="p", default=2, description="Slope below break", fixed=False, bounds=(FLOAT_EPSILON_FOR_POWER_LAW, None)) + p = Parameter( + name="p", default=2, description="Slope below break", fixed=False, bounds=(FLOAT_EPSILON_FOR_POWER_LAW, None) + ) - break_energy = Parameter(name="break_energy", default=100, unit=u.keV, description="Break Energy", fixed=True, bounds=(FLOAT_EPSILON, None)) + break_energy = Parameter( + name="break_energy", + default=100, + unit=u.keV, + description="Break Energy", + fixed=True, + bounds=(FLOAT_EPSILON, None), + ) - q = Parameter(name="q", default=5, description="Slope above break", fixed=True, bounds=(FLOAT_EPSILON_FOR_POWER_LAW, None)) + q = Parameter( + name="q", default=5, description="Slope above break", fixed=True, bounds=(FLOAT_EPSILON_FOR_POWER_LAW, None) + ) low_e_cutoff = Parameter( - name="low_e_cutoff", default=7, unit=u.keV, description="Low energy electron cut off", fixed=False, bounds=(FLOAT_EPSILON, None) + name="low_e_cutoff", + default=7, + unit=u.keV, + description="Low energy electron cut off", + fixed=False, + bounds=(FLOAT_EPSILON, None), ) high_e_cutoff = Parameter( - name="high_e_cutoff", default=1500, unit=u.keV, description="High energy electron cut off", fixed=True, bounds=(FLOAT_EPSILON, None) + name="high_e_cutoff", + default=1500, + unit=u.keV, + description="High energy electron cut off", + fixed=True, + bounds=(FLOAT_EPSILON, None), ) total_eflux = Parameter( - name="total_eflux", default=1.5, unit=norm_thick_target_eflux_units, description="Total electron flux", fixed=False, bounds=(0, None) + name="total_eflux", + default=1.5, + unit=norm_thick_target_eflux_units, + description="Total electron flux", + fixed=False, + bounds=(0, None), ) _input_units_allow_dimensionless = True @@ -130,23 +156,26 @@ def evaluate(self, energy_edges, p, break_energy, q, low_e_cutoff, high_e_cutoff low_e_cutoff <<= self.low_e_cutoff.unit high_e_cutoff <<= self.high_e_cutoff.unit total_eflux <<= norm_thick_target_eflux_units - - flux = bremsstrahlung_thick_target( - energy_centers.value, - p, - break_energy.value, - q, - low_e_cutoff.value, - high_e_cutoff.value, - self.integrator, - ) * total_eflux.decompose().value + + flux = ( + bremsstrahlung_thick_target( + energy_centers.value, + p, + break_energy.value, + q, + low_e_cutoff.value, + high_e_cutoff.value, + self.integrator, + ) + * total_eflux.decompose().value + ) return flux * self.return_units[self.outputs[0]] @property def input_units(self): # The units for the 'energy_edges' variable should be an energy (default keV) - return {self.inputs[0]: u.keV} + return {self.inputs[0]: u.keV} @property def return_units(self): @@ -160,6 +189,7 @@ def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): "total_eflux": norm_thick_target_eflux_units, } + class ThinTarget(FittableModel): r"""Calculates the thin-target bremsstrahlung radiation of a dual power-law electron distribution. @@ -224,7 +254,11 @@ class ThinTarget(FittableModel): ) total_eflux = Parameter( - name="total_eflux", default=1.5, unit=norm_thin_target_eflux_units, description="Total electron flux", fixed=True + name="total_eflux", + default=1.5, + unit=norm_thin_target_eflux_units, + description="Total electron flux", + fixed=True, ) _input_units_allow_dimensionless = True @@ -256,22 +290,25 @@ def __init__( def evaluate(self, energy_edges, p, break_energy, q, low_e_cutoff, high_e_cutoff, total_eflux): energy_centers = energy_edges[:-1] + 0.5 * np.diff(energy_edges) - + energy_centers <<= u.keV break_energy <<= self.break_energy.unit low_e_cutoff <<= self.low_e_cutoff.unit high_e_cutoff <<= self.high_e_cutoff.unit - total_eflux <<= norm_thin_target_eflux_units - - flux = bremsstrahlung_thin_target( - energy_centers.value, - p, - break_energy.value, - q, - low_e_cutoff.value, - high_e_cutoff.value, - integrator=self.integrator, - ) * total_eflux.to((u.electron * u.cm ** (-2) * u.s**-1)).value + total_eflux <<= norm_thin_target_eflux_units + + flux = ( + bremsstrahlung_thin_target( + energy_centers.value, + p, + break_energy.value, + q, + low_e_cutoff.value, + high_e_cutoff.value, + integrator=self.integrator, + ) + * total_eflux.to(u.electron * u.cm ** (-2) * u.s**-1).value + ) return flux * self.return_units[self.outputs[0]] diff --git a/sunkit_spex/models/physical/tests/test_nonthermal.py b/sunkit_spex/models/physical/tests/test_nonthermal.py index 28fc4413..87c01b33 100644 --- a/sunkit_spex/models/physical/tests/test_nonthermal.py +++ b/sunkit_spex/models/physical/tests/test_nonthermal.py @@ -146,6 +146,7 @@ def test_thin_target_against_ssw(ssw): expected_value = expected.to_value(output.unit) np.testing.assert_allclose(output.value, expected_value, rtol=0.035) + def test_thick_target_flux_scaling(): """Test thick target flux units being scaled.""" energy_edges = np.arange(2, 15, 0.1) << u.keV @@ -155,8 +156,11 @@ def test_thick_target_flux_scaling(): model = nonthermal.ThickTarget(total_eflux=eflux) s_model = nonthermal.ThickTarget(total_eflux=s_eflux) np.testing.assert_allclose(model(energy_edges).value, s_model(energy_edges).value) - np.testing.assert_allclose(model.evaluate(energy_edges, *model.parameters).value, - s_model.evaluate(energy_edges, *s_model.parameters).value) + np.testing.assert_allclose( + model.evaluate(energy_edges, *model.parameters).value, + s_model.evaluate(energy_edges, *s_model.parameters).value, + ) + def test_thin_target_flux_scaling(): """Test thin target flux units being scaled.""" @@ -167,25 +171,28 @@ def test_thin_target_flux_scaling(): model = nonthermal.ThinTarget(total_eflux=eflux) s_model = nonthermal.ThinTarget(total_eflux=s_eflux) np.testing.assert_allclose(model(energy_edges).value, s_model(energy_edges).value) - np.testing.assert_allclose(model.evaluate(energy_edges, *model.parameters).value, - s_model.evaluate(energy_edges, *s_model.parameters).value) + np.testing.assert_allclose( + model.evaluate(energy_edges, *model.parameters).value, + s_model.evaluate(energy_edges, *s_model.parameters).value, + ) + def test_thick_target_parameter_check(): """Test non-physical, erroneous values for parameters.""" energy_edges = np.arange(2, 15, 0.1) << u.keV model = nonthermal.ThickTarget() + # produce a division by zero when calculating the internal n0 + p = 1 with pytest.warns(RuntimeWarning): - # produce a division by zero when calculating the internal n0 - p = 1 model.evaluate(energy_edges, p, *model.parameters[1:]) + # produce a division by zero when calculating the internal n0 + q = 1 with pytest.warns(RuntimeWarning): - # produce a division by zero when calculating the internal n0 - q = 1 model.evaluate(energy_edges, *model.parameters[:2], q, *model.parameters[3:]) + # produce a division by zero when calculating the internal n0 + ec = 0 with pytest.warns(RuntimeWarning): - # produce a division by zero when calculating the internal n0 - ec = 0 model.evaluate(energy_edges, *model.parameters[:3], ec, *model.parameters[4:]) diff --git a/sunkit_spex/models/physical/tests/test_thermal.py b/sunkit_spex/models/physical/tests/test_thermal.py index a0d837ce..1a28aca4 100644 --- a/sunkit_spex/models/physical/tests/test_thermal.py +++ b/sunkit_spex/models/physical/tests/test_thermal.py @@ -5,8 +5,8 @@ import astropy.units as u -from sunkit_spex.models.scaling import norm_thermal_emission_measure_units from sunkit_spex.models.physical import thermal +from sunkit_spex.models.scaling import norm_thermal_emission_measure_units # Manually load file that was used to compile expected flux values. thermal.setup_continuum_parameters( @@ -515,6 +515,7 @@ def test_abundances_should_not_change(): after_models = thermal.DEFAULT_ABUNDANCES[thermal.DEFAULT_ABUNDANCE_TYPE].data assert np.allclose(after_models.data, orig.data) + def test_thermal_emission_measure_scaling(): """Test thermal emission measure units being scaled.""" energy_edges = np.arange(2, 15, 0.1) << u.keV @@ -524,7 +525,10 @@ def test_thermal_emission_measure_scaling(): model = thermal.ThermalEmission(emission_measure=em) s_model = thermal.ThermalEmission(emission_measure=s_em) np.testing.assert_allclose(model(energy_edges).value, s_model(energy_edges).value) - np.testing.assert_allclose(model.evaluate(energy_edges, *model.parameters).value, - s_model.evaluate(energy_edges, *s_model.parameters).value) - -test_thermal_emission_measure_scaling() \ No newline at end of file + np.testing.assert_allclose( + model.evaluate(energy_edges, *model.parameters).value, + s_model.evaluate(energy_edges, *s_model.parameters).value, + ) + + +test_thermal_emission_measure_scaling() diff --git a/sunkit_spex/models/physical/thermal.py b/sunkit_spex/models/physical/thermal.py index 7d98bd64..135104e2 100644 --- a/sunkit_spex/models/physical/thermal.py +++ b/sunkit_spex/models/physical/thermal.py @@ -145,7 +145,7 @@ class ThermalEmission(FittableModel): default=1, unit=norm_thermal_emission_measure_units, description="Emission measure of the observer", - fixed=False, + fixed=False, bounds=(0, None), ) @@ -186,9 +186,15 @@ def __init__( mg, al, si, s, ar, ca, fe = _initialize_abundances(DEFAULT_ABUNDANCES[abundance_type]) emission_measure <<= norm_thermal_emission_measure_units - self.temperature.bounds = ((np.min([CONTINUUM_GRID["temperature range K"][0], LINE_GRID["temperature range K"][0]])< Date: Wed, 29 Jul 2026 10:57:54 -0500 Subject: [PATCH 5/9] changelog added --- changelog/291.feature.rst | 1 + 1 file changed, 1 insertion(+) create mode 100644 changelog/291.feature.rst diff --git a/changelog/291.feature.rst b/changelog/291.feature.rst new file mode 100644 index 00000000..f33fefea --- /dev/null +++ b/changelog/291.feature.rst @@ -0,0 +1 @@ +The scalar parameters for ``ThermalEmission``, ``ThickTarget``, and ``ThinTarget`` in ``models.physical`` now use consistent custom scaled values and units instead of scaled values with non-scaled units. Parameters for the ``ThermalEmission`` and ``ThickTarget`` are given suitable parameter bounds. \ No newline at end of file From ceb016555c8af9e6b1341856c41d47ae292bfaa1 Mon Sep 17 00:00:00 2001 From: Kristopher Cooper Date: Wed, 29 Jul 2026 13:50:16 -0500 Subject: [PATCH 6/9] Made the new scaled units accessable to the Astropy methods like .find_equivalent_units() --- sunkit_spex/models/physical/nonthermal.py | 21 ++++++++++--------- .../models/physical/tests/test_nonthermal.py | 6 +++--- .../models/physical/tests/test_thermal.py | 4 ++-- sunkit_spex/models/physical/thermal.py | 18 ++++++++-------- sunkit_spex/models/scaling.py | 13 ++++++------ 5 files changed, 32 insertions(+), 30 deletions(-) diff --git a/sunkit_spex/models/physical/nonthermal.py b/sunkit_spex/models/physical/nonthermal.py index 46a561e8..42b4be3b 100644 --- a/sunkit_spex/models/physical/nonthermal.py +++ b/sunkit_spex/models/physical/nonthermal.py @@ -1,4 +1,5 @@ import logging +from collections.abc import Iterable import numpy as np @@ -8,7 +9,7 @@ from sunkit_spex.legacy import constants as const from sunkit_spex.legacy.integrate import gauss_legendre -from sunkit_spex.models.scaling import norm_thick_target_eflux_units, norm_thin_target_eflux_units +from sunkit_spex.models.scaling import scaled_thick_eflux_units, scaled_thin_eflux_units const = const.Constants() @@ -111,11 +112,11 @@ class ThickTarget(FittableModel): fixed=True, bounds=(FLOAT_EPSILON, None), ) - + total_eflux = Parameter( name="total_eflux", default=1.5, - unit=norm_thick_target_eflux_units, + unit=scaled_thick_eflux_units, description="Total electron flux", fixed=False, bounds=(0, None), @@ -136,7 +137,7 @@ def __init__( ): self.integrator = integrator - total_eflux <<= norm_thick_target_eflux_units + total_eflux <<= scaled_thick_eflux_units super().__init__( p=p, @@ -155,7 +156,7 @@ def evaluate(self, energy_edges, p, break_energy, q, low_e_cutoff, high_e_cutoff break_energy <<= self.break_energy.unit low_e_cutoff <<= self.low_e_cutoff.unit high_e_cutoff <<= self.high_e_cutoff.unit - total_eflux <<= norm_thick_target_eflux_units + total_eflux <<= scaled_thick_eflux_units flux = ( bremsstrahlung_thick_target( @@ -186,7 +187,7 @@ def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): "break_energy": u.keV, "low_e_cutoff": u.keV, "high_e_cutoff": u.keV, - "total_eflux": norm_thick_target_eflux_units, + "total_eflux": scaled_thick_eflux_units, } @@ -256,7 +257,7 @@ class ThinTarget(FittableModel): total_eflux = Parameter( name="total_eflux", default=1.5, - unit=norm_thin_target_eflux_units, + unit=scaled_thin_eflux_units, description="Total electron flux", fixed=True, ) @@ -276,7 +277,7 @@ def __init__( ): self.integrator = integrator - total_eflux <<= norm_thin_target_eflux_units + total_eflux <<= scaled_thin_eflux_units super().__init__( p=p, @@ -295,7 +296,7 @@ def evaluate(self, energy_edges, p, break_energy, q, low_e_cutoff, high_e_cutoff break_energy <<= self.break_energy.unit low_e_cutoff <<= self.low_e_cutoff.unit high_e_cutoff <<= self.high_e_cutoff.unit - total_eflux <<= norm_thin_target_eflux_units + total_eflux <<= scaled_thin_eflux_units flux = ( bremsstrahlung_thin_target( @@ -326,7 +327,7 @@ def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): "break_energy": u.keV, "low_e_cutoff": u.keV, "high_e_cutoff": u.keV, - "total_eflux": norm_thin_target_eflux_units, + "total_eflux": scaled_thin_eflux_units, # "total_eflux": u.electron * u.s**-1, } diff --git a/sunkit_spex/models/physical/tests/test_nonthermal.py b/sunkit_spex/models/physical/tests/test_nonthermal.py index 87c01b33..01d65208 100644 --- a/sunkit_spex/models/physical/tests/test_nonthermal.py +++ b/sunkit_spex/models/physical/tests/test_nonthermal.py @@ -4,7 +4,7 @@ import astropy.units as u from sunkit_spex.models.physical import nonthermal -from sunkit_spex.models.scaling import norm_thick_target_eflux_units, norm_thin_target_eflux_units +from sunkit_spex.models.scaling import scaled_thick_eflux_units, scaled_thin_eflux_units SSW_INTENSITY_UNIT = u.ph / u.cm**2 / u.s / u.keV @@ -152,7 +152,7 @@ def test_thick_target_flux_scaling(): energy_edges = np.arange(2, 15, 0.1) << u.keV for _flux in np.arange(1, 10, 0.5): eflux = (_flux * 1e35) << (u.electron * u.s**-1) - s_eflux = _flux << norm_thick_target_eflux_units + s_eflux = _flux << scaled_thick_eflux_units model = nonthermal.ThickTarget(total_eflux=eflux) s_model = nonthermal.ThickTarget(total_eflux=s_eflux) np.testing.assert_allclose(model(energy_edges).value, s_model(energy_edges).value) @@ -167,7 +167,7 @@ def test_thin_target_flux_scaling(): energy_edges = np.arange(2, 15, 0.1) << u.keV for _flux in np.arange(1, 10, 0.5): eflux = (_flux * 1e55) << (u.electron * u.cm**-2 * u.s**-1) - s_eflux = _flux << norm_thin_target_eflux_units + s_eflux = _flux << scaled_thin_eflux_units model = nonthermal.ThinTarget(total_eflux=eflux) s_model = nonthermal.ThinTarget(total_eflux=s_eflux) np.testing.assert_allclose(model(energy_edges).value, s_model(energy_edges).value) diff --git a/sunkit_spex/models/physical/tests/test_thermal.py b/sunkit_spex/models/physical/tests/test_thermal.py index 1a28aca4..a1eb6607 100644 --- a/sunkit_spex/models/physical/tests/test_thermal.py +++ b/sunkit_spex/models/physical/tests/test_thermal.py @@ -6,7 +6,7 @@ import astropy.units as u from sunkit_spex.models.physical import thermal -from sunkit_spex.models.scaling import norm_thermal_emission_measure_units +from sunkit_spex.models.scaling import scaled_em_units # Manually load file that was used to compile expected flux values. thermal.setup_continuum_parameters( @@ -521,7 +521,7 @@ def test_thermal_emission_measure_scaling(): energy_edges = np.arange(2, 15, 0.1) << u.keV for _em in np.arange(1, 10, 0.5): em = (_em * 1e49) << (u.cm**-3) - s_em = _em << norm_thermal_emission_measure_units + s_em = _em << scaled_em_units model = thermal.ThermalEmission(emission_measure=em) s_model = thermal.ThermalEmission(emission_measure=s_em) np.testing.assert_allclose(model(energy_edges).value, s_model(energy_edges).value) diff --git a/sunkit_spex/models/physical/thermal.py b/sunkit_spex/models/physical/thermal.py index 135104e2..d54f9625 100644 --- a/sunkit_spex/models/physical/thermal.py +++ b/sunkit_spex/models/physical/thermal.py @@ -15,7 +15,7 @@ load_chianti_lines_lite, load_xray_abundances, ) -from sunkit_spex.models.scaling import norm_thermal_emission_measure_units +from sunkit_spex.models.scaling import scaled_em_units # The default elemental abundance values correspond to coronal values DEFAULT_ABUNDANCE_TYPE = "sun_coronal_ext" @@ -143,7 +143,7 @@ class ThermalEmission(FittableModel): emission_measure = Parameter( name="emission_measure", default=1, - unit=norm_thermal_emission_measure_units, + unit=scaled_em_units, description="Emission measure of the observer", fixed=False, bounds=(0, None), @@ -185,7 +185,7 @@ def __init__( if abundance_type != DEFAULT_ABUNDANCE_TYPE: mg, al, si, s, ar, ca, fe = _initialize_abundances(DEFAULT_ABUNDANCES[abundance_type]) - emission_measure <<= norm_thermal_emission_measure_units + emission_measure <<= scaled_em_units self.temperature.bounds = ( (np.min([CONTINUUM_GRID["temperature range K"][0], LINE_GRID["temperature range K"][0]]) << u.K).to( temperature.unit @@ -288,7 +288,7 @@ def return_units(self): return {self.outputs[0]: u.ph / u.keV * u.s**-1} def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): - return {"temperature": u.MK, "emission_measure": norm_thermal_emission_measure_units} + return {"temperature": u.MK, "emission_measure": scaled_em_units} class ContinuumEmission(FittableModel): @@ -342,7 +342,7 @@ class ContinuumEmission(FittableModel): emission_measure = Parameter( name="emission_measure", default=1, - unit=norm_thermal_emission_measure_units, + unit=scaled_em_units, description="Emission measure of the observer", fixed=False, ) @@ -432,7 +432,7 @@ def return_units(self): return {self.outputs[0]: u.ph / u.keV * u.s**-1} def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): - return {"temperature": u.MK, "emission_measure": norm_thermal_emission_measure_units} + return {"temperature": u.MK, "emission_measure": scaled_em_units} class LineEmission(FittableModel): @@ -482,7 +482,7 @@ class LineEmission(FittableModel): emission_measure = Parameter( name="emission_measure", default=1, - unit=norm_thermal_emission_measure_units, + unit=scaled_em_units, description="Emission measure of the observer", fixed=False, ) @@ -572,7 +572,7 @@ def return_units(self): return {self.outputs[0]: u.ph / u.keV * u.s**-1} def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): - return {"temperature": u.MK, "emission_measure": norm_thermal_emission_measure_units} + return {"temperature": u.MK, "emission_measure": scaled_em_units} def setup_continuum_parameters(filename=None): @@ -1236,7 +1236,7 @@ def _sanitize_inputs(energy_edges, temperature, emission_measure): # If they were not already Quantities, the parameters get the default units. energy_edges <<= u.keV temperature <<= u.K - emission_measure <<= norm_thermal_emission_measure_units + emission_measure <<= scaled_em_units energy_edges_keV = energy_edges.to(u.keV) diff --git a/sunkit_spex/models/scaling.py b/sunkit_spex/models/scaling.py index 349ffff6..e1c6dd02 100644 --- a/sunkit_spex/models/scaling.py +++ b/sunkit_spex/models/scaling.py @@ -7,15 +7,16 @@ __all__ = [ "Constant", "InverseSquareFluxScaling", - "norm_thermal_emission_measure_units", - "norm_thick_target_eflux_units", - "norm_thin_target_eflux_units", + "scaled_em_units", + "scaled_thick_eflux_units", + "scaled_thin_eflux_units", ] -norm_thick_target_eflux_units = u.def_unit("scaled_thick_eflux_units", 1e35 * (u.electron * u.s**-1)) -norm_thin_target_eflux_units = u.def_unit("scaled_thin_eflux_units", 1e55 * (u.electron * u.cm ** (-2) * u.s**-1)) -norm_thermal_emission_measure_units = u.def_unit("scaled_em_units", 1e49 * (u.cm ** (-3))) +scaled_thick_eflux_units = u.def_unit("scaled_thick_eflux_units", 1e35 * (u.electron * u.s**-1)) +scaled_thin_eflux_units = u.def_unit("scaled_thin_eflux_units", 1e55 * (u.electron * u.cm ** (-2) * u.s**-1)) +scaled_em_units = u.def_unit("scaled_em_units", 1e49 * (u.cm ** (-3))) +u.add_enabled_units([scaled_thick_eflux_units, scaled_thin_eflux_units, scaled_em_units]) # lets astropy's methods find these units class InverseSquareFluxScaling(FittableModel): From 4a0c3645341e391b09aa2c590ebf58c31f77e678 Mon Sep 17 00:00:00 2001 From: Kristopher Cooper Date: Wed, 29 Jul 2026 13:51:46 -0500 Subject: [PATCH 7/9] pre-commit --- changelog/291.feature.rst | 2 +- sunkit_spex/models/physical/nonthermal.py | 3 +-- sunkit_spex/models/scaling.py | 4 +++- 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/changelog/291.feature.rst b/changelog/291.feature.rst index f33fefea..06932b28 100644 --- a/changelog/291.feature.rst +++ b/changelog/291.feature.rst @@ -1 +1 @@ -The scalar parameters for ``ThermalEmission``, ``ThickTarget``, and ``ThinTarget`` in ``models.physical`` now use consistent custom scaled values and units instead of scaled values with non-scaled units. Parameters for the ``ThermalEmission`` and ``ThickTarget`` are given suitable parameter bounds. \ No newline at end of file +The scalar parameters for ``ThermalEmission``, ``ThickTarget``, and ``ThinTarget`` in ``models.physical`` now use consistent custom scaled values and units instead of scaled values with non-scaled units. Parameters for the ``ThermalEmission`` and ``ThickTarget`` are given suitable parameter bounds. diff --git a/sunkit_spex/models/physical/nonthermal.py b/sunkit_spex/models/physical/nonthermal.py index 42b4be3b..94234fe8 100644 --- a/sunkit_spex/models/physical/nonthermal.py +++ b/sunkit_spex/models/physical/nonthermal.py @@ -1,5 +1,4 @@ import logging -from collections.abc import Iterable import numpy as np @@ -112,7 +111,7 @@ class ThickTarget(FittableModel): fixed=True, bounds=(FLOAT_EPSILON, None), ) - + total_eflux = Parameter( name="total_eflux", default=1.5, diff --git a/sunkit_spex/models/scaling.py b/sunkit_spex/models/scaling.py index e1c6dd02..e57ab61f 100644 --- a/sunkit_spex/models/scaling.py +++ b/sunkit_spex/models/scaling.py @@ -16,7 +16,9 @@ scaled_thick_eflux_units = u.def_unit("scaled_thick_eflux_units", 1e35 * (u.electron * u.s**-1)) scaled_thin_eflux_units = u.def_unit("scaled_thin_eflux_units", 1e55 * (u.electron * u.cm ** (-2) * u.s**-1)) scaled_em_units = u.def_unit("scaled_em_units", 1e49 * (u.cm ** (-3))) -u.add_enabled_units([scaled_thick_eflux_units, scaled_thin_eflux_units, scaled_em_units]) # lets astropy's methods find these units +u.add_enabled_units( + [scaled_thick_eflux_units, scaled_thin_eflux_units, scaled_em_units] +) # lets astropy's methods find these units class InverseSquareFluxScaling(FittableModel): From db1165f418088d5e46d39f582c3dad23b94fb4e9 Mon Sep 17 00:00:00 2001 From: Kristopher Cooper Date: Wed, 29 Jul 2026 13:55:54 -0500 Subject: [PATCH 8/9] Removed accidental call to test function I was checking. --- sunkit_spex/models/physical/tests/test_thermal.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/sunkit_spex/models/physical/tests/test_thermal.py b/sunkit_spex/models/physical/tests/test_thermal.py index a1eb6607..fe60728c 100644 --- a/sunkit_spex/models/physical/tests/test_thermal.py +++ b/sunkit_spex/models/physical/tests/test_thermal.py @@ -529,6 +529,3 @@ def test_thermal_emission_measure_scaling(): model.evaluate(energy_edges, *model.parameters).value, s_model.evaluate(energy_edges, *s_model.parameters).value, ) - - -test_thermal_emission_measure_scaling() From cd8d2ec4d313b47b0597f26abb0c03a23958f97c Mon Sep 17 00:00:00 2001 From: Kristopher Cooper Date: Wed, 29 Jul 2026 14:27:57 -0500 Subject: [PATCH 9/9] use <<= to assign or convert albedo theta input to degrees in evaluate. Similar to other functions. --- sunkit_spex/models/physical/albedo.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sunkit_spex/models/physical/albedo.py b/sunkit_spex/models/physical/albedo.py index 8fa6e758..87ded2f8 100644 --- a/sunkit_spex/models/physical/albedo.py +++ b/sunkit_spex/models/physical/albedo.py @@ -94,8 +94,8 @@ def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) def evaluate(self, spectrum, theta, anisotropy): - if not isinstance(theta, Quantity): - theta = theta * u.deg + + theta <<= u.deg albedo_matrix = get_albedo_matrix(self.energy_edges, theta, anisotropy)