Optimize for map-making calculation patterns.

See this discussion:
https://github.com/cosinekitty/astronomy/issues/150

For the case of calculating a map, where each pixel
on the map represents a different location on the Earth,
it is more efficient to factor out expensive calculation
of sidereal times, assuming the entire map represents
some phenomenon at a single moment in time.

For example, to determine whether the Moon is visible
at different places on the Earth, the following
functions can be calculated across thousands of
different (lat, lon) geographic coordinates around
the world:

    ObserverVector
    Rotation_EQD_HOR

Before iterating over the map pixels, a program
can call GeoMoon, then convert EQJ coordinates to EQD.

Then by passing the same time value in a loop to
ObserverVector and Rotation_EQD_HOR, the program
can calculate a vector from the observer to the Moon
in EQD coordinates, then convert EQD to HOR.
The z-coordinate of the horizontal coordinates
determines whether the Moon is above or below the
observer's horizon at that point on the Earth.

This calculation pattern performed redundant
sidereal time calculations for each pixel on the map.
I changed the code for all 4 languages to cache
sidereal time so that it only needs to be calculated
once.

In the C version of Astronomy Engine, this resulted
in a speedup factor of about 2.3 in the above use case.
(See the function MapPerformanceTest in generate/ctest.c.)
This commit is contained in:
Don Cross
2022-01-22 20:47:46 -05:00
parent 0bfdb359b1
commit 90a9839d18
21 changed files with 577 additions and 406 deletions

View File

@@ -560,7 +560,8 @@ class Time:
self.tt = _TerrestrialTime(ut)
else:
self.tt = tt
self.etilt = None
self._et = None # lazy-cache for earth tilt
self._st = None # lazy-cache for sidereal time
@staticmethod
def FromTerrestrialTime(tt):
@@ -717,9 +718,9 @@ class Time:
# Calculates precession and nutation of the Earth's axis.
# The calculations are very expensive, so lazy-evaluate and cache
# the result inside this Time object.
if self.etilt is None:
self.etilt = _e_tilt(self)
return self.etilt
if self._et is None:
self._et = _e_tilt(self)
return self._et
def __lt__(self, other):
return self.tt < other.tt
@@ -1565,20 +1566,22 @@ def _era(time): # Earth Rotation Angle
return theta
def _sidereal_time(time):
t = time.tt / 36525.0
eqeq = 15.0 * time._etilt().ee # Replace with eqeq=0 to get GMST instead of GAST (if we ever need it)
theta = _era(time)
st = (eqeq + 0.014506 +
(((( - 0.0000000368 * t
- 0.000029956 ) * t
- 0.00000044 ) * t
+ 1.3915817 ) * t
+ 4612.156534 ) * t)
gst = math.fmod((st/3600.0 + theta), 360.0) / 15.0
if gst < 0.0:
gst += 24.0
if time._st is None:
t = time.tt / 36525.0
eqeq = 15.0 * time._etilt().ee # Replace with eqeq=0 to get GMST instead of GAST (if we ever need it)
theta = _era(time)
st = (eqeq + 0.014506 +
(((( - 0.0000000368 * t
- 0.000029956 ) * t
- 0.00000044 ) * t
+ 1.3915817 ) * t
+ 4612.156534 ) * t)
gst = math.fmod((st/3600.0 + theta), 360.0) / 15.0
if gst < 0.0:
gst += 24.0
time._st = gst
# return sidereal hours in the half-open range [0, 24).
return gst
return time._st
def _inverse_terra(ovec, st):
# Convert from AU to kilometers