From b7842173d0fd9fa036cbda9de0897dcaaa0a2ebb Mon Sep 17 00:00:00 2001 From: Mike Jarvis Date: Thu, 15 Oct 2026 16:01:53 -0400 Subject: [PATCH 1/3] Fix an error in _mul_bandpass when reference wavelength is not nm (#1365) --- galsim/sed.py | 2 +- tests/test_sed.py | 27 ++++++++++++++++++++------- 2 files changed, 21 insertions(+), 8 deletions(-) diff --git a/galsim/sed.py b/galsim/sed.py index 704957d446..2393f75686 100644 --- a/galsim/sed.py +++ b/galsim/sed.py @@ -533,7 +533,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 diff --git a/tests/test_sed.py b/tests/test_sed.py index 7430fb0403..d04ab10732 100644 --- a/tests/test_sed.py +++ b/tests/test_sed.py @@ -517,10 +517,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 +565,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 +583,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 From ce437ad40e9d73c143fec361e18f3b9ebcbecdac Mon Sep 17 00:00:00 2001 From: Mike Jarvis Date: Thu, 15 Oct 2026 20:35:30 -0400 Subject: [PATCH 2/3] Fix unit handling for redshifted SED * SED (#1365) --- galsim/sed.py | 14 ++++++-------- tests/test_sed.py | 18 ++++++++++++++++++ 2 files changed, 24 insertions(+), 8 deletions(-) diff --git a/galsim/sed.py b/galsim/sed.py index 2393f75686..a55dc3e656 100644 --- a/galsim/sed.py +++ b/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)) diff --git a/tests/test_sed.py b/tests/test_sed.py index d04ab10732..401013c187 100644 --- a/tests/test_sed.py +++ b/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... From cc6b9df23da14b4408c504f0f410df90ef312ba8 Mon Sep 17 00:00:00 2001 From: Mike Jarvis Date: Thu, 15 Oct 2026 20:36:20 -0400 Subject: [PATCH 3/3] Fix unit handling of SED / SED (#1365) --- galsim/sed.py | 6 +++++- tests/test_sed.py | 18 ++++++++++++++++++ 2 files changed, 23 insertions(+), 1 deletion(-) diff --git a/galsim/sed.py b/galsim/sed.py index a55dc3e656..5e16baa4e4 100644 --- a/galsim/sed.py +++ b/galsim/sed.py @@ -641,12 +641,16 @@ 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, @@ -654,7 +658,7 @@ def __div__(self, other): 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) diff --git a/tests/test_sed.py b/tests/test_sed.py index 401013c187..d4fc99d801 100644 --- a/tests/test_sed.py +++ b/tests/test_sed.py @@ -374,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