Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Chromatic seeing when photon shooting #1225

Closed
wants to merge 5 commits into from
Closed

Chromatic seeing when photon shooting #1225

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
851e979
Add seeing_exp param to AtmosphericScreen
jmeyers314 Jul 10, 2023
fce671c
test seeing_exp with const wavelength
jmeyers314 Jul 14, 2023
2ddb484
Also test with array of wavelengths
jmeyers314 Jul 14, 2023
d69b572
Switch to new codecov uploader
rmjarvis Jul 16, 2023
73c2733
Refactor test of array of wavelengths
jmeyers314 Jul 17, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 8 additions & 12 deletions .github/workflows/ci.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand Down Expand Up @@ -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
Expand Down
19 changes: 12 additions & 7 deletions galsim/phase_psf.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -914,7 +914,7 @@
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:
Expand All @@ -927,6 +927,7 @@
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.
Expand All @@ -935,20 +936,20 @@
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:
return np.sum([layer._wavefront(u, v, t, theta) for layer in self], axis=0)
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
Expand Down Expand Up @@ -1579,6 +1580,10 @@
u = photons.pupil_u
v = photons.pupil_v
t = photons.time
if photons.hasAllocatedWavelengths():
w = photons.wavelength

Check warning on line 1584 in galsim/phase_psf.py

View check run for this annotation

Codecov / codecov/patch

galsim/phase_psf.py#L1584 

Added line #L1584 was not covered by tests
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Need a test that uses this feature. (Which is the feature request that triggered this PR...) Maybe one that compares this implementation to one using ChromaticAtmosphere?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, I had naïvely thought that adding a WavelengthSampler photon op to drawImage would trigger this, but I guess those ops only get executed after the base object is drawn. Looking into how to actually trigger this now...

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's a little bit of a design conundrum even: PhaseScreenPSF is currently a GSObject (and not a ChromaticObject). That's perfectly fine for the original usage, and would continue to be fine if seeing_exp = 0. But the profile is chromatic, of course, if seeing_exp != 0. Duck-typing might work here, but may require changing some of the code around https://github.com/GalSim-developers/GalSim/blob/main/galsim/gsobject.py#L419. For instance, the chromatic version of this profile is not separable.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should we just be using ChromaticAtmosphere then?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay, I figured out how to trigger this branch without a big rewrite. I had to use:

star = galsim.DeltaFunction()*sed
psf = atm.makePSF(...)
star.drawImage(..., photon_ops=[psf])

rather than

obj = galsim.Convolve(star, psf)
obj.drawImage(...)

I think it's a bit of a problem that the latter doesn't work though.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmmm... Maybe? Let me try.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay, that does appear to be close. The main difference is that ChromaticAtmosphere applies its shift to the entire base object at once. So we couldn't specify an exponent per screen. This is probably good enough for the atmospheric screens, but won't be correct for any OpticalScreens and I don't think will be correct for the second_kick component of the atmosphere (which has a limiting case of an Airy, e.g.).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we do something a little custom in ChromaticAtmosphere to make it work right when the base is a PhaseScreenPSF?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah. That or tweak https://github.com/LSSTDESC/imSim/blob/main/imsim/atmPSF.py#L83 to build 3 separate things (atm+opt+2kick) instead of embedding them all inside a single PhaseScreenPSF. Then we can use ChromaticAtmosphere to only add chromatic seeing to the right part.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The latter is probably simpler...

else:
w = self.lam

n_photons = len(photons)

Expand All @@ -1588,7 +1593,7 @@
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
Expand Down
42 changes: 30 additions & 12 deletions galsim/phase_screens.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We had talked about deriving this value from L0. Is that possible? If so, it would make it easier for users who may be using a VonKarman profile and want to get this to match up correctly.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, my memory for how these are related was faulty. There's a whole slide deck here showing some attempts to infer a relation, but I'm not confident this can or should be used directly. Safer to just let the user specify a value I think.

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
Expand All @@ -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)
Expand All @@ -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()
Expand Down Expand Up @@ -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:
Expand All @@ -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.
Expand All @@ -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)
Expand All @@ -738,19 +748,25 @@ 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()
v = v - t*self.vy
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 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)


Expand Down Expand Up @@ -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:
Expand All @@ -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:
Expand All @@ -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)
Expand Down Expand Up @@ -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:
Expand All @@ -1222,12 +1239,13 @@ 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:
Arrays dWdu and dWdv of wavefront lag or lead gradient in nm/m.
"""
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)
56 changes: 56 additions & 0 deletions tests/test_phase_psf.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1532,6 +1532,62 @@ 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)

# 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)]
for testfn in testfns:
Expand Down
Loading