Looking at _CalcVsop as the possible origin of numeric errors.

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 datetime
 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: