1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
|
#!/usr/bin/env python3
"""
https://github.com/geospace-code/pymap3d/issues/103
"""
from datetime import datetime
import sys
from astropy.time import Time
import astropy
import numpy as np
import pymap3d as pm
try:
from pymap3d.tests.matlab_engine import matlab_ecef2eci, matlab_engine, has_matmap3d
except ImportError:
pass
print("Python version:", sys.version)
print("AstroPy version:", astropy.__version__)
lat = 33.6 # deg
lon = 134.3 # deg
alt = 0 # m
dt = datetime(2020, 8, 14, 0, 0, 41)
# %% Astropy time for comparison
astropy_time = Time(dt, scale="utc")
print("----------------------------------------")
print("Astropy Time (UTC):", astropy_time.utc)
print("Julian Date (UTC):", astropy_time.utc.jd)
print("Julian Date (TT):", astropy_time.tt.jd)
print("GMST:", astropy_time.sidereal_time("mean", "greenwich"))
np.set_printoptions(precision=3, suppress=True)
# %% 1. Geodetic to ECEF
ecef = pm.geodetic2ecef(lat, lon, alt)
print("\nECEF Coordinates (meters):")
print(np.array(ecef))
# %% AstroPy ECEF to ECI (J2000)
astropy_eci = np.array(pm.ecef2eci(ecef[0], ecef[1], ecef[2], dt))
print("\nAstroPy: ECI Coordinates (meters):")
print(astropy_eci)
numpy_eci = np.array(pm.ecef2eci(ecef[0], ecef[1], ecef[2], dt, force_non_astropy=True))
print("\nNumpy: ECI Coordinates (meters):")
print(numpy_eci)
print("\nAstroPy - Numpy Difference (ECI meters):", astropy_eci - numpy_eci)
# %% Matlab comparison
if matlab_engine in sys.modules:
eng = matlab_engine()
eci_matlab_aerospace = matlab_ecef2eci(eng, False, dt, ecef)
print("\nMatlab Aerospace Toolbox: ECI Coordinates (meters):")
print(np.array(eci_matlab_aerospace))
print(
"\nAstroPy - Matlab Aerospace Toolbox Difference (ECI meters):",
astropy_eci - eci_matlab_aerospace,
)
print(
"Numpy - Matlab Aerospace Toolbox Difference (ECI meters):",
numpy_eci - eci_matlab_aerospace,
)
if has_matmap3d(eng):
eci_matmap3d = matlab_ecef2eci(eng, True, dt, ecef)
print("\nMatlab Matmap3D: ECI Coordinates (meters):")
print(np.array(eci_matmap3d))
print("\nAstroPy - Matlab Matmap3D Difference (ECI meters):", astropy_eci - eci_matmap3d)
print("Numpy - Matlab Matmap3D Difference (ECI meters):", numpy_eci - eci_matmap3d)
|
