add correct_parallax

src/kbmod/reprojection_utils.py

@@ -0,0 +1,57 @@
+import astropy.units as u
+import numpy as np
+from astropy import units as u
+from astropy.coordinates import GCRS, ICRS
+from scipy.optimize import minimize
+
+
+def correct_parallax(coord, obstime, point_on_earth, guess_distance):
+    """Calculate the parallax corrected postions for a given object at a given time and distance from Earth.
+
+    Attributes
+    ----------
+    coord : `astropy.coordinate.SkyCoord`
+        The coordinate to be corrected for.
+    obstime : `astropy.time.Time` or `string`
+        The observation time.
+    point_on_earth : `astropy.coordinate.EarthLocation`
+        The location on Earth of the observation.
+    guess_distance : `float`
+        The guess distance to the object from Earth.
+
+    Returns
+    ----------
+    An `astropy.coordinate.SkyCoord` containing the ra and dec of the pointin ICRS.
+
+    References
+    ----------
+    .. [1] `Jupyter Notebook <https://github.com/DinoBektesevic/region_search_example/blob/main/02_accounting_parallax.ipynb>`_
+    """
+    loc = (
+        point_on_earth.x.to(u.m).value,
+        point_on_earth.y.to(u.m).value,
+        point_on_earth.z.to(u.m).value,
+    ) * u.m
+
+    # line of sight from earth to the object,
+    # the object has an unknown distance from earth
+    los_earth_obj = coord.transform_to(GCRS(obstime=obstime, obsgeoloc=loc))
+
+    cost = lambda d: np.abs(
+        guess_distance
+        - GCRS(ra=los_earth_obj.ra, dec=los_earth_obj.dec, distance=d * u.AU, obstime=obstime, obsgeoloc=loc)
+        .transform_to(ICRS())
+        .distance.to(u.AU)
+        .value
+    )
+
+    fit = minimize(
+        cost,
+        (guess_distance,),
+    )
+
+    answer = GCRS(
+        ra=los_earth_obj.ra, dec=los_earth_obj.dec, distance=fit.x[0] * u.AU, obstime=obstime, obsgeoloc=loc
+    ).transform_to(ICRS())
+
+    return answer

tests/test_reprojection_utils.py

@@ -0,0 +1,47 @@
+import unittest
+
+import numpy.testing as npt
+from astropy.coordinates import EarthLocation, SkyCoord, solar_system_ephemeris
+from astropy.time import Time
+
+from kbmod.reprojection_utils import correct_parallax
+
+
+class test_reprojection_utils(unittest.TestCase):
+    def test_parallax_equinox(self):
+        icrs_ra1 = 88.74513571
+        icrs_dec1 = 23.43426475
+        time1 = Time("2023-03-20T16:00:00", format="isot", scale="utc")
+
+        icrs_ra2 = 91.24261107
+        icrs_dec2 = 23.43437467
+        time2 = Time("2023-09-24T04:00:00", format="isot", scale="utc")
+
+        sc1 = SkyCoord(ra=icrs_ra1, dec=icrs_dec1, unit="deg")
+        sc2 = SkyCoord(ra=icrs_ra2, dec=icrs_dec2, unit="deg")
+
+        with solar_system_ephemeris.set("de432s"):
+            loc = EarthLocation.of_site("ctio")
+
+        corrected_coord1 = correct_parallax(
+            coord=sc1,
+            obstime=time1,
+            point_on_earth=loc,
+            guess_distance=50.0,
+        )
+
+        expected_ra = 90.0
+        expected_dec = 23.43952556
+
+        npt.assert_almost_equal(corrected_coord1.ra.value, expected_ra)
+        npt.assert_almost_equal(corrected_coord1.dec.value, expected_dec)
+
+        corrected_coord2 = correct_parallax(
+            coord=sc2,
+            obstime=time2,
+            point_on_earth=loc,
+            guess_distance=50.0,
+        )
+
+        npt.assert_almost_equal(corrected_coord2.ra.value, expected_ra)
+        npt.assert_almost_equal(corrected_coord2.dec.value, expected_dec)