Skip to content

Bugfix/filter sed conversion #80

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

Merged
PabloCorcho merged 4 commits into from
Oct 10, 2025
Merged

Bugfix/filter sed conversion #80

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
e7bb4cc
filter accepts fnu and flambda
PabloCorcho Oct 10, 2025
dd049c4
test fnu and flam
PabloCorcho Oct 10, 2025
b9c9887
check wavelength exists
PabloCorcho Oct 10, 2025
5657583
update check_unit
PabloCorcho Oct 10, 2025
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
45 changes: 42 additions & 3 deletions src/pst/observables.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -292,10 +292,49 @@ def interpolate(self, wavelength=None):
return self.response

def _check_spectra(self, spectra, default_unit=u.Lsun / u.angstrom / u.cm**2):
Copy link
Copy Markdown
Owner

@paranoya paranoya Oct 10, 2025 •
edited
Loading

Choose a reason for hiding this comment

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

Perhaps _ensure_f_lambda, _ensure_f_lambda_units, or ensure_spectra_units would be more descriptive... I think I lean towards the second (perhaps renaming spectra as f_lambda?)

Copy link
Copy Markdown
Owner

@paranoya paranoya Oct 10, 2025

Choose a reason for hiding this comment

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

And unit for the argument

Copy link
Copy Markdown
Collaborator Author

@PabloCorcho PabloCorcho Oct 10, 2025

Choose a reason for hiding this comment

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

This is a private method. I don't know to what extent we need to be that careful given the docstrings... I don't understand the comment about the unit argument.

if spectra is not None and not isinstance(spectra, u.Quantity):
return spectra * default_unit
else:
"""
Ensure that ``spectra`` is an astropy Quantity and return it as flux density
per wavelength (F_lambda), using `self.wavelength` for conversions if
necessary.

Parameters
----------
spectra : array-like or ~astropy.units.Quantity
Spectrum values. If array-like (dimensionless), it will be multiplied by
`default_unit`. If a Quantity, it may be either F_nu (per frequency) or
F_lambda (per wavelength).
default_unit : ~astropy.units.Unit, optional
Target per-wavelength flux-density unit to return, e.g. Lsun / Å / cm²
(must be equivalent to power / area / length).

Returns
-------
astropy.units.Quantity
Spectrum expressed in `default_unit` (per-wavelength flux density).

Notes
-----
- Requires `self.wavelength` to be defined as an astropy Quantity with
Copy link
Copy Markdown
Owner

@paranoya paranoya Oct 10, 2025

Choose a reason for hiding this comment

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

Looks like we should implement _ensure_wavelength_units

Copy link
Copy Markdown
Collaborator Author

@PabloCorcho PabloCorcho Oct 10, 2025

Choose a reason for hiding this comment

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

I am not sure, self.wavelength cannot be anything else but a quantity with length dimensions or None. I think what I am missing here is checking whether the attribute is None.

length units.
- If `spectra` is already per wavelength, only a unit conversion is applied.
- If `spectra` is per frequency, the conversion uses
`equivalencies=u.spectral_density(self.wavelength)`.
"""
if spectra is None:
return None

if not isinstance(spectra, u.Quantity):
spectra = u.Quantity(spectra, default_unit)
# Check flux density units
if spectra.unit.is_equivalent(u.W / (u.m**2 * u.m)):
return spectra
Copy link
Copy Markdown
Owner

@paranoya paranoya Oct 10, 2025

Choose a reason for hiding this comment

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

Actually, what about the following helper function?

def ensure_unit(x, unit, just_equivalent=True, equivalence=None):
    if not isinstance(x, u.Quantity):
        return x << unit
    if x.unit.is_equivalent(unit):
        if just_equivalent:
            return x
        else:
            return x.to(unit)
    else:
            return x.to(unit, equivalencies=equivalence)

and then

ensure_unit(wavelength, u.AA)
ensure_unit(spectra, u.Lsun/u.angstrom, equivalence=u.spectral_density(wavelength))

Copy link
Copy Markdown
Collaborator Author

@PabloCorcho PabloCorcho Oct 10, 2025

Choose a reason for hiding this comment

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

I modified the check_unit method to also accept equivalencies, either taking arguments (spectral_density) or not ( spectral). The main difference is that the equivalency is used only when the input quantity is not equivalent to the target unit. This reduces the number of calls when using SSP data, which by default is in Flam.

elif spectra.unit.is_equivalent(u.W / (u.m**2 * u.Hz)):
return spectra.to(default_unit, equivalencies=u.spectral_density(self.wavelength))
else:
raise ValueError(
f"Cannot interpret spectral flux units {spectra.unit}. "
"Expected per-frequency or per-wavelength flux density."
)

def get_photons(self, spectra, spectra_err=None, mask_nan=True):
r"""Compute the photon flux from an input spectra.
Expand Down
20 changes: 16 additions & 4 deletions tests/test_observables.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,8 +12,10 @@ def setUpClass(self):
print("Setting SSP model")
self.dummy_wavelength = np.geomspace(100, 1e5, 3000) * u.angstrom
# Monocromatic SED
self.dummy_spectra = np.ones(self.dummy_wavelength.size
self.dummy_flam = np.ones(self.dummy_wavelength.size
) * constants.c / self.dummy_wavelength**2 * 3631 * u.Jy
self.dummy_fnu = np.ones(self.dummy_wavelength.size
) * 3631 * u.Jy

def test_default_dir(self):
self.assertTrue(
Expand Down Expand Up @@ -52,20 +54,30 @@ def test_filter(self):
0.32484839450189695),
"Unexpected effective transmission value")

# Interpolate filter to input wavelength array
filter.interpolate(self.dummy_wavelength)
flux, _ = filter.get_fnu(self.dummy_spectra)
# Use flam
flux, _ = filter.get_fnu(self.dummy_flam)
self.assertTrue(np.isclose(flux, 3631.0 * u.Jy),
f"Unexpected integrated flux value: {flux}")

mag, _ = filter.get_ab(self.dummy_spectra)
mag, _ = filter.get_ab(self.dummy_flam)
self.assertTrue(np.isclose(mag, 0.0, atol=1e-4),
f"Unexpected magnitude value: {mag}")
# Use fnu
flux, _ = filter.get_fnu(self.dummy_fnu)
self.assertTrue(np.isclose(flux, 3631.0 * u.Jy),
f"Unexpected integrated flux value: {flux}")

mag, _ = filter.get_ab(self.dummy_fnu)
self.assertTrue(np.isclose(mag, 0.0, atol=1e-4),
f"Unexpected magnitude value: {mag}")

fig = filter.plot(show=False)

def test_equivalent_width(self):
eqwidth = observables.EquivalentWidth.from_name("lick_ha")
ew, ew_err = eqwidth.compute_ew(self.dummy_wavelength, self.dummy_spectra)
ew, ew_err = eqwidth.compute_ew(self.dummy_wavelength, self.dummy_flam)
self.assertTrue(np.isfinite(ew), "Unexpected EW value")


Expand Down