|
9 | 9 |
|
10 | 10 | from astropy.time import Time |
11 | 11 | import astropy |
| 12 | +import numpy as np |
| 13 | + |
12 | 14 | import pymap3d as pm |
13 | 15 |
|
| 16 | +try: |
| 17 | + from pymap3d.tests.matlab_engine import matlab_ecef2eci, matlab_engine, has_matmap3d |
| 18 | +except ImportError: |
| 19 | + matlab_engine = None |
| 20 | + |
14 | 21 | print("Python version:", sys.version) |
15 | 22 | print("AstroPy version:", astropy.__version__) |
16 | 23 |
|
|
28 | 35 | print("Julian Date (TT):", astropy_time.tt.jd) |
29 | 36 | print("GMST:", astropy_time.sidereal_time("mean", "greenwich")) |
30 | 37 |
|
| 38 | +np.set_printoptions(precision=3, suppress=True) |
| 39 | + |
31 | 40 | # %% 1. Geodetic to ECEF |
32 | 41 | ecef = pm.geodetic2ecef(lat, lon, alt) |
33 | 42 | print("\nECEF Coordinates (meters):") |
34 | | -print(f"X: {ecef[0]:.8f}, Y: {ecef[1]:.8f}, Z: {ecef[2]:.8f}") |
| 43 | +print(np.array(ecef)) |
| 44 | + |
| 45 | +# %% AstroPy ECEF to ECI (J2000) |
| 46 | +astropy_eci = pm.ecef2eci(ecef[0], ecef[1], ecef[2], dt) |
| 47 | +astropy_eci = np.array(astropy_eci) |
| 48 | +print("\nAstroPy: ECI Coordinates (meters):") |
| 49 | +print(astropy_eci) |
| 50 | + |
| 51 | +numpy_eci = pm.ecef2eci(ecef[0], ecef[1], ecef[2], dt, force_non_astropy=True) |
| 52 | +numpy_eci = np.array(numpy_eci) |
| 53 | +print("\nNumpy: ECI Coordinates (meters):") |
| 54 | +print(numpy_eci) |
| 55 | + |
| 56 | +print("\nAstroPy - Numpy Difference (ECI meters):", astropy_eci - numpy_eci) |
| 57 | + |
| 58 | +# %% Matlab comparison |
| 59 | +if matlab_engine is not None: |
| 60 | + eng = matlab_engine() |
| 61 | + eci_matlab_aerospace = matlab_ecef2eci(eng, False, dt, ecef) |
| 62 | + print("\nMatlab Aerospace Toolbox: ECI Coordinates (meters):") |
| 63 | + print(np.array(eci_matlab_aerospace)) |
| 64 | + print( |
| 65 | + "\nAstroPy - Matlab Aerospace Toolbox Difference (ECI meters):", |
| 66 | + astropy_eci - eci_matlab_aerospace, |
| 67 | + ) |
| 68 | + print( |
| 69 | + "Numpy - Matlab Aerospace Toolbox Difference (ECI meters):", |
| 70 | + numpy_eci - eci_matlab_aerospace, |
| 71 | + ) |
35 | 72 |
|
36 | | -# %% 2. ECEF to ECI (J2000) |
37 | | -eci = pm.ecef2eci(ecef[0], ecef[1], ecef[2], dt) |
38 | | -print("\nECI Coordinates (meters):") |
39 | | -print(f"X: {eci[0]}, Y: {eci[1]}, Z: {eci[2]}") |
| 73 | + if has_matmap3d(eng): |
| 74 | + eci_matmap3d = matlab_ecef2eci(eng, True, dt, ecef) |
| 75 | + print("\nMatlab Matmap3D: ECI Coordinates (meters):") |
| 76 | + print(np.array(eci_matmap3d)) |
| 77 | + print("\nAstroPy - Matlab Matmap3D Difference (ECI meters):", astropy_eci - eci_matmap3d) |
| 78 | + print("Numpy - Matlab Matmap3D Difference (ECI meters):", numpy_eci - eci_matmap3d) |
0 commit comments