diff --git a/.gitattributes b/.gitattributes deleted file mode 100644 index b94f62d2..00000000 --- a/.gitattributes +++ /dev/null @@ -1,9 +0,0 @@ -*.png filter=lfs diff=lfs merge=lfs -text -*.jpg filter=lfs diff=lfs merge=lfs -text -*.egm* filter=lfs diff=lfs merge=lfs -text -gggrx_1200a_sha* filter=lfs diff=lfs merge=lfs -text -*.bpc filter=lfs diff=lfs merge=lfs -text -*.bsp filter=lfs diff=lfs merge=lfs -text -*.dbf filter=lfs diff=lfs merge=lfs -text -*.shp filter=lfs diff=lfs merge=lfs -text -*.shx filter=lfs diff=lfs merge=lfs -text diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index e9c9deea..c0d56216 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -23,7 +23,6 @@ concurrency: env: CI_PROJECT_DIR: ${{ github.workspace }} - GIT_LFS_SKIP_SMUDGE: 1 jobs: test: @@ -39,7 +38,6 @@ jobs: - uses: actions/checkout@v4 with: submodules: 'recursive' - lfs: false token: ${{ github.token }} - name: Set up Python ${{ matrix.python-version }} uses: actions/setup-python@v5 @@ -50,15 +48,11 @@ jobs: run: | # Install dependencies sudo apt-get update - sudo apt-get install build-essential git git-lfs python3 python3-setuptools python3-venv graphviz cmake + sudo apt-get install build-essential git python3 python3-setuptools python3-venv graphviz cmake # Initialize and update submodules git submodule update --init --recursive - # Install Git LFS and pull LFS objects - # git lfs install - # git lfs pull - # Upgrade pip and install Python dependencies python -m pip install --upgrade pip setuptools flake8 python -m pip install -r requirements.txt @@ -77,12 +71,6 @@ jobs: # Use the actions/setup-python toolchain and avoid Homebrew updates, # which can fail CI when the hosted runner's brew state is stale. python -m pip install --upgrade pip setuptools flake8 cmake - - # Make LFS usable if the runner already has it, but do not require - # Homebrew just for CI wheel/import checks. - if command -v git-lfs >/dev/null 2>&1; then - git lfs install --skip-repo - fi git config --global credential.helper "" # avoid Keychain prompts # Use the actions/setup-python toolchain consistently diff --git a/.github/workflows/python-publish.yml b/.github/workflows/python-publish.yml index db431709..e951b39d 100644 --- a/.github/workflows/python-publish.yml +++ b/.github/workflows/python-publish.yml @@ -25,9 +25,7 @@ jobs: - name: Install dependencies run: | sudo apt-get update - sudo apt-get install build-essential git git-lfs cmake - git lfs install - git lfs pull + sudo apt-get install build-essential git cmake python -m pip install --upgrade pip build twine python -m pip install -r requirements.txt diff --git a/MANIFEST.in b/MANIFEST.in index ab5c37cf..b5ce3fd2 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -11,4 +11,3 @@ recursive-include pybind11/include/ * include pybind11/CMakeLists.txt include pybind11/LICENSE recursive-include pybind11/tools * -recursive-include ssapy/data * diff --git a/README.rst b/README.rst index 2d938873..f8a614c5 100644 --- a/README.rst +++ b/README.rst @@ -29,6 +29,9 @@ wrappers around coordinate conversions, Lambertian magnitude / brightness calculations, and related extensions, see the companion project `SSAPy-Toolkit `__. +Optional packaged datasets are available with ``llnl-ssapy[data]``. Install +``llnl-ssapy[all-data]`` to install every split ``ssatk-data-*`` distribution. + SSAPy includes: - Ability to define satellite parameters (area, mass, radiation and drag coefficients, etc.) diff --git a/docs/source/api.rst b/docs/source/api.rst index f847b50b..b7e0087f 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -20,6 +20,7 @@ and/or heavy dependencies. Below lists the packages available in the ``ssapy`` ssapy.constants ssapy.correlate_tracks ssapy.ellipsoid + ssapy.ephemeris ssapy.gravity ssapy.io ssapy.linker diff --git a/docs/source/installation.rst b/docs/source/installation.rst index 15593cfe..a364e85b 100644 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -19,7 +19,7 @@ These requirements can be easily installed on most modern macOS and Linux system .. code-block:: console apt update - apt install build-essential git git-lfs python3 python3-distutils python3-venv graphviz + apt install build-essential git python3 python3-distutils python3-venv graphviz .. tab:: RHEL @@ -27,14 +27,14 @@ These requirements can be easily installed on most modern macOS and Linux system dnf install epel-release dnf group install "Development Tools" - dnf install git git-lfs gcc-gfortran python3 python3-pip python3-setuptools graphviz + dnf install git gcc-gfortran python3 python3-pip python3-setuptools graphviz .. tab:: macOS Brew .. code-block:: console brew update - brew install gcc git git-lfs python3 graphviz + brew install gcc git python3 graphviz Installation ------------ @@ -45,6 +45,30 @@ As the package has been published on `PyPI =0.2.0 astropy pyerfa emcee diff --git a/setup.py b/setup.py index d8cca2af..5e3da68e 100644 --- a/setup.py +++ b/setup.py @@ -45,7 +45,6 @@ def build_extension(self, ext): ext_modules=[CMakeExtension("ssapy._ssapy")], cmdclass={"build_ext": CMakeBuild}, packages=find_namespace_packages(include=["ssapy", "ssapy.*"]), - package_data={'ssapy': ['data/*']}, zip_safe=False, include_package_data=True, ) diff --git a/skills.md b/skills.md index 19e8ea36..80843591 100644 --- a/skills.md +++ b/skills.md @@ -115,11 +115,15 @@ r, v = rv(orbit, times, propagator=KeplerianPropagator()) ## Data and Binary Rules -SSAPy currently includes required package data under `ssapy/data`, and that data -is still part of the base package. Do not migrate or delete it unless the task -explicitly calls for a coordinated data migration. - -Use `ssapy.datadir` or `ssapy.utils.find_file()` when reading existing bundled +SSAPy's data (ephemerides, gravity models, textures) comes from the required +`llnl-ssapy-data` package, under its `ssapy/` tree; the repository no longer +stores data or uses Git LFS. Planetary ephemerides are handled by `ssapy.ephemeris`: +DE440 (`de440s.bsp`, 1849-2150) by default and a 1900-2150 DE430 excerpt on +request. Epochs outside those spans download the full-span kernel from NAIF +(checksummed, cached, with a warning); `SSAPY_EPHEMERIS_DOWNLOAD=0` disables +that. + +Use `ssapy.datadir` or `ssapy.utils.find_file()` when reading bundled resources. Some model loaders already resolve packaged gravity, SPICE, and Earth orientation files internally; prefer those public loaders over hard-coded paths. @@ -184,7 +188,7 @@ python -m twine check dist/* ## Common Pitfalls -- Do not assume all packaged resources are small; base SSAPy still carries existing required data. +- Do not add data files to this repository; add them to SSAPy-Data and release a new `llnl-ssapy-data`. - Do not add Git LFS as a default solution for new data in this repository. - Do not commit generated `_ssapy*.so`, `build/`, `dist/`, egg-info, docs build output, caches, or local environments. - Do not rely on unbuilt in-tree imports for code paths that require `ssapy._ssapy`. diff --git a/ssapy/__init__.py b/ssapy/__init__.py index 83d10fd0..d8ee4160 100644 --- a/ssapy/__init__.py +++ b/ssapy/__init__.py @@ -2,9 +2,23 @@ __version__ = "1.1.10" -datadir = os.path.join(os.path.dirname(__file__), "data") + + +def _datadir(): + """Directory holding SSAPy's data files: ssapy/ in llnl-ssapy-data.""" + try: + from importlib.resources import files + return os.fspath(files("ssapy_data") / "data" / "ssapy") + except (ImportError, TypeError): + # llnl-ssapy-data is a required dependency; without it find_file + # raises a FileNotFoundError that says how to install it. + return os.path.join(os.path.dirname(__file__), "data") + + +datadir = _datadir() from . import _ssapy +from . import ephemeris # noqa: F401 (ssapy.ephemeris.set_planetary_ephemeris) from .orbit import Orbit, EarthObserver, OrbitalObserver from .propagator import ( KeplerianPropagator, SeriesPropagator, RK4Propagator, SGP4Propagator, diff --git a/ssapy/accel.py b/ssapy/accel.py index 67e334d0..31b17814 100644 --- a/ssapy/accel.py +++ b/ssapy/accel.py @@ -365,12 +365,14 @@ class AccelDrag(Accel): ---------- recalc_threshold : float, optional Number of seconds past which the code will recompute the - precession/nutation matrix. Default: 86400*30 (30 days) + precession/nutation matrix. Default: 3600 (1 hour), matching + EarthOrientation; the previous 30-day default let the frame drift + 5 arcsec. defaultkw : dict default parameters for kwargs passed to __call__, (area, mass, CR) """ - def __init__(self, recalc_threshold=86400 * 30, **defaultkw): + def __init__(self, recalc_threshold=3600.0, **defaultkw): from . import _ssapy self.recalc_threshold = recalc_threshold @@ -412,7 +414,10 @@ def __call__(self, r, v, t, _T=None, **kwargs): self._t = t self._T = erfa.pnm80(2400000.5, mjd_tt) _T = self._T - r_sun = sunPos(t) + # The density model works in the true-of-date frame, so the Sun's + # right ascension and declination must be of date too; sunPos is + # GCRF, which put the diurnal bulge 0.37 deg off in 2026. + r_sun = _T @ sunPos(t) r_tod = _T @ r v_tod = _T @ v @@ -425,10 +430,9 @@ def __call__(self, r, v, t, _T=None, **kwargs): dec_sun ) if not np.isfinite(density): - print(f"r_tod = {r_tod}") - print(f"ra_sun = {ra_sun}") - print(f"dec_sun = {dec_sun}") - raise ValueError("non finite density") + raise ValueError( + f"non finite density at r_tod = {r_tod}, ra_sun = {ra_sun}, dec_sun = {dec_sun}" + ) a_tod = -0.5 * kw['CD'] * kw['area'] / kw['mass'] * density * v_rel * norm(v_rel) return _T.T @ a_tod diff --git a/ssapy/body.py b/ssapy/body.py index 756bce97..0cb0d779 100644 --- a/ssapy/body.py +++ b/ssapy/body.py @@ -15,7 +15,23 @@ def _close_if_possible(obj): close() +def _planetary_ephemeris_path(family=None): + """Path of the shipped kernel of the selected planetary ephemeris.""" + from .ephemeris import PlanetaryEphemeris + return PlanetaryEphemeris(family).shipped_path() + + class _KernelBacked: + _ephemeris = None + + def _positions(self, chain, mjd_tt): + """Sum of SPK segment chains at TT/TDB MJDs, in km, choosing a longer + kernel of the selected ephemeris for epochs the shipped one lacks.""" + if self._ephemeris is None: + from .ephemeris import PlanetaryEphemeris + self._ephemeris = PlanetaryEphemeris() + return self._ephemeris.compute(chain, 2400000.5, mjd_tt, self.kernel) + def close(self): kernel = getattr(self, "kernel", None) if kernel is not None: @@ -42,10 +58,14 @@ class EarthOrientation: Parameters ---------- recalc_threshold : float - Threshold for recomputing the orientation matrix. Default is 30 days. + Seconds after which the precession-nutation matrix and UT1 - TT are + recomputed; Earth rotation itself is evaluated at every call. Default + 1 hour, which bounds the frame error at 0.003 arcsec. The previous + 30-day default let the frame drift 1.0 arcsec in a week and 5.1 arcsec + in 29 days (172 m at 7000 km, 1 km at GEO). """ - def __init__(self, recalc_threshold=86400 * 30): + def __init__(self, recalc_threshold=3600.0): self.recalc_threshold = recalc_threshold self._t = None @@ -78,11 +98,10 @@ class MoonOrientation(_KernelBacked): orientation matrix at a given time. """ def __init__(self): - import os from jplephem.pck import PCK - from . import datadir + from .utils import find_file - fn = os.path.join(datadir, "moon_pa_de440_200625.bpc") + fn = find_file("moon_pa_de440_200625.bpc") self.kernel = PCK.open(fn) def __call__(self, t, _E=None): @@ -127,13 +146,17 @@ def _Rz(alpha): class MoonPosition(_KernelBacked): """Position of moon in GCRF. This is a callable class that returns the position vector at a given time. + + Uses the planetary ephemeris selected with + ``ssapy.ephemeris.set_planetary_ephemeris`` (DE440 by default), switching + to a longer kernel for epochs outside the shipped one. """ def __init__(self): - import os from jplephem.spk import SPK - from . import datadir + from .ephemeris import PlanetaryEphemeris - fn = os.path.join(datadir, "de430.bsp") # https://naif.jpl.nasa.gov/pub/naif/LUCY/kernels/spk/de430s.bsp.lbl + self._ephemeris = PlanetaryEphemeris() + fn = _planetary_ephemeris_path(self._ephemeris.family) self.kernel = SPK.open(fn) def __call__(self, t): @@ -150,8 +173,8 @@ def __call__(self, t): Position vector at time t in meters. """ mjd_tt = _gpsToTT(t) - pos = self.kernel[3, 301].compute(2400000.5, mjd_tt) # Earth-moon barycenter -> moon - pos -= self.kernel[3, 399].compute(2400000.5, mjd_tt) # Earth-moon barycenter -> earth + # Earth-moon barycenter -> moon, minus Earth-moon barycenter -> earth + pos = self._positions(((3, 301, 1), (3, 399, -1)), mjd_tt) return pos * 1e3 @@ -160,11 +183,11 @@ class SunPosition(_KernelBacked): position vector at a given time. """ def __init__(self): - import os from jplephem.spk import SPK - from . import datadir + from .ephemeris import PlanetaryEphemeris - fn = os.path.join(datadir, "de430.bsp") # https://naif.jpl.nasa.gov/pub/naif/LUCY/kernels/spk/de430s.bsp.lbl + self._ephemeris = PlanetaryEphemeris() + fn = _planetary_ephemeris_path(self._ephemeris.family) self.kernel = SPK.open(fn) def __call__(self, t): @@ -181,10 +204,8 @@ def __call__(self, t): Position vector at time t in meters. """ mjd_tt = _gpsToTT(t) - - pos = self.kernel[0, 10].compute(2400000.5, mjd_tt) # SS bary -> sun - pos -= self.kernel[0, 3].compute(2400000.5, mjd_tt) # SS bary -> Earth-moon bary - pos -= self.kernel[3, 399].compute(2400000.5, mjd_tt) # Earth-moon bary -> Earth + # SS bary -> sun, minus SS bary -> Earth-moon bary, minus Earth-moon bary -> Earth + pos = self._positions(((0, 10, 1), (0, 3, -1), (3, 399, -1)), mjd_tt) return pos * 1e3 @@ -193,11 +214,11 @@ class PlanetPosition(_KernelBacked): position vector at a given time. """ def __init__(self, planet_index): - import os from jplephem.spk import SPK - from . import datadir + from .ephemeris import PlanetaryEphemeris - fn = os.path.join(datadir, "de430.bsp") + self._ephemeris = PlanetaryEphemeris() + fn = _planetary_ephemeris_path(self._ephemeris.family) self.kernel = SPK.open(fn) self.planet_index = planet_index @@ -215,9 +236,9 @@ def __call__(self, t): Position vector at time t in meters. """ mjd_tt = _gpsToTT(t) - pos = self.kernel[0, self.planet_index].compute(2400000.5, mjd_tt) # SS bary -> Jupiter - pos -= self.kernel[0, 3].compute(2400000.5, mjd_tt) # SS bary -> Earth-moon bary - pos -= self.kernel[3, 399].compute(2400000.5, mjd_tt) # Earth-moon bary -> Earth + # SS bary -> planet barycenter, minus SS bary -> Earth-moon bary, minus Earth-moon bary -> Earth + chain = ((0, self.planet_index, 1), (0, 3, -1), (3, 399, -1)) + pos = self._positions(chain, mjd_tt) return pos * 1e3 diff --git a/ssapy/compute.py b/ssapy/compute.py index 523eb415..d70fea5a 100644 --- a/ssapy/compute.py +++ b/ssapy/compute.py @@ -853,18 +853,15 @@ def radecRate( Returns ------- - ra : array_like (n, m) - right ascension in radians raRate : array_like (n, m) Rate of change of right ascension*cos(dec) in radians per second. - dec : array_link (n, m) - declination in radians decRate : array_like (n, m) Rate of change of declination in radians per second. - slantRange : array_like (n, m) - Range in meters slantRangeRate : array_like (n, m) Slant range rate in meters per second. + + Only the three rates are returned (the docstring previously listed six + values); use ``radec(..., rate=True)`` for angles, range and rates. """ import warnings warnings.warn("This function is deprecated; use " diff --git a/ssapy/constants.py b/ssapy/constants.py index 561524f1..7d656870 100644 --- a/ssapy/constants.py +++ b/ssapy/constants.py @@ -237,7 +237,7 @@ SUN_MASS = 1.98847e+30 MOON_MASS = 7.348e22 MERCURY_MASS = 3.301e23 -VENUS_MASS = 4.687e24 +VENUS_MASS = 4.8675e24 EARTH_MASS = 5.9722e24 MARS_MASS = 6.417e23 JUPITER_MASS = 1.899e27 diff --git a/ssapy/correlate_tracks.py b/ssapy/correlate_tracks.py index fec1ffee..203eb47b 100644 --- a/ssapy/correlate_tracks.py +++ b/ssapy/correlate_tracks.py @@ -222,7 +222,9 @@ class VolumeDistancePrior: """ def __init__(self, scale=RGEO): - self.scale = RGEO + # The argument used to be ignored (self.scale = RGEO), so every prior + # peaked at 2 RGEO whatever scale was requested. + self.scale = scale def __call__(self, orbit, distance, chi=False): """Return log prior probability of given orbit. diff --git a/ssapy/data/Earth_graphics/ne_50m_ocean.README.html b/ssapy/data/Earth_graphics/ne_50m_ocean.README.html deleted file mode 100644 index 9f9b1de5..00000000 --- a/ssapy/data/Earth_graphics/ne_50m_ocean.README.html +++ /dev/null @@ -1,395 +0,0 @@ - - - - - - - - -Ocean | Natural Earth - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - -
-
- - - « 1:50m Physical Vectors  - « Downloads  -
-

Ocean

-
-
-
water_thumb
-
Ocean polygon split into contiguous pieces.

-
- Download ocean (450.91 KB) version 4.0.0 -
-

-
-
About

-

Derives from 10m ocean coastlines. Ocean polygon is the same as Bathymetry 0 (sea level).

-

(below) Yucatan peninsula, Cuba, and Hispaniola.

-

water_banner

-

Issues

-

Caspian’s Garabogaz Aylagy and other “lagoons” might better belong in the lake theme than oceans.

-

Dependencies

-

This theme (and it’s inverse cousin, Land) are used to clip the raster imagery.

-

Version History

-
    -
  • - 4.0.0 -
    • -
  • - 2.0.0 -
    • -
  • - 1.3.0 -
    • -
  • - 1.1.0 -
    • -
  • - 1.0.0 -
    • -
- -

The master changelog is available on Github » -

- - -
- -
- - -
- - - - -
-
- - - - - - \ No newline at end of file diff --git a/ssapy/data/Earth_graphics/ne_50m_ocean.VERSION.txt b/ssapy/data/Earth_graphics/ne_50m_ocean.VERSION.txt deleted file mode 100644 index 5dbaad6d..00000000 --- a/ssapy/data/Earth_graphics/ne_50m_ocean.VERSION.txt +++ /dev/null @@ -1 +0,0 @@ -4.1.0 diff --git a/ssapy/data/Earth_graphics/ne_50m_ocean.cpg b/ssapy/data/Earth_graphics/ne_50m_ocean.cpg deleted file mode 100644 index 3ad133c0..00000000 --- a/ssapy/data/Earth_graphics/ne_50m_ocean.cpg +++ /dev/null @@ -1 +0,0 @@ -UTF-8 \ No newline at end of file diff --git a/ssapy/data/Earth_graphics/ne_50m_ocean.dbf b/ssapy/data/Earth_graphics/ne_50m_ocean.dbf deleted file mode 100644 index cb8ca214..00000000 --- a/ssapy/data/Earth_graphics/ne_50m_ocean.dbf +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:162d515911b17e35d36e22c5757e75fd37bdd38db2bc0011b9c79888a28bb4f5 -size 176 diff --git a/ssapy/data/Earth_graphics/ne_50m_ocean.prj b/ssapy/data/Earth_graphics/ne_50m_ocean.prj deleted file mode 100644 index a30c00a5..00000000 --- a/ssapy/data/Earth_graphics/ne_50m_ocean.prj +++ /dev/null @@ -1 +0,0 @@ -GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] \ No newline at end of file diff --git a/ssapy/data/Earth_graphics/ne_50m_ocean.shp b/ssapy/data/Earth_graphics/ne_50m_ocean.shp deleted file mode 100644 index 014215f5..00000000 --- a/ssapy/data/Earth_graphics/ne_50m_ocean.shp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:6ff1190e8430e493cf8eed99393f021b7562b600b9a7a11aa23ce45b6035d8b8 -size 978340 diff --git a/ssapy/data/Earth_graphics/ne_50m_ocean.shx b/ssapy/data/Earth_graphics/ne_50m_ocean.shx deleted file mode 100644 index 099819b2..00000000 --- a/ssapy/data/Earth_graphics/ne_50m_ocean.shx +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:449353f638399c48c4cfe59e8f945e5c3025a4e412e9bcd0eeba64ba137ae7c1 -size 108 diff --git a/ssapy/data/de430.bsp b/ssapy/data/de430.bsp deleted file mode 100644 index 533d4b74..00000000 --- a/ssapy/data/de430.bsp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:6e1b277c5f07135a84950604b83e56b736be696a7f3560bcddb1d4aeb944fca1 -size 119741440 diff --git a/ssapy/data/earth.png b/ssapy/data/earth.png deleted file mode 100644 index 16d9db60..00000000 --- a/ssapy/data/earth.png +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:6e8a1296caef991dcdba963fdc72816e4095720d005e7af6688b060e5fb85aa4 -size 8825662 diff --git a/ssapy/data/egm2008.egm b/ssapy/data/egm2008.egm deleted file mode 100644 index 9442ced7..00000000 --- a/ssapy/data/egm2008.egm +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:031ecae647986faadf3f3a0177ba07ffcb73a2504b084c2e1c82b97cd9e61c75 -size 790 diff --git a/ssapy/data/egm2008.egm.cof b/ssapy/data/egm2008.egm.cof deleted file mode 100644 index 4268c37d..00000000 --- a/ssapy/data/egm2008.egm.cof +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:ebed5b09e3da7aeae37ac49e7b8c803fda0a525462bba374f762d6345624d132 -size 75755304 diff --git a/ssapy/data/egm84.egm b/ssapy/data/egm84.egm deleted file mode 100644 index a6512e8b..00000000 --- a/ssapy/data/egm84.egm +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:549e21625b71a0f0daec8477e0f4928338885ed2269ff70e351b15c3e6e0d43e -size 800 diff --git a/ssapy/data/egm84.egm.cof b/ssapy/data/egm84.egm.cof deleted file mode 100644 index 2728f31d..00000000 --- a/ssapy/data/egm84.egm.cof +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:f95fa2967091b9a0630ae7f34934446f17a02f73e5b648ea919db4174b047633 -size 262112 diff --git a/ssapy/data/egm96.egm b/ssapy/data/egm96.egm deleted file mode 100644 index d6f2e3b7..00000000 --- a/ssapy/data/egm96.egm +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:cdf0896bb34c1f8ad7ea1c4126cf949fe5069ad9df17c2a66cf3d81da096970e -size 753 diff --git a/ssapy/data/egm96.egm.cof b/ssapy/data/egm96.egm.cof deleted file mode 100644 index 8323c719..00000000 --- a/ssapy/data/egm96.egm.cof +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:fd427a88a944e2df46d3a3e1aefdb17a68e892e3decdf1382f3ca14de9b6f2ae -size 2085160 diff --git a/ssapy/data/gggrx_1200a_sha.lbl b/ssapy/data/gggrx_1200a_sha.lbl deleted file mode 100644 index bb784ae3..00000000 --- a/ssapy/data/gggrx_1200a_sha.lbl +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:8d46351c2e73db699d4669a98da1f47551086cbee321d8f9baa3fdb6b65d3d5f -size 10637 diff --git a/ssapy/data/gggrx_1200a_sha.tab b/ssapy/data/gggrx_1200a_sha.tab deleted file mode 100644 index 0d404d19..00000000 --- a/ssapy/data/gggrx_1200a_sha.tab +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:fa04c3dce9376948ad243f3df74144e2602f12d183ea4d179604ed0a79da7ded -size 88059844 diff --git a/ssapy/data/moon.png b/ssapy/data/moon.png deleted file mode 100644 index 754fbf09..00000000 --- a/ssapy/data/moon.png +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:4d9e6eebd180262aefe7555a0231927e5142954144cfcf2c4a10876b7c218027 -size 13728894 diff --git a/ssapy/data/moon_pa_de440_200625.bpc b/ssapy/data/moon_pa_de440_200625.bpc deleted file mode 100644 index ed2f32ac..00000000 --- a/ssapy/data/moon_pa_de440_200625.bpc +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:60cd55aa401ea2ea97360636f567554bfe4e37bb829f901b4460a455dfaf783f -size 12863488 diff --git a/ssapy/data/wgs84.egm b/ssapy/data/wgs84.egm deleted file mode 100644 index c27532ed..00000000 --- a/ssapy/data/wgs84.egm +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:594ff1c9dc2363e50d663003deb9baca2e144d096d70fc9b7a6c1c86f5b781f0 -size 724 diff --git a/ssapy/data/wgs84.egm.cof b/ssapy/data/wgs84.egm.cof deleted file mode 100644 index 0734d247..00000000 --- a/ssapy/data/wgs84.egm.cof +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:aeedf849678826d171b731eb7a783ab5d7e47552413da551418d21445fc3a2ec -size 192 diff --git a/ssapy/ephemeris.py b/ssapy/ephemeris.py new file mode 100644 index 00000000..ab5ec8b6 --- /dev/null +++ b/ssapy/ephemeris.py @@ -0,0 +1,292 @@ +""" +Planetary ephemeris selection and long-span kernels. + +SSAPy evaluates Sun, Moon and planet positions from JPL development +ephemerides (DE). ``llnl-ssapy-data`` ships two short kernels: + +========== ======================== ========================= +Ephemeris Shipped kernel Coverage +========== ======================== ========================= +``de440`` ``de440s.bsp`` 1849-12-26 to 2150-01-22 +``de430`` ``de430_1900_2150.bsp`` 1900-01-01 to 2150-01-01 +========== ======================== ========================= + +DE440 is the default. Select DE430, for example to reproduce results made +with SSAPy <= 1.1.10, with :func:`set_planetary_ephemeris` or the +``SSAPY_EPHEMERIS`` environment variable. + +For epochs outside the shipped kernel, SSAPy uses a longer kernel of the same +ephemeris (``de440.bsp`` or ``de430.bsp``, 1550-2650) and, for DE440 beyond +that, DE441 (-13200 to +17191). These are 120 MB to 1.65 GB, so they are not +shipped: SSAPy looks for them in the working directory, ``llnl-ssapy-data`` +and the cache directory, and otherwise downloads them from NAIF, checks their +SHA-256, and caches them, warning with :class:`EphemerisDownloadWarning` +before each download. Set ``SSAPY_EPHEMERIS_DOWNLOAD=0`` to forbid downloads +(e.g. on compute nodes without network access; prefetch on a login node with +:func:`fetch`) and ``SSAPY_DATA_CACHE`` to choose the cache directory. + +The shipped and long kernels of one ephemeris hold the same solution, so +positions are continuous where SSAPy switches between them. DE441 is a +different solution from DE440 (no lunar core-mantle damping; less accurate +in the current century), so using it raises :class:`EphemerisRangeWarning`. +""" +from __future__ import annotations + +import hashlib +import os +import shutil +import tempfile +import warnings +from dataclasses import dataclass +from typing import Optional + +import numpy as np + +__all__ = [ + "EphemerisDownloadWarning", + "EphemerisRangeWarning", + "EphemerisUnavailableError", + "set_planetary_ephemeris", + "get_planetary_ephemeris", + "fetch", + "cache_dir", +] + +NAIF_PLANETS_URL = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/" + + +class EphemerisDownloadWarning(UserWarning): + """A long-span ephemeris kernel is being downloaded.""" + + +class EphemerisRangeWarning(UserWarning): + """An epoch is evaluated with a different solution than the one selected.""" + + +class EphemerisUnavailableError(FileNotFoundError): + """No kernel of the selected ephemeris covers the requested epoch.""" + + +@dataclass(frozen=True) +class Kernel: + name: str + start_jd: float # TDB Julian dates covered by every segment SSAPy uses + end_jd: float + sha256: Optional[str] = None # None for kernels shipped in llnl-ssapy-data + nbytes: Optional[int] = None + solution: str = "" + + @property + def url(self): + return NAIF_PLANETS_URL + self.name + + def covers(self, jd): + return (jd >= self.start_jd) & (jd <= self.end_jd) + + +EPHEMERIDES = { + "de440": ( + Kernel("de440s.bsp", 2396752.5, 2506352.5, solution="DE440"), + Kernel("de440.bsp", 2287184.5, 2688976.5, + "a4ce9bf9b3282becc9f4b2ac3cebe03a2ae7599981aabd7265fd8482fff7c4b5", 119799808, "DE440"), + Kernel("de441_part-2.bsp", 2440400.5, 8000016.5, + "3abb17dae2d78dd34880377544aacb54892104a0d4462b322cb9f4454d4887f6", 1656830976, "DE441"), + Kernel("de441_part-1.bsp", -3100015.5, 2440432.5, + "13757827f5db41b835a24bbd637488636ce79a8ca754062fed17844f7d5b618e", 1651119104, "DE441"), + ), + "de430": ( + Kernel("de430_1900_2150.bsp", 2415020.5, 2506331.5, solution="DE430"), + Kernel("de430.bsp", 2287184.5, 2688976.5, + "6e1b277c5f07135a84950604b83e56b736be696a7f3560bcddb1d4aeb944fca1", 119741440, "DE430"), + ), +} + +_selected = None +_open_kernels = {} + + +def set_planetary_ephemeris(name): + """Select the planetary ephemeris, ``"de440"`` (default) or ``"de430"``. + + Affects Sun, Moon and planet positions created afterwards + (``get_body``, ``SunPosition``, ``MoonPosition``, ``PlanetPosition``). + """ + global _selected + key = str(name).lower() + if key not in EPHEMERIDES: + raise ValueError(f"Unknown planetary ephemeris {name!r}; choose from {sorted(EPHEMERIDES)}") + _selected = key + + +def get_planetary_ephemeris(): + """Name of the selected planetary ephemeris.""" + if _selected is not None: + return _selected + key = os.environ.get("SSAPY_EPHEMERIS", "de440").lower() + if key not in EPHEMERIDES: + raise ValueError(f"SSAPY_EPHEMERIS={key!r} is not one of {sorted(EPHEMERIDES)}") + return key + + +def cache_dir(): + """Directory where downloaded kernels are kept (``SSAPY_DATA_CACHE``).""" + return os.path.expanduser(os.environ.get("SSAPY_DATA_CACHE", os.path.join("~", ".cache", "ssapy"))) + + +def _downloads_allowed(): + return os.environ.get("SSAPY_EPHEMERIS_DOWNLOAD", "1").strip().lower() not in ("0", "false", "no", "off") + + +def _sha256(path): + digest = hashlib.sha256() + with open(path, "rb") as handle: + for chunk in iter(lambda: handle.read(1 << 20), b""): + digest.update(chunk) + return digest.hexdigest() + + +def _open_url(url): # separated so tests can replace the network + from urllib.request import urlopen + return urlopen(url, timeout=60) + + +def _download(kernel, destination): + os.makedirs(os.path.dirname(destination), exist_ok=True) + handle, partial = tempfile.mkstemp(dir=os.path.dirname(destination), suffix=".part") + try: + with os.fdopen(handle, "wb") as out, _open_url(kernel.url) as response: + shutil.copyfileobj(response, out, 1 << 20) + if os.path.getsize(partial) != kernel.nbytes or _sha256(partial) != kernel.sha256: + raise EphemerisUnavailableError( + f"{kernel.name} downloaded from {kernel.url} does not match its recorded SHA-256" + ) + os.replace(partial, destination) + finally: + if os.path.exists(partial): + os.remove(partial) + + +def _local_path(kernel): + from .utils import find_file + try: + return find_file(kernel.name) + except FileNotFoundError: + pass + cached = os.path.join(cache_dir(), kernel.name) + return cached if os.path.isfile(cached) else None + + +def fetch(name=None, include_de441=True): + """Download (if needed) the long-span kernels of an ephemeris and return + their paths. Use on a login node to fill ``SSAPY_DATA_CACHE`` before + running jobs without network access.""" + family = EPHEMERIDES[(name or get_planetary_ephemeris()).lower()] + return [ + _path(kernel, "requested with ssapy.ephemeris.fetch") + for kernel in family[1:] + if include_de441 or kernel.solution == family[0].solution + ] + + +def _path(kernel, reason): + path = _local_path(kernel) + if path is not None: + return path + if kernel.sha256 is None: + raise EphemerisUnavailableError( + f"{kernel.name} is missing; install 'llnl-ssapy-data>=0.2.0'." + ) + destination = os.path.join(cache_dir(), kernel.name) + manual = (f"Download {kernel.url} ({kernel.nbytes / 1e6:.0f} MB) into {cache_dir()} " + f"(or set SSAPY_DATA_CACHE), or place it in the working directory.") + if not _downloads_allowed(): + raise EphemerisUnavailableError(f"{kernel.name} is needed ({reason}) and downloads are disabled " + f"(SSAPY_EPHEMERIS_DOWNLOAD=0). {manual}") + warnings.warn( + f"Downloading {kernel.name} ({kernel.nbytes / 1e6:.0f} MB) from NAIF to {destination}: {reason}.", + EphemerisDownloadWarning, + stacklevel=6, + ) + try: + _download(kernel, destination) + except EphemerisUnavailableError: + raise + except Exception as exc: # no network, proxy, HTTP error, disk full + raise EphemerisUnavailableError(f"Could not download {kernel.name}: {exc}. {manual}") from exc + return destination + + +def _open(kernel, path): + from jplephem.spk import SPK + key = os.path.abspath(path) + if key not in _open_kernels: + _open_kernels[key] = SPK.open(path) + return _open_kernels[key] + + +def _jd_text(jd): + from astropy.time import Time + return Time(float(jd), format="jd", scale="tdb").iso[:10] + + +class PlanetaryEphemeris: + """Evaluate SPK segment chains, choosing a kernel per epoch. + + The shipped kernel of ``family`` (passed to :meth:`compute` as + ``primary`` once opened) is used for every epoch it covers; other epochs + go to the first longer kernel of the family that covers them. + """ + + def __init__(self, family=None): + self.family = (family or get_planetary_ephemeris()).lower() + self.kernels = EPHEMERIDES[self.family] + + def _kernel_for(self, kernel, reason, primary): + if kernel is self.kernels[0] and primary is not None: + return primary + return _open(kernel, _path(kernel, reason)) + + def shipped_path(self): + return _path(self.kernels[0], "shipped kernel") + + def compute(self, chain, jd1, jd2, primary=None): + """Sum ``sign * kernel[center, target].compute(jd1, jd2)`` over the + ``(center, target, sign)`` terms in ``chain``; returns km.""" + jd2 = np.asarray(jd2, dtype=float) + shipped = self.kernels[0] + if primary is not None and np.all(shipped.covers(jd1 + jd2)): + # Fast path (every call in normal use): one kernel, no masking. + total = 0.0 + for center, target, sign in chain: + total = total + sign * np.asarray(primary[center, target].compute(jd1, jd2)) + return total + scalar = jd2.ndim == 0 + jd2 = np.atleast_1d(jd2) + jd = jd1 + jd2 + out = np.zeros((3, jd2.size)) + pending = np.ones(jd2.size, dtype=bool) + for kernel in self.kernels: + mask = pending & kernel.covers(jd) + if not np.any(mask): + continue + reason = None + if kernel is not shipped: + reason = (f"epochs {_jd_text(jd[mask].min())} to {_jd_text(jd[mask].max())} are outside " + f"{shipped.name} ({_jd_text(shipped.start_jd)} to {_jd_text(shipped.end_jd)})") + spk = self._kernel_for(kernel, reason, primary) + if kernel.solution != self.kernels[0].solution: + warnings.warn( + f"Using {kernel.solution} ({kernel.name}) for {reason}; {kernel.solution} differs from " + f"{self.kernels[0].solution}, and the switch makes positions jump slightly at its boundary.", + EphemerisRangeWarning, + stacklevel=4, + ) + for center, target, sign in chain: + out[:, mask] += sign * np.asarray(spk[center, target].compute(jd1, jd2[mask])).reshape(3, -1) + pending &= ~mask + if np.any(pending): + raise EphemerisUnavailableError( + f"No {self.family.upper()} kernel covers {_jd_text(jd[pending].min())} to " + f"{_jd_text(jd[pending].max())}." + ) + return out[:, 0] if scalar else out diff --git a/ssapy/gravity.py b/ssapy/gravity.py index 06d9fe98..89c3aa9a 100644 --- a/ssapy/gravity.py +++ b/ssapy/gravity.py @@ -5,6 +5,7 @@ # from functools import lru_cache import numpy as np +import warnings from .accel import Accel as _Accel, _invnorm from .utils import find_file, norm @@ -279,8 +280,8 @@ def __init__(self, body, n_max=None, m_max=None): m_max = body.harmonics.m_max self.body = body if n_max > body.harmonics.n_max: - print( - f"WARNING::The provided degree ({n_max}) is higher than the maximum value allowed by the " + warnings.warn( + f"The provided degree ({n_max}) is higher than the maximum value allowed by the " f"{body.harmonics.name} model ({body.harmonics.n_max}).\nSetting the degree to " f"{body.harmonics.n_max}." ) @@ -288,8 +289,8 @@ def __init__(self, body, n_max=None, m_max=None): else: self.n_max = n_max if m_max > body.harmonics.m_max: - print( - f"WARNING::The provided order ({m_max}) is higher than the maximum value allowed by the " + warnings.warn( + f"The provided order ({m_max}) is higher than the maximum value allowed by the " f"{body.harmonics.name} model ({body.harmonics.m_max}).\nSetting the order to " f"{body.harmonics.m_max}." ) @@ -297,13 +298,17 @@ def __init__(self, body, n_max=None, m_max=None): else: self.m_max = m_max if self.m_max > self.n_max: - print( - f"WARNING::The provided order ({self.m_max}) is higher than the selected degree " + warnings.warn( + f"The provided order ({self.m_max}) is higher than the selected degree " f"({self.n_max}).\nSetting the order to {self.n_max}." ) self.m_max = self.n_max + # The coefficients are normalised to the model's own GM and reference + # radius, so both scale the expansion. Using the body's central GM + # instead put GRGM1200A lunar accelerations 2.2e-7 low (DE200 Moon GM) + # and EGM2008 Earth accelerations 7.5e-10 high. self._harmonic = _ssapy.AccelHarmonic( - self.body.mu, + getattr(self.body.harmonics, "MG", None) or self.body.mu, self.body.harmonics.radius, self.body.harmonics.CS.shape[0], self.body.harmonics.CS.ctypes.data diff --git a/ssapy/io.py b/ssapy/io.py index c52841a0..3b076ca6 100644 --- a/ssapy/io.py +++ b/ssapy/io.py @@ -202,7 +202,7 @@ def parse_tle(tle): return a, e, i, pa, raan, trueAnomaly, epoch.gps -def make_tle(a, e, i, pa, raan, trueAnomaly, t): +def make_tle(a, e, i, pa, raan, trueAnomaly, t, drag_fields=None): """ Create a TLE from Kozai mean orbital elements @@ -223,11 +223,16 @@ def make_tle(a, e, i, pa, raan, trueAnomaly, t): t : float or astropy.time.Time If float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC + drag_fields : str, optional + Line-1 columns 34-61 (ndot/2, nddot/6 and B*, 28 characters) to copy + verbatim, e.g. ``line1[33:61]`` of the TLE an orbit came from. By default + all three terms are written as zero. Notes ----- - Dynamic TLE terms, including the drag coefficient and ballistic coefficient, - are ignored in this function. + Unless ``drag_fields`` is given, the dynamic TLE terms are zero. They used + to be written with B* = '99999-0', which SGP4 reads as 0.99999 per Earth + radius: an enormous drag term for LEO (11,900 km of drift in 3 days). """ from .orbit import ( _ellipticalEccentricToMeanAnomaly, @@ -245,7 +250,12 @@ def make_tle(a, e, i, pa, raan, trueAnomaly, t): line1 += "{:02d}".format(int(year) % 100) line1 += "{:012.8f}".format(day) - line1 += " +.00000000 +00000-0 99999-0 0 0000" + # Columns 34-61: ndot/2, nddot/6 and B*, zero unless supplied. + if drag_fields is None: + drag_fields = "+.00000000 +00000-0 +00000-0" + if len(drag_fields) != 28: + raise ValueError("drag_fields must be the 28 characters of TLE line-1 columns 34-61.") + line1 += " " + drag_fields + " 0 0000" def checksum(s): check = 0 diff --git a/ssapy/orbit.py b/ssapy/orbit.py index dcaf671a..be2b51ff 100644 --- a/ssapy/orbit.py +++ b/ssapy/orbit.py @@ -1846,7 +1846,9 @@ def tle(self): """Two line element for this Orbit. """ from .io import make_tle - return make_tle(*self.kozaiMeanKeplerianElements, self.t) + source = getattr(self, "_tle", None) + drag_fields = source[0][33:61] if source is not None else None + return make_tle(*self.kozaiMeanKeplerianElements, self.t, drag_fields=drag_fields) class EarthObserver: diff --git a/ssapy/orbit_solver.py b/ssapy/orbit_solver.py index 682a2bd5..bc66ab9c 100644 --- a/ssapy/orbit_solver.py +++ b/ssapy/orbit_solver.py @@ -9,6 +9,7 @@ import abc import numpy as np +import warnings from astropy.time import Time from .orbit import Orbit @@ -132,6 +133,12 @@ class GaussTwoPosOrbitSolver(TwoPosOrbitSolver): and implements a method to compute the orbital parameter `p` based on Shefer's equations. + Gauss's fixed-point iteration converges only for short arcs: it recovers + the velocity to 1e-11 m/s for a 38 deg LEO arc but diverges for sweeps + near 90 deg and beyond (2.7 km/s error for a quarter GEO orbit). It warns + with a RuntimeWarning when it does not converge; SheferTwoPosOrbitSolver + handles long arcs. + Attributes: eps (float): Convergence tolerance for iterative calculations. maxiter (int): Maximum number of iterations allowed for convergence. @@ -167,10 +174,20 @@ def _getP(self): h = 2 * np.arcsinh(np.sqrt(-x)) X = (np.sinh(2 * h) - 2 * h) / np.sinh(h)**3 else: # -1e-15 1e-10 * max(1.0, np.abs(eta)): + warnings.warn( + "Gauss two-position iteration did not converge (|d_eta| = {:.3g} after {} iterations); " + "the solution is unreliable. Use SheferTwoPosOrbitSolver for long arcs.".format(abs(d_eta), niter), + RuntimeWarning, + stacklevel=3, + ) # Plug eta into (2) to obtain p p = (0.5 * eta * self.kappa * self.sigma / self.tau)**2 / self.mu return p @@ -225,43 +242,83 @@ def dXdg(g): """Compute dX(g)/dg from Shefer (12).""" return (2 * (1 - np.cos(2 * g)) - 3 * (2 * g - np.sin(2 * g)) / np.tan(g)) / (np.sin(g)**3) + def _eta_iteration(self): + """Newton iteration on eta (Shefer 14); returns eta.""" + with np.errstate(divide="ignore", invalid="ignore"): + eta = 0.5 * (np.sqrt(np.float64(self.m) / (self.ell + 1)) + np.sqrt(np.float64(self.m) / self.ell)) + if not np.isfinite(eta): + raise RuntimeError("Invalid x") + Geta = 1 + niter = 0 + while np.abs(Geta) > self.eps and niter < self.maxiter: + x = self.m / eta**2 - self.ell # Shefer (14) + if not np.isfinite(x) or (1 - 2 * x) > 1 or (1 - 2 * x) < -1: + raise RuntimeError("Invalid x") + g = np.arccos(1 - 2 * x) # Shefer (9) + dgdx = 2 / np.sin(g) # Shefer (10) + dxdeta = -2 * self.m / eta**3 # next few lines are Shefer (14ish) + dXdeta = self.dXdg(g) * dgdx * dxdeta + Geta = eta - 1 - (self.ell + x) * self.X(g) + dGdeta = 1 - self.X(g) * dxdeta - (self.ell + x) * dXdeta + eta -= Geta / dGdeta + niter += 1 + if np.abs(Geta) > 1e-9 * max(1.0, np.abs(eta)): + raise RuntimeError("Danchick eta iteration did not converge") + return eta + + def _x_iteration(self): + """Damped Newton iteration on x in (0, 1) (Shefer 15); returns eta.""" + x = 0.5 + Fx = 1 + niter = 0 + while np.abs(Fx) > self.eps and niter < self.maxiter: + g = np.arccos(1 - 2 * x) # Shefer (9) + dgdx = 2 / np.sin(g) # Shefer (10) + eta = 1 + (self.ell + x) * self.X(g) # Shefer (15) + detadx = (self.ell + x) * self.dXdg(g) * dgdx + self.X(g) + Fx = x - self.m / eta**2 + self.ell # Shefer (15ish) + dFdx = 1 + (2 * self.m / eta**3) * detadx + step = Fx / dFdx + if not np.isfinite(step): + raise RuntimeError("Invalid x") + # Halve the step until x stays strictly inside (0, 1). + for _ in range(60): + if 0.0 < x - step < 1.0: + break + step *= 0.5 + else: + raise RuntimeError("Invalid x") + x -= step + niter += 1 + if np.abs(Fx) > 1e-9: + raise RuntimeError("Danchick x iteration did not converge") + return 1 + (self.ell + x) * self.X(np.arccos(1 - 2 * x)) + def _getP(self): - if self.cos2f >= 0: - # TODO: rewrite to use newton_raphson function. - eta = 0.5 * (np.sqrt(self.m / (self.ell + 1)) + np.sqrt(self.m / self.ell)) - Geta = 1 - niter = 0 - while np.abs(Geta) > self.eps and niter < self.maxiter: - x = self.m / eta**2 - self.ell # Shefer (14) - if (1 - 2 * x) > 1 or (1 - 2 * x) < -1: - raise RuntimeError("Invalid x") - g = np.arccos(1 - 2 * x) # Shefer (9) - dgdx = 2 / np.sin(g) # Shefer (10) - dxdeta = -2 * self.m / eta**3 # next few lines are Shefer (14ish) - dXdeta = self.dXdg(g) * dgdx * dxdeta - Geta = eta - 1 - (self.ell + x) * self.X(g) - dGdeta = 1 - self.X(g) * dxdeta - (self.ell + x) * dXdeta - eta -= Geta / dGdeta - niter += 1 - x = self.m / eta**2 - self.ell # Shefer (9) - elif self.cos2f < 0: - x = 0.5 - Fx = 1 - niter = 0 - while np.abs(Fx) > self.eps and niter < self.maxiter: - if 1 - 2 * x > 1 or 1 - 2 * x < -1: - raise RuntimeError("Invalid x") - g = np.arccos(1 - 2 * x) # Shefer (9) - dgdx = 2 / np.sin(g) # Shefer (10) - eta = 1 + (self.ell + x) * self.X(g) # Shefer (15) - detadx = (self.ell + x) * self.dXdg(g) * dgdx + self.X(g) - Fx = x - self.m / eta**2 + self.ell # Shefer (15ish) - dFdx = 1 + (2 * self.m / eta**3) * detadx - x -= Fx / dFdx - niter += 1 - eta = 1 + (self.ell + x) * self.X(np.arccos(1 - 2 * x)) - else: + # Danchick's eta iteration suits short-way arcs with cos(2f) >= 0 and + # the x iteration the rest. Long-way arcs beyond 270 deg also have + # cos(2f) >= 0 but negative kappa, and fast perigee passages can + # have cos(2f) < 0 on short arcs where the x iteration's plain Newton + # step leaves 0 < x < 1; both used to raise "Invalid x". The + # preferred iteration is tried first, then the other one, and the x + # iteration damps any step that would leave its domain. + if not np.isfinite(self.cos2f): raise ValueError("Invalid value of cos2f") + order = (self._eta_iteration, self._x_iteration) + if not (self.cos2f >= 0 and self.kappa > 0): + order = order[::-1] + errors = [] + for iteration in order: + try: + eta = iteration() + except RuntimeError as exc: + errors.append(str(exc)) + continue + if np.isfinite(eta): + break + errors.append("non-finite eta") + else: + raise RuntimeError("Danchick two-position iteration failed: " + "; ".join(errors)) # Shefer (2) p = (0.5 * eta * self.kappa * self.sigma / self.tau)**2 / self.mu diff --git a/ssapy/propagator.py b/ssapy/propagator.py index 574176bf..721298f2 100644 --- a/ssapy/propagator.py +++ b/ssapy/propagator.py @@ -311,7 +311,11 @@ def _getRVOne(self, orbit, time): sat_epoch = orbit.t if self.truncate: - line1, line2 = make_tle(*orbit.kozaiMeanKeplerianElements, orbit.t) + # Keep the drag terms of the TLE the orbit came from, so truncation + # changes only the precision of the mean elements. + source = getattr(orbit, "_tle", None) + drag_fields = source[0][33:61] if source is not None else None + line1, line2 = make_tle(*orbit.kozaiMeanKeplerianElements, orbit.t, drag_fields=drag_fields) sat = Satrec.twoline2rv(line1, line2) elif getattr(orbit, "_sat", None) is not None: # Path A: the orbit was built from a TLE and still carries its @@ -576,6 +580,51 @@ def __eq__(self, rhs): ) +class _StateInterpolant: + """Piecewise quintic Hermite interpolant through cached RK states. + + Positions are interpolated from position, velocity and acceleration at + each step boundary, and velocities are the derivative of that quintic, so + the interpolation error is O(h^6) in position. The cubic spline through + positions and velocities that this replaces was O(h^4): an RK8 arc + accurate to 1e-8 m at its steps was only 2.5 cm accurate between them at + h = 40 s in LEO. Accelerations are evaluated once per step boundary and + reused when the cache grows. + """ + + def __init__(self, interpolant, accels): + self._interpolant = interpolant + self.accels = accels + + @classmethod + def build(cls, propagator, times, states, propkw, previous=None, fallback=None): + accels = dict(previous.accels) if isinstance(previous, cls) else {} + accel = getattr(propagator, "accel", None) + if accel is not None and len(times) >= 2: + try: + from scipy.interpolate import BPoly + derivs = np.empty((len(times), 3, 3), dtype=np.float64) + for k, (t, state) in enumerate(zip(times, states)): + key = float(t) + if key not in accels: + accels[key] = np.asarray( + _safe_accel(accel, state[0:3], state[3:6], t, propkw), dtype=np.float64 + ).reshape(3) + derivs[k, 0] = state[0:3] + derivs[k, 1] = state[3:6] + derivs[k, 2] = accels[key] + return cls(BPoly.from_derivatives(times, derivs), accels) + except Exception: # fall back to the cubic spline below + pass + k = min(3, len(times) - 1) + return fallback(times, states, k=k) + + def __call__(self, t): + position = self._interpolant(t) + velocity = self._interpolant.derivative()(t) + return np.concatenate([position, velocity], axis=-1) + + class RKPropagator(Propagator, ABC): """Abstract base class for Runge-Kutta-based orbit propagators. @@ -649,8 +698,10 @@ def _getRVOneWithMask(self, orbit, tQuery): return out[:, 0:3], out[:, 3:6], valid if remake_spline or spline is None: - k = min(3, len(times_arr) - 1) # cubic, stable - spline = make_interp_spline(times_arr, states_arr, k=k) + spline = _StateInterpolant.build( + self, times_arr, states_arr, orbit.propkw, + previous=spline, fallback=make_interp_spline, + ) container.clear() container.extend([times, states, h_pre, h_app, spline]) times_arr = np.asarray(times, dtype=np.float64) @@ -1174,14 +1225,19 @@ def _prop(self, times, states, h, tthresh, propkw): h1 = w1 * h h0 = w0 * h + # Substep epochs are offsets from t, and t itself advances by + # exactly one h per step. Chaining t + h1 + h0 + h1 rounds three + # times at the ulp of t (2.4e-7 s at GPS 1.4e9), and the bias + # accumulates: the stored epochs drifted from the integrated + # states and capped the method near 2 m over 6000 s. r, v = self._leapfrog_step(self.accel, r, v, t, h1, propkw) t1 = t + h1 r, v = self._leapfrog_step(self.accel, r, v, t1, h0, propkw) - t2 = t1 + h0 + t2 = t + (h1 + h0) r, v = self._leapfrog_step(self.accel, r, v, t2, h1, propkw) - t = t2 + h1 + t = t + h state = np.hstack([r, v]) diff --git a/ssapy/rvsampler.py b/ssapy/rvsampler.py index 0a2e8f06..10b0685d 100644 --- a/ssapy/rvsampler.py +++ b/ssapy/rvsampler.py @@ -530,16 +530,18 @@ def circular_guess(arc, indices=[0, 1]): Returns ------- - state: array_like, (6,) + state : array_like, (6,) state corresponding to orbit passing through points + epoch : astropy.time.Time + epoch of `state`: the midpoint of the two observations, or the + first selected observation when the arc carries rates ('pmra'). """ from astropy import units as u usepm = 'pmra' in arc.dtype.names - if not usepm: - observations = arc[indices] - else: - observations = arc + # With rates only the first selected row is needed (so a one-row arc + # works); it used to be row 0 of the whole arc whatever `indices` said. + observations = arc[[indices[0]]] if usepm else arc[indices] rStation = observations['rStation_GCRF'].to(u.m).value ra = observations['ra'] dec = observations['dec'] diff --git a/ssapy/utils.py b/ssapy/utils.py index 1a45d1b8..0cb6d1a4 100644 --- a/ssapy/utils.py +++ b/ssapy/utils.py @@ -5,6 +5,7 @@ import os import re import numpy as np +import warnings from astropy.time import Time as _Time import astropy.units as u from typing import Union, Tuple @@ -30,23 +31,62 @@ def _emcee_version_before_3(version_string): import astropy._erfa as erfa +_LFS_POINTER_PREFIX = b"version https://git-lfs.github.com/spec/v1" + + +def _is_lfs_pointer(path): + """True if ``path`` is an un-fetched git LFS pointer rather than data.""" + try: + if os.path.getsize(path) > 1024: + return False + with open(path, "rb") as f: + return f.read(len(_LFS_POINTER_PREFIX)) == _LFS_POINTER_PREFIX + except OSError: + return False + + +def _ssapy_data_files(): + """Paths of the files shipped by the optional ``llnl-ssapy-data`` package.""" + try: + import ssapy_data + except ImportError: + return [] + try: + return [os.fspath(path) for path in ssapy_data.iter_data_files()] + except Exception: + return [] + + def find_file(filename, ext=None): - """ Find a file in the current directory or the ssapy datadir. If ext is - not None, also try appending ext to the filename. - """ - candidates = [ - filename, - os.path.join(datadir, filename), - ] - if ext is not None: - candidates.extend([ - filename + ext, - os.path.join(datadir, filename + ext), - ]) + """ Find a file in the current directory or the ssapy datadir (the + ``ssapy/`` tree of the ``llnl-ssapy-data`` package). If ext is not None, + also try appending ext to the filename. + + Git LFS pointer files (left over from SSAPy versions that stored data + with Git LFS) are skipped, and the rest of ``llnl-ssapy-data`` is searched + by file name. + """ + names = [filename] if ext is None else [filename, filename + ext] + candidates = [] + for name in names: + candidates.extend([name, os.path.join(datadir, name)]) for candidate in candidates: - if os.path.isfile(candidate): + if os.path.isfile(candidate) and not _is_lfs_pointer(candidate): return candidate - raise FileNotFoundError(filename) + basenames = {os.path.basename(name) for name in names} + for path in _ssapy_data_files(): + if os.path.basename(path) in basenames and os.path.isfile(path) and not _is_lfs_pointer(path): + return path + for candidate in candidates: + if os.path.isfile(candidate): + raise FileNotFoundError( + f"{candidate} is a git LFS pointer, not the data file. SSAPy's data now " + "comes from the llnl-ssapy-data package: pip install 'llnl-ssapy-data>=0.2.0'." + ) + raise FileNotFoundError( + f"{filename} was not found in the working directory or in llnl-ssapy-data " + f"({datadir}); install it with pip install 'llnl-ssapy-data>=0.2.0'." + ) def _wrapToPi(angle): @@ -987,15 +1027,13 @@ def lb_to_tan(lb, b, mul=None, mub=None, lcen=None, bcen=None): dechat2 = np.cross(unit, rahat2) dechat2 /= np.sqrt(np.sum(dechat2**2, axis=1, keepdims=True)) vv = mul[:, None] * rahat2 + mub[:, None] * dechat2 + # x = rahat . unit and y = dechat . unit with fixed rahat/dechat, so + # their time derivatives are exactly rahat . (d unit/dt) and + # dechat . (d unit/dt). The radial 1/cos(rho) stretch that used to + # follow is the gnomonic derivative; applied to this orthographic + # projection it overstated the rates by 1/cos(rho). vx = np.sum(rahat * vv, axis=1) vy = np.sum(dechat * vv, axis=1) - rr = np.hypot(xx, yy) - m = np.abs(rr) > 1e-9 - vr = (vx[m] * xx[m] + vy[m] * yy[m]) / rr[m] - va = (vy[m] * xx[m] - vx[m] * yy[m]) / rr[m] - vr /= np.sum(unitcen[m, :] * unit[m, :], axis=1) - vx[m] = vr * xx[m] / rr[m] - va * yy[m] / rr[m] - vy[m] = vr * yy[m] / rr[m] + va * xx[m] / rr[m] res = res + (vx, vy) return res @@ -1196,7 +1234,12 @@ def sunPos(t, fast=True): # MG section 3.3.2 T = (_gpsToTT(t) - 51544.5) / 36525.0 M = 6.239998880168239 + 628.3019326367721 * T - lam = (4.938234585592756 + M + 0.03341335890206922 * np.sin(M) + 0.00034906585039886593 * np.sin(2 * M)) + # Omega + omega = 282.9400 deg at J2000, advancing 0.32327364 deg per + # Julian century in the J2000 frame (JPL approximate elements of the + # Earth-Moon barycentre). MG treat it as constant, which lets the + # longitude drift 11.6 arcsec/yr: 5.3 arcmin by 2026. + lam = (4.938234585592756 + 0.005642189402906841 * T + M + + 0.03341335890206922 * np.sin(M) + 0.00034906585039886593 * np.sin(2 * M)) rs = (149.619 - 2.499 * np.cos(M) - 0.021 * np.cos(2 * M)) * 1e9 obliquity = 0.40909280420293637 co, so = np.cos(obliquity), np.sin(obliquity) @@ -1863,20 +1906,15 @@ def rightascension_to_hourangle(right_ascension, local_time): Example: -------- - rightascension_to_hourangle("10:30:00", "12:45:00") -> "02:15:00" - rightascension_to_hourangle(157.5, 191.25) -> "02:15:00" + rightascension_to_hourangle("10:30:00", "12:45:00") -> "2:15:0" + rightascension_to_hourangle(157.5, 191.25) -> "2:15:0" """ - if type(right_ascension) is not str: - right_ascension = dd_to_hms(right_ascension) - if type(local_time) is not str: - local_time = dd_to_dms(local_time) - _ra = float(right_ascension.split(':')[0]) - _lt = float(local_time.split(':')[0]) - if _ra > _lt: - __ltm, __lts = local_time.split(':')[1:] - local_time = f'{24 + _lt}:{__ltm}:{__lts}' - - return dd_to_dms(hms_to_dd(local_time) - hms_to_dd(right_ascension)) + ra_deg = hms_to_dd(right_ascension) if isinstance(right_ascension, str) else float(right_ascension) + lst_deg = hms_to_dd(local_time) if isinstance(local_time, str) else float(local_time) + # Hour angle = local sidereal time - right ascension, wrapped to [0, 24) h + # and returned as hours. The previous code formatted the difference in + # degrees as D:M:S, so a 2 h hour angle came back as "30:0:0". + return dd_to_hms((lst_deg - ra_deg) % 360.0) def equatorial_to_horizontal(observer_latitude, declination, right_ascension=None, hour_angle=None, local_time=None, hms=False): @@ -1902,17 +1940,18 @@ def equatorial_to_horizontal(observer_latitude, declination, right_ascension=Non tuple Azimuth and altitude in degrees. """ - if right_ascension is not None and hour_angle is not None: - print('Both right_ascension and hour_angle parameters are provided.\nUsing hour_angle for calculations.') - if hms: + if hour_angle is not None: + if right_ascension is not None: + warnings.warn( + "Both right_ascension and hour_angle parameters are provided; using hour_angle.", + UserWarning, + stacklevel=2, + ) + if hms or isinstance(hour_angle, str): hour_angle = hms_to_dd(hour_angle) elif right_ascension is not None: - hour_angle = rightascension_to_hourangle(right_ascension, local_time) - if hms: - hour_angle = hms_to_dd(hour_angle) - elif hour_angle is not None: - if hms: - hour_angle = hms_to_dd(hour_angle) + # rightascension_to_hourangle always returns an HH:MM:SS string. + hour_angle = hms_to_dd(rightascension_to_hourangle(right_ascension, local_time)) else: raise ValueError('Either right_ascension or hour_angle must be provided.') @@ -1922,12 +1961,14 @@ def equatorial_to_horizontal(observer_latitude, declination, right_ascension=Non altitude = zenithangle_to_altitude(zenith_angle, deg=False) - _num = np.sin(declination) - np.sin(observer_latitude) * np.cos(zenith_angle) - _den = np.cos(observer_latitude) * np.sin(zenith_angle) - azimuth = np.arccos(_num / _den) - - if observer_latitude < 0: - azimuth = np.pi - azimuth + # Azimuth from north through east. arccos cannot tell east from west and + # put every object west of the meridian in the east; the sign of + # sin(hour angle) does (positive hour angle = west). + azimuth = np.mod(np.arctan2( + -np.cos(declination) * np.sin(hour_angle), + np.sin(declination) * np.cos(observer_latitude) + - np.cos(declination) * np.sin(observer_latitude) * np.cos(hour_angle), + ), 2 * np.pi) altitude, azimuth = np.degrees([altitude, azimuth]) return azimuth, altitude @@ -1954,18 +1995,20 @@ def horizontal_to_equatorial(observer_latitude, azimuth, altitude): altitude, azimuth, latitude = np.radians([altitude, azimuth, observer_latitude]) zenith_angle = altitude_to_zenithangle(altitude, deg=False) - zenith_angle = [-zenith_angle if latitude < 0 else zenith_angle][0] - - declination = np.sin(latitude) * np.cos(zenith_angle) - declination = declination + (np.cos(latitude) * np.sin(zenith_angle) * np.cos(azimuth)) - declination = np.arcsin(declination) - - _num = np.cos(zenith_angle) - np.sin(latitude) * np.sin(declination) - _den = np.cos(latitude) * np.cos(declination) - hour_angle = np.arccos(_num / _den) + declination = np.arcsin( + np.sin(latitude) * np.cos(zenith_angle) + + np.cos(latitude) * np.sin(zenith_angle) * np.cos(azimuth) + ) - if (latitude > 0 > declination) or (latitude < 0 < declination): - hour_angle = 2 * np.pi - hour_angle + # atan2 resolves the quadrant (east of the meridian, sin azimuth > 0, is a + # negative hour angle, i.e. above 12 h) and stays well conditioned on the + # meridian, where arccos of a value near 1 loses half its digits. The + # previous test on the signs of latitude and declination mirrored the hour + # angle in the wrong cases. + hour_angle = np.mod(np.arctan2( + -np.sin(azimuth) * np.sin(zenith_angle), + np.cos(zenith_angle) * np.cos(latitude) - np.sin(zenith_angle) * np.sin(latitude) * np.cos(azimuth), + ), 2 * np.pi) declination, hour_angle = np.degrees([declination, hour_angle]) @@ -2123,7 +2166,9 @@ def equatorial_to_ecliptic(right_ascension, declination, degrees=False): else: ra, dec = right_ascension, declination ec_latitude = np.arcsin(cos_ec * np.sin(dec) - sin_ec * np.cos(dec) * np.sin(ra)) - ec_longitude = np.arctan((cos_ec * np.cos(dec) * np.sin(ra) + sin_ec * np.sin(dec)) / (np.cos(dec) * np.cos(ra))) + # arctan2, not arctan of the ratio: arctan folds right ascensions between + # 90 and 270 deg onto the opposite side of the sky. + ec_longitude = np.arctan2(cos_ec * np.cos(dec) * np.sin(ra) + sin_ec * np.sin(dec), np.cos(dec) * np.cos(ra)) if degrees: return deg0to360(np.degrees(ec_longitude)), np.degrees(ec_latitude) else: @@ -2150,7 +2195,9 @@ def ecliptic_to_equatorial(lon, lat, degrees=False): lon, lat = np.radians(lon), np.radians(lat) else: lon, lat = lon, lat - ra = np.arctan((cos_ec * np.cos(lat) * np.sin(lon) - sin_ec * np.sin(lat)) / (np.cos(lat) * np.cos(lon))) + # arctan2 resolves the quadrant; arctan of the ratio returned right + # ascensions between 90 and 270 deg on the opposite side of the sky. + ra = rad0to2pi(np.arctan2(cos_ec * np.cos(lat) * np.sin(lon) - sin_ec * np.sin(lat), np.cos(lat) * np.cos(lon))) dec = np.arcsin(cos_ec * np.sin(lat) + sin_ec * np.cos(lat) * np.sin(lon)) if degrees: return np.degrees(ra), np.degrees(dec) @@ -2267,10 +2314,12 @@ def dms_to_dd(dms): # Degree minute second to Degree decimal """ dms, out = [[dms] if type(dms) is str else dms][0], [] for i in dms: - deg, minute, sec = [float(j) for j in i.split(':')] - if deg < 0: - minute, sec = float(f'-{minute}'), float(f'-{sec}') - out.append(deg + minute / 60 + sec / 3600) + text = i.strip() + # The sign belongs to the whole angle, so "-00:30:00" is -0.5 deg; a + # sign carried only by the degree field is lost when degrees are 0. + sign = -1.0 if text.startswith('-') else 1.0 + deg, minute, sec = [abs(float(j)) for j in text.lstrip('+-').split(':')] + out.append(sign * (deg + minute / 60 + sec / 3600)) return [out[0] if type(dms) is str or len(dms) == 1 else out][0] @@ -2393,9 +2442,8 @@ def dd_to_hms(degree_decimal): Notes: ------ - - Decimal degrees are divided by 15 to convert to hours. - - If the input DD value is negative, the function assumes the absolute value for conversion - and prints a warning message. + - Decimal degrees are wrapped into [0, 360) and divided by 15 to convert to hours, so + the result always lies in [0, 24) h. - Handles edge cases where seconds reach 60, incrementing minutes accordingly. - Returns seconds as an integer if the value is a whole number. @@ -2405,18 +2453,16 @@ def dd_to_hms(degree_decimal): '12:34:56' >>> dd_to_hms(-236.375) - '15:45:30' # Assumes positive value for conversion. + '8:14:30' # -236.375 deg is 123.625 deg. >>> dd_to_hms("12:34:56") # DMS string converted to DD first. '0:50:18.4' """ if type(degree_decimal) is str: degree_decimal = dms_to_dd(degree_decimal) - if degree_decimal < 0: - print('dd for HMS conversion cannot be negative, assuming positive.') - _dd = -degree_decimal / 15 - else: - _dd = degree_decimal / 15 + # Angles are wrapped into [0, 360) deg, so -15 deg is 23h; taking the + # absolute value would turn it into 1h. + _dd = float(np.mod(degree_decimal, 360.0)) / 15 _h = int(np.trunc(_dd)) __h = _dd - np.trunc(_dd) _m, __m = np.trunc(__h * 60), __h * 60 - np.trunc(__h * 60) @@ -2426,6 +2472,8 @@ def dd_to_hms(degree_decimal): _m, _s = _m + 1, '00' if _m == 60: _h, _m = _h + 1, 0 + if _h >= 24: + _h = 0 return f'{int(_h)}:{int(_m)}:{_s}' diff --git a/tests/test_accel.py b/tests/test_accel.py index 20de4ecb..33b19e22 100644 --- a/tests/test_accel.py +++ b/tests/test_accel.py @@ -4,7 +4,6 @@ import sys import subprocess from types import SimpleNamespace -from pathlib import Path import pytest import erfa @@ -18,20 +17,29 @@ from ssapy.constants import EARTH_MU, EARTH_RADIUS -def _is_lfs_pointer(path): + +def _has_data(name): + # SSAPy's data comes from llnl-ssapy-data; skip only if it is missing. + try: + ssapy.utils.find_file(name) + return True + except FileNotFoundError: + return False + + +def _has_ephemeris(): try: - with open(path, "rb") as f: - first = f.readline() - return first.startswith(b"version https://git-lfs.github.com/spec/v1") - except Exception: + ssapy.body._planetary_ephemeris_path() + return True + except FileNotFoundError: return False -HAS_EGM84 = not _is_lfs_pointer(Path(ssapy.datadir) / "egm84.egm.cof") -HAS_DE430 = not _is_lfs_pointer(Path(ssapy.datadir) / "de430.bsp") -HAS_MOON_PA = not _is_lfs_pointer(Path(ssapy.datadir) / "moon_pa_de440_200625.bpc") +HAS_EGM84 = _has_data("egm84.egm.cof") +HAS_EPHEMERIS = _has_ephemeris() +HAS_MOON_PA = _has_data("moon_pa_de440_200625.bpc") -HAS_BODY_DATA = HAS_EGM84 and HAS_DE430 and HAS_MOON_PA +HAS_BODY_DATA = HAS_EGM84 and HAS_EPHEMERIS and HAS_MOON_PA def test_harmonic_coefficients_error_and_identity_paths(tmp_path): @@ -86,7 +94,7 @@ def test_harmonic_coefficients_error_and_identity_paths(tmp_path): assert hc1 != object() -def test_accel_harmonic_clamps_requested_degree_and_order(capsys): +def test_accel_harmonic_clamps_requested_degree_and_order(): harmonics = SimpleNamespace( name="tiny", radius=1.0, @@ -100,12 +108,13 @@ def test_accel_harmonic_clamps_requested_degree_and_order(capsys): assert default_accel.n_max == 1 assert default_accel.m_max == 1 - clamped = ssapy.AccelHarmonic(body, n_max=5, m_max=6) + with pytest.warns(UserWarning) as record: + clamped = ssapy.AccelHarmonic(body, n_max=5, m_max=6) assert clamped.n_max == 1 assert clamped.m_max == 1 - out = capsys.readouterr().out - assert "provided degree" in out - assert "provided order" in out + messages = " ".join(str(w.message) for w in record) + assert "provided degree" in messages + assert "provided order" in messages def _make_body_dependent_objects(): @@ -142,7 +151,7 @@ def _make_body_dependent_objects(): earth = moon = sun = earth_MG = moon_MG = sun_MG = None -@pytest.mark.skipif(not HAS_EGM84, reason="EGM84 gravity file unavailable (Git LFS pointer)") +@pytest.mark.skipif(not HAS_EGM84, reason="EGM84 gravity file unavailable (llnl-ssapy-data not installed)") @timer def test_MG_3_1(): """Exercise 3.1 from Montenbruck and Gill diff --git a/tests/test_angle_conversion_references.py b/tests/test_angle_conversion_references.py new file mode 100644 index 00000000..e2dd0756 --- /dev/null +++ b/tests/test_angle_conversion_references.py @@ -0,0 +1,72 @@ +"""Sexagesimal, hour-angle, horizon, and ecliptic conversions against references.""" + +import numpy as np +import pytest +import astropy.units as u +from astropy.coordinates import Angle + +from ssapy import utils + + +def test_negative_sub_degree_dms_keeps_its_sign(): + # The sign belongs to the whole angle: -00:30:00 is -0.5 deg (exact). + for text, degrees in {"-00:30:00": -0.5, "-0:00:36": -0.01, "-10:30:36": -10.51, "+00:30:00": 0.5}.items(): + assert utils.dms_to_dd(text) == pytest.approx(degrees, abs=1e-12), text + + +@pytest.mark.parametrize("degrees", [-15.0, -236.375, 0.0, 15.0, 187.5, 359.99999999, 757.5]) +def test_dd_to_hms_wraps_into_one_day_like_astropy(degrees): + # astropy Angle.wrap_at(360 deg).hour gives the same instant within + # 1e-4 s of time, and hours lie in [0, 24). + hours, minutes, seconds = (float(part) for part in utils.dd_to_hms(degrees).split(":")) + assert 0 <= hours < 24 + expected = Angle(degrees * u.deg).wrap_at(360 * u.deg).hour * 3600.0 + actual = hours * 3600 + minutes * 60 + seconds + assert min(abs(actual - expected), 86400 - abs(actual - expected)) < 1e-4 + + +def test_hour_angle_is_sidereal_time_minus_right_ascension_in_hours(): + # HA = LST - RA, wrapped to [0, 24) h (exact for these inputs). + assert utils.rightascension_to_hourangle("10:30:00", "12:45:00") == "2:15:0" + assert utils.rightascension_to_hourangle(157.5, 191.25) == "2:15:0" + assert utils.rightascension_to_hourangle("23:00:00", "01:00:00") == "2:0:0" + + +@pytest.mark.parametrize("latitude", [37.68, -33.87]) +@pytest.mark.parametrize("hour_angle", [-60.0, 0.0, 45.0, 120.0]) +@pytest.mark.parametrize("declination", [-20.0, 20.0]) +def test_horizon_conversions_match_astropy_altaz_geometry(latitude, hour_angle, declination): + # The topocentric direction of (HA, dec) in the local east/north/up frame, + # built from rotation matrices, gives the reference azimuth (from north + # through east) and altitude; the forward and inverse transforms agree + # with it to 1e-9 deg in both hemispheres. + phi, h, d = np.radians([latitude, hour_angle, declination]) + east = -np.cos(d) * np.sin(h) + north = np.sin(d) * np.cos(phi) - np.cos(d) * np.sin(phi) * np.cos(h) + up = np.sin(d) * np.sin(phi) + np.cos(d) * np.cos(phi) * np.cos(h) + az_ref = np.degrees(np.arctan2(east, north)) % 360.0 + alt_ref = np.degrees(np.arcsin(up)) + az, alt = utils.equatorial_to_horizontal(latitude, declination, hour_angle=hour_angle) + assert alt == pytest.approx(alt_ref, abs=1e-9) + if abs(alt_ref) < 89.9: + assert ((az - az_ref + 180.0) % 360.0) - 180.0 == pytest.approx(0.0, abs=1e-9) + ha_back, dec_back = utils.horizontal_to_equatorial(latitude, az, alt) + assert dec_back == pytest.approx(declination, abs=1e-9) + assert ((ha_back - hour_angle + 180.0) % 360.0) - 180.0 == pytest.approx(0.0, abs=1e-9) + + +@pytest.mark.parametrize("ra, dec", [(30.0, 10.0), (120.0, 10.0), (200.0, -20.0), (300.0, 40.0)]) +def test_ecliptic_conversions_match_the_obliquity_rotation(ra, dec): + # Rotating the unit vector about +x by the J2000 obliquity (SSAPy's + # cos_ec/sin_ec) gives the reference ecliptic longitude and latitude in + # every quadrant; both directions agree with it to 1e-9 deg. + ce, se = utils.cos_ec, utils.sin_ec + ra_r, dec_r = np.radians([ra, dec]) + x = np.cos(dec_r) * np.cos(ra_r) + y = ce * np.cos(dec_r) * np.sin(ra_r) + se * np.sin(dec_r) + z = -se * np.cos(dec_r) * np.sin(ra_r) + ce * np.sin(dec_r) + lon_ref, lat_ref = np.degrees(np.arctan2(y, x)) % 360.0, np.degrees(np.arcsin(z)) + lon, lat = utils.equatorial_to_ecliptic(ra, dec, degrees=True) + assert lon == pytest.approx(lon_ref, abs=1e-9) and lat == pytest.approx(lat_ref, abs=1e-9) + ra_back, dec_back = utils.ecliptic_to_equatorial(lon_ref, lat_ref, degrees=True) + assert ra_back == pytest.approx(ra, abs=1e-9) and dec_back == pytest.approx(dec, abs=1e-9) diff --git a/tests/test_body.py b/tests/test_body.py index 5069e5f5..96bed091 100644 --- a/tests/test_body.py +++ b/tests/test_body.py @@ -1,3 +1,4 @@ +import os import numpy as np import pytest @@ -163,7 +164,7 @@ def fake_open(path): assert planet_position.kernel == "planet-kernel" assert planet_position.planet_index == 5 - assert opened[0].endswith("de430.bsp") + assert os.path.basename(opened[0]) == "de440s.bsp" def test_orientation_helpers(monkeypatch): diff --git a/tests/test_body_constants.py b/tests/test_body_constants.py new file mode 100644 index 00000000..bfc6319e --- /dev/null +++ b/tests/test_body_constants.py @@ -0,0 +1,15 @@ +import pytest + +from ssapy import constants + +G = 6.67430e-11 # CODATA 2018 [m^3 kg^-1 s^-2] + + +@pytest.mark.parametrize("body", ["SUN", "MOON", "MERCURY", "VENUS", "EARTH", "MARS", "JUPITER", "SATURN", "URANUS", "NEPTUNE"]) +def test_masses_are_consistent_with_gravitational_parameters(body): + # Each tabulated mass must equal its GM / G (CODATA 2018 G) to 0.2%; the + # masses and GMs come from different compilations, which differ by <0.05% + # here. VENUS_MASS was 4.687e24 kg against GM / G = 4.8673e24 kg (3.7%). + mu = getattr(constants, f"{body}_MU") + mass = getattr(constants, f"{body}_MASS") + assert mass == pytest.approx(mu / G, rel=2e-3) diff --git a/tests/test_data_file_lookup.py b/tests/test_data_file_lookup.py new file mode 100644 index 00000000..1472528d --- /dev/null +++ b/tests/test_data_file_lookup.py @@ -0,0 +1,58 @@ +import os + +import pytest + +from ssapy import utils + +POINTER = ( + "version https://git-lfs.github.com/spec/v1\n" + "oid sha256:0123456789abcdef0123456789abcdef0123456789abcdef0123456789abcdef\n" + "size 119741440\n" +) + + +def test_lfs_pointer_is_not_returned_as_data(tmp_path, monkeypatch): + # A clone without `git lfs pull` holds 130-byte pointer files in place of + # the kernels; find_file used to return them, and jplephem then failed on + # an unreadable "kernel". It must now raise a FileNotFoundError that says + # what the file is (git LFS spec v1 pointer format). + (tmp_path / "de440s.bsp").write_text(POINTER) + monkeypatch.setattr(utils, "datadir", str(tmp_path)) + monkeypatch.setattr(utils, "_ssapy_data_files", lambda: []) + with pytest.raises(FileNotFoundError, match="git LFS pointer.*llnl-ssapy-data"): + utils.find_file("de440s.bsp") + + +def test_ssapy_data_is_searched_when_the_package_data_is_missing(tmp_path, monkeypatch): + # With the package copy only a pointer, a real file of the same name + # shipped by llnl-ssapy-data is returned instead. + local = tmp_path / "local" + shipped = tmp_path / "shipped" / "ephemerides" + local.mkdir() + shipped.mkdir(parents=True) + (local / "de440s.bsp").write_text(POINTER) + real = shipped / "de440s.bsp" + real.write_bytes(b"DAF/SPK " + b"\0" * 2048) + monkeypatch.setattr(utils, "datadir", str(local)) + monkeypatch.setattr(utils, "_ssapy_data_files", lambda: [str(real)]) + assert os.path.samefile(utils.find_file("de440s.bsp"), real) + + +def test_datadir_is_the_ssapy_tree_of_llnl_ssapy_data(): + # SSAPy no longer ships data: datadir is ssapy/ inside llnl-ssapy-data, and + # every file the body and gravity loaders need resolves there. + import ssapy + import ssapy_data + from ssapy.body import _planetary_ephemeris_path + + assert os.path.samefile(ssapy.datadir, os.path.join(os.path.dirname(ssapy_data.__file__), "data", "ssapy")) + for name in ("moon_pa_de440_200625.bpc", "egm84.egm", "egm2008.egm.cof", "gggrx_1200a_sha.tab", "earth.png", "moon.png"): + assert os.path.dirname(utils.find_file(name)) == os.path.normpath(ssapy.datadir) + assert os.path.basename(_planetary_ephemeris_path()) == "de440s.bsp" + + +def test_missing_file_names_the_data_package(monkeypatch, tmp_path): + monkeypatch.setattr(utils, "datadir", str(tmp_path)) + monkeypatch.setattr(utils, "_ssapy_data_files", lambda: []) + with pytest.raises(FileNotFoundError, match="llnl-ssapy-data"): + utils.find_file("no_such_file.bsp") diff --git a/tests/test_drag_sun_frame.py b/tests/test_drag_sun_frame.py new file mode 100644 index 00000000..e55e73dc --- /dev/null +++ b/tests/test_drag_sun_frame.py @@ -0,0 +1,29 @@ +from types import SimpleNamespace + +import numpy as np +from astropy.coordinates import TETE, get_sun +from astropy.time import Time + +from ssapy.accel import AccelDrag + + +def test_drag_density_sees_the_sun_of_date(): + # Harris-Priester's diurnal bulge follows the Sun's right ascension and + # declination in the true-of-date frame the density is evaluated in. + # The values AccelDrag passes must match astropy's true-equator, + # true-equinox (TETE) Sun to 0.05 deg (fast sunPos is good to 0.6 arcmin); + # the GCRF values it passed before were 0.37 deg off in RA in 2026. + t = Time("2026-10-08T00:00:00", scale="utc") + drag = AccelDrag() + seen = {} + + def capture(x, y, z, ra_sun, dec_sun): + seen["ra"], seen["dec"] = ra_sun, dec_sun + return 1e-12 + + drag.atm = SimpleNamespace(density=capture) # the C++ model's methods are read-only + drag(np.array([6778e3, 0.0, 0.0]), np.array([0.0, 7670.0, 0.0]), t.gps, area=1.0, mass=100.0, CD=2.2) + sun = get_sun(t).transform_to(TETE(obstime=t)) + d_ra = (np.degrees(seen["ra"]) - sun.ra.deg + 180.0) % 360.0 - 180.0 + assert abs(d_ra) < 0.05 + assert abs(np.degrees(seen["dec"]) - sun.dec.deg) < 0.05 diff --git a/tests/test_earth_orientation_cache.py b/tests/test_earth_orientation_cache.py new file mode 100644 index 00000000..7f14d008 --- /dev/null +++ b/tests/test_earth_orientation_cache.py @@ -0,0 +1,32 @@ +import numpy as np +import pytest +import astropy.units as u +from astropy.coordinates import GCRS, ITRS, CartesianRepresentation +from astropy.time import Time + +from ssapy.body import EarthOrientation + + +def _angle_arcsec(a, b): + # Small-angle rotation between two frames; the Frobenius norm of + # (A B^T - I) / sqrt(2) avoids arccos's loss of precision near 1. + return np.degrees(np.linalg.norm(a @ b.T - np.eye(3)) / np.sqrt(2.0)) * 3600.0 + + +@pytest.mark.parametrize("days", [1.0, 7.0, 29.0]) +def test_cached_earth_orientation_matches_a_fresh_evaluation(days): + # A long-lived EarthOrientation (as held by Body / AccelHarmonic over a + # propagation) must agree with a freshly built one to 0.005 arcsec at any + # later time, and with astropy's GCRS -> ITRS rotation to 0.5 arcsec (the + # IAU 1980 / 2000A model difference; measured 0.27 arcsec). The old 30-day + # cache drifted 0.06, 1.0 and 5.1 arcsec after 1, 7 and 29 days. + t0 = Time("2026-10-08T00:00:00", scale="utc").gps + held = EarthOrientation() + held(t0) + t = t0 + days * 86400.0 + assert _angle_arcsec(held(t), EarthOrientation()(t)) < 0.005 + r = np.array([7000e3, 1000e3, 500e3]) + when = Time(t, format="gps") + expected = GCRS(CartesianRepresentation(r * u.m), obstime=when).transform_to(ITRS(obstime=when)).cartesian.xyz.to_value(u.m) + got = held(t) @ r + assert np.degrees(np.linalg.norm(np.cross(got, expected)) / (np.linalg.norm(got) * np.linalg.norm(expected))) * 3600.0 < 0.5 diff --git a/tests/test_ephemeris_selection.py b/tests/test_ephemeris_selection.py new file mode 100644 index 00000000..d4e1ce21 --- /dev/null +++ b/tests/test_ephemeris_selection.py @@ -0,0 +1,174 @@ +"""Planetary ephemeris selection, long-span downloads, and their failure modes. + +The download tests stand in a local file for NAIF, so they run offline: a +2000-2050 excerpt of the shipped de440s.bsp plays the shipped kernel and the +full de440s.bsp plays the long-span kernel "downloaded" for other epochs. +""" +import io +import os +import subprocess +import sys +import warnings + +import numpy as np +import pytest +from astropy.time import Time +from jplephem.spk import SPK + +import ssapy +from ssapy import ephemeris +from ssapy.body import MoonPosition, SunPosition +from ssapy.utils import find_file + + +@pytest.fixture(autouse=True) +def _reset_selection(monkeypatch, tmp_path): + monkeypatch.setattr(ephemeris, "_selected", None) + monkeypatch.delenv("SSAPY_EPHEMERIS", raising=False) + monkeypatch.delenv("SSAPY_EPHEMERIS_DOWNLOAD", raising=False) + monkeypatch.setenv("SSAPY_DATA_CACHE", str(tmp_path / "cache")) + monkeypatch.setattr(ephemeris, "_open_kernels", {}) + + +def _gps(iso): + # TT avoids ERFA's "dubious year" warnings for UTC far from the present. + return Time(iso, scale="tt").gps + + +def test_de440_and_de430_agree_to_their_published_difference(): + # DE440 is the default; DE430 is selectable. In 2026 the two solutions' + # geocentric Moon differs by <= 10 m and Sun by <= 0.5 km over 1975-2050 + # (measured from the NAIF kernels); hold them to 20 m and 1 km. + t = _gps("2026-10-08T00:00:00") + assert ephemeris.get_planetary_ephemeris() == "de440" + moon440, sun440 = MoonPosition()(t), SunPosition()(t) + ephemeris.set_planetary_ephemeris("de430") + moon430, sun430 = MoonPosition()(t), SunPosition()(t) + assert 0 < np.linalg.norm(moon440 - moon430) < 20.0 + assert np.linalg.norm(sun440 - sun430) < 1e3 + with pytest.raises(ValueError): + ephemeris.set_planetary_ephemeris("de405") + + +def test_environment_variable_selects_the_ephemeris(monkeypatch): + monkeypatch.setenv("SSAPY_EPHEMERIS", "DE430") + assert ephemeris.get_planetary_ephemeris() == "de430" + assert os.path.basename(ssapy.body._planetary_ephemeris_path()) == "de430_1900_2150.bsp" + + +def test_no_download_inside_the_shipped_span(monkeypatch): + def no_network(url): + raise AssertionError("no download expected") + + monkeypatch.setattr(ephemeris, "_open_url", no_network) + with warnings.catch_warnings(): + warnings.simplefilter("error") + MoonPosition()(_gps(["1900-01-01T00:00:00", "2026-10-08T00:00:00", "2149-12-31T00:00:00"])) + + +def _harness(monkeypatch, tmp_path, *, long_solution="DE440", served=None): + """Shipped kernel = 2000-2050 excerpt of de440s; long kernel = de440s.""" + shipped_full = find_file("de440s.bsp") + excerpt = tmp_path / "short_test.bsp" + subprocess.run([sys.executable, "-m", "jplephem", "excerpt", "2000/1/1", "2050/1/1", shipped_full, str(excerpt)], + check=True, stdout=subprocess.DEVNULL) + payload = open(shipped_full, "rb").read() if served is None else served + monkeypatch.chdir(tmp_path) + full = SPK.open(shipped_full) + start = max(seg.start_jd for seg in full.segments) + end = min(seg.end_jd for seg in full.segments) + family = ( + ephemeris.Kernel("short_test.bsp", 2451544.5, 2469807.5, solution="DE440"), + ephemeris.Kernel("long_test.bsp", start, end, ephemeris._sha256(shipped_full), + os.path.getsize(shipped_full), long_solution), + ) + monkeypatch.setitem(ephemeris.EPHEMERIDES, "de440", family) + calls = [] + + def fake_open_url(url): + calls.append(url) + return io.BytesIO(payload) + + monkeypatch.setattr(ephemeris, "_open_url", fake_open_url) + return full, calls + + +def _direct_moon(kernel, t): + mjd = ssapy.utils._gpsToTT(t) + return (kernel[3, 301].compute(2400000.5, mjd) - kernel[3, 399].compute(2400000.5, mjd)) * 1e3 + + +def test_epochs_outside_the_shipped_span_download_with_a_warning(monkeypatch, tmp_path): + # An epoch beyond the shipped span downloads the long kernel once, warns, + # caches it, and evaluates it exactly; epochs inside stay on the shipped + # kernel, and the two agree at the boundary (same solution). + full, calls = _harness(monkeypatch, tmp_path) + moon = MoonPosition() + times = _gps(["2026-10-08T00:00:00", "2100-06-01T00:00:00"]) + with pytest.warns(ephemeris.EphemerisDownloadWarning, match="long_test.bsp"): + got = moon(times) + assert len(calls) == 1 and calls[0].endswith("long_test.bsp") + assert os.path.isfile(tmp_path / "cache" / "long_test.bsp") + np.testing.assert_allclose(got[:, 1], _direct_moon(full, times[1]), rtol=0, atol=1e-6) + np.testing.assert_allclose(got[:, 0], _direct_moon(full, times[0]), rtol=0, atol=1e-3) + with warnings.catch_warnings(): + warnings.simplefilter("error") + moon(times) # cached: no second download or warning + assert len(calls) == 1 + + +def test_a_different_long_span_solution_warns(monkeypatch, tmp_path): + _harness(monkeypatch, tmp_path, long_solution="DE441") + with pytest.warns(ephemeris.EphemerisDownloadWarning): + with pytest.warns(ephemeris.EphemerisRangeWarning, match="DE441"): + MoonPosition()(_gps("2100-06-01T00:00:00")) + + +@pytest.mark.parametrize("iso, kernel", [ + ("1700-01-01T00:00:00", "de440.bsp"), + ("2500-01-01T00:00:00", "de440.bsp"), + ("1000-01-01T00:00:00", "de441_part-1.bsp"), + ("3000-01-01T00:00:00", "de441_part-2.bsp"), +]) +def test_long_span_epochs_route_to_full_de440_then_de441(monkeypatch, iso, kernel): + # With downloads disabled the error names the kernel that would be fetched: + # full DE440 for 1550-2650 outside de440s, DE441 beyond (NAIF coverage). + monkeypatch.setenv("SSAPY_EPHEMERIS_DOWNLOAD", "0") + with pytest.raises(ephemeris.EphemerisUnavailableError, match=kernel.replace(".", r"\.")): + MoonPosition()(_gps(iso)) + + +def test_offline_or_disabled_downloads_fail_clearly(monkeypatch, tmp_path): + _harness(monkeypatch, tmp_path) + t = _gps("2100-06-01T00:00:00") + + def offline(url): + raise OSError("network unreachable") + + monkeypatch.setattr(ephemeris, "_open_url", offline) + with pytest.warns(ephemeris.EphemerisDownloadWarning): + with pytest.raises(ephemeris.EphemerisUnavailableError, match="Could not download long_test.bsp"): + MoonPosition()(t) + assert not os.listdir(tmp_path / "cache") # no partial file left + + monkeypatch.setenv("SSAPY_EPHEMERIS_DOWNLOAD", "0") + with pytest.raises(ephemeris.EphemerisUnavailableError, match="downloads are disabled"): + MoonPosition()(t) + + +def test_a_corrupt_download_is_rejected(monkeypatch, tmp_path): + _harness(monkeypatch, tmp_path, served=b"not a kernel") + with pytest.warns(ephemeris.EphemerisDownloadWarning): + with pytest.raises(ephemeris.EphemerisUnavailableError, match="SHA-256"): + MoonPosition()(_gps("2100-06-01T00:00:00")) + assert not os.path.exists(tmp_path / "cache" / "long_test.bsp") + + +@pytest.mark.skipif(os.environ.get("SSAPY_TEST_NETWORK") != "1", reason="set SSAPY_TEST_NETWORK=1 to download from NAIF") +@pytest.mark.timeout(600) +def test_real_download_of_full_de440(): + # Downloads de440.bsp (119.8 MB) from NAIF, checks its SHA-256, and + # evaluates the Moon in 1700, outside de440s.bsp. + with pytest.warns(ephemeris.EphemerisDownloadWarning, match="de440.bsp"): + position = MoonPosition()(_gps("1700-01-01T00:00:00")) + assert 3.5e8 < np.linalg.norm(position) < 4.1e8 diff --git a/tests/test_harmonic_gravity_reference.py b/tests/test_harmonic_gravity_reference.py new file mode 100644 index 00000000..4444ca84 --- /dev/null +++ b/tests/test_harmonic_gravity_reference.py @@ -0,0 +1,101 @@ +"""AccelHarmonic against an independent spherical-harmonic potential. + +The coefficient files are parsed here without SSAPy's loaders, the potential +U = GM/r sum_n (R/r)^n sum_m Pbar_nm(sin phi) (C cos m lam + S sin m lam) is +built from scipy's associated Legendre functions with the 4-pi normalisation, +and its gradient is taken by central differences. The 20x20 accelerations +agree with AccelHarmonic to 1e-8 relative (finite-difference limited, +measured 1.3e-9 to 2.7e-9). Before AccelHarmonic used each model's own GM, +GRGM1200A lunar accelerations were 2.2e-7 low. +""" +from math import factorial + +import numpy as np +import pytest +from scipy.special import lpmv + +from ssapy.body import get_body +from ssapy.gravity import AccelHarmonic +from ssapy.utils import find_file + +DEGREE = 20 +T = 1.4e9 + + +def _read_egm(name): + egm = find_file(name, ext=".egm") + header = open(egm).read().splitlines() + radius = float(next(line for line in header if line.startswith("ModelRadius"))[11:]) + gm = float(next(line for line in header if line.startswith("ModelMass"))[9:]) + with open(egm + ".cof", "rb") as f: + f.read(8) + n_max = int.from_bytes(f.read(4), "little") + m_max = int.from_bytes(f.read(4), "little") + n_c = (m_max + 1) * (2 * n_max - m_max + 2) // 2 + n_s = m_max * (2 * n_max - m_max + 1) // 2 + c_values = np.frombuffer(f.read(8 * n_c)) + s_values = np.frombuffer(f.read(8 * n_s)) + c, s, k = {}, {}, 0 + for m in range(m_max + 1): + for n in range(m, n_max + 1): + c[(n, m)] = c_values[k] + k += 1 + k = 0 + for m in range(1, m_max + 1): + for n in range(m, n_max + 1): + s[(n, m)] = s_values[k] + k += 1 + return c, s, radius, gm + + +def _read_tab(name): + with open(find_file(name)) as f: + header = [float(x) for x in f.readline().replace(",", " ").split()] + c, s = {}, {} + for line in f: + fields = line.replace(",", " ").split() + n, m = int(fields[0]), int(fields[1]) + if n > DEGREE: + break + c[(n, m)], s[(n, m)] = float(fields[2]), float(fields[3]) + return c, s, header[0] * 1e3, header[1] * 1e9 + + +def _potential(r, c, s, radius, gm): + rr = np.linalg.norm(r) + x = r[2] / rr + lam = np.arctan2(r[1], r[0]) + total = 0.0 + for n in range(1, DEGREE + 1): + scale = (radius / rr) ** n + for m in range(n + 1): + norm = np.sqrt((2 - (m == 0)) * (2 * n + 1) * factorial(n - m) / factorial(n + m)) + pbar = norm * (-1) ** m * lpmv(m, n, x) # remove scipy's Condon-Shortley phase + total += scale * pbar * (c.get((n, m), 0.0) * np.cos(m * lam) + s.get((n, m), 0.0) * np.sin(m * lam)) + return gm / rr * total + + +def _gradient(r, *model, h=0.5): + return np.array([(_potential(r + h * e, *model) - _potential(r - h * e, *model)) / (2 * h) for e in np.eye(3)]) + + +@pytest.mark.parametrize("model", ["egm96", "egm2008"]) +def test_earth_harmonics_match_an_independent_potential(model): + earth = get_body("earth", model=model.upper()) + coefficients = _read_egm(model) + rotation = np.asarray(earth.orientation(T)) + for r_body in (np.array([7000e3, 1000e3, 500e3]), np.array([-2000e3, 3000e3, 6000e3])): + expected = _gradient(r_body, *coefficients) + got = rotation @ AccelHarmonic(earth, DEGREE, DEGREE)(rotation.T @ r_body, np.zeros(3), T) + assert np.linalg.norm(got - expected) / np.linalg.norm(expected) < 1e-8 + + +def test_moon_harmonics_match_an_independent_potential(): + moon = get_body("moon") + coefficients = _read_tab("gggrx_1200a_sha.tab") + rotation = np.asarray(moon.orientation(T)) + moon_position = np.asarray(moon.position(T)).reshape(3) + for r_body in (np.array([1900e3, 300e3, -200e3]), np.array([-500e3, 800e3, 1700e3])): + expected = _gradient(r_body, *coefficients) + got = rotation @ AccelHarmonic(moon, DEGREE, DEGREE)(moon_position + rotation.T @ r_body, np.zeros(3), T) + assert np.linalg.norm(got - expected) / np.linalg.norm(expected) < 1e-8 diff --git a/tests/test_io.py b/tests/test_io.py index 92fcf492..9c851f52 100644 --- a/tests/test_io.py +++ b/tests/test_io.py @@ -229,3 +229,37 @@ def put(start, text): pos = b3obs2pos("".join(line)) assert pos["satnum"] == 12345 assert pos["sensnum"] == 511 + + +def test_make_tle_writes_zero_drag_terms_that_sgp4_reads_back(): + # The TLE line-1 format puts ndot/2, nddot/6 and B* in columns 34-61. + # make_tle documents that dynamic terms are ignored, so sgp4 must read all + # three as 0 and the epoch to 1 ms; the checksum digit must match the + # mod-10 rule. B* was previously written as 0.99999 per Earth radius. + from astropy.time import Time + from sgp4.api import Satrec + + t = Time("2026-10-08T06:30:15.5", scale="utc") + line1, line2 = io.make_tle(7000e3, 0.001, np.radians(51.6), np.radians(30.0), np.radians(120.0), np.radians(75.0), t) + sat = Satrec.twoline2rv(line1, line2) + assert (sat.bstar, sat.ndot, sat.nddot) == (0.0, 0.0, 0.0) + assert abs((sat.jdsatepoch + sat.jdsatepochF) - t.utc.jd) * 86400.0 < 1e-3 + for line in (line1, line2): + digits = sum(int(c) if c.isdigit() else (1 if c == "-" else 0) for c in line[:68]) + assert line[68] == str(digits % 10) + + +def test_tle_from_a_tle_orbit_keeps_its_drag_terms(): + # Orbit.tle and SGP4Propagator(truncate=True) rebuild a TLE from Kozai + # elements; for an orbit read from a TLE they must keep its ndot, nddot and + # B* (here the ISS element set from the TLE format documentation, B* = + # -1.1606e-5), as sgp4 reads them, exactly. + import ssapy + from sgp4.api import Satrec + + line1 = "1 25544U 98067A 08264.51782528 -.00002182 00000-0 -11606-4 0 2927" + line2 = "2 25544 51.6416 247.4627 0006703 130.5360 325.0288 15.72125391563537" + orbit = ssapy.Orbit.fromTLETuple((line1, line2)) + rebuilt = Satrec.twoline2rv(*orbit.tle) + original = Satrec.twoline2rv(line1, line2) + assert (rebuilt.bstar, rebuilt.ndot, rebuilt.nddot) == (original.bstar, original.ndot, original.nddot) diff --git a/tests/test_leapfrog4_epoch.py b/tests/test_leapfrog4_epoch.py new file mode 100644 index 00000000..9a6c0b52 --- /dev/null +++ b/tests/test_leapfrog4_epoch.py @@ -0,0 +1,22 @@ +import numpy as np +import pytest + +import ssapy +from ssapy import AccelKepler + + +@pytest.mark.parametrize("t0", [0.0, 1.4e9, 3.8e9]) +def test_leapfrog4_is_fourth_order_at_real_gps_epochs(t0): + # Reference: the Keplerian propagator. Halving the step must cut the error + # by ~16x at any epoch (accept >= 12x) down to h = 5 s, where the error over + # 6000 s is under 0.05 m. Summing substep times onto a large GPS epoch used + # to stall the error near 2 m. + orbit = ssapy.Orbit.fromKeplerianElements(7500e3, 0.05, 0.9, 0.3, 1.2, 0.4, t=t0) + times = t0 + np.arange(0.0, 6000.0 + 1e-9, 600.0) + r_ref, _ = ssapy.rv(orbit, times, propagator=ssapy.KeplerianPropagator()) + errors = [] + for h in (20.0, 10.0, 5.0): + r, _ = ssapy.rv(orbit, times, propagator=ssapy.Leapfrog4Propagator(AccelKepler(), h=h)) + errors.append(np.max(np.linalg.norm(r - r_ref, axis=1))) + assert errors[0] / errors[1] > 12 and errors[1] / errors[2] > 12 + assert errors[2] < 0.05 diff --git a/tests/test_orbit.py b/tests/test_orbit.py index ad1f6c47..c0b7c478 100755 --- a/tests/test_orbit.py +++ b/tests/test_orbit.py @@ -18,13 +18,13 @@ from . import ssapy_test_helpers from .ssapy_test_helpers import timer, checkAngle, checkSphere, sample_orbit, sample_LEO_orbit, sample_GEO_orbit -from pathlib import Path -def _is_lfs_pointer(path): +def _has_ephemeris(): + # SSAPy's ephemeris comes from llnl-ssapy-data; skip only if it is missing. try: - with open(path, "rb") as f: - return f.readline().startswith(b"version https://git-lfs.github.com/spec/v1") - except Exception: + ssapy.body._planetary_ephemeris_path() + return True + except FileNotFoundError: return False try: @@ -1570,10 +1570,7 @@ def test_light_time_correction(): np.testing.assert_allclose(dircos, dc_linear, rtol=0, atol=1e-7) np.testing.assert_allclose(dircos, dc_exact, rtol=0, atol=1e-14) -@pytest.mark.skipif( - _is_lfs_pointer(Path(ssapy.datadir) / "de430.bsp"), - reason="de430.bsp is a Git LFS pointer, not real ephemeris data", -) +@pytest.mark.skipif(not _has_ephemeris(), reason="planetary ephemeris unavailable (llnl-ssapy-data not installed)") @pytest.mark.timeout(90) @timer def test_find_passes(): diff --git a/tests/test_orbit_solver.py b/tests/test_orbit_solver.py index 4420b874..f96875b3 100644 --- a/tests/test_orbit_solver.py +++ b/tests/test_orbit_solver.py @@ -312,7 +312,10 @@ def _getP(self): danchick.ell = 0.0 monkeypatch.setattr(ssapy.DanchickTwoPosOrbitSolver, 'X', staticmethod(lambda g: 1.0)) monkeypatch.setattr(ssapy.DanchickTwoPosOrbitSolver, 'dXdg', staticmethod(lambda g: 0.0)) - with pytest.raises(RuntimeError, match='Invalid x'): + # The x iteration now damps steps that would leave 0 < x < 1 and falls + # back to the eta iteration, so a non-converging case fails with a + # RuntimeError naming both attempts. + with pytest.raises(RuntimeError, match='Danchick two-position iteration failed'): ssapy.DanchickTwoPosOrbitSolver._getP(danchick) danchick.cos2f = np.nan diff --git a/tests/test_rk_interpolation.py b/tests/test_rk_interpolation.py new file mode 100644 index 00000000..2715aff8 --- /dev/null +++ b/tests/test_rk_interpolation.py @@ -0,0 +1,21 @@ +import numpy as np +import pytest + +import ssapy +from ssapy import AccelKepler + + +@pytest.mark.parametrize("h", [40.0, 20.0]) +def test_rk8_is_as_accurate_between_steps_as_at_them(h): + # Reference: the Keplerian propagator. RK8 with a fixed step is accurate + # to ~1e-8 m at its step boundaries over 6000 s in LEO; queried between + # boundaries it must stay within 1e-6 m in position and 1e-6 m/s in + # velocity (measured 1.1e-7 m at h = 40 s). The cubic spline used before + # gave 2.5e-2 m at h = 40 s. + t0 = 1.4e9 + orbit = ssapy.Orbit.fromKeplerianElements(7500e3, 0.05, 0.9, 0.3, 1.2, 0.4, t=t0) + times = t0 + np.array([0.0, 1234.5, 3333.3, 5999.9]) + r_ref, v_ref = ssapy.rv(orbit, times, propagator=ssapy.KeplerianPropagator()) + r, v = ssapy.rv(orbit, times, propagator=ssapy.RK8Propagator(AccelKepler(), h=h)) + assert np.max(np.linalg.norm(r - r_ref, axis=1)) < 1e-6 + assert np.max(np.linalg.norm(v - v_ref, axis=1)) < 1e-6 diff --git a/tests/test_sampler.py b/tests/test_sampler.py index 7493cc8d..1bd5757a 100644 --- a/tests/test_sampler.py +++ b/tests/test_sampler.py @@ -930,3 +930,25 @@ def test_mh_sampler(): import subprocess cmd = "gprof2dot -f pstats {} -n1 -e1 | dot -Tpng -o {}".format(args.prof_out, args.prof_png) subprocess.run(cmd, shell=True) + + +def test_circular_guess_with_rates_uses_the_selected_row(): + # With 'pmra'/'pmdec' columns only one observation is needed; the one + # used must be the first of `indices` (row 0 used to be taken regardless). + orbit, arc = _circular_guess_arc() + rows = [] + for k in range(2): + ra, dec, _, pmra, pmdec, _ = ssapy.radec( + orbit, arc['time'][k].gps, obsPos=arc['rStation_GCRF'][k].to(u.m).value, + obsVel=arc['vStation_GCRF'][k].to(u.m / u.s).value, rate=True) + rows.append((ra, dec, pmra, pmdec)) + arcpm = QTable() + arcpm['ra'] = [r[0] for r in rows] * u.rad + arcpm['dec'] = [r[1] for r in rows] * u.rad + arcpm['time'] = arc['time'][:2] + arcpm['rStation_GCRF'] = arc['rStation_GCRF'][:2] + arcpm['vStation_GCRF'] = arc['vStation_GCRF'][:2] + arcpm['pmra'] = [r[2] for r in rows] * u.rad / u.s + arcpm['pmdec'] = [r[3] for r in rows] * u.rad / u.s + _state, epoch = ssapy.circular_guess(arcpm, indices=[1, 0]) + assert abs((epoch - arcpm['time'][1]).to(u.s).value) < 1e-6 diff --git a/tests/test_sun_position_fast.py b/tests/test_sun_position_fast.py new file mode 100644 index 00000000..5babdd23 --- /dev/null +++ b/tests/test_sun_position_fast.py @@ -0,0 +1,24 @@ +import numpy as np +import pytest +import astropy.units as u +from astropy.coordinates import get_body_barycentric +from astropy.time import Time + +from ssapy import utils + + +def test_fast_sun_position_stays_within_an_arcminute_from_1980_to_2060(): + # Reference: astropy's geometric geocentric Sun (built-in ephemeris). + # Montenbruck & Gill's series with the Earth-Moon barycentre's perihelion + # advance (0.32327364 deg/century, JPL approximate elements) stays within + # 1 arcmin in direction (measured 0.58 arcmin) and 3e-4 in distance over + # 1980-2060; without the advance it reached 11.6 arcmin. + worst = 0.0 + for jd in np.linspace(Time("1980-01-01", scale="tt").jd, Time("2060-01-01", scale="tt").jd, 200): + t = Time(jd, format="jd", scale="tt") + expected = (get_body_barycentric("sun", t) - get_body_barycentric("earth", t)).xyz.to_value(u.m) + position = np.asarray(utils.sunPos(t.gps, fast=True), dtype=float).reshape(3) + cosine = np.dot(position, expected) / (np.linalg.norm(position) * np.linalg.norm(expected)) + worst = max(worst, np.degrees(np.arccos(np.clip(cosine, -1.0, 1.0))) * 60.0) + assert np.linalg.norm(position) == pytest.approx(np.linalg.norm(expected), rel=3e-4) + assert worst < 1.0 diff --git a/tests/test_tangent_plane_rates.py b/tests/test_tangent_plane_rates.py new file mode 100644 index 00000000..7b389993 --- /dev/null +++ b/tests/test_tangent_plane_rates.py @@ -0,0 +1,20 @@ +import numpy as np +import pytest + +from ssapy import utils + + +@pytest.mark.parametrize("offset", [0.0, 0.05, 0.33, 0.8]) +def test_tangent_plane_rates_are_the_derivative_of_the_projection(offset): + # The projected rates must be the time derivative of the projected + # position: a central difference of lb_to_tan for a point moving at + # (mul, mub) on the sphere (mul = cos(b) dl/dt) agrees to 1e-8. + lcen, bcen = np.array([1.0]), np.array([0.4]) + lb, b = np.array([1.0 + offset]), np.array([0.4 - 0.6 * offset]) + mul, mub = np.array([0.01]), np.array([-0.02]) + _x, _y, vx, vy = utils.lb_to_tan(lb, b, mul=mul, mub=mub, lcen=lcen, bcen=bcen) + h = 1e-6 + xp, yp = utils.lb_to_tan(lb + mul / np.cos(b) * h, b + mub * h, lcen=lcen, bcen=bcen) + xm, ym = utils.lb_to_tan(lb - mul / np.cos(b) * h, b - mub * h, lcen=lcen, bcen=bcen) + np.testing.assert_allclose(vx, (xp - xm) / (2 * h), rtol=1e-8) + np.testing.assert_allclose(vy, (yp - ym) / (2 * h), rtol=1e-8) diff --git a/tests/test_tle_drag.py b/tests/test_tle_drag.py index 066c00e5..18e102ca 100644 --- a/tests/test_tle_drag.py +++ b/tests/test_tle_drag.py @@ -1,5 +1,4 @@ import warnings -from pathlib import Path import numpy as np import pytest @@ -16,22 +15,31 @@ _sgp4_reference_arc, fit_drag) -def _is_lfs_pointer(path): +# The numerical drag propagator needs the gravity and ephemeris data from +# llnl-ssapy-data; skip those tests if it is not installed. +def _has_data(name): + # SSAPy's data comes from llnl-ssapy-data; skip only if it is missing. try: - with open(path, "rb") as f: - return f.readline().startswith(b"version https://git-lfs.github.com/spec/v1") - except Exception: + ssapy.utils.find_file(name) + return True + except FileNotFoundError: return False -# The numerical drag propagator needs the (Git LFS) gravity + ephemeris data, -# which is absent in CI; skip those tests there, as tests/test_accel.py does. -HAS_EGM84 = not _is_lfs_pointer(Path(ssapy.datadir) / "egm84.egm.cof") -HAS_DE430 = not _is_lfs_pointer(Path(ssapy.datadir) / "de430.bsp") -HAS_MOON_PA = not _is_lfs_pointer(Path(ssapy.datadir) / "moon_pa_de440_200625.bpc") -HAS_BODY_DATA = HAS_EGM84 and HAS_DE430 and HAS_MOON_PA +def _has_ephemeris(): + try: + ssapy.body._planetary_ephemeris_path() + return True + except FileNotFoundError: + return False + + +HAS_EGM84 = _has_data("egm84.egm.cof") +HAS_EPHEMERIS = _has_ephemeris() +HAS_MOON_PA = _has_data("moon_pa_de440_200625.bpc") +HAS_BODY_DATA = HAS_EGM84 and HAS_EPHEMERIS and HAS_MOON_PA needs_body_data = pytest.mark.skipif( - not HAS_BODY_DATA, reason="body data unavailable (Git LFS pointer)") + not HAS_BODY_DATA, reason="body data unavailable (llnl-ssapy-data not installed)") ISS = ("1 25544U 98067A 24015.54791435 .00016717 00000-0 30074-3 0 9993", "2 25544 51.6416 247.4627 0006703 130.5360 325.0288 15.49514637123456") diff --git a/tests/test_tracks.py b/tests/test_tracks.py index 709699b7..18f90f5c 100644 --- a/tests/test_tracks.py +++ b/tests/test_tracks.py @@ -1099,3 +1099,16 @@ def run(self, **kwargs): ) assert created['kwargs']['approximate'] is True assert newmht.run_kwargs == {'first': 2, 'last': 5} + + +def test_volume_distance_prior_uses_the_requested_scale(): + # lnp(d) = 2 ln(d/s) - d/s - (2 ln 2 - 2 + 1e-9) peaks at d = 2 s with + # lnp = -1e-9 (closed form, 1e-12). With s = 1000 km the peak is at 2000 km; + # the scale argument used to be ignored in favour of RGEO. + scale = 1000e3 + prior = VolumeDistancePrior(scale=scale) + orbit = sample_LEO_orbit(t=0) + assert prior(orbit, 2 * scale) == pytest.approx(-1e-9, abs=1e-12) + d = 500e3 + expected = 2 * np.log(d / scale) - d / scale - (2 * np.log(2) - 2 + 1e-9) + assert prior(orbit, d) == pytest.approx(expected, abs=1e-12) diff --git a/tests/test_two_position_solver_references.py b/tests/test_two_position_solver_references.py new file mode 100644 index 00000000..b741acef --- /dev/null +++ b/tests/test_two_position_solver_references.py @@ -0,0 +1,60 @@ +import numpy as np +import pytest + +import ssapy + +T1 = 1.4e9 +CASES = [ # (a [m], e, i [rad], arc [s]) + (7000e3, 0.01, 0.9, 600.0), # LEO, 38 deg + (7000e3, 0.01, 0.9, 4400.0), # LEO, 270 deg (long way) + (42164e3, 0.001, 0.1, 21600.0), # GEO, 90 deg + (20000e3, 0.4, 1.0, 9000.0), # eccentric MEO, 134 deg + (24400e3, 0.73, 0.5, 20000.0), # GTO, 166 deg +] + + +def _endpoints(a, e, inc, arc): + orbit = ssapy.Orbit.fromKeplerianElements(a, e, inc, 0.5, 1.0, 0.3, t=T1) + r1, v1 = ssapy.rv(orbit, T1) + r2, _ = ssapy.rv(orbit, T1 + arc) + kappa_sign = 1 if np.dot(np.cross(r1, r2), np.cross(r1, v1)) > 0 else -1 + return r1, r2, v1, kappa_sign + + +@pytest.mark.parametrize("a, e, inc, arc", CASES) +def test_shefer_recovers_the_keplerian_velocity_for_short_and_long_arcs(a, e, inc, arc): + # Reference: the Keplerian orbit that produced both positions. Shefer's + # method recovers the departure velocity to 1e-9 m/s on every arc + # (measured <= 6e-12 m/s). + r1, r2, v1, kappa_sign = _endpoints(a, e, inc, arc) + solved = ssapy.SheferTwoPosOrbitSolver(r1, r2, T1, T1 + arc, kappaSign=kappa_sign).solve() + assert np.linalg.norm(solved.v - v1) < 1e-9 + + +def test_gauss_warns_instead_of_returning_a_wrong_long_arc_orbit(): + # Gauss's iteration recovers a 38 deg LEO arc to 1e-9 m/s without a + # warning, but diverges for a quarter GEO orbit; it must now warn there + # (it used to return a velocity 2.7 km/s off silently). + r1, r2, v1, kappa_sign = _endpoints(*CASES[0]) + solved = ssapy.GaussTwoPosOrbitSolver(r1, r2, T1, T1 + CASES[0][3], kappaSign=kappa_sign).solve() + assert np.linalg.norm(solved.v - v1) < 1e-9 + r1, r2, v1, kappa_sign = _endpoints(*CASES[2]) + with pytest.warns(RuntimeWarning, match="did not converge"): + ssapy.GaussTwoPosOrbitSolver(r1, r2, T1, T1 + CASES[2][3], kappaSign=kappa_sign).solve() + + +@pytest.mark.parametrize("a, e", [(7000e3, 0.01), (26560e3, 0.3), (42164e3, 0.001), (24400e3, 0.73)]) +def test_danchick_and_shefer_cover_the_whole_revolution(a, e): + # Reference: Keplerian truth. Over arcs from 2% to 98% of a period (both + # short- and long-way, including fast perigee passages on a GTO) both + # solvers recover the departure velocity to 1e-7 m/s (measured <= 3.6e-9). + # Danchick's method used to raise "Invalid x" on 19 of these 100 arcs. + orbit = ssapy.Orbit.fromKeplerianElements(a, e, 0.9, 0.5, 1.0, 0.3, t=T1) + for fraction in np.linspace(0.02, 0.98, 25): + arc = fraction * orbit.period + r1, v1 = ssapy.rv(orbit, T1) + r2, _ = ssapy.rv(orbit, T1 + arc) + kappa_sign = 1 if np.dot(np.cross(r1, r2), np.cross(r1, v1)) > 0 else -1 + for solver in (ssapy.DanchickTwoPosOrbitSolver, ssapy.SheferTwoPosOrbitSolver): + solved = solver(r1, r2, T1, T1 + arc, kappaSign=kappa_sign).solve() + assert np.linalg.norm(solved.v - v1) < 1e-7, (solver.__name__, fraction) diff --git a/tests/test_utils.py b/tests/test_utils.py index 816d2ae3..9e7c2301 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -184,8 +184,9 @@ def test_small_coordinate_and_angle_helpers(capsys): assert utils.zenithangle_to_altitude(60.0) == 30.0 assert np.isclose(utils.altitude_to_zenithangle(np.pi / 6, deg=False), np.pi / 3) assert np.isclose(utils.zenithangle_to_altitude(np.pi / 3, deg=False), np.pi / 6) - assert utils.rightascension_to_hourangle(30.0, 2.0) == "0:0:0" - assert utils.rightascension_to_hourangle("02:00:00", "01:00:00").count(":") == 2 + # Numeric arguments are degrees (docstring): LST 2 deg - RA 30 deg = -28 deg = 22h08m. + assert utils.rightascension_to_hourangle(30.0, 2.0) == "22:8:00" + assert utils.rightascension_to_hourangle("02:00:00", "01:00:00") == "23:0:0" assert utils.dms_to_dd("12:30:00") == 12.5 assert utils.dms_to_dd(["12:30:00", "-12:30:00"]) == [12.5, -12.5] @@ -193,7 +194,6 @@ def test_small_coordinate_and_angle_helpers(capsys): assert utils.hms_to_dd("12:00:00") == 180.0 assert utils.hms_to_dd(["12:00:00", "01:30:00"]) == [180.0, 22.5] assert utils.dd_to_hms(-180.0) == "12:0:0" - assert "cannot be negative" in capsys.readouterr().out def test_sun_ra_dec_matches_astropy_solar_position(): @@ -469,17 +469,17 @@ def test_coordinate_angle_edge_cases_and_errors(capsys): assert np.isfinite(az) assert np.isfinite(alt) - az2, alt2 = utils.equatorial_to_horizontal( - observer_latitude=45.0, - declination=0.0, - right_ascension="01:00:00", - hour_angle="00:00:00", - local_time="02:00:00", - hms=True, - ) + with pytest.warns(UserWarning, match="using hour_angle"): + az2, alt2 = utils.equatorial_to_horizontal( + observer_latitude=45.0, + declination=0.0, + right_ascension="01:00:00", + hour_angle="00:00:00", + local_time="02:00:00", + hms=True, + ) assert np.isfinite(az2) assert np.isfinite(alt2) - assert "Using hour_angle" in capsys.readouterr().out az3, alt3 = utils.equatorial_to_horizontal( observer_latitude=45.0,