From 39a94dbee0da233a02cf94c1c2e2c659296f582f Mon Sep 17 00:00:00 2001 From: Kristopher Cooper Date: Wed, 29 Jul 2026 16:45:45 -0500 Subject: [PATCH 1/5] Attempt at warm thick target. Maybe just have the thermal model wrapped and make a thick tharget warm contribution model that should have its e_c param tied to a thick target param e_c? --- sunkit_spex/models/physical/nonthermal.py | 207 +++++++++++++++++++++- 1 file changed, 206 insertions(+), 1 deletion(-) diff --git a/sunkit_spex/models/physical/nonthermal.py b/sunkit_spex/models/physical/nonthermal.py index 1e95e29a..714d4378 100644 --- a/sunkit_spex/models/physical/nonthermal.py +++ b/sunkit_spex/models/physical/nonthermal.py @@ -27,7 +27,7 @@ """ -__all__ = ["ThickTarget", "ThinTarget"] +__all__ = ["ThickTarget", "ThinTarget", "WarmThickTarget"] class ThickTarget(FittableModel): @@ -1238,3 +1238,208 @@ def bremsstrahlung_thick_target(photon_energies, p, break_energy, q, low_e_cutof return (fcoeff / decoeff) * flux raise Warning("The photon energies are higher than the highest electron energy or not greater than zero") + + +from sunkit_spex.models.physical.thermal import ThermalEmission +from astropy.modeling.functional_models import FLOAT_EPSILON +class WarmThickTarget(FittableModel): + r"""Calculates the warm thick-target bremsstrahlung radiation of a single power-law electron distribution. + + [1] Kontar et al, ApJ 2015 (http://adsabs.harvard.edu/abs/2015arXiv150503733K) + [2] https://hesperia.gsfc.nasa.gov/ssw/packages/xray/idl/f_thick_warm.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. + """ + + FLOAT_EPSILON_FOR_POWER_LAW = 1 + FLOAT_EPSILON * 1e30 + scaled_thick_eflux_units = u.def_unit("scaled_thick_eflux_units", 1e35 * (u.electron * u.s**-1)) + scaled_warmthick_desnity_units = u.def_unit("scaled_warmthick_desnity_units", 1e10 * (u.cm**-3)) + + 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) + ) + + 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), + ) + + total_eflux = Parameter( + name="total_eflux", + default=1.5, + unit=scaled_thick_eflux_units, + description="Total electron flux", + fixed=False, + bounds=(0, None), + ) + + temperature = Parameter( + name="temperature", + default=10, + unit=u.MK, + description="Temperature of the plasma", + fixed=False, + bounds=(1, 100) + ) + + plasma_density = Parameter( + name="plasma_density", + default=1, + unit=scaled_warmthick_desnity_units, + description="Number density of the plasma", + fixed=False, + bounds=(FLOAT_EPSILON, None) + ) + + length = Parameter( + name="length", + default=1, + unit=u.Mm, + description="Plasma column length", + fixed=False, + bounds=(FLOAT_EPSILON, None) + ) + + mg = Parameter(name="Mg", default=8.15, min=6.15, max=10.15, description="Mg relative abundance", fixed=True) + + al = Parameter(name="Al", default=7.04, min=5.04, max=9.04, description="Al relative abundance", fixed=True) + + si = Parameter(name="Si", default=8.1, min=6.1, max=10.1, description="Si relative abundance", fixed=True) + + s = Parameter(name="S", default=7.27, min=5.27, max=9.27, description="S relative abundance", fixed=True) + + ar = Parameter(name="Ar", default=6.58, min=4.58, max=8.58, description="Ar relative abundance", fixed=True) + + ca = Parameter(name="Ca", default=6.93, min=4.93, max=8.93, description="Ca relative abundance", fixed=True) + + fe = Parameter(name="Fe", default=8.1, min=6.1, max=10.1, description="Fe relative abundance", fixed=True) + + _input_units_allow_dimensionless = True + + def __init__( + self, + p=p.default, + low_e_cutoff=u.Quantity(low_e_cutoff.default, low_e_cutoff.unit), + total_eflux=u.Quantity(total_eflux.default, total_eflux.unit), + temperature=u.Quantity(temperature.default, temperature.unit), + plasma_density=u.Quantity(plasma_density.default, plasma_density.unit), + length=u.Quantity(length.default, length.unit), + mg=mg.default, + al=al.default, + si=si.default, + s=s.default, + ar=ar.default, + ca=ca.default, + fe=fe.default, + abundance_type=DEFAULT_ABUNDANCE_TYPE, + **kwargs, + ): + self.integrator = integrator + + total_eflux <<= self.scaled_thick_eflux_units + plasma_density <<= self.scaled_warmthick_desnity_units + + self.line = ThermalEmission( + temperature=temperature, + emission_measure=emission_measure, + mg=mg, + al=al, + si=si, + s=s, + ar=ar, + ca=ca, + fe=fe, + abundance_type=abundance_type, + ) + + self.cont = ContinuumEmission( + temperature=temperature, + emission_measure=emission_measure, + mg=mg, + al=al, + si=si, + s=s, + ar=ar, + ca=ca, + fe=fe, + abundance_type=abundance_type, + ) + + super().__init__( + p=p, + low_e_cutoff=low_e_cutoff, + total_eflux=total_eflux, + temperature=temperature, + plasma_density=plasma_density, + length=length, + **kwargs, + ) + + def evaluate(self, energy_edges, p, low_e_cutoff, total_eflux, temperature, plasma_density, length): + energy_centers = energy_edges[:-1] + 0.5 * np.diff(energy_edges) + + energy_centers <<= u.keV + low_e_cutoff <<= self.low_e_cutoff.unit + total_eflux <<= self.scaled_thick_eflux_units + temperature <<= self.temperature.unit + plasma_density <<= self.scaled_warmthick_desnity_units + length <<= self.length.unit + + high_and_break = np.float64(energy_edges.value.max() * 10) + q = 20 + + return flux + + @property + def input_units(self): + # The units for the 'energy_edges' variable should be an energy (default keV) + return {self.inputs[0]: u.keV} + + @property + def return_units(self): + return {self.outputs[0]: u.ph * u.keV**-1 * u.s**-1} + + def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): + return { + "break_energy": u.keV, + "low_e_cutoff": u.keV, + "high_e_cutoff": u.keV, + "total_eflux": u.electron * u.s**-1, + } From e94acf2246e99a447df9c5845b4aa05883b49595 Mon Sep 17 00:00:00 2001 From: Kristopher Cooper Date: Thu, 30 Jul 2026 14:56:34 -0500 Subject: [PATCH 2/5] Added warm thick contribution model --- sunkit_spex/models/physical/nonthermal.py | 208 +---------------- sunkit_spex/models/physical/thermal.py | 267 +++++++++++++++++++++- 2 files changed, 268 insertions(+), 207 deletions(-) diff --git a/sunkit_spex/models/physical/nonthermal.py b/sunkit_spex/models/physical/nonthermal.py index 714d4378..96528b1e 100644 --- a/sunkit_spex/models/physical/nonthermal.py +++ b/sunkit_spex/models/physical/nonthermal.py @@ -27,7 +27,7 @@ """ -__all__ = ["ThickTarget", "ThinTarget", "WarmThickTarget"] +__all__ = ["ThickTarget", "ThinTarget"] class ThickTarget(FittableModel): @@ -1238,208 +1238,4 @@ def bremsstrahlung_thick_target(photon_energies, p, break_energy, q, low_e_cutof return (fcoeff / decoeff) * flux raise Warning("The photon energies are higher than the highest electron energy or not greater than zero") - - -from sunkit_spex.models.physical.thermal import ThermalEmission -from astropy.modeling.functional_models import FLOAT_EPSILON -class WarmThickTarget(FittableModel): - r"""Calculates the warm thick-target bremsstrahlung radiation of a single power-law electron distribution. - - [1] Kontar et al, ApJ 2015 (http://adsabs.harvard.edu/abs/2015arXiv150503733K) - [2] https://hesperia.gsfc.nasa.gov/ssw/packages/xray/idl/f_thick_warm.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. - """ - - FLOAT_EPSILON_FOR_POWER_LAW = 1 + FLOAT_EPSILON * 1e30 - scaled_thick_eflux_units = u.def_unit("scaled_thick_eflux_units", 1e35 * (u.electron * u.s**-1)) - scaled_warmthick_desnity_units = u.def_unit("scaled_warmthick_desnity_units", 1e10 * (u.cm**-3)) - - 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) - ) - - 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), - ) - - total_eflux = Parameter( - name="total_eflux", - default=1.5, - unit=scaled_thick_eflux_units, - description="Total electron flux", - fixed=False, - bounds=(0, None), - ) - - temperature = Parameter( - name="temperature", - default=10, - unit=u.MK, - description="Temperature of the plasma", - fixed=False, - bounds=(1, 100) - ) - - plasma_density = Parameter( - name="plasma_density", - default=1, - unit=scaled_warmthick_desnity_units, - description="Number density of the plasma", - fixed=False, - bounds=(FLOAT_EPSILON, None) - ) - - length = Parameter( - name="length", - default=1, - unit=u.Mm, - description="Plasma column length", - fixed=False, - bounds=(FLOAT_EPSILON, None) - ) - - mg = Parameter(name="Mg", default=8.15, min=6.15, max=10.15, description="Mg relative abundance", fixed=True) - - al = Parameter(name="Al", default=7.04, min=5.04, max=9.04, description="Al relative abundance", fixed=True) - - si = Parameter(name="Si", default=8.1, min=6.1, max=10.1, description="Si relative abundance", fixed=True) - - s = Parameter(name="S", default=7.27, min=5.27, max=9.27, description="S relative abundance", fixed=True) - - ar = Parameter(name="Ar", default=6.58, min=4.58, max=8.58, description="Ar relative abundance", fixed=True) - - ca = Parameter(name="Ca", default=6.93, min=4.93, max=8.93, description="Ca relative abundance", fixed=True) - - fe = Parameter(name="Fe", default=8.1, min=6.1, max=10.1, description="Fe relative abundance", fixed=True) - - _input_units_allow_dimensionless = True - - def __init__( - self, - p=p.default, - low_e_cutoff=u.Quantity(low_e_cutoff.default, low_e_cutoff.unit), - total_eflux=u.Quantity(total_eflux.default, total_eflux.unit), - temperature=u.Quantity(temperature.default, temperature.unit), - plasma_density=u.Quantity(plasma_density.default, plasma_density.unit), - length=u.Quantity(length.default, length.unit), - mg=mg.default, - al=al.default, - si=si.default, - s=s.default, - ar=ar.default, - ca=ca.default, - fe=fe.default, - abundance_type=DEFAULT_ABUNDANCE_TYPE, - **kwargs, - ): - self.integrator = integrator - - total_eflux <<= self.scaled_thick_eflux_units - plasma_density <<= self.scaled_warmthick_desnity_units - - self.line = ThermalEmission( - temperature=temperature, - emission_measure=emission_measure, - mg=mg, - al=al, - si=si, - s=s, - ar=ar, - ca=ca, - fe=fe, - abundance_type=abundance_type, - ) - - self.cont = ContinuumEmission( - temperature=temperature, - emission_measure=emission_measure, - mg=mg, - al=al, - si=si, - s=s, - ar=ar, - ca=ca, - fe=fe, - abundance_type=abundance_type, - ) - - super().__init__( - p=p, - low_e_cutoff=low_e_cutoff, - total_eflux=total_eflux, - temperature=temperature, - plasma_density=plasma_density, - length=length, - **kwargs, - ) - - def evaluate(self, energy_edges, p, low_e_cutoff, total_eflux, temperature, plasma_density, length): - energy_centers = energy_edges[:-1] + 0.5 * np.diff(energy_edges) - - energy_centers <<= u.keV - low_e_cutoff <<= self.low_e_cutoff.unit - total_eflux <<= self.scaled_thick_eflux_units - temperature <<= self.temperature.unit - plasma_density <<= self.scaled_warmthick_desnity_units - length <<= self.length.unit - - high_and_break = np.float64(energy_edges.value.max() * 10) - q = 20 - - return flux - - @property - def input_units(self): - # The units for the 'energy_edges' variable should be an energy (default keV) - return {self.inputs[0]: u.keV} - - @property - def return_units(self): - return {self.outputs[0]: u.ph * u.keV**-1 * u.s**-1} - - def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): - return { - "break_energy": u.keV, - "low_e_cutoff": u.keV, - "high_e_cutoff": u.keV, - "total_eflux": u.electron * u.s**-1, - } + \ No newline at end of file diff --git a/sunkit_spex/models/physical/thermal.py b/sunkit_spex/models/physical/thermal.py index 73d09897..5aa99f13 100644 --- a/sunkit_spex/models/physical/thermal.py +++ b/sunkit_spex/models/physical/thermal.py @@ -5,7 +5,9 @@ from scipy import interpolate, stats import astropy.units as u +from astropy import constants as const from astropy.modeling import FittableModel, Parameter +from astropy.modeling.functional_models import FLOAT_EPSILON from astropy.table.column import Column from sunpy.data import manager @@ -19,7 +21,7 @@ # The default elemental abundance values correspond to coronal values DEFAULT_ABUNDANCE_TYPE = "sun_coronal_ext" -__all__ = ["ContinuumEmission", "LineEmission", "ThermalEmission"] +__all__ = ["ContinuumEmission", "LineEmission", "ThermalEmission", "ThickTargetWarmContribution"] doc_string_params = """ Parameters @@ -1309,3 +1311,266 @@ def _calculate_abundances(abundance_type, mg, al, si, s, ar, ca, fe): abundances[25] = 10 ** (fe - 12) return abundances + + + +class ThickTargetWarmContribution(FittableModel): + r""" + Calculates the warm thick-target bremsstrahlung radiation contribution + of a single power-law electron distribution with a given low-energy + cut-off energy. + + [1] Kontar et al, ApJ 2015 (http://adsabs.harvard.edu/abs/2015arXiv150503733K) + [2] https://hesperia.gsfc.nasa.gov/ssw/packages/xray/idl/f_thick_warm.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. + + low_e_cutoff : int or float + Low-energy cut-off of the electron distribution in units of keV. + + temperature: `astropy.units.Quantity` + The temperature of the plasma. + Can be scalar or 1D of any length. If not scalar, the flux for each temperature + will be calculated. The first dimension of the output flux will correspond + to temperature. + + plasma_density: `astropy.units.Quantity` + The number density o the plasma. + + length: `astropy.units.Quantity` + The plasma column length. + + Returns + ------- + A 1d array of warm component from thick-target bremsstrahlung radiation + in units of ph s^-1 keV^-1. + """ + + scaled_warmthick_desnity_units = u.def_unit("scaled_warmthick_desnity_units", 1e10 * (u.cm**-3)) + scaled_thick_eflux_units = u.def_unit("scaled_thick_eflux_units", 1e35 * (u.electron * u.s**-1)) + scaled_em_units = u.def_unit("scaled_em_units", 1e49 * (u.cm ** (-3))) + + name = "ThickTargetWarmContribution" + n_inputs = 1 + n_outputs = 1 + + 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), + ) + + total_eflux = Parameter( + name="total_eflux", + default=1.5, + unit=scaled_thick_eflux_units, + description="Total electron flux", + fixed=False, + bounds=(0, None), + ) + + plasma_density = Parameter( + name="plasma_density", + default=1, + unit=scaled_warmthick_desnity_units, + description="Number density of the plasma", + fixed=False, + bounds=(FLOAT_EPSILON, None) + ) + + length = Parameter( + name="length", + default=10, + unit=u.Mm, + description="Plasma column length", + fixed=False, + bounds=(FLOAT_EPSILON, None) + ) + + + temperature = Parameter( + name="temperature", + default=10, + min=1, + max=100, + unit=u.MK, + description="Temperature of the plasma", + fixed=False, + ) + + mg = Parameter(name="Mg", default=8.15, min=6.15, max=10.15, description="Mg relative abundance", fixed=True) + + al = Parameter(name="Al", default=7.04, min=5.04, max=9.04, description="Al relative abundance", fixed=True) + + si = Parameter(name="Si", default=8.1, min=6.1, max=10.1, description="Si relative abundance", fixed=True) + + s = Parameter(name="S", default=7.27, min=5.27, max=9.27, description="S relative abundance", fixed=True) + + ar = Parameter(name="Ar", default=6.58, min=4.58, max=8.58, description="Ar relative abundance", fixed=True) + + ca = Parameter(name="Ca", default=6.93, min=4.93, max=8.93, description="Ca relative abundance", fixed=True) + + fe = Parameter(name="Fe", default=8.1, min=6.1, max=10.1, description="Fe relative abundance", fixed=True) + + _input_units_allow_dimensionless = True + + def __init__( + self, + low_e_cutoff=u.Quantity(low_e_cutoff.default, low_e_cutoff.unit), + total_eflux=u.Quantity(total_eflux.default, total_eflux.unit), + plasma_density=u.Quantity(plasma_density.default, plasma_density.unit), + length=u.Quantity(length.default, length.unit), + temperature=u.Quantity(temperature.default, temperature.unit), + mg=mg.default, + al=al.default, + si=si.default, + s=s.default, + ar=ar.default, + ca=ca.default, + fe=fe.default, + abundance_type=DEFAULT_ABUNDANCE_TYPE, + **kwargs, + ): + + total_eflux <<= self.scaled_thick_eflux_units + plasma_density <<= self.scaled_warmthick_desnity_units + length <<= u.Mm + + self.c_cgs = const.c.cgs + self.me_e = (const.m_e * const.c**2).to(u.keV) + self.mk2keV = const.k_B.to(u.keV/u.MK) + # collision parameter, units assigned by Kris to make the maths work later + self.col_param = 2.6e-18 << (u.cm**2 * u.keV**2) + + self.therm_model = ThermalEmission(temperature=temperature, + mg=mg, + al=al, + si=si, + s=s, + ar=ar, + ca=ca, + fe=fe, + abundance_type=abundance_type, + **kwargs, + ) + + super().__init__( + low_e_cutoff=low_e_cutoff, + total_eflux=total_eflux, + plasma_density=plasma_density, + length=length, + temperature=self.therm_model.temperature, + mg=self.therm_model.mg, + al=self.therm_model.al, + si=self.therm_model.si, + s=self.therm_model.s, + ar=self.therm_model.ar, + ca=self.therm_model.ca, + fe=self.therm_model.fe, + **kwargs, + ) + + def calculate_emission_measure(self, temperature, plasma_density, length, low_e_cutoff, total_eflux): + """ + Calculate the emission measure of the warm thick target contribution. + + Need some maths and units magic here. + """ + + # convert units + loop_temp_kev = (temperature*self.mk2keV).to(u.keV) + plasma_density_cm = plasma_density.to(u.cm**-3) + length_cm = length.to(u.cm) + eflux = total_eflux.to(u.electron/u.second) + + col_stop_dist = loop_temp_kev**2 / (2 * self.col_param * plasma_density_cm) # collisional stopping distance for electrons of Tloop energy + emin = loop_temp_kev * 3 * (5 * col_stop_dist / length_cm) ** 4 + + if emin > 0.1<0.1 keV. Fixing emin to 0.1." + ) + emin = 0.1< length_cm: + warnings.warn(f"Minimum length ({lmin.to(u.Mm)}) > loop length ({length})") + + emission_measure = 3 * np.pi / 2 / self.col_param / self.c_cgs * np.sqrt(self.me_e / 8.0) * loop_temp_kev**2 / np.sqrt(emin) * eflux + emission_measure /= u.electron # electron interactions but let's remove the unit here + + return emission_measure << self.scaled_em_units + + def evaluate(self, + energy_edges, + low_e_cutoff, + total_eflux, + plasma_density, + length, + temperature, + mg, + al, + si, + s, + ar, + ca, + fe, + ): + + low_e_cutoff <<= self.low_e_cutoff.unit + total_eflux <<= self.scaled_thick_eflux_units + temperature <<= self.temperature.unit + plasma_density <<= self.scaled_warmthick_desnity_units + length <<= self.length.unit + + emission_measure = self.calculate_emission_measure( + temperature, + plasma_density, + length, + low_e_cutoff, + total_eflux + ) + + # TODO: this version of Thermal Emission still uses the scaled value with the non-scaled units + # so need to scale the value here with the inconsistent units + # Once ixed, replace ``emission_measure.value< Date: Fri, 31 Jul 2026 09:42:46 -0500 Subject: [PATCH 3/5] Added function that does all the warm thick target emission. --- sunkit_spex/models/physical/nonthermal.py | 275 +++++++++++++++++++++- sunkit_spex/models/physical/thermal.py | 16 +- 2 files changed, 282 insertions(+), 9 deletions(-) diff --git a/sunkit_spex/models/physical/nonthermal.py b/sunkit_spex/models/physical/nonthermal.py index 96528b1e..5bad9479 100644 --- a/sunkit_spex/models/physical/nonthermal.py +++ b/sunkit_spex/models/physical/nonthermal.py @@ -1238,4 +1238,277 @@ def bremsstrahlung_thick_target(photon_energies, p, break_energy, q, low_e_cutof return (fcoeff / decoeff) * flux raise Warning("The photon energies are higher than the highest electron energy or not greater than zero") - \ No newline at end of file + +from astropy.modeling.functional_models import FLOAT_EPSILON +from sunkit_spex.models.physical.thermal import ThickTargetWarmContribution, DEFAULT_ABUNDANCE_TYPE +class WarmThickTarget(FittableModel): + r""" + Calculates the thick-target + thermal contribution bremsstrahlung + radiation contribution. + + [1] Kontar et al, ApJ 2015 (http://adsabs.harvard.edu/abs/2015arXiv150503733K) + [2] https://hesperia.gsfc.nasa.gov/ssw/packages/xray/idl/f_thick_warm.pro + [3] https://www.astro.gla.ac.uk/users/natasha/rhessi_wt_tutorial_2017.pdf + + 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. + + plasma_density: `astropy.units.Quantity` + The number density o the plasma. + + length: `astropy.units.Quantity` + The plasma column length. + + temperature: `astropy.units.Quantity` + The temperature of the plasma. + Can be scalar or 1D of any length. If not scalar, the flux for each temperature + will be calculated. The first dimension of the output flux will correspond + to temperature. + + Returns + ------- + A 1d array of warm component from thick-target bremsstrahlung radiation + in units of ph s^-1 keV^-1. + """ + + scaled_warmthick_desnity_units = u.def_unit("scaled_warmthick_desnity_units", 1e10 * (u.cm**-3)) + scaled_thick_eflux_units = u.def_unit("scaled_thick_eflux_units", 1e35 * (u.electron * u.s**-1)) + scaled_em_units = u.def_unit("scaled_em_units", 1e49 * (u.cm ** (-3))) + + name = "WarmThickTarget" + n_inputs = 1 + n_outputs = 1 + + p = Parameter(name="p", default=2, description="Slope below break", fixed=False) + + break_energy = Parameter(name="break_energy", default=100, unit=u.keV, description="Break Energy", fixed=False) + + q = Parameter(name="q", default=5, min=0.01, description="Slope above break", fixed=True) + + low_e_cutoff = Parameter( + name="low_e_cutoff", default=7, unit=u.keV, description="Low energy electron cut off", fixed=False + ) + + high_e_cutoff = Parameter( + name="high_e_cutoff", default=1500, unit=u.keV, description="High energy electron cut off", fixed=True + ) + + total_eflux = Parameter( + name="total_eflux", default=1.5, unit=u.electron * u.s**-1, description="Total electron flux", fixed=True + ) + + plasma_density = Parameter( + name="plasma_density", + default=1, + unit=scaled_warmthick_desnity_units, + description="Number density of the plasma", + fixed=False, + bounds=(FLOAT_EPSILON, None) + ) + + length = Parameter( + name="length", + default=10, + unit=u.Mm, + description="Plasma column length", + fixed=False, + bounds=(FLOAT_EPSILON, None) + ) + + + temperature = Parameter( + name="temperature", + default=10, + min=1, + max=100, + unit=u.MK, + description="Temperature of the plasma", + fixed=False, + ) + + mg = Parameter(name="Mg", default=8.15, min=6.15, max=10.15, description="Mg relative abundance", fixed=True) + + al = Parameter(name="Al", default=7.04, min=5.04, max=9.04, description="Al relative abundance", fixed=True) + + si = Parameter(name="Si", default=8.1, min=6.1, max=10.1, description="Si relative abundance", fixed=True) + + s = Parameter(name="S", default=7.27, min=5.27, max=9.27, description="S relative abundance", fixed=True) + + ar = Parameter(name="Ar", default=6.58, min=4.58, max=8.58, description="Ar relative abundance", fixed=True) + + ca = Parameter(name="Ca", default=6.93, min=4.93, max=8.93, description="Ca relative abundance", fixed=True) + + fe = Parameter(name="Fe", default=8.1, min=6.1, max=10.1, description="Fe relative abundance", fixed=True) + + _input_units_allow_dimensionless = True + + def __init__( + self, + p=p.default, + break_energy=u.Quantity(break_energy.default, break_energy.unit), + q=q.default, + low_e_cutoff=u.Quantity(low_e_cutoff.default, low_e_cutoff.unit), + high_e_cutoff=u.Quantity(high_e_cutoff.default, high_e_cutoff.unit), + total_eflux=u.Quantity(total_eflux.default, total_eflux.unit), + plasma_density=u.Quantity(plasma_density.default, plasma_density.unit), + length=u.Quantity(length.default, length.unit), + temperature=u.Quantity(temperature.default, temperature.unit), + mg=mg.default, + al=al.default, + si=si.default, + s=s.default, + ar=ar.default, + ca=ca.default, + fe=fe.default, + abundance_type=DEFAULT_ABUNDANCE_TYPE, + integrator=None, + **kwargs, + ): + self.integrator = integrator + + # TODO: this version of ThickTarget still uses the scaled value with the non-scaled units + # so need to scale the value here with the inconsistent units + # Once fixed, replace ``total_eflux.to(u.electron/u.second)*1e-35`` with ``total_eflux`` + self.thick_model = ThickTarget(p=p, + break_energy=break_energy, + q=q, + low_e_cutoff=low_e_cutoff, + high_e_cutoff=high_e_cutoff, + total_eflux=total_eflux.to(u.electron/u.second)*1e-35, + integrator=integrator, + **kwargs + ) + + self.warm_component = ThickTargetWarmContribution(low_e_cutoff=low_e_cutoff, + total_eflux=total_eflux, + plasma_density=plasma_density, + length=length, + temperature=temperature, + mg=mg, + al=al, + si=si, + s=s, + ar=ar, + ca=ca, + fe=fe, + abundance_type=abundance_type, + **kwargs + ) + + super().__init__( + p=p, + break_energy=break_energy, + q=q, + low_e_cutoff=low_e_cutoff, + high_e_cutoff=high_e_cutoff, + total_eflux=total_eflux, + plasma_density=plasma_density, + length=length, + temperature=temperature, + mg=mg, + al=al, + si=si, + s=s, + ar=ar, + ca=ca, + fe=fe, + **kwargs, + ) + + def evaluate(self, + energy_edges, + p, + break_energy, + q, + low_e_cutoff, + high_e_cutoff, + total_eflux, + plasma_density, + length, + temperature, + mg, + al, + si, + s, + ar, + ca, + fe): + + energy_edges <<= 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.scaled_thick_eflux_units + plasma_density <<= self.scaled_warmthick_desnity_units + length <<= self.length.unit + temperature <<= self.temperature.unit + + thick = self.thick_model.evaluate(energy_edges, + p, + break_energy, + q, + low_e_cutoff, + high_e_cutoff, + total_eflux) + + warm = self.warm_component.evaluate(energy_edges, + low_e_cutoff, + total_eflux, + plasma_density, + length, + temperature, + mg, + al, + si, + s, + ar, + ca, + fe) + + print(thick, warm) + # TODO: this version of ThickTarget won't return units + # eventually, remove the line ``thick <<= warm.unit`` + thick <<= warm.unit + return thick + warm + + @property + def input_units(self): + # The units for the 'energy_edges' variable should be an energy (default keV) + return {self.inputs[0]: u.keV} + + @property + def return_units(self): + return {self.outputs[0]: u.ph * u.keV**-1 * u.s**-1} + + def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): + return { + "low_e_cutoff": u.keV, + "break_energy": u.keV, + "high_e_cutoff": u.keV, + "total_eflux": self.scaled_thick_eflux_units, + "temperature": u.MK, + "plasma_density":self.scaled_warmthick_desnity_units, + "legnth": u.Mm, + } \ No newline at end of file diff --git a/sunkit_spex/models/physical/thermal.py b/sunkit_spex/models/physical/thermal.py index 5aa99f13..0ed156d5 100644 --- a/sunkit_spex/models/physical/thermal.py +++ b/sunkit_spex/models/physical/thermal.py @@ -1317,8 +1317,8 @@ def _calculate_abundances(abundance_type, mg, al, si, s, ar, ca, fe): class ThickTargetWarmContribution(FittableModel): r""" Calculates the warm thick-target bremsstrahlung radiation contribution - of a single power-law electron distribution with a given low-energy - cut-off energy. + of a power-law electron distribution with a given low-energy cut-off + energy. [1] Kontar et al, ApJ 2015 (http://adsabs.harvard.edu/abs/2015arXiv150503733K) [2] https://hesperia.gsfc.nasa.gov/ssw/packages/xray/idl/f_thick_warm.pro @@ -1336,18 +1336,18 @@ class ThickTargetWarmContribution(FittableModel): low_e_cutoff : int or float Low-energy cut-off of the electron distribution in units of keV. - temperature: `astropy.units.Quantity` - The temperature of the plasma. - Can be scalar or 1D of any length. If not scalar, the flux for each temperature - will be calculated. The first dimension of the output flux will correspond - to temperature. - plasma_density: `astropy.units.Quantity` The number density o the plasma. length: `astropy.units.Quantity` The plasma column length. + temperature: `astropy.units.Quantity` + The temperature of the plasma. + Can be scalar or 1D of any length. If not scalar, the flux for each temperature + will be calculated. The first dimension of the output flux will correspond + to temperature. + Returns ------- A 1d array of warm component from thick-target bremsstrahlung radiation From bc77df3ba108b71f1bc51d63c343952e9c6525d1 Mon Sep 17 00:00:00 2001 From: Kristopher Cooper Date: Fri, 31 Jul 2026 09:46:34 -0500 Subject: [PATCH 4/5] pre-commit and changelog --- changelog/293.feature.rst | 1 + sunkit_spex/models/physical/nonthermal.py | 163 ++++++++-------- sunkit_spex/models/physical/thermal.py | 223 +++++++++++----------- 3 files changed, 192 insertions(+), 195 deletions(-) create mode 100644 changelog/293.feature.rst diff --git a/changelog/293.feature.rst b/changelog/293.feature.rst new file mode 100644 index 00000000..eb9677c5 --- /dev/null +++ b/changelog/293.feature.rst @@ -0,0 +1 @@ +The ``models.physical`` module now contains the physical models ``thermal.ThickTargetWarmContribution`` and ``nonthermal.WarmThickTarget``. diff --git a/sunkit_spex/models/physical/nonthermal.py b/sunkit_spex/models/physical/nonthermal.py index 5bad9479..9f6daf84 100644 --- a/sunkit_spex/models/physical/nonthermal.py +++ b/sunkit_spex/models/physical/nonthermal.py @@ -4,9 +4,11 @@ import astropy.units as u from astropy.modeling import FittableModel, Parameter +from astropy.modeling.functional_models import FLOAT_EPSILON from sunkit_spex.legacy import constants as const from sunkit_spex.legacy.integrate import gauss_legendre +from sunkit_spex.models.physical.thermal import DEFAULT_ABUNDANCE_TYPE, ThickTargetWarmContribution const = const.Constants() @@ -1239,11 +1241,10 @@ def bremsstrahlung_thick_target(photon_energies, p, break_energy, q, low_e_cutof raise Warning("The photon energies are higher than the highest electron energy or not greater than zero") -from astropy.modeling.functional_models import FLOAT_EPSILON -from sunkit_spex.models.physical.thermal import ThickTargetWarmContribution, DEFAULT_ABUNDANCE_TYPE + class WarmThickTarget(FittableModel): r""" - Calculates the thick-target + thermal contribution bremsstrahlung + Calculates the thick-target + thermal contribution bremsstrahlung radiation contribution. [1] Kontar et al, ApJ 2015 (http://adsabs.harvard.edu/abs/2015arXiv150503733K) @@ -1289,7 +1290,7 @@ class WarmThickTarget(FittableModel): Returns ------- - A 1d array of warm component from thick-target bremsstrahlung radiation + A 1d array of warm component from thick-target bremsstrahlung radiation in units of ph s^-1 keV^-1. """ @@ -1320,24 +1321,23 @@ class WarmThickTarget(FittableModel): ) plasma_density = Parameter( - name="plasma_density", - default=1, - unit=scaled_warmthick_desnity_units, - description="Number density of the plasma", - fixed=False, - bounds=(FLOAT_EPSILON, None) - ) + name="plasma_density", + default=1, + unit=scaled_warmthick_desnity_units, + description="Number density of the plasma", + fixed=False, + bounds=(FLOAT_EPSILON, None), + ) length = Parameter( - name="length", - default=10, - unit=u.Mm, - description="Plasma column length", - fixed=False, - bounds=(FLOAT_EPSILON, None) - ) + name="length", + default=10, + unit=u.Mm, + description="Plasma column length", + fixed=False, + bounds=(FLOAT_EPSILON, None), + ) - temperature = Parameter( name="temperature", default=10, @@ -1391,31 +1391,33 @@ def __init__( # TODO: this version of ThickTarget still uses the scaled value with the non-scaled units # so need to scale the value here with the inconsistent units # Once fixed, replace ``total_eflux.to(u.electron/u.second)*1e-35`` with ``total_eflux`` - self.thick_model = ThickTarget(p=p, - break_energy=break_energy, - q=q, - low_e_cutoff=low_e_cutoff, - high_e_cutoff=high_e_cutoff, - total_eflux=total_eflux.to(u.electron/u.second)*1e-35, - integrator=integrator, - **kwargs - ) - - self.warm_component = ThickTargetWarmContribution(low_e_cutoff=low_e_cutoff, - total_eflux=total_eflux, - plasma_density=plasma_density, - length=length, - temperature=temperature, - mg=mg, - al=al, - si=si, - s=s, - ar=ar, - ca=ca, - fe=fe, - abundance_type=abundance_type, - **kwargs - ) + self.thick_model = ThickTarget( + p=p, + break_energy=break_energy, + q=q, + low_e_cutoff=low_e_cutoff, + high_e_cutoff=high_e_cutoff, + total_eflux=total_eflux.to(u.electron / u.second) * 1e-35, + integrator=integrator, + **kwargs, + ) + + self.warm_component = ThickTargetWarmContribution( + low_e_cutoff=low_e_cutoff, + total_eflux=total_eflux, + plasma_density=plasma_density, + length=length, + temperature=temperature, + mg=mg, + al=al, + si=si, + s=s, + ar=ar, + ca=ca, + fe=fe, + abundance_type=abundance_type, + **kwargs, + ) super().__init__( p=p, @@ -1437,24 +1439,26 @@ def __init__( **kwargs, ) - def evaluate(self, - energy_edges, - p, - break_energy, - q, - low_e_cutoff, - high_e_cutoff, - total_eflux, - plasma_density, - length, - temperature, - mg, - al, - si, - s, - ar, - ca, - fe): + def evaluate( + self, + energy_edges, + p, + break_energy, + q, + low_e_cutoff, + high_e_cutoff, + total_eflux, + plasma_density, + length, + temperature, + mg, + al, + si, + s, + ar, + ca, + fe, + ): energy_edges <<= u.keV break_energy <<= self.break_energy.unit @@ -1465,29 +1469,12 @@ def evaluate(self, length <<= self.length.unit temperature <<= self.temperature.unit - thick = self.thick_model.evaluate(energy_edges, - p, - break_energy, - q, - low_e_cutoff, - high_e_cutoff, - total_eflux) - - warm = self.warm_component.evaluate(energy_edges, - low_e_cutoff, - total_eflux, - plasma_density, - length, - temperature, - mg, - al, - si, - s, - ar, - ca, - fe) - - print(thick, warm) + thick = self.thick_model.evaluate(energy_edges, p, break_energy, q, low_e_cutoff, high_e_cutoff, total_eflux) + + warm = self.warm_component.evaluate( + energy_edges, low_e_cutoff, total_eflux, plasma_density, length, temperature, mg, al, si, s, ar, ca, fe + ) + # TODO: this version of ThickTarget won't return units # eventually, remove the line ``thick <<= warm.unit`` thick <<= warm.unit @@ -1509,6 +1496,6 @@ def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): "high_e_cutoff": u.keV, "total_eflux": self.scaled_thick_eflux_units, "temperature": u.MK, - "plasma_density":self.scaled_warmthick_desnity_units, - "legnth": u.Mm, - } \ No newline at end of file + "plasma_density": self.scaled_warmthick_desnity_units, + "length": u.Mm, + } diff --git a/sunkit_spex/models/physical/thermal.py b/sunkit_spex/models/physical/thermal.py index 0ed156d5..d00c7d30 100644 --- a/sunkit_spex/models/physical/thermal.py +++ b/sunkit_spex/models/physical/thermal.py @@ -1313,11 +1313,10 @@ def _calculate_abundances(abundance_type, mg, al, si, s, ar, ca, fe): return abundances - class ThickTargetWarmContribution(FittableModel): r""" - Calculates the warm thick-target bremsstrahlung radiation contribution - of a power-law electron distribution with a given low-energy cut-off + Calculates the warm thick-target bremsstrahlung radiation contribution + of a power-law electron distribution with a given low-energy cut-off energy. [1] Kontar et al, ApJ 2015 (http://adsabs.harvard.edu/abs/2015arXiv150503733K) @@ -1350,7 +1349,7 @@ class ThickTargetWarmContribution(FittableModel): Returns ------- - A 1d array of warm component from thick-target bremsstrahlung radiation + A 1d array of warm component from thick-target bremsstrahlung radiation in units of ph s^-1 keV^-1. """ @@ -1381,24 +1380,23 @@ class ThickTargetWarmContribution(FittableModel): ) plasma_density = Parameter( - name="plasma_density", - default=1, - unit=scaled_warmthick_desnity_units, - description="Number density of the plasma", - fixed=False, - bounds=(FLOAT_EPSILON, None) - ) + name="plasma_density", + default=1, + unit=scaled_warmthick_desnity_units, + description="Number density of the plasma", + fixed=False, + bounds=(FLOAT_EPSILON, None), + ) length = Parameter( - name="length", - default=10, - unit=u.Mm, - description="Plasma column length", - fixed=False, - bounds=(FLOAT_EPSILON, None) - ) + name="length", + default=10, + unit=u.Mm, + description="Plasma column length", + fixed=False, + bounds=(FLOAT_EPSILON, None), + ) - temperature = Parameter( name="temperature", default=10, @@ -1426,107 +1424,121 @@ class ThickTargetWarmContribution(FittableModel): _input_units_allow_dimensionless = True def __init__( - self, - low_e_cutoff=u.Quantity(low_e_cutoff.default, low_e_cutoff.unit), - total_eflux=u.Quantity(total_eflux.default, total_eflux.unit), - plasma_density=u.Quantity(plasma_density.default, plasma_density.unit), - length=u.Quantity(length.default, length.unit), - temperature=u.Quantity(temperature.default, temperature.unit), - mg=mg.default, - al=al.default, - si=si.default, - s=s.default, - ar=ar.default, - ca=ca.default, - fe=fe.default, - abundance_type=DEFAULT_ABUNDANCE_TYPE, - **kwargs, - ): - + self, + low_e_cutoff=u.Quantity(low_e_cutoff.default, low_e_cutoff.unit), + total_eflux=u.Quantity(total_eflux.default, total_eflux.unit), + plasma_density=u.Quantity(plasma_density.default, plasma_density.unit), + length=u.Quantity(length.default, length.unit), + temperature=u.Quantity(temperature.default, temperature.unit), + mg=mg.default, + al=al.default, + si=si.default, + s=s.default, + ar=ar.default, + ca=ca.default, + fe=fe.default, + abundance_type=DEFAULT_ABUNDANCE_TYPE, + **kwargs, + ): + total_eflux <<= self.scaled_thick_eflux_units plasma_density <<= self.scaled_warmthick_desnity_units length <<= u.Mm self.c_cgs = const.c.cgs self.me_e = (const.m_e * const.c**2).to(u.keV) - self.mk2keV = const.k_B.to(u.keV/u.MK) + self.mk2keV = const.k_B.to(u.keV / u.MK) # collision parameter, units assigned by Kris to make the maths work later - self.col_param = 2.6e-18 << (u.cm**2 * u.keV**2) - - self.therm_model = ThermalEmission(temperature=temperature, - mg=mg, - al=al, - si=si, - s=s, - ar=ar, - ca=ca, - fe=fe, - abundance_type=abundance_type, - **kwargs, - ) + self.col_param = 2.6e-18 << (u.cm**2 * u.keV**2) + + self.therm_model = ThermalEmission( + temperature=temperature, + mg=mg, + al=al, + si=si, + s=s, + ar=ar, + ca=ca, + fe=fe, + abundance_type=abundance_type, + **kwargs, + ) super().__init__( - low_e_cutoff=low_e_cutoff, - total_eflux=total_eflux, - plasma_density=plasma_density, - length=length, - temperature=self.therm_model.temperature, - mg=self.therm_model.mg, - al=self.therm_model.al, - si=self.therm_model.si, - s=self.therm_model.s, - ar=self.therm_model.ar, - ca=self.therm_model.ca, - fe=self.therm_model.fe, - **kwargs, - ) + low_e_cutoff=low_e_cutoff, + total_eflux=total_eflux, + plasma_density=plasma_density, + length=length, + temperature=self.therm_model.temperature, + mg=self.therm_model.mg, + al=self.therm_model.al, + si=self.therm_model.si, + s=self.therm_model.s, + ar=self.therm_model.ar, + ca=self.therm_model.ca, + fe=self.therm_model.fe, + **kwargs, + ) def calculate_emission_measure(self, temperature, plasma_density, length, low_e_cutoff, total_eflux): """ Calculate the emission measure of the warm thick target contribution. - + Need some maths and units magic here. """ # convert units - loop_temp_kev = (temperature*self.mk2keV).to(u.keV) + loop_temp_kev = (temperature * self.mk2keV).to(u.keV) plasma_density_cm = plasma_density.to(u.cm**-3) length_cm = length.to(u.cm) - eflux = total_eflux.to(u.electron/u.second) + eflux = total_eflux.to(u.electron / u.second) - col_stop_dist = loop_temp_kev**2 / (2 * self.col_param * plasma_density_cm) # collisional stopping distance for electrons of Tloop energy + col_stop_dist = loop_temp_kev**2 / ( + 2 * self.col_param * plasma_density_cm + ) # collisional stopping distance for electrons of Tloop energy emin = loop_temp_kev * 3 * (5 * col_stop_dist / length_cm) ** 4 - if emin > 0.1< 0.1 << u.keV: warnings.warn( f"The loop temperature ({temperature}), plasma number density ({plasma_density}), and loop length ({length}) make emin ({emin}) >0.1 keV. Fixing emin to 0.1." - ) - emin = 0.1< length_cm: warnings.warn(f"Minimum length ({lmin.to(u.Mm)}) > loop length ({length})") - emission_measure = 3 * np.pi / 2 / self.col_param / self.c_cgs * np.sqrt(self.me_e / 8.0) * loop_temp_kev**2 / np.sqrt(emin) * eflux - emission_measure /= u.electron # electron interactions but let's remove the unit here + emission_measure = ( + 3 + * np.pi + / 2 + / self.col_param + / self.c_cgs + * np.sqrt(self.me_e / 8.0) + * loop_temp_kev**2 + / np.sqrt(emin) + * eflux + ) + emission_measure /= u.electron # electron interactions but let's remove the unit here return emission_measure << self.scaled_em_units - def evaluate(self, - energy_edges, - low_e_cutoff, - total_eflux, - plasma_density, - length, - temperature, - mg, - al, - si, - s, - ar, - ca, - fe, - ): + def evaluate( + self, + energy_edges, + low_e_cutoff, + total_eflux, + plasma_density, + length, + temperature, + mg, + al, + si, + s, + ar, + ca, + fe, + ): low_e_cutoff <<= self.low_e_cutoff.unit total_eflux <<= self.scaled_thick_eflux_units @@ -1535,27 +1547,24 @@ def evaluate(self, length <<= self.length.unit emission_measure = self.calculate_emission_measure( - temperature, - plasma_density, - length, - low_e_cutoff, - total_eflux - ) + temperature, plasma_density, length, low_e_cutoff, total_eflux + ) # TODO: this version of Thermal Emission still uses the scaled value with the non-scaled units # so need to scale the value here with the inconsistent units # Once ixed, replace ``emission_measure.value< Date: Wed, 12 Aug 2026 13:14:09 -0500 Subject: [PATCH 5/5] tests --- .../models/physical/tests/test_nonthermal.py | 82 ++++++++++++++++++- .../models/physical/tests/test_thermal.py | 67 +++++++++++++++ 2 files changed, 148 insertions(+), 1 deletion(-) diff --git a/sunkit_spex/models/physical/tests/test_nonthermal.py b/sunkit_spex/models/physical/tests/test_nonthermal.py index ac5419f8..7d2a9b43 100644 --- a/sunkit_spex/models/physical/tests/test_nonthermal.py +++ b/sunkit_spex/models/physical/tests/test_nonthermal.py @@ -3,7 +3,7 @@ import astropy.units as u -from sunkit_spex.models.physical import nonthermal +from sunkit_spex.models.physical import nonthermal, thermal SSW_INTENSITY_UNIT = u.ph / u.cm**2 / u.s / u.keV @@ -128,6 +128,71 @@ def thin_target(): return energy_edges, ssw_output +def warm_thick_target(): + """ + Defines an output for the ``WarmThickTarget`` model to be tested against. + """ + energy_edges = np.arange(1.6, 15, 0.1) << u.keV + p = 3 + break_energy = 8 << u.keV + q = 9 + low_e_cutoff = 6.5 << u.keV + high_e_cutoff = 20 << u.keV + total_eflux = 1.8e35 << (u.electron / u.second) + plasma_density = 6e9 << u.cm**-3 + length = 15 << u.Mm + temperature = 10.2 << u.MK + abundance_type = thermal.DEFAULT_ABUNDANCE_TYPE + # fmt: off + mg, al, si, s, ar, ca, fe = 8.15, 7.04, 8.1, 7.27, 6.58, 6.93, 8.4010299956639812 + + inputs = ( + energy_edges, p, break_energy, q, low_e_cutoff, high_e_cutoff, total_eflux, plasma_density, length, temperature, mg, al, si, s, ar, ca, fe, abundance_type + ) + inputs_class = ( + p, break_energy, q, low_e_cutoff, high_e_cutoff, total_eflux, plasma_density, length, temperature, mg, al, si, s, ar, ca, fe, abundance_type + ) + + sunkit_spex_output = [6.98739369e+30, 1.01719200e+31, 3.68431997e+31, 1.17447066e+31, + 9.12468731e+30, 4.81844542e+30, 3.73891273e+30, 3.59670674e+30, + 5.55001931e+30, 2.18283128e+30, 1.59862179e+30, 1.34409593e+30, + 1.32657257e+30, 1.11334927e+30, 1.04347888e+30, 1.03138066e+30, + 6.94568496e+29, 6.11863342e+29, 5.34767348e+29, 4.75676980e+29, + 4.31267036e+29, 3.74106058e+29, 5.32413715e+29, 3.95184372e+29, + 2.59124198e+29, 2.31076149e+29, 2.05884602e+29, 1.83867478e+29, + 1.68376672e+29, 1.57670792e+29, 1.34947945e+29, 1.23144678e+29, + 1.09685712e+29, 9.62801755e+28, 8.61523188e+28, 7.77996397e+28, + 7.01093694e+28, 6.32097056e+28, 5.70100880e+28, 5.14284911e+28, + 4.63963309e+28, 4.18510744e+28, 3.77402867e+28, 3.40166211e+28, + 3.06399500e+28, 2.75742346e+28, 2.48309862e+28, 2.31939862e+28, + 3.57553130e+28, 5.83649787e+28, 3.18994718e+28, 1.50731140e+28, + 1.27731459e+28, 1.14202129e+28, 1.02083863e+28, 9.13630768e+27, + 8.16290062e+27, 7.28365497e+27, 6.50723666e+27, 5.81378169e+27, + 5.19611558e+27, 4.75909108e+27, 4.24162300e+27, 3.72071796e+27, + 3.32832621e+27, 3.02096699e+27, 2.69705740e+27, 2.42628793e+27, + 2.16854906e+27, 1.94839050e+27, 1.75521895e+27, 1.58304562e+27, + 1.42931365e+27, 1.29156780e+27, 1.16816563e+27, 1.05752595e+27, + 9.58206892e+26, 8.68987528e+26, 7.88769191e+26, 7.16554939e+26, + 6.51520166e+26, 5.92859787e+26, 5.39939208e+26, 4.92119876e+26, + 4.48900219e+26, 4.09778810e+26, 3.74358794e+26, 3.42240662e+26, + 3.13113291e+26, 2.86656812e+26, 2.62620732e+26, 2.40755579e+26, + 2.20854639e+26, 2.02725453e+26, 1.86194774e+26, 1.71115861e+26, + 1.57343681e+26, 1.44760910e+26, 1.33253400e+26, 1.22720708e+26, + 1.13077112e+26, 1.04237417e+26, 9.61307568e+25, 8.86919646e+25, + 8.18597166e+25, 7.55821022e+25, 6.98101580e+25, 6.44990152e+25, + 5.96101034e+25, 5.51069826e+25, 5.09562885e+25, 4.71289794e+25, + 4.35981392e+25, 4.03385834e+25, 3.73281810e+25, 3.45471638e+25, + 3.19762661e+25, 2.95986396e+25, 2.73992361e+25, 2.53637342e+25, + 2.34790933e+25, 2.17335568e+25, 2.01166829e+25, 1.86182095e+25, + 1.72290963e+25, 1.59410762e+25, 1.47466772e+25, 1.36386829e+25, + 1.26106646e+25, 1.16567102e+25, 1.07714863e+25, 9.94983441e+24, + 9.18713405e+24] * (u.ph / (u.keV * u.s)) + return inputs, inputs_class, energy_edges, sunkit_spex_output + + +warm_thick_target() + + @pytest.mark.parametrize("ssw", [thick_target]) def test_thick_target_against_ssw(ssw): energy_edges, expected = ssw() @@ -144,3 +209,18 @@ 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) + + +@pytest.mark.parametrize("sunkit_spex", [warm_thick_target]) +def test_thick_target_against_previous(sunkit_spex): + _, input_args_class, energy_edges, expected = sunkit_spex() + model_class = nonthermal.WarmThickTarget(*input_args_class) + output_class = model_class(energy_edges) + expected_value = expected.to_value(output_class.unit) + # check direct output + np.testing.assert_allclose(output_class.value, expected_value, rtol=0.05, atol=1e-30) + # now check that the warm thick target is warm+thick target + wc_model = thermal.ThickTargetWarmContribution(*(input_args_class[3], *input_args_class[5:])) + tt_model = nonthermal.ThickTarget(*input_args_class[:6]) + output_comb = (wc_model + tt_model)(energy_edges) + np.testing.assert_allclose(output_comb.value, expected_value, rtol=0.035) diff --git a/sunkit_spex/models/physical/tests/test_thermal.py b/sunkit_spex/models/physical/tests/test_thermal.py index 89fb2c53..df2a56fc 100644 --- a/sunkit_spex/models/physical/tests/test_thermal.py +++ b/sunkit_spex/models/physical/tests/test_thermal.py @@ -386,6 +386,64 @@ def chianti_kev_lines_Fe2(): return inputs, inputs_class, energy_edges, ssw_output +def thick_target_warm_contribution(): + """ + Defines an output for the ``ThickTargetWarmContribution`` model + to be tested against. + """ + energy_edges = np.arange(1.6, 15, 0.1) << u.keV + temperature = 10.2 << u.MK + plasma_density = 6e9 << u.cm**-3 + low_e_cutoff = 6.5 << u.keV + total_eflux = 1.8e35 << (u.electron / u.second) + length = 15 << u.Mm + abundance_type = DEFAULT_ABUNDANCE_TYPE + # fmt: off + mg, al, si, s, ar, ca, fe = 8.15, 7.04, 8.1, 7.27, 6.58, 6.93, 8.4010299956639812 + + inputs = ( + energy_edges, low_e_cutoff, total_eflux, plasma_density, length, temperature, mg, al, si, s, ar, ca, fe, abundance_type + ) + inputs_class = ( + low_e_cutoff, total_eflux, plasma_density, length, temperature, mg, al, si, s, ar, ca, fe, abundance_type + ) + sunkit_spex_output = [6.44949131e+30, 9.68510522e+30, 3.64009470e+31, 1.13415879e+31, + 8.75614449e+30, 4.48061840e+30, 3.42850740e+30, 3.31089228e+30, + 5.28634751e+30, 1.93917241e+30, 1.37311322e+30, 1.13510049e+30, + 1.13264467e+30, 9.33206989e+29, 8.75980923e+29, 8.75507150e+29, + 5.49404869e+29, 4.76586797e+29, 4.08635361e+29, 3.58017526e+29, + 3.21470232e+29, 2.71616979e+29, 4.36726197e+29, 3.05835653e+29, + 1.75690250e+29, 1.53167611e+29, 1.33143219e+29, 1.15962972e+29, + 1.05003990e+29, 9.85476704e+28, 7.98127510e+28, 7.17544932e+28, + 6.18146249e+28, 5.17177540e+28, 4.47022426e+28, 3.92784720e+28, + 3.43454494e+28, 3.00421486e+28, 2.62878786e+28, 2.30096583e+28, + 2.01471949e+28, 1.76455555e+28, 1.54592668e+28, 1.35473494e+28, + 1.18754844e+28, 1.04129031e+28, 9.17588886e+27, 8.95244823e+27, + 2.28384164e+28, 4.66833951e+28, 2.13455961e+28, 5.54442739e+27, + 4.17576618e+27, 3.66814685e+27, 3.22285512e+27, 2.84513673e+27, + 2.49996801e+27, 2.18839217e+27, 1.92417640e+27, 1.69214868e+27, + 1.48945428e+27, 1.42493811e+27, 1.24119455e+27, 1.01864256e+27, + 8.92414274e+26, 8.22404648e+26, 7.10419441e+26, 6.29127713e+26, + 5.41019869e+26, 4.72918762e+26, 4.16221692e+26, 3.66685749e+26, + 3.23287648e+26, 2.84921249e+26, 2.51132595e+26, 2.21389061e+26, + 1.95175083e+26, 1.72089201e+26, 1.51758643e+26, 1.33830307e+26, + 1.18051372e+26, 1.04125183e+26, 9.18697822e+25, 8.10500671e+25, + 7.15244883e+25, 6.31132611e+25, 5.57072025e+25, 4.91643710e+25, + 4.34045128e+25, 3.83149052e+25, 3.38311163e+25, 2.98709602e+25, + 2.63782790e+25, 2.32962460e+25, 2.05739148e+25, 1.81749313e+25, + 1.60539864e+25, 1.41839259e+25, 1.25322066e+25, 1.10723052e+25, + 9.78551010e+24, 8.64742902e+24, 7.64240622e+24, 6.75545270e+24, + 5.97088163e+24, 5.27832255e+24, 4.66665193e+24, 4.12549626e+24, + 3.64775453e+24, 3.22571621e+24, 2.85226405e+24, 2.52235232e+24, + 2.23101750e+24, 1.97317226e+24, 1.74509605e+24, 1.54389741e+24, + 1.36579076e+24, 1.20813628e+24, 1.06893488e+24, 9.45857108e+23, + 8.36888795e+23, 7.40418659e+23, 6.55344272e+23, 5.80006514e+23, + 5.13293421e+23, 4.54249083e+23, 4.02146850e+23, 3.55997428e+23, + 3.15123042e+23, 2.78926178e+23, 2.46986782e+23, 2.18691162e+23, + 1.93624991e+23] * (u.ph / (u.keV * u.s)) + return inputs, inputs_class, energy_edges, sunkit_spex_output + + @pytest.mark.parametrize("ssw", [fvth_simple, fvth_Fe2]) def test_thermal_emission_against_ssw(ssw): _, input_args_class, energy_edges, expected = ssw() @@ -513,3 +571,12 @@ def test_abundances_should_not_change(): after_models = thermal.DEFAULT_ABUNDANCES[thermal.DEFAULT_ABUNDANCE_TYPE].data assert np.allclose(after_models.data, orig.data) + + +@pytest.mark.parametrize("sunkit_spex", [thick_target_warm_contribution]) +def test_thick_target_warm_contribution_against_previous(sunkit_spex): + _, input_args_class, energy_edges, expected = sunkit_spex() + model_class = thermal.ThickTargetWarmContribution(*input_args_class) + output_class = model_class(energy_edges) + expected_value = expected.to_value(output_class.unit) + np.testing.assert_allclose(output_class.value, expected_value, rtol=0.05, atol=1e-30)