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 1e95e29a..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() @@ -1238,3 +1240,262 @@ 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") + + +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 + ) + + # 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, + "length": u.Mm, + } 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) diff --git a/sunkit_spex/models/physical/thermal.py b/sunkit_spex/models/physical/thermal.py index 73d09897..d00c7d30 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,275 @@ 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 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. + + 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 = "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 << 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 << u.keV + + lmin = low_e_cutoff**2 / (2 * self.col_param * plasma_density_cm) / 3 + if lmin > 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<