GalSim-developers / GalSim

The modular galaxy image simulation toolkit. Documentation:
http://galsim-developers.github.io/GalSim/
Other
224 stars 105 forks source link

SED flux type units lost in translation #1254

Closed sidneymau closed 11 months ago

sidneymau commented 11 months ago

When constructing an SED with a custom flux type via astropy units, the flux type appears to get improperly converted to fnu such that multiplying this SED mangles the units (e.g., when calling any SED.__mul__ method).

>>> disk_sed = galsim.SED(disk_sed_table, u.Angstrom, u.Lsun / u.Hz / u.Mpc**2, redshift=mock["redshift"][igal]))
>>> disk_sed.flux_type
'fnu'

This results in wrong fluxes

lsst_g = galsim.Bandpass("LSST_g.dat", wave_type="nm")
>>> (disk_sed).calculateFlux(lsst_g)
0.04530669166489878
>>> (1 * disk_sed).calculateFlux(lsst_g)
112691591447658.53

We can see that the flux type of the remade SED (even without applying any multiplicative factors) is the cause of this

>>> galsim.SED(galsim._LookupTable(disk_sed._spec.getArgs(), np.array(disk_sed._spec.getVals())), disk_sed.wave_type, disk_sed.flux_type, disk_sed.redshift).calculateFlux(lsst_g)
112691591447658.53

whereas forcing the original flux type produces the correct value

>>> galsim.SED(galsim._LookupTable(disk_sed._spec.getArgs(), np.array(disk_sed._spec.getVals())), disk_sed.wave_type, u.Lsun / u.Hz / u.Mpc**2, disk_sed.redshift).calculateFlux(lsst_g)
0.04530669166489878

Tested in GalSim Version 2.5.0

sidneymau commented 11 months ago

I noticed that the original flux type is preserved as SED._flux_type. Using this to reconstruct the SED in the above methods works:

>>> galsim.SED(galsim._LookupTable(disk_sed._spec.getArgs(), np.array(disk_sed._spec.getVals())), disk_sed.wave_type, disk_sed._flux_type, disk_sed.redshift).calculateFlux(lsst_g)
0.04530669166489878
sidneymau commented 11 months ago

Per the above, changing https://github.com/GalSim-developers/GalSim/blob/releases/2.5/galsim/sed.py#L530 from

flux_type = self.flux_type

to

flux_type = self._flux_type

fixes the particular problem I'm seeing. With the change,

>>> (disk_sed).calculateFlux(lsst_g)
0.04530669166489878
>>> (disk_sed * 1).calculateFlux(lsst_g)
0.04530669166489878
>>> (disk_sed * disk).calculateFlux(lsst_g)
0.04530669166489878