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

@@ -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 <!--MINIFIED_SIZE-->116424 bytes.)
+- Minified JavaScript code less than 120K. (The current size is <!--MINIFIED_SIZE-->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.
 

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.

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.

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.

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.

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. */

generate/template/astronomy.cs

@@ -4162,6 +4162,7 @@ private static Observer inverse_terra(AstroVector ovec)
                 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 @@ private static Observer inverse_terra(AstroVector ovec)
                     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.

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.

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.

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.

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)
+    );
 }
 
 

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. */

source/csharp/astronomy.cs

@@ -5296,6 +5296,7 @@ private static Observer inverse_terra(AstroVector ovec)
                 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 @@ private static Observer inverse_terra(AstroVector ovec)
                     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.

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.