diff --git a/skycatalogs/objects/base_object.py b/skycatalogs/objects/base_object.py index c319317d..deb3eb12 100644 --- a/skycatalogs/objects/base_object.py +++ b/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,7 +195,7 @@ def _get_sed(self, component=None, resolution=None, mjd=None): Returns ------- - A pair (sed, mag_norm) + galsim.SED object Must be implemented by subclass ''' @@ -206,7 +203,7 @@ def _get_sed(self, component=None, resolution=None, mjd=None): 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] diff --git a/skycatalogs/objects/diffsky_object.py b/skycatalogs/objects/diffsky_object.py index ae0af9a0..992d0d99 100644 --- a/skycatalogs/objects/diffsky_object.py +++ b/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 diff --git a/skycatalogs/objects/galaxy_object.py b/skycatalogs/objects/galaxy_object.py index c2d55a61..d8aa035f 100644 --- a/skycatalogs/objects/galaxy_object.py +++ b/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,7 +32,7 @@ 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') @@ -41,9 +40,7 @@ def _get_sed(self, component=None, resolution=None): 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) diff --git a/skycatalogs/objects/snana_object.py b/skycatalogs/objects/snana_object.py index 1648b3ed..d5ad78ab 100644 --- a/skycatalogs/objects/snana_object.py +++ b/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 diff --git a/skycatalogs/objects/sncosmo_object.py b/skycatalogs/objects/sncosmo_object.py index bbea5173..66c7ac94 100644 --- a/skycatalogs/objects/sncosmo_object.py +++ b/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 diff --git a/skycatalogs/objects/sso_object.py b/skycatalogs/objects/sso_object.py index df706a8a..041e9d46 100644 --- a/skycatalogs/objects/sso_object.py +++ b/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 diff --git a/skycatalogs/objects/star_object.py b/skycatalogs/objects/star_object.py index 78796a3f..5182b1a8 100644 --- a/skycatalogs/objects/star_object.py +++ b/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,7 +20,7 @@ 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: @@ -27,9 +28,7 @@ def get_gsobject_components(self, gsparams=None, rng=None): 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) diff --git a/skycatalogs/utils/sed_tools.py b/skycatalogs/utils/sed_tools.py index eb43a732..de453423 100644 --- a/skycatalogs/utils/sed_tools.py +++ b/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) diff --git a/tests/test_sed_tools.py b/tests/test_sed_tools.py new file mode 100644 index 00000000..0441d6f0 --- /dev/null +++ b/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()