Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Modernize #136

Merged
merged 4 commits into from
Dec 14, 2022
Merged

Modernize #136

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
4942382
Migrate to current Horizons web API
taldcroft Dec 13, 2022
da2929a
Allow any valid Horizons body in get_planet_chandra_horizons
taldcroft Dec 13, 2022
5dc97ff
Update docs and add type annotations
taldcroft Dec 13, 2022
fd8187b
Better exception handling and tests
taldcroft Dec 13, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
61 changes: 43 additions & 18 deletions chandra_aca/planets.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,12 +22,13 @@
"""
from datetime import datetime
from pathlib import Path
from typing import Optional, Union

import astropy.constants as const
import astropy.units as u
import numpy as np
from astropy.io import ascii
from cxotime import CxoTime
from cxotime import CxoTime, CxoTimeLike
from ska_helpers.utils import LazyVal

from chandra_aca.transform import eci_to_radec
Expand Down Expand Up @@ -93,7 +94,7 @@ def load_kernel():
("pluto", [(0, 9)]),
]
)
URL_HORIZONS = "https://ssd.jpl.nasa.gov/horizons_batch.cgi?"
URL_HORIZONS = "https://ssd.jpl.nasa.gov/api/horizons.api?"


def get_planet_angular_sep(
Expand Down Expand Up @@ -153,11 +154,14 @@ def get_planet_angular_sep(
return sep


def get_planet_barycentric(body, time=None):
def get_planet_barycentric(body: str, time: CxoTimeLike = None):
"""Get barycentric position for solar system ``body`` at ``time``.

This uses the built-in JPL ephemeris file DE432s and jplephem.

``body`` must be one of "sun", "mercury", "venus", "earth-moon-barycenter", "earth",
"moon", "mars", "jupiter", "saturn", "uranus", "neptune", or "pluto".

:param body: Body name (lower case planet name)
:param time: Time or times for returned position (default=NOW)
:returns: barycentric position (km) as (x, y, z) or N x (x, y, z)
Expand All @@ -178,13 +182,19 @@ def get_planet_barycentric(body, time=None):
return pos.transpose() # SPK returns (3, N) but we need (N, 3)


def get_planet_eci(body, time=None, pos_observer=None):
def get_planet_eci(
body: str, time: CxoTimeLike = None, pos_observer: Optional[str] = None
):
"""Get ECI apparent position for solar system ``body`` at ``time``.

This uses the built-in JPL ephemeris file DE432s and jplephem. The position
is computed at the supplied ``time`` minus the light-travel time from the
observer to ``body`` to generate the apparent position on Earth at ``time``.

``body`` and ``pos_observer`` must be one of "sun", "mercury", "venus",
"earth-moon-barycenter", "earth", "moon", "mars", "jupiter", "saturn", "uranus",
"neptune", or "pluto".

Estimated accuracy of planet coordinates (RA, Dec) is as follows, where the
JPL Horizons positions are used as the "truth". This assumes the observer
position is Earth (default).
Expand Down Expand Up @@ -214,7 +224,7 @@ def get_planet_eci(body, time=None, pos_observer=None):
return pos_planet - pos_observer


def get_planet_chandra(body, time=None):
def get_planet_chandra(body: str, time: CxoTimeLike = None):
"""Get position for solar system ``body`` at ``time`` relative to Chandra.

This uses the built-in JPL ephemeris file DE432s and jplephem, along with
Expand Down Expand Up @@ -267,18 +277,27 @@ def get_planet_chandra(body, time=None):
return planet_chandra


def get_planet_chandra_horizons(body, timestart, timestop, n_times=10, timeout=10):
"""Get body position, rate, mag, surface brightness, diameter from Horizons.
def get_planet_chandra_horizons(
body: Union[str, int],
timestart: CxoTimeLike,
timestop: CxoTimeLike,
n_times: int = 10,
timeout: float = 10,
):
"""Get body position and other info as seen from Chandra using JPL Horizons.

In addition to the planet names, the ``body`` argument can be any identifier that
Horizon supports, e.g. ``sun`` or ``geo`` (Earth geocenter).

This function queries the JPL Horizons site using the CGI batch interface
(See https://ssd.jpl.nasa.gov/horizons_batch.cgi for docs).
This function queries the JPL Horizons site using the web API interface
(See https://ssd-api.jpl.nasa.gov/doc/horizons.html for docs).

The return value is an astropy Table with columns: time, ra, dec, rate_ra,
rate_dec, mag, surf_brt, ang_diam. The units are included in the table
columns. The ``time`` column is a ``CxoTime`` object.

The returned Table has a meta key value ``response_text`` with the full text
of the Horizons response.
of the Horizons response and a ``response_json`` key with the parsed JSON.

Example::

Expand All @@ -297,7 +316,7 @@ def get_planet_chandra_horizons(body, timestart, timestop, n_times=10, timeout=1
2020:002:00:00:00.000 277.21454 -23.18970 34.69 2.51 -1.839 5.408 31.76

:param body: one of 'mercury', 'venus', 'mars', 'jupiter', 'saturn',
'uranus', 'neptune'
'uranus', 'neptune', or any other body that Horizons supports.
:param timestart: start time (any CxoTime-compatible time)
:param timestop: stop time (any CxoTime-compatible time)
:param n_times: number of time samples (inclusive, default=10)
Expand All @@ -319,11 +338,9 @@ def get_planet_chandra_horizons(body, timestart, timestop, n_times=10, timeout=1
"uranus": "799",
"neptune": "899",
}
if body not in planet_ids:
raise ValueError(f"body must be one of {tuple(planet_ids)}")

params = dict(
COMMAND=planet_ids[body],
COMMAND=planet_ids.get(body, str(body).lower()),
MAKE_EPHEM="YES",
CENTER="@-151",
TABLE_TYPE="OBSERVER",
Expand All @@ -336,18 +353,24 @@ def get_planet_chandra_horizons(body, timestart, timestop, n_times=10, timeout=1
)

# The HORIZONS web API seems to require all params to be quoted strings.
# See: https://ssd.jpl.nasa.gov/horizons_batch.cgi
# See: https://ssd-api.jpl.nasa.gov/doc/horizons.html
for key, val in params.items():
params[key] = repr(val)
params["batch"] = 1
resp = requests.get(URL_HORIZONS, params=params, timeout=timeout)

if resp.status_code != requests.codes["ok"]:
raise ValueError(
"request {resp.url} failed: {resp.reason} ({resp.status_code})"
f"request {resp.url} failed: {resp.reason} ({resp.status_code})"
)

lines = resp.text.splitlines()
resp_json: dict = resp.json()
result: str = resp_json["result"]
lines = result.splitlines()

if "$$SOE" not in lines:
msg = "problem with Horizons query:\n" + "\n".join(lines)
raise ValueError(msg)

idx0 = lines.index("$$SOE") + 1
idx1 = lines.index("$$EOE")
lines = lines[idx0:idx1]
Expand All @@ -368,6 +391,7 @@ def get_planet_chandra_horizons(body, timestart, timestop, n_times=10, timeout=1
"ang_diam",
"null3",
],
fill_values=[("n.a.", "0")],
)

times = [datetime.strptime(val[:20], "%Y-%b-%d %H:%M:%S") for val in dat["time"]]
Expand All @@ -390,6 +414,7 @@ def get_planet_chandra_horizons(body, timestart, timestop, n_times=10, timeout=1
dat["ang_diam"].info.format = ".2f"

dat.meta["response_text"] = resp.text
dat.meta["response_json"] = resp_json

del dat["null1"]
del dat["null2"]
Expand Down
21 changes: 19 additions & 2 deletions chandra_aca/tests/test_planets.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -114,8 +114,9 @@ def test_planet_positions_array():


@pytest.mark.skipif(not HAS_INTERNET, reason="Requires network access")
def test_get_chandra_planet_horizons():
dat = get_planet_chandra_horizons("jupiter", "2020:001", "2020:002", n_times=11)
@pytest.mark.parametrize("body", ["jupiter", 599])
def test_get_chandra_planet_horizons(body):
dat = get_planet_chandra_horizons(body, "2020:001", "2020:002", n_times=11)
exp = [
" time ra dec rate_ra rate_dec mag "
" surf_brt ang_diam",
Expand Down Expand Up @@ -150,6 +151,22 @@ def test_get_chandra_planet_horizons():
assert dat.pformat_all() == exp


@pytest.mark.skipif(not HAS_INTERNET, reason="Requires network access")
def test_get_chandra_planet_horizons_non_planet():
dat = get_planet_chandra_horizons(
"ACE (spacecraft)", "2020:001", "2020:002", n_times=2
)
exp = [
" ra dec rate_ra rate_dec mag surf_brt ang_diam",
" deg deg arcsec / h arcsec / h mag mag / arcsec2 arcsec ",
"--------- --------- ---------- ---------- --- ------------- --------",
"274.14524 -24.72559 10.75 -325.04 -- -- --",
"275.29375 -23.83915 349.43 659.97 -- -- --",
]
del dat["time"]
assert dat.pformat_all() == exp


@pytest.mark.skipif(not HAS_INTERNET, reason="Requires network access")
@pytest.mark.parametrize(
"obs_pos,exp_sep", [("chandra-horizons", 0.0), ("chandra", 0.74), ("earth", 23.02)]
Expand Down