Add galsim.EmissionLine

CHANGELOG.rst

@@ -37,6 +37,7 @@ New Features
   people might want to use to the installed galsim.utilities. (#1240)
 - Added `galsim.utilities.merge_sorted` which merges two or more sorted numpy arrays faster than
   the available numpy options. (#1243)
+- Added `galsim.EmissionLine` class to represent emission line SEDs. (#1249)
 
 
 Performance Improvements

galsim/sed.py

@@ -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
@@ -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
@@ -1095,5 +1097,133 @@ 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.  [default: 1.0]
+        flux:          The integrated flux of the line.  [default: 1.0]
+        wave_type:     The units of ``wavelength`` and ``fwhm`` above.  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=1.0,
+        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

tests/test_chromatic.py

@@ -530,14 +530,11 @@ def test_monochromatic_sed():
         obj_achrom.drawImage(image=im1)
         print('im1.max,sum = ', im1.array.max(), im1.array.sum())
 
-        # Next do this chromatically using an SED that is basically a single emission line
-        # at the given wavelength.
+        # Next do this chromatically using an emission line SED.
         psf_chrom = galsim.ChromaticOpticalPSF(lam=wave, diam=diam,
                                                aberrations=aberrations, obscuration=obscuration,
                                                nstruts=nstruts)
-        sed = galsim.SED(galsim.LookupTable([500, wave-1, wave, wave+1, 1000],
-                                            [0, 0, 1, 0, 0], 'linear'),
-                         wave_type='nm', flux_type='fphotons')
+        sed = galsim.EmissionLine(wave)
         gal_chrom = (gal_achrom * sed).withFlux(flux, bandpass=bandpass)
         obj_chrom = galsim.Convolve(gal_chrom, psf_chrom)
         im2 = obj_chrom.drawImage(bandpass, image=im1.copy())
@@ -2148,12 +2145,10 @@ def test_phot():
     """
     import time
 
-    # Use a fairly extreme SED to exaggerate the effect of the wavelength dependence of the
+    # Use an emission line SED to exaggerate the effect of the wavelength dependence of the
     # Airy profile.  Otherwise we need a lot of photons to notice the effect.
     bandpass = galsim.Bandpass(galsim.LookupTable([500,1000], [1,1], 'linear'), wave_type='nm')
-    sed = galsim.SED(galsim.LookupTable([500, 510,515,520, 970,975,980, 1000],
-                                        [0, 0,1,0, 0,1,0, 0], 'linear'),
-                     wave_type='nm', flux_type='fphotons')
+    sed = galsim.EmissionLine(510, fwhm=5) + galsim.EmissionLine(975, fwhm=5)
     flux = 1.e6
     gal_achrom = galsim.Sersic(n=2.8, half_light_radius=0.03, flux=flux)
     gal = (gal_achrom * sed).withFlux(flux, bandpass=bandpass)
@@ -2307,9 +2302,7 @@ def test_phot():
     #       one in InterpolatedChromaticObject._shoot.
     rng = galsim.BaseDeviate(1234)
     bp2 = galsim.Bandpass(galsim.LookupTable([400,601], [1,1], 'linear'), wave_type='nm')
-    sed2 = galsim.SED(galsim.LookupTable([400, 410,415,420, 510, 515, 520, 1000],
-                                        [0, 0,1,0, 0,2,0, 0], 'linear'),
-                      wave_type='nm', flux_type='fphotons')
+    sed2 = galsim.EmissionLine(415, fwhm=5) + 2*galsim.EmissionLine(515, fwhm=5)
     gal2 = (gal_achrom * sed2).withFlux(flux, bandpass=bp2)
     obj2 = galsim.Convolve(gal2, psf7)
     with assert_raises(galsim.GalSimRangeError):

tests/test_sed.py

@@ -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')
@@ -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)
 
 
@@ -1125,12 +1132,59 @@ def test_SED_calculateFlux_inf():
         wave_type='nm',
         flux_type='fphotons'
     )
+    sed3 = galsim.EmissionLine(622)
 
     bp = galsim.Bandpass("LSST_r.dat", 'nm')
     flux1 = sed1.calculateFlux(bp)
     flux2 = sed2.calculateFlux(bp)
-    print('flux = ', flux1, flux2)
+    flux3 = sed3.calculateFlux(bp)
+    print('flux = ', flux1, flux2, flux3)
     assert flux1 == flux2
+    assert flux1 == flux3
+
+
+@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)]