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