Skip to content

Commit

Permalink
Merge pull request #484 from welucas2/u/welucas2/correct-photon-op-tests
Browse files Browse the repository at this point in the history 
Ensure self-consistent use of sky_pos and local_wcs in photon op tests.
  • Loading branch information
@welucas2
welucas2 authored Jan 17, 2025
2 parents a4c2d97 + 4b3f933 commit a96345c
Showing 1 changed file with 81 additions and 28 deletions.
109 changes: 81 additions & 28 deletions tests/test_photon_ops.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,22 @@ def create_test_photon_array(t=0.0, n_photons=10000):
time=np.full(n_photons, t),
)

def create_test_img_wcs(boresight, rottelpos=np.pi / 3 * galsim.radians):
telescope = create_test_telescope(rottelpos)
wcs_factory = BatoidWCSFactory(
boresight,
obstime=Time.strptime("2022-08-06 06:50:59.337600", "%Y-%m-%d %H:%M:%S.%f"),
telescope=telescope,
wavelength=622.3195217611445, # nm
camera=get_camera(),
temperature=280.0,
pressure=72.7,
H2O_pressure=1.0,
)
det_name = "R22_S11"
camera = get_camera()[det_name]
wcs = wcs_factory.getWCS(det=camera)
return wcs

def create_test_rubin_optics(**kwargs):
return photon_ops.RubinOptics(**create_test_rubin_optics_kwargs(**kwargs))
Expand All @@ -73,17 +89,17 @@ def create_test_rubin_optics(**kwargs):
def create_test_rubin_optics_kwargs(
boresight=galsim.CelestialCoord(0.543 * galsim.radians, -0.174 * galsim.radians),
icrf_to_field=None,
sky_pos=galsim.CelestialCoord(0.543 * galsim.radians, -0.174 * galsim.radians),
image_pos=galsim.PositionD(809.6510740536025, 3432.6477953336625),
rottelpos=np.pi / 3 * galsim.radians,
):
det_name = "R22_S11"
if icrf_to_field is None:
icrf_to_field = create_test_icrf_to_field(boresight, det_name=det_name)
img_wcs = create_test_img_wcs(boresight, rottelpos)
return dict(
telescope=create_test_telescope(rottelpos),
boresight=boresight,
sky_pos=sky_pos,
sky_pos=img_wcs.toWorld(image_pos),
image_pos=image_pos,
icrf_to_field=icrf_to_field,
det_name=det_name,
Expand Down Expand Up @@ -140,7 +156,6 @@ def create_test_rubin_diffraction_optics(
boresight,
icrf_to_field,
image_pos=image_pos,
sky_pos=sky_pos,
rottelpos=rottelpos,
),
rubin_diffraction=rubin_diffraction,
Expand All @@ -153,17 +168,19 @@ def create_test_rng():

def test_rubin_optics() -> None:
"""Check that the image of a star is contained in a disc."""

rubin_optics = create_test_rubin_optics(rottelpos=0.0 * galsim.radians)
boresight = galsim.CelestialCoord(0.543 * galsim.radians, -0.174 * galsim.radians)
image_pos = galsim.PositionD(809.6510740536025, 3432.6477953336625)
rottelpos = 0.0 * galsim.radians
rubin_optics = create_test_rubin_optics(boresight=boresight, image_pos=image_pos, rottelpos=rottelpos)
photon_array = create_test_photon_array()
local_wcs = create_test_wcs()
local_wcs = create_test_img_wcs(boresight, rottelpos).local(image_pos)
u = photon_array.pupil_u.copy()
v = photon_array.pupil_v.copy()
rubin_optics.applyTo(photon_array, local_wcs=local_wcs, rng=create_test_rng())
np.testing.assert_array_equal(photon_array.pupil_u, u)
np.testing.assert_array_equal(photon_array.pupil_v, v)
expected_x_pic_center = 564.5
expected_y_pic_center = -1431.4
expected_x_pic_center = -960.28
expected_y_pic_center = -308.09
expected_r_pic_center = 20.0
np.testing.assert_array_less(
np.hypot(
Expand All @@ -176,20 +193,28 @@ def test_rubin_optics() -> None:

def test_rubin_diffraction_produces_spikes() -> None:
"""Checks that we have spike photons and that the spkies form a cross."""
boresight = galsim.CelestialCoord(0.543 * galsim.radians, -0.174 * galsim.radians)
image_pos = galsim.PositionD(809.6510740536025, 3432.6477953336625)
rottelpos = 0.0 * galsim.radians
img_wcs = create_test_img_wcs(boresight, rottelpos)
sky_pos = img_wcs.toWorld(image_pos)
rubin_diffraction_optics = create_test_rubin_diffraction_optics(
rottelpos=0.0 * galsim.radians
boresight=boresight,
sky_pos=sky_pos,
image_pos=image_pos,
rottelpos=rottelpos
)
photon_array = create_test_photon_array(n_photons=1000000)
local_wcs = create_test_wcs()
local_wcs = img_wcs.local(image_pos)
rubin_diffraction_optics.applyTo(
photon_array, local_wcs=local_wcs, rng=create_test_rng()
)

# The expected image is contained in a disc + spikes outside the disc:
spike_angles = extract_spike_angles(
photon_array,
x_center=564.5,
y_center=-1431.4,
x_center=-960.28,
y_center=-308.09,
r=20.0,
)

Expand Down Expand Up @@ -231,13 +256,29 @@ def test_rubin_diffraction_optics_is_same_as_diffraction_and_optics() -> None:
"""Checks that the result of applying RubinDiffraction and then RubinOptics
is the same as applying the combined photon op RubinDiffractionOptics."""
photon_array_combined = create_test_photon_array(n_photons=100000)
local_wcs = create_test_wcs()
rubin_diffraction_optics = create_test_rubin_diffraction_optics()
boresight = galsim.CelestialCoord(0.543 * galsim.radians, -0.174 * galsim.radians)
image_pos = galsim.PositionD(809.6510740536025, 3432.6477953336625)
rottelpos = np.pi / 3 * galsim.radians
img_wcs = create_test_img_wcs(boresight, rottelpos)
local_wcs = img_wcs.local(image_pos)
sky_pos = img_wcs.toWorld(image_pos)
rubin_diffraction_optics = create_test_rubin_diffraction_optics(
boresight=boresight,
image_pos=image_pos,
rottelpos=rottelpos,
sky_pos=sky_pos
)
rubin_diffraction_optics.applyTo(
photon_array_combined, local_wcs=local_wcs, rng=create_test_rng()
)
rubin_diffraction = create_test_rubin_diffraction()
rubin_optics = create_test_rubin_optics()
rubin_diffraction = create_test_rubin_diffraction(
sky_pos=sky_pos
)
rubin_optics = create_test_rubin_optics(
boresight=boresight,
image_pos=image_pos,
rottelpos=rottelpos
)
photon_array_modular = create_test_photon_array(n_photons=100000)
rubin_diffraction.applyTo(
photon_array_modular, local_wcs=local_wcs, rng=create_test_rng()
Expand Down Expand Up @@ -276,31 +317,43 @@ def test_rubin_diffraction_shows_field_rotation() -> None:
latitude = -30.24463 * degrees
azimuth = 45.0 * degrees
altitude = 89.9 * degrees
boresight = galsim.CelestialCoord(0.543 * galsim.radians, -0.174 * galsim.radians)
image_pos = galsim.PositionD(809.6510740536025, 3432.6477953336625)
rottelpos = 0.0 * galsim.radians
img_wcs = create_test_img_wcs(boresight, rottelpos)
sky_pos = img_wcs.toWorld(image_pos)
rubin_diffraction_optics = create_test_rubin_diffraction_optics(
latitude, azimuth, altitude, rottelpos=0.0 * galsim.radians
latitude,
azimuth,
altitude,
boresight=boresight,
sky_pos=sky_pos,
image_pos=image_pos,
rottelpos=rottelpos
)
dt = 1.0
photon_array_0 = create_test_photon_array(t=0.0, n_photons=1000000)
photon_array_1 = create_test_photon_array(t=dt, n_photons=1000000)
local_wcs = create_test_wcs()
local_wcs = img_wcs.local(image_pos)
rubin_diffraction_optics.applyTo(
photon_array_0, local_wcs=local_wcs, rng=create_test_rng()
)
rubin_diffraction_optics.applyTo(
photon_array_1, local_wcs=local_wcs, rng=create_test_rng()
)

expected_center = np.array([-960.28, -308.09])
# The expected image is contained in a disc + spikes outside the disc:
spike_angles_0 = extract_spike_angles(
photon_array_0,
x_center=564.5,
y_center=-1431.4,
x_center=expected_center[0],
y_center=expected_center[1],
r=20.0,
)
spike_angles_1 = extract_spike_angles(
photon_array_1,
x_center=564.5,
y_center=-1431.4,
x_center=expected_center[0],
y_center=expected_center[1],
r=20.0,
)

Expand Down Expand Up @@ -579,6 +632,7 @@ def test_config_rubin_optics():
"""Check the config interface to RubinOptics."""

image_pos = galsim.PositionD(3076.4462608524213, 1566.4896702703757)
boresight = galsim.CelestialCoord(1.1047934165124105 * galsim.radians, -0.5261230452954583 * galsim.radians)
config = {
**deepcopy(TEST_BASE_CONFIG),
"image_pos": image_pos, # This would get set appropriately during normal config processing.
Expand All @@ -588,20 +642,19 @@ def test_config_rubin_optics():
"type": "RubinOptics",
"camera": "LsstCam",
"det_name": "R22_S11",
"boresight": {
"type": "RADec",
"ra": "1.1047934165124105 radians",
"dec": "-0.5261230452954583 radians",
},
"boresight": boresight,
},
]
},
"sky_pos": create_test_img_wcs(
boresight=boresight,
rottelpos=np.pi / 3 * galsim.radians
).toWorld(image_pos),
}
galsim.config.ProcessInput(config)
galsim.config.input.SetupInputsForImage(config, None)
[photon_op] = galsim.config.BuildPhotonOps(config["stamp"], "photon_ops", config)
reference_op = create_test_rubin_optics(
sky_pos=TEST_BASE_CONFIG["sky_pos"],
image_pos=image_pos,
icrf_to_field=TEST_BASE_CONFIG["_icrf_to_field"],
boresight=photon_op.boresight,
Expand Down

0 comments on commit a96345c

Please sign in to comment.