From d01c4fd0a5cab98c73ae3340f803706c06fab4e0 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Mon, 1 Jul 2019 17:00:28 -0400 Subject: [PATCH] Python: Added tests for astronomy.Illumination(). I still need to add a test for SearchPeakMagnitude(). --- generate/template/astronomy.py | 13 ++--- generate/test.py | 97 ++++++++++++++++++++++++++++++++++ generate/unit_test_python | 1 + source/python/astronomy.py | 13 ++--- 4 files changed, 106 insertions(+), 18 deletions(-) diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index 0ec1aef7..b2090f53 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -1578,8 +1578,8 @@ def _SaturnMagnitude(phase, helio_dist, geo_dist, gc, time): # Find geocentric ecliptic coordinates of Saturn. eclip = Ecliptic(gc) - ir = DEG2RAD * 28.06 # tilt of Saturn's rings to the ecliptic, in radians - Nr = DEG2RAD * (169.51 + (3.82e-5 * time.tt)) # ascending node of Saturn's rings, in radians + ir = _DEG2RAD * 28.06 # tilt of Saturn's rings to the ecliptic, in radians + Nr = _DEG2RAD * (169.51 + (3.82e-5 * time.tt)) # ascending node of Saturn's rings, in radians # Find tilt of Saturn's rings, as seen from Earth. lat = _DEG2RAD * eclip.elat @@ -1627,6 +1627,7 @@ def Illumination(body, time): earth = _CalcEarth(time) if body == BODY_SUN: gc = Vector(-earth.x, -earth.y, -earth.z, time) + hc = Vector(0.0, 0.0, 0.0, time) phase = 0.0 # placeholder value; the Sun does not have a phase angle. else: if body == BODY_MOON: @@ -1957,7 +1958,7 @@ def SearchLunarApsis(startTime): def NextLunarApsis(apsis): skip = 11.0 # number of days to skip to start looking for next apsis event expected = APSIS_INVALID - apsis.time = apsis.time.AddDays(skip) + time = apsis.time.AddDays(skip) next = SearchLunarApsis(time) # Verify that we found the opposite apsis from the previous one. if apsis.kind == APSIS_APOCENTER: @@ -1971,12 +1972,6 @@ def NextLunarApsis(apsis): return next #================================================================================================== -# + Ecliptic -# + EclipticLongitude -# + AngleFromSun -# + SearchHourAngle -# + SearchRiseSet -# + Illumination # + SearchPeakMagnitude # + SearchLunarApsis # + NextLunarApsis diff --git a/generate/test.py b/generate/test.py index 0f461aca..67181d94 100644 --- a/generate/test.py +++ b/generate/test.py @@ -436,6 +436,101 @@ def Test_Elongation(): #----------------------------------------------------------------------------------------------------------- +def ParseJplHorizonsDateTime(line): + m = re.match(r'^\s*(\d{4})-(Jan|Feb|Mar|Apr|May|Jun|Jul|Aug|Sep|Oct|Nov|Dec)-(\d{2})\s(\d{2}):(\d{2})\s+(.*)$', line) + if not m: + return None, None + year = int(m.group(1)) + month = 1 + ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'].index(m.group(2)) + day = int(m.group(3)) + hour = int(m.group(4)) + minute = int(m.group(5)) + rest = m.group(6) + time = astronomy.Time.Make(year, month, day, hour, minute, 0) + return time, rest + +def CheckMagnitudeData(body, filename): + limit = 0.012 + sum_squared_diff = 0.0 + with open(filename, 'rt') as infile: + count = lnum = 0 + for line in infile: + lnum += 1 + line = line.strip() + (time, rest) = ParseJplHorizonsDateTime(line) + if (time is not None) and (rest is not None) and not ('n.a.' in rest): + data = [float(t) for t in rest.split()] + if len(data) != 7: + print('CheckMagnitudeData({} line {}): invalid data format'.format(filename, lnum)) + return 1 + (mag, sbrt, dist, rdot, delta, deldot, phase_angle) = data + illum = astronomy.Illumination(body, time) + diff = illum.mag - mag + if abs(diff) > limit: + print('CheckMagnitudeData({} line {}): EXCESSIVE ERROR: correct mag={}, calc mag={}'.format(filename, lnum, mag, illum.mag)) + return 1 + sum_squared_diff += diff * diff + if count == 0: + diff_lo = diff_hi = diff + else: + diff_lo = min(diff_lo, diff) + diff_hi = max(diff_hi, diff) + count += 1 + + if count == 0: + print('CheckMagnitudeData: Did not find any data in file: {}'.format(filename)) + return 1 + rms = math.sqrt(sum_squared_diff / count) + print('CheckMagnitudeData: {:<21s} {:5d} rows diff_lo={:0.4f} diff_hi={:0.4f} rms={:0.4f}'.format(filename, count, diff_lo, diff_hi, rms)) + return 0 + +def CheckSaturn(): + # JPL Horizons does not include Saturn's rings in its magnitude models. + # I still don't have authoritative test data for Saturn's magnitude. + # For now, I just test for consistency with Paul Schlyter's formulas at: + # http://www.stjarnhimlen.se/comp/ppcomp.html#15 + data = [ + ( "1972-01-01T00:00Z", -0.31904865, +24.50061220 ), + ( "1980-01-01T00:00Z", +0.85213663, -1.85761461 ), + ( "2009-09-04T00:00Z", +1.01626809, +0.08380716 ), + ( "2017-06-15T00:00Z", -0.12318790, -26.60871409 ), + ( "2019-05-01T00:00Z", +0.32954097, -23.53880802 ), + ( "2025-09-25T00:00Z", +0.51286575, +1.52327932 ), + ( "2032-05-15T00:00Z", -0.04652109, +26.95717765 ) + ] + for (dtext, mag, tilt) in data: + time = ParseDate(dtext) + illum = astronomy.Illumination(astronomy.BODY_SATURN, time) + print('Saturn: date={} calc mag={:12.8f} ring_tilt={:12.8f}'.format(dtext, illum.mag, illum.ring_tilt)) + mag_diff = abs(illum.mag - mag) + if mag_diff > 1.0e-8: + print('CheckSaturn: Excessive magnitude error {}'.format(mag_diff)) + return 1 + tilt_diff = abs(illum.ring_tilt - tilt) + if (tilt_diff > 1.0e-8): + print('CheckSaturn: Excessive ring tilt error {}'.format(tilt_diff)) + return 1 + return 0 + +def Test_Magnitude(): + nfailed = 0 + nfailed += CheckMagnitudeData(astronomy.BODY_SUN, 'magnitude/Sun.txt') + nfailed += CheckMagnitudeData(astronomy.BODY_MOON, 'magnitude/Moon.txt') + nfailed += CheckMagnitudeData(astronomy.BODY_MERCURY, 'magnitude/Mercury.txt') + nfailed += CheckMagnitudeData(astronomy.BODY_VENUS, 'magnitude/Venus.txt') + nfailed += CheckMagnitudeData(astronomy.BODY_MARS, 'magnitude/Mars.txt') + nfailed += CheckMagnitudeData(astronomy.BODY_JUPITER, 'magnitude/Jupiter.txt') + nfailed += CheckSaturn() + nfailed += CheckMagnitudeData(astronomy.BODY_URANUS, 'magnitude/Uranus.txt') + nfailed += CheckMagnitudeData(astronomy.BODY_NEPTUNE, 'magnitude/Neptune.txt') + nfailed += CheckMagnitudeData(astronomy.BODY_PLUTO, 'magnitude/Pluto.txt') + if nfailed > 0: + print('Test_Magnitude: failed {} test(s).'.format(nfailed)) + return 1 + return 0 + +#----------------------------------------------------------------------------------------------------------- + def Test_RiseSet(filename): sum_minutes = 0.0 max_minutes = 0.0 @@ -539,6 +634,8 @@ if len(sys.argv) == 2: sys.exit(0) if sys.argv[1] == 'elongation': sys.exit(Test_Elongation()) + if sys.argv[1] == 'magnitude': + sys.exit(Test_Magnitude()) if len(sys.argv) == 3: if sys.argv[1] == 'seasons': diff --git a/generate/unit_test_python b/generate/unit_test_python index 312678ff..d5900713 100755 --- a/generate/unit_test_python +++ b/generate/unit_test_python @@ -7,6 +7,7 @@ Fail() python3 --version || Fail "Cannot print python version" python3 test.py time || Fail "Failure reported by test.py (time)" +python3 test.py magnitude || Fail "Failed Python magnitude tests." python3 test.py moon || Fail "Failure reported by test.py (moon)" python3 test.py seasons seasons/seasons.txt || Fail "Failed Python seasons test." python3 test.py moonphase moonphase/moonphases.txt || Fail "Failed Python moon phase test." diff --git a/source/python/astronomy.py b/source/python/astronomy.py index 24609a6f..684c28af 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -2416,8 +2416,8 @@ def _SaturnMagnitude(phase, helio_dist, geo_dist, gc, time): # Find geocentric ecliptic coordinates of Saturn. eclip = Ecliptic(gc) - ir = DEG2RAD * 28.06 # tilt of Saturn's rings to the ecliptic, in radians - Nr = DEG2RAD * (169.51 + (3.82e-5 * time.tt)) # ascending node of Saturn's rings, in radians + ir = _DEG2RAD * 28.06 # tilt of Saturn's rings to the ecliptic, in radians + Nr = _DEG2RAD * (169.51 + (3.82e-5 * time.tt)) # ascending node of Saturn's rings, in radians # Find tilt of Saturn's rings, as seen from Earth. lat = _DEG2RAD * eclip.elat @@ -2465,6 +2465,7 @@ def Illumination(body, time): earth = _CalcEarth(time) if body == BODY_SUN: gc = Vector(-earth.x, -earth.y, -earth.z, time) + hc = Vector(0.0, 0.0, 0.0, time) phase = 0.0 # placeholder value; the Sun does not have a phase angle. else: if body == BODY_MOON: @@ -2795,7 +2796,7 @@ def SearchLunarApsis(startTime): def NextLunarApsis(apsis): skip = 11.0 # number of days to skip to start looking for next apsis event expected = APSIS_INVALID - apsis.time = apsis.time.AddDays(skip) + time = apsis.time.AddDays(skip) next = SearchLunarApsis(time) # Verify that we found the opposite apsis from the previous one. if apsis.kind == APSIS_APOCENTER: @@ -2809,12 +2810,6 @@ def NextLunarApsis(apsis): return next #================================================================================================== -# + Ecliptic -# + EclipticLongitude -# + AngleFromSun -# + SearchHourAngle -# + SearchRiseSet -# + Illumination # + SearchPeakMagnitude # + SearchLunarApsis # + NextLunarApsis