Skip to content

Commit

Permalink
Fix astrometric_utils tests.
Browse files Browse the repository at this point in the history
  • Loading branch information
@mairanteodoro
mairanteodoro committed Oct 11, 2023
1 parent 5bc4560 commit 75b7857
Showing 1 changed file with 189 additions and 59 deletions.
248 changes: 189 additions & 59 deletions romancal/tweakreg/tests/test_astrometric_utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,12 +24,160 @@
get_catalog,
)

ARAD = np.pi / 180.0


class MockConnectionError:
def __init__(self, *args, **kwargs):
raise requests.exceptions.ConnectionError


def get_earth_sun_orbit_elem(epoch):
amjd = (epoch - 2012.0) * 365.25 + 55927.0
imjd = int(amjd)
mjd = float(imjd)

time_argument = mjd - 51544.5
fmean_longitude = 280.460 + 0.9856474 * time_argument
fmean_anomaly = 357.528 + 0.9856003 * time_argument
iremainder = int(fmean_longitude / 360.0)
fmean_longitude = fmean_longitude - (iremainder * 360.0)
if fmean_longitude < 0:
fmean_longitude = fmean_longitude + 360.0
iremainder = int(fmean_anomaly / 360.0)
fmean_anomaly = fmean_anomaly - (iremainder * 360.0)
if fmean_anomaly < 0:
fmean_anomaly = fmean_anomaly + 360.0
ecliptic_longitude = (
fmean_longitude
+ 1.915 * np.sin(fmean_anomaly * ARAD)
+ 0.02 * np.sin(2.0 * fmean_anomaly * ARAD)
)
ecliptic_obliquity = 23.439 - 4.0 * 0.0000001 * time_argument

# radius vector (in AU)
# (approximation using an algorithm from USNO)
earth_sun_distance = (
1.00014
- 0.01671 * np.cos(fmean_anomaly * ARAD)
- 0.00014 * np.cos(2.0 * fmean_anomaly * ARAD)
)

earth_sun_distance = u.Quantity(earth_sun_distance, unit="AU")

return {
"earth_sun_distance": earth_sun_distance,
"ecliptic_longitude": ecliptic_longitude,
"ecliptic_obliquity": ecliptic_obliquity,
}


def get_parallax_correction_mast(epoch, gaia_ref_epoch_coords):
orbit_elem = get_earth_sun_orbit_elem(epoch)
earth_sun_distance = orbit_elem["earth_sun_distance"]
ecliptic_longitude = orbit_elem["ecliptic_longitude"]
ecliptic_obliquity = orbit_elem["ecliptic_obliquity"]

# cartesian coordinates of the sun
xsun = earth_sun_distance * np.cos(ecliptic_longitude * ARAD)
ysun = (
earth_sun_distance
* np.cos(ecliptic_obliquity * ARAD)
* np.sin(ecliptic_longitude * ARAD)
)
zsun = (
earth_sun_distance
* np.sin(ecliptic_obliquity * ARAD)
* np.sin(ecliptic_longitude * ARAD)
)

# angular displacement components
delta_ra = (
u.Quantity(gaia_ref_epoch_coords["parallax"], unit="mas")
/ np.cos(gaia_ref_epoch_coords["dec"] * ARAD)
* (
-xsun.value * np.sin(gaia_ref_epoch_coords["ra"] * ARAD)
+ ysun.value * np.cos(gaia_ref_epoch_coords["ra"] * ARAD)
)
).to("deg")
delta_dec = (
u.Quantity(gaia_ref_epoch_coords["parallax"], unit="mas")
* (
-xsun.value
* np.cos(gaia_ref_epoch_coords["ra"] * ARAD)
* np.sin(gaia_ref_epoch_coords["dec"] * ARAD)
- ysun.value
* np.sin(gaia_ref_epoch_coords["ra"] * ARAD)
* np.sin(gaia_ref_epoch_coords["dec"] * ARAD)
+ zsun.value * np.cos(gaia_ref_epoch_coords["dec"] * ARAD)
)
).to("deg")

return delta_ra, delta_dec


def get_parallax_correction_barycenter(epoch, gaia_ref_epoch_coords):
"""
Calculates the parallax correction in the Earth barycenter frame for a given epoch
and Gaia reference epoch coordinates (i.e. Gaia coordinates at the reference epoch).
Parameters
----------
epoch : float
The epoch for which the parallax correction is calculated.
gaia_ref_epoch_coords : dict
The Gaia reference epoch coordinates, including 'ra', 'dec', and 'parallax'.
Returns
-------
tuple
A tuple containing the delta_ra and delta_dec values of the parallax correction
in degrees.
Examples
--------
.. code-block :: python
epoch = 2022.5
gaia_coords = {'ra': 180.0, 'dec': 45.0, 'parallax': 10.0}
correction = get_parallax_correction_earth_barycenter(epoch, gaia_coords)
print(correction)
(0.001, -0.002)
""" # noqa: E501

obs_date = Time(epoch, format="decimalyear")
earths_center_barycentric_coords = coord.get_body_barycentric(
"earth", obs_date, ephemeris="builtin"
)
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)
delta_ra = (
u.Quantity(gaia_ref_epoch_coords["parallax"], unit="mas").to(u.rad)
* (1 / np.cos(gaia_ref_epoch_coords["dec"] * ARAD))
* (
earth_X.value * np.sin(gaia_ref_epoch_coords["ra"] * ARAD)
- earth_Y.value * np.cos(gaia_ref_epoch_coords["ra"] * ARAD)
)
).to("deg")
delta_dec = (
u.Quantity(gaia_ref_epoch_coords["parallax"], unit="mas").to(u.rad)
* (
earth_X.value
* np.cos(gaia_ref_epoch_coords["ra"] * ARAD)
* np.sin(gaia_ref_epoch_coords["dec"] * ARAD)
+ earth_Y.value
* np.sin(gaia_ref_epoch_coords["ra"] * ARAD)
* np.sin(gaia_ref_epoch_coords["dec"] * ARAD)
- earth_Z.value * np.cos(gaia_ref_epoch_coords["dec"] * ARAD)
)
).to("deg")

return delta_ra, delta_dec


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
Expand Down Expand Up @@ -88,62 +236,22 @@ def get_proper_motion_correction(epoch, gaia_ref_epoch_coords, gaia_ref_epoch):


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)
"""
# get parallax correction using textbook calculations (i.e. Earth's barycenter)
parallax_corr = get_parallax_correction_barycenter(
epoch=epoch, gaia_ref_epoch_coords=gaia_ref_epoch_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
# get parallax using the same coordinates frame as MAST
# (i.e. Sun's geocentric coordinates)
parallax_corr_mast = get_parallax_correction_mast(
epoch=epoch, gaia_ref_epoch_coords=gaia_ref_epoch_coords
)

# 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")
# add parallax corrections columns to the main table
gaia_ref_epoch_coords["parallax_delta_ra"] = parallax_corr[0]
gaia_ref_epoch_coords["parallax_delta_dec"] = parallax_corr[1]
gaia_ref_epoch_coords["parallax_delta_ra_mast"] = parallax_corr_mast[0]
gaia_ref_epoch_coords["parallax_delta_dec_mast"] = parallax_corr_mast[1]


def update_wcsinfo(input_dm):
Expand Down Expand Up @@ -198,6 +306,7 @@ def create_wcs_for_tweakreg_pipeline(input_dm, shift_1=0, shift_2=0):

# create required frames
detector = cf.Frame2D(name="detector", axes_order=(0, 1), unit=(u.pix, u.pix))
detector = cf.Frame2D(name="detector", axes_order=(0, 1), unit=(u.pix, u.pix))
v2v3 = cf.Frame2D(
name="v2v3",
axes_order=(0, 1),
Expand Down Expand Up @@ -455,6 +564,7 @@ def test_create_astrometric_catalog_write_results_to_disk(tmp_path, base_image):
("GAIADR3", None),
],
)
def test_create_astrometric_catalog_using_epoch(tmp_path, catalog, epoch, request):
def test_create_astrometric_catalog_using_epoch(tmp_path, catalog, epoch, request):
"""Test fetching data from supported catalogs for a specific epoch."""
output_filename = "ref_cat.ecsv"
Expand Down Expand Up @@ -591,12 +701,11 @@ def test_get_catalog_valid_parameters_but_no_sources_returned():
(0, 0, 2000),
(0, 0, 2010.3),
(0, 0, 2030),
(269.4521, 4.6933, 2030),
(89, 80, 2010),
],
)
def test_get_catalog_using_epoch(ra, dec, epoch):
"""Test that get_catalog returns coordinates corrected by proper motion."""
"""Test that get_catalog returns coordinates corrected by proper motion
and parallax."""

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

Expand All @@ -606,6 +715,8 @@ def test_get_catalog_using_epoch(ra, dec, epoch):
mask = result_all["parallax_over_error"] > 5

result = result_all[mask]

# updated coordinates at the provided epoch
returned_ra = np.array(result["ra"])
returned_dec = np.array(result["dec"])

Expand All @@ -623,9 +734,12 @@ def test_get_catalog_using_epoch(ra, dec, epoch):
)
# calculate parallax corrections
get_parallax_correction(epoch=epoch, gaia_ref_epoch_coords=gaia_ref_epoch_coords)
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

# textbook (barycentric frame)
expected_ra = (
gaia_ref_epoch_coords["ra"]
+ gaia_ref_epoch_coords["pm_delta_ra"]
Expand All @@ -637,10 +751,26 @@ def test_get_catalog_using_epoch(ra, dec, epoch):
+ gaia_ref_epoch_coords["parallax_delta_dec"]
)

# mast (geocentric frame)
expected_ra_mast = (
gaia_ref_epoch_coords["ra"]
+ gaia_ref_epoch_coords["pm_delta_ra"]
+ gaia_ref_epoch_coords["parallax_delta_ra_mast"]
)
expected_dec_mast = (
gaia_ref_epoch_coords["dec"]
+ gaia_ref_epoch_coords["pm_delta_dec"]
+ gaia_ref_epoch_coords["parallax_delta_dec_mast"]
)

assert len(result) > 0

assert np.isclose(returned_ra, expected_ra, atol=1e-5, rtol=0).all()
assert np.isclose(returned_dec, expected_dec, atol=1e-5, rtol=0).all()
# N.B.: atol=1e-8 (in deg) corresponds to a coordinate difference of ~ 40 uas
assert np.isclose(expected_ra, returned_ra, atol=1e-8, rtol=0).all()
assert np.isclose(expected_dec, returned_dec, atol=1e-8, rtol=0).all()

assert np.isclose(expected_ra_mast, returned_ra, atol=5e-10, rtol=0).all()
assert np.isclose(expected_dec_mast, returned_dec, atol=5e-10, rtol=0).all()


def test_get_catalog_timeout():
Expand Down

0 comments on commit 75b7857

Please sign in to comment.