From 8b8a0dbe3acac20537559b286993b3caf083779c Mon Sep 17 00:00:00 2001 From: Dar Dahlen Date: Fri, 5 Jul 2024 09:28:40 -0700 Subject: [PATCH] Remove SpiceKernels class, changing methods to functions (#69) Remove SpiceKernels class, changing meths to funcs See #68 for more details --- CHANGELOG.md | 2 + docs/tutorials/kona.rst | 2 +- docs/tutorials/propagation_plots.py | 8 +- src/examples/galactic_center.py | 2 +- src/examples/on_sky.py | 2 +- src/examples/plot_close_approach.py | 2 +- src/examples/plot_close_encounter.py | 7 +- src/examples/plot_elevation.py | 6 +- src/examples/plot_light_curve.py | 4 +- src/examples/plot_mpc_state.py | 4 +- src/examples/plot_observability.py | 8 +- src/examples/plot_phases.py | 2 +- src/examples/plot_trojans.py | 2 +- src/examples/plot_uncertainty.py | 4 +- src/neospy/__init__.py | 4 +- src/neospy/fov.py | 4 +- src/neospy/mpc.py | 4 +- src/neospy/propagation.py | 11 +- src/neospy/rust/time.rs | 24 +- src/neospy/spice.py | 637 +++++++++++++-------------- src/neospy/wise.py | 4 +- src/neospy/ztf.py | 4 +- src/tests/test_propagation.py | 16 +- src/tests/test_spice.py | 34 +- 24 files changed, 403 insertions(+), 394 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9f1706a..97639a9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,6 +20,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed +- Removed `SpiceKernel` as a class, lowering all its methods to the submodule level, + see #68 for more discussion. - Removed `Time` object which was a wrapper over astropy.Time, instead making a custom implementation of time which is ~3-400x faster than previous. - Improved orbital element conversion, leading to a 2-3x speedup in two body orbit diff --git a/docs/tutorials/kona.rst b/docs/tutorials/kona.rst index 96df0de..a4ea63d 100644 --- a/docs/tutorials/kona.rst +++ b/docs/tutorials/kona.rst @@ -35,7 +35,7 @@ calculated from the corners of the frame. frame = fits.open("data/01772a127-w3-int-1b.fits.gz")[0] # Here we compute a State of the observer, this could also be constructed - # using spice kernels using neospy.SpiceKernels.state + # using spice kernels using neospy.spice.state sc_pos = neospy.Vector([frame.header['SUN2SCX'], frame.header['SUN2SCY'], frame.header['SUN2SCZ']], diff --git a/docs/tutorials/propagation_plots.py b/docs/tutorials/propagation_plots.py index a28a1fc..19085c1 100644 --- a/docs/tutorials/propagation_plots.py +++ b/docs/tutorials/propagation_plots.py @@ -2,12 +2,12 @@ import matplotlib.pyplot as plt import neospy -neospy.SpiceKernels.kernel_reload(["./data/20000042.bsp"]) +neospy.spice.kernel_reload(["./data/20000042.bsp"]) jd_start = neospy.Time.from_ymd(1920, 1, 1).jd jd_end = neospy.Time.from_ymd(2020, 1, 1).jd -state = neospy.SpiceKernels.state("42", jd_end) +state = neospy.spice.state("42", jd_end) jds = np.logspace(np.log10(jd_end), np.log10(jd_end - 10), 1000) jds = np.concatenate( @@ -21,9 +21,9 @@ error_2body = [] n_body_no_asteroids = [] n_body_ast = [] -state = neospy.SpiceKernels.state("42", jd_end) +state = neospy.spice.state("42", jd_end) for jd in jds: - jpl_pos = neospy.SpiceKernels.state("42", jd).pos + jpl_pos = neospy.spice.state("42", jd).pos line = state.pos + state.vel * (jd - state.jd) error_line.append((jpl_pos - line).r * neospy.constants.AU_KM) diff --git a/src/examples/galactic_center.py b/src/examples/galactic_center.py index f9f9430..dbfa929 100644 --- a/src/examples/galactic_center.py +++ b/src/examples/galactic_center.py @@ -36,7 +36,7 @@ while states[0].jd < jd_end: states = neospy.propagate_two_body(states, states[0].jd + 1) - earth = neospy.SpiceKernels.state("Earth", states[0].jd) + earth = neospy.spice.state("Earth", states[0].jd) earth_to_obj = [(s.pos - earth.pos) for s in states] dist_to_galactic_center.append( diff --git a/src/examples/on_sky.py b/src/examples/on_sky.py index aaac74d..a891092 100644 --- a/src/examples/on_sky.py +++ b/src/examples/on_sky.py @@ -56,7 +56,7 @@ print(obs_info) # Load the position of the observer in the solar system: -observer = neospy.SpiceKernels.mpc_code_to_ecliptic("Palomar Mountain", jd) +observer = neospy.spice.mpc_code_to_ecliptic("Palomar Mountain", jd) # %% Propagation and Orbit calculations diff --git a/src/examples/plot_close_approach.py b/src/examples/plot_close_approach.py index 90564a2..bdfbfd4 100644 --- a/src/examples/plot_close_approach.py +++ b/src/examples/plot_close_approach.py @@ -40,7 +40,7 @@ rel_earth_state = cur_state.change_center(399) pos.append([rel_earth_state.pos.x, rel_earth_state.pos.y]) # gets earths position and record the distances. - moon = neospy.SpiceKernels.state("Moon", cur_state.jd, center=399).pos + moon = neospy.spice.state("Moon", cur_state.jd, center=399).pos dist_to_moon.append((moon - rel_earth_state.pos).r * neospy.constants.AU_KM) dist_to_earth.append(rel_earth_state.pos.r * neospy.constants.AU_KM) moon_pos.append([moon.x, moon.y]) diff --git a/src/examples/plot_close_encounter.py b/src/examples/plot_close_encounter.py index 84113c0..43fca11 100644 --- a/src/examples/plot_close_encounter.py +++ b/src/examples/plot_close_encounter.py @@ -29,7 +29,6 @@ # Sample the time from start to end in 1 day steps jds = np.arange(jd_start, jd_end) -eph = neospy.SpiceKernels() # Propagate the object to all of the state obj_pos = [] @@ -39,7 +38,7 @@ for jd in jds: last_state = neospy.propagate_two_body([last_state], jd)[0] obj_pos.append(last_state.pos) - earth_pos.append(eph.state("Earth", jd).pos) + earth_pos.append(neospy.spice.state("Earth", jd).pos) earth_r.append(neospy.Vector(obj_pos[-1] - earth_pos[-1]).r) # Find the time where the closest encounter with the earth @@ -56,10 +55,10 @@ ax.plot(pos[0], pos[1], pos[2], alpha=0.5, color="black") ax.scatter(pos_center.x, pos_center.y, pos_center.z, s=5, color="black") for i, planet in enumerate(["Mercury", "Venus", "Earth", "Mars", "Jupiter"]): - states = [eph.state(planet, jd).pos for jd in jds] + states = [neospy.spice.state(planet, jd).pos for jd in jds] pos = np.array(states).T ax.plot(pos[0], pos[1], pos[2], color=f"C{i}", alpha=0.5) - pos = eph.state(planet, jd_center).pos + pos = neospy.spice.state(planet, jd_center).pos ax.scatter(pos.x, pos.y, pos.z, color=f"C{i}", s=5) ax.scatter(0, 0, 0, color="red") ax.set_xticks([-zoom, 0, zoom]) diff --git a/src/examples/plot_elevation.py b/src/examples/plot_elevation.py index 3e827e6..f98efa2 100644 --- a/src/examples/plot_elevation.py +++ b/src/examples/plot_elevation.py @@ -59,10 +59,8 @@ # for each time step, calculate the elevation and sun elevation jd_step = jd_start + subrange approx_state = neospy.propagate_two_body(states, jd_step) - sun2obs = neospy.SpiceKernels.mpc_code_to_ecliptic(site, jd_step).pos - earth2obs = neospy.SpiceKernels.mpc_code_to_ecliptic( - site, jd_step, center="399" - ).pos + sun2obs = neospy.spice.mpc_code_to_ecliptic(site, jd_step).pos + earth2obs = neospy.spice.mpc_code_to_ecliptic(site, jd_step, center="399").pos sun_elevation.append(90 - earth2obs.angle_between(-sun2obs)) cur_elev = [] diff --git a/src/examples/plot_light_curve.py b/src/examples/plot_light_curve.py index 528061c..ce9436e 100644 --- a/src/examples/plot_light_curve.py +++ b/src/examples/plot_light_curve.py @@ -11,7 +11,7 @@ # Using ceres as a source for a state vector -state = neospy.SpiceKernels.state("ceres", 2460000.5) +state = neospy.spice.state("ceres", 2460000.5) # Various input values albedo = 0.1 @@ -39,7 +39,7 @@ for dt in dts: # Find the observer and object positions some time in the future. - earth_pos = neospy.SpiceKernels.state("Earth", jd + dt).pos + earth_pos = neospy.spice.state("Earth", jd + dt).pos final_pos = neospy.propagate_two_body([state], jd + dt, earth_pos)[0].pos obj2obs = final_pos - earth_pos diff --git a/src/examples/plot_mpc_state.py b/src/examples/plot_mpc_state.py index d39d460..4396713 100644 --- a/src/examples/plot_mpc_state.py +++ b/src/examples/plot_mpc_state.py @@ -39,10 +39,10 @@ # Plot the planets for 1 orbit each for i, planet in enumerate(["Mercury", "Venus", "Earth", "Mars", "Jupiter"]): jd = states[0].jd - plan = neospy.SpiceKernels().state(planet, jd) + plan = neospy.spice.state(planet, jd) ax.scatter(plan.pos.x, plan.pos.y, color=f"C{i}", s=10) jds = np.linspace(jd - plan.elements.orbital_period, jd, 100) - pos = np.array([neospy.SpiceKernels().state(planet, jd).pos for jd in jds]).T + pos = np.array([neospy.spice.state(planet, jd).pos for jd in jds]).T ax.plot(pos[0], pos[1], pos[2], color="black", alpha=0.2) for state in states: diff --git a/src/examples/plot_observability.py b/src/examples/plot_observability.py index 2daccd3..2b0c335 100644 --- a/src/examples/plot_observability.py +++ b/src/examples/plot_observability.py @@ -56,8 +56,8 @@ # For each time, compute the geometry, then compute the mag as well as all of the # various angle and distance values. state = neospy.propagate_n_body([state], t)[0] - sun2obs = neospy.SpiceKernels.mpc_code_to_ecliptic(site, t).pos - earth2obs = neospy.SpiceKernels.mpc_code_to_ecliptic(site, t, center="399").pos + sun2obs = neospy.spice.mpc_code_to_ecliptic(site, t).pos + earth2obs = neospy.spice.mpc_code_to_ecliptic(site, t, center="399").pos sun2obj = state.pos obs2obj = -sun2obs + sun2obj @@ -130,8 +130,8 @@ nights = [] for subrange in substeps: approx_state = neospy.propagate_two_body([state], t + subrange)[0] - sun2obs = neospy.SpiceKernels.mpc_code_to_ecliptic(site, t + subrange).pos - earth2obs = neospy.SpiceKernels.mpc_code_to_ecliptic( + sun2obs = neospy.spice.mpc_code_to_ecliptic(site, t + subrange).pos + earth2obs = neospy.spice.mpc_code_to_ecliptic( site, t + subrange, center="399" ).pos diff --git a/src/examples/plot_phases.py b/src/examples/plot_phases.py index f55d6c6..c52b856 100644 --- a/src/examples/plot_phases.py +++ b/src/examples/plot_phases.py @@ -27,7 +27,7 @@ # Move the entire population of asteroids to that time using 2-body # mechanics, this can be directly substituted with propagate_n_body if you # want more precision. -sun2earth = neospy.SpiceKernels.state("Earth", jd).pos +sun2earth = neospy.spice.state("Earth", jd).pos states = neospy.propagate_two_body(neos, jd) # Compute the expected V-mags for these objects at this time diff --git a/src/examples/plot_trojans.py b/src/examples/plot_trojans.py index 5534df5..f4336ad 100644 --- a/src/examples/plot_trojans.py +++ b/src/examples/plot_trojans.py @@ -16,7 +16,7 @@ states = neospy.propagate_n_body(states, states[0].jd) # Where is jupiter? -jupiter = neospy.SpiceKernels.state("Jupiter", states[0].jd) +jupiter = neospy.spice.state("Jupiter", states[0].jd) # found it! # Compute the positions, and relative longitudinal distance from jupiter diff --git a/src/examples/plot_uncertainty.py b/src/examples/plot_uncertainty.py index d30ccd8..0c3f3d4 100644 --- a/src/examples/plot_uncertainty.py +++ b/src/examples/plot_uncertainty.py @@ -50,7 +50,7 @@ mags = [] for jd in jds: states = neospy.propagate_n_body(states, jd) - earth = neospy.SpiceKernels.state("earth", jd) + earth = neospy.spice.state("earth", jd) m = [ neospy.flux.hg_apparent_mag( sun2obj=x.pos, sun2obs=earth.pos, h_mag=obj.h_mag, g_param=g @@ -65,7 +65,7 @@ # position at lowest mag states = neospy.propagate_n_body(states, brightest_jd) -earth = neospy.SpiceKernels.state("earth", brightest_jd) +earth = neospy.spice.state("earth", brightest_jd) vecs = [(s.pos - earth.pos).as_equatorial for s in states] ras = np.array([v.ra for v in vecs]) decs = np.array([v.dec for v in vecs]) diff --git a/src/neospy/__init__.py b/src/neospy/__init__.py index ab7e0e7..836a261 100644 --- a/src/neospy/__init__.py +++ b/src/neospy/__init__.py @@ -18,8 +18,8 @@ ztf, fov, shape, + spice, ) -from .spice import SpiceKernels from .propagation import ( propagate_n_body, propagate_two_body, @@ -64,12 +64,12 @@ "population", "Time", "set_logging", - "SpiceKernels", "flux", "fov", "wise", "neos", "mpc", + "spice", "SimultaneousStates", "propagate_n_body", "propagate_two_body", diff --git a/src/neospy/fov.py b/src/neospy/fov.py index 5470af1..2cfaf0e 100644 --- a/src/neospy/fov.py +++ b/src/neospy/fov.py @@ -16,7 +16,7 @@ fov_spk_check as _fov_spk_check, ) from .vector import State -from .spice import SpiceKernels +from . import spice __all__ = [ @@ -49,5 +49,5 @@ def fov_spice_check(desigs: list[str], fovs) -> list[State]: fov: A list of field of views from which to subselect objects which are visible. """ - obj_ids = [SpiceKernels.name_lookup(n)[1] for n in desigs] + obj_ids = [spice.name_lookup(n)[1] for n in desigs] return _fov_spk_check(obj_ids, fovs) diff --git a/src/neospy/mpc.py b/src/neospy/mpc.py index ed2058c..9e0820b 100644 --- a/src/neospy/mpc.py +++ b/src/neospy/mpc.py @@ -920,7 +920,7 @@ def _read_first_line(line): @staticmethod def _read_second_line(line, jd): - from .spice import SpiceKernels + from . import spice if line[14] != "s": raise ValueError("No second line of spacecraft observation found.") @@ -929,7 +929,7 @@ def _read_second_line(line, jd): y = float(line[46:57].replace(" ", "")) / constants.AU_KM z = float(line[58:69].replace(" ", "")) / constants.AU_KM earth2sc = Vector([x, y, z], Frames.Equatorial).as_ecliptic - sun2earth = SpiceKernels.state("Earth", jd).pos + sun2earth = spice.state("Earth", jd).pos sun2sc = sun2earth + earth2sc return list(sun2sc) diff --git a/src/neospy/propagation.py b/src/neospy/propagation.py index 1c6f45d..574b003 100644 --- a/src/neospy/propagation.py +++ b/src/neospy/propagation.py @@ -8,9 +8,8 @@ from scipy import optimize import numpy as np -from .spice import SpiceKernels from .vector import Vector, State -from . import _core +from . import _core, spice from ._core import NonGravModel @@ -133,10 +132,10 @@ def moid(state: State, other: Optional[State] = None): state: The state describing an object. other: - The state of the object to calculate the MOID for, if this is not provided, then - Earth is fetched from :class:`~neospy.spice.SpiceKernels` and is used in - the calculation. + The state of the object to calculate the MOID for, if this is not provided, + then Earth is fetched from :mod:`~neospy.spice` and is used in the + calculation. """ if other is None: - other = SpiceKernels.state("Earth", state.jd) + other = spice.state("Earth", state.jd) return _moid_single(state, other) diff --git a/src/neospy/rust/time.rs b/src/neospy/rust/time.rs index 1985110..b4c2016 100644 --- a/src/neospy/rust/time.rs +++ b/src/neospy/rust/time.rs @@ -101,12 +101,6 @@ impl PyTime { PyTime(Time::::from_year_month_day(year, month, day, frac_day).tdb()) } - /// Time in the UTC ISO time format. - #[getter] - pub fn iso(&self) -> PyResult { - Ok(self.0.utc().to_iso()?) - } - /// Time in the current time. #[staticmethod] pub fn now() -> Self { @@ -139,6 +133,24 @@ impl PyTime { self.0.mjd() } + /// Julian Date in UTC scaled time. + #[getter] + pub fn utc_jd(&self) -> f64 { + self.0.utc().jd + } + + /// Modified Julian Date in UTC scaled time. + #[getter] + pub fn utc_mjd(&self) -> f64 { + self.0.utc().jd + } + + /// Time in the UTC ISO time format. + #[getter] + pub fn iso(&self) -> PyResult { + Ok(self.0.utc().to_iso()?) + } + #[staticmethod] pub fn j2000() -> Self { PyTime(Time::::new(2451545.0)) diff --git a/src/neospy/spice.py b/src/neospy/spice.py index 5bd4e2c..8322459 100644 --- a/src/neospy/spice.py +++ b/src/neospy/spice.py @@ -11,7 +11,20 @@ from .cache import cached_file_download, cache_path from .vector import Frames, State -__all__ = ["SpiceKernels", "SpkInfo"] +__all__ = [ + "SpkInfo", + "state", + "name_lookup", + "loaded_objects", + "loaded_object_info", + "kernel_ls", + "kernel_fetch_from_url", + "kernel_reload", + "kernel_header_comments", + "mpc_code_to_ecliptic", + "earth_pos_to_ecliptic", + "moon_illumination_frac", +] SpkInfo = namedtuple("SpkInfo", "name, jd_start, jd_end, center, frame, spk_type") @@ -55,325 +68,311 @@ def _validate_time(time: Union[float, Time]) -> float: raise TypeError("Invalid jd type, use Time or float") from exc -class SpiceKernels: +_NAME_CACHE: dict = {} + + +def state( + target: Union[str, int], + jd: Union[float, Time], + center: str = "Sun", + frame: Frames = Frames.Ecliptic, +) -> State: + """ + Calculates the :class:`~neospy.State` of the target object at the + specified time `jd`. + + This defaults to the ecliptic heliocentric state, though other centers may be + chosen. + + Parameters + ---------- + target: + The names of the target object, this can include any object name listed in + :meth:`~neospy.spice.loaded_objects` + jd: + Julian time (TDB) of the desired record. + center: + The center point, this defaults to being heliocentric. + frame: + Coordinate frame of the state, defaults to ecliptic. + + Returns + ------- + State + Returns the ecliptic state of the target in AU and AU/days. + + Raises + ------ + ValueError + If the desired time is outside of the range of the source binary file. + """ + target, ids = name_lookup(target) + center, center_id = name_lookup(center) + jd = _validate_time(jd) + return _core.spk_state(ids, jd, center_id, frame) + + +def name_lookup(name: Union[int, str]) -> tuple[str, int]: + """ + Given the provided partial name or integer, find the full name contained within + the loaded SPICE kernels. + + >>> neospy.spice.name_lookup("jupi") + ('jupiter barycenter', 5) + + >>> neospy.spice.name_lookup(10) + ('sun', 10) + + If there are multiple names, but an exact match, the exact match is returned. In + the case of ``Earth``, there is also ``Earth Barycenter``, but asking for Earth + will return the exact match. Putting ``eart`` will raise an exception as there + are 2 partial matches. + + >>> neospy.spice.name_lookup("Earth") + ('earth', 399) + + >>> neospy.spice.name_lookup("Earth b") + ('earth barycenter', 3) + + Parameters + ---------- + name : + Name, partial name, or integer id value of the object. + + Returns + ------- + tuple : + Two elements in the tuple, the full name and the integer id value. + """ + if isinstance(name, str): + name = name.lower() + if name in _NAME_CACHE: + return _NAME_CACHE[name] + + # barycenter of the solar system is special + if name == 0: + return ("ssb", 0) + + try: + lookup_name = _core.spk_get_name_from_id(int(name)) + except ValueError: + lookup_name = name + lookup_name = lookup_name.lower() + + found = [] + for loaded in loaded_objects(): + loaded_lower = loaded[0].lower() + # If it is an exact match, finish early + if lookup_name == loaded_lower: + _NAME_CACHE[name] = loaded + return loaded + if lookup_name in loaded_lower: + found.append(loaded) + found = list(set(found)) + + if len(found) == 1: + _NAME_CACHE[name] = found[0] + return found[0] + elif len(found) > 1: + raise ValueError(f"Multiple objects match this name {found}") + raise ValueError(f"No loaded objects which match this name ({name})") + + +def loaded_objects() -> list[tuple[str, int]]: + """ + Return the name of all objects which are currently loaded in the SPICE kernels. + """ + objects = _core.spk_loaded() + return [(_core.spk_get_name_from_id(o), o) for o in objects] + + +def loaded_object_info(desig: Union[int, str]) -> list[SpkInfo]: + """ + Return the available SPK information for the target object. + + Parameters + ---------- + desig : + Name or integer id value of the object. + """ + name, naif = name_lookup(desig) + return [SpkInfo(name, *k) for k in _core.spk_available_info(naif)] + + +def kernel_ls(): + """ + List all files contained within the kernels cache folder. """ - Class for estimating the ephemeris for given bodies in the solar system using SPICE - and the DE440 dataset. + path = os.path.join(cache_path(), "kernels", "**") + return glob.glob(path) - This allows for the loading of additional spice kernels along with querying. + +def kernel_fetch_from_url(url, force_download: bool = False): + """ + Download the target url into the cache folder of spice kernels. + """ + cached_file_download(url, force_download=force_download, subfolder="kernels") + + +def kernel_reload(filenames: Optional[list[str]] = None, include_cache=False): + """ + Load the specified spice kernels into memory, this resets the currently loaded + kernels. + + If `include_cache` is true, this will reload the kernels contained within the + neospy cache folder as well. + + Parameters + ---------- + filenames : + Paths to the specified files to load, this must be a list of filenames. + include_cache: + This decides if all of the files contained within the neospy cache should + be loaded in addition to the specified files. + """ + _core.spk_reset() + _core.pck_reset() + + if include_cache: + cache_files = glob.glob(os.path.join(cache_path(), "kernels", "**.bsp")) + _core.spk_load(cache_files) + + cache_files = glob.glob(os.path.join(cache_path(), "kernels", "**.bpc")) + _core.pck_load(cache_files) + + if filenames: + _core.spk_load(filenames) + + +def kernel_header_comments(filename: str): + """ + Return the comments contained within the header of the provided DAF file, this + includes SPK and PCK files. + + This does not load the contents of the file into memory, it only prints the + header contents. + + Parameters + ---------- + filename : + Path to a DAF file. + """ + return _core.daf_header_comments(filename).replace("\x04", "\n").strip() + + +def mpc_code_to_ecliptic( + obs_code: str, jd: Union[float, Time], center: str = "Sun", full_name=False +) -> State: + """ + Load an MPC Observatory code as an ecliptic state. + + This only works for ground based observatories. + + Parameters + ---------- + obs_code: + MPC observatory code or name of observatory. + jd: + Julian time (TDB) of the desired state. + center: + The new center point, this defaults to being heliocentric. + full_name: + Should the final state include the full name of the observatory or just its + code. + + Returns + ------- + State + Returns the equatorial state of the observatory in AU and AU/days. + """ + from .mpc import find_obs_code + + jd = _validate_time(jd) + + obs = find_obs_code(obs_code) + return earth_pos_to_ecliptic( + jd, + geodetic_lat=obs[0], + geodetic_lon=obs[1], + height_above_surface=obs[2], + name=obs[3] if full_name else obs[4], + center=center, + ) + + +def earth_pos_to_ecliptic( + jd: Union[float, Time], + geodetic_lat: float, + geodetic_lon: float, + height_above_surface: float, + name: Optional[str] = None, + center: str = "Sun", +) -> State: + """ + Given a position in the frame of the Earth at a specific time, convert that to + Sun centered ecliptic state. + + This uses the WGS84 model of Earth's shape to compute state. This uses Geodetic + latitude and longitude, not geocentric. + + The frame conversion is done using a high precision PCK file from the NAIF/JPL + website. The file provided in neospy is valid from ~2000 to ~2024. This file + gives positional error on the scale of a few cm. Additionally a lower resolution + file is provided which is valid until 2099, and will be used automatically when + querying past the end of the high resolution file. + + Parameters + ---------- + jd: + Julian time (TDB) of the desired state. + geodetic_lat: + Latitude on Earth's surface in degrees. + geodetic_lon: + Latitude on Earth's surface in degrees. + height_above_surface: + Height of the observer above the surface of the Earth in km. + name : + Optional name of the position on Earth. + center: + The new center point, this defaults to being heliocentric. + + Returns + ------- + State + Returns the equatorial state of the target in AU and AU/days. + """ + jd = _validate_time(jd) + pos = _core.wgs_lat_lon_to_ecef(geodetic_lat, geodetic_lon, height_above_surface) + pos = np.array(pos) / AU_KM + _, center_id = name_lookup(center) + return _core.pck_earth_frame_to_ecliptic(pos, jd, center_id, name) + + +def moon_illumination_frac(jd: Union[float, Time], observer: str = "399"): + """ + Compute the fraction of the moon which is illuminated at the specified time. + + This is a simple approximation using basic spherical geometry, and defaults to + having the observer located at the geocenter of the Earth. + + >>> neospy.spice.moon_illumination_frac(Time.from_ymd(2024, 2, 24)) + 0.9964936478732302 + + Parameters + ---------- + jd: + Julian time (TDB) of the desired state. + observer: + NAIF ID of the observer location, defaults to Earth geocenter. + + Returns + ------- + State + Fraction between 0 and 1 of the moons visible surface which is illuminated. """ + jd = _validate_time(jd) - _name_cache: dict = {} - - @classmethod - def state( - cls, - target: Union[str, int], - jd: Union[float, Time], - center: str = "Sun", - frame: Frames = Frames.Ecliptic, - ) -> State: - """ - Calculates the :class:`~neospy.State` of the target object at the - specified time `jd`. - - This defaults to the ecliptic heliocentric state, though other centers may be - chosen. - - Parameters - ---------- - target: - The names of the target object, this can include any object name listed in - :meth:`~neospy.spice.SpiceKernels.loaded_objects` - jd: - Julian time (TDB) of the desired record. - center: - The center point, this defaults to being heliocentric. - frame: - Coordinate frame of the state, defaults to ecliptic. - - Returns - ------- - State - Returns the ecliptic state of the target in AU and AU/days. - - Raises - ------ - ValueError - If the desired time is outside of the range of the source binary file. - """ - target, ids = cls.name_lookup(target) - center, center_id = cls.name_lookup(center) - jd = _validate_time(jd) - return _core.spk_state(ids, jd, center_id, frame) - - @classmethod - def name_lookup(cls, name: Union[int, str]) -> tuple[str, int]: - """ - Given the provided partial name or integer, find the full name contained within - the loaded SPICE kernels. - - >>> neospy.SpiceKernels.name_lookup("jupi") - ('jupiter barycenter', 5) - - >>> neospy.SpiceKernels.name_lookup(10) - ('sun', 10) - - If there are multiple names, but an exact match, the exact match is returned. In - the case of ``Earth``, there is also ``Earth Barycenter``, but asking for Earth - will return the exact match. Putting ``eart`` will raise an exception as there - are 2 partial matches. - - >>> neospy.SpiceKernels.name_lookup("Earth") - ('earth', 399) - - >>> neospy.SpiceKernels.name_lookup("Earth b") - ('earth barycenter', 3) - - Parameters - ---------- - name : - Name, partial name, or integer id value of the object. - - Returns - ------- - tuple : - Two elements in the tuple, the full name and the integer id value. - """ - if isinstance(name, str): - name = name.lower() - if name in cls._name_cache: - return cls._name_cache[name] - - # barycenter of the solar system is special - if name == 0: - return ("ssb", 0) - - try: - lookup_name = _core.spk_get_name_from_id(int(name)) - except ValueError: - lookup_name = name - lookup_name = lookup_name.lower() - - found = [] - for loaded in cls.loaded_objects(): - loaded_lower = loaded[0].lower() - # If it is an exact match, finish early - if lookup_name == loaded_lower: - cls._name_cache[name] = loaded - return loaded - if lookup_name in loaded_lower: - found.append(loaded) - found = list(set(found)) - - if len(found) == 1: - cls._name_cache[name] = found[0] - return found[0] - elif len(found) > 1: - raise ValueError(f"Multiple objects match this name {found}") - raise ValueError(f"No loaded objects which match this name ({name})") - - @staticmethod - def loaded_objects() -> list[tuple[str, int]]: - """ - Return the name of all objects which are currently loaded in the SPICE kernels. - """ - objects = _core.spk_loaded() - return [(_core.spk_get_name_from_id(o), o) for o in objects] - - @classmethod - def loaded_object_info(cls, desig: Union[int, str]) -> list[SpkInfo]: - """ - Return the available SPK information for the target object. - - Parameters - ---------- - desig : - Name or integer id value of the object. - """ - name, naif = cls.name_lookup(desig) - return [SpkInfo(name, *k) for k in _core.spk_available_info(naif)] - - @staticmethod - def kernel_ls(): - """ - List all files contained within the kernels cache folder. - """ - path = os.path.join(cache_path(), "kernels", "**") - return glob.glob(path) - - @staticmethod - def kernel_fetch_from_url(url, force_download: bool = False): - """ - Download the target url into the cache folder of spice kernels. - """ - cached_file_download(url, force_download=force_download, subfolder="kernels") - - @staticmethod - def kernel_reload(filenames: Optional[list[str]] = None, include_cache=False): - """ - Load the specified spice kernels into memory, this resets the currently loaded - kernels. - - If `include_cache` is true, this will reload the kernels contained within the - neospy cache folder as well. - - Parameters - ---------- - filenames : - Paths to the specified files to load, this must be a list of filenames. - include_cache: - This decides if all of the files contained within the neospy cache should - be loaded in addition to the specified files. - """ - _core.spk_reset() - _core.pck_reset() - - if include_cache: - cache_files = glob.glob(os.path.join(cache_path(), "kernels", "**.bsp")) - _core.spk_load(cache_files) - - cache_files = glob.glob(os.path.join(cache_path(), "kernels", "**.bpc")) - _core.pck_load(cache_files) - - if filenames: - _core.spk_load(filenames) - - @staticmethod - def kernel_header_comments(filename: str): - """ - Return the comments contained within the header of the provided DAF file, this - includes SPK and PCK files. - - This does not load the contents of the file into memory, it only prints the - header contents. - - Parameters - ---------- - filename : - Path to a DAF file. - """ - return _core.daf_header_comments(filename).replace("\x04", "\n").strip() - - @classmethod - def mpc_code_to_ecliptic( - cls, obs_code: str, jd: Union[float, Time], center: str = "Sun", full_name=False - ) -> State: - """ - Load an MPC Observatory code as an ecliptic state. - - This only works for ground based observatories. - - Parameters - ---------- - obs_code: - MPC observatory code or name of observatory. - jd: - Julian time (TDB) of the desired state. - center: - The new center point, this defaults to being heliocentric. - full_name: - Should the final state include the full name of the observatory or just its - code. - - Returns - ------- - State - Returns the equatorial state of the observatory in AU and AU/days. - """ - from .mpc import find_obs_code - - jd = _validate_time(jd) - - obs = find_obs_code(obs_code) - return cls.earth_pos_to_ecliptic( - jd, - geodetic_lat=obs[0], - geodetic_lon=obs[1], - height_above_surface=obs[2], - name=obs[3] if full_name else obs[4], - center=center, - ) - - @classmethod - def earth_pos_to_ecliptic( - cls, - jd: Union[float, Time], - geodetic_lat: float, - geodetic_lon: float, - height_above_surface: float, - name: Optional[str] = None, - center: str = "Sun", - ) -> State: - """ - Given a position in the frame of the Earth at a specific time, convert that to - Sun centered ecliptic state. - - This uses the WGS84 model of Earth's shape to compute state. This uses Geodetic - latitude and longitude, not geocentric. - - The frame conversion is done using a high precision PCK file from the NAIF/JPL - website. The file provided in neospy is valid from ~2000 to ~2024. This file - gives positional error on the scale of a few cm. Additionally a lower resolution - file is provided which is valid until 2099, and will be used automatically when - querying past the end of the high resolution file. - - Parameters - ---------- - jd: - Julian time (TDB) of the desired state. - geodetic_lat: - Latitude on Earth's surface in degrees. - geodetic_lon: - Latitude on Earth's surface in degrees. - height_above_surface: - Height of the observer above the surface of the Earth in km. - name : - Optional name of the position on Earth. - center: - The new center point, this defaults to being heliocentric. - - Returns - ------- - State - Returns the equatorial state of the target in AU and AU/days. - """ - jd = _validate_time(jd) - pos = _core.wgs_lat_lon_to_ecef( - geodetic_lat, geodetic_lon, height_above_surface - ) - pos = np.array(pos) / AU_KM - _, center_id = cls.name_lookup(center) - return _core.pck_earth_frame_to_ecliptic(pos, jd, center_id, name) - - @classmethod - def moon_illumination_frac(cls, jd: Union[float, Time], observer: str = "399"): - """ - Compute the fraction of the moon which is illuminated at the specified time. - - This is a simple approximation using basic spherical geometry, and defaults to - having the observer located at the geocenter of the Earth. - - >>> SpiceKernels.moon_illumination_frac(Time.from_ymd(2024, 2, 24)) - 0.9964858477054903 - - Parameters - ---------- - jd: - Julian time (TDB) of the desired state. - observer: - NAIF ID of the observer location, defaults to Earth geocenter. - - Returns - ------- - State - Fraction between 0 and 1 of the moons visible surface which is illuminated. - """ - jd = _validate_time(jd) - - moon2sun = -cls.state("moon", jd).pos - moon2earth = -cls.state("moon", jd, center=observer).pos - perc = 1.0 - moon2sun.angle_between(moon2earth) / 180 - return 0.5 - np.cos(np.pi * perc) / 2 - - __call__ = state + moon2sun = -state("moon", jd).pos + moon2earth = -state("moon", jd, center=observer).pos + perc = 1.0 - moon2sun.angle_between(moon2earth) / 180 + return 0.5 - np.cos(np.pi * perc) / 2 diff --git a/src/neospy/wise.py b/src/neospy/wise.py index a1e988d..9393828 100644 --- a/src/neospy/wise.py +++ b/src/neospy/wise.py @@ -18,9 +18,9 @@ from astropy.coordinates import SkyCoord +from . import spice from .cache import cache_path from .time import Time -from .spice import SpiceKernels from .vector import Vector from .irsa import IRSA_URL, query_irsa_tap @@ -579,7 +579,7 @@ def fetch_WISE_fovs(phase): fovs = [] for row in res.itertuples(): - state = SpiceKernels.state("WISE", row.jd) + state = spice.state("WISE", row.jd) pointing = Vector.from_ra_dec(row.ra, row.dec).as_ecliptic fov = WiseCmos( diff --git a/src/neospy/ztf.py b/src/neospy/ztf.py index b3af021..8bfc68d 100644 --- a/src/neospy/ztf.py +++ b/src/neospy/ztf.py @@ -13,7 +13,7 @@ from .irsa import query_irsa_tap from .mpc import find_obs_code from .vector import Vector, State -from .spice import SpiceKernels +from . import spice __all__ = ["fetch_ZTF_file", "fetch_ZTF_fovs"] @@ -112,7 +112,7 @@ def fetch_ZTF_fovs(year: int): ra = getattr(row, f"ra{i + 1}") dec = getattr(row, f"dec{i + 1}") corners.append(Vector.from_ra_dec(ra, dec)) - observer = SpiceKernels.earth_pos_to_ecliptic(jd, *obs_info[:-1]) + observer = spice.earth_pos_to_ecliptic(jd, *obs_info[:-1]) observer = State("ZTF", observer.jd, observer.pos, observer.vel) fov = ZtfCcdQuad( diff --git a/src/tests/test_propagation.py b/src/tests/test_propagation.py index afd281c..e6eb9df 100644 --- a/src/tests/test_propagation.py +++ b/src/tests/test_propagation.py @@ -2,7 +2,7 @@ import numpy as np from neospy import ( constants, - SpiceKernels, + spice, propagate_two_body, propagate_n_body, moid, @@ -72,27 +72,27 @@ def test_a_terms(self): class TestTwoBodyPropagation: def test_propagation_single(self): """ - Propagate Venus using StateVector and compare the result to using SpiceKernels. + Propagate Venus using StateVector and compare the result to using spice. This is only over 5 days the 2 body approximation is accurate over that time. """ - state = SpiceKernels.state("Venus", 2461161.5) + state = spice.state("Venus", 2461161.5) for jd in range(-5, 5): jd = state.jd + jd vec_state = propagate_two_body([state], jd)[0] - jpl_state = SpiceKernels.state("Venus", jd) + jpl_state = spice.state("Venus", jd) assert vec_state.jd == jpl_state.jd assert np.allclose(vec_state.vel, jpl_state.vel) assert np.allclose(vec_state.pos, jpl_state.pos) def test_propagation_light_delay(self): """ - Propagate Venus using StateVector and compare the result to using SpiceKernels. + Propagate Venus using StateVector and compare the result to using spice. This is only over 5 days the 2 body approximation is accurate over that time. Place an observer X AU away from Venus and ensure that the delay is correct. """ - state = SpiceKernels.state("Venus", 2461161.5) + state = spice.state("Venus", 2461161.5) for au in range(0, 5): sun2obs = Vector(state.pos + [au, 0.0, 0.0]) @@ -108,9 +108,9 @@ def test_propagation_light_delay(self): def test_moid(planet): if planet is None: state = None - vs = SpiceKernels.state("Earth", 2461161.5) + vs = spice.state("Earth", 2461161.5) else: - state = SpiceKernels.state(planet, 2461161.5) + state = spice.state(planet, 2461161.5) vs = state assert np.isclose(moid(vs, state), 0) diff --git a/src/tests/test_spice.py b/src/tests/test_spice.py index a4a45ea..8e40e20 100644 --- a/src/tests/test_spice.py +++ b/src/tests/test_spice.py @@ -1,7 +1,7 @@ import numpy as np import pytest -from neospy import SpiceKernels, State, Time, Frames +from neospy import spice, State, Time, Frames from neospy.mpc import find_obs_code from neospy.spice import SpkInfo @@ -25,28 +25,28 @@ # fmt: on -class TestSpiceKernels: +class TestSpice: @pytest.mark.parametrize("expected", EXPECTED_EQUITORIAL) def test_ecliptic_state(self, expected): e_pos = expected[2:5] e_vel = expected[5:8] - state = SpiceKernels.state(expected[0], expected[1]) + state = spice.state(expected[0], expected[1]) assert isinstance(state, State) assert np.allclose(state.pos.as_equatorial, e_pos) assert np.allclose(state.vel.as_equatorial, e_vel) # Ensure that the same values are calculated if using an astropy time jd = Time(expected[1]) - state = SpiceKernels.state(expected[0], jd) + state = spice.state(expected[0], jd) assert isinstance(state, State) assert np.allclose(state.pos.as_equatorial, e_pos) assert np.allclose(state.vel.as_equatorial, e_vel) def test_loaded_objects(self): - assert len(SpiceKernels.loaded_objects()) > 8 + assert len(spice.loaded_objects()) > 8 def test_earth_to_ecliptic(self): - state = SpiceKernels.earth_pos_to_ecliptic(2454832.5, 0, 0, 0, center="399") + state = spice.earth_pos_to_ecliptic(2454832.5, 0, 0, 0, center="399") assert np.allclose( list(-state.pos), [7.686364937857906e-06, -3.847842319488030e-05, 1.667609317250255e-05], @@ -57,17 +57,17 @@ def test_earth_to_ecliptic(self): ) def test_name_lookup(self): - assert SpiceKernels.name_lookup("jupi") == ("jupiter barycenter", 5) - assert SpiceKernels.name_lookup(0) == ("ssb", 0) + assert spice.name_lookup("jupi") == ("jupiter barycenter", 5) + assert spice.name_lookup(0) == ("ssb", 0) with pytest.raises(ValueError, match="Multiple objects"): - SpiceKernels.name_lookup("ear") + spice.name_lookup("ear") with pytest.raises(ValueError, match="No loaded objects"): - SpiceKernels.name_lookup("banana") + spice.name_lookup("banana") def test_loaded_info(self): - info = SpiceKernels.loaded_object_info("ceres") + info = spice.loaded_object_info("ceres") expected = SpkInfo( name="ceres", jd_start=2415020.5, @@ -79,27 +79,27 @@ def test_loaded_info(self): assert info == [expected] def test_cache(self): - SpiceKernels.kernel_ls() + spice.kernel_ls() def test_reload(self): - SpiceKernels.kernel_reload([], include_cache=True) + spice.kernel_reload([], include_cache=True) def test_mpc_code(self): jd = Time.j2000().jd - state0 = SpiceKernels.mpc_code_to_ecliptic("Palomar Mountain", jd) + state0 = spice.mpc_code_to_ecliptic("Palomar Mountain", jd) code = find_obs_code("Palomar Mountain") - state1 = SpiceKernels.earth_pos_to_ecliptic(jd, *code[:-1]) + state1 = spice.earth_pos_to_ecliptic(jd, *code[:-1]) assert np.isclose((state0.pos - state1.pos).r, 0.0) assert np.isclose((state0.vel - state1.vel).r, 0.0) def test_moon_frac(self): jd = Time.j2000().jd - frac = SpiceKernels.moon_illumination_frac(jd) + frac = spice.moon_illumination_frac(jd) # matches JPL Horizons to 4 decimal places # Horizons reports: 0.230064 assert np.isclose(frac, 0.23018685933) # Horizons reports: 0.200455 - assert np.isclose(SpiceKernels.moon_illumination_frac(2451555), 0.2003318) + assert np.isclose(spice.moon_illumination_frac(2451555), 0.2003318)