Merge pull request #98 from LSSTDESC/u/jchiang/magnorm_fix

magnorm fix

skycatalogs/objects/base_object.py

@@ -89,9 +89,6 @@ class BaseObject(object):
     Likely need a variant for SSO.
     '''
 
-    _bp500 = galsim.Bandpass(galsim.LookupTable([499, 500, 501], [0, 1, 0]),
-                             wave_type='nm').withZeropoint('AB')
-
     def __init__(self, ra, dec, id, object_type, belongs_to, belongs_index):
         '''
         Save at least minimum info needed a fixed (not SSO) object to
@@ -198,15 +195,15 @@ def _get_sed(self, component=None, resolution=None, mjd=None):
 
         Returns
         -------
-        A pair (sed, mag_norm)
+        galsim.SED object
 
         Must be implemented by subclass
         '''
         raise NotImplementedError('Must be implemented by BaseObject subclass if needed')
 
     def write_sed(self, sed_file_path, component=None, resolution=None,
                   mjd=None):
-        sed, _ = self._get_sed(component=component, resolution=None, mjd=None)
+        sed = self._get_sed(component=component, resolution=None, mjd=None)
 
         wl = sed.wave_list
         flambda = [sed(w) for w in wl]

skycatalogs/objects/diffsky_object.py

@@ -25,7 +25,7 @@ def _get_sed(self, component=None, resolution=None):
 
         Returns
         -------
-        A pair (sed, mag_norm)
+        galsim.SED object
         '''
 
         if component not in ['disk', 'bulge', 'knots']:
@@ -39,8 +39,7 @@ def _get_sed(self, component=None, resolution=None):
             sky_cat = self._belongs_to._sky_catalog
             self._seds = sky_cat.observed_sed_factory.create(pixel, self.id, z_h, z)
 
-        magnorm = None
-        return self._seds[component], magnorm
+        return self._seds[component]
 
     def get_knot_size(self,z):
         """
@@ -147,7 +146,7 @@ def get_gsobject_components(self, gsparams=None, rng=None):
         return obj_dict
 
     def get_observer_sed_component(self, component, mjd=None, resolution=None):
-        sed, _ = self._get_sed(component)
+        sed = self._get_sed(component)
         if sed is not None:
             sed = self._apply_component_extinction(sed)
         return sed

skycatalogs/objects/galaxy_object.py

@@ -12,17 +12,16 @@ class GalaxyObject(BaseObject):
 
     def _get_sed(self, component=None, resolution=None):
         '''
-        Return sed and mag_norm for a galaxy component or for a star
+        Return sed for a galaxy component
         Parameters
         ----------
         component    one of 'bulge', 'disk', 'knots' for now. Other components
-                     may be supported.  Ignored for stars
-        resolution   desired resolution of lambda in nanometers. Ignored
-                     for stars.
+                     may be supported.
+        resolution   desired resolution of lambda in nanometers.
 
         Returns
         -------
-        A pair (sed, mag_norm)
+        galsim.SED object
         '''
         if component not in ['disk', 'bulge', 'knots']:
             raise ValueError(f'Cannot fetch SED for component type {component}')
@@ -33,17 +32,15 @@ def _get_sed(self, component=None, resolution=None):
 
         # if values are all zeros or nearly no point in trying to convert
         if max(th_val) < np.finfo('float').resolution:
-            return None, 0.0
+            return None
 
         z_h = self.get_native_attribute('redshift_hubble')
         z = self.get_native_attribute('redshift')
 
         sky_cat = self._belongs_to._sky_catalog
         sed = sky_cat.observed_sed_factory.create(th_val, z_h, z,
                                                   resolution=resolution)
-        magnorm = sky_cat.observed_sed_factory.magnorm(th_val, z_h)
-
-        return sed, magnorm
+        return sed
 
     def get_knot_size(self,z):
         """
@@ -146,7 +143,7 @@ def get_gsobject_components(self, gsparams=None, rng=None):
 
 
     def get_observer_sed_component(self, component, mjd=None, resolution=None):
-        sed, _ = self._get_sed(component=component, resolution=resolution)
+        sed = self._get_sed(component=component, resolution=resolution)
         if sed is not None:
             sed = self._apply_component_extinction(sed)
 

skycatalogs/objects/snana_object.py

@@ -32,9 +32,9 @@ def _get_sed(self, mjd=None):
         mjd_start = self.get_native_attribute('start_mjd')
         mjd_end = self.get_native_attribute('end_mjd')
         if mjd < mjd_start or mjd > mjd_end:
-            return None, 0.0
+            return None
 
-        return self._linear_interp_SED(mjd), 0.0
+        return self._linear_interp_SED(mjd)
 
     def _apply_flux_correction(self, flux, snana_band, mjd):
         def _flux_ratio(mag):
@@ -95,7 +95,7 @@ def get_observer_sed_component(self, component, mjd=None):
             txt = 'SnananObject.get_observer_sed_component: no mjd specified for this call\n'
             txt += 'nor when generating object list'
             raise ValueError(txt)
-        sed, _ = self._get_sed(mjd=mjd)
+        sed = self._get_sed(mjd=mjd)
         if sed is not None:
             sed = self._apply_component_extinction(sed)
         return sed

skycatalogs/objects/sncosmo_object.py

@@ -14,9 +14,8 @@ def _get_sed(self, mjd=None):
         sn = SncosmoModel(params=params)
 
         if mjd < sn.mintime() or mjd > sn.maxtime():
-            return None, 0.0
-        # Already normalized so magnorm is zero
-        return sn.get_sed(mjd), 0.0
+            return None
+        return sn.get_sed(mjd)
 
     def get_gsobject_components(self, gsparams=None, rng=None):
         if gsparams is not None:
@@ -30,7 +29,7 @@ def get_observer_sed_component(self, component, mjd=None):
             txt = 'SncosmoObject._get_sed: no mjd specified for this call\n'
             txt += 'nor when generating object list'
             raise ValueError(txt)
-        sed, _ = self._get_sed(mjd=mjd)
+        sed = self._get_sed(mjd=mjd)
         if sed is not None:
             sed = self._apply_component_extinction(sed)
         return sed

skycatalogs/objects/sso_object.py

@@ -4,6 +4,7 @@
 import galsim
 from .base_object import BaseObject, ObjectCollection
 from lsst.sphgeom import UnitVector3d, LonLat
+from ..utils import normalize_sed
 
 EXPOSURE_DEFAULT = 30.0         # seconds
 __all__ = ['SsoObject', 'SsoCollection', 'EXPOSURE_DEFAULT']
@@ -27,7 +28,7 @@ def mjd(self):
 
     def _get_sed(self, mjd=None):
         '''
-        returns a SED and magnorm
+        returns the SED
         mjd is required
         '''
         if SsoObject._solar_sed is None:
@@ -37,7 +38,8 @@ def _get_sed(self, mjd=None):
             # directly or are there other effects to be applied?
             # Have to find it by looking for entry for this id, this mjd
             # Do we look for specific entry or do we allow interpolation?
-        return SsoObject._solar_sed, self.get_native_attribute('trailed_source_mag')
+        magnorm = self.get_native_attribute('trailed_source_mag')
+        return normalize_sed(SsoObject._solar_sed, magnorm)
 
     def get_gsobject_components(self, gsparams=None, rng=None,
                                 streak=True):
@@ -91,8 +93,7 @@ def get_observer_sed_component(self, component, mjd=None):
             txt += 'nor when generating object list'
             raise ValueError(txt)
 
-        sed, magnorm = self._get_sed(mjd=mjd)
-        sed = sed.withMagnitude(magnorm, self._bp500)
+        sed = self._get_sed(mjd=mjd)
 
         # no extinction
         return sed

skycatalogs/objects/star_object.py

@@ -2,6 +2,7 @@
 import numpy as np
 import galsim
 from .base_object import BaseObject
+from ..utils import normalize_sed
 
 __all__ = ['StarObject']
 
@@ -19,17 +20,15 @@ def _get_sed(self, mjd=None, redshift=0):
         fpath = os.path.join(os.getenv('SIMS_SED_LIBRARY_DIR'),
                              self.get_native_attribute('sed_filepath'))
 
-        return self._get_sed_from_file(fpath), mag_norm
+        return normalize_sed(self._get_sed_from_file(fpath), mag_norm)
 
     def get_gsobject_components(self, gsparams=None, rng=None):
         if gsparams is not None:
             gsparams = galsim.GSParams(**gsparams)
         return {'this_object': galsim.DeltaFunction(gsparams=gsparams)}
 
     def get_observer_sed_component(self, component, mjd=None):
-        sed, magnorm = self._get_sed(mjd=mjd)
-
-        sed = sed.withMagnitude(magnorm, self._bp500)
+        sed = self._get_sed(mjd=mjd)
 
         if sed is not None:
             sed = self._apply_component_extinction(sed)

skycatalogs/utils/sed_tools.py

@@ -12,7 +12,8 @@
 import galsim
 
 __all__ = ['TophatSedFactory', 'DiffskySedFactory', 'SsoSedFactory',
-           'MilkyWayExtinction', 'get_star_sed_path', 'generate_sed_path']
+           'MilkyWayExtinction', 'get_star_sed_path', 'generate_sed_path',
+           'normalize_sed']
 
 
 class TophatSedFactory:
@@ -331,3 +332,35 @@ def generate_sed_path(ids, subdir, cmp):
     '''
     r = [f'{subdir}/{cmp}_{id}.txt' for id in ids]
     return r
+
+
+def normalize_sed(sed, magnorm, wl=500*u.nm):
+    """
+    Set the normalization of a GalSim SED object given a monochromatic
+    magnitude at a reference wavelength.
+
+    Parameters
+    ----------
+    sed : galsim.SED
+        The GalSim SED object.
+    magnorm : float
+        The monochromatic magnitude at the reference wavelength.
+    wl : astropy.units.nm
+        The reference wavelength.
+
+    Returns
+    -------
+    galsim.SED : The renormalized SED object.
+    """
+    # Compute the flux density from magnorm in units of erg/cm^2/s/nm.
+    fnu = (magnorm * u.ABmag).to_value(u.erg/u.s/u.cm**2/u.Hz)
+    flambda = fnu * (astropy.constants.c/wl**2).to_value(u.Hz/u.nm)
+
+    # GalSim expects the flux density in units of photons/cm^2/s/nm,
+    # so divide flambda by the photon energy at the reference
+    # wavelength.
+    hnu = (astropy.constants.h * astropy.constants.c / wl).to_value(u.erg)
+
+    flux_density = flambda/hnu
+
+    return sed.withFluxDensity(flux_density, wl)

tests/test_sed_tools.py

@@ -0,0 +1,44 @@
+import unittest
+import astropy.units as u
+import astropy.constants
+import numpy as np
+import galsim
+from skycatalogs.utils import normalize_sed
+
+
+class NormalizeSedTestCase(unittest.TestCase):
+    """TestCase class for normalizing SEDs with magnorm."""
+
+    def setUp(self):
+        """No set up needed."""
+        pass
+
+    def tearDown(self):
+        """Nothing to tear down."""
+        pass
+
+    def test_normalize_sed(self):
+        """Test the normalize_sed function"""
+        sed = galsim.SED(lambda wl: 1, "nm", "flambda", fast=False)
+        wl = 400 * u.nm
+
+        # Check that the ratio of differently normalized SEDs have the
+        # expected magnitude difference at the reference wavelength:
+        magnorm0 = 20.0
+        magnorm1 = 25.0
+        sed0 = normalize_sed(sed, magnorm0, wl=wl)
+        sed1 = normalize_sed(sed, magnorm1, wl=wl)
+
+        self.assertAlmostEqual(sed0(wl) / sed1(wl), 10 ** ((magnorm1 - magnorm0) / 2.5))
+
+        # Check the implied magnitude of the renormalized SED at the
+        # reference wavelength against the magnorm value.
+        hnu = (astropy.constants.h * astropy.constants.c / wl).to_value(u.erg)
+        fnu = sed0(wl) * hnu / (astropy.constants.c / wl**2).to_value(u.Hz / u.nm)
+        mag = -2.5 * np.log10(fnu) - 48.60
+
+        self.assertAlmostEqual(mag, magnorm0)
+
+
+if __name__ == "__main__":
+    unittest.main()