From 86674e3c73e85a38595af8301da9d9a586f4da32 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Sun, 30 Jun 2019 18:58:04 -0400 Subject: [PATCH] Python: implemented rise/set tests and fixed associated bugs. --- generate/ctest.c | 2 +- generate/template/astronomy.c | 2 +- generate/template/astronomy.py | 36 ++++++------- generate/test.py | 93 ++++++++++++++++++++++++++++++++++ generate/unit_test_python | 1 + source/c/astronomy.c | 2 +- source/python/astronomy.py | 36 ++++++------- 7 files changed, 127 insertions(+), 45 deletions(-) diff --git a/generate/ctest.c b/generate/ctest.c index 74071c27..59d313a9 100644 --- a/generate/ctest.c +++ b/generate/ctest.c @@ -1213,7 +1213,7 @@ static int RiseSet(const char *filename) if (error_minutes > max_minutes) max_minutes = error_minutes; - if (error_minutes > 2.0) + if (error_minutes > 0.56) { fprintf(stderr, "RiseSet(%s line %d): excessive prediction time error = %lg minutes.\n", filename, lnum, error_minutes); error = 1; diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index 70908dc6..26a20d26 100644 --- a/generate/template/astronomy.c +++ b/generate/template/astronomy.c @@ -3578,7 +3578,7 @@ astro_search_result_t Astronomy_SearchRiseSet( return result; } - /* If we didn't find the desired event, use time_after to find the next before-event. */ + /* If we didn't find the desired event, use evt_after.time to find the next before-event. */ evt_before = Astronomy_SearchHourAngle(body, observer, ha_before, evt_after.time); if (evt_before.status != ASTRO_SUCCESS) return SearchError(evt_before.status); diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index d8ec8ca8..0ec1aef7 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -1807,7 +1807,8 @@ def _peak_altitude(context, time): # We calculate altitude without refraction, then add fixed refraction near the horizon. # This gives us the time of rise/set without the extra work. hor = Horizon(time, context.observer, ofdate.ra, ofdate.dec, REFRACTION_NONE) - return context.direction*(hor.altitude + _RAD2DEG*(context.body_radius_au / ofdate.dist)) + _REFRACTION_NEAR_HORIZON + alt = hor.altitude + _RAD2DEG*(context.body_radius_au / ofdate.dist) + return context.direction * (alt + _REFRACTION_NEAR_HORIZON) def SearchRiseSet(body, observer, direction, startTime, limitDays): if body == BODY_EARTH: @@ -1841,38 +1842,31 @@ def SearchRiseSet(body, observer, direction, startTime, limitDays): alt_before = _peak_altitude(context, time_start) if alt_before > 0.0: # We are past the sought event, so we have to wait for the next "before" event (culm/bottom). - time_before = SearchHourAngle(body, observer, ha_before, time_start) - if time_before is None: - return None + evt_before = SearchHourAngle(body, observer, ha_before, time_start) + time_before = evt_before.time alt_before = _peak_altitude(context, time_before) else: # We are before or at the sought ebvent, so we find the next "after" event (bottom/culm), # and use the current time as the "before" event. time_before = time_start - time_after = SearchHourAngle(body, observer, ha_after, time_before) - if time_after is None: - return None - alt_after = _peak_altitude(context, time_after) + evt_after = SearchHourAngle(body, observer, ha_after, time_before) + alt_after = _peak_altitude(context, evt_after.time) while True: if alt_before <= 0.0 and alt_after > 0.0: - # Search between time_before and time_after for the desired event. - event_time = Search(_peak_altitude, context, time_before, time_after, 1.0) + # Search between the "before time" and the "after time" for the desired event. + event_time = Search(_peak_altitude, context, time_before, evt_after.time, 1.0) if event_time is not None: return event_time - - # We didn't find the desired event, so use time_after to find the next "before" event. - time_before = SearchHourAngle(body, observer, ha_before, time_after) - if time_before is None: + # We didn't find the desired event, so use the "after" time to find the next "before" event. + evt_before = SearchHourAngle(body, observer, ha_before, evt_after.time) + evt_after = SearchHourAngle(body, observer, ha_after, evt_before.time) + if evt_before.time.ut >= time_start.ut + limitDays: return None - - time_after = SearchHourAngle(body, observer, ha_after, time_before) - if time_after is None: - return None - - alt_before = _peak_altitude(context, time_before) - alt_after = _peak_altitude(context, time_after) + time_before = evt_before.time + alt_before = _peak_altitude(context, evt_before.time) + alt_after = _peak_altitude(context, evt_after.time) class SeasonInfo: def __init__(self, mar_equinox, jun_solstice, sep_equinox, dec_solstice): diff --git a/generate/test.py b/generate/test.py index 6de31074..0f461aca 100644 --- a/generate/test.py +++ b/generate/test.py @@ -436,6 +436,97 @@ def Test_Elongation(): #----------------------------------------------------------------------------------------------------------- +def Test_RiseSet(filename): + sum_minutes = 0.0 + max_minutes = 0.0 + nudge_days = 0.01 + observer = None + current_body = None + a_dir = 0 + b_dir = 0 + with open(filename, 'rt') as infile: + lnum = 0 + for line in infile: + lnum += 1 + line = line.strip() + # Moon 103 -61 1944-01-02T17:08Z s + # Moon 103 -61 1944-01-03T05:47Z r + m = re.match(r'^([A-Za-z]+)\s+(-?[0-9\.]+)\s+(-?[0-9\.]+)\s+(\d+)-(\d+)-(\d+)T(\d+):(\d+)Z\s+([sr])$', line) + if not m: + print('Test_RiseSet({} line {}): invalid data format'.format(filename, lnum)) + return 1 + name = m.group(1) + longitude = float(m.group(2)) + latitude = float(m.group(3)) + year = int(m.group(4)) + month = int(m.group(5)) + day = int(m.group(6)) + hour = int(m.group(7)) + minute = int(m.group(8)) + kind = m.group(9) + correct_time = astronomy.Time.Make(year, month, day, hour, minute, 0) + direction = astronomy.DIRECTION_RISE if kind == 'r' else astronomy.DIRECTION_SET + body = astronomy.BodyCode(name) + if body < 0: + print('Test_RiseSet({} line {}): invalid body name "{}"'.format(filename, lnum, name)) + return 1 + + # Every time we see a new geographic location, start a new iteration + # of finding all rise/set times for that UTC calendar year. + if (observer is None) or (observer.latitude != latitude) or (observer.longitude != longitude) or (current_body != body): + current_body = body + observer = astronomy.Observer(latitude, longitude, 0) + r_search_date = s_search_date = astronomy.Time.Make(year, 1, 1, 0, 0, 0) + b_evt = None + print('Test_RiseSet: {:<7s} lat={:0.1f} lon={:0.1f}'.format(name, latitude, longitude)) + + if b_evt is not None: + # Recycle the second event from the previous iteration as the first event. + a_evt = b_evt + a_dir = b_dir + b_evt = None + else: + r_evt = astronomy.SearchRiseSet(body, observer, astronomy.DIRECTION_RISE, r_search_date, 366.0) + if r_evt is None: + print('Test_RiseSet({} line {}): rise search failed'.format(filename, lnum)) + return 1 + s_evt = astronomy.SearchRiseSet(body, observer, astronomy.DIRECTION_SET, s_search_date, 366.0) + if s_evt is None: + print('Test_RiseSet({} line {}): set search failed'.format(filename, lnum)) + return 1 + # Expect the current event to match the earlier of the found times. + if r_evt.tt < s_evt.tt: + a_evt = r_evt + b_evt = s_evt + a_dir = astronomy.DIRECTION_RISE + b_dir = astronomy.DIRECTION_SET + else: + a_evt = s_evt + b_evt = r_evt + a_dir = astronomy.DIRECTION_SET + b_dir = astronomy.DIRECTION_RISE + # Nudge the event times forward a tiny amount. + r_search_date = r_evt.AddDays(nudge_days) + s_search_date = s_evt.AddDays(nudge_days) + + if a_dir != direction: + print('Test_RiseSet({} line {}): expected dir={} but found {}'.format(filename, lnum, a_dir, direction)) + return 1 + + error_minutes = (24.0 * 60.0) * abs(a_evt.tt - correct_time.tt) + sum_minutes += error_minutes ** 2 + max_minutes = max(max_minutes, error_minutes) + if error_minutes > 0.56: + print('Test_RiseSet({} line {}): excessive prediction time error = {} minutes.'.format(filename, lnum, error_minutes)) + print(' correct = {}, calculated = {}'.format(correct_time, a_evt)) + return 1 + + rms_minutes = math.sqrt(sum_minutes / lnum) + print('Test_RiseSet: passed {} lines: time errors in minutes: rms={:0.4f}, max={:0.4f}'.format(lnum, rms_minutes, max_minutes)) + return 0 + +#----------------------------------------------------------------------------------------------------------- + if len(sys.argv) == 2: if sys.argv[1] == 'time': Test_AstroTime() @@ -454,6 +545,8 @@ if len(sys.argv) == 3: sys.exit(Test_Seasons(sys.argv[2])) if sys.argv[1] == 'moonphase': sys.exit(Test_MoonPhase(sys.argv[2])) + if sys.argv[1] == 'riseset': + sys.exit(Test_RiseSet(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 ca65347e..312678ff 100755 --- a/generate/unit_test_python +++ b/generate/unit_test_python @@ -10,6 +10,7 @@ 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." python3 test.py moonphase moonphase/moonphases.txt || Fail "Failed Python moon phase test." +python3 test.py riseset riseset/riseset.txt || Fail "Failed Python rise/set tests." python3 test.py elongation || Fail "Failed Python elongation tests." for file in temp/py_longitude_*.txt; do diff --git a/source/c/astronomy.c b/source/c/astronomy.c index f435ca31..f87fd190 100644 --- a/source/c/astronomy.c +++ b/source/c/astronomy.c @@ -4634,7 +4634,7 @@ astro_search_result_t Astronomy_SearchRiseSet( return result; } - /* If we didn't find the desired event, use time_after to find the next before-event. */ + /* If we didn't find the desired event, use evt_after.time to find the next before-event. */ evt_before = Astronomy_SearchHourAngle(body, observer, ha_before, evt_after.time); if (evt_before.status != ASTRO_SUCCESS) return SearchError(evt_before.status); diff --git a/source/python/astronomy.py b/source/python/astronomy.py index 7536716e..24609a6f 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -2645,7 +2645,8 @@ def _peak_altitude(context, time): # We calculate altitude without refraction, then add fixed refraction near the horizon. # This gives us the time of rise/set without the extra work. hor = Horizon(time, context.observer, ofdate.ra, ofdate.dec, REFRACTION_NONE) - return context.direction*(hor.altitude + _RAD2DEG*(context.body_radius_au / ofdate.dist)) + _REFRACTION_NEAR_HORIZON + alt = hor.altitude + _RAD2DEG*(context.body_radius_au / ofdate.dist) + return context.direction * (alt + _REFRACTION_NEAR_HORIZON) def SearchRiseSet(body, observer, direction, startTime, limitDays): if body == BODY_EARTH: @@ -2679,38 +2680,31 @@ def SearchRiseSet(body, observer, direction, startTime, limitDays): alt_before = _peak_altitude(context, time_start) if alt_before > 0.0: # We are past the sought event, so we have to wait for the next "before" event (culm/bottom). - time_before = SearchHourAngle(body, observer, ha_before, time_start) - if time_before is None: - return None + evt_before = SearchHourAngle(body, observer, ha_before, time_start) + time_before = evt_before.time alt_before = _peak_altitude(context, time_before) else: # We are before or at the sought ebvent, so we find the next "after" event (bottom/culm), # and use the current time as the "before" event. time_before = time_start - time_after = SearchHourAngle(body, observer, ha_after, time_before) - if time_after is None: - return None - alt_after = _peak_altitude(context, time_after) + evt_after = SearchHourAngle(body, observer, ha_after, time_before) + alt_after = _peak_altitude(context, evt_after.time) while True: if alt_before <= 0.0 and alt_after > 0.0: - # Search between time_before and time_after for the desired event. - event_time = Search(_peak_altitude, context, time_before, time_after, 1.0) + # Search between the "before time" and the "after time" for the desired event. + event_time = Search(_peak_altitude, context, time_before, evt_after.time, 1.0) if event_time is not None: return event_time - - # We didn't find the desired event, so use time_after to find the next "before" event. - time_before = SearchHourAngle(body, observer, ha_before, time_after) - if time_before is None: + # We didn't find the desired event, so use the "after" time to find the next "before" event. + evt_before = SearchHourAngle(body, observer, ha_before, evt_after.time) + evt_after = SearchHourAngle(body, observer, ha_after, evt_before.time) + if evt_before.time.ut >= time_start.ut + limitDays: return None - - time_after = SearchHourAngle(body, observer, ha_after, time_before) - if time_after is None: - return None - - alt_before = _peak_altitude(context, time_before) - alt_after = _peak_altitude(context, time_after) + time_before = evt_before.time + alt_before = _peak_altitude(context, evt_before.time) + alt_after = _peak_altitude(context, evt_after.time) class SeasonInfo: def __init__(self, mar_equinox, jun_solstice, sep_equinox, dec_solstice):