From 9d304850df23b0a1ccdd0ec1d568914be7af867f Mon Sep 17 00:00:00 2001 From: Don Cross Date: Sat, 29 Jun 2019 23:02:27 -0400 Subject: [PATCH] Python: Added seasons test and fixed errors uncovered by it. --- generate/template/astronomy.py | 28 +++++++------- generate/test.py | 70 ++++++++++++++++++++++++++++++++++ generate/unit_test_python | 2 + source/python/astronomy.py | 28 +++++++------- 4 files changed, 98 insertions(+), 30 deletions(-) diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index 9a3b02ab..67be9ea1 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -1285,19 +1285,19 @@ def _RotateEquatorialToEcliptic(pos, obliq_radians): ez = -pos[1]*sin_ob + pos[2]*cos_ob xyproj = math.sqrt(ex*ex + ey*ey) if xyproj > 0.0: - elon = _RAD2DEG * math.atan(ey, ex) + elon = _RAD2DEG * math.atan2(ey, ex) if elon < 0.0: elon += 360.0 else: elon = 0.0 - elat = _RAD2DEG * math.atan(ez, xyproj) + elat = _RAD2DEG * math.atan2(ez, xyproj) return EclipticCoordinates(ex, ey, ez, elat, elon) def SunPosition(time): # Correct for light travel time from the Sun. # Otherwise season calculations (equinox, solstice) will all be early by about 8 minutes! adjusted_time = time.AddDays(-1.0 / _C_AUDAY) - earth2000 = CalcEarth(adjusted_time) + earth2000 = _CalcEarth(adjusted_time) sun2000 = [-earth2000.x, -earth2000.y, -earth2000.z] # Convert to equatorial Cartesian coordinates of date. @@ -1306,12 +1306,12 @@ def SunPosition(time): # Convert equatorial coordinates to ecliptic coordinates. true_obliq = _DEG2RAD * adjusted_time._etilt().tobl - return RotateEquatorialToEcliptic(sun_ofdate, true_obliq) + return _RotateEquatorialToEcliptic(sun_ofdate, true_obliq) def Ecliptic(equ): # Based on NOVAS functions equ2ecl() and equ2ecl_vec(). ob2000 = 0.40909260059599012 # mean obliquity of the J2000 ecliptic in radians - return RotateEquatorialToEcliptic([equ.x, equ.y, equ.z], ob2000) + return _RotateEquatorialToEcliptic([equ.x, equ.y, equ.z], ob2000) def EclipticLongitude(body, time): if body == BODY_SUN: @@ -1491,7 +1491,7 @@ def _sun_offset(targetLon, time): return _LongitudeOffset(ecl.elon - targetLon) def SearchSunLongitude(targetLon, startTime, limitDays): - t2 = startTime.Add(limitDays) + t2 = startTime.AddDays(limitDays) return Search(_sun_offset, targetLon, startTime, t2, 1.0) def MoonPhase(time): @@ -1877,7 +1877,7 @@ def SearchRiseSet(body, observer, direction, startTime, limitDays): alt_after = _peak_altitude(context, time_after) class SeasonInfo: - def __init__(mar_equinox, jun_solstice, sep_equinox, dec_solstice): + def __init__(self, mar_equinox, jun_solstice, sep_equinox, dec_solstice): self.mar_equinox = mar_equinox self.jun_solstice = jun_solstice self.sep_equinox = sep_equinox @@ -1885,14 +1885,18 @@ class SeasonInfo: def _FindSeasonChange(targetLon, year, month, day): startTime = Time.Make(year, month, day, 0, 0, 0) - return SearchSunLongitude(targetLon, startTime, 4.0) + time = SearchSunLongitude(targetLon, startTime, 4.0) + if time is None: + # We should always be able to find a season change. + raise InternalError() + return time def Seasons(year): mar_equinox = _FindSeasonChange(0, year, 3, 19) jun_solstice = _FindSeasonChange(90, year, 6, 19) sep_equinox = _FindSeasonChange(180, year, 9, 21) dec_solstice = _FindSeasonChange(270, year, 12, 20) - return SeasonInfo(mar_equinox, jun_solsitce, sep_equinox, dec_solstice) + return SeasonInfo(mar_equinox, jun_solstice, sep_equinox, dec_solstice) def _MoonDistance(time): return _CalcMoon(time).distance_au @@ -1975,10 +1979,6 @@ def NextLunarApsis(apsis): return next #================================================================================================== -# + SearchSunLongitude -# + _sun_offset -# + SunPosition -# + RotateEquatorialToEcliptic # + Ecliptic # + EclipticLongitude # + AngleFromSun @@ -1994,8 +1994,6 @@ def NextLunarApsis(apsis): # + NextMoonQuarter # + SearchHourAngle # + SearchRiseSet -# + Seasons -# + _FindSeasonChange # + Illumination # + SearchPeakMagnitude # + SearchLunarApsis diff --git a/generate/test.py b/generate/test.py index 9219f198..fd213937 100644 --- a/generate/test.py +++ b/generate/test.py @@ -1,6 +1,7 @@ #!/usr/bin/env python3 import sys import math +import re sys.path.append('../source/python') import astronomy @@ -82,6 +83,71 @@ def Test_AstroCheck(printflag): 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) +def Test_Seasons(filename): + with open(filename, 'rt') as infile: + lnum = 0 + current_year = 0 + mar_count = sep_count = jun_count = dec_count = 0 + max_minutes = 0.0 + for line in infile: + lnum += 1 + line = line.strip() + m = re.match(r'^(\d+)-(\d+)-(\d+)T(\d+):(\d+)Z\s+([A-Za-z]+)$', line) + if not m: + print('Test_Seasons: Invalid data on line {} of file {}'.format(lnum, filename)) + return 1 + year = int(m.group(1)) + month = int(m.group(2)) + day = int(m.group(3)) + hour = int(m.group(4)) + minute = int(m.group(5)) + name = m.group(6) + if year != current_year: + current_year = current_year + seasons = astronomy.Seasons(year) + correct_time = astronomy.Time.Make(year, month, day, hour, minute, 0) + if name == 'Equinox': + if month == 3: + calc_time = seasons.mar_equinox + mar_count += 1 + elif month == 9: + calc_time = seasons.sep_equinox + sep_count += 1 + else: + print('Test_Seasons: {} line {}: Invalid equinox date in test data'.format(filename, lnum)) + return 1 + elif name == 'Solstice': + if month == 6: + calc_time = seasons.jun_solstice + jun_count += 1 + elif month == 12: + calc_time = seasons.dec_solstice + dec_count += 1 + else: + print('Test_Seasons: {} line {}: Invalid solstice date in test data'.format(filename, lnum)) + return 1 + elif name == 'Aphelion': + continue # not yet calculated + elif name == 'Perihelion': + continue # not yet calculated + else: + print('Test_Seasons: {} line {}: unknown event type {}'.format(filename, lnum, name)) + return 1 + + # Verify that the calculated time matches the correct time for this event. + diff_minutes = (24.0 * 60.0) * abs(calc_time.tt - correct_time.tt) + if diff_minutes > max_minutes: + max_minutes = diff_minutes + + if diff_minutes > 1.7: + print('Test_Seasons: {} line {}: excessive error ({}): {} minutes.'.format(filename, lnum, name, diff_minutes)) + return 1 + print('Test_Seasons: verified {} lines from file {} : max error minutes = {:0.3f}'.format(lnum, filename, max_minutes)) + print('Test_Seasons: Event counts: mar={}, jun={}, sep={}, dec={}'.format(mar_count, jun_count, sep_count, dec_count)) + return 0 + +#----------------------------------------------------------------------------------------------------------- + if len(sys.argv) == 2: if sys.argv[1] == 'time': Test_AstroTime() @@ -95,5 +161,9 @@ if len(sys.argv) == 2: Test_AstroCheck(sys.argv[1] == 'astro_check') sys.exit(0) +if len(sys.argv) == 3: + if sys.argv[1] == 'seasons': + sys.exit(Test_Seasons(sys.argv[2])) + print('test.py: Invalid command line arguments.') sys.exit(1) diff --git a/generate/unit_test_python b/generate/unit_test_python index 4a9e71ef..37ac576d 100755 --- a/generate/unit_test_python +++ b/generate/unit_test_python @@ -8,6 +8,8 @@ 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 seasons seasons/seasons.txt || Fail "Failed Python seasons test." + 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." diff --git a/source/python/astronomy.py b/source/python/astronomy.py index c05d0709..f1e7d41d 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -2123,19 +2123,19 @@ def _RotateEquatorialToEcliptic(pos, obliq_radians): ez = -pos[1]*sin_ob + pos[2]*cos_ob xyproj = math.sqrt(ex*ex + ey*ey) if xyproj > 0.0: - elon = _RAD2DEG * math.atan(ey, ex) + elon = _RAD2DEG * math.atan2(ey, ex) if elon < 0.0: elon += 360.0 else: elon = 0.0 - elat = _RAD2DEG * math.atan(ez, xyproj) + elat = _RAD2DEG * math.atan2(ez, xyproj) return EclipticCoordinates(ex, ey, ez, elat, elon) def SunPosition(time): # Correct for light travel time from the Sun. # Otherwise season calculations (equinox, solstice) will all be early by about 8 minutes! adjusted_time = time.AddDays(-1.0 / _C_AUDAY) - earth2000 = CalcEarth(adjusted_time) + earth2000 = _CalcEarth(adjusted_time) sun2000 = [-earth2000.x, -earth2000.y, -earth2000.z] # Convert to equatorial Cartesian coordinates of date. @@ -2144,12 +2144,12 @@ def SunPosition(time): # Convert equatorial coordinates to ecliptic coordinates. true_obliq = _DEG2RAD * adjusted_time._etilt().tobl - return RotateEquatorialToEcliptic(sun_ofdate, true_obliq) + return _RotateEquatorialToEcliptic(sun_ofdate, true_obliq) def Ecliptic(equ): # Based on NOVAS functions equ2ecl() and equ2ecl_vec(). ob2000 = 0.40909260059599012 # mean obliquity of the J2000 ecliptic in radians - return RotateEquatorialToEcliptic([equ.x, equ.y, equ.z], ob2000) + return _RotateEquatorialToEcliptic([equ.x, equ.y, equ.z], ob2000) def EclipticLongitude(body, time): if body == BODY_SUN: @@ -2329,7 +2329,7 @@ def _sun_offset(targetLon, time): return _LongitudeOffset(ecl.elon - targetLon) def SearchSunLongitude(targetLon, startTime, limitDays): - t2 = startTime.Add(limitDays) + t2 = startTime.AddDays(limitDays) return Search(_sun_offset, targetLon, startTime, t2, 1.0) def MoonPhase(time): @@ -2715,7 +2715,7 @@ def SearchRiseSet(body, observer, direction, startTime, limitDays): alt_after = _peak_altitude(context, time_after) class SeasonInfo: - def __init__(mar_equinox, jun_solstice, sep_equinox, dec_solstice): + def __init__(self, mar_equinox, jun_solstice, sep_equinox, dec_solstice): self.mar_equinox = mar_equinox self.jun_solstice = jun_solstice self.sep_equinox = sep_equinox @@ -2723,14 +2723,18 @@ class SeasonInfo: def _FindSeasonChange(targetLon, year, month, day): startTime = Time.Make(year, month, day, 0, 0, 0) - return SearchSunLongitude(targetLon, startTime, 4.0) + time = SearchSunLongitude(targetLon, startTime, 4.0) + if time is None: + # We should always be able to find a season change. + raise InternalError() + return time def Seasons(year): mar_equinox = _FindSeasonChange(0, year, 3, 19) jun_solstice = _FindSeasonChange(90, year, 6, 19) sep_equinox = _FindSeasonChange(180, year, 9, 21) dec_solstice = _FindSeasonChange(270, year, 12, 20) - return SeasonInfo(mar_equinox, jun_solsitce, sep_equinox, dec_solstice) + return SeasonInfo(mar_equinox, jun_solstice, sep_equinox, dec_solstice) def _MoonDistance(time): return _CalcMoon(time).distance_au @@ -2813,10 +2817,6 @@ def NextLunarApsis(apsis): return next #================================================================================================== -# + SearchSunLongitude -# + _sun_offset -# + SunPosition -# + RotateEquatorialToEcliptic # + Ecliptic # + EclipticLongitude # + AngleFromSun @@ -2832,8 +2832,6 @@ def NextLunarApsis(apsis): # + NextMoonQuarter # + SearchHourAngle # + SearchRiseSet -# + Seasons -# + _FindSeasonChange # + Illumination # + SearchPeakMagnitude # + SearchLunarApsis