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

[samaloney]

Yea see #297 for potential method to not have to repeat the parameters would need to be tweaked a bit for here but similar approach should work.

[KriSun95]

Thank you so much for looking into this. I tried to get around it using inheritance but I kept running into meta-class stuff. This was so annoying to write so I'll edit using #297

[KriSun95]

@samaloney , looking at #297 , I think it would be better to merge that first then come back to this PR. The WarmThickTarget and ThickTargetWarmContribution models need all the thermal parameters anyway so I'd need to duplicate the work in #297 if it wasn't merged first.

@ianan , I'm also curious about your thoughts on how I've split up the warm thick target model, effectively just separating out the thermal contribution into its own separate model so we have both the ThickTargetWarmContribution mode and overall WarmThickTarget (equal to ThickTarget+ThickTargetWarmContribution) model. Additionally, any comments on the descriptions in the doc-strings (here and here ) or anything else would be greatly appreciated.

I'll move this off of a draft PR since I don't know if it being a "draft" stops folk editing/collaborating on it.