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

Use integer time throughout codebase #73

Merged
merged 15 commits into from
Oct 10, 2023
Merged

Use integer time throughout codebase #73

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
be47a92
Add utilities for converting to TDB and from MJD
spenczar Sep 26, 2023
940b220
Add Timestamp.to_numpy
spenczar Sep 27, 2023
d42b79a
Use integer Timestamps throughout codebase
spenczar Sep 27, 2023
89e21dc
Rewrite observers.tests.test_state
spenczar Oct 2, 2023
90007a8
Re-delete coordinates/io.py
spenczar Oct 2, 2023
a964474
Add Timestamp.unique
spenczar Oct 3, 2023
c9c9f79
Use integer time throughout SPICE and observer state code
spenczar Oct 3, 2023
19ac3dd
Replace trivial uses of astropy with direct construction of Timestamps
spenczar Oct 3, 2023
e5dcc63
Update benchmarks to use integer time
spenczar Oct 3, 2023
a82e444
Add Timestamp.rescale
spenczar Oct 4, 2023
0bf9f32
Use integer time in all interfaces, including propagators
spenczar Oct 5, 2023
24a8430
Remove straggling places where astropy time is used
spenczar Oct 10, 2023
dda3e18
Fix time import
spenczar Oct 10, 2023
382301c
Delete empty file
spenczar Oct 10, 2023
7b5a46f
Fix up mishmash that was missed in rebase
spenczar Oct 10, 2023
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
49 changes: 18 additions & 31 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,17 +14,13 @@

To define an orbit:
```python
from astropy.time import Time

from adam_core.coordinates import KeplerianCoordinates
from adam_core.coordinates import Times
from adam_core.coordinates import Origin
from adam_core.orbits import Orbits
from adam_core.time import Timestamp

keplerian_elements = KeplerianCoordinates.from_kwargs(
time=Times.from_astropy(
Time([59000.0], scale="tdb", format="mjd")
),
time=Timestamp.from_mjd([59000.0], scale="tdb"),
a=[1.0],
e=[0.002],
i=[10.],
Expand All @@ -47,17 +43,13 @@ can be done on demand by calling `to_cartesian()` on the coordinates object.
The underlying orbits class is 2 dimensional and can store elements and covariances for multiple orbits.

```python
from astropy.time import Time

from adam_core.coordinates import KeplerianCoordinates
from adam_core.coordinates import Times
from adam_core.coordinates import Origin
from adam_core.orbits import Orbits
from adam_core.time import Timestamp

keplerian_elements = KeplerianCoordinates.from_kwargs(
time=Times.from_astropy(
Time([59000.0, 60000.0], scale="tdb", format="mjd")
),
time=Timestamp.from_mjd([59000.0, 60000.0], scale="tdb"),
a=[1.0, 3.0],
e=[0.002, 0.0],
i=[10., 30.],
Expand Down Expand Up @@ -85,18 +77,14 @@ orbits.to_dataframe()
Orbits can also be defined with uncertainties.
```python
import numpy as np
from astropy.time import Time

from adam_core.coordinates import KeplerianCoordinates
from adam_core.coordinates import Times
from adam_core.coordinates import Origin
from adam_core.coordinates import CoordinateCovariances
from adam_core.orbits import Orbits
from adam_core.time import Timestamp

keplerian_elements = KeplerianCoordinates.from_kwargs(
time=Times.from_astropy(
Time([59000.0], scale="tdb", format="mjd")
),
time=Timestamp.from_mjd([59000.0], scale="tdb"),
a=[1.0],
e=[0.002],
i=[10.],
Expand All @@ -122,11 +110,10 @@ orbits.to_dataframe(sigmas=True)

To query orbits from JPL Horizons:
```python
from astropy.time import Time

from adam_core.orbits.query import query_horizons
from adam_core.time import Timestamp

times = Time([60000.0], scale="tdb", format="mjd")
times = Timestamp.from_mjd([60000.0], scale="tdb")
object_ids = ["Duende", "Eros", "Ceres"]
orbits = query_horizons(object_ids, times)
```
Expand Down Expand Up @@ -160,22 +147,22 @@ The propagator class in `adam_core` provides a generalized interface to the supp
To propagate orbits with PYOORB (here we grab some orbits from Horizons first):
```python
import numpy as np
from astropy.time import Time
from astropy import units as u

from adam_core.orbits.query import query_horizons
from adam_core.propagator import PYOORB
from adam_core.time import Timestamp

# Get orbits to propagate
initial_time = Time([60000.0], scale="tdb", format="mjd")
initial_time = Timestamp.from_mjd([60000.0], scale="tdb")
object_ids = ["Duende", "Eros", "Ceres"]
orbits = query_horizons(object_ids, initial_time)

# Make sure PYOORB is ready
propagator = PYOORB()

# Define propagation times
times = initial_time + np.arange(0, 100) * u.d
times = initial_time.from_mjd(initial_time.mjd() + np.arange(0, 100))

# Propagate orbits! This function supports multiprocessing for large
# propagation jobs.
Expand All @@ -192,23 +179,23 @@ The propagator class can also be used to generate ephemerides for a set of orbit

```python
import numpy as np
from astropy.time import Time
from astropy import units as u

from adam_core.orbits.query import query_horizons
from adam_core.propagator import PYOORB
from adam_core.observers import Observers
from adam_core.time import Timestamp

# Get orbits to propagate
initial_time = Time([60000.0], scale="tdb", format="mjd")
initial_time = Timestamp.from_mjd([60000.0], scale="tdb")
object_ids = ["Duende", "Eros", "Ceres"]
orbits = query_horizons(object_ids, initial_time)

# Make sure PYOORB is ready
propagator = PYOORB()

# Define a set of observers and observation times
times = initial_time + np.arange(0, 100) * u.d
times = Timestamp.from_mjd(initial_time.mjd() + np.arange(0, 100))
observers = Observers.from_code("I11", times)

# Generate ephemerides! This function supports multiprocessing for large
Expand All @@ -228,14 +215,14 @@ ephemeris = propagator.generate_ephemeris(
Getting the heliocentric ecliptic state vector of a DE440 body at a given set of times (in this case the barycenter of the Jovian system):
```python
import numpy as np
from astropy.time import Time

from adam_core.coordinates import OriginCodes
from adam_core.utils import get_perturber_state
from adam_core.time import Timestamp

states = get_perturber_state(
OriginCodes.JUPITER_BARYCENTER,
Time(np.arange(59000, 60000), format="mjd", scale="tdb"),
Timetamp.from_mjd(np.arange(59000, 60000), scale="tdb"),
frame="ecliptic",
origin=OriginCodes.SUN,
)
Expand All @@ -245,18 +232,18 @@ states = get_perturber_state(
`adam_core` also has 2-body propagation functionality. To propagate any orbit with 2-body dynamics:
```python
import numpy as np
from astropy.time import Time
from astropy import units as u

from adam_core.orbits.query import query_sbdb
from adam_core.dynamics import propagate_2body
from adam_core.time import Timestamp

# Get orbit to propagate
object_ids = ["Duende", "Eros", "Ceres"]
orbits = query_sbdb(object_ids)

# Define propagation times
times = Time(np.arange(59000, 60000), scale="tdb", format="mjd")
times = Timestamp.from_mjd(np.arange(59000, 60000), scale="tdb")

# Propagate orbits with 2-body dynamics
propagated_orbits = propagate_2body(
Expand Down
4 changes: 2 additions & 2 deletions adam_core/coordinates/cartesian.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,9 +5,9 @@
import quivr as qv
from astropy import units as u

from ..time import Timestamp
from .covariances import CoordinateCovariances
from .origin import Origin
from .times import Times

if TYPE_CHECKING:
from .cometary import CometaryCoordinates
Expand Down Expand Up @@ -37,7 +37,7 @@ class CartesianCoordinates(qv.Table):
vx = qv.Float64Column(nullable=True)
vy = qv.Float64Column(nullable=True)
vz = qv.Float64Column(nullable=True)
time = Times.as_column(nullable=True)
time = Timestamp.as_column(nullable=True)
covariance = CoordinateCovariances.as_column(nullable=True)
origin = Origin.as_column()
frame = qv.StringAttribute()
Expand Down
12 changes: 6 additions & 6 deletions adam_core/coordinates/cometary.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,10 +4,10 @@
import quivr as qv
from astropy import units as u

from ..time import Timestamp
from .cartesian import CartesianCoordinates
from .covariances import CoordinateCovariances, transform_covariances_jacobian
from .origin import Origin
from .times import Times

if TYPE_CHECKING:
from .keplerian import KeplerianCoordinates
Expand Down Expand Up @@ -44,7 +44,7 @@ class CometaryCoordinates(qv.Table):
raan = qv.Float64Column()
ap = qv.Float64Column()
tp = qv.Float64Column()
time = Times.as_column(nullable=True)
time = Timestamp.as_column(nullable=True)
covariance = CoordinateCovariances.as_column(nullable=True)
origin = Origin.as_column()
frame = qv.StringAttribute()
Expand Down Expand Up @@ -203,7 +203,7 @@ def to_cartesian(self) -> CartesianCoordinates:

coords_cartesian = cometary_to_cartesian(
self.values,
t0=self.time.to_astropy().tdb.mjd,
t0=self.time.to_numpy(),
mu=mu,
max_iter=100,
tol=1e-15,
Expand All @@ -218,7 +218,7 @@ def to_cartesian(self) -> CartesianCoordinates:
_cometary_to_cartesian,
in_axes=(0, 0, None, None, None),
out_axes=0,
t0=self.time.to_astropy().tdb.mjd,
t0=self.time.to_numpy(),
mu=mu,
max_iter=100,
tol=1e-15,
Expand Down Expand Up @@ -262,7 +262,7 @@ def from_cartesian(cls, cartesian: CartesianCoordinates) -> "CometaryCoordinates

coords_cometary = cartesian_to_cometary(
cartesian.values,
cartesian.time.to_astropy().tdb.mjd,
cartesian.time.to_numpy(),
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've decided to add to_numpy as a convenience method for the conversion to TDD-scaled, MJD-epoched, double-precisioned timestamps. Perhaps the name is opaque?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think to_numpy is a good addition. Though I agree that the name is a little opaque. Do you think it might make sense to add a scale kwarg that will return the numpy array in the desired scale (to_numpy(scale="tdb"))? The default can be set to its current scale?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Another option (now having read more of the code) might be to enable symmetry between .from_mjd by adding a to_mjd(scale="tdb"). The current mjd() could suffice as well with a scale parameter. The latter could just return a pyarrow array (the additional call to .to_numpy() from there doesn't seem super cumbersome from my perspective).

I suspect a lot of astronomers will expect an interface similar to astropy. Don't know how much we want to mimic it but that's a thought. This doesn't have to be addressed in this PR obviously. Here are some options:

import numpy as np
from adam_core.time import Timestamp

mjd_tdb = np.arange(59000, 60000)

times = Timestamp.from_mjd(mjd_tdb, scale="tdb") # this feels very nice and not much different than astropy
mjd_utc = times.utc().mjd().to_numpy() 
mjd_utc = times.mjd(scale="utc").to_numpy() 
mjd_utc = times.to_numpy(scale="utc")
mjd_utc = times.to_mjd(scale="utc") # the symmetry between .from_mjd() and .to_mjd() feels quite intuitive

mu=mu,
)
coords_cometary = np.array(coords_cometary)
Expand All @@ -275,7 +275,7 @@ def from_cartesian(cls, cartesian: CartesianCoordinates) -> "CometaryCoordinates
_cartesian_to_cometary,
in_axes=(0, 0, None),
out_axes=0,
t0=cartesian.time.to_astropy().tdb.mjd,
t0=cartesian.time.to_numpy(),
mu=mu,
)
else:
Expand Down
8 changes: 4 additions & 4 deletions adam_core/coordinates/keplerian.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,10 +4,10 @@
import quivr as qv
from astropy import units as u

from ..time import Timestamp
from .cartesian import CartesianCoordinates
from .covariances import CoordinateCovariances, transform_covariances_jacobian
from .origin import Origin
from .times import Times

if TYPE_CHECKING:
from .cometary import CometaryCoordinates
Expand Down Expand Up @@ -40,7 +40,7 @@ class KeplerianCoordinates(qv.Table):
raan = qv.Float64Column()
ap = qv.Float64Column()
M = qv.Float64Column()
time = Times.as_column(nullable=True)
time = Timestamp.as_column(nullable=True)
covariance = CoordinateCovariances.as_column(nullable=True)
origin = Origin.as_column()
frame = qv.StringAttribute()
Expand Down Expand Up @@ -240,7 +240,7 @@ def from_cartesian(cls, cartesian: CartesianCoordinates):

coords_keplerian = cartesian_to_keplerian(
cartesian.values,
cartesian.time.to_astropy().tdb.mjd,
cartesian.time.to_numpy(),
mu=mu,
)
coords_keplerian = np.array(coords_keplerian)
Expand All @@ -253,7 +253,7 @@ def from_cartesian(cls, cartesian: CartesianCoordinates):
_cartesian_to_keplerian6,
in_axes=(0, 0, None),
out_axes=0,
t0=cartesian.time.to_astropy().tdb.mjd,
t0=cartesian.time.to_numpy(),
mu=mu,
)
else:
Expand Down
4 changes: 2 additions & 2 deletions adam_core/coordinates/spherical.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,10 +4,10 @@
import quivr as qv
from astropy import units as u

from ..time import Timestamp
from .cartesian import CartesianCoordinates
from .covariances import CoordinateCovariances, transform_covariances_jacobian
from .origin import Origin
from .times import Times

if TYPE_CHECKING:
from .cometary import CometaryCoordinates
Expand Down Expand Up @@ -40,7 +40,7 @@ class SphericalCoordinates(qv.Table):
vrho = qv.Float64Column(nullable=True)
vlon = qv.Float64Column(nullable=True)
vlat = qv.Float64Column(nullable=True)
time = Times.as_column(nullable=True)
time = Timestamp.as_column(nullable=True)
covariance = CoordinateCovariances.as_column(nullable=True)
origin = Origin.as_column()
frame = qv.StringAttribute()
Expand Down
8 changes: 3 additions & 5 deletions adam_core/coordinates/tests/test_benchmarks.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,14 +1,13 @@
import numpy as np
import pytest
from astropy.time import Time

from ...time import Timestamp
from ..cartesian import CartesianCoordinates
from ..cometary import CometaryCoordinates
from ..covariances import CoordinateCovariances
from ..keplerian import KeplerianCoordinates
from ..origin import Origin, OriginCodes
from ..spherical import SphericalCoordinates
from ..times import Times
from ..transform import transform_coordinates


Expand Down Expand Up @@ -39,10 +38,9 @@ def test_benchmark_transform_cartesian_coordinates(
vx=np.array([1]),
vy=np.array([1]),
vz=np.array([1]),
time=Times.from_astropy(
Time([50000], format="mjd"),
),
time=Timestamp.from_mjd([50000]),
origin=Origin.from_kwargs(code=["SUN"]),
frame="ecliptic",
)
benchmark(
transform_coordinates,
Expand Down
30 changes: 10 additions & 20 deletions adam_core/coordinates/tests/test_transforms_rotation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,8 +1,6 @@
from astropy.time import Time

from ...time import Timestamp
from ..cartesian import CartesianCoordinates
from ..origin import Origin
from ..times import Times
from ..transform import cartesian_to_frame, transform_coordinates
from .test_transforms_translation import assert_coords_equal

Expand All @@ -18,9 +16,7 @@ def test_cartesian_to_frame(orbital_elements, orbital_elements_equatorial):
vx=orbital_elements["vx"].values,
vy=orbital_elements["vy"].values,
vz=orbital_elements["vz"].values,
time=Times.from_astropy(
Time(orbital_elements["mjd_tdb"].values, format="mjd", scale="tdb")
),
time=Timestamp.from_mjd(orbital_elements["mjd_tdb"].values, scale="tdb"),
frame="ecliptic",
origin=Origin.from_kwargs(code=["SUN" for i in range(len(orbital_elements))]),
)
Expand All @@ -32,12 +28,8 @@ def test_cartesian_to_frame(orbital_elements, orbital_elements_equatorial):
vx=orbital_elements_equatorial["vx"].values,
vy=orbital_elements_equatorial["vy"].values,
vz=orbital_elements_equatorial["vz"].values,
time=Times.from_astropy(
Time(
orbital_elements_equatorial["mjd_tdb"].values,
format="mjd",
scale="tdb",
)
time=Timestamp.from_mjd(
orbital_elements_equatorial["mjd_tdb"].values, scale="tdb"
),
frame="equatorial",
origin=Origin.from_kwargs(
Expand Down Expand Up @@ -81,8 +73,9 @@ def test_transform_coordinates_frame(orbital_elements, orbital_elements_equatori
vx=orbital_elements["vx"].values,
vy=orbital_elements["vy"].values,
vz=orbital_elements["vz"].values,
time=Times.from_astropy(
Time(orbital_elements["mjd_tdb"].values, format="mjd", scale="tdb")
time=Timestamp.from_mjd(
orbital_elements["mjd_tdb"].values,
scale="tdb",
),
frame="ecliptic",
origin=Origin.from_kwargs(code=["SUN" for i in range(len(orbital_elements))]),
Expand All @@ -95,12 +88,9 @@ def test_transform_coordinates_frame(orbital_elements, orbital_elements_equatori
vx=orbital_elements_equatorial["vx"].values,
vy=orbital_elements_equatorial["vy"].values,
vz=orbital_elements_equatorial["vz"].values,
time=Times.from_astropy(
Time(
orbital_elements_equatorial["mjd_tdb"].values,
format="mjd",
scale="tdb",
)
time=Timestamp.from_mjd(
orbital_elements_equatorial["mjd_tdb"].values,
scale="tdb",
),
frame="equatorial",
origin=Origin.from_kwargs(
Expand Down
Loading