mirror/astronomy
1
0
Fork 0
mirror of https://github.com/cosinekitty/astronomy.git synced 2026-04-07 16:17:25 -04:00
Code Issues Packages Projects Releases Wiki Activity

Simplified and optimized nutation formula.

Browse Source 
While trying to convert ecliptic coordinates from mean
equinox of date to true equinox of date, I ran into excessive
overhead from the IAU2000B nutation model. The fact that it
uses 77 trigonometric terms made the calculations a lot slower.

https://apps.dtic.mil/sti/pdfs/AD1112517.pdf
Page 4 in the above document mentions a shorter series
“NOD version 2” that has 13 terms instead of 77 as used in IAU2000B.
I had not noticed NOD2 before, because it appears only in
the FORTRAN version of NOVAS 3.x, not the C version.

After reading the FORTRAN code, I realized NOD2 is the same
as IAU2000B, only it keeps the first 13 of 77 terms.
The terms are already arranged in descending order of
significance, so it is easy to truncate the series.

Based on this discovery, I realized I could achieve all of
the required accuracy needed for Astronomy Engine by
keeping only the first 5 terms of the nutation series.
This tremendously speeds up nutation calculations while
sacrificing only a couple of arcseconds of accuracy.

It also makes the minified JavaScript code smaller:
Before: 119500 bytes.
After:  116653 bytes.

So that's what I did here. Most of the work was updating
unit tests for accepting slightly different calculation
results.

The nutation formula change did trigger detection of a
lurking bug in the inverse_terra functions, which convert
a geocentric vector into latitude, longitude, and elevation
(i.e. an Observer object). The Newton's Method loop in
this function was not always converging, resulting in
an infinite loop. I fixed that by increasing the
convergence threshold and throwing an exception
if the loop iterates more than 10 times.

I also fixed a couple of bugs in the `demotest` scripts.
This commit is contained in:
Don Cross
2022-12-04 10:31:15 -05:00
parent b2b426096f
commit 8a153315cf
57 changed files with 16550 additions and 18229 deletions
Download Patch File Download Diff File Expand all files Collapse all files

509
source/python/astronomy/astronomy.py
View File

@@ -1030,7 +1030,6 @@ class Spherical:
class _iau2000b:
def __init__(self, time):
t = time.tt / 36525.0
el = math.fmod((485868.249036 + t*1717915923.2178), _ASEC360) * _ASEC2RAD
elp = math.fmod((1287104.79305 + t*129596581.0481), _ASEC360) * _ASEC2RAD
f = math.fmod((335779.526232 + t*1739527262.8478), _ASEC360) * _ASEC2RAD
d = math.fmod((1072260.70369 + t*1602961601.2090), _ASEC360) * _ASEC2RAD
@@ -1072,508 +1071,6 @@ class _iau2000b:
de += (73871.0 - 184.0*t)*carg - 1924.0*sarg
arg = elp + 2.0*f - 2.0*d + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-516821.0 + 1226.0*t)*sarg - 524.0*carg
de += (224386.0 - 677.0*t)*carg - 174.0*sarg
sarg = math.sin(el)
carg = math.cos(el)
dp += (711159.0 + 73.0*t)*sarg - 872.0*carg
de += (-6750.0)*carg + 358.0*sarg
arg = 2.0*f + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-387298.0 - 367.0*t)*sarg + 380.0*carg
de += (200728.0 + 18.0*t)*carg + 318.0*sarg
arg = el + 2.0*f + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-301461.0 - 36.0*t)*sarg + 816.0*carg
de += (129025.0 - 63.0*t)*carg + 367.0*sarg
arg = -elp + 2.0*f - 2.0*d + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (215829.0 - 494.0*t)*sarg + 111.0*carg
de += (-95929.0 + 299.0*t)*carg + 132.0*sarg
arg = 2.0*f - 2.0*d + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (128227.0 + 137.0*t)*sarg + 181.0*carg
de += (-68982.0 - 9.0*t)*carg + 39.0*sarg
arg = -el + 2.0*f + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (123457.0 + 11.0*t)*sarg + 19.0*carg
de += (-53311.0 + 32.0*t)*carg - 4.0*sarg
arg = -el + 2.0*d
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (156994.0 + 10.0*t)*sarg - 168.0*carg
de += (-1235.0)*carg + 82.0*sarg
arg = el + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (63110.0 + 63.0*t)*sarg + 27.0*carg
de += (-33228.0)*carg - 9.0*sarg
arg = -el + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-57976.0 - 63.0*t)*sarg - 189.0*carg
de += (31429.0)*carg - 75.0*sarg
arg = -el + 2.0*f + 2.0*d + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-59641.0 - 11.0*t)*sarg + 149.0*carg
de += (25543.0 - 11.0*t)*carg + 66.0*sarg
arg = el + 2.0*f + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-51613.0 - 42.0*t)*sarg + 129.0*carg
de += (26366.0)*carg + 78.0*sarg
arg = -2.0*el + 2.0*f + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (45893.0 + 50.0*t)*sarg + 31.0*carg
de += (-24236.0 - 10.0*t)*carg + 20.0*sarg
arg = 2.0*d
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (63384.0 + 11.0*t)*sarg - 150.0*carg
de += (-1220.0)*carg + 29.0*sarg
arg = 2.0*f + 2.0*d + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-38571.0 - 1.0*t)*sarg + 158.0*carg
de += (16452.0 - 11.0*t)*carg + 68.0*sarg
arg = -2.0*elp + 2.0*f - 2.0*d + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (32481.0)*sarg
de += (-13870.0)*carg
arg = -2.0*el + 2.0*d
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-47722.0)*sarg - 18.0*carg
de += (477.0)*carg - 25.0*sarg
arg = 2.0*el + 2.0*f + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-31046.0 - 1.0*t)*sarg + 131.0*carg
de += (13238.0 - 11.0*t)*carg + 59.0*sarg
arg = el + 2.0*f - 2.0*d + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (28593.0)*sarg - carg
de += (-12338.0 + 10.0*t)*carg - 3.0*sarg
arg = -el + 2.0*f + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (20441.0 + 21.0*t)*sarg + 10.0*carg
de += (-10758.0)*carg - 3.0*sarg
arg = 2.0*el
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (29243.0)*sarg - 74.0*carg
de += (-609.0)*carg + 13.0*sarg
arg = 2.0*f
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (25887.0)*sarg - 66.0*carg
de += (-550.0)*carg + 11.0*sarg
arg = elp + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-14053.0 - 25.0*t)*sarg + 79.0*carg
de += (8551.0 - 2.0*t)*carg - 45.0*sarg
arg = -el + 2.0*d + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (15164.0 + 10.0*t)*sarg + 11.0*carg
de += (-8001.0)*carg - sarg
arg = 2.0*elp + 2.0*f - 2.0*d + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-15794.0 + 72.0*t)*sarg - 16.0*carg
de += (6850.0 - 42.0*t)*carg - 5.0*sarg
arg = -2.0*f + 2.0*d
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (21783.0)*sarg + 13.0*carg
de += (-167.0)*carg + 13.0*sarg
arg = el - 2.0*d + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-12873.0 - 10.0*t)*sarg - 37.0*carg
de += (6953.0)*carg - 14.0*sarg
arg = -elp + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-12654.0 + 11.0*t)*sarg + 63.0*carg
de += (6415.0)*carg + 26.0*sarg
arg = -el + 2.0*f + 2.0*d + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-10204.0)*sarg + 25.0*carg
de += (5222.0)*carg + 15.0*sarg
arg = 2.0*elp
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (16707.0 - 85.0*t)*sarg - 10.0*carg
de += (168.0 - 1.0*t)*carg + 10.0*sarg
arg = el + 2.0*f + 2.0*d + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-7691.0)*sarg + 44.0*carg
de += (3268.0)*carg + 19.0*sarg
arg = -2.0*el + 2.0*f
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-11024.0)*sarg - 14.0*carg
de += (104.0)*carg + 2.0*sarg
arg = elp + 2.0*f + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (7566.0 - 21.0*t)*sarg - 11.0*carg
de += (-3250.0)*carg - 5.0*sarg
arg = 2.0*f + 2.0*d + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-6637.0 - 11.0*t)*sarg + 25.0*carg
de += (3353.0)*carg + 14.0*sarg
arg = -elp + 2.0*f + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-7141.0 + 21.0*t)*sarg + 8.0*carg
de += (3070.0)*carg + 4.0*sarg
arg = 2.0*d + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-6302.0 - 11.0*t)*sarg + 2.0*carg
de += (3272.0)*carg + 4.0*sarg
arg = el + 2.0*f - 2.0*d + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (5800.0 + 10.0*t)*sarg + 2.0*carg
de += (-3045.0)*carg - sarg
arg = 2.0*el + 2.0*f - 2.0*d + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (6443.0)*sarg - 7.0*carg
de += (-2768.0)*carg - 4.0*sarg
arg = -2.0*el + 2.0*d + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-5774.0 - 11.0*t)*sarg - 15.0*carg
de += (3041.0)*carg - 5.0*sarg
arg = 2.0*el + 2.0*f + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-5350.0)*sarg + 21.0*carg
de += (2695.0)*carg + 12.0*sarg
arg = -elp + 2.0*f - 2.0*d + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-4752.0 - 11.0*t)*sarg - 3.0*carg
de += (2719.0)*carg - 3.0*sarg
arg = -2.0*d + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-4940.0 - 11.0*t)*sarg - 21.0*carg
de += (2720.0)*carg - 9.0*sarg
arg = -el - elp + 2.0*d
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (7350.0)*sarg - 8.0*carg
de += (-51.0)*carg + 4.0*sarg
arg = 2.0*el - 2.0*d + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (4065.0)*sarg + 6.0*carg
de += (-2206.0)*carg + sarg
arg = el + 2.0*d
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (6579.0)*sarg - 24.0*carg
de += (-199.0)*carg + 2.0*sarg
arg = elp + 2.0*f - 2.0*d + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (3579.0)*sarg + 5.0*carg
de += (-1900.0)*carg + sarg
arg = el - elp
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (4725.0)*sarg - 6.0*carg
de += (-41.0)*carg + 3.0*sarg
arg = -2.0*el + 2.0*f + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-3075.0)*sarg - 2.0*carg
de += (1313.0)*carg - sarg
arg = 3.0*el + 2.0*f + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-2904.0)*sarg + 15.0*carg
de += (1233.0)*carg + 7.0*sarg
arg = -elp + 2.0*d
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (4348.0)*sarg - 10.0*carg
de += (-81.0)*carg + 2.0*sarg
arg = el - elp + 2.0*f + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-2878.0)*sarg + 8.0*carg
de += (1232.0)*carg + 4.0*sarg
sarg = math.sin(d)
carg = math.cos(d)
dp += (-4230.0)*sarg + 5.0*carg
de += (-20.0)*carg - 2.0*sarg
arg = -el - elp + 2.0*f + 2.0*d + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-2819.0)*sarg + 7.0*carg
de += (1207.0)*carg + 3.0*sarg
arg = -el + 2.0*f
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-4056.0)*sarg + 5.0*carg
de += (40.0)*carg - 2.0*sarg
arg = -elp + 2.0*f + 2.0*d + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-2647.0)*sarg + 11.0*carg
de += (1129.0)*carg + 5.0*sarg
arg = -2.0*el + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-2294.0)*sarg - 10.0*carg
de += (1266.0)*carg - 4.0*sarg
arg = el + elp + 2.0*f + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (2481.0)*sarg - 7.0*carg
de += (-1062.0)*carg - 3.0*sarg
arg = 2.0*el + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (2179.0)*sarg - 2.0*carg
de += (-1129.0)*carg - 2.0*sarg
arg = -el + elp + d
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (3276.0)*sarg + carg
de += (-9.0)*carg
arg = el + elp
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-3389.0)*sarg + 5.0*carg
de += (35.0)*carg - 2.0*sarg
arg = el + 2.0*f
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (3339.0)*sarg - 13.0*carg
de += (-107.0)*carg + sarg
arg = -el + 2.0*f - 2.0*d + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-1987.0)*sarg - 6.0*carg
de += (1073.0)*carg - 2.0*sarg
arg = el + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-1981.0)*sarg
de += (854.0)*carg
arg = -el + d
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (4026.0)*sarg - 353.0*carg
de += (-553.0)*carg - 139.0*sarg
arg = 2.0*f + d + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (1660.0)*sarg - 5.0*carg
de += (-710.0)*carg - 2.0*sarg
arg = -el + 2.0*f + 4.0*d + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-1521.0)*sarg + 9.0*carg
de += (647.0)*carg + 4.0*sarg
arg = -el + elp + d + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (1314.0)*sarg
de += (-700.0)*carg
arg = -2.0*elp + 2.0*f - 2.0*d + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-1283.0)*sarg
de += (672.0)*carg
arg = el + 2.0*f + 2.0*d + om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (-1331.0)*sarg + 8.0*carg
de += (663.0)*carg + 4.0*sarg
arg = -2.0*el + 2.0*f + 2.0*d + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (1383.0)*sarg - 2.0*carg
de += (-594.0)*carg - 2.0*sarg
arg = -el + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (1405.0)*sarg + 4.0*carg
de += (-610.0)*carg + 2.0*sarg
arg = el + elp + 2.0*f - 2.0*d + 2.0*om
sarg = math.sin(arg)
carg = math.cos(arg)
dp += (1290.0)*sarg
de += (-556.0)*carg
self.dpsi = -0.000135 + (dp * 1.0e-7)
self.deps = +0.000388 + (de * 1.0e-7)
@@ -1870,7 +1367,11 @@ def _inverse_terra(ovec, st):
# Numerically solve for exact latitude, using Newton's Method.
# Start with initial latitude estimate, based on a spherical Earth.
lat = math.atan2(z, p)
count = 0
while True:
count += 1
if count > 10:
raise NoConvergeError()
# Calculate the error function W(lat).
# We try to find the root of W, meaning where the error is 0.
cos = math.cos(lat)
@@ -1881,7 +1382,7 @@ def _inverse_terra(ovec, st):
radicand = cos2 + _EARTH_FLATTENING_SQUARED*sin2
denom = math.sqrt(radicand)
W = (factor*sin*cos)/denom - z*cos + p*sin
if abs(W) < 1.0e-12:
if abs(W) < 1.0e-8:
# The error is now negligible
break
# Error is still too large. Find the next estimate.