From 16aa0b33130437f78aa03a7334b860cdf880ceab Mon Sep 17 00:00:00 2001 From: Don Cross Date: Sat, 29 Jun 2019 20:30:52 -0400 Subject: [PATCH] Python: Finished first pass of porting functions. The next step is to write all the unit tests. --- generate/template/astronomy.py | 106 ++++++++++++++++++++++++++++++++- source/python/astronomy.py | 106 ++++++++++++++++++++++++++++++++- 2 files changed, 206 insertions(+), 6 deletions(-) diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index b86e9e79..9a3b02ab 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -1876,6 +1876,103 @@ def SearchRiseSet(body, observer, direction, startTime, limitDays): alt_before = _peak_altitude(context, time_before) alt_after = _peak_altitude(context, time_after) +class SeasonInfo: + def __init__(mar_equinox, jun_solstice, sep_equinox, dec_solstice): + self.mar_equinox = mar_equinox + self.jun_solstice = jun_solstice + self.sep_equinox = sep_equinox + self.dec_solstice = dec_solstice + +def _FindSeasonChange(targetLon, year, month, day): + startTime = Time.Make(year, month, day, 0, 0, 0) + return SearchSunLongitude(targetLon, startTime, 4.0) + +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) + +def _MoonDistance(time): + return _CalcMoon(time).distance_au + +def _distance_slope(direction, time): + dt = 0.001 + t1 = time.AddDays(-dt/2.0) + t2 = time.AddDays(+dt/2.0) + dist1 = _MoonDistance(t1) + dist2 = _MoonDistance(t2) + return direction * (dist2 - dist1) / dt + +APSIS_PERICENTER = 0 +APSIS_APOCENTER = 1 +APSIS_INVALID = 2 + +class Apsis: + def __init__(self, time, kind, dist_au): + self.time = time + self.kind = kind + self.dist_au = dist_au + self.dist_km = dist_au * _KM_PER_AU + +def SearchLunarApsis(startTime): + increment = 5.0 # number of days to skip on each iteration + t1 = startTime + m1 = _distance_slope(+1, t1) + iter = 0 + while iter * increment < 2.0 * _MEAN_SYNODIC_MONTH: + t2 = t1.AddDays(increment) + m2 = _distance_slope(+1, t2) + if m1 * m2 <= 0.0: + # There is a change of slope polarity within the time range [t1, t2]. + # Therefore this time range contains an apsis. + # Figure out whether it is perigee or apogee. + if m1 < 0.0 or m2 > 0.0: + # We found a minimum-distance event: perigee. + # Search the time range for the time when the slope goes from negative to positive. + apsis_time = Search(_distance_slope, +1, t1, t2, 1.0) + kind = APSIS_PERICENTER + elif m1 > 0.0 or m2 < 0.0: + # We found a maximum-distance event: apogee. + # Search the time range for the time when the slope goes from positive to negative. + apsis_time = Search(_distance_slope, -1, t1, t2, 1.0) + kind = APSIS_APOCENTER + else: + # This should never happen. It should not be possible for both slopes to be zero. + raise InternalError() + + if apsis_time is None: + # We should always be able to find a lunar apsis! + raise InternalError() + + dist = _MoonDistance(apsis_time) + return Apsis(apsis_time, kind, dist) + + # We have not yet found a slope polarity change. Keep searching. + t1 = t2 + m1 = m2 + iter += 1 + + # It should not be possible to fail to find an apsis within 2 synodic months. + raise InternalError() + + +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) + next = SearchLunarApsis(time) + # Verify that we found the opposite apsis from the previous one. + if apsis.kind == APSIS_APOCENTER: + expected = APSIS_PERICENTER + elif apsis.kind == APSIS_PERICENTER: + expected = APSIS_APOCENTER + else: + raise Error('Parameter "apsis" contains an invalid "kind" value.') + if next.kind != expected: + raise InternalError() + return next #================================================================================================== # + SearchSunLongitude @@ -1897,8 +1994,11 @@ def SearchRiseSet(body, observer, direction, startTime, limitDays): # + NextMoonQuarter # + SearchHourAngle # + SearchRiseSet -# - Seasons +# + Seasons +# + _FindSeasonChange # + Illumination # + SearchPeakMagnitude -# - SearchLunarApsis -# - NextLunarApsis +# + SearchLunarApsis +# + _distance_slope +# + _MoonDistance +# + NextLunarApsis diff --git a/source/python/astronomy.py b/source/python/astronomy.py index 000e1357..c05d0709 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -2714,6 +2714,103 @@ def SearchRiseSet(body, observer, direction, startTime, limitDays): alt_before = _peak_altitude(context, time_before) alt_after = _peak_altitude(context, time_after) +class SeasonInfo: + def __init__(mar_equinox, jun_solstice, sep_equinox, dec_solstice): + self.mar_equinox = mar_equinox + self.jun_solstice = jun_solstice + self.sep_equinox = sep_equinox + self.dec_solstice = dec_solstice + +def _FindSeasonChange(targetLon, year, month, day): + startTime = Time.Make(year, month, day, 0, 0, 0) + return SearchSunLongitude(targetLon, startTime, 4.0) + +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) + +def _MoonDistance(time): + return _CalcMoon(time).distance_au + +def _distance_slope(direction, time): + dt = 0.001 + t1 = time.AddDays(-dt/2.0) + t2 = time.AddDays(+dt/2.0) + dist1 = _MoonDistance(t1) + dist2 = _MoonDistance(t2) + return direction * (dist2 - dist1) / dt + +APSIS_PERICENTER = 0 +APSIS_APOCENTER = 1 +APSIS_INVALID = 2 + +class Apsis: + def __init__(self, time, kind, dist_au): + self.time = time + self.kind = kind + self.dist_au = dist_au + self.dist_km = dist_au * _KM_PER_AU + +def SearchLunarApsis(startTime): + increment = 5.0 # number of days to skip on each iteration + t1 = startTime + m1 = _distance_slope(+1, t1) + iter = 0 + while iter * increment < 2.0 * _MEAN_SYNODIC_MONTH: + t2 = t1.AddDays(increment) + m2 = _distance_slope(+1, t2) + if m1 * m2 <= 0.0: + # There is a change of slope polarity within the time range [t1, t2]. + # Therefore this time range contains an apsis. + # Figure out whether it is perigee or apogee. + if m1 < 0.0 or m2 > 0.0: + # We found a minimum-distance event: perigee. + # Search the time range for the time when the slope goes from negative to positive. + apsis_time = Search(_distance_slope, +1, t1, t2, 1.0) + kind = APSIS_PERICENTER + elif m1 > 0.0 or m2 < 0.0: + # We found a maximum-distance event: apogee. + # Search the time range for the time when the slope goes from positive to negative. + apsis_time = Search(_distance_slope, -1, t1, t2, 1.0) + kind = APSIS_APOCENTER + else: + # This should never happen. It should not be possible for both slopes to be zero. + raise InternalError() + + if apsis_time is None: + # We should always be able to find a lunar apsis! + raise InternalError() + + dist = _MoonDistance(apsis_time) + return Apsis(apsis_time, kind, dist) + + # We have not yet found a slope polarity change. Keep searching. + t1 = t2 + m1 = m2 + iter += 1 + + # It should not be possible to fail to find an apsis within 2 synodic months. + raise InternalError() + + +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) + next = SearchLunarApsis(time) + # Verify that we found the opposite apsis from the previous one. + if apsis.kind == APSIS_APOCENTER: + expected = APSIS_PERICENTER + elif apsis.kind == APSIS_PERICENTER: + expected = APSIS_APOCENTER + else: + raise Error('Parameter "apsis" contains an invalid "kind" value.') + if next.kind != expected: + raise InternalError() + return next #================================================================================================== # + SearchSunLongitude @@ -2735,8 +2832,11 @@ def SearchRiseSet(body, observer, direction, startTime, limitDays): # + NextMoonQuarter # + SearchHourAngle # + SearchRiseSet -# - Seasons +# + Seasons +# + _FindSeasonChange # + Illumination # + SearchPeakMagnitude -# - SearchLunarApsis -# - NextLunarApsis +# + SearchLunarApsis +# + _distance_slope +# + _MoonDistance +# + NextLunarApsis