diff --git a/src/kbmod/reprojection_utils.py b/src/kbmod/reprojection_utils.py new file mode 100644 index 000000000..c118e5867 --- /dev/null +++ b/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 `_ + """ + 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 diff --git a/tests/test_reprojection_utils.py b/tests/test_reprojection_utils.py new file mode 100644 index 000000000..9cecbad45 --- /dev/null +++ b/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)