From f85025da310d1fd3d7e0e6fb9bc6931cc581b21f Mon Sep 17 00:00:00 2001 From: Don Cross Date: Wed, 26 Jun 2019 21:12:27 -0400 Subject: [PATCH] Implemented astro_check in Python. Exercises major parts of the code. There is still a slight discrepancy in calculations of altitude,azimuth that is larger than between JS and C. --- generate/template/astronomy.py | 121 +++++++++++++++++++++++++++++++-- generate/test.py | 9 ++- generate/unit_test_python | 8 ++- source/python/astronomy.py | 121 +++++++++++++++++++++++++++++++-- 4 files changed, 247 insertions(+), 12 deletions(-) diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index 405060ff..a25bd43e 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -363,11 +363,11 @@ _cls_t = [ class _iau2000b: def __init__(self, time): t = time.tt / 36525 - el = ((485868.249036 + t * 1717915923.2178) % ASEC360) * ASEC2RAD - elp = ((1287104.79305 + t * 129596581.0481) % ASEC360) * ASEC2RAD - f = ((335779.526232 + t * 1739527262.8478) % ASEC360) * ASEC2RAD - d = ((1072260.70369 + t * 1602961601.2090) % ASEC360) * ASEC2RAD - om = ((450160.398036 - t * 6962890.5431) % ASEC360) * ASEC2RAD + el = ((485868.249036 + t * 1717915923.2178) % _ASEC360) * _ASEC2RAD + elp = ((1287104.79305 + t * 129596581.0481) % _ASEC360) * _ASEC2RAD + f = ((335779.526232 + t * 1739527262.8478) % _ASEC360) * _ASEC2RAD + d = ((1072260.70369 + t * 1602961601.2090) % _ASEC360) * _ASEC2RAD + om = ((450160.398036 - t * 6962890.5431) % _ASEC360) * _ASEC2RAD dp = 0 de = 0 i = 76 @@ -1143,3 +1143,114 @@ def GeoVector(body, time, aberration): raise Error('Light-travel time solver did not converge: dt={}'.format(dt)) + +def Equator(body, time, observer, ofdate, aberration): + gc_observer = _geo_pos(time, observer) + gc = GeoVector(body, time, aberration) + j2000 = [ + gc.x - gc_observer[0], + gc.y - gc_observer[1], + gc.z - gc_observer[2] + ] + if not ofdate: + return _vector2radec(j2000) + temp = _precession(0, j2000, time.tt) + datevect = _nutation(time, 0, temp) + return _vector2radec(datevect) + +REFRACTION_NONE = 0 +REFRACTION_NORMAL = 1 +REFRACTION_JPLHOR = 2 + +class HorizontalCoordinates: + def __init__(self, azimuth, altitude, ra, dec): + self.azimuth = azimuth + self.altitude = altitude + self.ra = ra + self.dec = dec + +def Horizon(time, observer, ra, dec, refraction): + if not (REFRACTION_NONE <= refraction <= REFRACTION_JPLHOR): + raise Error('Invalid refraction type: ' + str(refraction)) + + sinlat = math.sin(observer.latitude * _DEG2RAD) + coslat = math.cos(observer.latitude * _DEG2RAD) + sinlon = math.sin(observer.longitude * _DEG2RAD) + coslon = math.cos(observer.longitude * _DEG2RAD) + sindc = math.sin(dec * _DEG2RAD) + cosdc = math.cos(dec * _DEG2RAD) + sinra = math.sin(ra * 15 * _DEG2RAD) + cosra = math.cos(ra * 15 * _DEG2RAD) + + uze = [coslat*coslon, coslat*sinlon, sinlat] + une = [-sinlat*coslon, -sinlat*sinlon, coslat] + uwe = [sinlon, -coslon, 0.0] + + uz = _ter2cel(time, uze) + un = _ter2cel(time, une) + uw = _ter2cel(time, uwe) + + p = [cosdc*cosra, cosdc*sinra, sindc] + + pz = p[0]*uz[0] + p[1]*uz[1] + p[2]*uz[2] + pn = p[0]*un[0] + p[1]*un[1] + p[2]*un[2] + pw = p[0]*uw[0] + p[1]*uw[1] + p[2]*uw[2] + + proj = math.sqrt(pn*pn + pw*pw) + az = 0.0 + if proj > 0.0: + az = -math.atan2(pw, pn) * _RAD2DEG + if az < 0: + az += 360 + if az >= 360: + az -= 360 + zd = math.atan2(proj, pz) * _RAD2DEG + hor_ra = ra + hor_dec = dec + + if refraction != REFRACTION_NONE: + zd0 = zd + + # http://extras.springer.com/1999/978-1-4471-0555-8/chap4/horizons/horizons.pdf + # JPL Horizons says it uses refraction algorithm from + # Meeus "Astronomical Algorithms", 1991, p. 101-102. + # I found the following Go implementation: + # https://github.com/soniakeys/meeus/blob/master/v3/refraction/refract.go + # This is a translation from the function "Saemundsson" there. + # I found experimentally that JPL Horizons clamps the angle to 1 degree below the horizon. + # This is important because the 'refr' formula below goes crazy near hd = -5.11. + hd = 90.0 - zd + if hd < -1.0: + hd = -1.0 + + refr = (1.02 / math.tan((hd+10.3/(hd+5.11))*_DEG2RAD)) / 60.0 + + if refraction == REFRACTION_NORMAL and zd > 91.0: + # In "normal" mode we gradually reduce refraction toward the nadir + # so that we never get an altitude angle less than -90 degrees. + # When horizon angle is -1 degrees, zd = 91, and the factor is exactly 1. + # As zd approaches 180 (the nadir), the fraction approaches 0 linearly. + refr *= (180.0 - zd) / 89.0 + + zd -= refr + + if refr > 0.0 and zd > 3.0e-4: + sinzd = math.sin(zd * _DEG2RAD) + coszd = math.cos(zd * _DEG2RAD) + sinzd0 = math.sin(zd0 * _DEG2RAD) + coszd0 = math.cos(zd0 * _DEG2RAD) + + pr = [(((p[j] - coszd0 * uz[j]) / sinzd0)*sinzd + uz[j]*coszd) for j in range(3)] + proj = math.sqrt(pr[0]*pr[0] + pr[1]*pr[1]) + if proj > 0: + hor_ra = math.atan2(pr[1], pr[0]) * _RAD2DEG / 15 + if hor_ra < 0: + hor_ra += 24 + if hor_ra >= 24: + hor_ra -= 24 + else: + hor_ra = 0 + hor_dec = math.atan2(pr[2], proj) * _RAD2DEG + + return HorizontalCoordinates(az, 90.0 - zd, hor_ra, hor_dec) + \ No newline at end of file diff --git a/generate/test.py b/generate/test.py index 5ea8be44..9124cbad 100644 --- a/generate/test.py +++ b/generate/test.py @@ -65,9 +65,16 @@ def Test_AstroCheck(): pos = astronomy.HelioVector(body, time) print('v {} {:0.16f} {:0.16f} {:0.16f} {:0.16f}'.format(name, pos.t.tt, pos.x, pos.y, pos.z)) if body != astronomy.BODY_EARTH: - pass + j2000 = astronomy.Equator(body, time, observer, False, False) + ofdate = astronomy.Equator(body, time, observer, True, True) + hor = astronomy.Horizon(time, observer, ofdate.ra, ofdate.dec, astronomy.REFRACTION_NONE) + print('s {} {:0.16f} {:0.16f} {:0.16f} {:0.16f} {:0.16f} {:0.16f} {:0.16f}'.format(name, time.tt, time.ut, j2000.ra, j2000.dec, j2000.dist, hor.azimuth, hor.altitude)) pos = astronomy.GeoMoon(time) print('v GM {:0.16f} {:0.16f} {:0.16f} {:0.16f}'.format(pos.t.tt, pos.x, pos.y, pos.z)) + j2000 = astronomy.Equator(astronomy.BODY_MOON, time, observer, False, False) + ofdate = astronomy.Equator(astronomy.BODY_MOON, time, observer, True, True) + hor = astronomy.Horizon(time, observer, ofdate.ra, ofdate.dec, astronomy.REFRACTION_NONE) + print('s GM {:0.16f} {:0.16f} {:0.16f} {:0.16f} {:0.16f} {:0.16f} {:0.16f}'.format(time.tt, time.ut, j2000.ra, j2000.dec, j2000.dist, hor.azimuth, hor.altitude)) time = time.AddDays(dt) if len(sys.argv) == 2: diff --git a/generate/unit_test_python b/generate/unit_test_python index 3c823e3d..4a9e71ef 100755 --- a/generate/unit_test_python +++ b/generate/unit_test_python @@ -8,5 +8,11 @@ Fail() python3 --version || Fail "Cannot print python version" python3 test.py time || Fail "Failure reported by test.py (time)" python3 test.py moon || Fail "Failure reported by test.py (moon)" -python3 test.py astro_check > temp/py_check.txt || Fail "Failure in Python astro_check" +echo "$0: Generating Python test output." +time python3 test.py astro_check > temp/py_check.txt || Fail "Failure in Python astro_check" +./generate check temp/py_check.txt || Fail "Verification failure for Python unit test output." +if false; then + # For some reason (alt, az) calculations are a tiny bit off (8e-11) in Python. + ./ctest diff temp/{py,js}_check.txt || Fail "Diff(py,js) failure." +fi exit 0 diff --git a/source/python/astronomy.py b/source/python/astronomy.py index 1cf7625d..54810607 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -454,11 +454,11 @@ _cls_t = [ class _iau2000b: def __init__(self, time): t = time.tt / 36525 - el = ((485868.249036 + t * 1717915923.2178) % ASEC360) * ASEC2RAD - elp = ((1287104.79305 + t * 129596581.0481) % ASEC360) * ASEC2RAD - f = ((335779.526232 + t * 1739527262.8478) % ASEC360) * ASEC2RAD - d = ((1072260.70369 + t * 1602961601.2090) % ASEC360) * ASEC2RAD - om = ((450160.398036 - t * 6962890.5431) % ASEC360) * ASEC2RAD + el = ((485868.249036 + t * 1717915923.2178) % _ASEC360) * _ASEC2RAD + elp = ((1287104.79305 + t * 129596581.0481) % _ASEC360) * _ASEC2RAD + f = ((335779.526232 + t * 1739527262.8478) % _ASEC360) * _ASEC2RAD + d = ((1072260.70369 + t * 1602961601.2090) % _ASEC360) * _ASEC2RAD + om = ((450160.398036 - t * 6962890.5431) % _ASEC360) * _ASEC2RAD dp = 0 de = 0 i = 76 @@ -1981,3 +1981,114 @@ def GeoVector(body, time, aberration): raise Error('Light-travel time solver did not converge: dt={}'.format(dt)) + +def Equator(body, time, observer, ofdate, aberration): + gc_observer = _geo_pos(time, observer) + gc = GeoVector(body, time, aberration) + j2000 = [ + gc.x - gc_observer[0], + gc.y - gc_observer[1], + gc.z - gc_observer[2] + ] + if not ofdate: + return _vector2radec(j2000) + temp = _precession(0, j2000, time.tt) + datevect = _nutation(time, 0, temp) + return _vector2radec(datevect) + +REFRACTION_NONE = 0 +REFRACTION_NORMAL = 1 +REFRACTION_JPLHOR = 2 + +class HorizontalCoordinates: + def __init__(self, azimuth, altitude, ra, dec): + self.azimuth = azimuth + self.altitude = altitude + self.ra = ra + self.dec = dec + +def Horizon(time, observer, ra, dec, refraction): + if not (REFRACTION_NONE <= refraction <= REFRACTION_JPLHOR): + raise Error('Invalid refraction type: ' + str(refraction)) + + sinlat = math.sin(observer.latitude * _DEG2RAD) + coslat = math.cos(observer.latitude * _DEG2RAD) + sinlon = math.sin(observer.longitude * _DEG2RAD) + coslon = math.cos(observer.longitude * _DEG2RAD) + sindc = math.sin(dec * _DEG2RAD) + cosdc = math.cos(dec * _DEG2RAD) + sinra = math.sin(ra * 15 * _DEG2RAD) + cosra = math.cos(ra * 15 * _DEG2RAD) + + uze = [coslat*coslon, coslat*sinlon, sinlat] + une = [-sinlat*coslon, -sinlat*sinlon, coslat] + uwe = [sinlon, -coslon, 0.0] + + uz = _ter2cel(time, uze) + un = _ter2cel(time, une) + uw = _ter2cel(time, uwe) + + p = [cosdc*cosra, cosdc*sinra, sindc] + + pz = p[0]*uz[0] + p[1]*uz[1] + p[2]*uz[2] + pn = p[0]*un[0] + p[1]*un[1] + p[2]*un[2] + pw = p[0]*uw[0] + p[1]*uw[1] + p[2]*uw[2] + + proj = math.sqrt(pn*pn + pw*pw) + az = 0.0 + if proj > 0.0: + az = -math.atan2(pw, pn) * _RAD2DEG + if az < 0: + az += 360 + if az >= 360: + az -= 360 + zd = math.atan2(proj, pz) * _RAD2DEG + hor_ra = ra + hor_dec = dec + + if refraction != REFRACTION_NONE: + zd0 = zd + + # http://extras.springer.com/1999/978-1-4471-0555-8/chap4/horizons/horizons.pdf + # JPL Horizons says it uses refraction algorithm from + # Meeus "Astronomical Algorithms", 1991, p. 101-102. + # I found the following Go implementation: + # https://github.com/soniakeys/meeus/blob/master/v3/refraction/refract.go + # This is a translation from the function "Saemundsson" there. + # I found experimentally that JPL Horizons clamps the angle to 1 degree below the horizon. + # This is important because the 'refr' formula below goes crazy near hd = -5.11. + hd = 90.0 - zd + if hd < -1.0: + hd = -1.0 + + refr = (1.02 / math.tan((hd+10.3/(hd+5.11))*_DEG2RAD)) / 60.0 + + if refraction == REFRACTION_NORMAL and zd > 91.0: + # In "normal" mode we gradually reduce refraction toward the nadir + # so that we never get an altitude angle less than -90 degrees. + # When horizon angle is -1 degrees, zd = 91, and the factor is exactly 1. + # As zd approaches 180 (the nadir), the fraction approaches 0 linearly. + refr *= (180.0 - zd) / 89.0 + + zd -= refr + + if refr > 0.0 and zd > 3.0e-4: + sinzd = math.sin(zd * _DEG2RAD) + coszd = math.cos(zd * _DEG2RAD) + sinzd0 = math.sin(zd0 * _DEG2RAD) + coszd0 = math.cos(zd0 * _DEG2RAD) + + pr = [(((p[j] - coszd0 * uz[j]) / sinzd0)*sinzd + uz[j]*coszd) for j in range(3)] + proj = math.sqrt(pr[0]*pr[0] + pr[1]*pr[1]) + if proj > 0: + hor_ra = math.atan2(pr[1], pr[0]) * _RAD2DEG / 15 + if hor_ra < 0: + hor_ra += 24 + if hor_ra >= 24: + hor_ra -= 24 + else: + hor_ra = 0 + hor_dec = math.atan2(pr[2], proj) * _RAD2DEG + + return HorizontalCoordinates(az, 90.0 - zd, hor_ra, hor_dec) + \ No newline at end of file