Skip to content

Commit

Permalink
ecef2eci: scalar output for scalar input with Numpy backend
Browse files Browse the repository at this point in the history
add test for Numpy and AstroPy backends
  • Loading branch information
scivision committed Jan 8, 2025
1 parent bf34608 commit 39e262d
Show file tree
Hide file tree
Showing 2 changed files with 69 additions and 64 deletions.
103 changes: 41 additions & 62 deletions src/pymap3d/eci.py
Original file line number Diff line number Diff line change
Expand Up @@ -54,25 +54,7 @@ def eci2ecef_astropy(x, y, z, t: datetime) -> tuple:
"""
eci2ecef using Astropy
Parameters
----------
x : float
ECI x-location [meters]
y : float
ECI y-location [meters]
z : float
ECI z-location [meters]
t : datetime.datetime
time of obsevation (UTC)
Results
-------
x_ecef : float
x ECEF coordinate
y_ecef : float
y ECEF coordinate
z_ecef : float
z ECEF coordinate
see eci2ecef() for description
"""

gcrs = GCRS(CartesianRepresentation(x * u.m, y * u.m, z * u.m), obstime=t)
Expand All @@ -89,27 +71,7 @@ def eci2ecef_numpy(x, y, z, t: datetime) -> tuple:
"""
eci2ecef using Numpy
less accurate than astropy, but may be good enough.
Parameters
----------
x : float
ECI x-location [meters]
y : float
ECI y-location [meters]
z : float
ECI z-location [meters]
t : datetime.datetime
time of obsevation (UTC)
Results
-------
x_ecef : float
x ECEF coordinate
y_ecef : float
y ECEF coordinate
z_ecef : float
z ECEF coordinate
see eci2ecef() for description
"""

x = numpy.atleast_1d(x)
Expand Down Expand Up @@ -165,33 +127,50 @@ def ecef2eci(x, y, z, time: datetime) -> tuple:
"""

try:
itrs = ITRS(CartesianRepresentation(x * u.m, y * u.m, z * u.m), obstime=time)
gcrs = itrs.transform_to(GCRS(obstime=time))
eci = EarthLocation(*gcrs.cartesian.xyz)

x_eci = eci.x.value
y_eci = eci.y.value
z_eci = eci.z.value
return ecef2eci_astropy(x, y, z, time)
except NameError:
x = numpy.atleast_1d(x)
y = numpy.atleast_1d(y)
z = numpy.atleast_1d(z)
gst = numpy.atleast_1d(greenwichsrt(juliandate(time)))
assert x.shape == y.shape == z.shape
assert x.size == gst.size

ecef = numpy.column_stack((x.ravel(), y.ravel(), z.ravel()))
eci = numpy.empty((x.size, 3))
for i in range(x.size):
eci[i, :] = R3(gst[i]).T @ ecef[i, :]

x_eci = eci[:, 0].reshape(x.shape)
y_eci = eci[:, 1].reshape(y.shape)
z_eci = eci[:, 2].reshape(z.shape)
return ecef2eci_numpy(x, y, z, time)


def ecef2eci_astropy(x, y, z, t: datetime) -> tuple:
"""ecef2eci using Astropy
see ecef2eci() for description
"""
itrs = ITRS(CartesianRepresentation(x * u.m, y * u.m, z * u.m), obstime=t)
gcrs = itrs.transform_to(GCRS(obstime=t))
eci = EarthLocation(*gcrs.cartesian.xyz)

x_eci = eci.x.value
y_eci = eci.y.value
z_eci = eci.z.value

return x_eci, y_eci, z_eci


def ecef2eci_numpy(x, y, z, t: datetime) -> tuple:
"""ecef2eci using Numpy
see ecef2eci() for description
"""

x = numpy.atleast_1d(x)
y = numpy.atleast_1d(y)
z = numpy.atleast_1d(z)
gst = numpy.atleast_1d(greenwichsrt(juliandate(t)))
assert x.shape == y.shape == z.shape
assert x.size == gst.size

ecef = numpy.column_stack((x.ravel(), y.ravel(), z.ravel()))
eci = numpy.empty((x.size, 3))
for i in range(x.size):
eci[i, :] = R3(gst[i]).T @ ecef[i, :]

x_eci = eci[:, 0].reshape(x.shape)
y_eci = eci[:, 1].reshape(y.shape)
z_eci = eci[:, 2].reshape(z.shape)

return x_eci.squeeze()[()], y_eci.squeeze()[()], z_eci.squeeze()[()]


def R3(x: float):
"""Rotation matrix for ECI"""
return numpy.array(
Expand Down
30 changes: 28 additions & 2 deletions src/pymap3d/tests/test_eci.py
Original file line number Diff line number Diff line change
Expand Up @@ -55,9 +55,35 @@ def test_ecef2eci():
# this example from Matlab ecef2eci docs
eci = pm.ecef2eci(*ECEF, UTC)

rel = 0.01 if astropy is None else 0.0001
assert isinstance(eci[0], float)
assert isinstance(eci[1], float)
assert isinstance(eci[2], float)


def test_ecef2eci_numpy():
pytest.importorskip("numpy")

eci = pm.eci.ecef2eci_numpy(*ECEF, UTC)

rel = 0.025

assert eci == approx(ECI, rel=rel)
assert isinstance(eci[0], float)
assert isinstance(eci[1], float)
assert isinstance(eci[2], float)


def test_ecef2eci_astropy():
pytest.importorskip("astropy")

eci = pm.eci.ecef2eci_astropy(*ECEF, UTC)

rel = 0.0001

assert eci == approx([-2981810.6, 5207039.5, 3161595.1], rel=rel)
assert eci == approx(ECI, rel=rel)
assert isinstance(eci[0], float)
assert isinstance(eci[1], float)
assert isinstance(eci[2], float)


def test_eci2geodetic():
Expand Down

0 comments on commit 39e262d

Please sign in to comment.