Skip to content

Commit

Permalink
Add galsim.EmissionLine
Browse files Browse the repository at this point in the history
  • Loading branch information
@jmeyers314
jmeyers314 committed Oct 3, 2023
1 parent 363120b commit 33b454d
Show file tree
Hide file tree
Showing 2 changed files with 223 additions and 44 deletions.
130 changes: 129 additions & 1 deletion galsim/sed.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@
# and/or other materials provided with the distribution.
#

__all__ = [ 'SED' ]
__all__ = [ 'SED', 'EmissionLine' ]

import numpy as np
from astropy import units
Expand Down Expand Up @@ -789,6 +789,8 @@ def calculateFlux(self, bandpass):
"""
if self.dimensionless:
raise GalSimSEDError("Cannot calculate flux of dimensionless SED.", self)
if bandpass is None: # compute bolometric flux
bandpass = Bandpass(lambda w: 1., 'nm', blue_limit=0., red_limit=1e100)
if len(bandpass.wave_list) > 0 or len(self.wave_list) > 0:
slop = 1e-6 # nm
if (self.blue_limit > bandpass.blue_limit + slop
Expand Down Expand Up @@ -1095,5 +1097,131 @@ def __setstate__(self, d):
self._initialize_spec()
self._setup_funcs()


class EmissionLine(SED):
"""Emission line SED.
Creates a simple triangle-shaped emission line with a given wavelength and
FWHM.
Parameters:
wavelength: The wavelength of the line.
fwhm: The full-width-half-max of the line.
flux: The integrated flux of the line. [default: 1.0]
wave_type: The units of wavelength used for this SED. See SED
constructor for options. [default: 'nm']
flux_type: The units of flux used for this SED. See SED constructor
for options. [default: 'fphotons']
**kwargs: Any additional keyword arguments to pass to the `SED`
constructor.
"""
def __init__(
self, wavelength, fwhm,
flux=1.0,
wave_type='nm',
flux_type='fphotons',
**kwargs
):
self.wavelength = wavelength
self.fwhm = fwhm
self.flux = flux
spec = LookupTable(
[0.0, wavelength-fwhm, wavelength, wavelength+fwhm, np.inf],
[0, 0, flux/fwhm, 0, 0],
interpolant='linear'
)
super().__init__(
spec,
wave_type=wave_type,
flux_type=flux_type,
**kwargs
)

def atRedshift(self, redshift):
"""Return a new `EmissionLine` with redshifted wavelengths.
Parameters:
redshift: The redshift for the returned `EmissionLine`
Returns:
the redshifted `EmissionLine`.
"""
if redshift == self.redshift:
return self
if redshift <= -1:
raise GalSimRangeError("Invalid redshift", redshift, -1.)
zfactor = (1.0 + redshift) / (1.0 + self.redshift)
return EmissionLine(
self.wavelength * zfactor,
self.fwhm * zfactor,
flux=self.flux,
wave_type=self.wave_type,
flux_type=self.flux_type,
redshift=redshift,
fast=self.fast
)

def __mul__(self, other):
if isinstance(other, (int, float)):
flux = self.flux * other
return EmissionLine(
self.wavelength,
self.fwhm,
flux=flux,
wave_type=self.wave_type,
flux_type=self.flux_type,
redshift=self.redshift,
fast=self.fast
)
else:
return super().__mul__(other)

def __div__(self, other):
if isinstance(other, (int, float)):
flux = self.flux / other
return EmissionLine(
self.wavelength,
self.fwhm,
flux=flux,
wave_type=self.wave_type,
flux_type=self.flux_type,
redshift=self.redshift,
fast=self.fast
)
else:
return super().__div__(other)

__truediv__ = __div__

def __eq__(self, other):
return (self is other or
(isinstance(other, EmissionLine) and
self.wavelength == other.wavelength and
self.fwhm == other.fwhm and
self.flux == other.flux and
self.wave_type == other.wave_type and
self.flux_type == other.flux_type and
self.redshift == other.redshift))

def __hash__(self):
return hash(("galsim.EmissionLine", self.wavelength, self.fwhm, self.flux))

def __repr__(self):
outstr = ('galsim.EmissionLine(%r, %r, flux=%r, wave_type=%r, flux_type=%r, redshift=%r,'
' fast=%r)')%(
self.wavelength, self.fwhm, self.flux, self.wave_type, self._flux_type,
self.redshift, self.fast)
return outstr

def __str__(self):
outstr = 'galsim.EmissionLine(wavelength=%s, fwhm=%s'%(self.wavelength, self.fwhm)
if self.flux != 1.0:
outstr += ', flux=%s'%self.flux
if self.redshift != 0.0:
outstr += ', redshift=%s'%self.redshift
outstr += ')'
return outstr


# Put this at the bottom to avoid circular import errors.
from .bandpass import Bandpass
137 changes: 94 additions & 43 deletions tests/test_sed.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -512,48 +512,53 @@ def test_SED_withFlux():
rband = galsim.Bandpass(os.path.join(bppath, 'LSST_r.dat'), 'nm')
for z in [0, 0.2, 0.4]:
for fast in [True, False]:
a = galsim.SED(os.path.join(sedpath, 'CWW_E_ext.sed'), wave_type='ang',
flux_type='flambda', fast=fast)
b = galsim.SED('wave', wave_type='nm', flux_type='fphotons')
if z != 0:
a = a.atRedshift(z)
b = b.atRedshift(z)
a = a.withFlux(1.0, rband)
b = b.withFlux(1.0, rband)
np.testing.assert_array_almost_equal(a.calculateFlux(rband), 1.0, 5,
"Setting SED flux failed.")
np.testing.assert_array_almost_equal(b.calculateFlux(rband), 1.0, 5,
"Setting SED flux failed.")

# Should be almost equivalent to multiplying an SED * Bandpass and computing the
# "bolometric" flux. The above is a bit more accurate, since it correctly does
# the integration of the product of two linear segments between each tabulated point.
ab = a * rband
bb = b * rband
bolo_bp = galsim.Bandpass('1', blue_limit=ab.blue_limit, red_limit=ab.red_limit,
wave_type='nm')
np.testing.assert_array_almost_equal(ab.calculateFlux(bolo_bp), 1.0, 3,
"Calculating SED flux from sed * bp failed.")
np.testing.assert_array_almost_equal(bb.calculateFlux(bolo_bp), 1.0, 3,
"Calculating SED flux from sed * bp failed.")

# If one or the other table has finer wavelength gridding, then the agreement
# will be better. Check with finer gridding for rband.
fine_wave = np.linspace(ab.blue_limit, ab.red_limit, 169101)
rband_fine = galsim.Bandpass(galsim.LookupTable(fine_wave, rband(fine_wave), 'linear'),
'nm')
ab = a * rband_fine
bb = b * rband_fine

np.testing.assert_array_almost_equal(ab.calculateFlux(bolo_bp), 1.0, 5,
"Calculating SED flux from sed * bp failed.")
np.testing.assert_array_almost_equal(bb.calculateFlux(bolo_bp), 1.0, 5,
"Calculating SED flux from sed * bp failed.")

# Multiplying in the other order also works.
ba = rband_fine * a
np.testing.assert_array_almost_equal(ba.calculateFlux(bolo_bp), 1.0, 5,
"Calculating SED flux from sed * bp failed.")
for sed in [
galsim.SED(os.path.join(sedpath, 'CWW_E_ext.sed'), wave_type='ang',
flux_type='flambda', fast=fast),
galsim.SED('wave', wave_type='nm', flux_type='fphotons'),
galsim.EmissionLine(620.0, 1.0) + 2*galsim.EmissionLine(450.0, 0.5)
]:
if z != 0:
sed = sed.atRedshift(z)
sed = sed.withFlux(1.0, rband)
np.testing.assert_array_almost_equal(
sed.calculateFlux(rband), 1.0, 5,
"Setting SED flux failed."
)

# Should be almost equivalent to multiplying an SED * Bandpass and computing the
# "bolometric" flux. The above is a bit more accurate, since it correctly does
# the integration of the product of two linear segments between each tabulated point.
ab = sed * rband
bolo_bp = galsim.Bandpass(
'1', blue_limit=ab.blue_limit, red_limit=ab.red_limit,
wave_type='nm'
)
np.testing.assert_array_almost_equal(
ab.calculateFlux(bolo_bp), 1.0, 3,
"Calculating SED flux from sed * bp failed."
)

# If one or the other table has finer wavelength gridding, then the agreement
# will be better. Check with finer gridding for rband.
fine_wave = np.linspace(ab.blue_limit, ab.red_limit, 169101)
rband_fine = galsim.Bandpass(
galsim.LookupTable(fine_wave, rband(fine_wave), 'linear'),
'nm'
)
ab = sed * rband_fine

np.testing.assert_array_almost_equal(
ab.calculateFlux(bolo_bp), 1.0, 5,
"Calculating SED flux from sed * bp failed."
)

# Multiplying in the other order also works.
ba = rband_fine * sed
np.testing.assert_array_almost_equal(
ba.calculateFlux(bolo_bp), 1.0, 5,
"Calculating SED flux from sed * bp failed."
)

# Invalid for dimensionless SED
flat = galsim.SED(2.0, 'nm', '1')
Expand Down Expand Up @@ -1003,7 +1008,9 @@ def test_ne():
galsim.SED(spec1, wave_type='nm', flux_type='fnu'),
galsim.SED(spec1, 'nm', 'flambda', redshift=1.0),
galsim.SED(spec2, 'nm', 'flambda'),
galsim.SED(spec3, 'nm', 'flambda')]
galsim.SED(spec3, 'nm', 'flambda'),
galsim.EmissionLine(500.0, 1.0),
]
check_all_diff(seds)


Expand Down Expand Up @@ -1132,6 +1139,50 @@ def test_SED_calculateFlux_inf():
print('flux = ', flux1, flux2)
assert flux1 == flux2


@timer
def test_emission_line():
spectral = galsim.SED("vega.txt", wave_type='nm', flux_type='flambda')
dimensionless = galsim.SED("1", wave_type='nm', flux_type='1')

for wavelength, fwhm in [
(500.0, 1.0),
(650.0, 0.3),
(700.0, 4.3)
]:
for sed in [
galsim.EmissionLine(wavelength, fwhm),
galsim.EmissionLine(wavelength*10, fwhm*10, wave_type='ang')
]:
np.testing.assert_allclose(sed.calculateFlux(None), 1.0)
np.testing.assert_allclose((sed*2).calculateFlux(None), 2.0)
np.testing.assert_allclose((3*sed).calculateFlux(None), 3.0)
np.testing.assert_allclose((sed/2).calculateFlux(None), 0.5)
np.testing.assert_allclose(sed(wavelength), 1./fwhm)
np.testing.assert_allclose(sed(wavelength+fwhm), 0.0, rtol=0, atol=1e-11)
np.testing.assert_allclose(sed(wavelength-fwhm), 0.0, rtol=0, atol=1e-11)
np.testing.assert_allclose(sed(wavelength+fwhm/2), 0.5/fwhm, rtol=0, atol=1e-11)
np.testing.assert_allclose(sed(wavelength-fwhm/2), 0.5/fwhm, rtol=0, atol=1e-11)
np.testing.assert_allclose((sed*dimensionless).calculateFlux(None), 1.0)

check_pickle(sed)
check_pickle(sed*2)
check_pickle(sed/3)
check_pickle(sed.atRedshift(0.3))

with np.testing.assert_raises(galsim.GalSimIncompatibleValuesError):
sed * spectral
with np.testing.assert_raises(galsim.GalSimSEDError):
sed / spectral
with np.testing.assert_raises(galsim.GalSimSEDError):
spectral / sed

sed = sed.atRedshift(1.1)
assert sed is sed.atRedshift(1.1)
with np.testing.assert_raises(galsim.GalSimRangeError):
sed.atRedshift(-2.0)


if __name__ == "__main__":
testfns = [v for k, v in vars().items() if k[:5] == 'test_' and callable(v)]
for testfn in testfns:
Expand Down

0 comments on commit 33b454d

Please sign in to comment.