Add warm thick

[KriSun95]

Opening this WIP PR that introduces the warm thick target model...at least the core part of it

The warm thick target is the addition of the thick target and a thermal component; that thermal component is the ThermalEmission model where the emission measure input is calculated from the thick target's low-energy cut-off and electron flux values, and new plasma number density and column length parameters.

I've made a new model called ThickTargetWarmContribution that, given these parameters, will calculate the thermal component o the whole warm thick target model we know rom IDL. That should mean that to re create the whole warm thick target model we would only need to do:

warm_thick = ThickTarget()+ThickTargetWarmContribution()

and make sure the low-energy cut-off and electron flux parameters are shared.

It is a work in progress, still need to add tests and comparisons to IDL, etc. but we'll want this physical model in eventually I'd imagine so might as well put it here for folk to see.

Here is a snippet of code I've been using it in:

import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np

from sunkit_spex.models.physical.nonthermal import ThickTarget
from sunkit_spex.models.physical.thermal import ThermalEmission,  ThickTargetWarmContribution

energy_edges = np.arange(1.6, 15, 0.1) << u.keV

energy_mids = (energy_edges[:-1]+energy_edges[1:])/2

## parameters from https://iopscience.iop.org/article/10.3847/2041-8213/ab7341 Table 1 Lindsay's paper
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) 
p = 6.3
length = 15 << u.Mm

model_thick = ThickTarget(total_eflux=total_eflux/1e35, p=p, low_e_cutoff=low_e_cutoff, break_energy=1500<< u.keV, q=20, high_e_cutoff=1500<< u.keV)
model_warmtherm2 =  ThickTargetWarmContribution(temperature=temperature,
                                                plasma_density=plasma_density,
                                                low_e_cutoff=low_e_cutoff,
                                                total_eflux=total_eflux,
                                                length=length,
                                                )
# print(model_warmtherm2.parameters)
flux_c2 = model_warmtherm2(energy_edges)
# print(model_warmtherm2.parameters)
flux_e2 = model_warmtherm2.evaluate(energy_edges, *model_warmtherm2.parameters)

plt.figure()
plt.plot(energy_mids, flux_c2, marker="o", label="warm called")
plt.plot(energy_mids, flux_e2, marker="x", ls=":", label="warm evaluated")
plt.plot(energy_mids, model_thick(energy_edges), marker="o", label="thick called")
plt.plot(energy_mids, model_thick(energy_edges)+flux_c2, marker="o", label="warm thick total")
plt.legend()
plt.yscale("log")
plt.ylabel(f"[{flux_c2.unit:latex}]")
plt.xlabel(f"Energy [{energy_mids.unit:latex}]")
plt.show()

[settwi]

Related to #235

[KriSun95]

Just another small addition. I've now made a model that combines ThickTarget and ThickTargetWarmContribution.

Still to add:

• Tests
• Examples

Here is a snippet of code and the plot showing it working:

import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np

from sunkit_spex.models.physical.nonthermal import ThickTarget, WarmThickTarget
from sunkit_spex.models.physical.thermal import ThermalEmission,  ThickTargetWarmContribution

energy_edges = np.arange(1.6, 15, 0.1) << u.keV

energy_mids = (energy_edges[:-1]+energy_edges[1:])/2

## parameters from https://iopscience.iop.org/article/10.3847/2041-8213/ab7341 Table 1 Lindsay's paper
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) 
p = 6.3
length = 15 << u.Mm

model_thick = ThickTarget(total_eflux=total_eflux/1e35, p=p, low_e_cutoff=low_e_cutoff, break_energy=1500<< u.keV, q=20, high_e_cutoff=1500<< u.keV)
model_warmtherm2 =  ThickTargetWarmContribution(temperature=temperature,
                                                plasma_density=plasma_density,
                                                low_e_cutoff=low_e_cutoff,
                                                total_eflux=total_eflux,
                                                length=length,
                                                )
# print(model_warmtherm2.parameters)
flux_c2 = model_warmtherm2(energy_edges)
# print(model_warmtherm2.parameters)
flux_e2 = model_warmtherm2.evaluate(energy_edges, *model_warmtherm2.parameters)

model_awt = WarmThickTarget(total_eflux=total_eflux, 
                            p=p, 
                            low_e_cutoff=low_e_cutoff, 
                            break_energy=1500<< u.keV, 
                            q=20, 
                            high_e_cutoff=1500<< u.keV,
                            temperature=temperature,
                            plasma_density=plasma_density,
                            length=length
                            )
flux_c3 = model_awt(energy_edges)
flux_e3 = model_awt.evaluate(energy_edges, *model_awt.parameters)

plt.figure()
plt.plot(energy_mids, flux_c2, marker="o", label="warm called")
plt.plot(energy_mids, flux_e2, marker="x", ls=":", label="warm evaluated")
plt.plot(energy_mids, model_thick(energy_edges), marker="o", label="thick called")
plt.plot(energy_mids, model_thick(energy_edges)+flux_c2, marker="o", label="warm thick total")
plt.plot(energy_mids, flux_c3, marker="o", ms=3, label="WARM THICK called")
plt.plot(energy_mids, flux_e3, marker="o", ms=2, label="WARM THICK eval")
plt.legend()
plt.yscale("log")
plt.ylabel(f"[{flux_c2.unit:latex}]")
plt.xlabel(f"Energy [{energy_mids.unit:latex}]")
plt.show()

[samaloney]

There must be a better way to do this without repeating the same parameter definition everywhere