Skip to content

Commit

Permalink
eci2ecef: broadcast_to scalar time
Browse files Browse the repository at this point in the history
  • Loading branch information
scivision committed Nov 20, 2023
1 parent 5f27d75 commit a90b7b2
Show file tree
Hide file tree
Showing 2 changed files with 28 additions and 23 deletions.
13 changes: 7 additions & 6 deletions scripts/compare/compare_ecef2eci.py.orig
Original file line number Diff line number Diff line change
Expand Up @@ -26,12 +26,13 @@ else:
eci_matlab = eng.matmap3d.ecef2eci(utc_matlab, *ecef, nargout=3)

<<<<<<< HEAD
eci_good = np.isclose(eci_py, eci_matlab, rtol=0.01).all()
=======
eci_m = np.array(eci_matlab.tomemoryview()).reshape(eci_matlab.size, order="F")

eci_good = np.isclose(eci_py, eci_m, rtol=0.01).all()
>>>>>>> 3b1c510 (f)
=======
assert eci_matlab.size == 3
eci_good = np.isclose(eci_py, eci_matlab, rtol=0.01).all()
>>>>>>> 249a291 (eci2ecef: broadcast_to scalar time)

if eci_good:
print("OK: PyMap3d ecef2eci vs. Matlab ecef2eci")
Expand All @@ -42,10 +43,10 @@ ecef_py = eci2ecef(eci_m[0], eci_m[1], eci_m[2], utc)

if matlab_aerospace(eng):
<<<<<<< HEAD
ecef_matlab = eng.eci2ecef(utc_matlab, *eci_matlab, nargout=3)
=======
ecef_matlab = eng.eci2ecef(utc_matlab, eci_m, nargout=1)
>>>>>>> 3b1c510 (f)
=======
ecef_matlab = eng.eci2ecef(utc_matlab, eci_matlab, nargout=3)
>>>>>>> 249a291 (eci2ecef: broadcast_to scalar time)
else:
ecef_matlab = eng.matmap3d.eci2ecef(utc_matlab, *eci_m, nargout=3)

Expand Down
38 changes: 21 additions & 17 deletions src/pymap3d/eci.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@

from datetime import datetime

from numpy import array, atleast_1d, column_stack, cos, empty, sin
import numpy as np

try:
import astropy.units as u
Expand Down Expand Up @@ -52,15 +52,19 @@ def eci2ecef(x, y, z, time: datetime) -> tuple:
y_ecef = itrs.y.value
z_ecef = itrs.z.value
except NameError:
x = atleast_1d(x)
y = atleast_1d(y)
z = atleast_1d(z)
gst = atleast_1d(greenwichsrt(juliandate(time)))
assert x.shape == y.shape == z.shape, f"shape mismatch: x: ${x.shape} y: {y.shape} z: {z.shape}"
assert x.size == gst.size

eci = column_stack((x.ravel(), y.ravel(), z.ravel()))
ecef = empty((x.size, 3))
x = np.atleast_1d(x)
y = np.atleast_1d(y)
z = np.atleast_1d(z)
gst = np.atleast_1d(greenwichsrt(juliandate(time)))
assert (
x.shape == y.shape == z.shape
), f"shape mismatch: x: ${x.shape} y: {y.shape} z: {z.shape}"
if gst.size == 1 and x.size != 1:
gst = np.broadcast_to(gst, x.shape[0])
assert x.size == gst.size, f"shape mismatch: x: {x.shape} gst: {gst.shape}"

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

Expand Down Expand Up @@ -108,15 +112,15 @@ def ecef2eci(x, y, z, time: datetime) -> tuple:
y_eci = eci.y.value
z_eci = eci.z.value
except NameError:
x = atleast_1d(x)
y = atleast_1d(y)
z = atleast_1d(z)
gst = atleast_1d(greenwichsrt(juliandate(time)))
x = np.atleast_1d(x)
y = np.atleast_1d(y)
z = np.atleast_1d(z)
gst = np.atleast_1d(greenwichsrt(juliandate(time)))
assert x.shape == y.shape == z.shape
assert x.size == gst.size

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

Expand All @@ -129,4 +133,4 @@ def ecef2eci(x, y, z, time: datetime) -> tuple:

def R3(x: float):
"""Rotation matrix for ECI"""
return array([[cos(x), sin(x), 0], [-sin(x), cos(x), 0], [0, 0, 1]])
return np.array([[np.cos(x), np.sin(x), 0], [-np.sin(x), np.cos(x), 0], [0, 0, 1]])

0 comments on commit a90b7b2

Please sign in to comment.