From d3621e7206fd57bf377990696d2fee5c3537c3d0 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Thu, 23 Sep 2021 14:27:56 -0400 Subject: [PATCH] Implemented Python function SearchAltitude. --- demo/python/astronomy.py | 152 +++++++++++++++++++++++---------- generate/template/astronomy.c | 4 +- generate/template/astronomy.py | 152 +++++++++++++++++++++++---------- generate/test.py | 39 +++++++++ pydown/py_prefix.md | 1 + source/c/README.md | 4 +- source/c/astronomy.c | 4 +- source/python/README.md | 32 +++++++ source/python/astronomy.py | 152 +++++++++++++++++++++++---------- 9 files changed, 405 insertions(+), 135 deletions(-) diff --git a/demo/python/astronomy.py b/demo/python/astronomy.py index 1ae73465..c52af2d0 100644 --- a/demo/python/astronomy.py +++ b/demo/python/astronomy.py @@ -5805,6 +5805,55 @@ class Direction(enum.Enum): Rise = +1 Set = -1 + +def _InternalSearchAltitude(body, observer, direction, startTime, limitDays, altitude_error_func, altitude_error_context): + if direction == Direction.Rise: + ha_before = 12.0 # minimum altitude (bottom) happens BEFORE the body rises. + ha_after = 0.0 # maximum altitude (culmination) happens AFTER the body rises. + elif direction == Direction.Set: + ha_before = 0.0 # culmination happens BEFORE the body sets. + ha_after = 12.0 # bottom happens AFTER the body sets. + else: + raise Error('Invalid value for direction parameter') + + # See if the body is currently above/below the horizon. + # If we are looking for next rise time and the body is below the horizon, + # we use the current time as the lower time bound and the next culmination + # as the upper bound. + # If the body is above the horizon, we search for the next bottom and use it + # as the lower bound and the next culmination after that bottom as the upper bound. + # The same logic applies for finding set times, only we swap the hour angles. + time_start = startTime + alt_before = altitude_error_func(altitude_error_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). + evt_before = SearchHourAngle(body, observer, ha_before, time_start) + time_before = evt_before.time + alt_before = altitude_error_func(altitude_error_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 + + evt_after = SearchHourAngle(body, observer, ha_after, time_before) + alt_after = altitude_error_func(altitude_error_context, evt_after.time) + + while True: + if alt_before <= 0.0 and alt_after > 0.0: + # Search between the "before time" and the "after time" for the desired event. + event_time = Search(altitude_error_func, altitude_error_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 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_before = evt_before.time + alt_before = altitude_error_func(altitude_error_context, evt_before.time) + alt_after = altitude_error_func(altitude_error_context, evt_after.time) + + class _peak_altitude_context: def __init__(self, body, direction, observer, body_radius_au): self.body = body @@ -5812,6 +5861,7 @@ class _peak_altitude_context: self.observer = observer self.body_radius_au = body_radius_au + def _peak_altitude(context, time): # Return the angular altitude above or below the horizon # of the highest part (the peak) of the given object. @@ -5829,6 +5879,7 @@ def _peak_altitude(context, time): alt = hor.altitude + math.degrees(context.body_radius_au / ofdate.dist) return context.direction.value * (alt + _REFRACTION_NEAR_HORIZON) + def SearchRiseSet(body, observer, direction, startTime, limitDays): """Searches for the next time a celestial body rises or sets as seen by an observer on the Earth. @@ -5881,53 +5932,68 @@ def SearchRiseSet(body, observer, direction, startTime, limitDays): else: body_radius = 0.0 - if direction == Direction.Rise: - ha_before = 12.0 # minimum altitude (bottom) happens BEFORE the body rises. - ha_after = 0.0 # maximum altitude (culmination) happens AFTER the body rises. - elif direction == Direction.Set: - ha_before = 0.0 # culmination happens BEFORE the body sets. - ha_after = 12.0 # bottom happens AFTER the body sets. - else: - raise Error('Invalid value for direction parameter') - context = _peak_altitude_context(body, direction, observer, body_radius) + return _InternalSearchAltitude(body, observer, direction, startTime, limitDays, _peak_altitude, context) - # See if the body is currently above/below the horizon. - # If we are looking for next rise time and the body is below the horizon, - # we use the current time as the lower time bound and the next culmination - # as the upper bound. - # If the body is above the horizon, we search for the next bottom and use it - # as the lower bound and the next culmination after that bottom as the upper bound. - # The same logic applies for finding set times, only we swap the hour angles. - time_start = startTime - 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). - 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 - evt_after = SearchHourAngle(body, observer, ha_after, time_before) - alt_after = _peak_altitude(context, evt_after.time) +class _altitude_error_context: + def __init__(self, body, direction, observer, altitude): + self.body = body + self.direction = direction + self.observer = observer + self.altitude = altitude - while True: - if alt_before <= 0.0 and alt_after > 0.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 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_before = evt_before.time - alt_before = _peak_altitude(context, evt_before.time) - alt_after = _peak_altitude(context, evt_after.time) +def _altitude_error_func(context, time): + ofdate = Equator(context.body, time, context.observer, True, True) + hor = Horizon(time, context.observer, ofdate.ra, ofdate.dec, Refraction.Airless) + return context.direction.value * (hor.altitude - context.altitude) + + +def SearchAltitude(body, observer, direction, dateStart, limitDays, altitude): + """Finds the next time a body reaches a given altitude. + + Finds when the given body ascends or descends through a given + altitude angle, as seen by an observer at the specified location on the Earth. + By using the appropriate combination of `direction` and `altitude` parameters, + this function can be used to find when civil, nautical, or astronomical twilight + begins (dawn) or ends (dusk). + + Civil dawn begins before sunrise when the Sun ascends through 6 degrees below + the horizon. To find civil dawn, pass `Direction.Rise` for `direction` and -6 for `altitude`. + + Civil dusk ends after sunset when the Sun descends through 6 degrees below the horizon. + To find civil dusk, pass `Direction.Set` for `direction` and -6 for `altitude`. + + Nautical twilight is similar to civil twilight, only the `altitude` value should be -12 degrees. + + Astronomical twilight uses -18 degrees as the `altitude` value. + + Parameters + ---------- + body : Body + The Sun, Moon, or any planet other than the Earth. + observer : Observer + The location where observation takes place. + direction : Direction + Either `Direction.Rise` to find an ascending altitude event + or `Direction.Set` to find a descending altitude event. + startTime : Time + The date and time at which to start the search. + limitDays : float + The fractional number of days after `dateStart` that limits + when the altitude event is to be found. Must be a positive number. + + Returns + ------- + Time or `None` + If the altitude event time is found within the specified time window, + this function returns that time. Otherwise, it returns `None`. + """ + if not (-90.0 <= altitude <= +90.0): + raise Error('Invalid altitude: {}'.format(altitude)) + + context = _altitude_error_context(body, direction, observer, altitude) + return _InternalSearchAltitude(body, observer, direction, dateStart, limitDays, _altitude_error_func, context) class SeasonInfo: """The dates and times of changes of season for a given calendar year. diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index 28605aad..30acd1ba 100644 --- a/generate/template/astronomy.c +++ b/generate/template/astronomy.c @@ -4957,10 +4957,10 @@ astro_search_result_t Astronomy_SearchRiseSet( * begins (dawn) or ends (dusk). * * Civil dawn begins before sunrise when the Sun ascends through 6 degrees below - * the horizon. To find civil dawn, pass +1 for `direction` and -6 for `altitude`. + * the horizon. To find civil dawn, pass `DIRECTION_RISE` for `direction` and -6 for `altitude`. * * Civil dusk ends after sunset when the Sun descends through 6 degrees below the horizon. - * To find civil dusk, pass -1 for `direction` and -6 for `altitude`. + * To find civil dusk, pass `DIRECTION_SET` for `direction` and -6 for `altitude`. * * Nautical twilight is similar to civil twilight, only the `altitude` value should be -12 degrees. * diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index e2956393..edb63444 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -3775,6 +3775,55 @@ class Direction(enum.Enum): Rise = +1 Set = -1 + +def _InternalSearchAltitude(body, observer, direction, startTime, limitDays, altitude_error_func, altitude_error_context): + if direction == Direction.Rise: + ha_before = 12.0 # minimum altitude (bottom) happens BEFORE the body rises. + ha_after = 0.0 # maximum altitude (culmination) happens AFTER the body rises. + elif direction == Direction.Set: + ha_before = 0.0 # culmination happens BEFORE the body sets. + ha_after = 12.0 # bottom happens AFTER the body sets. + else: + raise Error('Invalid value for direction parameter') + + # See if the body is currently above/below the horizon. + # If we are looking for next rise time and the body is below the horizon, + # we use the current time as the lower time bound and the next culmination + # as the upper bound. + # If the body is above the horizon, we search for the next bottom and use it + # as the lower bound and the next culmination after that bottom as the upper bound. + # The same logic applies for finding set times, only we swap the hour angles. + time_start = startTime + alt_before = altitude_error_func(altitude_error_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). + evt_before = SearchHourAngle(body, observer, ha_before, time_start) + time_before = evt_before.time + alt_before = altitude_error_func(altitude_error_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 + + evt_after = SearchHourAngle(body, observer, ha_after, time_before) + alt_after = altitude_error_func(altitude_error_context, evt_after.time) + + while True: + if alt_before <= 0.0 and alt_after > 0.0: + # Search between the "before time" and the "after time" for the desired event. + event_time = Search(altitude_error_func, altitude_error_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 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_before = evt_before.time + alt_before = altitude_error_func(altitude_error_context, evt_before.time) + alt_after = altitude_error_func(altitude_error_context, evt_after.time) + + class _peak_altitude_context: def __init__(self, body, direction, observer, body_radius_au): self.body = body @@ -3782,6 +3831,7 @@ class _peak_altitude_context: self.observer = observer self.body_radius_au = body_radius_au + def _peak_altitude(context, time): # Return the angular altitude above or below the horizon # of the highest part (the peak) of the given object. @@ -3799,6 +3849,7 @@ def _peak_altitude(context, time): alt = hor.altitude + math.degrees(context.body_radius_au / ofdate.dist) return context.direction.value * (alt + _REFRACTION_NEAR_HORIZON) + def SearchRiseSet(body, observer, direction, startTime, limitDays): """Searches for the next time a celestial body rises or sets as seen by an observer on the Earth. @@ -3851,53 +3902,68 @@ def SearchRiseSet(body, observer, direction, startTime, limitDays): else: body_radius = 0.0 - if direction == Direction.Rise: - ha_before = 12.0 # minimum altitude (bottom) happens BEFORE the body rises. - ha_after = 0.0 # maximum altitude (culmination) happens AFTER the body rises. - elif direction == Direction.Set: - ha_before = 0.0 # culmination happens BEFORE the body sets. - ha_after = 12.0 # bottom happens AFTER the body sets. - else: - raise Error('Invalid value for direction parameter') - context = _peak_altitude_context(body, direction, observer, body_radius) + return _InternalSearchAltitude(body, observer, direction, startTime, limitDays, _peak_altitude, context) - # See if the body is currently above/below the horizon. - # If we are looking for next rise time and the body is below the horizon, - # we use the current time as the lower time bound and the next culmination - # as the upper bound. - # If the body is above the horizon, we search for the next bottom and use it - # as the lower bound and the next culmination after that bottom as the upper bound. - # The same logic applies for finding set times, only we swap the hour angles. - time_start = startTime - 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). - 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 - evt_after = SearchHourAngle(body, observer, ha_after, time_before) - alt_after = _peak_altitude(context, evt_after.time) +class _altitude_error_context: + def __init__(self, body, direction, observer, altitude): + self.body = body + self.direction = direction + self.observer = observer + self.altitude = altitude - while True: - if alt_before <= 0.0 and alt_after > 0.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 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_before = evt_before.time - alt_before = _peak_altitude(context, evt_before.time) - alt_after = _peak_altitude(context, evt_after.time) +def _altitude_error_func(context, time): + ofdate = Equator(context.body, time, context.observer, True, True) + hor = Horizon(time, context.observer, ofdate.ra, ofdate.dec, Refraction.Airless) + return context.direction.value * (hor.altitude - context.altitude) + + +def SearchAltitude(body, observer, direction, dateStart, limitDays, altitude): + """Finds the next time a body reaches a given altitude. + + Finds when the given body ascends or descends through a given + altitude angle, as seen by an observer at the specified location on the Earth. + By using the appropriate combination of `direction` and `altitude` parameters, + this function can be used to find when civil, nautical, or astronomical twilight + begins (dawn) or ends (dusk). + + Civil dawn begins before sunrise when the Sun ascends through 6 degrees below + the horizon. To find civil dawn, pass `Direction.Rise` for `direction` and -6 for `altitude`. + + Civil dusk ends after sunset when the Sun descends through 6 degrees below the horizon. + To find civil dusk, pass `Direction.Set` for `direction` and -6 for `altitude`. + + Nautical twilight is similar to civil twilight, only the `altitude` value should be -12 degrees. + + Astronomical twilight uses -18 degrees as the `altitude` value. + + Parameters + ---------- + body : Body + The Sun, Moon, or any planet other than the Earth. + observer : Observer + The location where observation takes place. + direction : Direction + Either `Direction.Rise` to find an ascending altitude event + or `Direction.Set` to find a descending altitude event. + startTime : Time + The date and time at which to start the search. + limitDays : float + The fractional number of days after `dateStart` that limits + when the altitude event is to be found. Must be a positive number. + + Returns + ------- + Time or `None` + If the altitude event time is found within the specified time window, + this function returns that time. Otherwise, it returns `None`. + """ + if not (-90.0 <= altitude <= +90.0): + raise Error('Invalid altitude: {}'.format(altitude)) + + context = _altitude_error_context(body, direction, observer, altitude) + return _InternalSearchAltitude(body, observer, direction, dateStart, limitDays, _altitude_error_func, context) class SeasonInfo: """The dates and times of changes of season for a given calendar year. diff --git a/generate/test.py b/generate/test.py index 63162c32..7bc15309 100755 --- a/generate/test.py +++ b/generate/test.py @@ -2051,6 +2051,44 @@ def Aberration(): #----------------------------------------------------------------------------------------------------------- +def Twilight(): + tolerance_seconds = 60.0 + max_diff = 0.0 + filename = 'riseset/twilight.txt' + with open(filename, 'rt') as infile: + lnum = 0 + for line in infile: + lnum += 1 + tokens = line.split() + if len(tokens) != 9: + print('PY Twilight({} line {}): incorrect number of tokens = {}'.format(filename, lnum, len(tokens))) + return 1 + observer = astronomy.Observer(float(tokens[0]), float(tokens[1]), 0.0) + searchDate = astronomy.Time.Parse(tokens[2]) + correctTimes = [astronomy.Time.Parse(t) for t in tokens[3:]] + calcTimes = [ + astronomy.SearchAltitude(astronomy.Body.Sun, observer, astronomy.Direction.Rise, searchDate, 1.0, -18.0), # astronomical dawn + astronomy.SearchAltitude(astronomy.Body.Sun, observer, astronomy.Direction.Rise, searchDate, 1.0, -12.0), # nautical dawn + astronomy.SearchAltitude(astronomy.Body.Sun, observer, astronomy.Direction.Rise, searchDate, 1.0, -6.0), # civil dawn + astronomy.SearchAltitude(astronomy.Body.Sun, observer, astronomy.Direction.Set, searchDate, 1.0, -6.0), # civil dusk + astronomy.SearchAltitude(astronomy.Body.Sun, observer, astronomy.Direction.Set, searchDate, 1.0, -12.0), # nautical dusk + astronomy.SearchAltitude(astronomy.Body.Sun, observer, astronomy.Direction.Set, searchDate, 1.0, -18.0) # astronomical dusk + ] + for i in range(6): + correct = correctTimes[i] + calc = calcTimes[i] + diff = 86400.0 * vabs(calc.ut - correct.ut) + if diff > tolerance_seconds: + print('PY Twilight({} line {}): EXCESSIVE ERROR = {} seconds for case {}'.format(filename, lnum, diff, i)) + return 1 + if diff > max_diff: + max_diff = diff + + print('PY Twilight: PASS ({} test cases, max error = {} seconds)'.format(lnum, max_diff)) + return 0 + +#----------------------------------------------------------------------------------------------------------- + UnitTests = { 'aberration': Aberration, 'barystate': BaryState, @@ -2075,6 +2113,7 @@ UnitTests = { 'seasons': Seasons, 'time': AstroTime, 'transit': Transit, + 'twilight': Twilight, } #----------------------------------------------------------------------------------------------------------- diff --git a/pydown/py_prefix.md b/pydown/py_prefix.md index 203dee81..849bda76 100644 --- a/pydown/py_prefix.md +++ b/pydown/py_prefix.md @@ -50,6 +50,7 @@ To get started quickly, here are some [examples](../../demo/python/). | Function | Description | | -------- | ----------- | | [SearchRiseSet](#SearchRiseSet) | Finds time of rise or set for a body as seen by an observer on the Earth. | +| [SearchAltitude](#SearchAltitude) | Finds time when a body reaches a given altitude above or below the horizon. Useful for finding civil, nautical, or astronomical twilight. | | [SearchHourAngle](#SearchHourAngle) | Finds when body reaches a given hour angle for an observer on the Earth. Hour angle = 0 finds culmination, the highest point in the sky. | ### Moon phases diff --git a/source/c/README.md b/source/c/README.md index d1c9cffc..0462939c 100644 --- a/source/c/README.md +++ b/source/c/README.md @@ -1739,9 +1739,9 @@ If the search does not converge within 20 iterations, it will fail with status c Finds when the given body ascends or descends through a given altitude angle, as seen by an observer at the specified location on the Earth. By using the appropriate combination of `direction` and `altitude` parameters, this function can be used to find when civil, nautical, or astronomical twilight begins (dawn) or ends (dusk). -Civil dawn begins before sunrise when the Sun ascends through 6 degrees below the horizon. To find civil dawn, pass +1 for `direction` and -6 for `altitude`. +Civil dawn begins before sunrise when the Sun ascends through 6 degrees below the horizon. To find civil dawn, pass `DIRECTION_RISE` for `direction` and -6 for `altitude`. -Civil dusk ends after sunset when the Sun descends through 6 degrees below the horizon. To find civil dusk, pass -1 for `direction` and -6 for `altitude`. +Civil dusk ends after sunset when the Sun descends through 6 degrees below the horizon. To find civil dusk, pass `DIRECTION_SET` for `direction` and -6 for `altitude`. Nautical twilight is similar to civil twilight, only the `altitude` value should be -12 degrees. diff --git a/source/c/astronomy.c b/source/c/astronomy.c index 82fa0c51..d0c5566d 100644 --- a/source/c/astronomy.c +++ b/source/c/astronomy.c @@ -6183,10 +6183,10 @@ astro_search_result_t Astronomy_SearchRiseSet( * begins (dawn) or ends (dusk). * * Civil dawn begins before sunrise when the Sun ascends through 6 degrees below - * the horizon. To find civil dawn, pass +1 for `direction` and -6 for `altitude`. + * the horizon. To find civil dawn, pass `DIRECTION_RISE` for `direction` and -6 for `altitude`. * * Civil dusk ends after sunset when the Sun descends through 6 degrees below the horizon. - * To find civil dusk, pass -1 for `direction` and -6 for `altitude`. + * To find civil dusk, pass `DIRECTION_SET` for `direction` and -6 for `altitude`. * * Nautical twilight is similar to civil twilight, only the `altitude` value should be -12 degrees. * diff --git a/source/python/README.md b/source/python/README.md index c93c6292..22e9a4ef 100644 --- a/source/python/README.md +++ b/source/python/README.md @@ -50,6 +50,7 @@ To get started quickly, here are some [examples](../../demo/python/). | Function | Description | | -------- | ----------- | | [SearchRiseSet](#SearchRiseSet) | Finds time of rise or set for a body as seen by an observer on the Earth. | +| [SearchAltitude](#SearchAltitude) | Finds time when a body reaches a given altitude above or below the horizon. Useful for finding civil, nautical, or astronomical twilight. | | [SearchHourAngle](#SearchHourAngle) | Finds when body reaches a given hour angle for an observer on the Earth. Hour angle = 0 finds culmination, the highest point in the sky. | ### Moon phases @@ -2142,6 +2143,37 @@ the function returns `None`. --- + +### SearchAltitude(body, observer, direction, dateStart, limitDays, altitude) + +**Finds the next time a body reaches a given altitude.** + +Finds when the given body ascends or descends through a given +altitude angle, as seen by an observer at the specified location on the Earth. +By using the appropriate combination of `direction` and `altitude` parameters, +this function can be used to find when civil, nautical, or astronomical twilight +begins (dawn) or ends (dusk). +Civil dawn begins before sunrise when the Sun ascends through 6 degrees below +the horizon. To find civil dawn, pass `Direction.Rise` for `direction` and -6 for `altitude`. +Civil dusk ends after sunset when the Sun descends through 6 degrees below the horizon. +To find civil dusk, pass `Direction.Set` for `direction` and -6 for `altitude`. +Nautical twilight is similar to civil twilight, only the `altitude` value should be -12 degrees. +Astronomical twilight uses -18 degrees as the `altitude` value. + +| Type | Parameter | Description | +| --- | --- | --- | +| [`Body`](#Body) | `body` | The Sun, Moon, or any planet other than the Earth. | +| [`Observer`](#Observer) | `observer` | The location where observation takes place. | +| [`Direction`](#Direction) | `direction` | Either `Direction.Rise` to find an ascending altitude event or `Direction.Set` to find a descending altitude event. | +| [`Time`](#Time) | `startTime` | The date and time at which to start the search. | +| `float` | `limitDays` | The fractional number of days after `dateStart` that limits when the altitude event is to be found. Must be a positive number. | + +### Returns: [`Time`](#Time) or `None` +If the altitude event time is found within the specified time window, +this function returns that time. Otherwise, it returns `None`. + +--- + ### SearchGlobalSolarEclipse(startTime) diff --git a/source/python/astronomy.py b/source/python/astronomy.py index 1ae73465..c52af2d0 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -5805,6 +5805,55 @@ class Direction(enum.Enum): Rise = +1 Set = -1 + +def _InternalSearchAltitude(body, observer, direction, startTime, limitDays, altitude_error_func, altitude_error_context): + if direction == Direction.Rise: + ha_before = 12.0 # minimum altitude (bottom) happens BEFORE the body rises. + ha_after = 0.0 # maximum altitude (culmination) happens AFTER the body rises. + elif direction == Direction.Set: + ha_before = 0.0 # culmination happens BEFORE the body sets. + ha_after = 12.0 # bottom happens AFTER the body sets. + else: + raise Error('Invalid value for direction parameter') + + # See if the body is currently above/below the horizon. + # If we are looking for next rise time and the body is below the horizon, + # we use the current time as the lower time bound and the next culmination + # as the upper bound. + # If the body is above the horizon, we search for the next bottom and use it + # as the lower bound and the next culmination after that bottom as the upper bound. + # The same logic applies for finding set times, only we swap the hour angles. + time_start = startTime + alt_before = altitude_error_func(altitude_error_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). + evt_before = SearchHourAngle(body, observer, ha_before, time_start) + time_before = evt_before.time + alt_before = altitude_error_func(altitude_error_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 + + evt_after = SearchHourAngle(body, observer, ha_after, time_before) + alt_after = altitude_error_func(altitude_error_context, evt_after.time) + + while True: + if alt_before <= 0.0 and alt_after > 0.0: + # Search between the "before time" and the "after time" for the desired event. + event_time = Search(altitude_error_func, altitude_error_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 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_before = evt_before.time + alt_before = altitude_error_func(altitude_error_context, evt_before.time) + alt_after = altitude_error_func(altitude_error_context, evt_after.time) + + class _peak_altitude_context: def __init__(self, body, direction, observer, body_radius_au): self.body = body @@ -5812,6 +5861,7 @@ class _peak_altitude_context: self.observer = observer self.body_radius_au = body_radius_au + def _peak_altitude(context, time): # Return the angular altitude above or below the horizon # of the highest part (the peak) of the given object. @@ -5829,6 +5879,7 @@ def _peak_altitude(context, time): alt = hor.altitude + math.degrees(context.body_radius_au / ofdate.dist) return context.direction.value * (alt + _REFRACTION_NEAR_HORIZON) + def SearchRiseSet(body, observer, direction, startTime, limitDays): """Searches for the next time a celestial body rises or sets as seen by an observer on the Earth. @@ -5881,53 +5932,68 @@ def SearchRiseSet(body, observer, direction, startTime, limitDays): else: body_radius = 0.0 - if direction == Direction.Rise: - ha_before = 12.0 # minimum altitude (bottom) happens BEFORE the body rises. - ha_after = 0.0 # maximum altitude (culmination) happens AFTER the body rises. - elif direction == Direction.Set: - ha_before = 0.0 # culmination happens BEFORE the body sets. - ha_after = 12.0 # bottom happens AFTER the body sets. - else: - raise Error('Invalid value for direction parameter') - context = _peak_altitude_context(body, direction, observer, body_radius) + return _InternalSearchAltitude(body, observer, direction, startTime, limitDays, _peak_altitude, context) - # See if the body is currently above/below the horizon. - # If we are looking for next rise time and the body is below the horizon, - # we use the current time as the lower time bound and the next culmination - # as the upper bound. - # If the body is above the horizon, we search for the next bottom and use it - # as the lower bound and the next culmination after that bottom as the upper bound. - # The same logic applies for finding set times, only we swap the hour angles. - time_start = startTime - 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). - 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 - evt_after = SearchHourAngle(body, observer, ha_after, time_before) - alt_after = _peak_altitude(context, evt_after.time) +class _altitude_error_context: + def __init__(self, body, direction, observer, altitude): + self.body = body + self.direction = direction + self.observer = observer + self.altitude = altitude - while True: - if alt_before <= 0.0 and alt_after > 0.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 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_before = evt_before.time - alt_before = _peak_altitude(context, evt_before.time) - alt_after = _peak_altitude(context, evt_after.time) +def _altitude_error_func(context, time): + ofdate = Equator(context.body, time, context.observer, True, True) + hor = Horizon(time, context.observer, ofdate.ra, ofdate.dec, Refraction.Airless) + return context.direction.value * (hor.altitude - context.altitude) + + +def SearchAltitude(body, observer, direction, dateStart, limitDays, altitude): + """Finds the next time a body reaches a given altitude. + + Finds when the given body ascends or descends through a given + altitude angle, as seen by an observer at the specified location on the Earth. + By using the appropriate combination of `direction` and `altitude` parameters, + this function can be used to find when civil, nautical, or astronomical twilight + begins (dawn) or ends (dusk). + + Civil dawn begins before sunrise when the Sun ascends through 6 degrees below + the horizon. To find civil dawn, pass `Direction.Rise` for `direction` and -6 for `altitude`. + + Civil dusk ends after sunset when the Sun descends through 6 degrees below the horizon. + To find civil dusk, pass `Direction.Set` for `direction` and -6 for `altitude`. + + Nautical twilight is similar to civil twilight, only the `altitude` value should be -12 degrees. + + Astronomical twilight uses -18 degrees as the `altitude` value. + + Parameters + ---------- + body : Body + The Sun, Moon, or any planet other than the Earth. + observer : Observer + The location where observation takes place. + direction : Direction + Either `Direction.Rise` to find an ascending altitude event + or `Direction.Set` to find a descending altitude event. + startTime : Time + The date and time at which to start the search. + limitDays : float + The fractional number of days after `dateStart` that limits + when the altitude event is to be found. Must be a positive number. + + Returns + ------- + Time or `None` + If the altitude event time is found within the specified time window, + this function returns that time. Otherwise, it returns `None`. + """ + if not (-90.0 <= altitude <= +90.0): + raise Error('Invalid altitude: {}'.format(altitude)) + + context = _altitude_error_context(body, direction, observer, altitude) + return _InternalSearchAltitude(body, observer, direction, dateStart, limitDays, _altitude_error_func, context) class SeasonInfo: """The dates and times of changes of season for a given calendar year.