#1365 Fix an error in SED * Bandpass when not using nm for the wavelengths

galsim/sed.py

@@ -502,26 +502,24 @@ def _mul_sed(self, other):
 
         wave_list, blue_limit, red_limit = utilities.combine_wave_list(self, other)
         if fast:
+            zfactor1 = (1.+redshift) / (1.+self.redshift)
+            zfactor2 = (1.+redshift) / (1.+other.redshift)
             if (isinstance(self._fast_spec, LookupTable)
                     and not self._fast_spec.x_log
                     and not self._fast_spec.f_log):
-                x = wave_list / (1.0 + self.redshift)
+                x = wave_list / (1.0 + redshift)
                 # Add in 500 uniformly spaced values to help improve accuracy.
                 x = utilities.merge_sorted([x, np.linspace(x[0], x[-1], 500)])
-                zfactor2 = (1.+redshift) / (1.+other.redshift)
-                f = self._fast_spec(x) * other._fast_spec(x*zfactor2)
+                f = self._fast_spec(x*zfactor1) * other._fast_spec(x*zfactor2)
                 spec = _LookupTable(x, f, self._fast_spec.interpolant)
             elif (isinstance(other._fast_spec, LookupTable)
                     and not other._fast_spec.x_log
                     and not other._fast_spec.f_log):
-                x = wave_list / (1.0 + other.redshift)
+                x = wave_list / (1.0 + redshift)
                 x = utilities.merge_sorted([x, np.linspace(x[0], x[-1], 500)])
-                zfactor1 = (1.+redshift) / (1.+other.redshift)
-                f = self._fast_spec(x*zfactor1) * other._fast_spec(x)
+                f = self._fast_spec(x*zfactor1) * other._fast_spec(x*zfactor2)
                 spec = _LookupTable(x, f, other._fast_spec.interpolant)
             else:
-                zfactor1 = (1.+redshift) / (1.+self.redshift)
-                zfactor2 = (1.+redshift) / (1.+other.redshift)
                 spec = lambda w: self._fast_spec(w * zfactor1) * other._fast_spec(w * zfactor2)
         else:
             spec = lambda w: self(w * (1.+redshift)) * other(w * (1.+redshift))
@@ -533,7 +531,7 @@ def _mul_sed(self, other):
     def _mul_bandpass(self, other):
         """Equivalent to self * other when other is a Bandpass"""
         wave_list, blue_limit, red_limit = utilities.combine_wave_list(self, other)
-        zfactor = (1.0+self.redshift) / other.wave_factor
+        zfactor = (1.0+self.redshift) * other.wave_factor
         if self.fast:
             if (isinstance(self._fast_spec, LookupTable)
                     and not self._fast_spec.x_log
@@ -643,20 +641,24 @@ def __div__(self, other):
         # Enable division by scalars or dimensionless callables (including dimensionless SEDs.)
         if isinstance(other, SED) and other.spectral:
             raise GalSimSEDError("Cannot divide by spectral SED.", other)
+        wave_type = 'nm'
+        flux_type = 'fphotons' if self.spectral else '1'
         if hasattr(other, '__call__'):
             spec = lambda w: self(w * (1.0 + self.redshift)) / other(w * (1.0 + self.redshift))
         elif isinstance(self._spec, LookupTable):
             # If other is not a function, then there is no loss of accuracy by applying the
             # factor directly to the LookupTable, if that's what we are using.
             # Make sure to keep the same properties about the table, flux_type, wave_type.
+            wave_type = self.wave_type
+            flux_type = self.flux_type
             x = self._spec.getArgs()
             f = [ val / other for val in self._spec.getVals() ]
             spec = _LookupTable(x, f, x_log=self._spec.x_log, f_log=self._spec.f_log,
                                 interpolant=self._spec.interpolant)
         else:
             spec = lambda w: self(w * (1.0 + self.redshift)) / other
 
-        return SED(spec, flux_type=self.flux_type, wave_type=self.wave_type,
+        return SED(spec, flux_type=flux_type, wave_type=wave_type,
                    redshift=self.redshift, fast=self.fast,
                    _wave_list=self.wave_list,
                    _blue_limit=self.blue_limit, _red_limit=self.red_limit)

tests/test_sed.py

@@ -320,6 +320,24 @@ def test_SED_mul():
     np.testing.assert_almost_equal((sed2*sed4)(10), 8.0, 10, "Found wrong value in SED.__mul__")
     np.testing.assert_almost_equal((sed3*sed4)(30), 12.0, 10, "Found wrong value in SED.__mul__")
 
+    # Mix callable and tabulated SEDs, with the spectral SED carrying the redshift.
+    # Both multiplication orders must evaluate each operand at the observed wavelength.
+    wave = np.array([1., 2., 3., 4., 5.])
+    x = np.array([2.5, 3.5])
+    for wave_type, factor in [('nm', 1.), ('Angstrom', 10.), (units.m, 1.e-9)]:
+        specs = [lambda w: w / factor,
+                 galsim.LookupTable(wave * factor, wave, 'linear')]
+        for fast in [True, False]:
+            for z in [0, 0.2, 0.4]:
+                for spectral_spec in specs:
+                    a = galsim.SED(spectral_spec, wave_type, 'fphotons', fast=fast).atRedshift(z)
+                    for dimensionless_spec in specs:
+                        b = galsim.SED(dimensionless_spec, wave_type, '1', fast=fast)
+                        expected = a(x) * b(x)
+                        # A tabulated product approximates the quadratic between sample points.
+                        np.testing.assert_allclose((a*b)(x), expected, rtol=1.e-5)
+                        np.testing.assert_allclose((b*a)(x), expected, rtol=1.e-5)
+
 @timer
 def test_SED_div():
     """Check that SEDs divide like I think they should...
@@ -356,6 +374,24 @@ def test_SED_div():
             np.testing.assert_almost_equal(d(x), a(x)/4.2/2/e(x), 10,
                                            err_msg="Found wrong value in SED.__div__")
 
+    # Division must preserve the evaluated photon density regardless of the input units.
+    wave = np.array([1., 2., 3., 4., 5.])
+    for wave_type, factor in [('nm', 1.), ('Angstrom', 10.), (units.m, 1.e-9)]:
+        specs = [lambda w: w / factor,
+                 galsim.LookupTable(wave * factor, wave, 'linear')]
+        for flux_type in ['fphotons', 'flambda', 'fnu', '1']:
+            for fast in [True, False]:
+                for z in [0, 0.2, 0.4]:
+                    # Use table nodes to avoid interpolation error in the quotient.
+                    x = np.array([2., 3.]) * (1.+z)
+                    for spec in specs:
+                        a = galsim.SED(spec, wave_type, flux_type, fast=fast).atRedshift(z)
+                        np.testing.assert_allclose((a/4.2)(x), a(x)/4.2, rtol=1.e-12)
+                        b = lambda w: w**2
+                        np.testing.assert_allclose((a/b)(x), a(x)/b(x), rtol=1.e-12)
+                        e = galsim.SED('wave', 'nm', '1', fast=fast)
+                        np.testing.assert_allclose((a/e)(x), a(x)/e(x), rtol=1.e-12)
+
     # Can't divide by spectral SED
     with assert_raises(galsim.GalSimSEDError):
         a0_lt / a0_fn
@@ -517,10 +553,24 @@ def test_SED_withFlux():
     """ Check that setting the flux works.
     """
     rband = galsim.Bandpass(os.path.join(bppath, 'LSST_r.dat'), 'nm')
+
+    # Make some equivalent bandpasses made in different ways.
+    # This one checks giving the file name using Path. Should be exactly equivalent.
     rband2 = galsim.Bandpass(Path(rband._orig_tp), 'nm')
     np.testing.assert_array_equal(rband.wave_list, rband2.wave_list)
     np.testing.assert_array_equal(rband(rband.wave_list), rband2(rband2.wave_list))
 
+    # These two are functionally the same, but uses Angstrom and meter for their wavelengths.
+    rband_A = galsim.Bandpass(
+        galsim.LookupTable(rband.wave_list * 10, rband(rband.wave_list), 'linear'), 'Angstrom')
+    rband_m = galsim.Bandpass(
+        galsim.LookupTable(rband.wave_list * 1.e-9, rband(rband.wave_list), 'linear'), units.m)
+
+    # This one has a finer binning for the wavelengths, but follows the same shape.
+    fine_wave = np.linspace(rband.blue_limit, rband.red_limit, 169101)
+    rband_fine = galsim.Bandpass(
+        galsim.LookupTable(fine_wave, rband(fine_wave), 'linear'), 'nm')
+
     for z in [0, 0.2, 0.4]:
         for fast in [True, False]:
             for sed in [
@@ -551,16 +601,15 @@ def test_SED_withFlux():
                     ab.calculateFlux(bolo_bp), 1.0, 3,
                     "Calculating SED flux from sed * bp failed."
                 )
+                for bp in [rband2, rband_A, rband_m]:
+                    filtered = sed * bp
+                    np.testing.assert_allclose(filtered([600., 650.]), ab([600., 650.]))
+                    np.testing.assert_allclose(filtered.calculateFlux(bolo_bp),
+                                               ab.calculateFlux(bolo_bp))
 
                 # 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."
@@ -570,7 +619,7 @@ def test_SED_withFlux():
                 ba = rband_fine * sed
                 np.testing.assert_array_almost_equal(
                     ba.calculateFlux(bolo_bp), 1.0, 5,
-                    "Calculating SED flux from sed * bp failed."
+                    "Calculating SED flux from bp * sed failed."
                 )
 
     # Invalid for dimensionless SED