Skip to content

Commit

Permalink
add basic unit handling, and update test
Browse files Browse the repository at this point in the history
  • Loading branch information
@mperrin
mperrin committed Nov 25, 2024
1 parent ef0024c commit 65f7747
Show file tree
Hide file tree
Showing 2 changed files with 60 additions and 43 deletions.
24 changes: 16 additions & 8 deletions poppy/poppy_core.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3402,15 +3402,23 @@ def __init__(self, pixelscale=1 * (u.arcsec / u.pixel), fov_pixels=None, fov_arc

if offset is not None:
if len(offset) != 2:
raise ValueError("If a detector offset is specified, it must be a tuple or list with 2 elements, giving the (X, Y) offsets.")
raise ValueError("If a detector offset is specified, it must be a tuple or list with 2 elements, "
"giving the (X, Y) offsets.")
# The offset is specified in pixels, so this can have units of pixels,
# or else if an integer or float, that's considered as implicitly a number of pixels
if isinstance(offset, u.Quantity):
try:
offset = offset.to_value(u.pixel)
except u.UnitConversionError:
raise(ValueError(f"A detector offset must be specified in units of detector pixels, not '{offset.unit}'"))
offset = np.asarray(offset) # ensure it's an ndarray, not just a list or tuple
# todo, do something sensible with units. These should have units consistent with the detector pixescale without the /pixel part
# note the detector pixelscale can be in arcsec/pix or meters/pix, depending on context
# a note on sign convention for detector offset. This is regrettably confusing.
# The implementation in matrixDFT has the sense of "how much should the source be offset", i.e. an offset of +5 pix moves the source by +5 pix.
# However, physically we would like the opposite sign convention. Moving the detector by +5 pix should move the source by -5 pix.
# This is implemented by a sign flip multplication by -1 which is applied in the _propagate_mft methods.
# That could just be a hard-coded -1, but we choose to implement as a named variable to help make this logic clear later to readers of the code.
# A note on sign convention for detector offset: (This is regrettably confusing.)
# The implementation in matrixDFT has the sense of "how much should the source be offset",
# i.e. an offset of +5 pix moves the source by +5 pix.
# However, physically we would like the opposite sign convention: Moving the detector by +5 pix
# should move the source by -5 pix. This is implemented by a sign flip multplication by -1
# which is applied in the _propagate_mft methods. That could just be a hard-coded -1,
# but we choose to implement as a named variable to help make this logic clear later to readers of this code:
self.offset = offset
self._offset_sign = -1

Expand Down
79 changes: 44 additions & 35 deletions poppy/tests/test_core.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -719,43 +719,52 @@ def test_detector_offsets(plot=False, pixscale=0.01, fov_pixels=100):
In other words, shifting a source by (+dX,+dY) should look the same as
shifting the detector by (-dX, -dY)
And you can specify the detector offsets in units of pixels or just as floats.
"""
source_offset_r = .1
for offset_theta in [0, 45, 90, 180]:

# Compute offsets from radial to cartesian coords.
# recall astronomy convention is PA=0 is +Y, increasing CCW
source_offset_x = -source_offset_r * np.sin(np.deg2rad(offset_theta)) # arcsec
source_offset_y = source_offset_r * np.cos(np.deg2rad(offset_theta))
print(f"offset theta {offset_theta} is x = {source_offset_x}, y = {source_offset_y}")

# Create a PSF with a shifted source
offset_source_sys = poppy_core.OpticalSystem(npix=1024, oversample=1)
offset_source_sys.add_pupil(optics.ParityTestAperture())
offset_source_sys.add_detector(pixelscale=pixscale, fov_pixels=fov_pixels, oversample=1) #, offset=(pixscale/2, pixscale/2))
# This interface only has r, theta offsets available. Can't use _x, _y offsets here
offset_source_sys.source_offset_r = source_offset_r
offset_source_sys.source_offset_theta = offset_theta
offset_source_psf = offset_source_sys.calc_psf()

# Create a PSF with a shifted detector, the other way
offset_det_sys = poppy_core.OpticalSystem(npix=1024, oversample=1)
offset_det_sys.add_pupil(optics.ParityTestAperture())
offset_det_sys.add_detector(pixelscale=pixscale, fov_pixels=fov_pixels, oversample=1,
offset=(-source_offset_y/pixscale, -source_offset_x/pixscale)) # Y, X, in pixels
offset_det_psf = offset_det_sys.calc_psf()

if plot:
fig, axes = plt.subplots(figsize=(16,9), ncols=3)
poppy.display_psf(offset_source_psf, ax=axes[0], crosshairs=True, colorbar=False,
title=f'Offset Source: {source_offset_x:.3f} arcsec, {source_offset_y:.3f}')
poppy.display_psf(offset_det_psf, ax=axes[1], crosshairs=True, colorbar=False,
title=f'Offset Det: {offset_det_sys.planes[-1].offset[1]:.2f}, {offset_det_sys.planes[-1].offset[0]:.2f} pix')
poppy.display_psf_difference(offset_source_psf, offset_det_psf, ax=axes[2],
title='difference', colorbar=False)

# Check the equality of the two results from the two above
assert np.allclose(offset_source_psf[0].data, offset_det_psf[0].data), "Offset source and offset detector the opposite way should be equivalent"
for with_units in [True, False]:
for offset_theta in [0, 45, 90, 180]:

# Compute offsets from radial to cartesian coords.
# recall astronomy convention is PA=0 is +Y, increasing CCW
source_offset_x = -source_offset_r * np.sin(np.deg2rad(offset_theta)) # arcsec
source_offset_y = source_offset_r * np.cos(np.deg2rad(offset_theta))
print(f"offset theta {offset_theta} is x = {source_offset_x}, y = {source_offset_y}")

# Create a PSF with a shifted source
offset_source_sys = poppy_core.OpticalSystem(npix=1024, oversample=1)
offset_source_sys.add_pupil(optics.ParityTestAperture())
offset_source_sys.add_detector(pixelscale=pixscale, fov_pixels=fov_pixels, oversample=1) #, offset=(pixscale/2, pixscale/2))
# This interface only has r, theta offsets available. Can't use _x, _y offsets here
offset_source_sys.source_offset_r = source_offset_r
offset_source_sys.source_offset_theta = offset_theta
offset_source_psf = offset_source_sys.calc_psf()

# Create a PSF with a shifted detector, the other way
offset_det_sys = poppy_core.OpticalSystem(npix=1024, oversample=1)
offset_det_sys.add_pupil(optics.ParityTestAperture())

det_offset = (-source_offset_y/pixscale, -source_offset_x/pixscale) # Y, X in pixels
if with_units:
det_offset = np.asarray(det_offset) * u.pixel
offset_det_sys.add_detector(pixelscale=pixscale, fov_pixels=fov_pixels, oversample=1,
offset=det_offset)
offset_det_psf = offset_det_sys.calc_psf()

if plot:
fig, axes = plt.subplots(figsize=(16,9), ncols=3)
poppy.display_psf(offset_source_psf, ax=axes[0], crosshairs=True, colorbar=False,
title=f'Offset Source: {source_offset_x:.3f} arcsec, {source_offset_y:.3f}')
poppy.display_psf(offset_det_psf, ax=axes[1], crosshairs=True, colorbar=False,
title=f'Offset Det: {offset_det_sys.planes[-1].offset[1]:.2f}, {offset_det_sys.planes[-1].offset[0]:.2f} pix')
poppy.display_psf_difference(offset_source_psf, offset_det_psf, ax=axes[2],
title='difference', colorbar=False)

# Check the equality of the two results from the two above
assert np.allclose(offset_source_psf[0].data, offset_det_psf[0].data), "Offset source and offset detector the opposite way should be equivalent"




# Tests for CompoundOpticalSystem
Expand Down

0 comments on commit 65f7747

Please sign in to comment.