Looking at _CalcVsop as the possible origin of numeric errors.

demo/python/astronomy.py

@@ -36,6 +36,7 @@
 import enum
 import re
 import abc
+from os import getenv
 from typing import Any, List, Tuple, Optional, Union, Callable, Dict
 
 def _cbrt(x: float) -> float:
@@ -3180,12 +3181,16 @@ def _VsopSphereToRect(lon: float, lat: float, rad: float) -> _TerseVector:
         rad * math.sin(lat)
     )
 
+_ProblemDebug:bool = ('1' == getenv('ASTRONOMY_ENGINE_PYTHON_PROBLEM'))
+
 def _CalcVsop(model: _vsop_model_t, time: Time) -> Vector:
     t = time.tt / _DAYS_PER_MILLENNIUM
     lon = _VsopFormula(model.lon, t, True)
     lat = _VsopFormula(model.lat, t, False)
     rad = _VsopFormula(model.rad, t, False)
     eclip = _VsopSphereToRect(lon, lat, rad)
+    if _ProblemDebug:
+        print('_CalcVsop: lon={:0.16g}, lat={:0.16g}, rad={:0.16g}, eclip=({:0.16g}, {:0.16g}, {:0.16g})'.format(lon, lat, rad, eclip.x, eclip.y, eclip.z))
     return _VsopRotate(eclip).ToAstroVector(time)
 
 class _body_state_t:

generate/.gitignore

@@ -13,3 +13,4 @@ obj/
 profile/
 private/
 node_modules/
+precision_problem.txt

generate/commit_hook

@@ -1,6 +1,6 @@
 #/bin/bash
 [[ -d generate ]] && cd generate
-python3 test.py precision || exit $?
+ASTRONOMY_ENGINE_PYTHON_PROBLEM=1 python3 test.py precision || exit $?
 rm -f output/vsop*.txt output/*.eph output/jupiter_moons.txt
 ./run || exit $?
 ./verify_clean || exit $?

generate/commit_hook.bat

@@ -5,9 +5,10 @@ REM     Astronomy Engine - GitHub Actions steps for Windows.
 REM     This batch file is executed on every push to GitHub.
 REM --------------------------------------------------------------------------------
 
-REM **** HACK **** Just test the Python precision problem and quit.
 cd %~dp0
+set ASTRONOMY_ENGINE_PYTHON_PROBLEM=1
 py test.py precision || exit /b 1
+set ASTRONOMY_ENGINE_PYTHON_PROBLEM=
 
 REM Change to project/repo root directory.
 cd %~dp0\..

generate/template/astronomy.py

@@ -36,6 +36,7 @@
 import enum
 import re
 import abc
+from os import getenv
 from typing import Any, List, Tuple, Optional, Union, Callable, Dict
 
 def _cbrt(x: float) -> float:
@@ -1865,12 +1866,16 @@ def _VsopSphereToRect(lon: float, lat: float, rad: float) -> _TerseVector:
         rad * math.sin(lat)
     )
 
+_ProblemDebug:bool = ('1' == getenv('ASTRONOMY_ENGINE_PYTHON_PROBLEM'))
+
 def _CalcVsop(model: _vsop_model_t, time: Time) -> Vector:
     t = time.tt / _DAYS_PER_MILLENNIUM
     lon = _VsopFormula(model.lon, t, True)
     lat = _VsopFormula(model.lat, t, False)
     rad = _VsopFormula(model.rad, t, False)
     eclip = _VsopSphereToRect(lon, lat, rad)
+    if _ProblemDebug:
+        print('_CalcVsop: lon={:0.16g}, lat={:0.16g}, rad={:0.16g}, eclip=({:0.16g}, {:0.16g}, {:0.16g})'.format(lon, lat, rad, eclip.x, eclip.y, eclip.z))
     return _VsopRotate(eclip).ToAstroVector(time)
 
 class _body_state_t:

source/python/README.md

@@ -3567,3 +3567,12 @@ latitude, longitude, and elevation for that vector.
 The geographic latitude, longitude, and elevation above sea level
 that corresponds to the given equatorial vector.
 
+---
+
+<a name="getenv"></a>
+### getenv(key, default=None)
+
+Get an environment variable, return None if it doesn't exist.
+The optional second argument can specify an alternate default.
+key, default and the result are str.
+

source/python/astronomy/astronomy.py

@@ -36,6 +36,7 @@
 import enum
 import re
 import abc
+from os import getenv
 from typing import Any, List, Tuple, Optional, Union, Callable, Dict
 
 def _cbrt(x: float) -> float:
@@ -3180,12 +3181,16 @@ def _VsopSphereToRect(lon: float, lat: float, rad: float) -> _TerseVector:
         rad * math.sin(lat)
     )
 
+_ProblemDebug:bool = ('1' == getenv('ASTRONOMY_ENGINE_PYTHON_PROBLEM'))
+
 def _CalcVsop(model: _vsop_model_t, time: Time) -> Vector:
     t = time.tt / _DAYS_PER_MILLENNIUM
     lon = _VsopFormula(model.lon, t, True)
     lat = _VsopFormula(model.lat, t, False)
     rad = _VsopFormula(model.rad, t, False)
     eclip = _VsopSphereToRect(lon, lat, rad)
+    if _ProblemDebug:
+        print('_CalcVsop: lon={:0.16g}, lat={:0.16g}, rad={:0.16g}, eclip=({:0.16g}, {:0.16g}, {:0.16g})'.format(lon, lat, rad, eclip.x, eclip.y, eclip.z))
     return _VsopRotate(eclip).ToAstroVector(time)
 
 class _body_state_t: