diff --git a/README.md b/README.md
index c46e0c5c..edcf6634 100644
--- a/README.md
+++ b/README.md
@@ -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.],
@@ -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.],
@@ -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.],
@@ -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)
 ```
@@ -160,14 +147,14 @@ 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)
 
@@ -175,7 +162,7 @@ orbits = query_horizons(object_ids, initial_time)
 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.
@@ -192,15 +179,15 @@ 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)
 
@@ -208,7 +195,7 @@ orbits = query_horizons(object_ids, initial_time)
 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
@@ -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,
 )
@@ -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(
diff --git a/adam_core/dynamics/tests/test_ephemeris.py b/adam_core/dynamics/tests/test_ephemeris.py
index 62be0523..97c30581 100644
--- a/adam_core/dynamics/tests/test_ephemeris.py
+++ b/adam_core/dynamics/tests/test_ephemeris.py
@@ -2,10 +2,10 @@
 import pytest
 import quivr as qv
 from astropy import units as u
-from astropy.time import Time
 
 from ...observers import Observers
 from ..ephemeris import generate_ephemeris_2body
+from ..time import Timestamp
 
 OBJECT_IDS = [
     "594913 'Aylo'chaxnim (2020 AV2)",
@@ -81,9 +81,8 @@ def test_generate_ephemeris_2body(object_id, propagated_orbits, ephemeris):
         observatory_mask = ephemeris_orbit["observatory_code"] == observatory_code
         observer_i = Observers.from_code(
             observatory_code,
-            Time(
+            Timestamp.from_mjd(
                 ephemeris_orbit[observatory_mask]["mjd_utc"].values,
-                format="mjd",
                 scale="utc",
             ),
         )
diff --git a/adam_core/observers/tests/testdata/README.md b/adam_core/observers/tests/testdata/README.md
index 6a8d0166..9f3784f7 100644
--- a/adam_core/observers/tests/testdata/README.md
+++ b/adam_core/observers/tests/testdata/README.md
@@ -5,13 +5,12 @@ generated with the following code:
 
 ```python
 import numpy as np
-from astropy.time import Time
 from astroquery.jplhorizons import Horizons
 
 from adam_core.coordinates.cartesian import CartesianCoordinates
 from adam_core.coordinates.covariances import CoordinateCovariances
 from adam_core.coordinates.origin import Origin, OriginCodes
-from adam_core.coordinates.times import Times
+from adam_core.time import Timestamp
 
 observatory_codes = ["I41", "X05", "F51", "W84", "000", "500"]
 times = Time(
@@ -32,7 +31,7 @@ for code in observatory_codes:
 
         # Flip the signs of the state to get the state of the observer
         states = CartesianCoordinates.from_kwargs(
-            time=Times.from_astropy(Time(result["datetime_jd"].values, format="jd", scale="tdb")),
+            time=Timestamp.from_jd(result["datetime_jd"].values, scale="tdb"),
             x=-result["x"].values,
             y=-result["y"].values,
             z=-result["z"].values,
diff --git a/adam_core/orbits/query/horizons.py b/adam_core/orbits/query/horizons.py
index 0d31f88e..df883a55 100644
--- a/adam_core/orbits/query/horizons.py
+++ b/adam_core/orbits/query/horizons.py
@@ -31,8 +31,8 @@ def _get_horizons_vectors(
     ----------
     object_ids : Union[List, `~numpy.ndarray`] (N)
         Object IDs / designations recognizable by HORIZONS.
-    times : `~astropy.core.time.Time` (M)
-        Astropy time object at which to gather state vectors.
+    times : Timestamp (M)
+        Time at which to gather state vectors.
     location : str, optional
         Location of the origin typically a NAIF code.
         ('0' or '@ssb' for solar system barycenter, '10' or '@sun' for heliocenter)
@@ -89,8 +89,8 @@ def _get_horizons_elements(
     ----------
     object_ids : Union[List, `~numpy.ndarray`] (N)
         Object IDs / designations recognizable by HORIZONS.
-    times : `~astropy.core.time.Time`
-        Astropy time object at which to gather state vectors.
+    times : Timestamp
+        Time at which to gather state vectors.
     location : str, optional
         Location of the origin typically a NAIF code.
         ('0' or '@ssb' for solar system barycenter, '10' or '@sun' for heliocenter)
@@ -144,8 +144,8 @@ def _get_horizons_ephemeris(
     ----------
     object_ids : Union[List, `~numpy.ndarray`] (N)
         Object IDs / designations recognizable by HORIZONS.
-    times : `~astropy.core.time.Time`
-        Astropy time object at which to gather state vectors.
+    times : Timestamp
+        Time at which to gather state vectors.
     location : str, optional
         Location of the origin typically a NAIF code or MPC observatory code
     id_type : {'majorbody', 'smallbody', 'designation',
@@ -242,7 +242,7 @@ def query_horizons(
     object_ids : npt.ArrayLike (N)
         Object IDs / designations recognizable by HORIZONS.
     times : Timestamp (M)
-        Astropy time object at which to gather state vectors.
+        Time at which to gather state vectors.
     coordinate_type : {'cartesian', 'keplerian', 'cometary'}
         Type of orbital elements to return.
     location : str, optional
diff --git a/adam_core/utils/helpers/data/get_test_data.py b/adam_core/utils/helpers/data/get_test_data.py
index ba73c6da..ff0513be 100644
--- a/adam_core/utils/helpers/data/get_test_data.py
+++ b/adam_core/utils/helpers/data/get_test_data.py
@@ -5,7 +5,6 @@
 import pandas as pd
 import pyarrow as pa
 import quivr as qv
-from astropy.time import Time
 
 from adam_core.observers import Observers
 from adam_core.orbits import Orbits
@@ -14,10 +13,11 @@
     _get_horizons_elements,
     _get_horizons_vectors,
 )
+from adam_core.time import Timestamp
 
 
 def _get_orbital_elements(
-    object_ids: List[str], time: Time, location: str, refplane: str = "ecliptic"
+    object_ids: List[str], time: Timestamp, location: str, refplane: str = "ecliptic"
 ) -> pd.DataFrame:
     """
     Get orbital elements as Cartesian, Cometary, and Keplerian representations from JPL Horizons.
@@ -26,7 +26,7 @@ def _get_orbital_elements(
     ----------
     object_id : List[str]
         Object IDs to
-    epoch : `~astropy.time.core.Time`
+    epoch : Timestamp
         Epoch at which to query orbital elements.
     location : str
         Location of the observer (in this case typically "@sun" or "@ssb"
@@ -78,12 +78,12 @@ def _get_orbital_elements(
     horizons_elements = vectors_df.merge(elements_df, on=["targetname", "datetime_jd"])
 
     horizons_elements.insert(
-        1, "mjd_tdb", Time(horizons_elements["datetime_jd"], format="jd").mjd
+        1, "mjd_tdb", Timestamp.from_jd(horizons_elements["datetime_jd"]).mjd()
     )
     horizons_elements.insert(
         len(horizons_elements.columns),
         "tp_mjd",
-        Time(horizons_elements["Tp_jd"], format="jd").mjd,
+        Timestamp.from_jd(horizons_elements["Tp_jd"]).mjd(),
     )
     horizons_elements.drop(columns=["datetime_jd", "Tp_jd"], inplace=True)
     return horizons_elements
@@ -134,7 +134,7 @@ def _get_orbital_elements(
 
                 horizons_elements_df = _get_orbital_elements(
                     [prov_designation],
-                    orbit.coordinates.time.to_astropy(),
+                    orbit.coordinates.time,
                     location=origin,
                     refplane=refplane,
                 )
diff --git a/adam_core/utils/spice.py b/adam_core/utils/spice.py
index 3bad0489..4c0ecf1b 100644
--- a/adam_core/utils/spice.py
+++ b/adam_core/utils/spice.py
@@ -31,7 +31,7 @@
 
 def _jd_tdb_to_et(jd_tdb: np.ndarray) -> np.ndarray:
     """
-    Convert an astropy Time object to an ephemeris time (ET) in seconds.
+    Convert TDB-scaled JD times to an ephemeris time (ET) in seconds.
 
     Parameters
     ----------