diff --git a/.github/workflows/testsuite.yaml b/.github/workflows/testsuite.yaml index 3fadf55..fb28529 100644 --- a/.github/workflows/testsuite.yaml +++ b/.github/workflows/testsuite.yaml @@ -27,8 +27,12 @@ jobs: pre-commit install pre-commit run -a datacache: - name: Cache Data - runs-on: ubuntu-latest + name: Cache Data (${{ matrix.os }}) + strategy: + fail-fast: false + matrix: + os: [ubuntu-latest, macos-latest] + runs-on: ${{ matrix.os }} steps: - uses: actions/checkout@v3 with: @@ -37,19 +41,17 @@ jobs: uses: actions/setup-python@v4 with: python-version: 3.9 - - name: Get current date - id: date - run: echo "::set-output name=date::$(date +'%Y-%m-%d')" - name: cache-kernels + id: cache-kernels uses: actions/cache@v3 with: - path: ~/ap_cache - key: kernel_cache-${{ steps.date.outputs.date }} - - name: Download and cache - run: | - pip install astropy - python ./ci/cache_kernels.py save + path: ~/.astropy/cache/download + key: kernel-de430-${{ runner.os }} + - name: Download kernels if: steps.cache-kernels.outputs.cache-hit != 'true' + run: | + pip install . + python -c "from lunarsky.spice_utils import download_big_kernels; download_big_kernels(show_progress=False)" tests: env: @@ -72,18 +74,11 @@ jobs: uses: actions/setup-python@v4 with: python-version: ${{ matrix.python-version }} - - name: Get current date - id: date - run: echo "::set-output name=date::$(date +'%Y-%m-%d')" - name: load-cache uses: actions/cache@v3 with: - path: ~/ap_cache - key: kernel_cache-${{ steps.date.outputs.date }} - - name: load kernels - run: | - pip install astropy - python ./ci/cache_kernels.py load + path: ~/.astropy/cache/download + key: kernel-de430-${{ runner.os }} - name: Install run: | pip install .[test] pytest-cov diff --git a/CHANGELOG.md b/CHANGELOG.md index c88945a..e7a3faf 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,14 @@ # Changelog +## [Unreleased] + +## Changed +- spice_utils now has a PCK file reader to read the moon's planetary constants +- Transforms are now defined in Python, not using SPICE + +## Deprecated +- Removed dependency on spiceypy for coordinate definitions. Kernels are now handled with jplephem + ## [0.2.6] -- 2024-12-12 ## Fixed - Avoid deprecation warning with numpy.broadcast_shape diff --git a/README.md b/README.md index ee96430..f7ef392 100644 --- a/README.md +++ b/README.md @@ -9,15 +9,13 @@ An extension to `astropy`, providing selenocentric and topocentric reference fra for the Moon and transformations of star positions from the ICRS system to these frames. This is to describe the sky as observed from the surface of the Moon. -Non-relativistic transformations are calculated using the SPICE toolkit. Relativistic -corrections (stellar aberration) will be added. +Non-relativistic transformations are calculated using lunar orientation data from JPL SPICE kernel files, read with `jplephem`. ## Dependencies * `numpy` -* `astropy>=3.0` +* `astropy>=6.0.0` * `jplephem` -* `spiceypy` ## Installation @@ -33,6 +31,14 @@ git clone https://github.com/aelanman/lunarsky python setup.py install ``` +The first time lunarsky is used, it will automatically download the JPL DE430 ephemeris +(~150 MB) into astropy's data cache (`~/.astropy/cache/`). Subsequent uses read from the +cache. To pre-fetch the kernel (e.g. before going offline, or in CI), run: +``` +from lunarsky import spice_utils +spice_utils.download_big_kernels() +``` + ## Usage ![mcmf_coords](./docs/figure.png) @@ -53,7 +59,7 @@ The cartesian axes of the selenocentric system are those of the MCMF frame. Sele ## Credit -This package makes use of the ``spiceypy`` wrapper [2] for the JPL SPICE Toolkit, produced by the NASA Navigation and Ancillary Information Facility (NAIF) [3] [4]. The transformations are defined using data in kernel files ``pck/moon_pa_de421_1900-2050.bpc``, ``moon_080317.tf``, and ``moon_assoc_me.tf``. These may be found at [the NAIF website](https://naif.jpl.nasa.gov/pub/naif/generic_kernels), and were produced by Nat Bachman (NAIF/JPL) in March 2008. Further information may be found in the comments in these files in the `data` directory. +This package uses lunar orientation data from SPICE kernel files produced by the NASA Navigation and Ancillary Information Facility (NAIF) [3] [4]. The transformations are defined using data in kernel files ``pck/moon_pa_de421_1900-2050.bpc``, ``moon_080317.tf``, and ``moon_assoc_me.tf``. These may be found at [the NAIF website](https://naif.jpl.nasa.gov/pub/naif/generic_kernels), and were produced by Nat Bachman (NAIF/JPL) in March 2008. Binary kernel files are read using ``jplephem`` [5]. Further information may be found in the comments in the kernel files in the `data` directory. ## References [1]: Ye, Hanlin, et al. "Looking Vector Direction Analysis for the Moon-Based Earth Observation Optical Sensor." IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 11, no. 11, Nov. 2018, pp. 4488–99. IEEE Xplore, doi:10.1109/JSTARS.2018.2870247. @@ -63,3 +69,5 @@ This package makes use of the ``spiceypy`` wrapper [2] for the JPL SPICE Toolkit [3]: Acton, C.H.; "Ancillary Data Services of NASA's Navigation and Ancillary Information Facility;" Planetary and Space Science, Vol. 44, No. 1, pp. 65-70, 1996. [4]: Charles Acton, Nathaniel Bachman, Boris Semenov, Edward Wright; A look toward the future in the handling of space science mission geometry; Planetary and Space Science (2017); DOI 10.1016/j.pss.2017.02.013; https://doi.org/10.1016/j.pss.2017.02.013 + +[5]: Rhodes, Brandon. jplephem: Python module for evaluating JPL ephemerides. https://github.com/brandon-rhodes/python-jplephem diff --git a/ci/cache_kernels.py b/ci/cache_kernels.py deleted file mode 100644 index 2241661..0000000 --- a/ci/cache_kernels.py +++ /dev/null @@ -1,27 +0,0 @@ -from astropy.utils.data import ( - download_files_in_parallel, - export_download_cache, - import_download_cache, -) -import os -import sys - -home = os.path.expanduser("~") -path = os.path.join(home, "ap_cache") -cache_file = os.path.join(path, "astropy_cache.zip") -os.makedirs(path, exist_ok=True) - -if sys.argv[1] == "save": - if not os.path.exists(cache_file): - knames = [ - "lsk/naif0012.tls", - "spk/planets/de430.bsp", - "pck/earth_latest_high_prec.bpc", - ] - _naif_kernel_url = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/" - kurls = [_naif_kernel_url + kname for kname in knames] - paths = download_files_in_parallel(kurls, cache=True, show_progress=False) - export_download_cache(cache_file, urls=kurls, overwrite=True) - -if sys.argv[1] == "load": - import_download_cache(cache_file) diff --git a/lunarsky/__init__.py b/lunarsky/__init__.py index 35c308e..9eeef94 100644 --- a/lunarsky/__init__.py +++ b/lunarsky/__init__.py @@ -1,11 +1,3 @@ -import sys - -if sys.platform == "darwin": - # On MacOS, need to set multiprocessing to fork by default - import multiprocessing - - multiprocessing.set_start_method("fork", force=True) - from .moon import * # noqa from .mcmf import * # noqa from .topo import * # noqa diff --git a/lunarsky/mcmf.py b/lunarsky/mcmf.py index e462f6a..9371b10 100644 --- a/lunarsky/mcmf.py +++ b/lunarsky/mcmf.py @@ -17,8 +17,6 @@ from astropy.coordinates.transformations import FunctionTransformWithFiniteDifference from astropy.coordinates.builtin_frames.icrs import ICRS -import spiceypy as spice - __all__ = ["MCMF"] _J2000 = Time("J2000", scale="tt") @@ -56,6 +54,7 @@ def moon_location(self): def make_transform(coo, toframe): + from .spice_utils import j2000_to_moon_me, body_position ap_to_spice = {"icrs": ("J2000", 0), "mcmf": ("MOON_ME", 301)} @@ -67,12 +66,12 @@ def make_transform(coo, toframe): obstime = toframe.obstime if to_mcmf else coo.obstime # Make arrays - ets = np.atleast_1d((obstime - _J2000).sec) + ets = np.atleast_1d((obstime.tdb - _J2000).sec) shape_out = np.broadcast_shapes(coo.shape, ets.shape) coo_cart = coo.cartesian - mats = np.stack([spice.pxform("J2000", "MOON_ME", et) for et in ets], axis=0) + mats = j2000_to_moon_me(ets) if not to_mcmf: mats = np.linalg.inv(mats) @@ -85,11 +84,8 @@ def make_transform(coo, toframe): # If not unitspherical, shift by origin vector before rotating. if not is_unitspherical: # Make origin vector(s) in coo's frame. - orig_posvel = ( - np.asarray([spice.spkgeo(to_id, et, from_name, from_id)[0] for et in ets]) - * un.km - ) - coo_cart -= CartesianRepresentation((orig_posvel.T)[:3]) + orig_pos = body_position(to_id, ets, from_name, from_id) * un.km + coo_cart -= CartesianRepresentation(orig_pos.T) newrepr = coo_cart.transform(mats).reshape(shape_out) diff --git a/lunarsky/moon.py b/lunarsky/moon.py index da0b1ec..fa7ffe8 100644 --- a/lunarsky/moon.py +++ b/lunarsky/moon.py @@ -12,22 +12,15 @@ from astropy.coordinates.attributes import Attribute from astropy.utils import minversion -from .spice_utils import remove_topo - COPY_IF_NEEDED = None if minversion(np, "2.0") else False -LUNAR_RADIUS = 1737.1e3 # m - __all__ = ["MoonLocation", "MoonLocationAttribute"] class SPHERESelenodeticRepresentation(BaseGeodeticRepresentation): - """Lunar ellipsoid as a sphere - - Radius defined by lunarsky.spice_utils.LUNAR_RADIUS - """ + """Lunar ellipsoid as a sphere""" - _equatorial_radius = LUNAR_RADIUS * u.m + _equatorial_radius = 1737.1e3 * u.m _flattening = 0.0 @@ -231,16 +224,6 @@ class MoonLocation(u.Quantity): _location_dtype = np.dtype({"names": ["x", "y", "z"], "formats": [np.float64] * 3}) _array_dtype = np.dtype((np.float64, (3,))) - # Manage the set of defined ephemerides. - # Class attributes only - _inuse_stat_ids = [] - _avail_stat_ids = None - _existing_locs = [] - _ref_count = [] - - # This instance's station id(s) - station_ids = [] - info = MoonLocationInfo() def __new__(cls, *args, **kwargs): @@ -284,56 +267,6 @@ def of_site(cls, site_name: str): lon_lat = sites[site_name.lower()] return cls.from_selenodetic(*lon_lat, ellipsoid="IAU2000") - @classmethod - def _set_site_id(cls, inst): - """ - Set the station ID number and manage registry of defined stations. - """ - if inst.isscalar: - llh_arr = [ - ( - inst.lon.deg.item(), - inst.lat.deg.item(), - inst.height.to_value("km").item(), - ) - ] - ncrds = 1 - else: - llh_arr = zip(inst.lon.deg, inst.lat.deg, inst.height.to_value("km")) - ncrds = inst.lon.size - - statids = [] - if cls._avail_stat_ids is None: - cls._avail_stat_ids = list(range(999, 0, -1)) - - if len(cls._avail_stat_ids) < ncrds: - raise ValueError("Too many unique MoonLocation objects open at once.") - - for llh in llh_arr: - lonlatheight = "_".join( - ["{:.4f}".format(ll) for ll in llh] + [inst._ellipsoid] - ) - if lonlatheight not in cls._existing_locs: - new_stat_id = cls._avail_stat_ids.pop() - cls._existing_locs.append(lonlatheight) - statids.append(new_stat_id) - cls._inuse_stat_ids.append(new_stat_id) - cls._ref_count.append(1) - else: - ind = cls._existing_locs.index(lonlatheight) - cls._ref_count[ind] += 1 - statids.append(cls._inuse_stat_ids[ind]) - inst.station_ids = statids - return inst - - def _set_station_id(self): - """ - Run classmethod for setting station IDs. - - Convenience function used for testing mostly - """ - self.__class__._set_site_id(self) - @classmethod def from_selenocentric(cls, x, y, z, unit=None): """ @@ -457,35 +390,6 @@ def from_selenodetic(cls, lon, lat, height=0.0, ellipsoid=None): def __str__(self): return self.__repr__() - def copy(self): - # Necessary to preserve station_ids list - return self.__copy__() - - def __copy__(self): - # Ensure that the station_ids are copied as well under shallow copy - obj = super().copy() - obj.station_ids = self.station_ids - return obj - - def __del__(self): - # Remove this MoonLocation's station_ids from _inuse_stat_ids and - # locations from _existing_locs. - # Also clear the corresponding frames from spice variable pool. - for si, stat_id in enumerate(self.station_ids): - try: - ind = self.__class__._inuse_stat_ids.index(stat_id) - except ValueError: - continue - count = self.__class__._ref_count[ind] - if count <= 1: - sid = self.__class__._inuse_stat_ids.pop(ind) - self.__class__._existing_locs.pop(ind) - self.__class__._avail_stat_ids.insert(0, sid) - remove_topo(stat_id) - self.__class__._ref_count.pop(ind) - else: - self.__class__._ref_count[ind] -= 1 - @property def selenodetic(self): """Convert to selenodetic coordinates.""" @@ -532,7 +436,7 @@ def lon(self): @property def lat(self): - """Longitude of the location""" + """Latitude of the location""" return self.selenodetic[1] @property @@ -573,7 +477,7 @@ def get_mcmf(self, obstime=None): mcmf = property( get_mcmf, - doc="""An `~astropy.coordinates.MCMF` object with + doc="""An `~astropy.coordinates.MCMF` object for the location of this object at the default ``obstime``.""", ) diff --git a/lunarsky/spice_utils.py b/lunarsky/spice_utils.py index 254f2d5..d312ab4 100644 --- a/lunarsky/spice_utils.py +++ b/lunarsky/spice_utils.py @@ -1,230 +1,323 @@ import numpy as np import os -import tempfile -from astropy.utils.data import download_files_in_parallel +from astropy.utils.data import download_file from astropy.coordinates.matrix_utilities import rotation_matrix from astropy.time import Time import astropy.units as unit -import spiceypy as spice - from .mcmf import MCMF from .data import DATA_PATH _J2000 = Time("J2000") -TEMPORARY_KERNEL_DIR = tempfile.TemporaryDirectory() +# --------------------- +# Binary PCK reader +# --------------------- + +_PCK_DATA = None -def check_is_loaded(search): +def _read_bpc(filepath): """ - Search the kernel pool variable names for a given string. + Read a DAF-format binary PCK file (Type 2 Chebyshev). + + Returns a dict with segment metadata and coefficient arrays. """ - try: - spice.gnpool(search, 0, 100) - except (spice.support_types.SpiceyError): - return False - return True + from jplephem.daf import DAF + with open(filepath, "rb") as f: + daf = DAF(f) -def list_kernels(): - """ - List loaded kernels. + segments = list(daf.summaries()) + if len(segments) != 1: + raise ValueError(f"Expected 1 segment in binary PCK, got {len(segments)}") - Returns - ------- - list of str - Kernel names (file paths) - list of str - Corresponding kernel types - """ - knames, ktypes = [], [] - for typ in ["spk", "fk", "tk", "pck", "lsk"]: - for ii in range(spice.ktotal(typ)): - dat = spice.kdata(ii, typ) - knames.append(dat[0]) - ktypes.append(dat[1]) - return knames, ktypes - - -def furnish_kernels(): - kernel_names = [ - "pck/moon_pa_de421_1900-2050.bpc", - "fk/satellites/moon_080317.tf", - "fk/satellites/moon_assoc_me.tf", - ] - kernel_paths = [os.path.join(DATA_PATH, kn) for kn in kernel_names] - for kp in kernel_paths: - spice.furnsh(kp) + name, values = segments[0] + data_type = int(values[4]) + start_addr = int(values[5]) + end_addr = int(values[6]) - # LSK and DE430 Kernels - knames = ["lsk/naif0012.tls", "spk/planets/de430.bsp"] - _naif_kernel_url = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/" - kurls = [_naif_kernel_url + kname for kname in knames] - paths = download_files_in_parallel(kurls, cache=True, show_progress=False) - for kp in paths: - kernel_paths.append(kp) - spice.furnsh(kp) + if data_type != 2: + raise ValueError(f"Unsupported binary PCK data type {data_type}") - return kernel_paths + arr = daf.read_array(start_addr, end_addr) + init_epoch = arr[-4] + interval = arr[-3] + rsize = int(arr[-2]) + n_records = int(arr[-1]) + n_coeffs = (rsize - 2) // 3 -def lunar_surface_ephem(pos_x, pos_y, pos_z, station_num=98): - """ - Make an SPK for the point on the lunar surface + data = arr[: n_records * rsize].reshape(n_records, rsize) + + return { + "init_epoch": init_epoch, + "interval": interval, + "n_coeffs": n_coeffs, + "n_records": n_records, + "data": data, + } + + +def _ensure_pck(): + global _PCK_DATA + if _PCK_DATA is None: + _PCK_DATA = _read_bpc( + os.path.join(DATA_PATH, "pck", "moon_pa_de421_1900-2050.bpc") + ) + return _PCK_DATA + + +def _cheby_eval(coeffs, tau): + """Evaluate Chebyshev polynomial at tau in [-1, 1]. coeffs shape: (..., n_coeffs).""" + n = coeffs.shape[-1] + if n == 0: + return np.zeros(coeffs.shape[:-1]) + T_prev = np.ones(coeffs.shape[:-1]) + if n == 1: + return coeffs[..., 0] * T_prev + T_curr = tau + result = coeffs[..., 0] * T_prev + coeffs[..., 1] * T_curr + for i in range(2, n): + T_prev, T_curr = T_curr, 2 * tau * T_curr - T_prev + result = result + coeffs[..., i] * T_curr + return result + + +def _R1(a): + c, s = np.cos(a), np.sin(a) + return np.array([[1, 0, 0], [0, c, s], [0, -s, c]]) + + +def _R2(a): + c, s = np.cos(a), np.sin(a) + return np.array([[c, 0, -s], [0, 1, 0], [s, 0, c]]) + + +def _R3(a): + c, s = np.cos(a), np.sin(a) + return np.array([[c, s, 0], [-s, c, 0], [0, 0, 1]]) + + +# PA to ME constant rotation from moon_080317.tf (FRAME_31007 = MOON_ME_DE421). +# These angles are specific to DE421; they will differ for other ephemerides. +# TKFRAME_31007_ANGLES = (67.92, 78.56, 0.30) arcseconds, AXES = (3, 2, 1) +# SPICE convention: R_PA_to_ME = R1(-a3) @ R2(-a2) @ R3(-a1) +_a1 = np.radians(67.92 / 3600) +_a2 = np.radians(78.56 / 3600) +_a3 = np.radians(0.30 / 3600) +_PA_TO_ME = _R1(-_a3) @ _R2(-_a2) @ _R3(-_a1) - Creates a temporary file and furnishes from that. + +def j2000_to_moon_me(ets): + """ + Compute J2000 -> MOON_ME rotation matrices. Parameters ---------- - pos_x, pos_y, pos_z: float - MCMF frame cartesian position in km - station_num: int - Station number + ets : float or array_like + TDB seconds past J2000. Returns ------- - int: - Ephemeris ID number + mats : ndarray, shape (..., 3, 3) """ - point_id = 301000 + station_num + pck = _ensure_pck() + ets = np.atleast_1d(np.asarray(ets, dtype=float)) - ets = np.array([spice.str2et("1950-01-01"), spice.str2et("2150-01-01")]) + record_idx = ((ets - pck["init_epoch"]) / pck["interval"]).astype(int) + record_idx = np.clip(record_idx, 0, pck["n_records"] - 1) - states = np.zeros((len(ets), 6)) - states[:, :3] = np.repeat([[pos_x], [pos_y], [pos_z]], len(ets), axis=1).T + records = pck["data"][record_idx] + mid = records[:, 0] + half = records[:, 1] + nc = pck["n_coeffs"] - center = 301 - frame = "MOON_ME" - degree = 1 + tau = (ets - mid) / half - fname = os.path.join(TEMPORARY_KERNEL_DIR.name, "lunar_points.bsp") - if os.path.exists(fname): - spice.unload(fname) - handle = spice.spkopa(fname) - else: - handle = spice.spkopn(fname, "SPK_FILE", 0) - spice.spkw09( - handle, - point_id, - center, - frame, - ets[0], - ets[-1], - "0", - degree, - len(ets), - states.tolist(), - ets.tolist(), - ) - spice.spkcls(handle) - spice.furnsh(fname) - - return point_id - - -def topo_frame_def(latitude, longitude, station_num=98, moon=True): - """ - Make a list of strings defining a topocentric frame. This can then be loaded - with spiceypy.lmpool. - """ - if moon: - idnum = 1301000 - station_name = f"LUNAR-TOPO-{station_num}" - relative = "MOON_ME" - else: - # Used in tests only - idnum = 1399000 - station_name = f"EARTH-TOPO-{station_num}" - relative = "ITRF93" + ra_coeffs = records[:, 2 : 2 + nc] + dec_coeffs = records[:, 2 + nc : 2 + 2 * nc] + w_coeffs = records[:, 2 + 2 * nc : 2 + 3 * nc] - # The DSS stations are built into SPICE, and they number up to 66. - # We will call this station number 98. - idnum += station_num - fm_center_id = idnum - 1000000 + ra = _cheby_eval(ra_coeffs, tau) + dec = _cheby_eval(dec_coeffs, tau) + w = _cheby_eval(w_coeffs, tau) - ecef_to_enu = np.matmul( - rotation_matrix(-longitude, "z", unit="deg"), - rotation_matrix(latitude, "y", unit="deg"), - ).T - # Reorder the axes so that X,Y,Z = E,N,U - ecef_to_enu = ecef_to_enu[[2, 1, 0]] + # J2000 -> MOON_PA: R3(W) @ R1(Dec) @ R3(RA) + # Then apply PA -> ME constant rotation + mats = np.zeros((len(ets), 3, 3)) + for i in range(len(ets)): + j2000_to_pa = _R3(w[i]) @ _R1(dec[i]) @ _R3(ra[i]) + mats[i] = _PA_TO_ME @ j2000_to_pa + + return mats + + +def moon_me_to_j2000(ets): + """Transpose of j2000_to_moon_me.""" + return np.swapaxes(j2000_to_moon_me(ets), -2, -1) - mat = " ".join(map("{:.7f}".format, ecef_to_enu.flatten())) - - fmt_strs = [ - "FRAME_{1} = {0}", - "FRAME_{0}_NAME = '{1}'", - "FRAME_{0}_CLASS = 4", - "FRAME_{0}_CLASS_ID = {0}", - "FRAME_{0}_CENTER = {2}", - "OBJECT_{2}_FRAME = '{1}'", - "TKFRAME_{0}_RELATIVE = '{3}'", - "TKFRAME_{0}_SPEC = 'MATRIX'", - "TKFRAME_{0}_MATRIX = {4}", - ] - frame_specs = [ - s.format(idnum, station_name, fm_center_id, relative, mat) for s in fmt_strs +# --------------------- +# Ephemeris (jplephem) +# --------------------- + +_SPK = None +_BIG_KERNELS = ["spk/planets/de430.bsp"] +_NAIF_KERNEL_URL = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/" + + +def download_big_kernels(show_progress=True): + """ + Pre-fetch large (~150 MB) SPICE kernel files into astropy's data cache. + + Calling this is optional — the kernels are auto-downloaded on first use. + It is useful for offline preparation or to control when the download + happens (e.g. in CI, or before going offline). + """ + return [ + download_file(_NAIF_KERNEL_URL + kern, cache=True, show_progress=show_progress) + for kern in _BIG_KERNELS ] - frame_dict = {} - for spec in frame_specs: - k, v = map(str.strip, spec.split("=")) - frame_dict[k] = v +def _ensure_spk(): + global _SPK + if _SPK is None: + from jplephem.spk import SPK - latlon = ["{:.4f}".format(coord) for coord in [latitude, longitude]] + spk_path = download_file( + _NAIF_KERNEL_URL + "spk/planets/de430.bsp", + cache=True, + show_progress=True, + ) + _SPK = SPK.open(spk_path) + return _SPK - return station_name, idnum, frame_dict, latlon +def _et_to_jd(ets): + """Convert ET (TDB seconds past J2000) to Julian date.""" + return np.asarray(ets) / 86400.0 + 2451545.0 -def remove_topo(station_num): - """Remove a lunar station, by number, from variable pool.""" - idnum = 1301000 + station_num - fm_center_id = idnum - 1000000 - station_name = f"LUNAR-TOPO-{station_num}" +def _pos_ssb_j2000(body_id, ets): + """ + Get position of a solar system body relative to SSB in J2000. - fmt_vars = [ - "FRAME_{1}", - "FRAME_{0}_NAME", - "FRAME_{0}_CLASS", - "FRAME_{0}_CLASS_ID", - "FRAME_{0}_CENTER", - "OBJECT_{2}_FRAME", - "TKFRAME_{0}_RELATIVE", - "TKFRAME_{0}_SPEC", - "TKFRAME_{0}_MATRIX", - ] + Parameters + ---------- + body_id : int + 0=SSB, 301=Moon, 399=Earth + ets : ndarray + ET values - frame_vars = [s.format(idnum, station_name, fm_center_id) for s in fmt_vars] + Returns + ------- + pos : ndarray, shape (N, 3) + Position in km. + """ + if body_id == 0: + return np.zeros((len(ets), 3)) + + spk = _ensure_spk() + jd = _et_to_jd(ets) + + if body_id == 301: + emb = np.asarray(spk[0, 3].compute(jd)).T # spk keys are tuples (center, target) + moon_emb = np.asarray(spk[3, 301].compute(jd)).T + return emb + moon_emb + elif body_id == 399: + emb = np.asarray(spk[0, 3].compute(jd)).T + earth_emb = np.asarray(spk[3, 399].compute(jd)).T + return emb + earth_emb + else: + raise ValueError(f"Unsupported body ID {body_id}") - # Handle a bug in spiceypy for older versions of numpy - if np.str_ is None: - return - for var in frame_vars: - spice.dvpool(var) - # Ideally, one would also remove the ephemeris data for this station from - # the lunar_points.bsp file. This doesn't seem to be possible. However, - # the ephemeris _should_ be overwritten if the station_id is reused. +def body_position(target_id, ets, frame, observer_id): + """ + Get position of target relative to observer in given frame. + + Parameters + ---------- + target_id : int + NAIF body ID of target (0=SSB, 301=Moon, 399=Earth) + ets : array_like + ET values (TDB seconds past J2000) + frame : str + "J2000" or "MOON_ME" + observer_id : int + NAIF body ID of observer + + Returns + ------- + pos : ndarray, shape (N, 3) + Position in km. + """ + ets = np.atleast_1d(np.asarray(ets, dtype=float)) + + pos_j2000 = _pos_ssb_j2000(target_id, ets) - _pos_ssb_j2000(observer_id, ets) + + if frame == "J2000": + return pos_j2000 + elif frame == "MOON_ME": + mats = j2000_to_moon_me(ets) + return np.einsum("nij,nj->ni", mats, pos_j2000) + else: + raise ValueError(f"Unsupported frame {frame}") + + +def station_pos_ssb_j2000(pos_me_km, ets): + """ + Get position of a lunar surface station relative to SSB in J2000. + + Parameters + ---------- + pos_me_km : ndarray, shape (3,) + Station position in MOON_ME frame, in km. + ets : ndarray + ET values. + + Returns + ------- + pos : ndarray, shape (N, 3) + Position in km. + """ + ets = np.atleast_1d(np.asarray(ets, dtype=float)) + moon_ssb = _pos_ssb_j2000(301, ets) + me_to_j2000 = moon_me_to_j2000(ets) + station_j2000 = np.einsum("nij,j->ni", me_to_j2000, pos_me_km) + return moon_ssb + station_j2000 + + +def topo_rotation_matrix(lat_deg, lon_deg): + """ + Compute the MOON_ME -> topocentric (E/N/U) rotation matrix for a surface location. + + Parameters + ---------- + lat_deg, lon_deg : float + Selenodetic latitude and longitude in degrees. + + Returns + ------- + matrix : ndarray, shape (3, 3) + """ + ecef_to_enu = np.matmul( + rotation_matrix(-lon_deg, "z", unit="deg"), + rotation_matrix(lat_deg, "y", unit="deg"), + ).T + ecef_to_enu = ecef_to_enu[[2, 1, 0]] + return np.asarray(ecef_to_enu) def earth_pos_mcmf(obstimes): """ Get the position of the Earth in the MCMF frame. - Using SPICE. - Used for tests. """ - ets = (obstimes - Time("J2000")).sec - earthpos = np.stack( - [spice.spkpos("399", et, "MOON_ME", "None", "301")[0] for et in ets] - ) + ets = (obstimes.tdb - Time("J2000")).sec + earthpos = body_position(399, ets, "MOON_ME", 301) earthpos = unit.Quantity(earthpos.T, "km") return MCMF(*earthpos, obstime=obstimes) - - -KERNEL_PATHS = furnish_kernels() diff --git a/lunarsky/tests/conftest.py b/lunarsky/tests/conftest.py index e79153d..89e3608 100644 --- a/lunarsky/tests/conftest.py +++ b/lunarsky/tests/conftest.py @@ -3,11 +3,7 @@ import pytest import warnings -from lunarsky.spice_utils import remove_topo -from lunarsky.moon import MoonLocation - -# Ignore deprecation warnings from spiceypy @pytest.fixture(autouse=True) def ignore_representation_deprecation(): warnings.filterwarnings( @@ -30,17 +26,3 @@ def grcat(): stars = SkyCoord(ra=ras, dec=decs, unit="deg", frame="icrs") return stars - - -@pytest.fixture -def cleanup_moonlocs(): - # Reset the existing set of station IDs and corresponding kernel vars. - yield - - for sid in MoonLocation._inuse_stat_ids: - remove_topo(sid) - - MoonLocation._inuse_stat_ids = [] - MoonLocation._avail_stat_ids = None - MoonLocation._existing_locs = [] - MoonLocation._ref_count = [] diff --git a/lunarsky/tests/data/spice_fixtures.npz b/lunarsky/tests/data/spice_fixtures.npz new file mode 100644 index 0000000..1d3e4e8 Binary files /dev/null and b/lunarsky/tests/data/spice_fixtures.npz differ diff --git a/lunarsky/tests/test_moon.py b/lunarsky/tests/test_moon.py index 08b4cbd..fd605da 100644 --- a/lunarsky/tests/test_moon.py +++ b/lunarsky/tests/test_moon.py @@ -1,5 +1,4 @@ import pytest -import copy import numpy as np import astropy.units as unit from astropy.coordinates import ( @@ -8,7 +7,6 @@ ) from astropy.tests.helper import quantity_allclose from lunarsky import MoonLocation, MoonLocationAttribute, MCMF -from lunarsky.moon import SELENOIDS class TestsWithObject: @@ -107,151 +105,3 @@ def test_moonlocation_attribute(): attr2, boo = mlattr.convert_input(moonloc.mcmf) assert np.all(attr2 == moonloc) - - -def test_moonlocation_copy(): - # Check that station_ids are copied properly - loc0 = MoonLocation.from_selenodetic(lat=["-15d", "25d"], lon=["97d", "0d"]) - loc0._set_station_id() - lcop = loc0.copy() - assert lcop.station_ids == loc0.station_ids - lcop2 = copy.copy(loc0) - assert lcop2.station_ids == loc0.station_ids - - -def test_moonlocation_delete(): - # Check that making multiple instances of the same location raises the _ref_count - # appropriately, and deleting them removes them correctly. - before = copy.copy(MoonLocation._inuse_stat_ids) - if MoonLocation._avail_stat_ids is not None: - exp = before + MoonLocation._avail_stat_ids[-2:][::-1] - else: - exp = before + [1, 2] - - locs = [] - for ii in range(5): - locs.append(MoonLocation.from_selenodetic(lat=["-15d", "25d"], lon=["97d", "0d"])) - locs[-1]._set_station_id() - - check0 = copy.copy(MoonLocation._inuse_stat_ids) - for ii in range(5, 0, -1): - assert MoonLocation._ref_count[-1] == ii - assert MoonLocation._inuse_stat_ids[-1] == check0[-1] - locs.pop(ii - 1) - - check1 = copy.copy(MoonLocation._inuse_stat_ids) - - assert exp == check0 - assert check1 == before - - -@pytest.mark.parametrize("ell", SELENOIDS) -def test_station_ids(ell): - # Check that when a MoonLocations are made, the appropriate station_ids are assigned. - - # Get whatever are already recorded in the class. - orig_statids = copy.copy(MoonLocation._inuse_stat_ids) - - # Random positions with some repeats, in five groups - lonlatheights = [ - [ - (294.67, -67.68, 15.23), - (78.31, 51.60, 16.59), - (335.27, 57.45, 27.32), - ], - [ - (335.27, 57.45, 27.32), - (133.29, -80.63, 25.92), - (323.28, -2.46, 4.74), - (45.94, -25.17, 20.22), - (326.40, 29.21, 9.67), - (197.94, 3.92, 16.34), - ], - [(197.94, 3.92, 16.34)], - [ - (197.94, 3.92, 16.34), - (242.63, -20.86, 4.90), - (70.30, -40.22, 29.91), - (187.70, -85.44, 9.35), - (129.55, 14.11, 0.24), - (140.46, 4.19, 28.03), - (232.22, 56.89, 15.22), - (51.31, -30.92, 11.53), - (55.59, -61.90, 3.98), - ], - [ - (307.76, 32.31, 29.88), - (248.07, 2.89, 22.41), - (349.94, -77.10, 14.27), - (315.41, -72.71, 13.26), - (173.39, 14.43, 3.37), - (183.19, -83.01, 16.13), - (258.53, -56.66, 29.75), - (115.36, -44.42, 11.34), - (270.61, 57.85, 26.16), - (320.07, -75.06, 21.15), - (216.74, 54.75, 1.92), - ], - ] - - all_pos = np.array(sum(lonlatheights, [])) - n_unique = np.unique(all_pos, axis=0).shape[0] - - locs = [] - locstrs = [] - - for gp in lonlatheights: - lons, lats, heights = np.array(gp).T - locs.append( - MoonLocation.from_selenodetic( - lat=lats, lon=lons, height=heights, ellipsoid=ell - ) - ) - locs[-1]._set_station_id() - - # Check that only unique positions got added - added = len(MoonLocation._inuse_stat_ids) - len(orig_statids) - assert added == n_unique - - statids = [loc.station_ids for loc in locs] - - locstrs = [] - for inst in locs: - llhs = [] - if inst.isscalar: - llh_arr = [ - ( - inst.lon.deg.item(), - inst.lat.deg.item(), - inst.height.to_value("km").item(), - ) - ] - else: - llh_arr = zip(inst.lon.deg, inst.lat.deg, inst.height.to_value("km")) - for llh in llh_arr: - llhs.append("_".join(["{:.4f}".format(ll) for ll in llh] + [ell])) - locstrs.append(llhs) - - # Check that saved location strings correspond with station IDs in each instance - for gi, gp in enumerate(statids): - for sti, sid in enumerate(gp): - ind = MoonLocation._inuse_stat_ids.index(sid) - assert locstrs[gi][sti] == MoonLocation._existing_locs[ind] - - -def test_statid_management(cleanup_moonlocs): - lats = np.random.uniform(-90, 90, 1000) - lons = np.random.uniform(0, 360, 1000) - - for ii in range(2): - loc = MoonLocation.from_selenodetic( - lon=lons[ii * 500 : (ii + 1) * 500], - lat=lats[ii * 500 : (ii + 1) * 500], - ) - loc._set_station_id() - assert len(loc.station_ids) == 500 - del loc - - bigloc = MoonLocation.from_selenodetic(lon=lons, lat=lats) - with pytest.raises(ValueError, match="Too many unique MoonLocation"): - bigloc._set_station_id() diff --git a/lunarsky/tests/test_spice.py b/lunarsky/tests/test_spice.py index fcdfc51..2c90c5d 100644 --- a/lunarsky/tests/test_spice.py +++ b/lunarsky/tests/test_spice.py @@ -1,120 +1,88 @@ -from astropy.coordinates.baseframe import frame_transform_graph -from astropy.coordinates.transformations import FunctionTransformWithFiniteDifference -from astropy.coordinates import AltAz, ICRS, EarthLocation, Angle -from astropy.utils.data import download_files_in_parallel -import lunarsky +import numpy as np +import os +import pytest import lunarsky.spice_utils as spice_utils -from lunarsky.time import Time -import spiceypy as spice +# spice_fixtures.npz contains saved results of the earlier version of lunarsky that was based +# on the SPICE toolkit. Tests here are meant to verify that the new code is consistent with +# older results. +FIXTURE_PATH = os.path.join(os.path.dirname(__file__), "data", "spice_fixtures.npz") -def test_topo_frame_setup(): - # Check that the frame setup puts all the correct values in the kernel pool. +@pytest.fixture +def fixtures(): + return np.load(FIXTURE_PATH) - latitude, longitude = 30, 25 - name, idnum, frame_dict, latlon = spice_utils.topo_frame_def(latitude, longitude) - frame_strs = ["{}={}".format(k, v) for (k, v) in frame_dict.items()] - spice.lmpool(frame_strs) - for k, v in frame_dict.items(): - N, typecode = spice.dtpool(k) - if typecode == "N": - res = spice.gdpool(k, 0, 100) - if len(res) == 1: - res = [int(res[0])] - if N > 1: - v = [float(it) for it in v.split(" ")] - res = res.tolist() - else: - res = res[0] - v = int(v) - assert v == res - else: - res = spice.gcpool(k, 0, 100)[0] - v = v.replace("'", "") - assert v == res +def test_kernels_available(): + """ + Ensure the large SPK kernel is downloadable and cached. + On a fresh environment this hits the JPL server; on subsequent runs + (and in CI with a populated astropy cache) it is a no-op. + """ + paths = spice_utils.download_big_kernels(show_progress=False) + for p in paths: + assert os.path.exists(p) -def test_kernel_paths(): - # Check that the correct kernel files are downloaded - # Need to unhash the file names - assert len(lunarsky.spice_utils.KERNEL_PATHS) == 5 +def test_j2000_to_moon_me(fixtures): + ets = fixtures["ets"] + expected = fixtures["j2000_to_me"] + computed = spice_utils.j2000_to_moon_me(ets) + np.testing.assert_allclose(computed, expected, atol=1e-10) -def test_spice_earth(grcat): - # Replace the ICRS->AltAz transform in astropy with one using SPICE. - # Confirm that star positions are the same as with the original transform - # to within the error due to relativistic aberration (~ 21 arcsec) +def test_body_position_moon_earth(fixtures): + ets = fixtures["ets"] + expected = fixtures["moon_earth_j2000"][:, :3] + computed = spice_utils.body_position(301, ets, "J2000", 399) + np.testing.assert_allclose(computed, expected, atol=1e-6) - # DSS-15, 399015 - lat, lon, height = 35.42185, -116.8871951, 973.211 - loc = EarthLocation.from_geodetic(lon, lat, height) - t0 = Time.now() - _J2000 = Time("J2000") +def test_body_position_moon_ssb(fixtures): + ets = fixtures["ets"] + expected = fixtures["moon_ssb_j2000"][:, :3] + computed = spice_utils.body_position(301, ets, "J2000", 0) + np.testing.assert_allclose(computed, expected, atol=1e-6) - altaz = grcat.transform_to(AltAz(location=loc, obstime=t0)) - # Make the Earth topo frame in spice. - framename, idnum, frame_dict, latlon = lunarsky.spice_utils.topo_frame_def( - lat, lon, moon=False - ) - - # One more kernel is needed for the ITRF93 frame. - kname = "pck/earth_latest_high_prec.bpc" - _naif_kernel_url = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels" - kurl = [_naif_kernel_url + "/" + kname] - kernpath = download_files_in_parallel( - kurl, - cache=True, - show_progress=False, - ) - spice.furnsh(kernpath) - - frame_strs = ["{}={}".format(k, v) for (k, v) in frame_dict.items()] - spice.lmpool(frame_strs) - et = (t0 - _J2000).sec +def test_body_position_ssb_moon_me(fixtures): + ets = fixtures["ets"] + expected = fixtures["ssb_moon_me"][:, :3] + computed = spice_utils.body_position(0, ets, "MOON_ME", 301) + np.testing.assert_allclose(computed, expected, atol=1e-4) - @frame_transform_graph.transform(FunctionTransformWithFiniteDifference, ICRS, AltAz) - def icrs_to_ecef(icrs_coo, ecef_frame): - mat = spice.pxform("J2000", framename, et) - newrepr = icrs_coo.cartesian.transform(mat) - return ecef_frame.realize_frame(newrepr) +def test_earth_pos_mcmf(fixtures): + ets = fixtures["ets"] + expected = fixtures["earth_me"] + computed = spice_utils.body_position(399, ets, "MOON_ME", 301) + np.testing.assert_allclose(computed, expected, atol=1e-4) - tr = frame_transform_graph.get_transform(ICRS, AltAz).transforms - assert tr[0].func.__name__ == "icrs_to_ecef" # Confirm new transform added correctly - - altaz2 = grcat.transform_to(AltAz(location=loc, obstime=t0)) - - # Having done the transform, remove the spice transform from the graph - frame_transform_graph.remove_transform(ICRS, AltAz, None) - - # Compare astropy topocentric coordinate transform results to - assert all(altaz.separation(altaz2) < Angle("25arcsec")) - # NOTE Large deviations in azimuth high altitudes. +def test_topo_rotation_matrix(fixtures): + expected = fixtures["topo_to_me"] + topo_matrix = spice_utils.topo_rotation_matrix( + float(fixtures["topo_station_lat_deg"]), + float(fixtures["topo_station_lon_deg"]), + ) + # TOPO -> MOON_ME should be topo_matrix.T (constant) + for i in range(len(expected)): + np.testing.assert_allclose(topo_matrix.T, expected[i], atol=1e-7) -def test_topo_kernel_setup(): - # Tests a single function in topo.py - # Need to clear the kernel pool first + # Verify full precision matrix is orthogonal + np.testing.assert_allclose(topo_matrix @ topo_matrix.T, np.eye(3), atol=1e-14) - spice.clpool() - # Confirm no variables are loaded. - assert not lunarsky.spice_utils.check_is_loaded("*") - for filepath in lunarsky.spice_utils.KERNEL_PATHS: - spice.furnsh(filepath) - loc = lunarsky.MoonLocation.from_selenodetic(lat="30d", lon="20d") - station_name, idnum, frame_specs, latlon = lunarsky.spice_utils.topo_frame_def( - loc.lat.deg, loc.lon.deg, moon=True +def test_topo_to_j2000(fixtures): + ets = fixtures["ets"] + expected_to_j2000 = fixtures["topo_to_j2000"] + topo_matrix = spice_utils.topo_rotation_matrix( + float(fixtures["topo_station_lat_deg"]), + float(fixtures["topo_station_lon_deg"]), ) - statnum = idnum - 1301000 - lunarsky.topo._spice_setup(loc, [statnum]) - try: - assert lunarsky.spice_utils.check_is_loaded("*{}*".format(idnum)) - finally: - lunarsky.spice_utils.remove_topo(statnum) + me_to_j2000 = spice_utils.moon_me_to_j2000(ets) + computed_to_j2000 = np.einsum("nij,jk->nik", me_to_j2000, topo_matrix.T) + np.testing.assert_allclose(computed_to_j2000, expected_to_j2000, atol=1e-10) diff --git a/lunarsky/time.py b/lunarsky/time.py index 1372a77..23b46a4 100644 --- a/lunarsky/time.py +++ b/lunarsky/time.py @@ -6,7 +6,7 @@ from astropy.coordinates import EarthLocation, Longitude from .moon import MoonLocation, COPY_IF_NEEDED -import spiceypy as spice +from .spice_utils import moon_me_to_j2000 __all__ = ["Time", "TimeDelta"] @@ -70,8 +70,8 @@ def sidereal_time(self, kind, longitude=None, model=None): # to the selenodetic coordinate system. "self.location" must # be defined in order to get here. - et = np.atleast_1d((self - Time("J2000")).sec) - mats = np.array([spice.pxform("MOON_ME", "J2000", t) for t in et]) + et = np.atleast_1d((self.tdb - Time("J2000")).sec) + mats = moon_me_to_j2000(et) # Zenith vector zvec = self.location.mcmf.cartesian.xyz diff --git a/lunarsky/topo.py b/lunarsky/topo.py index 3146481..2fa473f 100644 --- a/lunarsky/topo.py +++ b/lunarsky/topo.py @@ -22,9 +22,12 @@ from .mcmf import MCMF from .moon import MoonLocationAttribute -from .spice_utils import check_is_loaded, topo_frame_def, lunar_surface_ephem - -import spiceypy as spice +from .spice_utils import ( + j2000_to_moon_me, + moon_me_to_j2000, + station_pos_ssb_j2000, + topo_rotation_matrix, +) _90DEG = 90 * un.deg __all__ = ["LunarTopo"] @@ -88,22 +91,24 @@ def zen(self): # ----------------- -def _spice_setup(locations, station_ids): - for li, loc in enumerate(np.atleast_1d(locations)): - sid = int(station_ids[li]) - frameloaded = check_is_loaded(f"FRAME_LUNAR-TOPO-{sid}") - if not frameloaded: - lunar_surface_ephem( - loc.x.to_value("km"), - loc.y.to_value("km"), - loc.z.to_value("km"), - station_num=sid, - ) # Furnishes SPK for lunar surface point - station_name, idnum, frame_specs, latlon = topo_frame_def( - loc.lat.deg, loc.lon.deg, moon=True, station_num=sid - ) - frame_strs = ["{}={}".format(k, v) for (k, v) in frame_specs.items()] - spice.lmpool(frame_strs) +def _get_topo_data(location): + """ + Compute topo rotation matrices and MOON_ME positions for each location element. + + Returns + ------- + topo_mats : ndarray, shape (M, 3, 3) + ME-to-topo rotation matrices. + pos_mes : ndarray, shape (M, 3) + Station positions in MOON_ME, in km. + """ + locations = np.atleast_1d(location) + topo_mats = [] + pos_mes = [] + for loc in locations: + topo_mats.append(topo_rotation_matrix(loc.lat.deg, loc.lon.deg)) + pos_mes.append([loc.x.to_value("km"), loc.y.to_value("km"), loc.z.to_value("km")]) + return np.array(topo_mats), np.array(pos_mes) def make_transform(coo, toframe): @@ -123,29 +128,28 @@ def make_transform(coo, toframe): if location is None: raise ValueError("location must be defined for LunarTopo transformations") - # Initialize station_ids if not defined. - if location.station_ids == []: - location.__class__._set_site_id(location) - # Make arrays - ets = (obstime - _J2000).sec - stat_ids = np.asarray(location.station_ids) + ets = (obstime.tdb - _J2000).sec + topo_mats, pos_mes = _get_topo_data(location) + loc_indices = np.arange(len(topo_mats)) shape_out = np.broadcast_shapes(coo.shape, ets.shape, location.shape) - # Set up SPICE ephemerides and frame details - _spice_setup(location, stat_ids) - - ets_ids = np.atleast_2d(np.stack(np.broadcast_arrays(ets, stat_ids)).T) + ets_locs = np.atleast_2d(np.stack(np.broadcast_arrays(ets, loc_indices)).T) coo_cart = coo.cartesian # Make rotation matrices - mats = np.asarray( - [ - spice.pxform(f"LUNAR-TOPO-{int(stat_id)}", frame_spice_name, et) - for (et, stat_id) in ets_ids - ] - ) + # TOPO->ME is topo_matrix.T (constant); TOPO->J2000 = R(ME->J2000) @ topo_matrix.T + mats_list = [] + for et, loc_idx in ets_locs: + loc_idx = int(loc_idx) + topo_mat = topo_mats[loc_idx] + if frame_spice_name == "MOON_ME": + mats_list.append(topo_mat.T) + else: + me_to_j2000 = moon_me_to_j2000(np.array([et]))[0] + mats_list.append(me_to_j2000 @ topo_mat.T) + mats = np.asarray(mats_list) if totopo: mats = np.linalg.inv(mats) @@ -157,29 +161,33 @@ def make_transform(coo, toframe): # If not unitspherical, shift by origin vector before rotating. if not is_unitspherical: - # Make origin vector(s) in coo's frame. - if totopo: - origin_id = lambda n: int(frame_id) # MCMF or ICRS frame origin - target_id = lambda n: int(n) + 301000 # Station ID - frame_name = lambda n: frame_spice_name - - else: - origin_id = lambda n: int(n) + 301000 - target_id = lambda n: int(frame_id) - frame_name = lambda n: f"LUNAR-TOPO-{int(n)}" - - orig_posvel = ( - np.asarray( - [ - spice.spkgeo( - target_id(stat_id), et, frame_name(stat_id), origin_id(stat_id) - )[0] - for (et, stat_id) in ets_ids - ] - ) - * un.km - ) - coo_cart -= CartesianRepresentation((orig_posvel.T)[:3]) + orig_pos_list = [] + for et, loc_idx in ets_locs: + loc_idx = int(loc_idx) + pos_me = pos_mes[loc_idx] + topo_mat = topo_mats[loc_idx] + et_arr = np.array([et]) + + if totopo: + # Origin = station position relative to frame origin, in coo's frame + if frame_id == 0: # SSB / ICRS + orig_pos_list.append(station_pos_ssb_j2000(pos_me, et_arr)[0]) + else: # Moon center / MCMF + orig_pos_list.append(pos_me) + else: + # Origin = frame origin relative to station, in coo's frame (topo) + if frame_id == 0: # SSB + station_ssb = station_pos_ssb_j2000(pos_me, et_arr)[0] + ssb_station_j2000 = -station_ssb + # Rotate to topo + j2000_me = j2000_to_moon_me(et_arr)[0] + ssb_station_me = j2000_me @ ssb_station_j2000 + orig_pos_list.append(topo_mat @ ssb_station_me) + else: # Moon center + orig_pos_list.append(-(topo_mat @ pos_me)) + + orig_pos = np.asarray(orig_pos_list) * un.km + coo_cart -= CartesianRepresentation(orig_pos.T) newrepr = coo_cart.transform(mats).reshape(shape_out) diff --git a/pyproject.toml b/pyproject.toml index eaf220e..57c857c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,7 +24,6 @@ keywords = [ dependencies = [ "numpy>=1.15", "astropy>=6.0.0", - "spiceypy", "jplephem", "fastkml[lxml]", ]