From 7c475fcadaffdbec81e8294f2733e4f38d0cda44 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Mon, 27 May 2024 17:07:30 -0400 Subject: [PATCH] Expanded the fix for issue #347. I tried more distant objects like Jupiter ... Neptune. This revealed that at increasing distances, the convergence threshold in inverse_terra needed to increased also. So now I use 1 AU as a baseline, and scale up linearly for more distant objects. --- README.md | 2 +- demo/browser/astronomy.browser.js | 11 +++-- demo/nodejs/astronomy.js | 11 +++-- demo/nodejs/calendar/astronomy.ts | 11 +++-- demo/python/astronomy.py | 3 +- generate/template/astronomy.c | 7 +++- generate/template/astronomy.cs | 3 +- generate/template/astronomy.kt | 3 +- generate/template/astronomy.py | 3 +- generate/template/astronomy.ts | 11 +++-- generate/test.js | 30 +++++++++----- source/c/astronomy.c | 7 +++- source/csharp/astronomy.cs | 3 +- source/js/astronomy.browser.js | 11 +++-- source/js/astronomy.browser.min.js | 38 +++++++++--------- source/js/astronomy.js | 11 +++-- source/js/astronomy.min.js | 40 +++++++++---------- source/js/astronomy.ts | 11 +++-- source/js/esm/astronomy.js | 11 +++-- .../github/cosinekitty/astronomy/astronomy.kt | 3 +- source/python/astronomy/astronomy.py | 3 +- 21 files changed, 140 insertions(+), 93 deletions(-) diff --git a/README.md b/README.md index b7c5b7b8..7e9b8410 100644 --- a/README.md +++ b/README.md @@ -189,7 +189,7 @@ of complexity. So I decided to create Astronomy Engine with the following engine - Support JavaScript, C, C#, and Python with the same algorithms, and verify them to produce identical results. (Kotlin support was added in 2022.) - No external dependencies! The code must not require anything outside the standard library for each language. -- Minified JavaScript code less than 120K. (The current size is 116424 bytes.) +- Minified JavaScript code less than 120K. (The current size is 116485 bytes.) - Accuracy always within 1 arcminute of results from NOVAS. - It would be well documented, relatively easy to use, and support a wide variety of common use cases. diff --git a/demo/browser/astronomy.browser.js b/demo/browser/astronomy.browser.js index b9ac051e..3fb1f874 100644 --- a/demo/browser/astronomy.browser.js +++ b/demo/browser/astronomy.browser.js @@ -1938,9 +1938,12 @@ function inverse_terra(ovec, st) { let lat = Math.atan2(z, p); let cos, sin, denom; let count = 0; + let W = 0; + const distanceAu = Math.max(1, Math.hypot(ovec[0], ovec[1], ovec[2])); for (;;) { - if (++count > 10) - throw `inverse_terra failed to converge.`; + if (++count > 10) { + throw `inverse_terra failed to converge: W=${W}, distanceAu=${distanceAu}`; + } // Calculate the error function W(lat). // We try to find the root of W, meaning where the error is 0. cos = Math.cos(lat); @@ -1950,8 +1953,8 @@ function inverse_terra(ovec, st) { const sin2 = sin * sin; const radicand = cos2 + EARTH_FLATTENING_SQUARED * sin2; denom = Math.sqrt(radicand); - const W = (factor * sin * cos) / denom - z * cos + p * sin; - if (Math.abs(W) < 2.0e-8) + W = (factor * sin * cos) / denom - z * cos + p * sin; + if (Math.abs(W) < distanceAu * 2.0e-8) break; // The error is now negligible // Error is still too large. Find the next estimate. // Calculate D = the derivative of W with respect to lat. diff --git a/demo/nodejs/astronomy.js b/demo/nodejs/astronomy.js index 94dd87ad..1e36c36b 100644 --- a/demo/nodejs/astronomy.js +++ b/demo/nodejs/astronomy.js @@ -1937,9 +1937,12 @@ function inverse_terra(ovec, st) { let lat = Math.atan2(z, p); let cos, sin, denom; let count = 0; + let W = 0; + const distanceAu = Math.max(1, Math.hypot(ovec[0], ovec[1], ovec[2])); for (;;) { - if (++count > 10) - throw `inverse_terra failed to converge.`; + if (++count > 10) { + throw `inverse_terra failed to converge: W=${W}, distanceAu=${distanceAu}`; + } // Calculate the error function W(lat). // We try to find the root of W, meaning where the error is 0. cos = Math.cos(lat); @@ -1949,8 +1952,8 @@ function inverse_terra(ovec, st) { const sin2 = sin * sin; const radicand = cos2 + EARTH_FLATTENING_SQUARED * sin2; denom = Math.sqrt(radicand); - const W = (factor * sin * cos) / denom - z * cos + p * sin; - if (Math.abs(W) < 2.0e-8) + W = (factor * sin * cos) / denom - z * cos + p * sin; + if (Math.abs(W) < distanceAu * 2.0e-8) break; // The error is now negligible // Error is still too large. Find the next estimate. // Calculate D = the derivative of W with respect to lat. diff --git a/demo/nodejs/calendar/astronomy.ts b/demo/nodejs/calendar/astronomy.ts index 52c0970d..6a1b8156 100644 --- a/demo/nodejs/calendar/astronomy.ts +++ b/demo/nodejs/calendar/astronomy.ts @@ -2108,9 +2108,12 @@ function inverse_terra(ovec: ArrayVector, st: number): Observer { let lat = Math.atan2(z, p); let cos:number, sin:number, denom:number; let count = 0; + let W = 0; + const distanceAu = Math.max(1, Math.hypot(ovec[0], ovec[1], ovec[2])); for(;;) { - if (++count > 10) - throw `inverse_terra failed to converge.`; + if (++count > 10) { + throw `inverse_terra failed to converge: W=${W}, distanceAu=${distanceAu}`; + } // Calculate the error function W(lat). // We try to find the root of W, meaning where the error is 0. cos = Math.cos(lat); @@ -2120,8 +2123,8 @@ function inverse_terra(ovec: ArrayVector, st: number): Observer { const sin2 = sin*sin; const radicand = cos2 + EARTH_FLATTENING_SQUARED*sin2; denom = Math.sqrt(radicand); - const W = (factor*sin*cos)/denom - z*cos + p*sin; - if (Math.abs(W) < 2.0e-8) + W = (factor*sin*cos)/denom - z*cos + p*sin; + if (Math.abs(W) < distanceAu * 2.0e-8) break; // The error is now negligible // Error is still too large. Find the next estimate. // Calculate D = the derivative of W with respect to lat. diff --git a/demo/python/astronomy.py b/demo/python/astronomy.py index d20ae486..e4b922eb 100644 --- a/demo/python/astronomy.py +++ b/demo/python/astronomy.py @@ -1435,6 +1435,7 @@ def _inverse_terra(ovec: List[float], st: float) -> Observer: # Start with initial latitude estimate, based on a spherical Earth. lat = math.atan2(z, p) count = 0 + distanceAu = max(1.0, math.hypot(ovec[0], ovec[1], ovec[2])) while True: count += 1 if count > 10: @@ -1449,7 +1450,7 @@ def _inverse_terra(ovec: List[float], st: float) -> Observer: radicand = cos2 + _EARTH_FLATTENING_SQUARED*sin2 denom = math.sqrt(radicand) W = (factor*sin*cos)/denom - z*cos + p*sin - if abs(W) < 2.0e-8: + if abs(W) < distanceAu * 2.0e-8: # The error is now negligible break # Error is still too large. Find the next estimate. diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index 0fe7bc3e..616d924e 100644 --- a/generate/template/astronomy.c +++ b/generate/template/astronomy.c @@ -1809,7 +1809,7 @@ static astro_observer_t inverse_terra(const double ovec[3], double st) { double x, y, z, p, F, W, D, c, s, c2, s2; double lon_deg, lat_deg, lat, radicand, factor, denom, adjust; - double height_km, stlocl; + double height_km, stlocl, distance_au; astro_observer_t observer; int count; @@ -1841,6 +1841,9 @@ static astro_observer_t inverse_terra(const double ovec[3], double st) F = EARTH_FLATTENING * EARTH_FLATTENING; /* Start with initial latitude estimate, based on a spherical Earth. */ lat = atan2(z, p); + distance_au = sqrt(ovec[0]*ovec[0] + ovec[1]*ovec[1] + ovec[2]*ovec[2]); + if (distance_au < 1.0) + distance_au = 1.0; for (count = 0; ; ++count) { if (count > 10) @@ -1859,7 +1862,7 @@ static astro_observer_t inverse_terra(const double ovec[3], double st) radicand = c2 + F*s2; denom = sqrt(radicand); W = (factor*s*c)/denom - z*c + p*s; - if (fabs(W) < 2.0e-8) + if (fabs(W) < distance_au * 2.0e-8) break; /* The error is now negligible. */ /* Error is still too large. Find the next estimate. */ /* Calculate D = the derivative of W with respect to lat. */ diff --git a/generate/template/astronomy.cs b/generate/template/astronomy.cs index 667fc395..9917e279 100644 --- a/generate/template/astronomy.cs +++ b/generate/template/astronomy.cs @@ -4162,6 +4162,7 @@ namespace CosineKitty double lat = Math.Atan2(z, p); double c, s, denom; int count = 0; + double distanceAu = Math.Max(1.0, ovec.Length()); for(;;) { if (++count > 10) @@ -4176,7 +4177,7 @@ namespace CosineKitty double radicand = c2 + F*s2; denom = Math.Sqrt(radicand); double W = (factor*s*c)/denom - z*c + p*s; - if (Math.Abs(W) < 2.0e-8) + if (Math.Abs(W) < distanceAu * 2.0e-8) break; // The error is now negligible. // Error is still too large. Find the next estimate. // Calculate D = the derivative of W with respect to lat. diff --git a/generate/template/astronomy.kt b/generate/template/astronomy.kt index fbebfc66..b4b00e35 100644 --- a/generate/template/astronomy.kt +++ b/generate/template/astronomy.kt @@ -4370,6 +4370,7 @@ private fun inverseTerra(ovec: Vector): Observer { var s: Double var denom: Double var count = 0 + val distanceAu = max(1.0, ovec.length()) while (true) { ++count if (count > 10) @@ -4384,7 +4385,7 @@ private fun inverseTerra(ovec: Vector): Observer { val radicand = c2 + F*s2 denom = sqrt(radicand) val W = ((factor * s * c) / denom) - (z * c) + (p * s) - if (W.absoluteValue < 2.0e-8) + if (W.absoluteValue < distanceAu * 2.0e-8) break // The error is now negligible. // Error is still too large. Find the next estimate. // Calculate D = the derivative of W with respect to lat. diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index eccf13be..12040796 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -1435,6 +1435,7 @@ def _inverse_terra(ovec: List[float], st: float) -> Observer: # Start with initial latitude estimate, based on a spherical Earth. lat = math.atan2(z, p) count = 0 + distanceAu = max(1.0, math.hypot(ovec[0], ovec[1], ovec[2])) while True: count += 1 if count > 10: @@ -1449,7 +1450,7 @@ def _inverse_terra(ovec: List[float], st: float) -> Observer: radicand = cos2 + _EARTH_FLATTENING_SQUARED*sin2 denom = math.sqrt(radicand) W = (factor*sin*cos)/denom - z*cos + p*sin - if abs(W) < 2.0e-8: + if abs(W) < distanceAu * 2.0e-8: # The error is now negligible break # Error is still too large. Find the next estimate. diff --git a/generate/template/astronomy.ts b/generate/template/astronomy.ts index 2ce04931..7b5e2eb1 100644 --- a/generate/template/astronomy.ts +++ b/generate/template/astronomy.ts @@ -1375,9 +1375,12 @@ function inverse_terra(ovec: ArrayVector, st: number): Observer { let lat = Math.atan2(z, p); let cos:number, sin:number, denom:number; let count = 0; + let W = 0; + const distanceAu = Math.max(1, Math.hypot(ovec[0], ovec[1], ovec[2])); for(;;) { - if (++count > 10) - throw `inverse_terra failed to converge.`; + if (++count > 10) { + throw `inverse_terra failed to converge: W=${W}, distanceAu=${distanceAu}`; + } // Calculate the error function W(lat). // We try to find the root of W, meaning where the error is 0. cos = Math.cos(lat); @@ -1387,8 +1390,8 @@ function inverse_terra(ovec: ArrayVector, st: number): Observer { const sin2 = sin*sin; const radicand = cos2 + EARTH_FLATTENING_SQUARED*sin2; denom = Math.sqrt(radicand); - const W = (factor*sin*cos)/denom - z*cos + p*sin; - if (Math.abs(W) < 2.0e-8) + W = (factor*sin*cos)/denom - z*cos + p*sin; + if (Math.abs(W) < distanceAu * 2.0e-8) break; // The error is now negligible // Error is still too large. Find the next estimate. // Calculate D = the derivative of W with respect to lat. diff --git a/generate/test.js b/generate/test.js index 2d18f345..ba48b700 100644 --- a/generate/test.js +++ b/generate/test.js @@ -2847,18 +2847,28 @@ function TwilightTest() { } +function VectorObserverCase(body) { + for (let ts = 1717780096000; ts <= 1717780096200; ts++) { + var date = new Date(ts); + var vect = Astronomy.GeoVector(body, date, true) + //console.log('ts, date, vector', ts, date.valueOf(), vect) + var point = Astronomy.VectorObserver(vect) + //console.log('point', point) + //console.log('*') + } + return Pass(`VectorObserverCase(${body})`); + } + + function VectorObserverIssue347() { // https://github.com/cosinekitty/astronomy/issues/347 - const body = Astronomy.Body.Sun - for (let ts = 1717780096000; ts <= 1717780096200; ts++) { - var date = new Date(ts); - var vect = Astronomy.GeoVector(body, date, true) - //console.log('ts, date, vector', ts, date.valueOf(), vect) - var point = Astronomy.VectorObserver(vect) - //console.log('point', point) - //console.log('*') - } - return Pass('VectorObserverIssue347'); + return ( + VectorObserverCase(Astronomy.Body.Sun) || + VectorObserverCase(Astronomy.Body.Jupiter) || + VectorObserverCase(Astronomy.Body.Saturn) || + VectorObserverCase(Astronomy.Body.Uranus) || + VectorObserverCase(Astronomy.Body.Neptune) + ); } diff --git a/source/c/astronomy.c b/source/c/astronomy.c index 291534a8..bb9119d6 100644 --- a/source/c/astronomy.c +++ b/source/c/astronomy.c @@ -1815,7 +1815,7 @@ static astro_observer_t inverse_terra(const double ovec[3], double st) { double x, y, z, p, F, W, D, c, s, c2, s2; double lon_deg, lat_deg, lat, radicand, factor, denom, adjust; - double height_km, stlocl; + double height_km, stlocl, distance_au; astro_observer_t observer; int count; @@ -1847,6 +1847,9 @@ static astro_observer_t inverse_terra(const double ovec[3], double st) F = EARTH_FLATTENING * EARTH_FLATTENING; /* Start with initial latitude estimate, based on a spherical Earth. */ lat = atan2(z, p); + distance_au = sqrt(ovec[0]*ovec[0] + ovec[1]*ovec[1] + ovec[2]*ovec[2]); + if (distance_au < 1.0) + distance_au = 1.0; for (count = 0; ; ++count) { if (count > 10) @@ -1865,7 +1868,7 @@ static astro_observer_t inverse_terra(const double ovec[3], double st) radicand = c2 + F*s2; denom = sqrt(radicand); W = (factor*s*c)/denom - z*c + p*s; - if (fabs(W) < 2.0e-8) + if (fabs(W) < distance_au * 2.0e-8) break; /* The error is now negligible. */ /* Error is still too large. Find the next estimate. */ /* Calculate D = the derivative of W with respect to lat. */ diff --git a/source/csharp/astronomy.cs b/source/csharp/astronomy.cs index b886cf28..f01a2a49 100644 --- a/source/csharp/astronomy.cs +++ b/source/csharp/astronomy.cs @@ -5296,6 +5296,7 @@ namespace CosineKitty double lat = Math.Atan2(z, p); double c, s, denom; int count = 0; + double distanceAu = Math.Max(1.0, ovec.Length()); for(;;) { if (++count > 10) @@ -5310,7 +5311,7 @@ namespace CosineKitty double radicand = c2 + F*s2; denom = Math.Sqrt(radicand); double W = (factor*s*c)/denom - z*c + p*s; - if (Math.Abs(W) < 2.0e-8) + if (Math.Abs(W) < distanceAu * 2.0e-8) break; // The error is now negligible. // Error is still too large. Find the next estimate. // Calculate D = the derivative of W with respect to lat. diff --git a/source/js/astronomy.browser.js b/source/js/astronomy.browser.js index b9ac051e..3fb1f874 100644 --- a/source/js/astronomy.browser.js +++ b/source/js/astronomy.browser.js @@ -1938,9 +1938,12 @@ function inverse_terra(ovec, st) { let lat = Math.atan2(z, p); let cos, sin, denom; let count = 0; + let W = 0; + const distanceAu = Math.max(1, Math.hypot(ovec[0], ovec[1], ovec[2])); for (;;) { - if (++count > 10) - throw `inverse_terra failed to converge.`; + if (++count > 10) { + throw `inverse_terra failed to converge: W=${W}, distanceAu=${distanceAu}`; + } // Calculate the error function W(lat). // We try to find the root of W, meaning where the error is 0. cos = Math.cos(lat); @@ -1950,8 +1953,8 @@ function inverse_terra(ovec, st) { const sin2 = sin * sin; const radicand = cos2 + EARTH_FLATTENING_SQUARED * sin2; denom = Math.sqrt(radicand); - const W = (factor * sin * cos) / denom - z * cos + p * sin; - if (Math.abs(W) < 2.0e-8) + W = (factor * sin * cos) / denom - z * cos + p * sin; + if (Math.abs(W) < distanceAu * 2.0e-8) break; // The error is now negligible // Error is still too large. Find the next estimate. // Calculate D = the derivative of W with respect to lat. diff --git a/source/js/astronomy.browser.min.js b/source/js/astronomy.browser.min.js index e86dcdfa..eda17747 100644 --- a/source/js/astronomy.browser.min.js +++ b/source/js/astronomy.browser.min.js @@ -175,25 +175,25 @@ k));l=Math.asin(-Math.sin(t)*Math.cos(h)*Math.sin(l)-Math.sin(h)*Math.cos(l));t= (72.56+20.186*a))+-1.4E-4*Math.sin(2*d)+1.4E-4*Math.cos(2*q-2*k)+-1.2E-4*Math.sin(q-2*k)+-1.2E-4*Math.sin(2*q)+1.1E-4*Math.sin(2*q-2*p-2*d))+(t*Math.cos(f)+x*Math.sin(f))*Math.tan(l)),e.RAD2DEG*h,e.RAD2DEG*b,g,2*e.RAD2DEG*Math.atan(1737.4/Math.sqrt(g*g-1737.4*1737.4)))};var eb;e.SiderealTime=bc;var E=function(a,b,c,d){this.x=a;this.y=b;this.z=c;this.t=d};E.prototype.Length=function(){return Math.hypot(this.x,this.y,this.z)};e.Vector=E;var I=function(a,b,c,d,f,g,h){this.x=a;this.y=b;this.z=c;this.vx= d;this.vy=f;this.vz=g;this.t=h};e.StateVector=I;var Ba=function(a,b,c){this.lat=w(a);this.lon=w(b);this.dist=w(c)};e.Spherical=Ba;var gb=function(a,b,c,d){this.ra=w(a);this.dec=w(b);this.dist=w(c);this.vec=d};e.EquatorialCoordinates=gb;var J=function(a){this.rot=a};e.RotationMatrix=J;e.MakeRotation=function(a){if(!qd(a))throw"Argument must be a [3][3] array of numbers";return new J(a)};var ec=function(a,b,c,d){this.azimuth=w(a);this.altitude=w(b);this.ra=w(c);this.dec=w(d)};e.HorizontalCoordinates= ec;var gc=function(a,b,c){this.vec=a;this.elat=w(b);this.elon=w(c)};e.EclipticCoordinates=gc;e.Horizon=hb;var zb=function(a,b,c){this.latitude=a;this.longitude=b;this.height=c;za(this)};e.Observer=zb;e.SunPosition=fc;e.Equator=Pa;e.ObserverVector=function(a,b,c){a=v(a);var d=ca(a);b=xb(b,d).pos;c||(b=fb(b,a,F.Into2000));return yb(b,a)};e.ObserverState=function(a,b,c){a=v(a);var d=ca(a);b=xb(b,d);b=new I(b.pos[0],b.pos[1],b.pos[2],b.vel[0],b.vel[1],b.vel[2],a);return c?b:(c=F.Into2000,c===F.Into2000? -(d=Na(a,c),b=Ea(d,b),a=Ma(a,c),a=Ea(a,b)):(d=Ma(a,c),b=Ea(d,b),a=Na(a,c),a=Ea(a,b)),a)};e.VectorObserver=function(a,b){var c=ca(a.t),d=[a.x,a.y,a.z];b||(d=wa(d,a.t,F.From2000),d=xa(d,a.t,F.From2000));b=d[0]*e.KM_PER_AU;var f=d[1]*e.KM_PER_AU;d=d[2]*e.KM_PER_AU;a=Math.hypot(b,f);if(1E-6>a){c=0;var g=0=c;)c+=360;for(;180Math.abs(q))break;g-=q/(-42.69778487239616*((k-n)/h-n*k*-.006694397995865464/(-42.69778487239616*p))+d*f+a*b)}g*=e.RAD2DEG;h=6378.1366/h;d=Math.abs(f)>Math.abs(b)?d/f-.9933056020041345*h:a/b-h}return new zb(g,c,1E3*d)};e.ObserverGravity=function(a,b){a=Math.sin(a*e.DEG2RAD);a*=a;return 9.7803253359*(1+.00193185265241*a)/Math.sqrt(1-.00669437999013*a)*(1-(3.15704E-7-2.10269E-9* -a)*b+7.37452E-14*b*b)};e.Ecliptic=Qa;e.GeoMoon=Y;e.EclipticGeoMoon=ib;e.GeoMoonState=Ra;e.GeoEmbState=Bb;var na=[[-73E4,[-26.118207232108,-14.376168177825,3.384402515299],[.0016339372163656,-.0027861699588508,-.0013585880229445]],[-700800,[41.974905202127,-.448502952929,-12.770351505989],[7.3458569351457E-4,.0022785014891658,4.8619778602049E-4]],[-671600,[14.706930780744,44.269110540027,9.353698474772],[-.00210001479998,2.2295915939915E-4,7.0143443551414E-4]],[-642400,[-29.441003929957,-6.43016153057, -6.858481011305],[8.4495803960544E-4,-.0030783914758711,-.0012106305981192]],[-613200,[39.444396946234,-6.557989760571,-13.913760296463],[.0011480029005873,.0022400006880665,3.5168075922288E-4]],[-584E3,[20.2303809507,43.266966657189,7.382966091923],[-.0019754081700585,5.3457141292226E-4,7.5929169129793E-4]],[-554800,[-30.65832536462,2.093818874552,9.880531138071],[6.1010603013347E-5,-.0031326500935382,-9.9346125151067E-4]],[-525600,[35.737703251673,-12.587706024764,-14.677847247563],[.0015802939375649, -.0021347678412429,1.9074436384343E-4]],[-496400,[25.466295188546,41.367478338417,5.216476873382],[-.0018054401046468,8.328308359951E-4,8.0260156912107E-4]],[-467200,[-29.847174904071,10.636426313081,12.297904180106],[-6.3257063052907E-4,-.0029969577578221,-7.4476074151596E-4]],[-438E3,[30.774692107687,-18.236637015304,-14.945535879896],[.0020113162005465,.0019353827024189,-2.0937793168297E-6]],[-408800,[30.243153324028,38.656267888503,2.938501750218],[-.0016052508674468,.0011183495337525,8.3333973416824E-4]], -[-379600,[-27.288984772533,18.643162147874,14.023633623329],[-.0011856388898191,-.0027170609282181,-4.9015526126399E-4]],[-350400,[24.519605196774,-23.245756064727,-14.626862367368],[.0024322321483154,.0016062008146048,-2.3369181613312E-4]],[-321200,[34.505274805875,35.125338586954,.557361475637],[-.0013824391637782,.0013833397561817,8.4823598806262E-4]],[-292E3,[-23.275363915119,25.818514298769,15.055381588598],[-.0016062295460975,-.0023395961498533,-2.4377362639479E-4]],[-262800,[17.050384798092, --27.180376290126,-13.608963321694],[.0028175521080578,.0011358749093955,-4.9548725258825E-4]],[-233600,[38.093671910285,30.880588383337,-1.843688067413],[-.0011317697153459,.0016128814698472,8.4177586176055E-4]],[-204400,[-18.197852930878,31.932869934309,15.438294826279],[-.0019117272501813,-.0019146495909842,-1.9657304369835E-5]],[-175200,[8.528924039997,-29.618422200048,-11.805400994258],[.0031034370787005,5.139363329243E-4,-7.7293066202546E-4]],[-146E3,[40.94685725864,25.904973592021,-4.256336240499], -[-8.3652705194051E-4,.0018129497136404,8.156422827306E-4]],[-116800,[-12.326958895325,36.881883446292,15.217158258711],[-.0021166103705038,-.001481442003599,1.7401209844705E-4]],[-87600,[-.633258375909,-30.018759794709,-9.17193287495],[.0032016994581737,-2.5279858672148E-4,-.0010411088271861]],[-58400,[42.936048423883,20.344685584452,-6.588027007912],[-5.0525450073192E-4,.0019910074335507,7.7440196540269E-4]],[-29200,[-5.975910552974,40.61180995846,14.470131723673],[-.0022184202156107,-.0010562361130164, -3.3652250216211E-4]],[0,[-9.875369580774,-27.978926224737,-5.753711824704],[.0030287533248818,-.0011276087003636,-.0012651326732361]],[29200,[43.958831986165,14.214147973292,-8.808306227163],[-1.4717608981871E-4,.0021404187242141,7.1486567806614E-4]],[58400,[.67813676352,43.094461639362,13.243238780721],[-.0022358226110718,-6.3233636090933E-4,4.7664798895648E-4]],[87600,[-18.282602096834,-23.30503958666,-1.766620508028],[.0025567245263557,-.0019902940754171,-.0013943491701082]],[116800,[43.873338744526, -7.700705617215,-10.814273666425],[2.3174803055677E-4,.0022402163127924,6.2988756452032E-4]],[146E3,[7.392949027906,44.382678951534,11.629500214854],[-.002193281545383,-2.1751799585364E-4,5.9556516201114E-4]],[175200,[-24.981690229261,-16.204012851426,2.466457544298],[.001819398914958,-.0026765419531201,-.0013848283502247]],[204400,[42.530187039511,.845935508021,-12.554907527683],[6.5059779150669E-4,.0022725657282262,5.1133743202822E-4]],[233600,[13.999526486822,44.462363044894,9.669418486465],[-.0021079296569252, -1.7533423831993E-4,6.9128485798076E-4]],[262800,[-29.184024803031,-7.371243995762,6.493275957928],[9.3581363109681E-4,-.0030610357109184,-.0012364201089345]],[292E3,[39.831980671753,-6.078405766765,-13.909815358656],[.0011117769689167,.0022362097830152,3.6230548231153E-4]],[321200,[20.294955108476,43.417190420251,7.450091985932],[-.0019742157451535,5.3102050468554E-4,7.5938408813008E-4]],[350400,[-30.66999230216,2.318743558955,9.973480913858],[4.5605107450676E-5,-.0031308219926928,-9.9066533301924E-4]], -[379600,[35.626122155983,-12.897647509224,-14.777586508444],[.0016015684949743,.0021171931182284,1.8002516202204E-4]],[408800,[26.133186148561,41.232139187599,5.00640132622],[-.0017857704419579,8.6046232702817E-4,8.0614690298954E-4]],[438E3,[-29.57674022923,11.863535943587,12.631323039872],[-7.2292830060955E-4,-.0029587820140709,-7.08242964503E-4]],[467200,[29.910805787391,-19.159019294,-15.013363865194],[.0020871080437997,.0018848372554514,-3.8528655083926E-5]],[496400,[31.375957451819,38.050372720763, -2.433138343754],[-.0015546055556611,.0011699815465629,8.3565439266001E-4]],[525600,[-26.360071336928,20.662505904952,14.414696258958],[-.0013142373118349,-.0026236647854842,-4.2542017598193E-4]],[554800,[22.599441488648,-24.508879898306,-14.484045731468],[.0025454108304806,.0014917058755191,-3.0243665086079E-4]],[584E3,[35.877864013014,33.894226366071,-.224524636277],[-.0012941245730845,.0014560427668319,8.4762160640137E-4]],[613200,[-21.538149762417,28.204068269761,15.321973799534],[-.001731211740901, --.0021939631314577,-1.631691327518E-4]],[642400,[13.971521374415,-28.339941764789,-13.083792871886],[.0029334630526035,9.1860931752944E-4,-5.9939422488627E-4]],[671600,[39.526942044143,28.93989736011,-2.872799527539],[-.0010068481658095,.001702113288809,8.3578230511981E-4]],[700800,[-15.576200701394,34.399412961275,15.466033737854],[-.0020098814612884,-.0017191109825989,7.0414782780416E-5]],[73E4,[4.24325283709,-30.118201690825,-10.707441231349],[.0031725847067411,1.609846120227E-4,-9.0672150593868E-4]]], -B=function(a,b,c){this.x=a;this.y=b;this.z=c};B.prototype.clone=function(){return new B(this.x,this.y,this.z)};B.prototype.ToAstroVector=function(a){return new E(this.x,this.y,this.z,a)};B.zero=function(){return new B(0,0,0)};B.prototype.quadrature=function(){return this.x*this.x+this.y*this.y+this.z*this.z};B.prototype.add=function(a){return new B(this.x+a.x,this.y+a.y,this.z+a.z)};B.prototype.sub=function(a){return new B(this.x-a.x,this.y-a.y,this.z-a.z)};B.prototype.incr=function(a){this.x+=a.x; -this.y+=a.y;this.z+=a.z};B.prototype.decr=function(a){this.x-=a.x;this.y-=a.y;this.z-=a.z};B.prototype.mul=function(a){return new B(a*this.x,a*this.y,a*this.z)};B.prototype.div=function(a){return new B(this.x/a,this.y/a,this.z/a)};B.prototype.mean=function(a){return new B((this.x+a.x)/2,(this.y+a.y)/2,(this.z+a.z)/2)};B.prototype.neg=function(){return new B(-this.x,-this.y,-this.z)};var X=function(a,b,c){this.tt=a;this.r=b;this.v=c};X.prototype.clone=function(){return new X(this.tt,this.r,this.v)}; -X.prototype.sub=function(a){return new X(this.tt,this.r.sub(a.r),this.v.sub(a.v))};var Da=function(a){var b=new X(a,new B(0,0,0),new B(0,0,0));this.Jupiter=R(b,a,m.Jupiter,2.825345909524226E-7);this.Saturn=R(b,a,m.Saturn,8.459715185680659E-8);this.Uranus=R(b,a,m.Uranus,1.292024916781969E-8);this.Neptune=R(b,a,m.Neptune,1.524358900784276E-8);this.Jupiter.r.decr(b.r);this.Jupiter.v.decr(b.v);this.Saturn.r.decr(b.r);this.Saturn.v.decr(b.v);this.Uranus.r.decr(b.r);this.Uranus.v.decr(b.v);this.Neptune.r.decr(b.r); -this.Neptune.v.decr(b.v);this.Sun=new X(a,b.r.mul(-1),b.v.mul(-1))};Da.prototype.Acceleration=function(a){var b=Ta(a,2.959122082855911E-4,this.Sun.r);b.incr(Ta(a,2.825345909524226E-7,this.Jupiter.r));b.incr(Ta(a,8.459715185680659E-8,this.Saturn.r));b.incr(Ta(a,1.292024916781969E-8,this.Uranus.r));b.incr(Ta(a,1.524358900784276E-8,this.Neptune.r));return b};var Ua=function(a,b,c,d){this.tt=a;this.r=b;this.v=c;this.a=d};Ua.prototype.clone=function(){return new Ua(this.tt,this.r.clone(),this.v.clone(), -this.a.clone())};var ic=function(a,b){this.bary=a;this.grav=b},Ib=[],rd=new J([[.999432765338654,-.0336771074697641,0],[.0303959428906285,.902057912352809,.430543388542295],[-.0144994559663353,-.430299169409101,.902569881273754]]),tb=[{mu:2.82489428433814E-7,al:[1.446213296021224,3.5515522861824],a:[[.0028210960212903,0,0]],l:[[-1.925258348666E-4,4.9369589722645,.01358483658305],[-9.70803596076E-5,4.3188796477322,.01303413843243],[-8.988174165E-5,1.9080016428617,.00305064867158],[-5.53101050262E-5, +(d=Na(a,c),b=Ea(d,b),a=Ma(a,c),a=Ea(a,b)):(d=Ma(a,c),b=Ea(d,b),a=Na(a,c),a=Ea(a,b)),a)};e.VectorObserver=function(a,b){var c=ca(a.t),d=[a.x,a.y,a.z];b||(d=wa(d,a.t,F.From2000),d=xa(d,a.t,F.From2000));var f=d[0]*e.KM_PER_AU,g=d[1]*e.KM_PER_AU;a=d[2]*e.KM_PER_AU;b=Math.hypot(f,g);if(1E-6>b)c=0,g=0=c;)c+=360;for(;180a&&(a+=360);return a}var sidereal_time_cache; +k,p=-h*e-g*l*d,q=-h*g*c+l*d*e*c-a*l*k;h=-h*g*a+l*d*e*a+c*l*k;g*=k;l=-k*e*c-a*d;a=-k*e*a+d*c;if(b===PrecessDirection.Into2000)return new RotationMatrix([[f,m,n],[p,q,h],[g,l,a]]);if(b===PrecessDirection.From2000)return new RotationMatrix([[f,p,g],[m,q,l],[n,h,a]]);throw"Invalid precess direction";}function era(a){a=(.779057273264+.00273781191135448*a.ut+a.ut%1)%1*360;0>a&&(a+=360);return a}var sidereal_time_cache; function sidereal_time(a){if(!sidereal_time_cache||sidereal_time_cache.tt!==a.tt){var b=a.tt/36525,c=15*e_tilt(a).ee,d=era(a);b=((c+.014506+((((-3.68E-8*b-2.9956E-5)*b-4.4E-7)*b+1.3915817)*b+4612.156534)*b)/3600+d)%360/15;0>b&&(b+=24);sidereal_time_cache={tt:a.tt,st:b}}return sidereal_time_cache.st}function SiderealTime(a){a=MakeTime(a);return sidereal_time(a)}exports.SiderealTime=SiderealTime; -function inverse_terra(a,b){var c=a[0]*exports.KM_PER_AU,d=a[1]*exports.KM_PER_AU;a=a[2]*exports.KM_PER_AU;var e=Math.hypot(c,d);if(1E-6>e){b=0;var f=0=b;)b+=360;for(;180Math.abs(p))break;f-=p/(h*((l-m)/g-m*l*(EARTH_FLATTENING_SQUARED-1)/(h*n))+a*d+e*c)}f*=exports.RAD2DEG;g=EARTH_EQUATORIAL_RADIUS_KM/g;a=Math.abs(d)>Math.abs(c)?a/d-EARTH_FLATTENING_SQUARED*g:e/c-g}return new Observer(f,b,1E3*a)} +function inverse_terra(a,b){var c=a[0]*exports.KM_PER_AU,d=a[1]*exports.KM_PER_AU,e=a[2]*exports.KM_PER_AU,f=Math.hypot(c,d);if(1E-6>f)b=0,d=0=b;)b+=360;for(;180d[1]&&(d[1]+=PI2);e=$jscomp.makeIterator(b.z);for(f=e.next();!f.done;f= e.next())g=$jscomp.makeIterator(f.value),f=g.next().value,k=g.next().value,g=g.next().value,k+=c*g,d[2]+=f*Math.cos(k),d[3]+=f*Math.sin(k);e=$jscomp.makeIterator(b.zeta);for(f=e.next();!f.done;f=e.next())g=$jscomp.makeIterator(f.value),f=g.next().value,k=g.next().value,g=g.next().value,k+=c*g,d[4]+=f*Math.cos(k),d[5]+=f*Math.sin(k);a=JupiterMoon_elem2pv(a,b.mu,d);return RotateState(Rotation_JUP_EQJ,a)} @@ -204,8 +204,8 @@ GeoMoonState(b):GeoEmbState(b);return new StateVector(a.x+c.Sun.r.x+d.r.x,a.y+c. function HelioState(a,b){b=MakeTime(b);switch(a){case Body.Sun:return new StateVector(0,0,0,0,0,0,b);case Body.SSB:return a=new major_bodies_t(b.tt),new StateVector(-a.Sun.r.x,-a.Sun.r.y,-a.Sun.r.z,-a.Sun.v.x,-a.Sun.v.y,-a.Sun.v.z,b);case Body.Mercury:case Body.Venus:case Body.Earth:case Body.Mars:case Body.Jupiter:case Body.Saturn:case Body.Uranus:case Body.Neptune:return a=CalcVsopPosVel(vsop[a],b.tt),ExportState(a,b);case Body.Pluto:return CalcPluto(b,!0);case Body.Moon:case Body.EMB:var c=CalcVsopPosVel(vsop.Earth, b.tt);a=a==Body.Moon?GeoMoonState(b):GeoEmbState(b);return new StateVector(a.x+c.r.x,a.y+c.r.y,a.z+c.r.z,a.vx+c.v.x,a.vy+c.v.y,a.vz+c.v.z,b);default:if(UserDefinedStar(a))return a=HelioVector(a,b),new StateVector(a.x,a.y,a.z,0,0,0,b);throw'HelioState: Unsupported body "'+a+'"';}}exports.HelioState=HelioState; function QuadInterp(a,b,c,d,e){var f=(e+c)/2-d;c=(e-c)/2;if(0==f){if(0==c)return null;d=-d/c;if(-1>d||1=d)return null;e=Math.sqrt(d);d=(-c+e)/(2*f);e=(-c-e)/(2*f);if(-1<=d&&1>=d){if(-1<=e&&1>=e)return null}else if(-1<=e&&1>=e)d=e;else return null}return{t:a+d*b,df_dt:(2*f*d+c)/b}} -function Search(a,b,c,d){var e=VerifyNumber(d&&d.dt_tolerance_seconds||1);e=Math.abs(e/SECONDS_PER_DAY);var f=d&&d.init_f1||a(b),g=d&&d.init_f2||a(c),k=NaN,h=0;d=d&&d.iter_limit||20;for(var l=!0;;){if(++h>d)throw"Excessive iteration in Search()";var m=InterpolateTime(b,c,.5),n=m.ut-b.ut;if(Math.abs(n)(n.ut-b.ut)*(n.ut-c.ut)&&0>(r.ut-b.ut)*(r.ut-c.ut))){p=a(n);var t=a(r);if(0>p&&0<=t){f=p;g=t;b=n;c=r;k=q;l=!1;continue}}}}if(0>f&&0<=k)c=m,g=k;else if(0>k&&0<=g)b=m,f=k;else return null}}exports.Search=Search;function LongitudeOffset(a){for(;-180>=a;)a+=360;for(;180a;)a+=360;for(;360<=a;)a-=360;return a} +function Search(a,b,c,d){var e=VerifyNumber(d&&d.dt_tolerance_seconds||1);e=Math.abs(e/SECONDS_PER_DAY);var f=d&&d.init_f1||a(b),g=d&&d.init_f2||a(c),k=NaN,h=0;d=d&&d.iter_limit||20;for(var l=!0;;){if(++h>d)throw"Excessive iteration in Search()";var m=InterpolateTime(b,c,.5),n=m.ut-b.ut;if(Math.abs(n)(n.ut-b.ut)*(n.ut-c.ut)&&0>(q.ut-b.ut)*(q.ut-c.ut))){p=a(n);var t=a(q);if(0>p&&0<=t){f=p;g=t;b=n;c=q;k=r;l=!1;continue}}}}if(0>f&&0<=k)c=m,g=k;else if(0>k&&0<=g)b=m,f=k;else return null}}exports.Search=Search;function LongitudeOffset(a){for(;-180>=a;)a+=360;for(;180a;)a+=360;for(;360<=a;)a-=360;return a} function SearchSunLongitude(a,b,c){VerifyNumber(a);VerifyNumber(c);b=MakeTime(b);c=b.AddDays(c);return Search(function(d){d=SunPosition(d);return LongitudeOffset(d.elon-a)},b,c,{dt_tolerance_seconds:.01})}exports.SearchSunLongitude=SearchSunLongitude; function PairLongitude(a,b,c){if(a===Body.Earth||b===Body.Earth)throw"The Earth does not have a longitude as seen from itself.";c=MakeTime(c);a=GeoVector(a,c,!1);a=Ecliptic(a);b=GeoVector(b,c,!1);b=Ecliptic(b);return NormalizeLongitude(a.elon-b.elon)}exports.PairLongitude=PairLongitude;function AngleFromSun(a,b){if(a==Body.Earth)throw"The Earth does not have an angle as seen from itself.";var c=MakeTime(b);b=GeoVector(Body.Sun,c,!0);a=GeoVector(a,c,!0);return AngleBetween(b,a)} exports.AngleFromSun=AngleFromSun;function EclipticLongitude(a,b){if(a===Body.Sun)throw"Cannot calculate heliocentric longitude of the Sun.";a=HelioVector(a,b);return Ecliptic(a).elon}exports.EclipticLongitude=EclipticLongitude; @@ -225,8 +225,8 @@ function SearchAltitude(a,b,c,d,e,f){if(!Number.isFinite(f)||-90>f||90f&&0<=g)return new AscentInfo(d,e,f,g);if(0<=f&&0>g)return null;if(17k*SECONDS_PER_DAY||Math.min(Math.abs(f),Math.abs(g))>k/2*c)return null;k=new AstroTime((d.ut+e.ut)/2);var h=b(k);return FindAscent(1+a,b,c,d,k,f,h)||FindAscent(1+a,b,c,k,e,h,g)} function MaxAltitudeSlope(a,b){if(-90>b||90g||90e?(l=m.AddDays(-.42),n=k(l)):(m=l.AddDays(.42),p=k(m));var r=FindAscent(0,k,h,l,m,n,p);if(r){if(h=Search(k,r.tx,r.ty,{dt_tolerance_seconds:.1, -init_f1:r.ax,init_f2:r.ay})){if(0>e){if(h.utd.ut+e)return null;return h}throw"Rise/set search failed after finding ascent: t1="+l+", t2="+m+", a1="+n+", a2="+p;}if(0>e){if(l.utd.ut+e)return null;l=m;n=p}}}var HourAngleEvent=function(a,b){this.time=a;this.hor=b};exports.HourAngleEvent=HourAngleEvent; +function InternalSearchAltitude(a,b,c,d,e,f,g){function k(r){var t=Equator(a,r,b,!0,!0);r=Horizon(r,b,t.ra,t.dec).altitude+exports.RAD2DEG*Math.asin(f/t.dist);return c*(r-g)}VerifyObserver(b);VerifyNumber(e);VerifyNumber(f);VerifyNumber(g);if(-90>g||90e?(l=m.AddDays(-.42),n=k(l)):(m=l.AddDays(.42),p=k(m));var q=FindAscent(0,k,h,l,m,n,p);if(q){if(h=Search(k,q.tx,q.ty,{dt_tolerance_seconds:.1, +init_f1:q.ax,init_f2:q.ay})){if(0>e){if(h.utd.ut+e)return null;return h}throw"Rise/set search failed after finding ascent: t1="+l+", t2="+m+", a1="+n+", a2="+p;}if(0>e){if(l.utd.ut+e)return null;l=m;n=p}}}var HourAngleEvent=function(a,b){this.time=a;this.hor=b};exports.HourAngleEvent=HourAngleEvent; function SearchHourAngle(a,b,c,d,e){e=void 0===e?1:e;VerifyObserver(b);d=MakeTime(d);var f=0;if(a===Body.Earth)throw"Cannot search for hour angle of the Earth.";VerifyNumber(c);if(0>c||24<=c)throw"Invalid hour angle "+c;VerifyNumber(e);if(0===e)throw"Direction must be positive or negative.";for(;;){++f;var g=sidereal_time(d),k=Equator(a,d,b,!0,!0);g=(c+k.ra-b.longitude/15-g)%24;1===f?0g&&(g+=24):0g?g+=24:123600*Math.abs(g))return a=Horizon(d,b,k.ra,k.dec,"normal"), new HourAngleEvent(d,a);d=d.AddDays(g/24*SOLAR_DAYS_PER_SIDEREAL_DAY)}}exports.SearchHourAngle=SearchHourAngle;function HourAngle(a,b,c){var d=MakeTime(b);b=SiderealTime(d);a=Equator(a,d,c,!0,!0);c=(c.longitude/15+b-a.ra)%24;0>c&&(c+=24);return c}exports.HourAngle=HourAngle;var SeasonInfo=function(a,b,c,d){this.mar_equinox=a;this.jun_solstice=b;this.sep_equinox=c;this.dec_solstice=d};exports.SeasonInfo=SeasonInfo; function Seasons(a){function b(g,k,h){k=new Date(Date.UTC(a,k-1,h));g=SearchSunLongitude(g,k,20);if(!g)throw"Cannot find season change near "+k.toISOString();return g}a instanceof Date&&Number.isFinite(a.getTime())&&(a=a.getUTCFullYear());if(!Number.isSafeInteger(a))throw"Cannot calculate seasons because year argument "+a+" is neither a Date nor a safe integer.";var c=b(0,3,10),d=b(90,6,10),e=b(180,9,10),f=b(270,12,10);return new SeasonInfo(c,d,e,f)}exports.Seasons=Seasons; @@ -314,17 +314,17 @@ d;c=e}}exports.SearchMoonNode=SearchMoonNode; function NextMoonNode(a){var b=a.time.AddDays(MoonNodeStepDays);b=SearchMoonNode(b);switch(a.kind){case NodeEventKind.Ascending:if(b.kind!==NodeEventKind.Descending)throw"Internal error: previous node was ascending, but this node was: "+b.kind;break;case NodeEventKind.Descending:if(b.kind!==NodeEventKind.Ascending)throw"Internal error: previous node was descending, but this node was: "+b.kind;break;default:throw"Previous node has an invalid node kind: "+a.kind;}return b}exports.NextMoonNode=NextMoonNode; var AxisInfo=function(a,b,c,d){this.ra=a;this.dec=b;this.spin=c;this.north=d};exports.AxisInfo=AxisInfo;function EarthRotationAxis(a){var b=nutation([0,0,1],a,PrecessDirection.Into2000);b=precession(b,a,PrecessDirection.Into2000);b=new Vector(b[0],b[1],b[2],a);var c=EquatorFromVector(b);return new AxisInfo(c.ra,c.dec,190.41375788700253+360.9856122880876*a.ut,b)} function RotationAxis(a,b){b=MakeTime(b);var c=b.tt,d=c/36525;switch(a){case Body.Sun:a=286.13;var e=63.87;c=84.176+14.1844*c;break;case Body.Mercury:a=281.0103-.0328*d;e=61.4155-.0049*d;c=329.5988+6.1385108*c+.01067257*Math.sin(exports.DEG2RAD*(174.7910857+4.092335*c))-.00112309*Math.sin(exports.DEG2RAD*(349.5821714+8.18467*c))-1.104E-4*Math.sin(exports.DEG2RAD*(164.3732571+12.277005*c))-2.539E-5*Math.sin(exports.DEG2RAD*(339.1643429+16.36934*c))-5.71E-6*Math.sin(exports.DEG2RAD*(153.9554286+20.461675* -c));break;case Body.Venus:a=272.76;e=67.16;c=160.2-1.4813688*c;break;case Body.Earth:return EarthRotationAxis(b);case Body.Moon:var f=exports.DEG2RAD*(125.045-.0529921*c),g=exports.DEG2RAD*(250.089-.1059842*c),k=exports.DEG2RAD*(260.008+13.0120009*c),h=exports.DEG2RAD*(176.625+13.3407154*c),l=exports.DEG2RAD*(357.529+.9856003*c),m=exports.DEG2RAD*(311.589+26.4057084*c),n=exports.DEG2RAD*(134.963+13.064993*c),p=exports.DEG2RAD*(276.617+.3287146*c),r=exports.DEG2RAD*(34.226+1.7484877*c),q=exports.DEG2RAD* -(15.134-.1589763*c),t=exports.DEG2RAD*(119.743+.0036096*c),H=exports.DEG2RAD*(239.961+.1643573*c),F=exports.DEG2RAD*(25.053+12.9590088*c);a=269.9949+.0031*d-3.8787*Math.sin(f)-.1204*Math.sin(g)+.07*Math.sin(k)-.0172*Math.sin(h)+.0072*Math.sin(m)-.0052*Math.sin(q)+.0043*Math.sin(F);e=66.5392+.013*d+1.5419*Math.cos(f)+.0239*Math.cos(g)-.0278*Math.cos(k)+.0068*Math.cos(h)-.0029*Math.cos(m)+9E-4*Math.cos(n)+8E-4*Math.cos(q)-9E-4*Math.cos(F);c=38.3213+(13.17635815-1.4E-12*c)*c+3.561*Math.sin(f)+.1208* -Math.sin(g)-.0642*Math.sin(k)+.0158*Math.sin(h)+.0252*Math.sin(l)-.0066*Math.sin(m)-.0047*Math.sin(n)-.0046*Math.sin(p)+.0028*Math.sin(r)+.0052*Math.sin(q)+.004*Math.sin(t)+.0019*Math.sin(H)-.0044*Math.sin(F);break;case Body.Mars:a=317.269202-.10927547*d+6.8E-5*Math.sin(exports.DEG2RAD*(198.991226+19139.4819985*d))+2.38E-4*Math.sin(exports.DEG2RAD*(226.292679+38280.8511281*d))+5.2E-5*Math.sin(exports.DEG2RAD*(249.663391+57420.7251593*d))+9E-6*Math.sin(exports.DEG2RAD*(266.18351+76560.636795*d))+.419057* +c));break;case Body.Venus:a=272.76;e=67.16;c=160.2-1.4813688*c;break;case Body.Earth:return EarthRotationAxis(b);case Body.Moon:var f=exports.DEG2RAD*(125.045-.0529921*c),g=exports.DEG2RAD*(250.089-.1059842*c),k=exports.DEG2RAD*(260.008+13.0120009*c),h=exports.DEG2RAD*(176.625+13.3407154*c),l=exports.DEG2RAD*(357.529+.9856003*c),m=exports.DEG2RAD*(311.589+26.4057084*c),n=exports.DEG2RAD*(134.963+13.064993*c),p=exports.DEG2RAD*(276.617+.3287146*c),q=exports.DEG2RAD*(34.226+1.7484877*c),r=exports.DEG2RAD* +(15.134-.1589763*c),t=exports.DEG2RAD*(119.743+.0036096*c),H=exports.DEG2RAD*(239.961+.1643573*c),F=exports.DEG2RAD*(25.053+12.9590088*c);a=269.9949+.0031*d-3.8787*Math.sin(f)-.1204*Math.sin(g)+.07*Math.sin(k)-.0172*Math.sin(h)+.0072*Math.sin(m)-.0052*Math.sin(r)+.0043*Math.sin(F);e=66.5392+.013*d+1.5419*Math.cos(f)+.0239*Math.cos(g)-.0278*Math.cos(k)+.0068*Math.cos(h)-.0029*Math.cos(m)+9E-4*Math.cos(n)+8E-4*Math.cos(r)-9E-4*Math.cos(F);c=38.3213+(13.17635815-1.4E-12*c)*c+3.561*Math.sin(f)+.1208* +Math.sin(g)-.0642*Math.sin(k)+.0158*Math.sin(h)+.0252*Math.sin(l)-.0066*Math.sin(m)-.0047*Math.sin(n)-.0046*Math.sin(p)+.0028*Math.sin(q)+.0052*Math.sin(r)+.004*Math.sin(t)+.0019*Math.sin(H)-.0044*Math.sin(F);break;case Body.Mars:a=317.269202-.10927547*d+6.8E-5*Math.sin(exports.DEG2RAD*(198.991226+19139.4819985*d))+2.38E-4*Math.sin(exports.DEG2RAD*(226.292679+38280.8511281*d))+5.2E-5*Math.sin(exports.DEG2RAD*(249.663391+57420.7251593*d))+9E-6*Math.sin(exports.DEG2RAD*(266.18351+76560.636795*d))+.419057* Math.sin(exports.DEG2RAD*(79.398797+.5042615*d));e=54.432516-.05827105*d+5.1E-5*Math.cos(exports.DEG2RAD*(122.433576+19139.9407476*d))+1.41E-4*Math.cos(exports.DEG2RAD*(43.058401+38280.8753272*d))+3.1E-5*Math.cos(exports.DEG2RAD*(57.663379+57420.7517205*d))+5E-6*Math.cos(exports.DEG2RAD*(79.476401+76560.6495004*d))+1.591274*Math.cos(exports.DEG2RAD*(166.325722+.5042615*d));c=176.049863+350.891982443297*c+1.45E-4*Math.sin(exports.DEG2RAD*(129.071773+19140.0328244*d))+1.57E-4*Math.sin(exports.DEG2RAD* (36.352167+38281.0473591*d))+4E-5*Math.sin(exports.DEG2RAD*(56.668646+57420.929536*d))+1E-6*Math.sin(exports.DEG2RAD*(67.364003+76560.2552215*d))+1E-6*Math.sin(exports.DEG2RAD*(104.79268+95700.4387578*d))+.584542*Math.sin(exports.DEG2RAD*(95.391654+.5042615*d));break;case Body.Jupiter:e=exports.DEG2RAD*(99.360714+4850.4046*d);f=exports.DEG2RAD*(175.895369+1191.9605*d);g=exports.DEG2RAD*(300.323162+262.5475*d);k=exports.DEG2RAD*(114.012305+6070.2476*d);h=exports.DEG2RAD*(49.511251+64.3*d);a=268.056595- .006499*d+1.17E-4*Math.sin(e)+9.38E-4*Math.sin(f)+.001432*Math.sin(g)+3E-5*Math.sin(k)+.00215*Math.sin(h);e=64.495303+.002413*d+5E-5*Math.cos(e)+4.04E-4*Math.cos(f)+6.17E-4*Math.cos(g)-1.3E-5*Math.cos(k)+9.26E-4*Math.cos(h);c=284.95+870.536*c;break;case Body.Saturn:a=40.589-.036*d;e=83.537-.004*d;c=38.9+810.7939024*c;break;case Body.Uranus:a=257.311;e=-15.175;c=203.81-501.1600928*c;break;case Body.Neptune:d=exports.DEG2RAD*(357.85+52.316*d);a=299.36+.7*Math.sin(d);e=43.46-.51*Math.cos(d);c=249.978+ 541.1397757*c-.48*Math.sin(d);break;case Body.Pluto:a=132.993;e=-6.163;c=302.695+56.3625225*c;break;default:throw"Invalid body: "+a;}d=e*exports.DEG2RAD;f=a*exports.DEG2RAD;g=Math.cos(d);b=new Vector(g*Math.cos(f),g*Math.sin(f),Math.sin(d),b);return new AxisInfo(a/15,e,c,b)}exports.RotationAxis=RotationAxis; function LagrangePoint(a,b,c,d){var e=MakeTime(b);b=MassProduct(c);var f=MassProduct(d);c===Body.Earth&&d===Body.Moon?(c=new StateVector(0,0,0,0,0,0,e),d=GeoMoonState(e)):(c=HelioState(c,e),d=HelioState(d,e));return LagrangePointFast(a,c,b,d,f)}exports.LagrangePoint=LagrangePoint; -function LagrangePointFast(a,b,c,d,e){if(1>a||5=c)throw"Major mass must be a positive number.";if(!Number.isFinite(e)||0>=e)throw"Minor mass must be a negative number.";var f=d.x-b.x,g=d.y-b.y,k=d.z-b.z,h=f*f+g*g+k*k,l=Math.sqrt(h),m=d.vx-b.vx,n=d.vy-b.vy;d=d.vz-b.vz;if(4===a||5===a){h=g*d-k*n;c=k*m-f*d;var p=f*n-g*m,r=c*k-p*g;p=p*f-h*k;h=h*g-c*f;c=Math.sqrt(r*r+p*p+h*h);r/=c;p/=c;h/=c;f/=l;g/=l;k/=l;a=4==a?.8660254037844386:-.8660254037844386; -c=.5*f+a*r;e=.5*g+a*p;var q=.5*k+a*h,t=m*f+n*g+d*k;m=m*r+n*p+d*h;b=new StateVector(l*c,l*e,l*q,t*c+m*(.5*r-a*f),t*e+m*(.5*p-a*g),t*q+m*(.5*h-a*k),b.t)}else{r=e/(c+e)*-l;p=c/(c+e)*+l;h=(c+e)/(h*l);if(1===a||2===a)q=c/(c+e)*Math.cbrt(e/(3*c)),c=-c,1==a?(q=1-q,a=+e):(q=1+q,a=-e);else if(3===a)q=(7/12*e-c)/(e+c),c=+c,a=+e;else throw"Invalid Langrage point "+a+". Must be an integer 1..5.";e=l*q-r;do q=e-r,t=e-p,q=(h*e+c/(q*q)+a/(t*t))/(h-2*c/(q*q*q)-2*a/(t*t*t)),e-=q;while(1E-14a||5=c)throw"Major mass must be a positive number.";if(!Number.isFinite(e)||0>=e)throw"Minor mass must be a negative number.";var f=d.x-b.x,g=d.y-b.y,k=d.z-b.z,h=f*f+g*g+k*k,l=Math.sqrt(h),m=d.vx-b.vx,n=d.vy-b.vy;d=d.vz-b.vz;if(4===a||5===a){h=g*d-k*n;c=k*m-f*d;var p=f*n-g*m,q=c*k-p*g;p=p*f-h*k;h=h*g-c*f;c=Math.sqrt(q*q+p*p+h*h);q/=c;p/=c;h/=c;f/=l;g/=l;k/=l;a=4==a?.8660254037844386:-.8660254037844386; +c=.5*f+a*q;e=.5*g+a*p;var r=.5*k+a*h,t=m*f+n*g+d*k;m=m*q+n*p+d*h;b=new StateVector(l*c,l*e,l*r,t*c+m*(.5*q-a*f),t*e+m*(.5*p-a*g),t*r+m*(.5*h-a*k),b.t)}else{q=e/(c+e)*-l;p=c/(c+e)*+l;h=(c+e)/(h*l);if(1===a||2===a)r=c/(c+e)*Math.cbrt(e/(3*c)),c=-c,1==a?(r=1-r,a=+e):(r=1+r,a=-e);else if(3===a)r=(7/12*e-c)/(e+c),c=+c,a=+e;else throw"Invalid Langrage point "+a+". Must be an integer 1..5.";e=l*r-q;do r=e-q,t=e-p,r=(h*e+c/(r*r)+a/(t*t))/(h-2*c/(r*r*r)-2*a/(t*t*t)),e-=r;while(1E-14 10) - throw `inverse_terra failed to converge.`; + if (++count > 10) { + throw `inverse_terra failed to converge: W=${W}, distanceAu=${distanceAu}`; + } // Calculate the error function W(lat). // We try to find the root of W, meaning where the error is 0. cos = Math.cos(lat); @@ -2120,8 +2123,8 @@ function inverse_terra(ovec: ArrayVector, st: number): Observer { const sin2 = sin*sin; const radicand = cos2 + EARTH_FLATTENING_SQUARED*sin2; denom = Math.sqrt(radicand); - const W = (factor*sin*cos)/denom - z*cos + p*sin; - if (Math.abs(W) < 2.0e-8) + W = (factor*sin*cos)/denom - z*cos + p*sin; + if (Math.abs(W) < distanceAu * 2.0e-8) break; // The error is now negligible // Error is still too large. Find the next estimate. // Calculate D = the derivative of W with respect to lat. diff --git a/source/js/esm/astronomy.js b/source/js/esm/astronomy.js index d095213c..2722857d 100644 --- a/source/js/esm/astronomy.js +++ b/source/js/esm/astronomy.js @@ -1920,9 +1920,12 @@ function inverse_terra(ovec, st) { let lat = Math.atan2(z, p); let cos, sin, denom; let count = 0; + let W = 0; + const distanceAu = Math.max(1, Math.hypot(ovec[0], ovec[1], ovec[2])); for (;;) { - if (++count > 10) - throw `inverse_terra failed to converge.`; + if (++count > 10) { + throw `inverse_terra failed to converge: W=${W}, distanceAu=${distanceAu}`; + } // Calculate the error function W(lat). // We try to find the root of W, meaning where the error is 0. cos = Math.cos(lat); @@ -1932,8 +1935,8 @@ function inverse_terra(ovec, st) { const sin2 = sin * sin; const radicand = cos2 + EARTH_FLATTENING_SQUARED * sin2; denom = Math.sqrt(radicand); - const W = (factor * sin * cos) / denom - z * cos + p * sin; - if (Math.abs(W) < 2.0e-8) + W = (factor * sin * cos) / denom - z * cos + p * sin; + if (Math.abs(W) < distanceAu * 2.0e-8) break; // The error is now negligible // Error is still too large. Find the next estimate. // Calculate D = the derivative of W with respect to lat. diff --git a/source/kotlin/src/main/kotlin/io/github/cosinekitty/astronomy/astronomy.kt b/source/kotlin/src/main/kotlin/io/github/cosinekitty/astronomy/astronomy.kt index a78fd16a..0ff7d11d 100644 --- a/source/kotlin/src/main/kotlin/io/github/cosinekitty/astronomy/astronomy.kt +++ b/source/kotlin/src/main/kotlin/io/github/cosinekitty/astronomy/astronomy.kt @@ -4370,6 +4370,7 @@ private fun inverseTerra(ovec: Vector): Observer { var s: Double var denom: Double var count = 0 + val distanceAu = max(1.0, ovec.length()) while (true) { ++count if (count > 10) @@ -4384,7 +4385,7 @@ private fun inverseTerra(ovec: Vector): Observer { val radicand = c2 + F*s2 denom = sqrt(radicand) val W = ((factor * s * c) / denom) - (z * c) + (p * s) - if (W.absoluteValue < 2.0e-8) + if (W.absoluteValue < distanceAu * 2.0e-8) break // The error is now negligible. // Error is still too large. Find the next estimate. // Calculate D = the derivative of W with respect to lat. diff --git a/source/python/astronomy/astronomy.py b/source/python/astronomy/astronomy.py index d20ae486..e4b922eb 100644 --- a/source/python/astronomy/astronomy.py +++ b/source/python/astronomy/astronomy.py @@ -1435,6 +1435,7 @@ def _inverse_terra(ovec: List[float], st: float) -> Observer: # Start with initial latitude estimate, based on a spherical Earth. lat = math.atan2(z, p) count = 0 + distanceAu = max(1.0, math.hypot(ovec[0], ovec[1], ovec[2])) while True: count += 1 if count > 10: @@ -1449,7 +1450,7 @@ def _inverse_terra(ovec: List[float], st: float) -> Observer: radicand = cos2 + _EARTH_FLATTENING_SQUARED*sin2 denom = math.sqrt(radicand) W = (factor*sin*cos)/denom - z*cos + p*sin - if abs(W) < 2.0e-8: + if abs(W) < distanceAu * 2.0e-8: # The error is now negligible break # Error is still too large. Find the next estimate.