From c85881c25b752abe2f47c4f4264d168bd143f684 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Mon, 13 Mar 2023 20:42:46 -0400 Subject: [PATCH] Python: added `metersAboveGround` parameter to SearchRiseSet. --- demo/python/astronomy.py | 56 +++++++++++++++++++++++- generate/template/astronomy.py | 56 +++++++++++++++++++++++- generate/test.py | 65 ++++++++++++++++++++++++++++ source/python/README.md | 6 ++- source/python/astronomy/astronomy.py | 56 +++++++++++++++++++++++- 5 files changed, 232 insertions(+), 7 deletions(-) diff --git a/demo/python/astronomy.py b/demo/python/astronomy.py index ab3d5451..b413f9ee 100644 --- a/demo/python/astronomy.py +++ b/demo/python/astronomy.py @@ -4958,6 +4958,9 @@ def RefractionAngle(refraction: Refraction, altitude: float) -> float: the amount of "lift" caused by atmospheric refraction. This is the number of degrees higher in the sky an object appears due to lensing of the Earth's atmosphere. + This function works best near sea level. + To correct for higher elevations, call #Atmosphere for that + elevation and multiply the refraction angle by the resulting relative density. Parameters ---------- @@ -6470,7 +6473,34 @@ def _InternalSearchAltitude(body: Body, observer: Observer, direction: Direction a1 = a2 -def SearchRiseSet(body: Body, observer: Observer, direction: Direction, startTime: Time, limitDays: float) -> Optional[Time]: +def _HorizonDipAngle(observer: Observer, metersAboveGround: float) -> float: + # Calculate the effective radius of the Earth at ground level below the observer. + # Correct for the Earth's oblateness. + phi = math.radians(observer.latitude) + sinphi = math.sin(phi) + cosphi = math.cos(phi) + c = 1.0 / math.hypot(cosphi, sinphi*_EARTH_FLATTENING) + s = c * (_EARTH_FLATTENING * _EARTH_FLATTENING) + ht_km = (observer.height - metersAboveGround) / 1000.0 # height of ground above sea level + ach = _EARTH_EQUATORIAL_RADIUS_KM*c + ht_km + ash = _EARTH_EQUATORIAL_RADIUS_KM*s + ht_km + radius_m = 1000.0 * math.hypot(ach*cosphi, ash*sinphi) + + # Correct refraction of a ray of light traveling tangent to the Earth's surface. + # Based on: https://www.largeformatphotography.info/sunmooncalc/SMCalc.js + # which in turn derives from: + # Sweer, John. 1938. The Path of a Ray of Light Tangent to the Surface of the Earth. + # Journal of the Optical Society of America 28 (September):327-329. + + # k = refraction index + k = 0.175 * (1.0 - (6.5e-3/283.15)*(observer.height - (2.0/3.0)*metersAboveGround))**3.256 + + # Calculate how far below the observer's horizontal plane the observed horizon dips. + return math.degrees(-(math.sqrt(2*(1 - k)*metersAboveGround / radius_m) / (1 - k))) + + + +def SearchRiseSet(body: Body, observer: Observer, direction: Direction, startTime: Time, limitDays: float, metersAboveGround: float = 0.0) -> Optional[Time]: """Searches for the next time a celestial body rises or sets as seen by an observer on the Earth. This function finds the next rise or set time of the Sun, Moon, or planet other than the Earth. @@ -6515,6 +6545,13 @@ def SearchRiseSet(body: Body, observer: Observer, direction: Direction, startTim In cases where you want to find the next rise or set time no matter how far in the future (for example, for an observer near the south pole), you can pass in a larger value like 365. + metersAboveGround : float + Default value = 0.0. + Usually the observer is located at ground level. Then this parameter + should be zero. But if the observer is significantly higher than ground + level, for example in an airplane, this parameter should be a positive + number indicating how far above the ground the observer is. + An exception occurs if `metersAboveGround` is negative. Returns ------- @@ -6522,13 +6559,28 @@ def SearchRiseSet(body: Body, observer: Observer, direction: Direction, startTim If the rise or set time is found within the specified time window, this function returns that time. Otherwise, it returns `None`. """ + if not math.isfinite(metersAboveGround) or metersAboveGround < 0.0: + raise Error('Invalid value for metersAboveGround: {}'.format(metersAboveGround)) + + # Determine the radius of the body to be observed. if body == Body.Sun: bodyRadiusAu = _SUN_RADIUS_AU elif body == Body.Moon: bodyRadiusAu = _MOON_EQUATORIAL_RADIUS_AU else: bodyRadiusAu = 0.0 - return _InternalSearchAltitude(body, observer, direction, startTime, limitDays, bodyRadiusAu, -_REFRACTION_NEAR_HORIZON) + + # Calculate atmospheric density at ground level. + atmos = Atmosphere(observer.height - metersAboveGround) + + # Calculate the apparent angular dip of the horizon. + dip = _HorizonDipAngle(observer, metersAboveGround) + + # Correct refraction for objects near the horizon, using atmospheric density at the ground. + altitude = dip - (_REFRACTION_NEAR_HORIZON * atmos.density) + + # Search for the top of the body crossing the corrected altitude angle. + return _InternalSearchAltitude(body, observer, direction, startTime, limitDays, bodyRadiusAu, altitude) def SearchAltitude(body: Body, observer: Observer, direction: Direction, startTime: Time, limitDays: float, altitude: float) -> Optional[Time]: diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index a511ccbf..baae3c15 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -3452,6 +3452,9 @@ def RefractionAngle(refraction: Refraction, altitude: float) -> float: the amount of "lift" caused by atmospheric refraction. This is the number of degrees higher in the sky an object appears due to lensing of the Earth's atmosphere. + This function works best near sea level. + To correct for higher elevations, call #Atmosphere for that + elevation and multiply the refraction angle by the resulting relative density. Parameters ---------- @@ -4964,7 +4967,34 @@ def _InternalSearchAltitude(body: Body, observer: Observer, direction: Direction a1 = a2 -def SearchRiseSet(body: Body, observer: Observer, direction: Direction, startTime: Time, limitDays: float) -> Optional[Time]: +def _HorizonDipAngle(observer: Observer, metersAboveGround: float) -> float: + # Calculate the effective radius of the Earth at ground level below the observer. + # Correct for the Earth's oblateness. + phi = math.radians(observer.latitude) + sinphi = math.sin(phi) + cosphi = math.cos(phi) + c = 1.0 / math.hypot(cosphi, sinphi*_EARTH_FLATTENING) + s = c * (_EARTH_FLATTENING * _EARTH_FLATTENING) + ht_km = (observer.height - metersAboveGround) / 1000.0 # height of ground above sea level + ach = _EARTH_EQUATORIAL_RADIUS_KM*c + ht_km + ash = _EARTH_EQUATORIAL_RADIUS_KM*s + ht_km + radius_m = 1000.0 * math.hypot(ach*cosphi, ash*sinphi) + + # Correct refraction of a ray of light traveling tangent to the Earth's surface. + # Based on: https://www.largeformatphotography.info/sunmooncalc/SMCalc.js + # which in turn derives from: + # Sweer, John. 1938. The Path of a Ray of Light Tangent to the Surface of the Earth. + # Journal of the Optical Society of America 28 (September):327-329. + + # k = refraction index + k = 0.175 * (1.0 - (6.5e-3/283.15)*(observer.height - (2.0/3.0)*metersAboveGround))**3.256 + + # Calculate how far below the observer's horizontal plane the observed horizon dips. + return math.degrees(-(math.sqrt(2*(1 - k)*metersAboveGround / radius_m) / (1 - k))) + + + +def SearchRiseSet(body: Body, observer: Observer, direction: Direction, startTime: Time, limitDays: float, metersAboveGround: float = 0.0) -> Optional[Time]: """Searches for the next time a celestial body rises or sets as seen by an observer on the Earth. This function finds the next rise or set time of the Sun, Moon, or planet other than the Earth. @@ -5009,6 +5039,13 @@ def SearchRiseSet(body: Body, observer: Observer, direction: Direction, startTim In cases where you want to find the next rise or set time no matter how far in the future (for example, for an observer near the south pole), you can pass in a larger value like 365. + metersAboveGround : float + Default value = 0.0. + Usually the observer is located at ground level. Then this parameter + should be zero. But if the observer is significantly higher than ground + level, for example in an airplane, this parameter should be a positive + number indicating how far above the ground the observer is. + An exception occurs if `metersAboveGround` is negative. Returns ------- @@ -5016,13 +5053,28 @@ def SearchRiseSet(body: Body, observer: Observer, direction: Direction, startTim If the rise or set time is found within the specified time window, this function returns that time. Otherwise, it returns `None`. """ + if not math.isfinite(metersAboveGround) or metersAboveGround < 0.0: + raise Error('Invalid value for metersAboveGround: {}'.format(metersAboveGround)) + + # Determine the radius of the body to be observed. if body == Body.Sun: bodyRadiusAu = _SUN_RADIUS_AU elif body == Body.Moon: bodyRadiusAu = _MOON_EQUATORIAL_RADIUS_AU else: bodyRadiusAu = 0.0 - return _InternalSearchAltitude(body, observer, direction, startTime, limitDays, bodyRadiusAu, -_REFRACTION_NEAR_HORIZON) + + # Calculate atmospheric density at ground level. + atmos = Atmosphere(observer.height - metersAboveGround) + + # Calculate the apparent angular dip of the horizon. + dip = _HorizonDipAngle(observer, metersAboveGround) + + # Correct refraction for objects near the horizon, using atmospheric density at the ground. + altitude = dip - (_REFRACTION_NEAR_HORIZON * atmos.density) + + # Search for the top of the body crossing the corrected altitude angle. + return _InternalSearchAltitude(body, observer, direction, startTime, limitDays, bodyRadiusAu, altitude) def SearchAltitude(body: Body, observer: Observer, direction: Direction, startTime: Time, limitDays: float, altitude: float) -> Optional[Time]: diff --git a/generate/test.py b/generate/test.py index 5ff16208..23d7846f 100755 --- a/generate/test.py +++ b/generate/test.py @@ -3264,6 +3264,70 @@ def Atmosphere(): #----------------------------------------------------------------------------------------------------------- +def RiseSetElevationBodyCase(body, observer, direction, metersAboveGround, startTime, eventOffsetDays): + time = astronomy.SearchRiseSet(body, observer, direction, startTime, 2.0, metersAboveGround) + if not time: + return Fail('RiseSetElevationBodyCase {} {}: search failed.'.format(body, direction)) + diff = v(time.ut - (startTime.ut + eventOffsetDays)) + if diff > 0.5: + diff -= 1.0 # assume event actually takes place on the next day + diff = vabs(MINUTES_PER_DAY * diff) # convert signed days to absolute minutes + if diff > 0.5: + return Fail('RiseSetElevationBodyCase {} {}: EXCESSIVE diff = {}.'.format(body, direction, diff)) + return 0 + + +def RiseSetElevation(): + regex = re.compile(r'^(\d+)-(\d+)-(\d+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\d+):(\d+)\s+(\d+):(\d+)\s+(\d+):(\d+)\s+(\d+):(\d+)\s+(\S+)\s*$') + filename = 'riseset/elevation.txt' + with open(filename, 'rt') as infile: + lnum = 0 + for line in infile: + lnum += 1 + if line.startswith('#'): + continue + m = regex.match(line) + if not m: + return Fail('RiseSetElevation({} line {})'.format(filename, lnum), 'Invalid data format') + year = int(m.group(1)) + month = int(m.group(2)) + day = int(m.group(3)) + latitude = v(float(m.group(4))) + longitude = v(float(m.group(5))) + height = v(float(m.group(6))) + metersAboveGround = v(float(m.group(7))) + srh = int(m.group( 8)) + srm = int(m.group( 9)) + ssh = int(m.group(10)) + ssm = int(m.group(11)) + mrh = int(m.group(12)) + mrm = int(m.group(13)) + msh = int(m.group(14)) + msm = int(m.group(15)) + + # Get search origin time + time = astronomy.Time.Make(year, month, day, 0, 0, 0.0) + + # Convert scanned values into sunrise, sunset, moonrise, moonset day offsets. + sr = (srh + srm/60.0) / 24.0 + ss = (ssh + ssm/60.0) / 24.0 + mr = (mrh + mrm/60.0) / 24.0 + ms = (msh + msm/60.0) / 24.0 + + observer = astronomy.Observer(latitude, longitude, height) + + if (0 != RiseSetElevationBodyCase(astronomy.Body.Sun, observer, astronomy.Direction.Rise, metersAboveGround, time, sr) or + 0 != RiseSetElevationBodyCase(astronomy.Body.Sun, observer, astronomy.Direction.Set, metersAboveGround, time, ss) or + 0 != RiseSetElevationBodyCase(astronomy.Body.Moon, observer, astronomy.Direction.Rise, metersAboveGround, time, mr) or + 0 != RiseSetElevationBodyCase(astronomy.Body.Moon, observer, astronomy.Direction.Set, metersAboveGround, time, ms)): + return 1 + + + return Pass('RiseSetElevation') + +#----------------------------------------------------------------------------------------------------------- + + UnitTests = { 'aberration': Aberration, 'atmosphere': Atmosphere, @@ -3297,6 +3361,7 @@ UnitTests = { 'refraction': Refraction, 'repr': Repr, 'riseset': RiseSet, + 'riseset_elevation': RiseSetElevation, 'riseset_reverse': RiseSetReverse, 'rotation': Rotation, 'seasons': Seasons, diff --git a/source/python/README.md b/source/python/README.md index 26d49857..c184d81d 100644 --- a/source/python/README.md +++ b/source/python/README.md @@ -2494,6 +2494,9 @@ Given an altitude angle and a refraction option, calculates the amount of "lift" caused by atmospheric refraction. This is the number of degrees higher in the sky an object appears due to lensing of the Earth's atmosphere. +This function works best near sea level. +To correct for higher elevations, call [`Atmosphere`](#Atmosphere) for that +elevation and multiply the refraction angle by the resulting relative density. | Type | Parameter | Description | | --- | --- | --- | @@ -3324,7 +3327,7 @@ The date and time of the relative longitude event. --- -### SearchRiseSet(body: [`Body`](#Body), observer: [`Observer`](#Observer), direction: [`Direction`](#Direction), startTime: [`Time`](#Time), limitDays: float) → Optional\[[`Time`](#Time)\] +### SearchRiseSet(body: [`Body`](#Body), observer: [`Observer`](#Observer), direction: [`Direction`](#Direction), startTime: [`Time`](#Time), limitDays: float, metersAboveGround: float = 0.0) → Optional\[[`Time`](#Time)\] **Searches for the next time a celestial body rises or sets as seen by an observer on the Earth.** @@ -3355,6 +3358,7 @@ Therefore callers must not assume that the function will always succeed. | [`Direction`](#Direction) | `direction` | Either `Direction.Rise` to find a rise time or `Direction.Set` to find a set time. | | [`Time`](#Time) | `startTime` | The date and time at which to start the search. | | `float` | `limitDays` | Limits how many days to search for a rise or set time, and defines the direction in time to search. When `limitDays` is positive, the search is performed into the future, after `startTime`. When negative, the search is performed into the past, before `startTime`. To limit a rise or set time to the same day, you can use a value of 1 day. In cases where you want to find the next rise or set time no matter how far in the future (for example, for an observer near the south pole), you can pass in a larger value like 365. | +| `float` | `metersAboveGround` | Default value = 0.0. Usually the observer is located at ground level. Then this parameter should be zero. But if the observer is significantly higher than ground level, for example in an airplane, this parameter should be a positive number indicating how far above the ground the observer is. An exception occurs if `metersAboveGround` is negative. | **Returns**: [`Time`](#Time) or `None` If the rise or set time is found within the specified time window, diff --git a/source/python/astronomy/astronomy.py b/source/python/astronomy/astronomy.py index ab3d5451..b413f9ee 100644 --- a/source/python/astronomy/astronomy.py +++ b/source/python/astronomy/astronomy.py @@ -4958,6 +4958,9 @@ def RefractionAngle(refraction: Refraction, altitude: float) -> float: the amount of "lift" caused by atmospheric refraction. This is the number of degrees higher in the sky an object appears due to lensing of the Earth's atmosphere. + This function works best near sea level. + To correct for higher elevations, call #Atmosphere for that + elevation and multiply the refraction angle by the resulting relative density. Parameters ---------- @@ -6470,7 +6473,34 @@ def _InternalSearchAltitude(body: Body, observer: Observer, direction: Direction a1 = a2 -def SearchRiseSet(body: Body, observer: Observer, direction: Direction, startTime: Time, limitDays: float) -> Optional[Time]: +def _HorizonDipAngle(observer: Observer, metersAboveGround: float) -> float: + # Calculate the effective radius of the Earth at ground level below the observer. + # Correct for the Earth's oblateness. + phi = math.radians(observer.latitude) + sinphi = math.sin(phi) + cosphi = math.cos(phi) + c = 1.0 / math.hypot(cosphi, sinphi*_EARTH_FLATTENING) + s = c * (_EARTH_FLATTENING * _EARTH_FLATTENING) + ht_km = (observer.height - metersAboveGround) / 1000.0 # height of ground above sea level + ach = _EARTH_EQUATORIAL_RADIUS_KM*c + ht_km + ash = _EARTH_EQUATORIAL_RADIUS_KM*s + ht_km + radius_m = 1000.0 * math.hypot(ach*cosphi, ash*sinphi) + + # Correct refraction of a ray of light traveling tangent to the Earth's surface. + # Based on: https://www.largeformatphotography.info/sunmooncalc/SMCalc.js + # which in turn derives from: + # Sweer, John. 1938. The Path of a Ray of Light Tangent to the Surface of the Earth. + # Journal of the Optical Society of America 28 (September):327-329. + + # k = refraction index + k = 0.175 * (1.0 - (6.5e-3/283.15)*(observer.height - (2.0/3.0)*metersAboveGround))**3.256 + + # Calculate how far below the observer's horizontal plane the observed horizon dips. + return math.degrees(-(math.sqrt(2*(1 - k)*metersAboveGround / radius_m) / (1 - k))) + + + +def SearchRiseSet(body: Body, observer: Observer, direction: Direction, startTime: Time, limitDays: float, metersAboveGround: float = 0.0) -> Optional[Time]: """Searches for the next time a celestial body rises or sets as seen by an observer on the Earth. This function finds the next rise or set time of the Sun, Moon, or planet other than the Earth. @@ -6515,6 +6545,13 @@ def SearchRiseSet(body: Body, observer: Observer, direction: Direction, startTim In cases where you want to find the next rise or set time no matter how far in the future (for example, for an observer near the south pole), you can pass in a larger value like 365. + metersAboveGround : float + Default value = 0.0. + Usually the observer is located at ground level. Then this parameter + should be zero. But if the observer is significantly higher than ground + level, for example in an airplane, this parameter should be a positive + number indicating how far above the ground the observer is. + An exception occurs if `metersAboveGround` is negative. Returns ------- @@ -6522,13 +6559,28 @@ def SearchRiseSet(body: Body, observer: Observer, direction: Direction, startTim If the rise or set time is found within the specified time window, this function returns that time. Otherwise, it returns `None`. """ + if not math.isfinite(metersAboveGround) or metersAboveGround < 0.0: + raise Error('Invalid value for metersAboveGround: {}'.format(metersAboveGround)) + + # Determine the radius of the body to be observed. if body == Body.Sun: bodyRadiusAu = _SUN_RADIUS_AU elif body == Body.Moon: bodyRadiusAu = _MOON_EQUATORIAL_RADIUS_AU else: bodyRadiusAu = 0.0 - return _InternalSearchAltitude(body, observer, direction, startTime, limitDays, bodyRadiusAu, -_REFRACTION_NEAR_HORIZON) + + # Calculate atmospheric density at ground level. + atmos = Atmosphere(observer.height - metersAboveGround) + + # Calculate the apparent angular dip of the horizon. + dip = _HorizonDipAngle(observer, metersAboveGround) + + # Correct refraction for objects near the horizon, using atmospheric density at the ground. + altitude = dip - (_REFRACTION_NEAR_HORIZON * atmos.density) + + # Search for the top of the body crossing the corrected altitude angle. + return _InternalSearchAltitude(body, observer, direction, startTime, limitDays, bodyRadiusAu, altitude) def SearchAltitude(body: Body, observer: Observer, direction: Direction, startTime: Time, limitDays: float, altitude: float) -> Optional[Time]: