From 787812287c7dce46f92cd29a7c4acabc0faf9f44 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Mon, 24 Jun 2019 15:25:28 -0400 Subject: [PATCH] Python: calling GeoMoon() without errors, but calculations not yet validated. --- generate/template/astronomy.py | 55 ++++++++++++++++++++++++---------- generate/test.py | 11 +++++++ generate/unit_test_python | 3 +- source/python/astronomy.py | 55 ++++++++++++++++++++++++---------- 4 files changed, 93 insertions(+), 31 deletions(-) diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index 76e6bae4..ff36f0b8 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -42,8 +42,15 @@ def _NormalizeLongitude(lon): lon -= 360.0 return lon -def VectorLength(vector): - return math.sqrt(vector.x**2 + vector.y**2 + vector.z**2) +class Vector: + def __init__(self, x, y, z, t): + self.x = x + self.y = y + self.z = z + self.t = t + + def Length(self): + return math.sqrt(self.x**2 + self.y**2 + self.z**2) BODY_INVALID = -1 BODY_MERCURY = 0 @@ -116,7 +123,7 @@ def _SynodicPeriod(body): return abs(_EARTH_ORBITAL_PERIOD / (_EARTH_ORBITAL_PERIOD/_PlanetOrbitalPeriod[body] - 1.0)) def _AngleBetween(a, b): - r = VectorLength(a) * VectorLength(b) + r = a.Length() * b.Length() if r < 1.0e-8: return BadVectorError() dot = (a.x*b.x + a.y*b.y + a.z*b.z) / r @@ -395,7 +402,7 @@ class _e_tilt: self.tt = time.tt self.ee = e.dpsi * math.cos(self.mobl * _DEG2RAD) / 15.0 -def _ecl2equ_vec(time, ecl, equ): +def _ecl2equ_vec(time, ecl): obl = _mean_obliq(time.tt) * _DEG2RAD cos_obl = math.cos(obl) sin_obl = math.sin(obl) @@ -603,7 +610,7 @@ def _ter2cel(time, vec1): # BEGIN CalcMoon class _Array1: - def __init__(xmin, xmax): + def __init__(self, xmin, xmax): self.min = xmin self.array = [0] * (xmax - xmin + 1) @@ -614,7 +621,7 @@ class _Array1: self.array[key - self.min] = value class _Array2: - def __init__(xmin, xmax, ymin, ymax): + def __init__(self, xmin, xmax, ymin, ymax): self.min = xmin self.array = [_Array1(ymin, ymax) for i in range(xmax - xmin + 1)] @@ -665,11 +672,11 @@ def _CalcMoon(time): -539E-9 * Sine(0.35498-5.37899*T) -64E-9 * Sine(0.39943-5.37511*T))) - L0 = PI2*Frac(0.60643382+1336.85522467*T-0.00000313*T2) + DL0/ARC - L = PI2*Frac(0.37489701+1325.55240982*T+0.00002565*T2) + DL /ARC - LS = PI2*Frac(0.99312619+ 99.99735956*T-0.00000044*T2) + DLS/ARC - F = PI2*Frac(0.25909118+1342.22782980*T-0.00000892*T2) + DF /ARC - D = PI2*Frac(0.82736186+1236.85308708*T-0.00000397*T2) + DD /ARC + L0 = _PI2*Frac(0.60643382+1336.85522467*T-0.00000313*T2) + DL0/_ARC + L = _PI2*Frac(0.37489701+1325.55240982*T+0.00002565*T2) + DL /_ARC + LS = _PI2*Frac(0.99312619+ 99.99735956*T-0.00000044*T2) + DLS/_ARC + F = _PI2*Frac(0.25909118+1342.22782980*T-0.00000892*T2) + DF /_ARC + D = _PI2*Frac(0.82736186+1236.85308708*T-0.00000397*T2) + DD /_ARC I = 1 while I <= 4: @@ -702,7 +709,7 @@ def _CalcMoon(time): def Term(p, q, r, s): result = (1, 0) - I = [null, p, q, r, s] + I = [None, p, q, r, s] k = 1 while k <= 4: if I[k] != 0: @@ -848,14 +855,32 @@ def _CalcMoon(time): +0.24*Sine(0.2275 -5.7374*T)+0.28*Sine(0.2965 +2.6929*T) +0.33*Sine(0.3132 +6.3368*T) ) - S = F + DS/ARC + S = F + DS/_ARC lat_seconds = (1.000002708 + 139.978*DGAM)*(18518.511+1.189+GAM1C)*math.sin(S) - 6.24*math.sin(3*S) + N return _moonpos( - _PI2 * Frac((L0+DLAM/ARC) / _PI2), + _PI2 * Frac((L0+DLAM/_ARC) / _PI2), (math.pi / (180 * 3600)) * lat_seconds, - (ARC * (ERAD / AU)) / (0.999953253 * SINPI) + (_ARC * (_ERAD / _AU)) / (0.999953253 * SINPI) ) # END CalcMoon #---------------------------------------------------------------------------- +def GeoMoon(time): + m = _CalcMoon(time) + + # Convert geocentric ecliptic spherical coordinates to Cartesian coordinates. + dist_cos_lat = m.distance_au * math.cos(m.geo_eclip_lat) + gepos = [ + dist_cos_lat * math.cos(m.geo_eclip_lon), + dist_cos_lat * math.sin(m.geo_eclip_lon), + m.distance_au * math.sin(m.geo_eclip_lat) + ] + + # Convert ecliptic coordinates to equatorial coordinates, both in mean equinox of date. + mpos1 = _ecl2equ_vec(time, gepos) + + # Convert from mean equinox of date to J2000. + mpos2 = _precession(time.tt, mpos1, 0) + return Vector(mpos2[0], mpos2[1], mpos2[2], time) + diff --git a/generate/test.py b/generate/test.py index b824ea37..1e7c9de6 100644 --- a/generate/test.py +++ b/generate/test.py @@ -32,10 +32,21 @@ def Test_AstroTime(): sys.exit(1) print('Current time =', astronomy.Time.Now()) + +def Test_GeoMoon(): + time = astronomy.Time.Make(2019, 6, 24, 15, 45, 37) + vec = astronomy.GeoMoon(time) + print('Test_GeoMoon: vec = {:0.10f}, {:0.10f}, {:0.10f}'.format(vec.x, vec.y, vec.z)) + + if len(sys.argv) == 2: if sys.argv[1] == 'time': Test_AstroTime() sys.exit(0) + if sys.argv[1] == 'moon': + Test_GeoMoon() + sys.exit(0) + print('test.py: Invalid command line arguments.') sys.exit(1) diff --git a/generate/unit_test_python b/generate/unit_test_python index 8b1f602e..be265c01 100755 --- a/generate/unit_test_python +++ b/generate/unit_test_python @@ -6,5 +6,6 @@ Fail() } python3 --version || Fail "Cannot print python version" -python3 test.py time || Fail "Failure reported by test.py" +python3 test.py time || Fail "Failure reported by test.py (time)" +python3 test.py moon || Fail "Failure reported by test.py (moon)" exit 0 diff --git a/source/python/astronomy.py b/source/python/astronomy.py index 03c1a92a..41c05a0d 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -42,8 +42,15 @@ def _NormalizeLongitude(lon): lon -= 360.0 return lon -def VectorLength(vector): - return math.sqrt(vector.x**2 + vector.y**2 + vector.z**2) +class Vector: + def __init__(self, x, y, z, t): + self.x = x + self.y = y + self.z = z + self.t = t + + def Length(self): + return math.sqrt(self.x**2 + self.y**2 + self.z**2) BODY_INVALID = -1 BODY_MERCURY = 0 @@ -116,7 +123,7 @@ def _SynodicPeriod(body): return abs(_EARTH_ORBITAL_PERIOD / (_EARTH_ORBITAL_PERIOD/_PlanetOrbitalPeriod[body] - 1.0)) def _AngleBetween(a, b): - r = VectorLength(a) * VectorLength(b) + r = a.Length() * b.Length() if r < 1.0e-8: return BadVectorError() dot = (a.x*b.x + a.y*b.y + a.z*b.z) / r @@ -486,7 +493,7 @@ class _e_tilt: self.tt = time.tt self.ee = e.dpsi * math.cos(self.mobl * _DEG2RAD) / 15.0 -def _ecl2equ_vec(time, ecl, equ): +def _ecl2equ_vec(time, ecl): obl = _mean_obliq(time.tt) * _DEG2RAD cos_obl = math.cos(obl) sin_obl = math.sin(obl) @@ -694,7 +701,7 @@ def _ter2cel(time, vec1): # BEGIN CalcMoon class _Array1: - def __init__(xmin, xmax): + def __init__(self, xmin, xmax): self.min = xmin self.array = [0] * (xmax - xmin + 1) @@ -705,7 +712,7 @@ class _Array1: self.array[key - self.min] = value class _Array2: - def __init__(xmin, xmax, ymin, ymax): + def __init__(self, xmin, xmax, ymin, ymax): self.min = xmin self.array = [_Array1(ymin, ymax) for i in range(xmax - xmin + 1)] @@ -756,11 +763,11 @@ def _CalcMoon(time): -539E-9 * Sine(0.35498-5.37899*T) -64E-9 * Sine(0.39943-5.37511*T))) - L0 = PI2*Frac(0.60643382+1336.85522467*T-0.00000313*T2) + DL0/ARC - L = PI2*Frac(0.37489701+1325.55240982*T+0.00002565*T2) + DL /ARC - LS = PI2*Frac(0.99312619+ 99.99735956*T-0.00000044*T2) + DLS/ARC - F = PI2*Frac(0.25909118+1342.22782980*T-0.00000892*T2) + DF /ARC - D = PI2*Frac(0.82736186+1236.85308708*T-0.00000397*T2) + DD /ARC + L0 = _PI2*Frac(0.60643382+1336.85522467*T-0.00000313*T2) + DL0/_ARC + L = _PI2*Frac(0.37489701+1325.55240982*T+0.00002565*T2) + DL /_ARC + LS = _PI2*Frac(0.99312619+ 99.99735956*T-0.00000044*T2) + DLS/_ARC + F = _PI2*Frac(0.25909118+1342.22782980*T-0.00000892*T2) + DF /_ARC + D = _PI2*Frac(0.82736186+1236.85308708*T-0.00000397*T2) + DD /_ARC I = 1 while I <= 4: @@ -793,7 +800,7 @@ def _CalcMoon(time): def Term(p, q, r, s): result = (1, 0) - I = [null, p, q, r, s] + I = [None, p, q, r, s] k = 1 while k <= 4: if I[k] != 0: @@ -939,14 +946,32 @@ def _CalcMoon(time): +0.24*Sine(0.2275 -5.7374*T)+0.28*Sine(0.2965 +2.6929*T) +0.33*Sine(0.3132 +6.3368*T) ) - S = F + DS/ARC + S = F + DS/_ARC lat_seconds = (1.000002708 + 139.978*DGAM)*(18518.511+1.189+GAM1C)*math.sin(S) - 6.24*math.sin(3*S) + N return _moonpos( - _PI2 * Frac((L0+DLAM/ARC) / _PI2), + _PI2 * Frac((L0+DLAM/_ARC) / _PI2), (math.pi / (180 * 3600)) * lat_seconds, - (ARC * (ERAD / AU)) / (0.999953253 * SINPI) + (_ARC * (_ERAD / _AU)) / (0.999953253 * SINPI) ) # END CalcMoon #---------------------------------------------------------------------------- +def GeoMoon(time): + m = _CalcMoon(time) + + # Convert geocentric ecliptic spherical coordinates to Cartesian coordinates. + dist_cos_lat = m.distance_au * math.cos(m.geo_eclip_lat) + gepos = [ + dist_cos_lat * math.cos(m.geo_eclip_lon), + dist_cos_lat * math.sin(m.geo_eclip_lon), + m.distance_au * math.sin(m.geo_eclip_lat) + ] + + # Convert ecliptic coordinates to equatorial coordinates, both in mean equinox of date. + mpos1 = _ecl2equ_vec(time, gepos) + + # Convert from mean equinox of date to J2000. + mpos2 = _precession(time.tt, mpos1, 0) + return Vector(mpos2[0], mpos2[1], mpos2[2], time) +