From 851e979ad253ca577e0fdd83fc4830cd3b68497e Mon Sep 17 00:00:00 2001 From: Josh Meyers Date: Mon, 10 Jul 2023 14:56:49 -0700 Subject: [PATCH 1/5] Add seeing_exp param to AtmosphericScreen --- galsim/phase_psf.py | 19 ++++++++++++------- galsim/phase_screens.py | 42 +++++++++++++++++++++++++++++------------ 2 files changed, 42 insertions(+), 19 deletions(-) diff --git a/galsim/phase_psf.py b/galsim/phase_psf.py index be16b845eea..76c6c8e88ba 100644 --- a/galsim/phase_psf.py +++ b/galsim/phase_psf.py @@ -914,7 +914,7 @@ def wavefront(self, u, v, t, theta=(0.0*radians, 0.0*radians)): else: return self._layers[0].wavefront(u, v, t, theta) - def wavefront_gradient(self, u, v, t, theta=(0.0*radians, 0.0*radians)): + def wavefront_gradient(self, u, v, t, w=None, theta=(0.0*radians, 0.0*radians)): """ Compute cumulative wavefront gradient due to all phase screens in `PhaseScreenList`. Parameters: @@ -927,6 +927,7 @@ def wavefront_gradient(self, u, v, t, theta=(0.0*radians, 0.0*radians)): will be used for all u, v. If scalar, then the size will be broadcast up to match that of u and v. If iterable, then the shape must match the shapes of u and v. + w: Wavelengths (in nanometers) at which to evaluate wavefront gradient. theta: Field angle at which to evaluate wavefront, as a 2-tuple of `galsim.Angle` instances. [default: (0.0*galsim.arcmin, 0.0*galsim.arcmin)] Only a single theta is permitted. @@ -935,9 +936,9 @@ def wavefront_gradient(self, u, v, t, theta=(0.0*radians, 0.0*radians)): Arrays dWdu and dWdv of wavefront lag or lead gradient in nm/m. """ if len(self._layers) > 1: - return np.sum([layer.wavefront_gradient(u, v, t, theta) for layer in self], axis=0) + return np.sum([layer.wavefront_gradient(u, v, t, w, theta) for layer in self], axis=0) else: - return self._layers[0].wavefront_gradient(u, v, t, theta) + return self._layers[0].wavefront_gradient(u, v, t, w, theta) def _wavefront(self, u, v, t, theta): if len(self._layers) > 1: @@ -945,10 +946,10 @@ def _wavefront(self, u, v, t, theta): else: return self._layers[0]._wavefront(u, v, t, theta) - def _wavefront_gradient(self, u, v, t, theta): - gradx, grady = self._layers[0]._wavefront_gradient(u, v, t, theta) + def _wavefront_gradient(self, u, v, t, w, theta): + gradx, grady = self._layers[0]._wavefront_gradient(u, v, t, w, theta) for layer in self._layers[1:]: - gx, gy = layer._wavefront_gradient(u, v, t, theta) + gx, gy = layer._wavefront_gradient(u, v, t, w, theta) gradx += gx grady += gy return gradx, grady @@ -1579,6 +1580,10 @@ def _shoot(self, photons, rng): u = photons.pupil_u v = photons.pupil_v t = photons.time + if photons.hasAllocatedWavelengths(): + w = photons.wavelength + else: + w = None n_photons = len(photons) @@ -1588,7 +1593,7 @@ def _shoot(self, photons, rng): nm_to_arcsec = 1.e-9 * radians / arcsec if self._fft_sign == '+': nm_to_arcsec *= -1 - photons.x, photons.y = self._screen_list._wavefront_gradient(u, v, t, self.theta) + photons.x, photons.y = self._screen_list._wavefront_gradient(u, v, t, w, self.theta) photons.x *= nm_to_arcsec photons.y *= nm_to_arcsec photons.flux = self._flux / n_photons diff --git a/galsim/phase_screens.py b/galsim/phase_screens.py index b85fa369aed..bc93f07715c 100644 --- a/galsim/phase_screens.py +++ b/galsim/phase_screens.py @@ -221,6 +221,11 @@ def work(i, atm): significantly faster than computing PSFs with non-frozen-flow atmospheres. If ``alpha`` != 1.0, then it is required that a ``time_step`` is also specified. [default: 1.0] + seeing_exp: Power law index for wavelength dependent seeing when using geometric + photon shooting. For pure Kolmogorov turbulence, the correct value is + -0.2, but values in the range [-0.2, -0.4] are expected for Von Karman + turbulence with finite outer scales. This value is ignored if drawing + in FFT mode. [default: 0.0]. time_step: Time interval between phase boiling updates. Note that this is distinct from the time interval used to integrate the PSF over time, which is set by the ``time_step`` keyword argument to `PhaseScreenPSF` or @@ -246,8 +251,8 @@ def work(i, atm): September 2014 """ def __init__(self, screen_size, screen_scale=None, altitude=0.0, r0_500=0.2, L0=25.0, - vx=0.0, vy=0.0, alpha=1.0, time_step=None, rng=None, suppress_warning=False, - mp_context=None): + vx=0.0, vy=0.0, alpha=1.0, seeing_exp=0.0, time_step=None, rng=None, + suppress_warning=False, mp_context=None): if (alpha != 1.0 and time_step is None): raise GalSimIncompatibleValuesError( "No time_step provided when alpha != 1.0", alpha=alpha, time_step=time_step) @@ -274,6 +279,7 @@ def __init__(self, screen_size, screen_scale=None, altitude=0.0, r0_500=0.2, L0= self.vx = vx self.vy = vy self.alpha = alpha + self.seeing_exp = seeing_exp if rng is None: rng = BaseDeviate() @@ -689,7 +695,7 @@ def _wavefront(self, u, v, t, theta): v += self._altitude*theta[1].tan() return self._tab2d._call_wrap(u.ravel(), v.ravel()).reshape(u.shape) - def wavefront_gradient(self, u, v, t=None, theta=(0.0*radians, 0.0*radians)): + def wavefront_gradient(self, u, v, t=None, w=None, theta=(0.0*radians, 0.0*radians)): """ Compute gradient of wavefront due to atmospheric phase screen. Parameters: @@ -702,6 +708,10 @@ def wavefront_gradient(self, u, v, t=None, theta=(0.0*radians, 0.0*radians)): will be used for all u, v. If scalar, then the size will be broadcast up to match that of u and v. If iterable, then the shape must match the shapes of u and v. [default: None] + w: Wavelength in nanometers to use to compute gradient. This is a bit of a + hack to implement chromatic seeing in photon-shooting mode. Returned gradients + are multiplied by (wavelength/500)^(self.seeing_exp). If not present, then + no chromatic seeing correction is applied. [default: None] theta: Field angle at which to evaluate wavefront, as a 2-tuple of `galsim.Angle` instances. [default: (0.0*galsim.arcmin, 0.0*galsim.arcmin)] Only a single theta is permitted. @@ -728,7 +738,7 @@ def wavefront_gradient(self, u, v, t=None, theta=(0.0*radians, 0.0*radians)): self.instantiate() # noop if already instantiated if self.reversible: - return self._wavefront_gradient(u, v, t, theta) + return self._wavefront_gradient(u, v, t, w, theta) else: dwdu = np.empty_like(u, dtype=np.float64) dwdv = np.empty_like(u, dtype=np.float64) @@ -738,12 +748,15 @@ def wavefront_gradient(self, u, v, t=None, theta=(0.0*radians, 0.0*radians)): while tt <= tmax: self._seek(tt) here = ((tt <= t) & (t < tt+self.time_step)) - dwdu[here], dwdv[here] = self._wavefront_gradient(u[here], v[here], t[here], theta) + w_here = None if w is None else w[here] + dwdu[here], dwdv[here] = self._wavefront_gradient( + u[here], v[here], t[here], w_here, theta + ) tt += self.time_step return dwdu, dwdv - def _wavefront_gradient(self, u, v, t, theta): - # Same as wavefront(), but no argument checking and no boiling updates. + def _wavefront_gradient(self, u, v, t, w, theta): + # Same as wavefront_gradient(), but no argument checking and no boiling updates. u = u - t*self.vx if theta[0].rad != 0: u += self._altitude*theta[0].tan() @@ -751,6 +764,9 @@ def _wavefront_gradient(self, u, v, t, theta): if theta[1].rad != 0: v += self._altitude*theta[1].tan() dfdx, dfdy = self._tab2d._gradient_wrap(u.ravel(), v.ravel()) + if w is not None: + dfdx *= (w/500.)**self.seeing_exp + dfdy *= (w/500.)**self.seeing_exp return dfdx.reshape(u.shape), dfdy.reshape(u.shape) @@ -1087,7 +1103,7 @@ def _wavefront(self, u, v, t, theta): # Note, this phase screen is actually independent of time and theta. return self._zernike.evalCartesian(u, v) * self.lam_0 - def wavefront_gradient(self, u, v, t=None, theta=None): + def wavefront_gradient(self, u, v, t=None, w=None, theta=None): """ Compute gradient of wavefront due to optical phase screen. Parameters: @@ -1096,6 +1112,7 @@ def wavefront_gradient(self, u, v, t=None, theta=None): v: Vertical pupil coordinate (in meters) at which to evaluate wavefront. Can be a scalar or an iterable. The shapes of u and v must match. t: Ignored for `OpticalScreen`. + w: Ignored for `OpticalScreen`. theta: Ignored for `OpticalScreen`. Returns: @@ -1105,10 +1122,10 @@ def wavefront_gradient(self, u, v, t=None, theta=None): v = np.array(v, dtype=float, copy=False) if u.shape != v.shape: raise GalSimIncompatibleValuesError("u.shape not equal to v.shape", u=u, v=v) - return self._wavefront_gradient(u, v, t, theta) + return self._wavefront_gradient(u, v, t, w, theta) - def _wavefront_gradient(self, u, v, t, theta): + def _wavefront_gradient(self, u, v, t, w, theta): # Same as wavefront(), but no argument checking. # Note, this phase screen is actually independent of time and theta. gradx, grady = self._zernike.evalCartesianGrad(u, v) @@ -1213,7 +1230,7 @@ def wavefront(self, u, v, t=None, theta=None): def _wavefront(self, u, v, t, theta): return self.table(u, v) - def wavefront_gradient(self, u, v, t=None, theta=None): + def wavefront_gradient(self, u, v, t=None, w=None, theta=None): """ Evaluate gradient of wavefront from lookup table. Parameters: @@ -1222,6 +1239,7 @@ def wavefront_gradient(self, u, v, t=None, theta=None): v: Vertical pupil coordinate (in meters) at which to evaluate wavefront. Can be a scalar or an iterable. The shapes of u and v must match. t: Ignored for `UserScreen`. + w: Ignored for `UserScreen`. theta: Ignored for `UserScreen`. Returns: @@ -1229,5 +1247,5 @@ def wavefront_gradient(self, u, v, t=None, theta=None): """ return self.table.gradient(u, v) - def _wavefront_gradient(self, u, v, t, theta): + def _wavefront_gradient(self, u, v, t, w, theta): return self.table.gradient(u, v) From fce671c2e1848b99366a5279f67224dd6336f843 Mon Sep 17 00:00:00 2001 From: Josh Meyers Date: Thu, 13 Jul 2023 23:06:38 -0700 Subject: [PATCH 2/5] test seeing_exp with const wavelength --- galsim/phase_psf.py | 2 +- galsim/phase_screens.py | 2 +- tests/test_phase_psf.py | 38 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 40 insertions(+), 2 deletions(-) diff --git a/galsim/phase_psf.py b/galsim/phase_psf.py index 76c6c8e88ba..d7900335b1e 100644 --- a/galsim/phase_psf.py +++ b/galsim/phase_psf.py @@ -1583,7 +1583,7 @@ def _shoot(self, photons, rng): if photons.hasAllocatedWavelengths(): w = photons.wavelength else: - w = None + w = self.lam n_photons = len(photons) diff --git a/galsim/phase_screens.py b/galsim/phase_screens.py index bc93f07715c..a90805b64e5 100644 --- a/galsim/phase_screens.py +++ b/galsim/phase_screens.py @@ -764,7 +764,7 @@ def _wavefront_gradient(self, u, v, t, w, theta): if theta[1].rad != 0: v += self._altitude*theta[1].tan() dfdx, dfdy = self._tab2d._gradient_wrap(u.ravel(), v.ravel()) - if w is not None: + if w is not None and self.seeing_exp is not None: dfdx *= (w/500.)**self.seeing_exp dfdy *= (w/500.)**self.seeing_exp return dfdx.reshape(u.shape), dfdy.reshape(u.shape) diff --git a/tests/test_phase_psf.py b/tests/test_phase_psf.py index 908e0b2c5c6..118f0d2ee53 100644 --- a/tests/test_phase_psf.py +++ b/tests/test_phase_psf.py @@ -1532,6 +1532,44 @@ def test_t_persistence(): assert np.max(photons.time) < 25.0 +def test_seeing_exp(): + rng = galsim.BaseDeviate(10) + atm1 = galsim.Atmosphere( + screen_size=10.0, + altitude=[0,1,2,3], + r0_500=0.2, + seeing_exp=-0.3, + rng=rng.duplicate() + ) + atm2 = galsim.Atmosphere( + screen_size=10.0, + altitude=[0,1,2,3], + r0_500=0.2, + seeing_exp=None, + rng=rng.duplicate() + ) + # When all photons are 500nm, seeing_exp has no effect. + psf1 = atm1.makePSF(exptime=15.0, lam=500.0, diam=1.0) + psf2 = atm2.makePSF(exptime=15.0, lam=500.0, diam=1.0) + nphot = 100_000 + photons1 = psf1.drawImage(save_photons=True, method='phot', n_photons=nphot, rng=rng.duplicate()).photons + photons2 = psf2.drawImage(save_photons=True, method='phot', n_photons=nphot, rng=rng.duplicate()).photons + np.testing.assert_allclose(photons1.x, photons2.x, rtol=0, atol=1e-15) + np.testing.assert_allclose(photons1.y, photons2.y, rtol=0, atol=1e-15) + + # But if lam != 500, we should see a scaling + psf1 = atm1.makePSF(exptime=15.0, lam=700.0, diam=1.0, second_kick=False) # second kick is achromatic + psf2 = atm2.makePSF(exptime=15.0, lam=700.0, diam=1.0, second_kick=False) # so turn it off + nphot = 100_000 + + img = galsim.Image(11, 11, scale=0.2) # odd makes center at (0,0) + photons1 = psf1.drawImage(image=img, save_photons=True, method='phot', n_photons=nphot, rng=rng.duplicate()).photons + photons2 = psf2.drawImage(image=img, save_photons=True, method='phot', n_photons=nphot, rng=rng.duplicate()).photons + + np.testing.assert_allclose(photons1.x, photons2.x*(700/500)**-0.3, rtol=0, atol=1e-15) + np.testing.assert_allclose(photons1.y, photons2.y*(700/500)**-0.3, rtol=0, atol=1e-15) + + if __name__ == "__main__": testfns = [v for k, v in vars().items() if k[:5] == 'test_' and callable(v)] for testfn in testfns: From 2ddb484965ec97feaf0c3681a436cc4b4293f9d0 Mon Sep 17 00:00:00 2001 From: Josh Meyers Date: Thu, 13 Jul 2023 23:29:56 -0700 Subject: [PATCH 3/5] Also test with array of wavelengths --- tests/test_phase_psf.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/tests/test_phase_psf.py b/tests/test_phase_psf.py index 118f0d2ee53..8fc0b3c6a3c 100644 --- a/tests/test_phase_psf.py +++ b/tests/test_phase_psf.py @@ -1569,6 +1569,24 @@ def test_seeing_exp(): np.testing.assert_allclose(photons1.x, photons2.x*(700/500)**-0.3, rtol=0, atol=1e-15) np.testing.assert_allclose(photons1.y, photons2.y*(700/500)**-0.3, rtol=0, atol=1e-15) + # Also try when using a range of wavelengths + sed = galsim.SED('vega.txt', 'nm', 'flambda') + op1 = galsim.WavelengthSampler(sed) + op2 = galsim.WavelengthSampler(sed) + + photons1 = psf1.drawImage( + image=img, save_photons=True, method='phot', n_photons=nphot, rng=rng.duplicate(), + photon_ops=[op1] + ).photons + photons2 = psf2.drawImage( + image=img, save_photons=True, method='phot', n_photons=nphot, rng=rng.duplicate(), + photon_ops=[op2] + ).photons + np.testing.assert_equal(photons1.wavelength, photons2.wavelength) + assert np.std(photons1.wavelength) > 0 + np.testing.assert_allclose(photons1.x, photons2.x*(700/500)**-0.3, rtol=0, atol=1e-15) + np.testing.assert_allclose(photons1.y, photons2.y*(700/500)**-0.3, rtol=0, atol=1e-15) + if __name__ == "__main__": testfns = [v for k, v in vars().items() if k[:5] == 'test_' and callable(v)] From d69b572ce53452989ee7647122864f6777426043 Mon Sep 17 00:00:00 2001 From: Mike Jarvis Date: Sat, 15 Jul 2023 21:53:33 -0400 Subject: [PATCH 4/5] Switch to new codecov uploader --- .github/workflows/ci.yml | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index ae759fcb88a..31863c1b1fd 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -181,6 +181,8 @@ jobs: cd tests coverage run -m pytest -v pytest -v run_examples.py + coverage combine + coverage xml cd .. # N.B. This seems to happen automatically if omitted. # Less confusing to include it explicitly. @@ -209,18 +211,12 @@ jobs: - name: Upload coverage to codecov if: matrix.py != 'pypy-3.7' - run: | - cd tests - pwd -P - ls -la - coverage combine || true # (Not necessary I think, but just in case.) - coverage report - ls -la - #codecov # This didn't work. - # cf. https://community.codecov.io/t/github-not-getting-codecov-report-after-switching-from-travis-to-github-actions-for-ci/ - # The solution was to switch to the bash uploader line instead. - bash <(curl -s https://codecov.io/bash) - cd .. + uses: codecov/codecov-action@v3 + with: + token: ${{ secrets.CODECOV_TOKEN }} + files: tests/coverage.xml + fail_ci_if_error: false + verbose: true - name: Pre-cache cleanup continue-on-error: true From 73c273375c29193769a94ba4ab6de9c1542ac8a2 Mon Sep 17 00:00:00 2001 From: Josh Meyers Date: Mon, 17 Jul 2023 09:29:18 -0700 Subject: [PATCH 5/5] Refactor test of array of wavelengths --- tests/test_phase_psf.py | 62 ++++++++++++++++++++++++++++++----------- 1 file changed, 45 insertions(+), 17 deletions(-) diff --git a/tests/test_phase_psf.py b/tests/test_phase_psf.py index 8fc0b3c6a3c..158762dcebc 100644 --- a/tests/test_phase_psf.py +++ b/tests/test_phase_psf.py @@ -1552,40 +1552,68 @@ def test_seeing_exp(): psf1 = atm1.makePSF(exptime=15.0, lam=500.0, diam=1.0) psf2 = atm2.makePSF(exptime=15.0, lam=500.0, diam=1.0) nphot = 100_000 - photons1 = psf1.drawImage(save_photons=True, method='phot', n_photons=nphot, rng=rng.duplicate()).photons - photons2 = psf2.drawImage(save_photons=True, method='phot', n_photons=nphot, rng=rng.duplicate()).photons + photons1 = psf1.drawImage( + save_photons=True, method='phot', n_photons=nphot, rng=rng.duplicate() + ).photons + photons2 = psf2.drawImage( + save_photons=True, method='phot', n_photons=nphot, rng=rng.duplicate() + ).photons np.testing.assert_allclose(photons1.x, photons2.x, rtol=0, atol=1e-15) np.testing.assert_allclose(photons1.y, photons2.y, rtol=0, atol=1e-15) # But if lam != 500, we should see a scaling - psf1 = atm1.makePSF(exptime=15.0, lam=700.0, diam=1.0, second_kick=False) # second kick is achromatic - psf2 = atm2.makePSF(exptime=15.0, lam=700.0, diam=1.0, second_kick=False) # so turn it off + # also, second kick is achromatic, so turn it off + psf1 = atm1.makePSF(exptime=15.0, lam=700.0, diam=1.0, second_kick=False) + psf2 = atm2.makePSF(exptime=15.0, lam=700.0, diam=1.0, second_kick=False) nphot = 100_000 img = galsim.Image(11, 11, scale=0.2) # odd makes center at (0,0) - photons1 = psf1.drawImage(image=img, save_photons=True, method='phot', n_photons=nphot, rng=rng.duplicate()).photons - photons2 = psf2.drawImage(image=img, save_photons=True, method='phot', n_photons=nphot, rng=rng.duplicate()).photons + photons1 = psf1.drawImage( + image=img, save_photons=True, method='phot', n_photons=nphot, rng=rng.duplicate() + ).photons + photons2 = psf2.drawImage( + image=img, save_photons=True, method='phot', n_photons=nphot, rng=rng.duplicate() + ).photons np.testing.assert_allclose(photons1.x, photons2.x*(700/500)**-0.3, rtol=0, atol=1e-15) np.testing.assert_allclose(photons1.y, photons2.y*(700/500)**-0.3, rtol=0, atol=1e-15) # Also try when using a range of wavelengths sed = galsim.SED('vega.txt', 'nm', 'flambda') - op1 = galsim.WavelengthSampler(sed) - op2 = galsim.WavelengthSampler(sed) - - photons1 = psf1.drawImage( - image=img, save_photons=True, method='phot', n_photons=nphot, rng=rng.duplicate(), - photon_ops=[op1] + bandpass = galsim.Bandpass('LSST_r.dat', 'nm') + star = (galsim.DeltaFunction()*sed).withFlux(nphot, bandpass) + + photons1 = star.drawImage( + bandpass, + image=img, + method='phot', + n_photons=nphot, + rng=rng.duplicate(), + save_photons=True, + photon_ops=[psf1] ).photons - photons2 = psf2.drawImage( - image=img, save_photons=True, method='phot', n_photons=nphot, rng=rng.duplicate(), - photon_ops=[op2] + photons2 = star.drawImage( + bandpass, + image=img, + method='phot', + n_photons=nphot, + rng=rng.duplicate(), + save_photons=True, + photon_ops=[psf2] ).photons + np.testing.assert_equal(photons1.wavelength, photons2.wavelength) assert np.std(photons1.wavelength) > 0 - np.testing.assert_allclose(photons1.x, photons2.x*(700/500)**-0.3, rtol=0, atol=1e-15) - np.testing.assert_allclose(photons1.y, photons2.y*(700/500)**-0.3, rtol=0, atol=1e-15) + np.testing.assert_allclose( + photons1.x, + photons2.x*(photons2.wavelength/500)**-0.3, + rtol=0, atol=1e-15 + ) + np.testing.assert_allclose( + photons1.y, + photons2.y*(photons2.wavelength/500)**-0.3, + rtol=0, atol=1e-15 + ) if __name__ == "__main__":