Add warm thick

changelog/293.feature.rst

@@ -0,0 +1 @@
+The ``models.physical`` module now contains the physical models ``thermal.ThickTargetWarmContribution`` and ``nonthermal.WarmThickTarget``.

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,
+        }

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)