Skip to content

Commit c85881c

Browse files
cosinekittycosinekitty
committed
Python: added metersAboveGround parameter to SearchRiseSet.
1 parent 6c59d14 commit c85881c

File tree

5 files changed

+232
-7
lines changed

5 files changed

+232
-7
lines changed

Diff for: ‎demo/python/astronomy.py

+54-2
Original file line numberDiff line numberDiff line change
@@ -4958,6 +4958,9 @@ def RefractionAngle(refraction: Refraction, altitude: float) -> float:
49584958
the amount of "lift" caused by atmospheric refraction.
49594959
This is the number of degrees higher in the sky an object appears
49604960
due to lensing of the Earth's atmosphere.
4961+
This function works best near sea level.
4962+
To correct for higher elevations, call #Atmosphere for that
4963+
elevation and multiply the refraction angle by the resulting relative density.
49614964
49624965
Parameters
49634966
----------
@@ -6470,7 +6473,34 @@ def _InternalSearchAltitude(body: Body, observer: Observer, direction: Direction
64706473
a1 = a2
64716474

64726475

6473-
def SearchRiseSet(body: Body, observer: Observer, direction: Direction, startTime: Time, limitDays: float) -> Optional[Time]:
6476+
def _HorizonDipAngle(observer: Observer, metersAboveGround: float) -> float:
6477+
# Calculate the effective radius of the Earth at ground level below the observer.
6478+
# Correct for the Earth's oblateness.
6479+
phi = math.radians(observer.latitude)
6480+
sinphi = math.sin(phi)
6481+
cosphi = math.cos(phi)
6482+
c = 1.0 / math.hypot(cosphi, sinphi*_EARTH_FLATTENING)
6483+
s = c * (_EARTH_FLATTENING * _EARTH_FLATTENING)
6484+
ht_km = (observer.height - metersAboveGround) / 1000.0 # height of ground above sea level
6485+
ach = _EARTH_EQUATORIAL_RADIUS_KM*c + ht_km
6486+
ash = _EARTH_EQUATORIAL_RADIUS_KM*s + ht_km
6487+
radius_m = 1000.0 * math.hypot(ach*cosphi, ash*sinphi)
6488+
6489+
# Correct refraction of a ray of light traveling tangent to the Earth's surface.
6490+
# Based on: https://www.largeformatphotography.info/sunmooncalc/SMCalc.js
6491+
# which in turn derives from:
6492+
# Sweer, John. 1938. The Path of a Ray of Light Tangent to the Surface of the Earth.
6493+
# Journal of the Optical Society of America 28 (September):327-329.
6494+
6495+
# k = refraction index
6496+
k = 0.175 * (1.0 - (6.5e-3/283.15)*(observer.height - (2.0/3.0)*metersAboveGround))**3.256
6497+
6498+
# Calculate how far below the observer's horizontal plane the observed horizon dips.
6499+
return math.degrees(-(math.sqrt(2*(1 - k)*metersAboveGround / radius_m) / (1 - k)))
6500+
6501+
6502+
6503+
def SearchRiseSet(body: Body, observer: Observer, direction: Direction, startTime: Time, limitDays: float, metersAboveGround: float = 0.0) -> Optional[Time]:
64746504
"""Searches for the next time a celestial body rises or sets as seen by an observer on the Earth.
64756505
64766506
This function finds the next rise or set time of the Sun, Moon, or planet other than the Earth.
@@ -6515,20 +6545,42 @@ def SearchRiseSet(body: Body, observer: Observer, direction: Direction, startTim
65156545
In cases where you want to find the next rise or set time no matter how far
65166546
in the future (for example, for an observer near the south pole), you can
65176547
pass in a larger value like 365.
6548+
metersAboveGround : float
6549+
Default value = 0.0.
6550+
Usually the observer is located at ground level. Then this parameter
6551+
should be zero. But if the observer is significantly higher than ground
6552+
level, for example in an airplane, this parameter should be a positive
6553+
number indicating how far above the ground the observer is.
6554+
An exception occurs if `metersAboveGround` is negative.
65186555
65196556
Returns
65206557
-------
65216558
Time or `None`
65226559
If the rise or set time is found within the specified time window,
65236560
this function returns that time. Otherwise, it returns `None`.
65246561
"""
6562+
if not math.isfinite(metersAboveGround) or metersAboveGround < 0.0:
6563+
raise Error('Invalid value for metersAboveGround: {}'.format(metersAboveGround))
6564+
6565+
# Determine the radius of the body to be observed.
65256566
if body == Body.Sun:
65266567
bodyRadiusAu = _SUN_RADIUS_AU
65276568
elif body == Body.Moon:
65286569
bodyRadiusAu = _MOON_EQUATORIAL_RADIUS_AU
65296570
else:
65306571
bodyRadiusAu = 0.0
6531-
return _InternalSearchAltitude(body, observer, direction, startTime, limitDays, bodyRadiusAu, -_REFRACTION_NEAR_HORIZON)
6572+
6573+
# Calculate atmospheric density at ground level.
6574+
atmos = Atmosphere(observer.height - metersAboveGround)
6575+
6576+
# Calculate the apparent angular dip of the horizon.
6577+
dip = _HorizonDipAngle(observer, metersAboveGround)
6578+
6579+
# Correct refraction for objects near the horizon, using atmospheric density at the ground.
6580+
altitude = dip - (_REFRACTION_NEAR_HORIZON * atmos.density)
6581+
6582+
# Search for the top of the body crossing the corrected altitude angle.
6583+
return _InternalSearchAltitude(body, observer, direction, startTime, limitDays, bodyRadiusAu, altitude)
65326584

65336585

65346586
def SearchAltitude(body: Body, observer: Observer, direction: Direction, startTime: Time, limitDays: float, altitude: float) -> Optional[Time]:

Diff for: ‎generate/template/astronomy.py

+54-2
Original file line numberDiff line numberDiff line change
@@ -3452,6 +3452,9 @@ def RefractionAngle(refraction: Refraction, altitude: float) -> float:
34523452
the amount of "lift" caused by atmospheric refraction.
34533453
This is the number of degrees higher in the sky an object appears
34543454
due to lensing of the Earth's atmosphere.
3455+
This function works best near sea level.
3456+
To correct for higher elevations, call #Atmosphere for that
3457+
elevation and multiply the refraction angle by the resulting relative density.
34553458
34563459
Parameters
34573460
----------
@@ -4964,7 +4967,34 @@ def _InternalSearchAltitude(body: Body, observer: Observer, direction: Direction
49644967
a1 = a2
49654968

49664969

4967-
def SearchRiseSet(body: Body, observer: Observer, direction: Direction, startTime: Time, limitDays: float) -> Optional[Time]:
4970+
def _HorizonDipAngle(observer: Observer, metersAboveGround: float) -> float:
4971+
# Calculate the effective radius of the Earth at ground level below the observer.
4972+
# Correct for the Earth's oblateness.
4973+
phi = math.radians(observer.latitude)
4974+
sinphi = math.sin(phi)
4975+
cosphi = math.cos(phi)
4976+
c = 1.0 / math.hypot(cosphi, sinphi*_EARTH_FLATTENING)
4977+
s = c * (_EARTH_FLATTENING * _EARTH_FLATTENING)
4978+
ht_km = (observer.height - metersAboveGround) / 1000.0 # height of ground above sea level
4979+
ach = _EARTH_EQUATORIAL_RADIUS_KM*c + ht_km
4980+
ash = _EARTH_EQUATORIAL_RADIUS_KM*s + ht_km
4981+
radius_m = 1000.0 * math.hypot(ach*cosphi, ash*sinphi)
4982+
4983+
# Correct refraction of a ray of light traveling tangent to the Earth's surface.
4984+
# Based on: https://www.largeformatphotography.info/sunmooncalc/SMCalc.js
4985+
# which in turn derives from:
4986+
# Sweer, John. 1938. The Path of a Ray of Light Tangent to the Surface of the Earth.
4987+
# Journal of the Optical Society of America 28 (September):327-329.
4988+
4989+
# k = refraction index
4990+
k = 0.175 * (1.0 - (6.5e-3/283.15)*(observer.height - (2.0/3.0)*metersAboveGround))**3.256
4991+
4992+
# Calculate how far below the observer's horizontal plane the observed horizon dips.
4993+
return math.degrees(-(math.sqrt(2*(1 - k)*metersAboveGround / radius_m) / (1 - k)))
4994+
4995+
4996+
4997+
def SearchRiseSet(body: Body, observer: Observer, direction: Direction, startTime: Time, limitDays: float, metersAboveGround: float = 0.0) -> Optional[Time]:
49684998
"""Searches for the next time a celestial body rises or sets as seen by an observer on the Earth.
49694999
49705000
This function finds the next rise or set time of the Sun, Moon, or planet other than the Earth.
@@ -5009,20 +5039,42 @@ def SearchRiseSet(body: Body, observer: Observer, direction: Direction, startTim
50095039
In cases where you want to find the next rise or set time no matter how far
50105040
in the future (for example, for an observer near the south pole), you can
50115041
pass in a larger value like 365.
5042+
metersAboveGround : float
5043+
Default value = 0.0.
5044+
Usually the observer is located at ground level. Then this parameter
5045+
should be zero. But if the observer is significantly higher than ground
5046+
level, for example in an airplane, this parameter should be a positive
5047+
number indicating how far above the ground the observer is.
5048+
An exception occurs if `metersAboveGround` is negative.
50125049
50135050
Returns
50145051
-------
50155052
Time or `None`
50165053
If the rise or set time is found within the specified time window,
50175054
this function returns that time. Otherwise, it returns `None`.
50185055
"""
5056+
if not math.isfinite(metersAboveGround) or metersAboveGround < 0.0:
5057+
raise Error('Invalid value for metersAboveGround: {}'.format(metersAboveGround))
5058+
5059+
# Determine the radius of the body to be observed.
50195060
if body == Body.Sun:
50205061
bodyRadiusAu = _SUN_RADIUS_AU
50215062
elif body == Body.Moon:
50225063
bodyRadiusAu = _MOON_EQUATORIAL_RADIUS_AU
50235064
else:
50245065
bodyRadiusAu = 0.0
5025-
return _InternalSearchAltitude(body, observer, direction, startTime, limitDays, bodyRadiusAu, -_REFRACTION_NEAR_HORIZON)
5066+
5067+
# Calculate atmospheric density at ground level.
5068+
atmos = Atmosphere(observer.height - metersAboveGround)
5069+
5070+
# Calculate the apparent angular dip of the horizon.
5071+
dip = _HorizonDipAngle(observer, metersAboveGround)
5072+
5073+
# Correct refraction for objects near the horizon, using atmospheric density at the ground.
5074+
altitude = dip - (_REFRACTION_NEAR_HORIZON * atmos.density)
5075+
5076+
# Search for the top of the body crossing the corrected altitude angle.
5077+
return _InternalSearchAltitude(body, observer, direction, startTime, limitDays, bodyRadiusAu, altitude)
50265078

50275079

50285080
def SearchAltitude(body: Body, observer: Observer, direction: Direction, startTime: Time, limitDays: float, altitude: float) -> Optional[Time]:

Diff for: ‎generate/test.py

+65
Original file line numberDiff line numberDiff line change
@@ -3264,6 +3264,70 @@ def Atmosphere():
32643264

32653265
#-----------------------------------------------------------------------------------------------------------
32663266

3267+
def RiseSetElevationBodyCase(body, observer, direction, metersAboveGround, startTime, eventOffsetDays):
3268+
time = astronomy.SearchRiseSet(body, observer, direction, startTime, 2.0, metersAboveGround)
3269+
if not time:
3270+
return Fail('RiseSetElevationBodyCase {} {}: search failed.'.format(body, direction))
3271+
diff = v(time.ut - (startTime.ut + eventOffsetDays))
3272+
if diff > 0.5:
3273+
diff -= 1.0 # assume event actually takes place on the next day
3274+
diff = vabs(MINUTES_PER_DAY * diff) # convert signed days to absolute minutes
3275+
if diff > 0.5:
3276+
return Fail('RiseSetElevationBodyCase {} {}: EXCESSIVE diff = {}.'.format(body, direction, diff))
3277+
return 0
3278+
3279+
3280+
def RiseSetElevation():
3281+
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*$')
3282+
filename = 'riseset/elevation.txt'
3283+
with open(filename, 'rt') as infile:
3284+
lnum = 0
3285+
for line in infile:
3286+
lnum += 1
3287+
if line.startswith('#'):
3288+
continue
3289+
m = regex.match(line)
3290+
if not m:
3291+
return Fail('RiseSetElevation({} line {})'.format(filename, lnum), 'Invalid data format')
3292+
year = int(m.group(1))
3293+
month = int(m.group(2))
3294+
day = int(m.group(3))
3295+
latitude = v(float(m.group(4)))
3296+
longitude = v(float(m.group(5)))
3297+
height = v(float(m.group(6)))
3298+
metersAboveGround = v(float(m.group(7)))
3299+
srh = int(m.group( 8))
3300+
srm = int(m.group( 9))
3301+
ssh = int(m.group(10))
3302+
ssm = int(m.group(11))
3303+
mrh = int(m.group(12))
3304+
mrm = int(m.group(13))
3305+
msh = int(m.group(14))
3306+
msm = int(m.group(15))
3307+
3308+
# Get search origin time
3309+
time = astronomy.Time.Make(year, month, day, 0, 0, 0.0)
3310+
3311+
# Convert scanned values into sunrise, sunset, moonrise, moonset day offsets.
3312+
sr = (srh + srm/60.0) / 24.0
3313+
ss = (ssh + ssm/60.0) / 24.0
3314+
mr = (mrh + mrm/60.0) / 24.0
3315+
ms = (msh + msm/60.0) / 24.0
3316+
3317+
observer = astronomy.Observer(latitude, longitude, height)
3318+
3319+
if (0 != RiseSetElevationBodyCase(astronomy.Body.Sun, observer, astronomy.Direction.Rise, metersAboveGround, time, sr) or
3320+
0 != RiseSetElevationBodyCase(astronomy.Body.Sun, observer, astronomy.Direction.Set, metersAboveGround, time, ss) or
3321+
0 != RiseSetElevationBodyCase(astronomy.Body.Moon, observer, astronomy.Direction.Rise, metersAboveGround, time, mr) or
3322+
0 != RiseSetElevationBodyCase(astronomy.Body.Moon, observer, astronomy.Direction.Set, metersAboveGround, time, ms)):
3323+
return 1
3324+
3325+
3326+
return Pass('RiseSetElevation')
3327+
3328+
#-----------------------------------------------------------------------------------------------------------
3329+
3330+
32673331
UnitTests = {
32683332
'aberration': Aberration,
32693333
'atmosphere': Atmosphere,
@@ -3297,6 +3361,7 @@ def Atmosphere():
32973361
'refraction': Refraction,
32983362
'repr': Repr,
32993363
'riseset': RiseSet,
3364+
'riseset_elevation': RiseSetElevation,
33003365
'riseset_reverse': RiseSetReverse,
33013366
'rotation': Rotation,
33023367
'seasons': Seasons,

Diff for: ‎source/python/README.md

+5-1
Original file line numberDiff line numberDiff line change
@@ -2494,6 +2494,9 @@ Given an altitude angle and a refraction option, calculates
24942494
the amount of "lift" caused by atmospheric refraction.
24952495
This is the number of degrees higher in the sky an object appears
24962496
due to lensing of the Earth's atmosphere.
2497+
This function works best near sea level.
2498+
To correct for higher elevations, call [`Atmosphere`](#Atmosphere) for that
2499+
elevation and multiply the refraction angle by the resulting relative density.
24972500

24982501
| Type | Parameter | Description |
24992502
| --- | --- | --- |
@@ -3324,7 +3327,7 @@ The date and time of the relative longitude event.
33243327
---
33253328

33263329
<a name="SearchRiseSet"></a>
3327-
### SearchRiseSet(body: [`Body`](#Body), observer: [`Observer`](#Observer), direction: [`Direction`](#Direction), startTime: [`Time`](#Time), limitDays: float) &#8594; Optional\[[`Time`](#Time)\]
3330+
### SearchRiseSet(body: [`Body`](#Body), observer: [`Observer`](#Observer), direction: [`Direction`](#Direction), startTime: [`Time`](#Time), limitDays: float, metersAboveGround: float = 0.0) &#8594; Optional\[[`Time`](#Time)\]
33283331

33293332
**Searches for the next time a celestial body rises or sets as seen by an observer on the Earth.**
33303333

@@ -3355,6 +3358,7 @@ Therefore callers must not assume that the function will always succeed.
33553358
| [`Direction`](#Direction) | `direction` | Either `Direction.Rise` to find a rise time or `Direction.Set` to find a set time. |
33563359
| [`Time`](#Time) | `startTime` | The date and time at which to start the search. |
33573360
| `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. |
3361+
| `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. |
33583362

33593363
**Returns**: [`Time`](#Time) or `None`
33603364
If the rise or set time is found within the specified time window,

Diff for: ‎source/python/astronomy/astronomy.py

+54-2
Original file line numberDiff line numberDiff line change
@@ -4958,6 +4958,9 @@ def RefractionAngle(refraction: Refraction, altitude: float) -> float:
49584958
the amount of "lift" caused by atmospheric refraction.
49594959
This is the number of degrees higher in the sky an object appears
49604960
due to lensing of the Earth's atmosphere.
4961+
This function works best near sea level.
4962+
To correct for higher elevations, call #Atmosphere for that
4963+
elevation and multiply the refraction angle by the resulting relative density.
49614964
49624965
Parameters
49634966
----------
@@ -6470,7 +6473,34 @@ def _InternalSearchAltitude(body: Body, observer: Observer, direction: Direction
64706473
a1 = a2
64716474

64726475

6473-
def SearchRiseSet(body: Body, observer: Observer, direction: Direction, startTime: Time, limitDays: float) -> Optional[Time]:
6476+
def _HorizonDipAngle(observer: Observer, metersAboveGround: float) -> float:
6477+
# Calculate the effective radius of the Earth at ground level below the observer.
6478+
# Correct for the Earth's oblateness.
6479+
phi = math.radians(observer.latitude)
6480+
sinphi = math.sin(phi)
6481+
cosphi = math.cos(phi)
6482+
c = 1.0 / math.hypot(cosphi, sinphi*_EARTH_FLATTENING)
6483+
s = c * (_EARTH_FLATTENING * _EARTH_FLATTENING)
6484+
ht_km = (observer.height - metersAboveGround) / 1000.0 # height of ground above sea level
6485+
ach = _EARTH_EQUATORIAL_RADIUS_KM*c + ht_km
6486+
ash = _EARTH_EQUATORIAL_RADIUS_KM*s + ht_km
6487+
radius_m = 1000.0 * math.hypot(ach*cosphi, ash*sinphi)
6488+
6489+
# Correct refraction of a ray of light traveling tangent to the Earth's surface.
6490+
# Based on: https://www.largeformatphotography.info/sunmooncalc/SMCalc.js
6491+
# which in turn derives from:
6492+
# Sweer, John. 1938. The Path of a Ray of Light Tangent to the Surface of the Earth.
6493+
# Journal of the Optical Society of America 28 (September):327-329.
6494+
6495+
# k = refraction index
6496+
k = 0.175 * (1.0 - (6.5e-3/283.15)*(observer.height - (2.0/3.0)*metersAboveGround))**3.256
6497+
6498+
# Calculate how far below the observer's horizontal plane the observed horizon dips.
6499+
return math.degrees(-(math.sqrt(2*(1 - k)*metersAboveGround / radius_m) / (1 - k)))
6500+
6501+
6502+
6503+
def SearchRiseSet(body: Body, observer: Observer, direction: Direction, startTime: Time, limitDays: float, metersAboveGround: float = 0.0) -> Optional[Time]:
64746504
"""Searches for the next time a celestial body rises or sets as seen by an observer on the Earth.
64756505
64766506
This function finds the next rise or set time of the Sun, Moon, or planet other than the Earth.
@@ -6515,20 +6545,42 @@ def SearchRiseSet(body: Body, observer: Observer, direction: Direction, startTim
65156545
In cases where you want to find the next rise or set time no matter how far
65166546
in the future (for example, for an observer near the south pole), you can
65176547
pass in a larger value like 365.
6548+
metersAboveGround : float
6549+
Default value = 0.0.
6550+
Usually the observer is located at ground level. Then this parameter
6551+
should be zero. But if the observer is significantly higher than ground
6552+
level, for example in an airplane, this parameter should be a positive
6553+
number indicating how far above the ground the observer is.
6554+
An exception occurs if `metersAboveGround` is negative.
65186555
65196556
Returns
65206557
-------
65216558
Time or `None`
65226559
If the rise or set time is found within the specified time window,
65236560
this function returns that time. Otherwise, it returns `None`.
65246561
"""
6562+
if not math.isfinite(metersAboveGround) or metersAboveGround < 0.0:
6563+
raise Error('Invalid value for metersAboveGround: {}'.format(metersAboveGround))
6564+
6565+
# Determine the radius of the body to be observed.
65256566
if body == Body.Sun:
65266567
bodyRadiusAu = _SUN_RADIUS_AU
65276568
elif body == Body.Moon:
65286569
bodyRadiusAu = _MOON_EQUATORIAL_RADIUS_AU
65296570
else:
65306571
bodyRadiusAu = 0.0
6531-
return _InternalSearchAltitude(body, observer, direction, startTime, limitDays, bodyRadiusAu, -_REFRACTION_NEAR_HORIZON)
6572+
6573+
# Calculate atmospheric density at ground level.
6574+
atmos = Atmosphere(observer.height - metersAboveGround)
6575+
6576+
# Calculate the apparent angular dip of the horizon.
6577+
dip = _HorizonDipAngle(observer, metersAboveGround)
6578+
6579+
# Correct refraction for objects near the horizon, using atmospheric density at the ground.
6580+
altitude = dip - (_REFRACTION_NEAR_HORIZON * atmos.density)
6581+
6582+
# Search for the top of the body crossing the corrected altitude angle.
6583+
return _InternalSearchAltitude(body, observer, direction, startTime, limitDays, bodyRadiusAu, altitude)
65326584

65336585

65346586
def SearchAltitude(body: Body, observer: Observer, direction: Direction, startTime: Time, limitDays: float, altitude: float) -> Optional[Time]:

0 commit comments

Comments
 (0)