Skip to content

Thick Target (and thermal) Update #291

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
KriSun95 wants to merge 10 commits into
base: main
Choose a base branch
from
Open

Thick Target (and thermal) Update #291

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
cd3a615
Trying to fix scaling of models
KriSun95 Jul 27, 2026
a949b0f
Units and norm scaling is fixed in the thermal, thick, thin models
KriSun95 Jul 27, 2026
e0972f9
Found reason why the thick target model would return NaNs
KriSun95 Jul 28, 2026
1f57e2e
pre-commit
KriSun95 Jul 28, 2026
9df211e
changelog added
KriSun95 Jul 29, 2026
ceb0165
Made the new scaled units accessable to the Astropy methods like .fin…
KriSun95 Jul 29, 2026
4a0c364
pre-commit
KriSun95 Jul 29, 2026
db1165f
Removed accidental call to test function I was checking.
KriSun95 Jul 29, 2026
cd8d2ec
use <<= to assign or convert albedo theta input to degrees in evaluat…
KriSun95 Jul 29, 2026
8c682a4
Merge branch 'sunpy:main' into tt-fix
KriSun95 Jul 29, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions changelog/291.feature.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
The scalar parameters for ``ThermalEmission``, ``ThickTarget``, and ``ThinTarget`` in ``models.physical`` now use consistent custom scaled values and units instead of scaled values with non-scaled units. Parameters for the ``ThermalEmission`` and ``ThickTarget`` are given suitable parameter bounds.
4 changes: 2 additions & 2 deletions sunkit_spex/models/physical/albedo.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -94,8 +94,8 @@ def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)

def evaluate(self, spectrum, theta, anisotropy):
if not isinstance(theta, Quantity):
theta = theta * u.deg

theta <<= u.deg

albedo_matrix = get_albedo_matrix(self.energy_edges, theta, anisotropy)

Expand Down
229 changes: 70 additions & 159 deletions sunkit_spex/models/physical/nonthermal.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.scaling import scaled_thick_eflux_units, scaled_thin_eflux_units

const = const.Constants()

Expand All @@ -29,6 +31,8 @@

__all__ = ["ThickTarget", "ThinTarget"]

FLOAT_EPSILON_FOR_POWER_LAW = 1 + FLOAT_EPSILON * 1e30

@samaloney samaloney Aug 6, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Where does the FLOAT_EPSILON * 1e30 come from, trail and error or something else?

@KriSun95 KriSun95 Aug 7, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Trial and error unfortunately. I needed a value as close to 1 as I could get but any of the Python epsilons added to 1 just returned 1 annoyingly. The FLOAT_EPSILON * 1e30 was the lowest value I could find that would return 1+(a little something).

I.e., 1 - (1 + FLOAT_EPSILON)==0 when evaluated in the model but 1 - (1 + FLOAT_EPSILON * 1e30)==FLOAT_EPSILON * 1e30.

I'm more than open to suggestions on using something better but I couldn't find anything :/



class ThickTarget(FittableModel):
r"""Calculates the thick-target bremsstrahlung radiation of a dual power-law electron distribution.
Expand Down Expand Up @@ -73,22 +77,48 @@ class ThickTarget(FittableModel):
n_inputs = 1
n_outputs = 1

p = Parameter(name="p", default=2, description="Slope below break", fixed=False)
p = Parameter(
name="p", default=2, description="Slope below break", fixed=False, bounds=(FLOAT_EPSILON_FOR_POWER_LAW, None)
)

break_energy = Parameter(name="break_energy", default=100, unit=u.keV, description="Break Energy", fixed=False)
break_energy = Parameter(
name="break_energy",
default=100,
unit=u.keV,
description="Break Energy",
fixed=True,
bounds=(FLOAT_EPSILON, None),
)

q = Parameter(name="q", default=5, min=0.01, description="Slope above break", fixed=True)
q = Parameter(
name="q", default=5, description="Slope above break", fixed=True, 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
name="low_e_cutoff",
default=7,
unit=u.keV,
description="Low energy electron cut off",
fixed=False,
bounds=(FLOAT_EPSILON, None),
)

high_e_cutoff = Parameter(
name="high_e_cutoff", default=1500, unit=u.keV, description="High energy electron cut off", fixed=True
name="high_e_cutoff",
default=1500,
unit=u.keV,
description="High energy electron cut off",
fixed=True,
bounds=(FLOAT_EPSILON, None),
)

total_eflux = Parameter(
name="total_eflux", default=1.5, unit=u.electron * u.s**-1, description="Total electron flux", fixed=True
name="total_eflux",
default=1.5,
unit=scaled_thick_eflux_units,
description="Total electron flux",
fixed=False,
bounds=(0, None),
)

_input_units_allow_dimensionless = True
Expand All @@ -106,6 +136,8 @@ def __init__(
):
self.integrator = integrator

total_eflux <<= scaled_thick_eflux_units

super().__init__(
p=p,
break_energy=break_energy,
Expand All @@ -119,29 +151,26 @@ def __init__(
def evaluate(self, energy_edges, p, break_energy, q, low_e_cutoff, high_e_cutoff, total_eflux):
energy_centers = energy_edges[:-1] + 0.5 * np.diff(energy_edges)

if (
hasattr(break_energy, "unit")
or hasattr(energy_centers, "unit")
or hasattr(low_e_cutoff, "unit")
or hasattr(high_e_cutoff, "unit")
or hasattr(total_eflux, "unit")
):
flux = thick_fn(
energy_centers <<= 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 <<= scaled_thick_eflux_units

flux = (
bremsstrahlung_thick_target(
energy_centers.value,
p,
break_energy.value,
q,
low_e_cutoff.value,
high_e_cutoff.value,
total_eflux.value,
self.integrator,
)
else:
flux = thick_fn(
energy_centers, p, break_energy, q, low_e_cutoff, high_e_cutoff, total_eflux, self.integrator
)
* total_eflux.decompose().value
)

return flux
return flux * self.return_units[self.outputs[0]]

@property
def input_units(self):
Expand All @@ -157,7 +186,7 @@ def _parameter_units_for_data_units(self, inputs_unit, outputs_unit):
"break_energy": u.keV,
"low_e_cutoff": u.keV,
"high_e_cutoff": u.keV,
"total_eflux": u.electron * u.s**-1,
"total_eflux": scaled_thick_eflux_units,
}


Expand Down Expand Up @@ -225,7 +254,11 @@ class ThinTarget(FittableModel):
)

total_eflux = Parameter(
name="total_eflux", default=1.5, unit=u.s**-1 * u.cm**-2, description="Total electron flux", fixed=True
name="total_eflux",
default=1.5,
unit=scaled_thin_eflux_units,
description="Total electron flux",
fixed=True,
)

_input_units_allow_dimensionless = True
Expand All @@ -243,6 +276,8 @@ def __init__(
):
self.integrator = integrator

total_eflux <<= scaled_thin_eflux_units

super().__init__(
p=p,
break_energy=break_energy,
Expand All @@ -256,29 +291,26 @@ def __init__(
def evaluate(self, energy_edges, p, break_energy, q, low_e_cutoff, high_e_cutoff, total_eflux):
energy_centers = energy_edges[:-1] + 0.5 * np.diff(energy_edges)

if (
hasattr(break_energy, "unit")
or hasattr(energy_centers, "unit")
or hasattr(low_e_cutoff, "unit")
or hasattr(high_e_cutoff, "unit")
or hasattr(total_eflux, "unit")
):
flux = thin_fn(
energy_centers <<= 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 <<= scaled_thin_eflux_units

flux = (
bremsstrahlung_thin_target(
energy_centers.value,
p,
break_energy.value,
q,
low_e_cutoff.value,
high_e_cutoff.value,
total_eflux.value,
self.integrator,
)
else:
flux = thin_fn(
energy_centers, p, break_energy, q, low_e_cutoff, high_e_cutoff, total_eflux, self.integrator
integrator=self.integrator,
)
* total_eflux.to(u.electron * u.cm ** (-2) * u.s**-1).value
)

return flux
return flux * self.return_units[self.outputs[0]]

@property
def input_units(self):
Expand All @@ -294,132 +326,11 @@ def _parameter_units_for_data_units(self, inputs_unit, outputs_unit):
"break_energy": u.keV,
"low_e_cutoff": u.keV,
"high_e_cutoff": u.keV,
"total_eflux": u.s**-1 * u.cm**-2,
"total_eflux": scaled_thin_eflux_units,
# "total_eflux": u.electron * u.s**-1,
}


def thick_fn(energy_centers, p, break_energy, q, low_e_cutoff, high_e_cutoff, total_eflux, integrator):
"""Calculates the thick-target bremsstrahlung radiation of a dual power-law electron distribution.

[1] Brown, Solar Physics 18, 489 (1971) (https://link.springer.com/article/10.1007/BF00149070)
[2] https://hesperia.gsfc.nasa.gov/ssw/packages/xray/doc/brm_thick_doc.pdf
[3] https://hesperia.gsfc.nasa.gov/ssw/packages/xray/idl/brm2/brm2_thicktarget.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.
"""

# hack = np.round([p, break_energy, q, low_e_cutoff, high_e_cutoff, total_eflux], 15)
# p, break_energy, q, low_e_cutoff, high_e_cutoff, total_eflux = hack[0], hack[1], hack[2], hack[3], hack[4], hack[5]

# energies = np.mean(energies, axis=1) # since energy bins are given, use midpoints though

# we want a single power law electron distribution,
# so set break_energy == high_e_cutoff at a high value.
# we don't care about q at E > break_energy.
# high_break = energies.max() * 10

output = bremsstrahlung_thick_target(energy_centers, p, break_energy, q, low_e_cutoff, high_e_cutoff, integrator)

output[np.isnan(output)] = 0
output[~np.isfinite(output)] = 0

# convert to 1e35 e-/s
return output * total_eflux * 1e35


# def thin_fn(total_eflux, index, e_c, energies=None):
def thin_fn(energy_centers, p, break_energy, q, low_e_cutoff, high_e_cutoff, total_eflux, integrator):
"""Calculates the thin-target bremsstrahlung radiation of a dual power-law electron distribution.

[1] Brown, Solar Physics 18, 489 (1971) (https://link.springer.com/article/10.1007/BF00149070)
[2] https://hesperia.gsfc.nasa.gov/ssw/packages/xray/doc/brm_thick_doc.pdf
[3] https://hesperia.gsfc.nasa.gov/ssw/packages/xray/idl/brm2/brm2_thicktarget.pro

Parameters
----------
energy_edges : 1d array
Edges of energy bins in units of keV.

total_eflux : int or float
normalization factor in units of 1.0d55 cm-2 sec-1,
i.e. plasma density * volume of source * integrated nonthermal electron flux density
Need to take care here as the model returns units of cm-2 sec-1 as the scaling factor of 1e55 is hidden.
So actual units are 1.0d55 cm-2 sec-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 thin-target bremsstrahlung radiation in units
of ph s^-1 keV^-1.
"""

# hack = np.round([total_eflux, index, e_c], 15)
# total_eflux, index, e_c = hack[0], hack[1], hack[2]

# energies = np.mean(energies, axis=1) # since energy bins are given, use midpoints though
# energies = energy_centers
# we want a single power law electron distribution,
# so set break_energy == high_e_cutoff at a high value.
# we don't care about q at E > break_energy.
# high_break = energies.max() * 10
output = bremsstrahlung_thin_target(
energy_centers, p, break_energy, q, low_e_cutoff, high_e_cutoff, total_eflux, integrator
)

output[np.isnan(output)] = 0
output[~np.isfinite(output)] = 0

# convert to 1e35 e-/s
return output * total_eflux * 1e55


class BrokenPowerLawElectronDistribution:
"""
A broken or double power law electron flux distribution and integral.
Expand Down
Loading
Loading