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.

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

Include parallax correction in unit tests. #925

Closed
Changes from all commits
Commits
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
173 changes: 151 additions & 22 deletions romancal/tweakreg/tests/test_astrometric_utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
from astropy import units as u
from astropy.modeling import models
from astropy.modeling.models import RotationSequence3D, Scale, Shift
from astropy.time import Time
from gwcs import coordinate_frames as cf
from gwcs import wcs
from gwcs.geometry import CartesianToSpherical, SphericalToCartesian
Expand All @@ -29,6 +30,122 @@ def __init__(self, *args, **kwargs):
raise requests.exceptions.ConnectionError


def get_proper_motion_correction(epoch, gaia_ref_epoch_coords, gaia_ref_epoch):
"""
Calculates the proper motion correction for a given epoch and Gaia reference epoch
coordinates.

Parameters
----------
epoch : float
The epoch for which the proper motion correction is calculated.
gaia_ref_epoch_coords : dict
A dictionary containing Gaia reference epoch coordinates.
gaia_ref_epoch : float
The Gaia reference epoch.

Returns
-------
None

Examples
--------
.. code-block:: python
epoch = 2022.5
gaia_coords = {
"ra": 180.0,
"dec": 45.0,
"pmra": 2.0,
"pmdec": 1.5
}
gaia_ref_epoch = 2020.0
get_proper_motion_correction(epoch, gaia_coords, gaia_ref_epoch)
"""

expected_new_dec = (
np.array(
gaia_ref_epoch_coords["dec"] * 3600
+ (epoch - gaia_ref_epoch) * gaia_ref_epoch_coords["pmdec"] / 1000
)
/ 3600
)
average_dec = np.array(
[
np.mean([new, old])
for new, old in zip(expected_new_dec, gaia_ref_epoch_coords["dec"])
]
)
pmra = gaia_ref_epoch_coords["pmra"] / np.cos(np.deg2rad(average_dec))

# angular displacement components
gaia_ref_epoch_coords["pm_delta_dec"] = u.Quantity(
(epoch - gaia_ref_epoch) * gaia_ref_epoch_coords["pmdec"] / 1000,
unit=u.arcsec,
).to(u.deg)
gaia_ref_epoch_coords["pm_delta_ra"] = u.Quantity(
(epoch - gaia_ref_epoch) * (pmra / 1000), unit=u.arcsec
).to(u.deg)


def get_parallax_correction(epoch, gaia_ref_epoch_coords):
"""
Calculates the parallax correction for a given epoch and Gaia reference epoch
coordinates.
Calculations based on Chapter 8 of "Spherical Astronomy" by Robin M. Green.

Parameters
----------
epoch : float
The epoch for which the parallax correction is calculated.
gaia_ref_epoch_coords : dict
A dictionary containing Gaia reference epoch coordinates.

Returns
-------
None

Examples
--------
.. code-block :: python
epoch = 2022.5
gaia_coords = {
"ra": 180.0,
"dec": 45.0,
"parallax": 10.0
}
get_parallax_correction(epoch, gaia_coords)
"""

obs_date = Time(epoch, format="decimalyear")
earths_center_barycentric_coords = coord.get_body_barycentric("earth", obs_date)
earth_X = earths_center_barycentric_coords.x
earth_Y = earths_center_barycentric_coords.y
earth_Z = earths_center_barycentric_coords.z

# angular displacement components
# (see eq. 8.15 of "Spherical Astronomy" by Robert M. Green)
gaia_ref_epoch_coords["parallax_delta_ra"] = (
u.Quantity(gaia_ref_epoch_coords["parallax"], unit="mas").to(u.rad)
* (1 / np.cos(gaia_ref_epoch_coords["dec"] / 180 * np.pi))
* (
earth_X.value * np.sin(gaia_ref_epoch_coords["ra"] / 180 * np.pi)
- earth_Y.value * np.cos(gaia_ref_epoch_coords["ra"] / 180 * np.pi)
)
).to("deg")
gaia_ref_epoch_coords["parallax_delta_dec"] = (
u.Quantity(gaia_ref_epoch_coords["parallax"], unit="mas").to(u.rad)
* (
earth_X.value
* np.cos(gaia_ref_epoch_coords["ra"] / 180 * np.pi)
* np.sin(gaia_ref_epoch_coords["dec"] / 180 * np.pi)
+ earth_Y.value
* np.sin(gaia_ref_epoch_coords["ra"] / 180 * np.pi)
* np.sin(gaia_ref_epoch_coords["dec"] / 180 * np.pi)
- earth_Z.value * np.cos(gaia_ref_epoch_coords["dec"] / 180 * np.pi)
)
).to("deg")


def update_wcsinfo(input_dm):
"""
Update WCSInfo with realistic data (i.e. obtained from romanisim simulations).
Expand Down Expand Up @@ -481,39 +598,51 @@ def test_get_catalog_valid_parameters_but_no_sources_returned():
def test_get_catalog_using_epoch(ra, dec, epoch):
"""Test that get_catalog returns coordinates corrected by proper motion."""

result = get_catalog(ra, dec, epoch=epoch)
result_all = get_catalog(ra, dec, epoch=epoch)

# select sources with reliable astrometric solutions based on the
# parallax_over_error parameter as discussed in Fabricius et al. 2021
# (https://www.aanda.org/articles/aa/full_html/2021/05/aa39834-20/aa39834-20.html)
mask = result_all["parallax_over_error"] > 5

result = result_all[mask]
returned_ra = np.array(result["ra"])
returned_dec = np.array(result["dec"])

# get GAIA data and update coords to requested epoch using pm measurements
gaia_ref_epoch = 2016.0
gaia_ref_epoch_coords = get_catalog(ra, dec, epoch=gaia_ref_epoch)
gaia_ref_epoch_coords_all = get_catalog(ra, dec, epoch=gaia_ref_epoch)

expected_new_dec = (
np.array(
gaia_ref_epoch_coords["dec"] * 3600
+ (epoch - gaia_ref_epoch) * gaia_ref_epoch_coords["pmdec"] / 1000
)
/ 3600
gaia_ref_epoch_coords = gaia_ref_epoch_coords_all[mask]

# calculate proper motion corrections
get_proper_motion_correction(
epoch=epoch,
gaia_ref_epoch_coords=gaia_ref_epoch_coords,
gaia_ref_epoch=gaia_ref_epoch,
)
average_dec = np.array(
[
np.mean([new, old])
for new, old in zip(expected_new_dec, gaia_ref_epoch_coords["dec"])
]
# calculate parallax corrections
get_parallax_correction(epoch=epoch, gaia_ref_epoch_coords=gaia_ref_epoch_coords)

# calculate the expected coordinates value after corrections have been applied to
# Gaia's reference epoch coordinates
expected_ra = (
gaia_ref_epoch_coords["ra"]
+ gaia_ref_epoch_coords["pm_delta_ra"]
+ gaia_ref_epoch_coords["parallax_delta_ra"]
)
pmra = gaia_ref_epoch_coords["pmra"] / np.cos(np.deg2rad(average_dec))
expected_new_ra = (
np.array(
gaia_ref_epoch_coords["ra"] * 3600
+ (epoch - gaia_ref_epoch) * (pmra / 1000)
)
/ 3600
expected_dec = (
gaia_ref_epoch_coords["dec"]
+ gaia_ref_epoch_coords["pm_delta_dec"]
+ gaia_ref_epoch_coords["parallax_delta_dec"]
)

assert len(result) > 0
assert np.isclose(returned_ra, expected_new_ra, atol=1e-10, rtol=1e-9).all()
assert np.isclose(returned_dec, expected_new_dec, atol=1e-10, rtol=1e-9).all()
# choosing atol=0 corresponds to abs(a - b) / abs(b) <= rtol,
# where a is the returned value from the VO API and b is the expected
# value from applying the calculated PM and parallax corrections.
assert np.isclose(returned_ra, expected_ra, atol=0, rtol=1e-5).all()
assert np.isclose(returned_dec, expected_dec, atol=0, rtol=1e-5).all()


def test_get_catalog_timeout():
Expand Down