From 4249239a9d038db62a80a168c71f97862c208eb8 Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Fri, 26 Dec 2025 11:28:36 -0500
Subject: [PATCH 01/34] New 48P syndynes example, and explain how to use
 ephemeris from Horizons.

---
 docs/sbpy/dynamics.rst | 209 +++++++++++++++++++++++------------------
 1 file changed, 116 insertions(+), 93 deletions(-)

diff --git a/docs/sbpy/dynamics.rst b/docs/sbpy/dynamics.rst
index 60549fc2..a95253b0 100644
--- a/docs/sbpy/dynamics.rst
+++ b/docs/sbpy/dynamics.rst
@@ -293,6 +293,24 @@ To support generating synchrones at specific times, the :func:`~sbpy.dynamics.sy
    <Time object: scale='utc' format='iso' value=2023-10-01 00:00:00.000>
 
 
+Using Horizons
+--------------
+
+`Ephem.from_horizons` returns equatorial coordinates in the ICRF reference frame, which has its origin at the Solar System barycenter.  For `State` to correctly convert the ephemeris object to vectors, we need to set the Horizons observer to the Solar System barycenter (``"@ssb"``).  However, dynamical integrations are done in a heliocentric reference frame, so we transform the result to ``"heliocentriceclipticiau76"``:
+
+.. doctest-requires:: astroquery
+
+   >>> eph = Ephem.from_horizons(
+   ...     "48P",
+   ...     id_type="designation",
+   ...     closest_apparition=True,
+   ...     epochs=obstime,
+   ...     location="@ssb",
+   ... )
+   >>> comet = State.from_ephem(eph, frame="icrs")
+   >>> comet = comet.transform_to("heliocentriceclipticiau76")
+
+
 Adding ejection velocity
 ------------------------
 
@@ -388,110 +406,115 @@ In this example, we compute the syndynes of a comet orbiting β Pic (1.8 solar m
 Plotting syndynes and synchrones
 --------------------------------
 
-Generally, we are interested in plotting syndynes and synchrones on an image of a comet.  The accuracy of the coordinates object depends on the the comet and observer states, but also on whether or not light travel time is accounted for, and the accuracy of the orbit integrator.  The `sbpy` testing suite shows that arcsecond-level accuracy is possible, but this is generally not accurate enough for direct comparison to typical images of comets.  Instead, it helps to compute the positions of the syndynes and synchrone coordinate objects relative to the comet, and plot the results.
+Generally, we are interested in plotting syndynes and synchrones on an image of a comet.  The accuracy of the coordinates object depends on the the comet and observer states, but also on whether or not light travel time is accounted for, and the accuracy of the orbit integrator.  The `sbpy` testing suite shows that arcsecond-level accuracy is possible, but this is generally not accurate enough for direct comparison to typical images of comets, which need sub-arcsecond alignment.  Instead, it helps plot results relative to the comet's calculated position.
+
+The following example compares syndynes to a Spitzer Space Telesocpe image of comet 48P/Johnson (`Reach et al. 2007 <https://scixplorer.org/abs/2007Icar..191..298R/abstract>`_).  The FITS world coordinate system is used to account for the image orientation and scale.  To precisely align the syndynes with the comet nucleus, we update the world coordinate system to use our calculated comet coordinates:
 
-.. doctest-requires:: scipy,matplotlib
+.. doctest-requires:: scipy,astroquery
 
+   >>> from astropy.io import fits
+   >>> from astropy.wcs import WCS
+   >>>
+   >>> image, header = fits.getdata("48p-spitzer-reach07.fits", header=True)
+   >>> obstime = Time(header["DATE_OBS"])
+   >>>
+   >>> # get the comet state
+   >>> eph = Ephem.from_horizons(
+   ...     "48P",
+   ...     id_type="designation",
+   ...     closest_apparition=True,
+   ...     epochs=obstime,
+   ...     location="@ssb",
+   ... )
+   >>> comet = State.from_ephem(eph, frame="icrs")
+   >>> comet = comet.transform_to("heliocentriceclipticiau76")
+   >>> 
+   >>> # get the Spitzer Space Telescope state
+   >>> eph = Ephem.from_horizons("-79", id_type=None, epochs=obstime, location="@ssb")
+   >>> observer = State.from_ephem(eph, frame="icrs")
+   >>> 
+   >>> # set up the world coordinate system object and update the origin to align with
+   >>> # the calculated position of the comet
+   >>> wcs = WCS(header)
+   >>> coords0 = observer.observe(comet)[0].unmasked
+   >>> wcs.wcs.crval = coords0.ra.deg, coords0.dec.deg
+   >>> wcs.wcs.crpix = 209, 99
+
+
+.. plot::
+
+   >>> import numpy as np
    >>> import matplotlib.pyplot as plt
    >>> 
-   >>> coords0 = observer.observe(comet)
-   >>> def plot(ax, coords, **kwargs):
-   ...     dRA = coords.ra - coords0.ra
-   ...     dDec = coords.dec - coords0.dec
-   ...     ax.plot(dRA.arcsec, dDec.arcsec, **kwargs)
+   >>> import astropy.units as u
+   >>> from astropy.time import Time
+   >>> from astropy.io import fits
+   >>> from astropy.wcs import WCS
    >>> 
-   >>> fig, ax = plt.subplots()
+   >>> from sbpy.dynamics import State
+   >>> from sbpy.data import Ephem
+   >>> from sbpy.dynamics import SynGenerator
    >>> 
-   >>> # plot all but the last (beta=0) syndyne
-   >>> for syndyne in dust.syndynes()[:-1]:
-   ...     plot(ax, syndyne.coords, label=f"$\\beta={syndyne.beta:.2g}$")
+   >>> image, header = fits.getdata("48p-spitzer-reach07.fits", header=True)
+   >>> obstime = Time(header["DATE_OBS"])
    >>> 
-   >>> # plot every 5th synchrone
-   >>> for synchrone in dust.synchrones()[4::5]:
-   ...     plot(
-   ...         ax,
-   ...         synchrone.coords,
-   ...         ls="--",
-   ...         label=f"$\\Delta t={synchrone.age.to(u.d):.2g}$"
-   ...     )
+   >>> # Ephem.from_horizons returns equatorial coordinates in the ICRF reference
+   >>> # frame, which has its origin at the Solar System barycenter.  For State to
+   >>> # correctly convert the ephemeris to vectors, we need to set the Horizons
+   >>> # observer to the Solar System barycenter: @ssb
    >>> 
-   >>> # and plot the orbit
-   >>> dt = np.linspace(-2, 2) * u.d
+   >>> # get the position of the comet and transform to a heliocentric frame for
+   >>> # integration
+   >>> eph = Ephem.from_horizons(
+   ...     "48P",
+   ...     id_type="designation",
+   ...     closest_apparition=True,
+   ...     epochs=obstime,
+   ...     location="@ssb",
+   ... )
+   >>> comet = State.from_ephem(eph, frame="icrs")
+   >>> comet = comet.transform_to("heliocentriceclipticiau76")
+   >>> 
+   >>> # get the position of the Spitzer Space Telescope
+   >>> eph = Ephem.from_horizons("-79", id_type=None, epochs=obstime, location="@ssb")
+   >>> observer = State.from_ephem(eph, frame="icrs")
+   >>> 
+   >>> # generate the syndynes
+   >>> betas = [1, 0.1, 0.01, 0.001]
+   >>> ages = np.linspace(0, 365, 51) * u.day
+   >>> dust = SynGenerator(comet[0], betas, ages, observer=observer[0])
+   >>> 
+   >>> # Set up the world coordinate system object and update the origin to align with
+   >>> # the calculated position of the comet.
+   >>> wcs = WCS(header)
+   >>> coords0 = observer.observe(comet)[0].unmasked
+   >>> wcs.wcs.crval = coords0.ra.deg, coords0.dec.deg
+   >>> wcs.wcs.crpix = 209, 99
+   >>> 
+   >>> # plot
+   >>> fig, ax = plt.subplots(num=1, clear=True, figsize=(6.5, 3.25))
+   >>> 
+   >>> ax.imshow(image, origin="lower", vmin=14, vmax=70, cmap="gray_r")
+   >>> 
+   >>> # save xlim and ylim for later
+   >>> xlim = ax.get_xlim()
+   >>> ylim = ax.get_ylim()
+   >>> 
+   >>> # plot the syndynes
+   >>> for syndyne in dust.syndynes():
+   ...     x, y = wcs.world_to_pixel(syndyne.coords.unmasked)
+   ...     ax.plot(x, y, label=f"$\\beta={syndyne.beta:.2g}$")
+   >>> 
+   >>> # plot orbit
+   >>> dt = np.linspace(-1, 1) * 3 * u.d
    >>> states, coords = dust.source_orbit(dt)
-   >>> plot(ax, coords, color="k", ls=":", label="Orbit")
-   >>>
-   >>> ax.invert_xaxis()
-   >>> plt.setp(ax,
-   ...          xlim=[100, -10],
-   ...          ylim=[-10, 100],
-   ...          xlabel="$\\Delta$RA (arcsec)",
-   ...          ylabel="$\\Delta$Dec (arcsec)",
-   ... ) # doctest: +SKIP
-   >>> plt.legend() # doctest: +SKIP
+   >>> x, y = wcs.world_to_pixel(coords.unmasked)
+   >>> ax.plot(x, y, color="tab:olive", ls="--", label="Orbit")
+   >>> 
+   >>> plt.setp(ax, xlim=xlim, ylim=ylim)
+   >>> plt.legend()
    >>> plt.tight_layout()
 
-.. plot::
-
-   import numpy as np
-   import matplotlib.pyplot as plt
-
-   import astropy.units as u
-   from astropy.time import Time
-   from sbpy.dynamics import State, SynGenerator
-
-   r = [2, 0, 0] * u.au
-   v = [0, 30, 0] * u.km / u.s
-   t = Time("2023-12-08")
-   frame = "heliocentriceclipticiau76"
-   comet = State(r, v, t, frame=frame)
-
-   betas = [1, 0.1, 0.01, 0]
-   ages = np.linspace(0, 100, 25) * u.day
-   observer = State(
-       r=[0, 2, 2] * u.au,
-       v=[0, 0, 0] * u.km / u.s,
-       t=comet.t,
-       frame="icrs",
-   )
-   dust = SynGenerator(comet, betas, ages, observer=observer)
-
-   coords0 = observer.observe(comet)
-   def plot(ax, coords, **kwargs):
-       dRA = coords.ra - coords0.ra
-       dDec = coords.dec - coords0.dec
-       ax.plot(dRA.arcsec, dDec.arcsec, **kwargs)
-   
-   fig, ax = plt.subplots()
-   
-   for syndyne in dust.syndynes():
-       # don't draw the beta = 0 syndyne
-       if syndyne.beta == 0:
-           continue
-       plot(ax, syndyne.coords, label=f"$\\beta={syndyne.beta:.2g}$")
-   
-   # plot every 5th synchrone
-   for synchrone in dust.synchrones()[4::5]:
-       plot(
-           ax,
-           synchrone.coords,
-           ls="--",
-           label=f"$\\Delta t={synchrone.age.to(u.d):.2g}$",
-       )
-   
-   # and plot the orbit
-   dt = np.linspace(-2, 2) * u.d
-   states, coords = dust.source_orbit(dt)
-   plot(ax, coords, color="k", ls=":", label="Orbit")
-
-   ax.invert_xaxis()
-   plt.setp(ax,
-            xlim=[100, -10],
-            ylim=[-10, 100],
-            xlabel="$\\Delta$RA (arcsec)",
-            ylabel="$\\Delta$Dec (arcsec)",
-   )
-   plt.legend()
-   plt.tight_layout()
 
 
 Reference/API

From 3046647483d58c6c0f1dfd3e1eed9c52167b4eb9 Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Sat, 27 Dec 2025 13:55:51 -0500
Subject: [PATCH 02/34] Update file location.

---
 docs/sbpy/dynamics.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/sbpy/dynamics.rst b/docs/sbpy/dynamics.rst
index a95253b0..ac1da032 100644
--- a/docs/sbpy/dynamics.rst
+++ b/docs/sbpy/dynamics.rst
@@ -415,7 +415,7 @@ The following example compares syndynes to a Spitzer Space Telesocpe image of co
    >>> from astropy.io import fits
    >>> from astropy.wcs import WCS
    >>>
-   >>> image, header = fits.getdata("48p-spitzer-reach07.fits", header=True)
+   >>> image, header = fits.getdata("https://sbpy.org/data/48p-spitzer-reach07.fits", header=True)
    >>> obstime = Time(header["DATE_OBS"])
    >>>
    >>> # get the comet state

From 9fcf5ae2f2337f5fc1bc332d8cfb0fa80d2a2444 Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Sat, 27 Dec 2025 17:23:00 -0500
Subject: [PATCH 03/34] New unmasked() function.

---
 sbpy/utils/core.py            | 25 ++++++++++++++++++++++++-
 sbpy/utils/tests/test_core.py | 15 ++++++++++++++-
 2 files changed, 38 insertions(+), 2 deletions(-)

diff --git a/sbpy/utils/core.py b/sbpy/utils/core.py
index caec9205..b0ab6a44 100644
--- a/sbpy/utils/core.py
+++ b/sbpy/utils/core.py
@@ -6,10 +6,20 @@
 
 """
 
-__all__ = ["required_packages", "optional_packages"]
+__all__ = ["required_packages", "optional_packages", "unmasked"]
 
 from importlib import import_module
 from warnings import warn
+from packaging.version import Version
+
+import numpy as np
+import astropy
+
+if Version(astropy.__version__) >= Version("7.0"):
+    from astropy.utils.masked import get_data_and_mask
+else:
+    get_data_and_mask = None
+
 from ..exceptions import RequiredPackageUnavailable, OptionalPackageUnavailable
 
 
@@ -100,3 +110,16 @@ def optional_packages(*packages, message=None):
             )
             return False
     return True
+
+
+def unmasked(array):
+    """Return an unmasked version of the array."""
+
+    if get_data_and_mask is None:
+        if hasattr(array, "unmasked"):
+            return array.unmasked
+        elif isinstance(array, np.ma.MaskedArray):
+            return array._data
+        return array
+    else:
+        return get_data_and_mask(array)[0]
diff --git a/sbpy/utils/tests/test_core.py b/sbpy/utils/tests/test_core.py
index 272b0cad..2020dab9 100644
--- a/sbpy/utils/tests/test_core.py
+++ b/sbpy/utils/tests/test_core.py
@@ -1,5 +1,7 @@
 import pytest
-from ..core import required_packages, optional_packages
+import numpy as np
+from astropy.utils.masked import Masked
+from ..core import required_packages, optional_packages, unmasked
 from ...exceptions import RequiredPackageUnavailable, OptionalPackageUnavailable
 
 
@@ -27,3 +29,14 @@ def test_optional_message():
     message = "Using linear interpolation."
     with pytest.warns(OptionalPackageUnavailable, match=message) as record:
         optional_packages("unavailable_package", message=message)
+
+
+def test_unmasked():
+    array = [1, 2, 3]
+    mask = [False, False, True]
+    ma1 = Masked([1.0, 2.0, 3.0], mask=[False, False, True])
+    ma2 = np.ma.MaskedArray(array, mask=mask)
+
+    assert unmasked(array) is array
+    assert (unmasked(ma1) == array).all()
+    assert (unmasked(ma2) == array).all()

From 83f942755182cba5dbd1477684b63fcb831668fa Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Sat, 27 Dec 2025 20:09:25 -0500
Subject: [PATCH 04/34] New syndynes example, plot methods.

---
 CHANGES.rst                          |  22 +++
 docs/sbpy/dynamics.rst               | 254 ++++++++++++++++++++-------
 sbpy/dynamics/syndynes.py            | 171 +++++++++++++++++-
 sbpy/dynamics/tests/test_syndynes.py |  12 ++
 4 files changed, 391 insertions(+), 68 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index f5bc53d4..c6a1722b 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -8,12 +8,34 @@ Update supported versions [#427]:
 New Features
 ------------
 
+sbpy.dynamics
+^^^^^^^^^^^^^
+
+- Plot syndynes and synchrones with the new ``plot()`` methods on
+  `sbpy.dynamics.syndynes.Syndyne`, `sbpy.dynamics.syndynes.Syndynes`,
+  `sbpy.dynamics.syndynes.Synchrone`, and `sbpy.dynamics.syndynes.Synchrones`.
+  [#434]
+
 sbpy.photometry
 ^^^^^^^^^^^^^^^
 
 - Added Rubin Observatory's LSSTCam filter set to `bandpass()`. [#431]
 
 
+API Changes
+-----------
+
+sbpy.dynamics
+^^^^^^^^^^^^^
+
+- Indexing `sbpy.dynamics.syndynes.Syndynes` or
+  `sbpy.dynamics.syndynes.Synchrones` with an integer, e.g., ``syndynes[0]``,
+  will return a single `sbpy.dynamics.syndynes.Syndyne` or
+  `sbpy.dynamics.syndynes.Synchrone`.  Indexing with a tuple or slice, e.g.,
+  ``syndynes[:2]`` will return `sbpy.dynamics.syndynes.Syndynes` or
+  `sbpy.dynamics.syndynes.Synchrones`.  Previously a list was returned.
+
+
 v0.6.0 (2025-12-02)
 ===================
 
diff --git a/docs/sbpy/dynamics.rst b/docs/sbpy/dynamics.rst
index ac1da032..0acb2ee8 100644
--- a/docs/sbpy/dynamics.rst
+++ b/docs/sbpy/dynamics.rst
@@ -406,11 +406,122 @@ In this example, we compute the syndynes of a comet orbiting β Pic (1.8 solar m
 Plotting syndynes and synchrones
 --------------------------------
 
-Generally, we are interested in plotting syndynes and synchrones on an image of a comet.  The accuracy of the coordinates object depends on the the comet and observer states, but also on whether or not light travel time is accounted for, and the accuracy of the orbit integrator.  The `sbpy` testing suite shows that arcsecond-level accuracy is possible, but this is generally not accurate enough for direct comparison to typical images of comets, which need sub-arcsecond alignment.  Instead, it helps plot results relative to the comet's calculated position.
+Generally, we are interested in visualizing syndynes and synchrones for an observer.  We can plot the coordinates themselves, or plot them relative to the comet.
+
+.. note::
+   The accuracy of the coordinates object depends on the the comet and observer states, but also on whether or not light travel time is accounted for, and the accuracy of the orbit integrator.  The `sbpy` testing suite shows that arcsecond-level accuracy is possible, but this is generally not enough for direct comparison to typical images of comets, which need sub-arcsecond alignment.  For small offsets, plotting results relative to the comet's calculated position is sufficient.
+
+Here is a simple example that plots the syndynes and synchrones as offsets from the comet's coordinates:
+
+.. doctest-requires:: scipy,matplotlib
+
+   >>> import matplotlib.pyplot as plt
+   >>>
+   >>> # comet's caluclated coordinates
+   >>> coords0 = observer.observe(comet)
+   >>>
+   >>> def plot(ax, coords, **kwargs):
+   ...     "plot coordinates relative to comet"
+   ...     dRA = coords.ra - coords0.ra
+   ...     dDec = coords.dec - coords0.dec
+   ...     ax.plot(dRA.arcsec, dDec.arcsec, **kwargs)
+   >>>
+   >>> fig, ax = plt.subplots()
+   >>>
+   >>> # plot all but the last (beta=0) syndyne
+   >>> for syndyne in dust.syndynes()[:-1]:
+   ...     plot(ax, syndyne.coords, label=f"$\\beta={syndyne.beta:.2g}$")
+   >>>
+   >>> # plot every 5th synchrone
+   >>> for synchrone in dust.synchrones()[4::5]:
+   ...     plot(
+   ...         ax,
+   ...         synchrone.coords,
+   ...         ls="--",
+   ...         label=f"$\\Delta t={synchrone.age.to(u.d):.2g}$"
+   ...     )
+   >>>
+   >>> # and plot the orbit
+   >>> dt = np.linspace(-2, 2) * u.d
+   >>> states, coords = dust.source_orbit(dt)
+   >>> plot(ax, coords, color="k", ls=":", label="Orbit")
+   >>>
+   >>> ax.invert_xaxis()
+   >>> plt.setp(ax,
+   ...          xlim=[100, -10],
+   ...          ylim=[-10, 100],
+   ...          xlabel="$\\Delta$RA (arcsec)",
+   ...          ylabel="$\\Delta$Dec (arcsec)",
+   ... ) # doctest: +SKIP
+   >>> plt.legend() # doctest: +SKIP
+   >>> plt.tight_layout()
+
+.. plot::
+
+   import numpy as np
+   import matplotlib.pyplot as plt
+
+   import astropy.units as u
+   from astropy.time import Time
+   from sbpy.dynamics import State, SynGenerator
+
+   r = [2, 0, 0] * u.au
+   v = [0, 30, 0] * u.km / u.s
+   t = Time("2023-12-08")
+   frame = "heliocentriceclipticiau76"
+   comet = State(r, v, t, frame=frame)
+
+   betas = [1, 0.1, 0.01, 0]
+   ages = np.linspace(0, 100, 25) * u.day
+   observer = State(
+       r=[0, 2, 2] * u.au,
+       v=[0, 0, 0] * u.km / u.s,
+       t=comet.t,
+       frame="icrs",
+   )
+   dust = SynGenerator(comet, betas, ages, observer=observer)
+
+   coords0 = observer.observe(comet)
+   def plot(ax, coords, **kwargs):
+       dRA = coords.ra - coords0.ra
+       dDec = coords.dec - coords0.dec
+       ax.plot(dRA.arcsec, dDec.arcsec, **kwargs)
+
+   fig, ax = plt.subplots()
+
+   for syndyne in dust.syndynes():
+       # don't draw the beta = 0 syndyne
+       if syndyne.beta == 0:
+           continue
+       plot(ax, syndyne.coords, label=f"$\\beta={syndyne.beta:.2g}$")
+
+   # plot every 5th synchrone
+   for synchrone in dust.synchrones()[4::5]:
+       plot(
+           ax,
+           synchrone.coords,
+           ls="--",
+           label=f"$\Delta t={synchrone.age.to(u.d):.2g}$",
+       )
+
+   # and plot the orbit
+   dt = np.linspace(-2, 2) * u.d
+   states, coords = dust.source_orbit(dt)
+   plot(ax, coords, color="k", ls=":", label="Orbit")
+
+   ax.invert_xaxis()
+   plt.setp(ax,
+            xlim=[100, -10],
+            ylim=[-10, 100],
+            xlabel="$\Delta$RA (arcsec)",
+            ylabel="$\Delta$Dec (arcsec)",
+   )
+   plt.legend()
+   plt.tight_layout()
 
 The following example compares syndynes to a Spitzer Space Telesocpe image of comet 48P/Johnson (`Reach et al. 2007 <https://scixplorer.org/abs/2007Icar..191..298R/abstract>`_).  The FITS world coordinate system is used to account for the image orientation and scale.  To precisely align the syndynes with the comet nucleus, we update the world coordinate system to use our calculated comet coordinates:
 
-.. doctest-requires:: scipy,astroquery
+.. doctest-requires:: scipy,astroquery,matplotlib
 
    >>> from astropy.io import fits
    >>> from astropy.wcs import WCS
@@ -439,81 +550,96 @@ The following example compares syndynes to a Spitzer Space Telesocpe image of co
    >>> coords0 = observer.observe(comet)[0].unmasked
    >>> wcs.wcs.crval = coords0.ra.deg, coords0.dec.deg
    >>> wcs.wcs.crpix = 209, 99
-
-
-.. plot::
-
-   >>> import numpy as np
-   >>> import matplotlib.pyplot as plt
-   >>> 
-   >>> import astropy.units as u
-   >>> from astropy.time import Time
-   >>> from astropy.io import fits
-   >>> from astropy.wcs import WCS
-   >>> 
-   >>> from sbpy.dynamics import State
-   >>> from sbpy.data import Ephem
-   >>> from sbpy.dynamics import SynGenerator
-   >>> 
-   >>> image, header = fits.getdata("48p-spitzer-reach07.fits", header=True)
-   >>> obstime = Time(header["DATE_OBS"])
-   >>> 
-   >>> # Ephem.from_horizons returns equatorial coordinates in the ICRF reference
-   >>> # frame, which has its origin at the Solar System barycenter.  For State to
-   >>> # correctly convert the ephemeris to vectors, we need to set the Horizons
-   >>> # observer to the Solar System barycenter: @ssb
-   >>> 
-   >>> # get the position of the comet and transform to a heliocentric frame for
-   >>> # integration
-   >>> eph = Ephem.from_horizons(
-   ...     "48P",
-   ...     id_type="designation",
-   ...     closest_apparition=True,
-   ...     epochs=obstime,
-   ...     location="@ssb",
-   ... )
-   >>> comet = State.from_ephem(eph, frame="icrs")
-   >>> comet = comet.transform_to("heliocentriceclipticiau76")
-   >>> 
-   >>> # get the position of the Spitzer Space Telescope
-   >>> eph = Ephem.from_horizons("-79", id_type=None, epochs=obstime, location="@ssb")
-   >>> observer = State.from_ephem(eph, frame="icrs")
-   >>> 
+   >>>
    >>> # generate the syndynes
    >>> betas = [1, 0.1, 0.01, 0.001]
    >>> ages = np.linspace(0, 365, 51) * u.day
    >>> dust = SynGenerator(comet[0], betas, ages, observer=observer[0])
    >>> 
-   >>> # Set up the world coordinate system object and update the origin to align with
-   >>> # the calculated position of the comet.
-   >>> wcs = WCS(header)
-   >>> coords0 = observer.observe(comet)[0].unmasked
-   >>> wcs.wcs.crval = coords0.ra.deg, coords0.dec.deg
-   >>> wcs.wcs.crpix = 209, 99
-   >>> 
-   >>> # plot
+   >>> # plot the image and syndynes
    >>> fig, ax = plt.subplots(num=1, clear=True, figsize=(6.5, 3.25))
    >>> 
    >>> ax.imshow(image, origin="lower", vmin=14, vmax=70, cmap="gray_r")
    >>> 
-   >>> # save xlim and ylim for later
-   >>> xlim = ax.get_xlim()
-   >>> ylim = ax.get_ylim()
-   >>> 
    >>> # plot the syndynes
    >>> for syndyne in dust.syndynes():
    ...     x, y = wcs.world_to_pixel(syndyne.coords.unmasked)
    ...     ax.plot(x, y, label=f"$\\beta={syndyne.beta:.2g}$")
-   >>> 
-   >>> # plot orbit
-   >>> dt = np.linspace(-1, 1) * 3 * u.d
-   >>> states, coords = dust.source_orbit(dt)
-   >>> x, y = wcs.world_to_pixel(coords.unmasked)
-   >>> ax.plot(x, y, color="tab:olive", ls="--", label="Orbit")
-   >>> 
-   >>> plt.setp(ax, xlim=xlim, ylim=ylim)
-   >>> plt.legend()
-   >>> plt.tight_layout()
+
+.. plot::
+   :show-source-link:
+
+   import numpy as np
+   import matplotlib.pyplot as plt
+   
+   import astropy.units as u
+   from astropy.time import Time
+   from astropy.io import fits
+   from astropy.wcs import WCS
+   
+   from sbpy.dynamics import State
+   from sbpy.data import Ephem
+   from sbpy.dynamics import SynGenerator
+   
+   image, header = fits.getdata("https://sbpy.org/data/48p-spitzer-reach07.fits", header=True)
+   obstime = Time(header["DATE_OBS"])
+   
+   # Ephem.from_horizons returns equatorial coordinates in the ICRF reference
+   # frame, which has its origin at the Solar System barycenter.  For State to
+   # correctly convert the ephemeris to vectors, we need to set the Horizons
+   # observer to the Solar System barycenter: @ssb
+   
+   # get the position of the comet and transform to a heliocentric frame for
+   # integration
+   eph = Ephem.from_horizons(
+       "48P",
+       id_type="designation",
+       closest_apparition=True,
+       epochs=obstime,
+       location="@ssb",
+   )
+   comet = State.from_ephem(eph, frame="icrs")
+   comet = comet.transform_to("heliocentriceclipticiau76")
+   
+   # get the position of the Spitzer Space Telescope
+   eph = Ephem.from_horizons("-79", id_type=None, epochs=obstime, location="@ssb")
+   observer = State.from_ephem(eph, frame="icrs")
+   
+   # generate the syndynes
+   betas = [1, 0.1, 0.01, 0.001]
+   ages = np.linspace(0, 365, 51) * u.day
+   dust = SynGenerator(comet[0], betas, ages, observer=observer[0])
+   
+   # Set up the world coordinate system object and update the origin to align with
+   # the calculated position of the comet.
+   wcs = WCS(header)
+   coords0 = observer.observe(comet)[0].unmasked
+   wcs.wcs.crval = coords0.ra.deg, coords0.dec.deg
+   wcs.wcs.crpix = 209, 99
+   
+   # plot
+   fig, ax = plt.subplots(num=1, clear=True, figsize=(6.5, 3.25))
+   
+   ax.imshow(image, origin="lower", vmin=14, vmax=70, cmap="gray_r")
+   
+   # save xlim and ylim for later
+   xlim = ax.get_xlim()
+   ylim = ax.get_ylim()
+   
+   # plot the syndynes
+   for syndyne in dust.syndynes():
+       x, y = wcs.world_to_pixel(syndyne.coords.unmasked)
+       ax.plot(x, y, label=f"$\\beta={syndyne.beta:.2g}$")
+   
+   # plot orbit
+   dt = np.linspace(-1, 1) * 3 * u.d
+   states, coords = dust.source_orbit(dt)
+   x, y = wcs.world_to_pixel(coords.unmasked)
+   ax.plot(x, y, color="tab:olive", ls="--", label="Orbit")
+   
+   plt.setp(ax, xlim=xlim, ylim=ylim)
+   plt.legend()
+   plt.tight_layout()
 
 
 
diff --git a/sbpy/dynamics/syndynes.py b/sbpy/dynamics/syndynes.py
index a198fbcb..cbbeb3c7 100644
--- a/sbpy/dynamics/syndynes.py
+++ b/sbpy/dynamics/syndynes.py
@@ -22,17 +22,27 @@
 import logging
 from typing import Iterable, Optional, TypeVar
 
+try:
+    import matplotlib as mpl
+    import matplotlib.pyplot as plt
+except ImportError:
+    pass
+
 import numpy as np
 import astropy.units as u
+from astropy.wcs import WCS
 from astropy.time import Time
 from astropy.table import vstack
 from astropy.coordinates import SkyCoord
 
 from ..data import Ephem
+from ..utils.decorators import requires
+from ..utils import unmasked
 from .models import DynamicalModel, SolarGravityAndRadiationPressure
-from .state import StateBase, StateBaseType, State
+from .state import StateBase, State
 
 SynGeneratorType = TypeVar("SynGeneratorType", bound="SynGenerator")
+SynCollectionType = TypeVar("SynCollectionType", bound="SynCollection")
 
 
 class SynStates(StateBase, abc.ABC):
@@ -172,6 +182,56 @@ def to_ephem(self) -> Ephem:
 
         return Ephem.from_dict(data, meta=meta)
 
+    @requires("matplotlib")
+    def plot(
+        self,
+        ax: mpl.axes.Axes | None = None,
+        *,
+        wcs: WCS | None = None,
+        unit: u.Unit | str = "arcsec",
+        **kwargs,
+    ) -> None:
+        """Plot the coordinates.
+
+        Requires `self.observer`.
+
+
+        Parameters
+        ----------
+        ax : `~matplotlib.axes.Axes`, optional
+            Plot to this axis object.
+
+        wcs : `~astropy.wcs.WCS`, optional
+            Transform coordinates to the image plane using this world coordinate
+            system object.
+
+        unit : `~astropy.unit.Unit` or str, optional
+            Plot in these angular units.
+
+        **kwargs : dict
+            Keyword arguments are passed along to ``ax.plot``.
+
+        """
+
+        if wcs is None:
+            # plot offsets from the source in the tangent plane
+
+            if self.observer is None:
+                raise ValueError("observer is not defined")
+
+            coords0 = self.observer.observe(self.source)
+
+            dRA = self.coords.ra - coords0.ra
+            dDec = self.coords.dec - coords0.dec
+            x = dRA.to(unit) / np.cos(self.coords.dec)
+            y = dDec.to(unit)
+        else:
+            # convert coordinates to plot units with the WCS object (avoid
+            # passing a masked object, or wcs will complain)
+            x, y = wcs.world_to_pixel(unmasked(self.coords))
+
+        ax.plot(x, y, **kwargs)
+
 
 class Syndyne(SynStates):
     """Collection of particle states that make up a syndyne.
@@ -179,7 +239,7 @@ class Syndyne(SynStates):
 
     Parameters
     ----------
-    source : State
+    source : `State`
         The source of the syndyne dust.
 
     beta : float
@@ -197,6 +257,9 @@ class Syndyne(SynStates):
     t : `~astropy.time.Time` or `~astropy.units.Quantity`
         Time of observation.
 
+    initial : `State`
+        The initial states for each particle.
+
     observer : `~sbpy.dynamics.State`, optional
         The observer, used to generate sky coordinates.
 
@@ -245,6 +308,9 @@ class Synchrone(SynStates):
     t : `~astropy.time.Time` or `~astropy.units.Quantity`
         Time of observation.
 
+    initial : `State`
+        The initial states for each particle.
+
     observer : `~sbpy.dynamics.State`, optional
         The observer, used to generate sky coordinates.
 
@@ -274,6 +340,48 @@ def epoch(self) -> Time:
         return self.source.t - self.ages[0]
 
 
+class SourceOrbit(SynStates):
+    """Collection of states that make up an orbit.
+
+    Parameters
+    ----------
+    source : State
+        The orbit is for this source state.
+
+    ages : ~astropy.units.Quantity
+        Array of ages, shape = (N,).
+
+    r : `~astropy.units.Quantity`
+        Position (x, y, z), shape = (N, 3).  Same coordinate frame as ``source``.
+
+    v : `~astropy.units.Quantity`
+        Velocity (x, y, z), shape = (N, 3).  Same coordinate frame as ``source``.
+
+    t : `~astropy.time.Time` or `~astropy.units.Quantity`
+        Time of observation.
+
+    observer : `~sbpy.dynamics.State`, optional
+        The observer, used to generate sky coordinates.
+
+    """
+
+    def __init__(
+        self,
+        source: State,
+        ages: u.Quantity,
+        r: u.Quantity,
+        v: u.Quantity,
+        t: Time,
+        observer: Optional[State] = None,
+    ) -> None:
+        super().__init__(source, 0, ages, r, v, t, source, observer=observer)
+
+    @property
+    def epoch(self) -> Time:
+        """Epoch of each orbit point."""
+        return self.source.t - self.ages[0]
+
+
 class SynCollection:
     """Immutable collection of syndynes or synchrones.
 
@@ -310,8 +418,13 @@ def __len__(self) -> int:
         """Number of items in the container."""
         return len(self._data)
 
-    def __getitem__(self, k: int | tuple | slice) -> SynStates:
-        return self._data[k]
+    def __getitem__(self, k: int | tuple | slice) -> SynStates | SynCollectionType:
+        # return a SynStates object
+        if isinstance(k, int):
+            return self._data[k]
+
+        # return a SynCollection
+        return type(self)(self._data[k])
 
     def to_ephem(self) -> Ephem:
         """Convert to an sbpy ephemeris object.
@@ -330,6 +443,56 @@ def to_ephem(self) -> Ephem:
             meta=tables[0].meta,
         )
 
+    @requires("matplotlib")
+    def plot(
+        self,
+        ax: mpl.axes.Axes | None = None,
+        *,
+        wcs: WCS | None = None,
+        unit: u.Unit | str = "arcsec",
+        label_format: str | None = None,
+        **kwargs,
+    ) -> None:
+        """Plot the collection.
+
+
+        Parameters
+        ----------
+        ax : `~matplotlib.axes.Axes`, optional
+            Plot to this axis object.
+
+        wcs : `~astropy.wcs.WCS`, optional
+            Transform coordinates to the image plane using this world coordinate
+            system object.
+
+        unit : `~astropy.unit.Unit` or str, optional
+            Plot in these angular units.
+
+        label_format : str, optional
+            A string defining the label format.  ``{beta}`` or ``{age}`` will be
+            replaced as appropriate for the syndyne or synchrone.
+
+        **kwargs : dict
+            Keyword arguments are passed along to ``ax.plot``.
+
+        """
+
+        ax = plt.gca() if ax is None else ax
+
+        match (label_format, self):
+            case (None, Syndynes()):
+                label_format = "$\\beta={beta}$"
+            case (None, Synchrones()):
+                label_format = "$\\Delta t={age}$"
+            case (None, _):
+                label_format = ""
+
+        for syn in self:
+            label = label_format.format(
+                beta=getattr(syn, "beta", None), age=getattr(syn, "age", None)
+            )
+            syn.plot(ax=ax, wcs=wcs, unit=unit, label=label, **kwargs)
+
 
 class Syndynes(SynCollection, data_type=Syndyne):
     """Collection of syndynes.
diff --git a/sbpy/dynamics/tests/test_syndynes.py b/sbpy/dynamics/tests/test_syndynes.py
index 8e239d48..005aac61 100644
--- a/sbpy/dynamics/tests/test_syndynes.py
+++ b/sbpy/dynamics/tests/test_syndynes.py
@@ -212,6 +212,12 @@ def test_init_1d(self, example_syndynes):
                 initial_states,
             )
 
+    def test_getitem(self, example_syndynes):
+        _, _, _, dust, _ = example_syndynes
+        assert isinstance(dust.syndynes(), Syndynes)
+        assert isinstance(dust.syndynes()[0], Syndyne)
+        assert isinstance(dust.syndynes()[:1], Syndynes)
+
 
 class TestSynchrone:
     def test_init(self, example_syndynes):
@@ -286,6 +292,12 @@ def test_to_ephem(self, example_syndynes):
         assert np.allclose(eph["beta_rad"], betas)
         assert u.allclose(eph["age"], age)
 
+    def test_getitem(self, example_syndynes):
+        _, _, _, dust, _ = example_syndynes
+        assert isinstance(dust.synchrones(), Synchrones)
+        assert isinstance(dust.synchrones()[0], Synchrone)
+        assert isinstance(dust.synchrones()[:1], Synchrones)
+
 
 def random_state(length=1):
     return State(

From 90af062409bef9e2c1376a07102f08a76a43a64a Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Sun, 28 Dec 2025 10:16:45 -0500
Subject: [PATCH 05/34] Add missing to_ephem

---
 sbpy/dynamics/state.py            | 24 ++++++++++++++++++++++++
 sbpy/dynamics/tests/test_state.py | 23 ++++++++++++++++++++++-
 2 files changed, 46 insertions(+), 1 deletion(-)

diff --git a/sbpy/dynamics/state.py b/sbpy/dynamics/state.py
index a41aad8a..bff66a2c 100644
--- a/sbpy/dynamics/state.py
+++ b/sbpy/dynamics/state.py
@@ -435,6 +435,30 @@ def from_ephem(
             " velocity fields."
         )
 
+    def to_ephem(self) -> Ephem:
+        """Convert to an ephemeris object.
+
+
+        Returns
+        -------
+        eph : `~sbpy.data.ephem.Ephem`
+
+        """
+
+        time_field = "dt" if self.arbitrary_time else "date"
+
+        eph = {
+            time_field: self.t,
+            "x": self.x,
+            "y": self.y,
+            "z": self.z,
+            "vx": self.v_x,
+            "vy": self.v_y,
+            "vz": self.v_z,
+        }
+
+        return Ephem.from_dict(eph, meta={"frame": self.frame})
+
     def observe(self, target: StateType) -> SkyCoord:
         """Project a target's position onto the sky.
 
diff --git a/sbpy/dynamics/tests/test_state.py b/sbpy/dynamics/tests/test_state.py
index 085c8647..7ae638de 100644
--- a/sbpy/dynamics/tests/test_state.py
+++ b/sbpy/dynamics/tests/test_state.py
@@ -5,7 +5,7 @@
 import numpy as np
 import astropy.units as u
 from astropy.time import Time
-from astropy.coordinates import HeliocentricEclipticIAU76, SkyCoord
+from astropy.coordinates import HeliocentricEclipticIAU76, SkyCoord, ICRS
 
 from ...data import Ephem
 from ..state import ArbitraryFrame, State
@@ -561,3 +561,24 @@ def test_from_ephem(self):
             del incomplete.table[k]
             with pytest.raises(ValueError):
                 State.from_ephem(incomplete)
+
+    def test_to_ephm(self):
+        r = np.arange(1, 7).reshape([2, 3]) * u.au
+        v = r.value / 10 * u.km / u.s
+        t = [0, 1] * u.day
+        state = State(r, v, t)
+        eph = state.to_ephem()
+
+        assert u.allclose(eph["x"], [1, 4] * u.au)
+        assert u.allclose(eph["y"], [2, 5] * u.au)
+        assert u.allclose(eph["z"], [3, 6] * u.au)
+        assert u.allclose(eph["vx"], [0.1, 0.4] * u.km / u.s)
+        assert u.allclose(eph["vy"], [0.2, 0.5] * u.km / u.s)
+        assert u.allclose(eph["vz"], [0.3, 0.6] * u.km / u.s)
+        assert np.allclose(eph["dt"], [0, 1] * u.day)
+        assert isinstance(eph.meta["frame"], ArbitraryFrame)
+
+        state = State(r, v, Time("2000-01-01") + t, frame="icrs")
+        eph = state.to_ephem()
+        assert np.allclose(eph["date"].mjd, [51544.0, 51545.0])
+        assert isinstance(eph.meta["frame"], ICRS)

From 2b0c298d65376519b642d41b9e71d13d6ecc5f49 Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Sun, 28 Dec 2025 10:24:34 -0500
Subject: [PATCH 06/34] Use relative time in SourceOrbit

---
 sbpy/dynamics/syndynes.py | 25 +++++++++++--------------
 1 file changed, 11 insertions(+), 14 deletions(-)

diff --git a/sbpy/dynamics/syndynes.py b/sbpy/dynamics/syndynes.py
index cbbeb3c7..7f0c012f 100644
--- a/sbpy/dynamics/syndynes.py
+++ b/sbpy/dynamics/syndynes.py
@@ -8,6 +8,7 @@
 """
 
 __all__ = [
+    "SourceOrbit",
     "SynGenerator",
     "SynStates",
     "SynCollection",
@@ -348,8 +349,8 @@ class SourceOrbit(SynStates):
     source : State
         The orbit is for this source state.
 
-    ages : ~astropy.units.Quantity
-        Array of ages, shape = (N,).
+    dt : ~astropy.units.Quantity
+        Each point's time relative to the observation time, shape = (N,).
 
     r : `~astropy.units.Quantity`
         Position (x, y, z), shape = (N, 3).  Same coordinate frame as ``source``.
@@ -368,13 +369,17 @@ class SourceOrbit(SynStates):
     def __init__(
         self,
         source: State,
-        ages: u.Quantity,
+        dt: u.Quantity,
         r: u.Quantity,
         v: u.Quantity,
         t: Time,
         observer: Optional[State] = None,
     ) -> None:
-        super().__init__(source, 0, ages, r, v, t, source, observer=observer)
+        super().__init__(source, 0, -dt, r, v, t, source, observer=observer)
+
+    @property
+    def dt(self) -> u.Quantity:
+        return -self.ages
 
     @property
     def epoch(self) -> Time:
@@ -777,13 +782,9 @@ def source_orbit(self, dt: u.Quantity) -> State | tuple[State, SkyCoord]:
         Returns
         -------
 
-        orbit : State
+        orbit : SourceOrbit
             The orbital states.
 
-        coords : SkyCoord, optional
-            The observed coordinates.  Only returned when ``.observer`` is
-            defined.
-
         """
 
         states: list[State] = []
@@ -792,8 +793,4 @@ def source_orbit(self, dt: u.Quantity) -> State | tuple[State, SkyCoord]:
             states.append(self.solver.solve(self.source, t, 0))
         states: State = State.from_states(states)
 
-        if self.observer is None:
-            return states
-
-        coords: SkyCoord = self.observer.observe(states)
-        return states, coords
+        return SourceOrbit(self.source, dt, states.r, states.v, states.t, self.observer)

From d258eaba5cac1342b364c37ca369efa18ea8c6bc Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Sun, 28 Dec 2025 10:24:55 -0500
Subject: [PATCH 07/34] Fix test_orbit

---
 sbpy/dynamics/tests/test_syndynes.py | 7 +------
 1 file changed, 1 insertion(+), 6 deletions(-)

diff --git a/sbpy/dynamics/tests/test_syndynes.py b/sbpy/dynamics/tests/test_syndynes.py
index 005aac61..6b43feab 100644
--- a/sbpy/dynamics/tests/test_syndynes.py
+++ b/sbpy/dynamics/tests/test_syndynes.py
@@ -497,7 +497,7 @@ def test_orbit(self, example_syndynes):
         comet, betas, ages, dust, observer = example_syndynes
 
         dt = [-1, 0, 1] * u.d
-        orbit, coords = dust.source_orbit(dt)
+        orbit = dust.source_orbit(dt)
         assert u.allclose((orbit.t - comet.t).to("s"), dt)
 
         solver = SolarGravity()
@@ -505,8 +505,3 @@ def test_orbit(self, example_syndynes):
             expected = solver.solve(comet, comet.t + dt[i])
             assert u.allclose(orbit[i].r, expected.r, rtol=1e-11, atol=1 * u.cm)
             assert u.allclose(orbit[i].v, expected.v, rtol=1e-11, atol=1 * u.um / u.s)
-
-        # test without observer
-        dust.observer = None
-        orbit2 = dust.source_orbit(dt)
-        assert np.allclose(orbit2.r, orbit.r)

From d9b1d08664a6079bcab7b96e300db6756499ffb9 Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Sun, 28 Dec 2025 10:25:59 -0500
Subject: [PATCH 08/34] Update examples, show how to do more complex plots

---
 docs/sbpy/dynamics.rst | 185 +++++++++++++++++++++--------------------
 1 file changed, 95 insertions(+), 90 deletions(-)

diff --git a/docs/sbpy/dynamics.rst b/docs/sbpy/dynamics.rst
index 0acb2ee8..fcfef001 100644
--- a/docs/sbpy/dynamics.rst
+++ b/docs/sbpy/dynamics.rst
@@ -168,6 +168,24 @@ Or, get the sky coordinates of the comet as seen from the Earth:
 The result, a `~astropy.coordinates.SkyCoord` object, is expressed in the reference frame of the observer.
 
 
+Fetching states from Horizons
+-----------------------------
+
+`Ephem.from_horizons` returns equatorial coordinates in the ICRF reference frame, which has its origin at the Solar System barycenter.  For `State` to correctly convert the ephemeris object to vectors, we need to set the Horizons observer to the Solar System barycenter (``"@ssb"``).  However, dynamical integrations are done in a heliocentric reference frame, so we transform the result to ``"heliocentriceclipticiau76"``:
+
+.. doctest-requires:: astroquery
+
+   >>> eph = Ephem.from_horizons(
+   ...     "48P",
+   ...     id_type="designation",
+   ...     closest_apparition=True,
+   ...     epochs=obstime,
+   ...     location="@ssb",
+   ... )
+   >>> comet = State.from_ephem(eph, frame="icrs")
+   >>> comet = comet.transform_to("heliocentriceclipticiau76")
+
+
 Dynamical integrators
 =====================
 
@@ -293,24 +311,6 @@ To support generating synchrones at specific times, the :func:`~sbpy.dynamics.sy
    <Time object: scale='utc' format='iso' value=2023-10-01 00:00:00.000>
 
 
-Using Horizons
---------------
-
-`Ephem.from_horizons` returns equatorial coordinates in the ICRF reference frame, which has its origin at the Solar System barycenter.  For `State` to correctly convert the ephemeris object to vectors, we need to set the Horizons observer to the Solar System barycenter (``"@ssb"``).  However, dynamical integrations are done in a heliocentric reference frame, so we transform the result to ``"heliocentriceclipticiau76"``:
-
-.. doctest-requires:: astroquery
-
-   >>> eph = Ephem.from_horizons(
-   ...     "48P",
-   ...     id_type="designation",
-   ...     closest_apparition=True,
-   ...     epochs=obstime,
-   ...     location="@ssb",
-   ... )
-   >>> comet = State.from_ephem(eph, frame="icrs")
-   >>> comet = comet.transform_to("heliocentriceclipticiau76")
-
-
 Adding ejection velocity
 ------------------------
 
@@ -381,7 +381,8 @@ Calculating the positions of the projected orbit of the source object may be hel
 .. doctest-requires:: scipy
 
    >>> dt = np.linspace(-2, 2) * u.d
-   >>> orbit, coords = dust.source_orbit(dt)
+   >>> orbit = dust.source_orbit(dt)
+   >>> orbit.coords
 
 
 Using other dynamical models
@@ -406,45 +407,25 @@ In this example, we compute the syndynes of a comet orbiting β Pic (1.8 solar m
 Plotting syndynes and synchrones
 --------------------------------
 
-Generally, we are interested in visualizing syndynes and synchrones for an observer.  We can plot the coordinates themselves, or plot them relative to the comet.
-
-.. note::
-   The accuracy of the coordinates object depends on the the comet and observer states, but also on whether or not light travel time is accounted for, and the accuracy of the orbit integrator.  The `sbpy` testing suite shows that arcsecond-level accuracy is possible, but this is generally not enough for direct comparison to typical images of comets, which need sub-arcsecond alignment.  For small offsets, plotting results relative to the comet's calculated position is sufficient.
+Generally, we are interested in visualizing the syndynes and synchrones for an observer.  `Syndynes`, `Synchrones`, and `SourceOrbit` have ``plot()`` methods to assist with this.  They can plot the coordinates relative to the comet with a simple tangent plane projection, or projected onto the image plane with an `astropy.wcs.WCS` object.
 
-Here is a simple example that plots the syndynes and synchrones as offsets from the comet's coordinates:
+Here is a simple example that plots the syndynes and synchrones from above as offsets from the comet's coordinates:
 
 .. doctest-requires:: scipy,matplotlib
 
    >>> import matplotlib.pyplot as plt
    >>>
-   >>> # comet's caluclated coordinates
-   >>> coords0 = observer.observe(comet)
-   >>>
-   >>> def plot(ax, coords, **kwargs):
-   ...     "plot coordinates relative to comet"
-   ...     dRA = coords.ra - coords0.ra
-   ...     dDec = coords.dec - coords0.dec
-   ...     ax.plot(dRA.arcsec, dDec.arcsec, **kwargs)
-   >>>
    >>> fig, ax = plt.subplots()
    >>>
    >>> # plot all but the last (beta=0) syndyne
-   >>> for syndyne in dust.syndynes()[:-1]:
-   ...     plot(ax, syndyne.coords, label=f"$\\beta={syndyne.beta:.2g}$")
+   >>> dust.syndynes()[:-1].plot(ax)
    >>>
    >>> # plot every 5th synchrone
-   >>> for synchrone in dust.synchrones()[4::5]:
-   ...     plot(
-   ...         ax,
-   ...         synchrone.coords,
-   ...         ls="--",
-   ...         label=f"$\\Delta t={synchrone.age.to(u.d):.2g}$"
-   ...     )
+   >>> dust.synchrones()[4::5].plot(ax)
    >>>
-   >>> # and plot the orbit
+   >>> # plot the orbit
    >>> dt = np.linspace(-2, 2) * u.d
-   >>> states, coords = dust.source_orbit(dt)
-   >>> plot(ax, coords, color="k", ls=":", label="Orbit")
+   >>> dust.source_orbit(dt).plot(ax, color="k", ls=":", label="Orbit")
    >>>
    >>> ax.invert_xaxis()
    >>> plt.setp(ax,
@@ -457,6 +438,7 @@ Here is a simple example that plots the syndynes and synchrones as offsets from
    >>> plt.tight_layout()
 
 .. plot::
+   :context:
 
    import numpy as np
    import matplotlib.pyplot as plt
@@ -481,45 +463,34 @@ Here is a simple example that plots the syndynes and synchrones as offsets from
    )
    dust = SynGenerator(comet, betas, ages, observer=observer)
 
-   coords0 = observer.observe(comet)
-   def plot(ax, coords, **kwargs):
-       dRA = coords.ra - coords0.ra
-       dDec = coords.dec - coords0.dec
-       ax.plot(dRA.arcsec, dDec.arcsec, **kwargs)
-
+   # plot
    fig, ax = plt.subplots()
 
-   for syndyne in dust.syndynes():
-       # don't draw the beta = 0 syndyne
-       if syndyne.beta == 0:
-           continue
-       plot(ax, syndyne.coords, label=f"$\\beta={syndyne.beta:.2g}$")
+   # plot all but the last (beta=0) syndyne
+   dust.syndynes()[:-1].plot(ax)
 
    # plot every 5th synchrone
-   for synchrone in dust.synchrones()[4::5]:
-       plot(
-           ax,
-           synchrone.coords,
-           ls="--",
-           label=f"$\Delta t={synchrone.age.to(u.d):.2g}$",
-       )
-
-   # and plot the orbit
+   dust.synchrones()[4::5].plot(ax, ls="--", lw=1)
+
+   # plot the orbit
    dt = np.linspace(-2, 2) * u.d
-   states, coords = dust.source_orbit(dt)
-   plot(ax, coords, color="k", ls=":", label="Orbit")
+   dust.source_orbit(dt).plot(ax, color="k", ls=":", label="Orbit")
 
    ax.invert_xaxis()
    plt.setp(ax,
             xlim=[100, -10],
             ylim=[-10, 100],
-            xlabel="$\Delta$RA (arcsec)",
-            ylabel="$\Delta$Dec (arcsec)",
+            xlabel="$\\Delta$RA (arcsec)",
+            ylabel="$\\Delta$Dec (arcsec)",
    )
    plt.legend()
    plt.tight_layout()
 
-The following example compares syndynes to a Spitzer Space Telesocpe image of comet 48P/Johnson (`Reach et al. 2007 <https://scixplorer.org/abs/2007Icar..191..298R/abstract>`_).  The FITS world coordinate system is used to account for the image orientation and scale.  To precisely align the syndynes with the comet nucleus, we update the world coordinate system to use our calculated comet coordinates:
+
+The following example compares syndynes to a Spitzer Space Telesocpe image of comet 48P/Johnson (`Reach et al. 2007 <https://scixplorer.org/abs/2007Icar..191..298R/abstract>`_).  The FITS world coordinate system is used to account for the image orientation and scale.  To precisely align the syndynes with the comet nucleus, we update the world coordinate system to use our calculated comet coordinates.
+
+.. note::
+   The `sbpy` testing suite shows that arcsecond-level accuracy is possible, but this is generally not enough for direct comparison to typical images of comets, which need sub-arcsecond alignment.  The accuracy of the coordinates object depends on the the comet and observer states, but also on whether or not light travel time is accounted for, and the accuracy of the orbit integrator.
 
 .. doctest-requires:: scipy,astroquery,matplotlib
 
@@ -559,12 +530,21 @@ The following example compares syndynes to a Spitzer Space Telesocpe image of co
    >>> # plot the image and syndynes
    >>> fig, ax = plt.subplots(num=1, clear=True, figsize=(6.5, 3.25))
    >>> 
-   >>> ax.imshow(image, origin="lower", vmin=14, vmax=70, cmap="gray_r")
+   >>> ax.imshow(image, origin="lower", vmin=49.1, vmax=49.5, cmap="gray_r")
    >>> 
-   >>> # plot the syndynes
-   >>> for syndyne in dust.syndynes():
-   ...     x, y = wcs.world_to_pixel(syndyne.coords.unmasked)
-   ...     ax.plot(x, y, label=f"$\\beta={syndyne.beta:.2g}$")
+   >>> # save xlim and ylim for later
+   >>> xlim = ax.get_xlim()
+   >>> ylim = ax.get_ylim()
+   >>> 
+   >>> # plot syndynes
+   >>> dust.syndynes().plot(ax, wcs=wcs)
+   >>> 
+   >>> # plot the orbit
+   >>> dt = np.linspace(-1, 1) * u.d
+   >>> dust.source_orbit(dt).plot(ax, wcs=wcs, color="tab:cyan", lw=1, label="Orbit")
+   >>> 
+   >>> plt.setp(ax, xlim=xlim, ylim=ylim)
+   >>> plt.legend()
 
 .. plot::
    :show-source-link:
@@ -619,28 +599,53 @@ The following example compares syndynes to a Spitzer Space Telesocpe image of co
    
    # plot
    fig, ax = plt.subplots(num=1, clear=True, figsize=(6.5, 3.25))
-   
-   ax.imshow(image, origin="lower", vmin=14, vmax=70, cmap="gray_r")
-   
+      
+   ax.imshow(image, origin="lower", vmin=49.1, vmax=49.5, cmap="gray_r")
+
    # save xlim and ylim for later
    xlim = ax.get_xlim()
    ylim = ax.get_ylim()
-   
-   # plot the syndynes
-   for syndyne in dust.syndynes():
-       x, y = wcs.world_to_pixel(syndyne.coords.unmasked)
-       ax.plot(x, y, label=f"$\\beta={syndyne.beta:.2g}$")
-   
-   # plot orbit
-   dt = np.linspace(-1, 1) * 3 * u.d
-   states, coords = dust.source_orbit(dt)
-   x, y = wcs.world_to_pixel(coords.unmasked)
-   ax.plot(x, y, color="tab:olive", ls="--", label="Orbit")
-   
+
+   # plot syndynes
+   dust.syndynes().plot(ax, wcs=wcs)
+
+   # plot the orbit
+   dt = np.linspace(-1, 1) * u.d
+   dust.source_orbit(dt).plot(ax, wcs=wcs, color="tab:cyan", lw=1, label="Orbit")
+
    plt.setp(ax, xlim=xlim, ylim=ylim)
    plt.legend()
    plt.tight_layout()
 
+For more complex plot logic, e.g., to use specific line colors and styles, we can use the plot methods of the individual syndynes/synchrones:
+
+.. doctest-requires:: scipy,matplotlib
+
+   >>> ls = ["-", "--", "-."]
+   >>> syndynes = dust.syndynes()
+   >>> for i in range(3):
+   ...     syndynes[i].plot(ax, color="k", ls=ls[i])
+
+Note that legend labels are not automatically generated in this case.
+
+.. plot::
+   :context:
+
+   fig, ax = plt.subplots()
+
+   ls = ["-", "--", "-."]
+   syndynes = dust.syndynes()
+   for i in range(3):
+      syndynes[i].plot(ax, color="k", ls=ls[i])
+      ax.invert_xaxis()
+
+   plt.setp(ax,
+            xlim=[100, -10],
+            ylim=[-10, 100],
+            xlabel="$\\Delta$RA (arcsec)",
+            ylabel="$\\Delta$Dec (arcsec)",
+   )
+   plt.tight_layout()
 
 
 Reference/API

From fa6e71a8bc118d71641d9b7798f176ddcb39cb7b Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Sun, 28 Dec 2025 10:26:10 -0500
Subject: [PATCH 09/34] Update change log

---
 CHANGES.rst | 18 +++++++++++++++++-
 1 file changed, 17 insertions(+), 1 deletion(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index c6a1722b..3f8ec0e1 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -21,6 +21,9 @@ sbpy.photometry
 
 - Added Rubin Observatory's LSSTCam filter set to `bandpass()`. [#431]
 
+- New `sbpy.dynamics.syndynes.SourceOrbit` class to encapsulate a collection of
+  points along an object's orbit. [#...]
+
 
 API Changes
 -----------
@@ -33,7 +36,20 @@ sbpy.dynamics
   will return a single `sbpy.dynamics.syndynes.Syndyne` or
   `sbpy.dynamics.syndynes.Synchrone`.  Indexing with a tuple or slice, e.g.,
   ``syndynes[:2]`` will return `sbpy.dynamics.syndynes.Syndynes` or
-  `sbpy.dynamics.syndynes.Synchrones`.  Previously a list was returned.
+  `sbpy.dynamics.syndynes.Synchrones`.  Previously a list was returned. [#...]
+
+- `sbpy.dynamics.syndynes.SynGenerator.source_orbit` now returns a
+  `sbpy.dynamics.syndynes.SourceOrbit` object.  Previously the states and
+  projected coordinates were returned as separate objects. [#...]
+
+
+Bug Fixes
+---------
+
+sbpy.dynamics
+^^^^^^^^^^^^^
+
+- Added missing `sbpy.dynamics.state.State.to_ephem`. [#...]
 
 
 v0.6.0 (2025-12-02)

From e3a76e6f20d00164d791e76d6e43b665aeb4b0f0 Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Sun, 28 Dec 2025 10:35:11 -0500
Subject: [PATCH 10/34] Debug tests

---
 docs/sbpy/dynamics.rst | 15 +++++++++------
 1 file changed, 9 insertions(+), 6 deletions(-)

diff --git a/docs/sbpy/dynamics.rst b/docs/sbpy/dynamics.rst
index fcfef001..ea11a0ae 100644
--- a/docs/sbpy/dynamics.rst
+++ b/docs/sbpy/dynamics.rst
@@ -179,7 +179,7 @@ Fetching states from Horizons
    ...     "48P",
    ...     id_type="designation",
    ...     closest_apparition=True,
-   ...     epochs=obstime,
+   ...     epochs=Time("2004-10-13T21:08:23.894"),
    ...     location="@ssb",
    ... )
    >>> comet = State.from_ephem(eph, frame="icrs")
@@ -380,9 +380,12 @@ Calculating the positions of the projected orbit of the source object may be hel
 
 .. doctest-requires:: scipy
 
-   >>> dt = np.linspace(-2, 2) * u.d
+   >>> dt = np.linspace(-2, 2, 3) * u.d
    >>> orbit = dust.source_orbit(dt)
-   >>> orbit.coords
+   >>> print("\n".join(orbit.coords.to_string("hmsdms", precision=0)))
+   20h57m34s -35d17m04s
+   20h59m23s -35d18m48s
+   21h01m13s -35d20m36s
 
 
 Using other dynamical models
@@ -530,7 +533,7 @@ The following example compares syndynes to a Spitzer Space Telesocpe image of co
    >>> # plot the image and syndynes
    >>> fig, ax = plt.subplots(num=1, clear=True, figsize=(6.5, 3.25))
    >>> 
-   >>> ax.imshow(image, origin="lower", vmin=49.1, vmax=49.5, cmap="gray_r")
+   >>> ax.imshow(image, origin="lower", vmin=49.1, vmax=49.5, cmap="gray_r")  # doctest: +SKIP
    >>> 
    >>> # save xlim and ylim for later
    >>> xlim = ax.get_xlim()
@@ -543,8 +546,8 @@ The following example compares syndynes to a Spitzer Space Telesocpe image of co
    >>> dt = np.linspace(-1, 1) * u.d
    >>> dust.source_orbit(dt).plot(ax, wcs=wcs, color="tab:cyan", lw=1, label="Orbit")
    >>> 
-   >>> plt.setp(ax, xlim=xlim, ylim=ylim)
-   >>> plt.legend()
+   >>> plt.setp(ax, xlim=xlim, ylim=ylim)  # doctest: +SKIP
+   >>> plt.legend()  # doctest: +SKIP
 
 .. plot::
    :show-source-link:

From 4cb0afc225f537ba1be40c36571c461915f1ea6f Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Sun, 28 Dec 2025 10:36:19 -0500
Subject: [PATCH 11/34] Remove note

---
 docs/sbpy/dynamics.rst | 2 --
 1 file changed, 2 deletions(-)

diff --git a/docs/sbpy/dynamics.rst b/docs/sbpy/dynamics.rst
index ea11a0ae..940f75d2 100644
--- a/docs/sbpy/dynamics.rst
+++ b/docs/sbpy/dynamics.rst
@@ -629,8 +629,6 @@ For more complex plot logic, e.g., to use specific line colors and styles, we ca
    >>> for i in range(3):
    ...     syndynes[i].plot(ax, color="k", ls=ls[i])
 
-Note that legend labels are not automatically generated in this case.
-
 .. plot::
    :context:
 

From 71b151fcf144a5479955a2702db6571dfa8be1db Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Sun, 28 Dec 2025 14:25:50 -0500
Subject: [PATCH 12/34] unmasked should be a private method.

---
 sbpy/utils/core.py            | 4 ++--
 sbpy/utils/tests/test_core.py | 8 ++++----
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/sbpy/utils/core.py b/sbpy/utils/core.py
index b0ab6a44..458d5d2f 100644
--- a/sbpy/utils/core.py
+++ b/sbpy/utils/core.py
@@ -6,7 +6,7 @@
 
 """
 
-__all__ = ["required_packages", "optional_packages", "unmasked"]
+__all__ = ["required_packages", "optional_packages"]
 
 from importlib import import_module
 from warnings import warn
@@ -112,7 +112,7 @@ def optional_packages(*packages, message=None):
     return True
 
 
-def unmasked(array):
+def _unmasked(array):
     """Return an unmasked version of the array."""
 
     if get_data_and_mask is None:
diff --git a/sbpy/utils/tests/test_core.py b/sbpy/utils/tests/test_core.py
index 2020dab9..d316aa83 100644
--- a/sbpy/utils/tests/test_core.py
+++ b/sbpy/utils/tests/test_core.py
@@ -1,7 +1,7 @@
 import pytest
 import numpy as np
 from astropy.utils.masked import Masked
-from ..core import required_packages, optional_packages, unmasked
+from ..core import required_packages, optional_packages, _unmasked
 from ...exceptions import RequiredPackageUnavailable, OptionalPackageUnavailable
 
 
@@ -37,6 +37,6 @@ def test_unmasked():
     ma1 = Masked([1.0, 2.0, 3.0], mask=[False, False, True])
     ma2 = np.ma.MaskedArray(array, mask=mask)
 
-    assert unmasked(array) is array
-    assert (unmasked(ma1) == array).all()
-    assert (unmasked(ma2) == array).all()
+    assert _unmasked(array) is array
+    assert (_unmasked(ma1) == array).all()
+    assert (_unmasked(ma2) == array).all()

From de73b9747ba6a04089881ee0651a4da928e4a74d Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Sun, 28 Dec 2025 16:22:54 -0500
Subject: [PATCH 13/34] Fix usage of _unmasked

---
 sbpy/dynamics/syndynes.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/sbpy/dynamics/syndynes.py b/sbpy/dynamics/syndynes.py
index 7f0c012f..51489dbe 100644
--- a/sbpy/dynamics/syndynes.py
+++ b/sbpy/dynamics/syndynes.py
@@ -7,6 +7,8 @@
 
 """
 
+from __future__ import annotations
+
 __all__ = [
     "SourceOrbit",
     "SynGenerator",
@@ -38,7 +40,7 @@
 
 from ..data import Ephem
 from ..utils.decorators import requires
-from ..utils import unmasked
+from ..utils.core import _unmasked
 from .models import DynamicalModel, SolarGravityAndRadiationPressure
 from .state import StateBase, State
 
@@ -229,7 +231,7 @@ def plot(
         else:
             # convert coordinates to plot units with the WCS object (avoid
             # passing a masked object, or wcs will complain)
-            x, y = wcs.world_to_pixel(unmasked(self.coords))
+            x, y = wcs.world_to_pixel(_unmasked(self.coords))
 
         ax.plot(x, y, **kwargs)
 

From 60e60b17cc8c1c1cd21b649e65cea59120aadb88 Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Wed, 31 Dec 2025 09:28:54 -0500
Subject: [PATCH 14/34] Add missing lower case version

---
 sbpy/data/__init__.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/sbpy/data/__init__.py b/sbpy/data/__init__.py
index 0a5a471b..310a8403 100644
--- a/sbpy/data/__init__.py
+++ b/sbpy/data/__init__.py
@@ -274,6 +274,7 @@ class Conf:
             "description": "RA*cos(Dec) Rate",
             "fieldnames": [
                 "RA*cos(Dec)_rate",
+                "ra*cos(dec)_rate",
                 "dra cos(dec)",
                 "dRA cos(Dec)",
                 "dra*cos(dec)",

From 7ce267b544e0d7798878dd6bd15087bceac20fc7 Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Wed, 31 Dec 2025 09:32:50 -0500
Subject: [PATCH 15/34] Update to_ephem

* Move to_ephem to StateBase
* compute coordinates if an observer is given
* put frame in Ephem metadata
---
 sbpy/dynamics/state.py | 121 +++++++++++++++++++++++++++++++++--------
 1 file changed, 97 insertions(+), 24 deletions(-)

diff --git a/sbpy/dynamics/state.py b/sbpy/dynamics/state.py
index bff66a2c..91f937b8 100644
--- a/sbpy/dynamics/state.py
+++ b/sbpy/dynamics/state.py
@@ -294,6 +294,103 @@ def to_skycoord(self) -> SkyCoord:
         else:
             return SkyCoord(self._data, obstime=self.t, representation_type="cartesian")
 
+    def to_ephem(self, observer: StateType | None = None) -> Ephem:
+        r"""Convert to an sbpy ephemeris object.
+
+
+        Parameters
+        ----------
+        observer : `State`, optional
+            Calculate RA/longitude, Dec/latitude, distance, etc for this
+            observer.
+
+
+        Returns
+        -------
+        eph : Ephem
+
+
+        Notes
+        -----
+
+        Observer state is stored in the `Ephem.meta` attribute. ``coords`` is
+        calculated for the given observer.
+
+        =========================   ====================
+        Attribute or quantity       ``Ephem`` field name
+        =========================   ====================
+        t, as `Time`                date
+        t, as `Quantity`            t_relative
+        :math:`|r|`                 r
+        :math:`|v \cdot \hat{r}|`   rdot
+        coords.ra                   ra
+        coords.dec                  dec
+        coords.pm_ra_cosdec         ra*cos(dec)_rate
+        coords.pm_dec               dec_rate
+        coords.lon                  lon
+        coords.lat                  lat
+        coords.pm_lon_coslat        lon*cos(lat)_rate
+        coords.pm_lat               lat_rate
+        coords.distance             delta
+        coords.radial_velocity      deltadot
+        x                           x
+        y                           y
+        z                           z
+        v_x                         vx
+        v_y                         vy
+        v_z                         vz
+        =========================   ====================
+
+        """
+
+        data: dict = {}
+
+        if isinstance(self.t, Time):
+            data["date"] = self.t
+        else:
+            data["t_relative"] = self.t
+
+        data["r"] = abs(self)[0]
+        data["rdot"] = np.sum(self.r * self.v, 1) / np.sqrt(np.sum(self.r * self.r, 1))
+
+        data["x"] = self.x
+        data["y"] = self.y
+        data["z"] = self.z
+        data["vx"] = self.v_x
+        data["vy"] = self.v_y
+        data["vz"] = self.v_z
+
+        meta = {"frame": self.frame}
+
+        if observer is not None:
+            meta["observer"] = {
+                "r": self.observer.r,
+                "v": self.observer.v,
+                "t": self.observer.t,
+                "frame": self.observer.frame,
+            }
+
+            coords = observer.observe(self)
+
+            # use SkyCoords's to_table() method, which will account for when
+            # RA/Dec vs lon/lat are used.
+            tab = coords.to_table()
+
+            # convert SkyCoords's column names to Ephem's field names
+            skycoords_to_ephem = {
+                "pm_ra_cosdec": "ra*cos(dec)_rate",
+                "pm_lon_coslat": "lon*cos(lat)_rate",
+                "pm_dec": "dec_rate",
+                "pm_lat": "lat_rate",
+                "distance": "delta",
+                "radial_velocity": "deltadot",
+            }
+            for col in tab.colnames:
+                field = skycoords_to_ephem.get(col, col)
+                data[field] = tab[col]
+
+        return Ephem.from_dict(data, meta=meta)
+
     def transform_to(self, frame: FrameInputTypes) -> StateBaseType:
         """Transform state into another reference frame.
 
@@ -435,30 +532,6 @@ def from_ephem(
             " velocity fields."
         )
 
-    def to_ephem(self) -> Ephem:
-        """Convert to an ephemeris object.
-
-
-        Returns
-        -------
-        eph : `~sbpy.data.ephem.Ephem`
-
-        """
-
-        time_field = "dt" if self.arbitrary_time else "date"
-
-        eph = {
-            time_field: self.t,
-            "x": self.x,
-            "y": self.y,
-            "z": self.z,
-            "vx": self.v_x,
-            "vy": self.v_y,
-            "vz": self.v_z,
-        }
-
-        return Ephem.from_dict(eph, meta={"frame": self.frame})
-
     def observe(self, target: StateType) -> SkyCoord:
         """Project a target's position onto the sky.
 

From 61573d1c14a053b8c32cc0dbbdd012f76f1a95e4 Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Wed, 31 Dec 2025 09:41:19 -0500
Subject: [PATCH 16/34] Improve Syndyne(s) Synchrone(s)

to_ephem:
* Remove coords column from to_ephem due to problem vstacking SkyCoords columns when frame contains obstime
* Use StateBase's to_ephem for most of the work.

plot()
* Fix cosine correction
* Allow spherical coords other than RA/Dec

Fix SourceOrbit.epoch

SynCollection:
* __getitem__ should return SynCollections for tuples and slices
* Use Ephem's vstack in to_ephem
---
 sbpy/dynamics/syndynes.py | 96 ++++++++++++++++++---------------------
 1 file changed, 44 insertions(+), 52 deletions(-)

diff --git a/sbpy/dynamics/syndynes.py b/sbpy/dynamics/syndynes.py
index 51489dbe..dd7edcd5 100644
--- a/sbpy/dynamics/syndynes.py
+++ b/sbpy/dynamics/syndynes.py
@@ -35,7 +35,6 @@
 import astropy.units as u
 from astropy.wcs import WCS
 from astropy.time import Time
-from astropy.table import vstack
 from astropy.coordinates import SkyCoord
 
 from ..data import Ephem
@@ -107,11 +106,14 @@ def to_ephem(self) -> Ephem:
         t, as `Quantity`            t_relative
         :math:`|r|`                 r
         :math:`|v \cdot \hat{r}|`   rdot
-        coords                      coords
         coords.ra                   ra
         coords.dec                  dec
+        coords.pm_ra_cosdec         ra*cos(dec)_rate
+        coords.pm_dec               dec_rate
         coords.lon                  lon
         coords.lat                  lat
+        coords.pm_lon_coslat        lon*cos(lat)_rate
+        coords.pm_lat               lat_rate
         coords.distance             delta
         coords.radial_velocity      deltadot
         x                           x
@@ -131,59 +133,36 @@ def to_ephem(self) -> Ephem:
 
         """
 
-        data: dict = {}
-        data["beta_rad"] = self.betas
-        data["age"] = self.ages
+        eph = super().to_ephem(observer=self.observer)
 
-        if isinstance(self.t, Time):
-            data["date"] = self.t
-        else:
-            data["t_relative"] = self.t
-
-        data["r"] = abs(self)[0]
-        data["rdot"] = np.sum(self.r * self.v, 1) / np.sqrt(np.sum(self.r * self.r, 1))
-
-        if self.observer is not None:
-            for k, v in self.coords.representation_component_names.items():
-                if v in ("lon", "lat"):
-                    data[k] = getattr(self.coords, k)
-                elif v == "distance":
-                    data["delta"] = getattr(self.coords, k)
-            data["deltadot"] = self.coords.radial_velocity
-            data["coords"] = self.coords
-
-        data["x"] = self.x
-        data["y"] = self.y
-        data["z"] = self.z
-        data["vx"] = self.v_x
-        data["vy"] = self.v_y
-        data["vz"] = self.v_z
-        data["x initial"] = self.initial.x
-        data["y initial"] = self.initial.y
-        data["z initial"] = self.initial.z
-        data["vx initial"] = self.initial.v_x
-        data["vy initial"] = self.initial.v_y
-        data["vz initial"] = self.initial.v_z
-        data["t initial"] = self.initial.t
-
-        meta: dict = {}
-        meta["source"] = {
+        eph["beta_rad"] = self.betas
+        eph["age"] = self.ages
+
+        eph["x initial"] = self.initial.x
+        eph["y initial"] = self.initial.y
+        eph["z initial"] = self.initial.z
+        eph["vx initial"] = self.initial.v_x
+        eph["vy initial"] = self.initial.v_y
+        eph["vz initial"] = self.initial.v_z
+        eph["t initial"] = self.initial.t
+
+        eph.meta["source"] = {
             "r": self.source.r,
             "v": self.source.v,
             "t": self.source.t,
             "frame": self.source.frame,
         }
         if self.observer is None:
-            meta["observer"] = None
+            eph.meta["observer"] = None
         else:
-            meta["observer"] = {
+            eph.meta["observer"] = {
                 "r": self.observer.r,
                 "v": self.observer.v,
                 "t": self.observer.t,
                 "frame": self.observer.frame,
             }
 
-        return Ephem.from_dict(data, meta=meta)
+        return eph
 
     @requires("matplotlib")
     def plot(
@@ -224,10 +203,17 @@ def plot(
 
             coords0 = self.observer.observe(self.source)
 
-            dRA = self.coords.ra - coords0.ra
-            dDec = self.coords.dec - coords0.dec
-            x = dRA.to(unit) / np.cos(self.coords.dec)
-            y = dDec.to(unit)
+            # lon could be coords.ra or coords.lon
+            # lat could be coords.dec or coords.lat
+            for name, component in coords0.get_representation_component_names().items():
+                if component == "lon":
+                    dlon = getattr(self.coords, name) - getattr(coords0, name)
+                elif component == "lat":
+                    lat0 = getattr(coords0, name)
+                    dlat = getattr(self.coords, name) - lat0
+
+            x = dlon.to(unit) * np.cos(lat0)
+            y = dlat.to(unit)
         else:
             # convert coordinates to plot units with the WCS object (avoid
             # passing a masked object, or wcs will complain)
@@ -386,7 +372,7 @@ def dt(self) -> u.Quantity:
     @property
     def epoch(self) -> Time:
         """Epoch of each orbit point."""
-        return self.source.t - self.ages[0]
+        return self.source.t + self.dt
 
 
 class SynCollection:
@@ -426,9 +412,12 @@ def __len__(self) -> int:
         return len(self._data)
 
     def __getitem__(self, k: int | tuple | slice) -> SynStates | SynCollectionType:
-        # return a SynStates object
         if isinstance(k, int):
+            # return a SynStates object
             return self._data[k]
+        elif isinstance(k, tuple):
+            # return a SynCollection
+            return type(self)([self._data[i] for i in k])
 
         # return a SynCollection
         return type(self)(self._data[k])
@@ -443,12 +432,15 @@ def to_ephem(self) -> Ephem:
         if len(self) == 0:
             return Ephem()
 
-        tables: list[Ephem] = [s.to_ephem().table for s in self]
+        result = self[0].to_ephem()
+        for syn in self[1:]:
+            # remove metadata or else it will be appended to the tables[0]'s
+            # metadata
+            eph = syn.to_ephem()
+            eph.table.meta = {}
+            result.vstack(eph)
 
-        return Ephem.from_table(
-            vstack(tables, metadata_conflicts="error"),
-            meta=tables[0].meta,
-        )
+        return result
 
     @requires("matplotlib")
     def plot(

From 8dd3acb37f1e257d84e3e8d657fb65d724da05b1 Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Wed, 31 Dec 2025 09:42:32 -0500
Subject: [PATCH 17/34] Update syndyne tests for full coverage and more control

Replace random_syn* with fake_syn*
---
 sbpy/dynamics/tests/test_syndynes.py | 336 ++++++++++++++++++++-------
 1 file changed, 257 insertions(+), 79 deletions(-)

diff --git a/sbpy/dynamics/tests/test_syndynes.py b/sbpy/dynamics/tests/test_syndynes.py
index 6b43feab..dd92f415 100644
--- a/sbpy/dynamics/tests/test_syndynes.py
+++ b/sbpy/dynamics/tests/test_syndynes.py
@@ -5,7 +5,16 @@
 import astropy.units as u
 from astropy.coordinates.errors import ConvertError
 from astropy.time import Time
-from ..syndynes import Syndyne, Synchrone, Syndynes, Synchrones, SynGenerator
+from astropy.wcs import WCS
+from ..syndynes import (
+    Syndyne,
+    Syndynes,
+    Synchrone,
+    Synchrones,
+    SourceOrbit,
+    SynCollection,
+    SynGenerator,
+)
 from ..state import State
 from ..models import SolarGravity, SolarGravityAndRadiationPressure
 
@@ -13,6 +22,94 @@
 pytest.importorskip("scipy")
 
 
+def fake_state(source, betas, ages):
+    """Fake dynamics: particles are displaced beta * 1e4 km from the source for
+    a little more control on the tests."""
+    r0 = 10000 * u.km
+    v0 = 0.1 * u.km / u.s
+    length = len(betas)
+    return State(
+        (betas * r0)[:, np.newaxis] + source.r,
+        (betas * v0)[:, np.newaxis] + source.v,
+        source.t - ages,
+        frame=source.frame,
+    )
+
+
+def fake_syndynes():
+    observer = State(
+        [0, 1, 0] * u.au,
+        [30, 0, 0] * u.km / u.s,
+        Time("2025-12-29"),
+        frame="icrs",
+    )
+    comet = State(
+        [1, 0, 0] * u.au,
+        [0, 30, 0] * u.km / u.s,
+        observer.t,
+        frame="heliocentriceclipticiau76",
+    )
+
+    betas = 10.0 ** -np.arange(4)
+    ages = np.arange(10) * u.day
+
+    items = []
+    for beta in betas:
+        particles = fake_state(comet, np.repeat(beta, 10), ages * 0)
+        initial = fake_state(comet, np.zeros(10), ages)
+        items.append(
+            Syndyne(
+                comet,
+                beta,
+                ages,
+                particles.r,
+                particles.v,
+                particles.t,
+                initial,
+                observer=observer,
+            )
+        )
+
+    return comet, observer, betas, ages, items
+
+
+def fake_synchrones():
+    observer = State(
+        [0, 1, 0] * u.au,
+        [30, 0, 0] * u.km / u.s,
+        Time("2025-12-29"),
+        frame="icrs",
+    )
+    comet = State(
+        [1, 0, 0] * u.au,
+        [0, 30, 0] * u.km / u.s,
+        observer.t,
+        frame="heliocentriceclipticiau76",
+    )
+
+    betas = 10.0 ** -np.arange(4)
+    ages = np.arange(10) * u.day
+
+    items = []
+    for age in ages:
+        particles = fake_state(comet, betas, age)
+        initial = fake_state(comet, np.zeros(10), ages)
+        items.append(
+            Synchrone(
+                comet,
+                betas,
+                age,
+                particles.r,
+                particles.v,
+                particles.t,
+                initial,
+                observer=observer,
+            )
+        )
+
+    return comet, observer, betas, ages, items
+
+
 @pytest.fixture
 def example_syndynes():
     comet = State(
@@ -115,8 +212,6 @@ def test_to_ephem(self, example_syndynes):
         assert u.allclose(eph["deltadot"], deltadot, atol=1e-12 * u.km / u.s)
         assert u.allclose(eph["lon"], coords.lon)
         assert u.allclose(eph["lat"], coords.lat)
-        assert u.allclose(eph["coords"].lon, coords.lon)
-        assert u.allclose(eph["coords"].lat, coords.lat)
         assert u.allclose(eph["x"], particles.x)
         assert u.allclose(eph["y"], particles.y)
         assert u.allclose(eph["z"], particles.z)
@@ -161,10 +256,7 @@ def test_to_ephem(self, example_syndynes):
         assert np.all(eph["t_relative"] == 0 * u.s)
         assert eph.meta["observer"] is None
         assert np.all(
-            [
-                k not in eph
-                for k in ("delta", "deltadot", "ra", "dec", "lon", "lat", "coords")
-            ]
+            [k not in eph for k in ("delta", "deltadot", "ra", "dec", "lon", "lat")]
         )
 
     def test_init_1d(self, example_syndynes):
@@ -219,6 +311,24 @@ def test_getitem(self, example_syndynes):
         assert isinstance(dust.syndynes()[:1], Syndynes)
 
 
+class TestSynCollection:
+    def test_plot(self):
+        """Test plot's label_format when plotting a mix of syndynes and synchrones"""
+
+        plt = pytest.importorskip("matplotlib.pyplot")
+
+        _, _, _, _, syndynes = fake_syndynes()
+        _, _, _, _, synchrones = fake_synchrones()
+        collection = SynCollection(syndynes + synchrones)
+
+        _, ax = plt.subplots()
+        collection.plot()
+
+        for i, line in enumerate(ax.get_lines()):
+            # when label="", mpl sets it to _child0, _child1, ...
+            assert line.get_label().startswith("_child")
+
+
 class TestSynchrone:
     def test_init(self, example_syndynes):
         comet, betas, ages, dust, observer = example_syndynes
@@ -287,8 +397,8 @@ def test_to_ephem(self, example_syndynes):
 
         eph = synchrone.to_ephem()
 
-        # most everything is tested in TestSyndyne.test_to_ephem
-        # we just need to make sure the betas and ages are right
+        # most everything is tested in TestSyndyne.test_to_ephem let's just
+        # check the betas and ages
         assert np.allclose(eph["beta_rad"], betas)
         assert u.allclose(eph["age"], age)
 
@@ -299,90 +409,158 @@ def test_getitem(self, example_syndynes):
         assert isinstance(dust.synchrones()[:1], Synchrones)
 
 
-def random_state(length=1):
-    return State(
-        np.random.rand(length, 3).squeeze() * u.au,
-        np.random.rand(length, 3).squeeze() * u.km / u.s,
-        np.random.rand() * u.day,
-    )
+class TestSyndynes:
+    def test_init(self):
+        comet, _, betas, ages, items = fake_syndynes()
 
+        # create from list
+        syndynes = Syndynes(items)
 
-def test_syndynes():
-    comet = random_state()
-    betas = np.linspace(0, 1, 3)
-    ages = np.arange(10) * u.day
+        # create from syndynes object
+        syndynes2 = Syndynes(syndynes)
+        assert np.all(syndynes2[0].r == syndynes[0].r)
 
-    items = []
-    for beta in betas:
-        particles = random_state(len(ages))
-        initial = random_state(len(ages))
-        items.append(
-            Syndyne(
-                comet,
-                beta,
-                ages,
-                particles.r,
-                particles.v,
-                particles.t,
-                initial,
-            )
+        # create from tuple
+        syndynes3 = Syndynes(tuple(items))
+        assert np.all(syndynes3[0].r == syndynes[0].r)
+
+        # cannot mix syndynes and synchrones
+        _, _, _, _, items2 = fake_synchrones()
+        with pytest.raises(TypeError):
+            Syndynes(items + items2)
+
+        # test empty syndynes objects
+        assert len(Syndynes([])) == 0
+
+    def test_getitem(self):
+        _, _, _, _, items = fake_syndynes()
+        syndynes = Syndynes(items)
+
+        # get a single Syndyne
+        assert isinstance(syndynes[2], Syndyne)
+        assert syndynes[2] is items[2]
+
+        # get Syndynes from a slice
+        assert len(syndynes[:3]) == 3
+        assert isinstance(syndynes[:3], Syndynes)
+
+        # get Syndynes from a tuple
+        assert isinstance(syndynes[(0, 1)], Syndynes)
+
+    def test_repr(self):
+        _, _, _, _, items = fake_syndynes()
+        syndynes = Syndynes(items)
+        assert repr(syndynes) == "<Syndynes: betas=[1.    0.1   0.01  0.001]>"
+
+    def test_to_ephem(self):
+        comet, _, _, _, items = fake_syndynes()
+        syndynes = Syndynes(items)
+
+        eph = syndynes.to_ephem()
+        assert len(eph) == 40
+        assert all(eph.meta["source"]["r"] == comet.r)
+
+        # [1] is 2nd index of 1st syndyne
+        assert eph["x"][1] == items[0].x[1]
+
+        # [13] is 4th index of 2nd syndyne
+        assert eph["vz"][13] == items[1].v_z[3]
+
+        assert len(Syndynes([]).to_ephem()) == 0
+
+    def test_plot(self):
+        plt = pytest.importorskip("matplotlib.pyplot")
+
+        comet, observer, betas, ages, items = fake_syndynes()
+        syndynes = Syndynes(items)
+
+        _, ax = plt.subplots()
+        syndynes.plot()
+
+        # observer is sqrt(2) au away, particles are up to 1e-4 au = 14959 km =
+        # 15"
+        for i, line in enumerate(ax.get_lines()):
+            assert line.get_label() == f"$\\beta={betas[i]}$"
+            assert all(line.get_xdata() <= 15 * betas[i] * u.arcsec)
+            assert all(line.get_ydata() <= 15 * betas[i] * u.arcsec)
+
+        coords = observer.observe(comet)
+
+        # 1 deg/pix
+        wcs = WCS()
+        wcs.wcs.ctype = "RA---TAN", "DEC--TAN"
+        wcs.wcs.crval = coords.ra.deg, coords.dec.deg
+
+        syndynes.plot(wcs=wcs)
+
+        # 15" / (1 deg/pix) = 0.00417 pix
+        for i, line in enumerate(ax.get_lines()[len(syndynes) :]):
+            assert all(line.get_xdata() <= 4.17e-3 * betas[i])
+            assert all(line.get_ydata() <= 4.17e-3 * betas[i])
+
+        # cannot plot relative to source when observer is None
+        state = fake_state(comet, [1], ages)
+        syndyne = Syndyne(
+            comet, [1], ages, state.r, state.v, state.t, comet, observer=None
         )
+        with pytest.raises(ValueError):
+            syndyne.plot()
 
-    syndynes = Syndynes(items)
-    assert syndynes[2] is items[2]
-    assert len(syndynes[:3]) == 3
-    assert repr(syndynes) == "<Syndynes: betas=[0.  0.5 1. ]>"
 
-    eph = syndynes.to_ephem()
-    assert len(eph) == 30
-    assert all(eph.meta["source"]["r"] == comet.r)
+class TestSynchrones:
+    # most of the relevant tests are covered by TestSyndynes
 
-    syndynes2 = Syndynes(syndynes)
-    assert np.all(syndynes2[0].r == syndynes[0].r)
-    syndynes3 = Syndynes(tuple(items))
-    assert np.all(syndynes3[0].r == syndynes[0].r)
+    def test_init(self):
+        _, _, _, _, items = fake_synchrones()
+        synchrones = Synchrones(items)
 
-    particles = random_state(len(betas))
-    synchrone = Synchrones(
-        [
-            Synchrone(
-                comet, betas, ages[0], particles.r, particles.v, particles.t, particles
-            )
-        ]
-    )
-    with pytest.raises(TypeError):
-        Syndynes(items + [synchrone])
+        assert synchrones[2] is items[2]
+        assert len(synchrones) == len(items)
 
-    # test emptpy syndynes objects
-    assert len(Syndynes([])) == 0
-    assert len(Syndynes([]).to_ephem()) == 0
+    def test_getitem(self):
+        _, _, _, _, items = fake_synchrones()
+        synchrones = Synchrones(items)
 
+        # index with integer returns one synchrone
+        assert synchrones[2] is items[2]
+        assert isinstance(synchrones[0], Synchrone)
 
-def test_synchrones():
-    comet = random_state()
-    betas = np.linspace(0, 1, 3)
-    ages = np.arange(10) * u.day
+        # index with slice returns synchrone collection
+        assert isinstance(synchrones[:3], Synchrones)
 
-    items = []
-    for age in ages:
-        particles = random_state(len(ages))
-        initial = random_state(len(ages))
-        items.append(
-            Synchrone(
-                comet,
-                betas,
-                age,
-                particles.r,
-                particles.v,
-                particles.t,
-                initial,
-            )
+        # index with tuple returns synchrone collection
+        assert isinstance(synchrones[(0, 1)], Synchrones)
+
+    def test_repr(self):
+        _, _, _, _, items = fake_synchrones()
+        synchrones = Synchrones(items)
+        assert (
+            repr(synchrones) == "<Synchrones: ages=[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.] d>"
         )
 
-    synchrones = Synchrones(items)
-    assert synchrones[2] is items[2]
-    assert len(synchrones[:3]) == 3
-    assert repr(synchrones) == "<Synchrones: ages=[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.] d>"
+    def test_plot(self):
+        plt = pytest.importorskip("matplotlib.pyplot")
+
+        comet, observer, betas, ages, items = fake_synchrones()
+        synchrones = Synchrones(items)
+
+        _, ax = plt.subplots()
+        synchrones.plot()
+
+        for i, line in enumerate(ax.get_lines()):
+            assert line.get_label() == f"$\\Delta t={ages[i]}$"
+
+
+class TestSourceOrbit:
+    def test_dt(self):
+        comet = State([1, 0, 0] * u.au, [0, 30, 0] * u.km / u.s, Time("2023-12-07"))
+        dt = [-1, 0, 1] * u.day
+        ages = -dt
+        state = fake_state(comet, [0, 1, 2], ages)
+        orbit = SourceOrbit(comet, dt, state.r, state.v, comet.t + dt)
+
+        assert all(orbit.dt == dt)
+        assert all(orbit.epoch == comet.t + dt)
 
 
 class TestSynGenerator:

From d50c28088e249ee191ca3b4e29b0007ecb6e56e2 Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Wed, 31 Dec 2025 10:21:50 -0500
Subject: [PATCH 18/34] to_ephem was in the wrong spot

---
 sbpy/dynamics/state.py            | 33 ++++++++++++++++++++-----------
 sbpy/dynamics/syndynes.py         |  2 +-
 sbpy/dynamics/tests/test_state.py | 15 +++++++++++++-
 3 files changed, 36 insertions(+), 14 deletions(-)

diff --git a/sbpy/dynamics/state.py b/sbpy/dynamics/state.py
index 91f937b8..a1c3fdae 100644
--- a/sbpy/dynamics/state.py
+++ b/sbpy/dynamics/state.py
@@ -294,15 +294,21 @@ def to_skycoord(self) -> SkyCoord:
         else:
             return SkyCoord(self._data, obstime=self.t, representation_type="cartesian")
 
-    def to_ephem(self, observer: StateType | None = None) -> Ephem:
+    def to_ephem(
+        self, observer: StateType | None = None, coords: SkyCoord | None = None
+    ) -> Ephem:
         r"""Convert to an sbpy ephemeris object.
 
 
         Parameters
         ----------
         observer : `State`, optional
-            Calculate RA/longitude, Dec/latitude, distance, etc for this
-            observer.
+            Include this observer in the metadata.  If ``coords`` is ``None``,
+            then also calculate RA/longitude, Dec/latitude, distance, etc for
+            this observer.
+
+        coords : `~astropy.coordinates.SkyCoord`, optional
+            Include these observer-based coordinates in the result.
 
 
         Returns
@@ -364,20 +370,23 @@ def to_ephem(self, observer: StateType | None = None) -> Ephem:
 
         if observer is not None:
             meta["observer"] = {
-                "r": self.observer.r,
-                "v": self.observer.v,
-                "t": self.observer.t,
-                "frame": self.observer.frame,
+                "r": observer.r,
+                "v": observer.v,
+                "t": observer.t,
+                "frame": observer.frame,
             }
 
-            coords = observer.observe(self)
+            if coords is None:
+                coords = observer.observe(self)
 
-            # use SkyCoords's to_table() method, which will account for when
+            # use SkyCoord's to_table() method, which will account for when
             # RA/Dec vs lon/lat are used.
             tab = coords.to_table()
 
-            # convert SkyCoords's column names to Ephem's field names
-            skycoords_to_ephem = {
+            # convert SkyCoord's column names to Ephem's field names
+            skycoord_to_ephem = {
+                "pm_ra": "ra_rate",
+                "pm_lon": "lon_rate",
                 "pm_ra_cosdec": "ra*cos(dec)_rate",
                 "pm_lon_coslat": "lon*cos(lat)_rate",
                 "pm_dec": "dec_rate",
@@ -386,7 +395,7 @@ def to_ephem(self, observer: StateType | None = None) -> Ephem:
                 "radial_velocity": "deltadot",
             }
             for col in tab.colnames:
-                field = skycoords_to_ephem.get(col, col)
+                field = skycoord_to_ephem.get(col, col)
                 data[field] = tab[col]
 
         return Ephem.from_dict(data, meta=meta)
diff --git a/sbpy/dynamics/syndynes.py b/sbpy/dynamics/syndynes.py
index dd7edcd5..52ffefcf 100644
--- a/sbpy/dynamics/syndynes.py
+++ b/sbpy/dynamics/syndynes.py
@@ -133,7 +133,7 @@ def to_ephem(self) -> Ephem:
 
         """
 
-        eph = super().to_ephem(observer=self.observer)
+        eph = super().to_ephem(observer=self.observer, coords=self.coords)
 
         eph["beta_rad"] = self.betas
         eph["age"] = self.ages
diff --git a/sbpy/dynamics/tests/test_state.py b/sbpy/dynamics/tests/test_state.py
index 7ae638de..d21069c7 100644
--- a/sbpy/dynamics/tests/test_state.py
+++ b/sbpy/dynamics/tests/test_state.py
@@ -562,7 +562,7 @@ def test_from_ephem(self):
             with pytest.raises(ValueError):
                 State.from_ephem(incomplete)
 
-    def test_to_ephm(self):
+    def test_to_ephem(self):
         r = np.arange(1, 7).reshape([2, 3]) * u.au
         v = r.value / 10 * u.km / u.s
         t = [0, 1] * u.day
@@ -582,3 +582,16 @@ def test_to_ephm(self):
         eph = state.to_ephem()
         assert np.allclose(eph["date"].mjd, [51544.0, 51545.0])
         assert isinstance(eph.meta["frame"], ICRS)
+
+        # observer with an observer
+        observer = State(
+            [0, 0, 1] * u.au, [0, 0, 0] * u.km / u.s, state.t[0], frame="icrs"
+        )
+        eph = state.to_ephem(observer=observer)
+        coords = observer.observe(state)
+        assert all(eph["ra"] == coords.ra)
+        assert all(eph["dec"] == coords.dec)
+        assert all(eph["ra_rate"] == coords.pm_ra)
+        assert all(eph["dec_rate"] == coords.pm_dec)
+        assert all(eph["delta"] == coords.distance)
+        assert all(eph["deltadot"] == coords.radial_velocity)

From 259d89611659bc5c7d28cc91cdcb8d91ba89e816 Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Wed, 31 Dec 2025 10:26:15 -0500
Subject: [PATCH 19/34] Update change log

---
 CHANGES.rst | 17 +++++++++++------
 1 file changed, 11 insertions(+), 6 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index 3f8ec0e1..aa720976 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -16,14 +16,14 @@ sbpy.dynamics
   `sbpy.dynamics.syndynes.Synchrone`, and `sbpy.dynamics.syndynes.Synchrones`.
   [#434]
 
+- New `sbpy.dynamics.syndynes.SourceOrbit` class to encapsulate a collection of
+  points along an object's orbit. [#434]
+
 sbpy.photometry
 ^^^^^^^^^^^^^^^
 
 - Added Rubin Observatory's LSSTCam filter set to `bandpass()`. [#431]
 
-- New `sbpy.dynamics.syndynes.SourceOrbit` class to encapsulate a collection of
-  points along an object's orbit. [#...]
-
 
 API Changes
 -----------
@@ -36,11 +36,16 @@ sbpy.dynamics
   will return a single `sbpy.dynamics.syndynes.Syndyne` or
   `sbpy.dynamics.syndynes.Synchrone`.  Indexing with a tuple or slice, e.g.,
   ``syndynes[:2]`` will return `sbpy.dynamics.syndynes.Syndynes` or
-  `sbpy.dynamics.syndynes.Synchrones`.  Previously a list was returned. [#...]
+  `sbpy.dynamics.syndynes.Synchrones`.  Previously a list was returned. [#434]
 
 - `sbpy.dynamics.syndynes.SynGenerator.source_orbit` now returns a
   `sbpy.dynamics.syndynes.SourceOrbit` object.  Previously the states and
-  projected coordinates were returned as separate objects. [#...]
+  projected coordinates were returned as separate objects. [#434]
+
+- Converting `sbpy.dynamics.syndynes.Syndynes` and
+  `sbpy.dynamics.syndynes.Synchrones` to ``Ephem`` objects no longer include the
+  ``SkyCoord`` column "coords", but all the data (RA, Dec, etc.) should still be
+  present as separate columns. [#434]
 
 
 Bug Fixes
@@ -49,7 +54,7 @@ Bug Fixes
 sbpy.dynamics
 ^^^^^^^^^^^^^
 
-- Added missing `sbpy.dynamics.state.State.to_ephem`. [#...]
+- Added missing `sbpy.dynamics.state.State.to_ephem`. [#434]
 
 
 v0.6.0 (2025-12-02)

From 3f368cc3c23440cec285bdcdb6840694eca2e6ab Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Wed, 31 Dec 2025 19:34:58 -0500
Subject: [PATCH 20/34] Mark tests with remote-data and require
 pytest-doctestplus 1.6 to do so.

---
 docs/sbpy/dynamics.rst | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/docs/sbpy/dynamics.rst b/docs/sbpy/dynamics.rst
index 940f75d2..8da3fce7 100644
--- a/docs/sbpy/dynamics.rst
+++ b/docs/sbpy/dynamics.rst
@@ -56,6 +56,7 @@ Convert to/from `Ephem` and `SkyCoord`
 State objects may be initialized from ephemeris objects (`~sbpy.data.Ephem`), provided they contain time, and 3D position and velocity:
 
 .. doctest-requires:: astroquery
+.. doctest-remote-data::
 
     >>> from sbpy.data import Ephem
     >>> eph = Ephem.from_horizons("9P",
@@ -67,6 +68,7 @@ State objects may be initialized from ephemeris objects (`~sbpy.data.Ephem`), pr
 And `State` may be converted to an `Ephem` object:
 
 .. doctest-requires:: astroquery
+.. doctest-remote-data::
 
     >>> eph = tempel1.to_ephem()  # doctest: +REMOTE_DATA
 
@@ -174,6 +176,7 @@ Fetching states from Horizons
 `Ephem.from_horizons` returns equatorial coordinates in the ICRF reference frame, which has its origin at the Solar System barycenter.  For `State` to correctly convert the ephemeris object to vectors, we need to set the Horizons observer to the Solar System barycenter (``"@ssb"``).  However, dynamical integrations are done in a heliocentric reference frame, so we transform the result to ``"heliocentriceclipticiau76"``:
 
 .. doctest-requires:: astroquery
+.. doctest-remote-data::
 
    >>> eph = Ephem.from_horizons(
    ...     "48P",
@@ -496,6 +499,7 @@ The following example compares syndynes to a Spitzer Space Telesocpe image of co
    The `sbpy` testing suite shows that arcsecond-level accuracy is possible, but this is generally not enough for direct comparison to typical images of comets, which need sub-arcsecond alignment.  The accuracy of the coordinates object depends on the the comet and observer states, but also on whether or not light travel time is accounted for, and the accuracy of the orbit integrator.
 
 .. doctest-requires:: scipy,astroquery,matplotlib
+.. doctest-remote-data::
 
    >>> from astropy.io import fits
    >>> from astropy.wcs import WCS

From 8e53a4178649b5b4439b08f9e8661858574e637e Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Wed, 31 Dec 2025 19:45:08 -0500
Subject: [PATCH 21/34] Age should not be used.

The t is observation time.
---
 sbpy/dynamics/tests/test_syndynes.py | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/sbpy/dynamics/tests/test_syndynes.py b/sbpy/dynamics/tests/test_syndynes.py
index dd92f415..c19de24a 100644
--- a/sbpy/dynamics/tests/test_syndynes.py
+++ b/sbpy/dynamics/tests/test_syndynes.py
@@ -22,7 +22,7 @@
 pytest.importorskip("scipy")
 
 
-def fake_state(source, betas, ages):
+def fake_state(source, betas):
     """Fake dynamics: particles are displaced beta * 1e4 km from the source for
     a little more control on the tests."""
     r0 = 10000 * u.km
@@ -31,7 +31,7 @@ def fake_state(source, betas, ages):
     return State(
         (betas * r0)[:, np.newaxis] + source.r,
         (betas * v0)[:, np.newaxis] + source.v,
-        source.t - ages,
+        source.t,
         frame=source.frame,
     )
 
@@ -55,8 +55,8 @@ def fake_syndynes():
 
     items = []
     for beta in betas:
-        particles = fake_state(comet, np.repeat(beta, 10), ages * 0)
-        initial = fake_state(comet, np.zeros(10), ages)
+        particles = fake_state(comet, np.repeat(beta, 10))
+        initial = fake_state(comet, np.zeros(10))
         items.append(
             Syndyne(
                 comet,
@@ -92,8 +92,8 @@ def fake_synchrones():
 
     items = []
     for age in ages:
-        particles = fake_state(comet, betas, age)
-        initial = fake_state(comet, np.zeros(10), ages)
+        particles = fake_state(comet, betas)
+        initial = fake_state(comet, np.zeros(10))
         items.append(
             Synchrone(
                 comet,
@@ -499,7 +499,7 @@ def test_plot(self):
             assert all(line.get_ydata() <= 4.17e-3 * betas[i])
 
         # cannot plot relative to source when observer is None
-        state = fake_state(comet, [1], ages)
+        state = fake_state(comet, [1])
         syndyne = Syndyne(
             comet, [1], ages, state.r, state.v, state.t, comet, observer=None
         )
@@ -556,7 +556,7 @@ def test_dt(self):
         comet = State([1, 0, 0] * u.au, [0, 30, 0] * u.km / u.s, Time("2023-12-07"))
         dt = [-1, 0, 1] * u.day
         ages = -dt
-        state = fake_state(comet, [0, 1, 2], ages)
+        state = fake_state(comet, [0, 1, 2])
         orbit = SourceOrbit(comet, dt, state.r, state.v, comet.t + dt)
 
         assert all(orbit.dt == dt)

From edfb03abcb7e947f057e24b8362d905ee91fd84d Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Fri, 2 Jan 2026 09:34:00 -0500
Subject: [PATCH 22/34] Mark doctests with REMOTE_DATA

---
 docs/sbpy/dynamics.rst | 62 ++++++++++++++++++++----------------------
 1 file changed, 29 insertions(+), 33 deletions(-)

diff --git a/docs/sbpy/dynamics.rst b/docs/sbpy/dynamics.rst
index 8da3fce7..c63e02a8 100644
--- a/docs/sbpy/dynamics.rst
+++ b/docs/sbpy/dynamics.rst
@@ -56,7 +56,6 @@ Convert to/from `Ephem` and `SkyCoord`
 State objects may be initialized from ephemeris objects (`~sbpy.data.Ephem`), provided they contain time, and 3D position and velocity:
 
 .. doctest-requires:: astroquery
-.. doctest-remote-data::
 
     >>> from sbpy.data import Ephem
     >>> eph = Ephem.from_horizons("9P",
@@ -68,7 +67,6 @@ State objects may be initialized from ephemeris objects (`~sbpy.data.Ephem`), pr
 And `State` may be converted to an `Ephem` object:
 
 .. doctest-requires:: astroquery
-.. doctest-remote-data::
 
     >>> eph = tempel1.to_ephem()  # doctest: +REMOTE_DATA
 
@@ -176,7 +174,6 @@ Fetching states from Horizons
 `Ephem.from_horizons` returns equatorial coordinates in the ICRF reference frame, which has its origin at the Solar System barycenter.  For `State` to correctly convert the ephemeris object to vectors, we need to set the Horizons observer to the Solar System barycenter (``"@ssb"``).  However, dynamical integrations are done in a heliocentric reference frame, so we transform the result to ``"heliocentriceclipticiau76"``:
 
 .. doctest-requires:: astroquery
-.. doctest-remote-data::
 
    >>> eph = Ephem.from_horizons(
    ...     "48P",
@@ -184,9 +181,9 @@ Fetching states from Horizons
    ...     closest_apparition=True,
    ...     epochs=Time("2004-10-13T21:08:23.894"),
    ...     location="@ssb",
-   ... )
-   >>> comet = State.from_ephem(eph, frame="icrs")
-   >>> comet = comet.transform_to("heliocentriceclipticiau76")
+   ... )  # doctest: +REMOTE_DATA
+   >>> comet = State.from_ephem(eph, frame="icrs")  # doctest: +REMOTE_DATA
+   >>> comet = comet.transform_to("heliocentriceclipticiau76")  # doctest: +REMOTE_DATA
 
 
 Dynamical integrators
@@ -499,13 +496,12 @@ The following example compares syndynes to a Spitzer Space Telesocpe image of co
    The `sbpy` testing suite shows that arcsecond-level accuracy is possible, but this is generally not enough for direct comparison to typical images of comets, which need sub-arcsecond alignment.  The accuracy of the coordinates object depends on the the comet and observer states, but also on whether or not light travel time is accounted for, and the accuracy of the orbit integrator.
 
 .. doctest-requires:: scipy,astroquery,matplotlib
-.. doctest-remote-data::
 
    >>> from astropy.io import fits
    >>> from astropy.wcs import WCS
    >>>
-   >>> image, header = fits.getdata("https://sbpy.org/data/48p-spitzer-reach07.fits", header=True)
-   >>> obstime = Time(header["DATE_OBS"])
+   >>> image, header = fits.getdata("https://sbpy.org/data/48p-spitzer-reach07.fits", header=True)  # doctest: +REMOTE_DATA
+   >>> obstime = Time(header["DATE_OBS"])  # doctest: +REMOTE_DATA
    >>>
    >>> # get the comet state
    >>> eph = Ephem.from_horizons(
@@ -514,44 +510,44 @@ The following example compares syndynes to a Spitzer Space Telesocpe image of co
    ...     closest_apparition=True,
    ...     epochs=obstime,
    ...     location="@ssb",
-   ... )
-   >>> comet = State.from_ephem(eph, frame="icrs")
-   >>> comet = comet.transform_to("heliocentriceclipticiau76")
+   ... )  # doctest: +REMOTE_DATA
+   >>> comet = State.from_ephem(eph, frame="icrs")  # doctest: +REMOTE_DATA
+   >>> comet = comet.transform_to("heliocentriceclipticiau76")  # doctest: +REMOTE_DATA
    >>> 
    >>> # get the Spitzer Space Telescope state
-   >>> eph = Ephem.from_horizons("-79", id_type=None, epochs=obstime, location="@ssb")
-   >>> observer = State.from_ephem(eph, frame="icrs")
+   >>> eph = Ephem.from_horizons("-79", id_type=None, epochs=obstime, location="@ssb")  # doctest: +REMOTE_DATA
+   >>> observer = State.from_ephem(eph, frame="icrs")  # doctest: +REMOTE_DATA
    >>> 
    >>> # set up the world coordinate system object and update the origin to align with
    >>> # the calculated position of the comet
-   >>> wcs = WCS(header)
-   >>> coords0 = observer.observe(comet)[0].unmasked
-   >>> wcs.wcs.crval = coords0.ra.deg, coords0.dec.deg
-   >>> wcs.wcs.crpix = 209, 99
+   >>> wcs = WCS(header)  # doctest: +REMOTE_DATA
+   >>> coords0 = observer.observe(comet)[0].unmasked  # doctest: +REMOTE_DATA
+   >>> wcs.wcs.crval = coords0.ra.deg, coords0.dec.deg  # doctest: +REMOTE_DATA
+   >>> wcs.wcs.crpix = 209, 99  # doctest: +REMOTE_DATA
    >>>
    >>> # generate the syndynes
-   >>> betas = [1, 0.1, 0.01, 0.001]
-   >>> ages = np.linspace(0, 365, 51) * u.day
-   >>> dust = SynGenerator(comet[0], betas, ages, observer=observer[0])
+   >>> betas = [1, 0.1, 0.01, 0.001]  # doctest: +REMOTE_DATA
+   >>> ages = np.linspace(0, 365, 51) * u.day  # doctest: +REMOTE_DATA
+   >>> dust = SynGenerator(comet[0], betas, ages, observer=observer[0])  # doctest: +REMOTE_DATA
    >>> 
    >>> # plot the image and syndynes
-   >>> fig, ax = plt.subplots(num=1, clear=True, figsize=(6.5, 3.25))
+   >>> fig, ax = plt.subplots(num=1, clear=True, figsize=(6.5, 3.25))  # doctest: +REMOTE_DATA
    >>> 
-   >>> ax.imshow(image, origin="lower", vmin=49.1, vmax=49.5, cmap="gray_r")  # doctest: +SKIP
+   >>> ax.imshow(image, origin="lower", vmin=49.1, vmax=49.5, cmap="gray_r")  # doctest: +SKIP +REMOTE_DATA
    >>> 
    >>> # save xlim and ylim for later
-   >>> xlim = ax.get_xlim()
-   >>> ylim = ax.get_ylim()
+   >>> xlim = ax.get_xlim()  # doctest: +REMOTE_DATA
+   >>> ylim = ax.get_ylim()  # doctest: +REMOTE_DATA
    >>> 
    >>> # plot syndynes
-   >>> dust.syndynes().plot(ax, wcs=wcs)
+   >>> dust.syndynes().plot(ax, wcs=wcs)  # doctest: +REMOTE_DATA
    >>> 
    >>> # plot the orbit
-   >>> dt = np.linspace(-1, 1) * u.d
-   >>> dust.source_orbit(dt).plot(ax, wcs=wcs, color="tab:cyan", lw=1, label="Orbit")
+   >>> dt = np.linspace(-1, 1) * u.d  # doctest: +REMOTE_DATA
+   >>> dust.source_orbit(dt).plot(ax, wcs=wcs, color="tab:cyan", lw=1, label="Orbit")  # doctest: +REMOTE_DATA
    >>> 
-   >>> plt.setp(ax, xlim=xlim, ylim=ylim)  # doctest: +SKIP
-   >>> plt.legend()  # doctest: +SKIP
+   >>> plt.setp(ax, xlim=xlim, ylim=ylim)  # doctest: +SKIP +REMOTE_DATA
+   >>> plt.legend()  # doctest: +SKIP  # doctest: +REMOTE_DATA
 
 .. plot::
    :show-source-link:
@@ -626,12 +622,12 @@ The following example compares syndynes to a Spitzer Space Telesocpe image of co
 
 For more complex plot logic, e.g., to use specific line colors and styles, we can use the plot methods of the individual syndynes/synchrones:
 
-.. doctest-requires:: scipy,matplotlib
+.. doctest-requires:: scipy,astroquery,matplotlib
 
    >>> ls = ["-", "--", "-."]
-   >>> syndynes = dust.syndynes()
+   >>> syndynes = dust.syndynes()  # doctest: +REMOTE_DATA
    >>> for i in range(3):
-   ...     syndynes[i].plot(ax, color="k", ls=ls[i])
+   ...     syndynes[i].plot(ax, color="k", ls=ls[i])  # doctest: +REMOTE_DATA
 
 .. plot::
    :context:

From d33b41533ce5942fe2db4808223d58b87672e1b7 Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Fri, 2 Jan 2026 09:42:44 -0500
Subject: [PATCH 23/34] Fix copy-paste mistake

---
 docs/sbpy/dynamics.rst | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/sbpy/dynamics.rst b/docs/sbpy/dynamics.rst
index c63e02a8..4b027d6d 100644
--- a/docs/sbpy/dynamics.rst
+++ b/docs/sbpy/dynamics.rst
@@ -546,8 +546,8 @@ The following example compares syndynes to a Spitzer Space Telesocpe image of co
    >>> dt = np.linspace(-1, 1) * u.d  # doctest: +REMOTE_DATA
    >>> dust.source_orbit(dt).plot(ax, wcs=wcs, color="tab:cyan", lw=1, label="Orbit")  # doctest: +REMOTE_DATA
    >>> 
-   >>> plt.setp(ax, xlim=xlim, ylim=ylim)  # doctest: +SKIP +REMOTE_DATA
-   >>> plt.legend()  # doctest: +SKIP  # doctest: +REMOTE_DATA
+   >>> plt.setp(ax, xlim=xlim, ylim=ylim)  # doctest: +SKIP
+   >>> plt.legend()  # doctest: +SKIP
 
 .. plot::
    :show-source-link:

From be708b6affc74f58c1fec3899b0511fa9f4be58d Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Thu, 8 Jan 2026 08:53:01 -0500
Subject: [PATCH 24/34] Back to using doctest-remote-data

---
 docs/sbpy/dynamics.rst | 63 +++++++++++++++++++++++-------------------
 1 file changed, 34 insertions(+), 29 deletions(-)

diff --git a/docs/sbpy/dynamics.rst b/docs/sbpy/dynamics.rst
index 4b027d6d..c1f847cd 100644
--- a/docs/sbpy/dynamics.rst
+++ b/docs/sbpy/dynamics.rst
@@ -56,19 +56,21 @@ Convert to/from `Ephem` and `SkyCoord`
 State objects may be initialized from ephemeris objects (`~sbpy.data.Ephem`), provided they contain time, and 3D position and velocity:
 
 .. doctest-requires:: astroquery
+.. doctest-remote-data::
 
     >>> from sbpy.data import Ephem
     >>> eph = Ephem.from_horizons("9P",
     ...                           epochs=Time("2005-07-04"),
     ...                           id_type="designation",
-    ...                           closest_apparition=True)  # doctest: +REMOTE_DATA
-    >>> tempel1 = State.from_ephem(eph)                     # doctest: +REMOTE_DATA
+    ...                           closest_apparition=True)
+    >>> tempel1 = State.from_ephem(eph)
 
 And `State` may be converted to an `Ephem` object:
 
 .. doctest-requires:: astroquery
+.. doctest-remote-data::
 
-    >>> eph = tempel1.to_ephem()  # doctest: +REMOTE_DATA
+    >>> eph = tempel1.to_ephem()
 
 `astropy`'s `~astropy.coordinates.SkyCoord` objects may also be used, assuming the time and 3D vectors are fully defined:
 
@@ -174,6 +176,7 @@ Fetching states from Horizons
 `Ephem.from_horizons` returns equatorial coordinates in the ICRF reference frame, which has its origin at the Solar System barycenter.  For `State` to correctly convert the ephemeris object to vectors, we need to set the Horizons observer to the Solar System barycenter (``"@ssb"``).  However, dynamical integrations are done in a heliocentric reference frame, so we transform the result to ``"heliocentriceclipticiau76"``:
 
 .. doctest-requires:: astroquery
+.. doctest-remote-data::
 
    >>> eph = Ephem.from_horizons(
    ...     "48P",
@@ -181,9 +184,9 @@ Fetching states from Horizons
    ...     closest_apparition=True,
    ...     epochs=Time("2004-10-13T21:08:23.894"),
    ...     location="@ssb",
-   ... )  # doctest: +REMOTE_DATA
-   >>> comet = State.from_ephem(eph, frame="icrs")  # doctest: +REMOTE_DATA
-   >>> comet = comet.transform_to("heliocentriceclipticiau76")  # doctest: +REMOTE_DATA
+   ... )
+   >>> comet = State.from_ephem(eph, frame="icrs")
+   >>> comet = comet.transform_to("heliocentriceclipticiau76")
 
 
 Dynamical integrators
@@ -496,12 +499,13 @@ The following example compares syndynes to a Spitzer Space Telesocpe image of co
    The `sbpy` testing suite shows that arcsecond-level accuracy is possible, but this is generally not enough for direct comparison to typical images of comets, which need sub-arcsecond alignment.  The accuracy of the coordinates object depends on the the comet and observer states, but also on whether or not light travel time is accounted for, and the accuracy of the orbit integrator.
 
 .. doctest-requires:: scipy,astroquery,matplotlib
+.. doctest-remote-data::
 
    >>> from astropy.io import fits
    >>> from astropy.wcs import WCS
    >>>
-   >>> image, header = fits.getdata("https://sbpy.org/data/48p-spitzer-reach07.fits", header=True)  # doctest: +REMOTE_DATA
-   >>> obstime = Time(header["DATE_OBS"])  # doctest: +REMOTE_DATA
+   >>> image, header = fits.getdata("https://sbpy.org/data/48p-spitzer-reach07.fits", header=True)
+   >>> obstime = Time(header["DATE_OBS"])
    >>>
    >>> # get the comet state
    >>> eph = Ephem.from_horizons(
@@ -510,41 +514,41 @@ The following example compares syndynes to a Spitzer Space Telesocpe image of co
    ...     closest_apparition=True,
    ...     epochs=obstime,
    ...     location="@ssb",
-   ... )  # doctest: +REMOTE_DATA
-   >>> comet = State.from_ephem(eph, frame="icrs")  # doctest: +REMOTE_DATA
-   >>> comet = comet.transform_to("heliocentriceclipticiau76")  # doctest: +REMOTE_DATA
+   ... )
+   >>> comet = State.from_ephem(eph, frame="icrs")
+   >>> comet = comet.transform_to("heliocentriceclipticiau76")
    >>> 
    >>> # get the Spitzer Space Telescope state
-   >>> eph = Ephem.from_horizons("-79", id_type=None, epochs=obstime, location="@ssb")  # doctest: +REMOTE_DATA
-   >>> observer = State.from_ephem(eph, frame="icrs")  # doctest: +REMOTE_DATA
+   >>> eph = Ephem.from_horizons("-79", id_type=None, epochs=obstime, location="@ssb")
+   >>> observer = State.from_ephem(eph, frame="icrs")
    >>> 
    >>> # set up the world coordinate system object and update the origin to align with
    >>> # the calculated position of the comet
-   >>> wcs = WCS(header)  # doctest: +REMOTE_DATA
-   >>> coords0 = observer.observe(comet)[0].unmasked  # doctest: +REMOTE_DATA
-   >>> wcs.wcs.crval = coords0.ra.deg, coords0.dec.deg  # doctest: +REMOTE_DATA
-   >>> wcs.wcs.crpix = 209, 99  # doctest: +REMOTE_DATA
+   >>> wcs = WCS(header)
+   >>> coords0 = observer.observe(comet)[0].unmasked
+   >>> wcs.wcs.crval = coords0.ra.deg, coords0.dec.deg
+   >>> wcs.wcs.crpix = 209, 99
    >>>
    >>> # generate the syndynes
-   >>> betas = [1, 0.1, 0.01, 0.001]  # doctest: +REMOTE_DATA
-   >>> ages = np.linspace(0, 365, 51) * u.day  # doctest: +REMOTE_DATA
-   >>> dust = SynGenerator(comet[0], betas, ages, observer=observer[0])  # doctest: +REMOTE_DATA
+   >>> betas = [1, 0.1, 0.01, 0.001]
+   >>> ages = np.linspace(0, 365, 51) * u.day
+   >>> dust = SynGenerator(comet[0], betas, ages, observer=observer[0])
    >>> 
    >>> # plot the image and syndynes
-   >>> fig, ax = plt.subplots(num=1, clear=True, figsize=(6.5, 3.25))  # doctest: +REMOTE_DATA
+   >>> fig, ax = plt.subplots(num=1, clear=True, figsize=(6.5, 3.25))
    >>> 
-   >>> ax.imshow(image, origin="lower", vmin=49.1, vmax=49.5, cmap="gray_r")  # doctest: +SKIP +REMOTE_DATA
+   >>> ax.imshow(image, origin="lower", vmin=49.1, vmax=49.5, cmap="gray_r")  # docte
    >>> 
    >>> # save xlim and ylim for later
-   >>> xlim = ax.get_xlim()  # doctest: +REMOTE_DATA
-   >>> ylim = ax.get_ylim()  # doctest: +REMOTE_DATA
+   >>> xlim = ax.get_xlim()
+   >>> ylim = ax.get_ylim()
    >>> 
    >>> # plot syndynes
-   >>> dust.syndynes().plot(ax, wcs=wcs)  # doctest: +REMOTE_DATA
+   >>> dust.syndynes().plot(ax, wcs=wcs)
    >>> 
    >>> # plot the orbit
-   >>> dt = np.linspace(-1, 1) * u.d  # doctest: +REMOTE_DATA
-   >>> dust.source_orbit(dt).plot(ax, wcs=wcs, color="tab:cyan", lw=1, label="Orbit")  # doctest: +REMOTE_DATA
+   >>> dt = np.linspace(-1, 1) * u.d
+   >>> dust.source_orbit(dt).plot(ax, wcs=wcs, color="tab:cyan", lw=1, label="Orbit")
    >>> 
    >>> plt.setp(ax, xlim=xlim, ylim=ylim)  # doctest: +SKIP
    >>> plt.legend()  # doctest: +SKIP
@@ -623,11 +627,12 @@ The following example compares syndynes to a Spitzer Space Telesocpe image of co
 For more complex plot logic, e.g., to use specific line colors and styles, we can use the plot methods of the individual syndynes/synchrones:
 
 .. doctest-requires:: scipy,astroquery,matplotlib
+.. doctest-remote-data::
 
    >>> ls = ["-", "--", "-."]
-   >>> syndynes = dust.syndynes()  # doctest: +REMOTE_DATA
+   >>> syndynes = dust.syndynes()
    >>> for i in range(3):
-   ...     syndynes[i].plot(ax, color="k", ls=ls[i])  # doctest: +REMOTE_DATA
+   ...     syndynes[i].plot(ax, color="k", ls=ls[i])
 
 .. plot::
    :context:

From 8a47462ef5f12e8709f2514e8b6600b1d502c8f3 Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Thu, 7 May 2026 09:07:50 -0400
Subject: [PATCH 25/34] Move orbit section next to syndynes and synchrones

---
 docs/sbpy/dynamics.rst | 37 ++++++++++++++++++++++---------------
 1 file changed, 22 insertions(+), 15 deletions(-)

diff --git a/docs/sbpy/dynamics.rst b/docs/sbpy/dynamics.rst
index c1f847cd..623e61b1 100644
--- a/docs/sbpy/dynamics.rst
+++ b/docs/sbpy/dynamics.rst
@@ -314,6 +314,28 @@ To support generating synchrones at specific times, the :func:`~sbpy.dynamics.sy
    <Time object: scale='utc' format='iso' value=2023-10-01 00:00:00.000>
 
 
+Source object's orbit
+---------------------
+
+The `SynGenerator` can also generate points along the source object's oribit.  Calculating the positions of the projected orbit of the source object may be helpful for comparison to syndynes, or a comet's dust trail.  They are calculated with the :func:`~sbpy.dynamics.syndynes.SynGenerator.source_orbit` method.  To produce 3 points at -2, 0, and +2 days from the object's current position:
+
+.. doctest-requires:: scipy
+
+   >>> dt = np.linspace(-2, 2, 3) * u.d
+   >>> orbit = dust.source_orbit(dt)
+
+This method produced a `SourceOrbit` object that is functionally the same as the `Syndynes` and `Synchrones` objects above:
+
+.. doctest-requires:: scipy
+
+   >>> type(orbit)
+   sbpy.dynamics.syndynes.SourceOrbit
+   >>> orbit.ages
+   <Quantity [ 2., -0., -2.] d>
+   >>> np.linalg.norm(orbit[0].r - orbit[1].r)
+   <Quantity 5183919.42549707 km>
+
+
 Adding ejection velocity
 ------------------------
 
@@ -376,21 +398,6 @@ With the observer and coordinate frames defined, the generated syndynes and sync
    21h05m59s -33d30m26s
 
 
-Source object's projected orbit
--------------------------------
-
-Calculating the positions of the projected orbit of the source object may be helpful for interpreting an observation or a set of syndynes.  They are calculated with the :func:`~sbpy.dynamics.syndynes.SynGenerator.source_orbit` method:
-
-.. doctest-requires:: scipy
-
-   >>> dt = np.linspace(-2, 2, 3) * u.d
-   >>> orbit = dust.source_orbit(dt)
-   >>> print("\n".join(orbit.coords.to_string("hmsdms", precision=0)))
-   20h57m34s -35d17m04s
-   20h59m23s -35d18m48s
-   21h01m13s -35d20m36s
-
-
 Using other dynamical models
 ----------------------------
 

From 9d3e84d47acaec32f7974aa4de6e73aa5274a160 Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Thu, 7 May 2026 09:33:32 -0400
Subject: [PATCH 26/34] Expand preamble to describe
 SynGenerator->syndynes()->Syndynes relationship

---
 docs/sbpy/dynamics.rst | 43 ++++++++++++++++++++++++++++++------------
 1 file changed, 31 insertions(+), 12 deletions(-)

diff --git a/docs/sbpy/dynamics.rst b/docs/sbpy/dynamics.rst
index 623e61b1..f9403b41 100644
--- a/docs/sbpy/dynamics.rst
+++ b/docs/sbpy/dynamics.rst
@@ -224,13 +224,11 @@ Dust syndynes and synchrones
 
 Syndynes are lines in space connecting particles that are experiencing the same forces.  A syndyne is parameterized by :math:`\beta`, the ratio of the force from solar radiation to the force from solar gravity, :math:`F_r / F_g`, and age (or time of release).  Thus, all particles in a syndyne have a constant :math:`\beta` but variable age.  Similarly, synchrones are lines of constant particle age, but variable :math:`\beta`.
 
+sbpy has the `sbpy.dynamics.syndynes.SynGenerator` class to assist the user in the generation of syndynes, synchrones, and points along a source object's orbit.  Use the `~sbpy.dynamics.syndynes.SynGenerator`'s :func:`~sbpy.dynamics.syndynes.SynGenerator.syndynes`, :func:`~sbpy.dynamics.syndynes.SynGenerator.synchrones`, and :func:`~sbpy.dynamices.syndynes.SynGenerator.source_orbit` methods to produce `~sbpy.dynamics.syndynes.Syndynes`, `~sbpy.dynamics.syndynes.Synchrones`, and `~sbpy.dynamics.syndynes.SourceOrbit` objects, described in further detail below.
 
-Syndynes
---------
+First, setup a `~sbpy.dynamics.syndynes.SynGenerator` instance.  This requires a dust source described by a `~sbpy.dynamics.state.State` object, :math:`\beta` values, and particle ages from which to generate the syndynes/synchrones.
 
-Zero-ejection velocity syndynes are generated with the `~sbpy.dynamics.syndynes.SynGenerator` class.  The class requires a dust source described by a `~sbpy.dynamics.state.State` object, :math:`\beta` values, and particle ages from which to generate the syndynes.
-
-First, define the source of the syndynes, a comet at 2 au from the Sun:
+First, define the dust source, a comet at 2 au from the Sun:
 
 .. doctest::
 
@@ -243,7 +241,7 @@ First, define the source of the syndynes, a comet at 2 au from the Sun:
    >>> t = Time("2023-12-08")
    >>> comet = State(r, v, t)
 
-Next, initialize the syndyne object:
+Next, initialize the syn-generator object.  Here, we request dust with :math:`\beta` values from 0 to 1, generated over a 100-day period:
 
 .. doctest-requires:: scipy
 
@@ -268,7 +266,28 @@ The computed particle positions are saved in the :attr:`~sbpy.dynamics.syndynes.
    >>> print(len(dust.particles))
    104
 
-Get the results with the :func:`~sbpy.dynamics.syndynes.SynGenerator.syndynes` method, which returns a list-like collection of `~sbpy.dynamics.syndynes.Syndyne` objects.  `Syndyne` objects are specialized `State` objects.  For example, we can compute the linear distance from the comet to farthest particle in each syndyne:
+
+Syndynes
+--------
+
+Zero-ejection velocity syndynes are generated with the `~sbpy.dynamics.syndynes.SynGenerator` class, and this example uses the generator above.  Get the syndynes with the :func:`~sbpy.dynamics.syndynes.SynGenerator.syndynes` method, which returns `~sbpy.dynamics.syndynes.Syndynes` a list-like collection of `~sbpy.dynamics.syndynes.Syndyne` objects (note the use of plural and singular forms).
+
+.. doctest-requires:: scipy
+
+   >>> syndynes = dust.syndynes()
+   >>> syndynes
+   <Syndynes: betas=[1.   0.1  0.01 0.  ]>
+   >>> syndynes[0]
+   <Syndyne (<ArbitraryFrame Frame>):
+    r
+     [[ 2.99195741e+08  0.00000000e+00  0.00000000e+00]
+    [ 2.99284224e+08 -2.04410384e+03  0.00000000e+00]
+    [ 2.99549029e+08 -1.63231106e+04  0.00000000e+00]
+    [ 2.99988249e+08 -5.49241301e+04  0.00000000e+00]
+    [ 3.00598735e+08 -1.29642418e+05  0.00000000e+00]
+   ...
+
+`Syndyne` objects are specialized `State` objects.  For example, we can compute the linear distance from the comet to farthest particle in each syndyne:
 
 .. doctest-requires:: scipy
 
@@ -280,7 +299,7 @@ Get the results with the :func:`~sbpy.dynamics.syndynes.SynGenerator.syndynes` m
    0.01, 0.003 AU
    0.0, 0.000 AU
 
-Individual syndynes may be retrieved with the :func:`~sbpy.dynamics.syndynes.SynGenerator.syndyne` method and a syndyne index.  The index for the syndyne matches the index of the `betas` array, i.e., to get the :math:`\beta=0.1` syndyne from our example:
+Individual syndynes may also be retrieved directly from the `~sbpy.dynamics.syndynes.SynGenerator` object using the :func:`~sbpy.dynamics.syndynes.SynGenerator.syndyne` method and a syndyne index.  The index for the syndyne matches the index of the `betas` array, i.e., to get the :math:`\beta=0.1` syndyne from our example:
 
 .. doctest-requires:: scipy
 
@@ -417,12 +436,12 @@ In this example, we compute the syndynes of a comet orbiting β Pic (1.8 solar m
    >>> betapic_dust = SynGenerator(comet, [1, 0], [0, 100] * u.d, solver=solver)
 
 
-Plotting syndynes and synchrones
---------------------------------
+Plotting syndynes, synchrones, and orbits
+-----------------------------------------
 
-Generally, we are interested in visualizing the syndynes and synchrones for an observer.  `Syndynes`, `Synchrones`, and `SourceOrbit` have ``plot()`` methods to assist with this.  They can plot the coordinates relative to the comet with a simple tangent plane projection, or projected onto the image plane with an `astropy.wcs.WCS` object.
+Generally, we are interested in visualizing the syndynes and synchrones for an observer.  `Syndynes`, `Synchrones`, and `SourceOrbit` objects have ``plot()`` methods to assist with this.  They can plot the coordinates relative to the comet with a simple tangent plane projection, or projected onto the image plane with an `astropy.wcs.WCS` object.
 
-Here is a simple example that plots the syndynes and synchrones from above as offsets from the comet's coordinates:
+Here is an example that plots the syndynes and synchrones from above as offsets from the comet's coordinates:
 
 .. doctest-requires:: scipy,matplotlib
 

From 756a62ce2e486133f4d294548394af77b0890a37 Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Thu, 7 May 2026 09:40:11 -0400
Subject: [PATCH 27/34] simplify plot examples

---
 docs/sbpy/dynamics.rst | 30 ++++++++++++++++++++----------
 1 file changed, 20 insertions(+), 10 deletions(-)

diff --git a/docs/sbpy/dynamics.rst b/docs/sbpy/dynamics.rst
index f9403b41..baa645ce 100644
--- a/docs/sbpy/dynamics.rst
+++ b/docs/sbpy/dynamics.rst
@@ -450,14 +450,17 @@ Here is an example that plots the syndynes and synchrones from above as offsets
    >>> fig, ax = plt.subplots()
    >>>
    >>> # plot all but the last (beta=0) syndyne
-   >>> dust.syndynes()[:-1].plot(ax)
+   >>> syndynes = dust.syndynes()
+   >>> syndynes[:-1].plot(ax)
    >>>
    >>> # plot every 5th synchrone
-   >>> dust.synchrones()[4::5].plot(ax)
+   >>> synchrones = dust.synchrones()
+   >>> synchrones[4::5].plot(ax, ls="--")
    >>>
    >>> # plot the orbit
    >>> dt = np.linspace(-2, 2) * u.d
-   >>> dust.source_orbit(dt).plot(ax, color="k", ls=":", label="Orbit")
+   >>> orbit = dust.source_orbit(dt)
+   >>> orbit.plot(ax, color="k", ls=":", label="Orbit")
    >>>
    >>> ax.invert_xaxis()
    >>> plt.setp(ax,
@@ -499,14 +502,17 @@ Here is an example that plots the syndynes and synchrones from above as offsets
    fig, ax = plt.subplots()
 
    # plot all but the last (beta=0) syndyne
-   dust.syndynes()[:-1].plot(ax)
+   syndynes = dust.syndynes()
+   syndynes[:-1].plot(ax)
 
    # plot every 5th synchrone
-   dust.synchrones()[4::5].plot(ax, ls="--", lw=1)
+   synchrones = dust.synchrones()
+   synchrones[4::5].plot(ax, ls="--", lw=1)
 
    # plot the orbit
    dt = np.linspace(-2, 2) * u.d
-   dust.source_orbit(dt).plot(ax, color="k", ls=":", label="Orbit")
+   orbit = dust.source_orbit(dt)
+   orbit.plot(ax, color="k", ls=":", label="Orbit")
 
    ax.invert_xaxis()
    plt.setp(ax,
@@ -570,11 +576,13 @@ The following example compares syndynes to a Spitzer Space Telesocpe image of co
    >>> ylim = ax.get_ylim()
    >>> 
    >>> # plot syndynes
-   >>> dust.syndynes().plot(ax, wcs=wcs)
+   >>> syndynes = dust.syndynes()
+   >>> syndynes.plot(ax, wcs=wcs)
    >>> 
    >>> # plot the orbit
    >>> dt = np.linspace(-1, 1) * u.d
-   >>> dust.source_orbit(dt).plot(ax, wcs=wcs, color="tab:cyan", lw=1, label="Orbit")
+   >>> orbit = dust.source_orbit(dt)
+   >>> orbit.plot(ax, wcs=wcs, color="tab:cyan", lw=1, label="Orbit")
    >>> 
    >>> plt.setp(ax, xlim=xlim, ylim=ylim)  # doctest: +SKIP
    >>> plt.legend()  # doctest: +SKIP
@@ -640,11 +648,13 @@ The following example compares syndynes to a Spitzer Space Telesocpe image of co
    ylim = ax.get_ylim()
 
    # plot syndynes
-   dust.syndynes().plot(ax, wcs=wcs)
+   syndynes = dust.syndynes()
+   syndynes.plot(ax, wcs=wcs)
 
    # plot the orbit
    dt = np.linspace(-1, 1) * u.d
-   dust.source_orbit(dt).plot(ax, wcs=wcs, color="tab:cyan", lw=1, label="Orbit")
+   orbit = dust.source_orbit(dt)
+   orbit.plot(ax, wcs=wcs, color="tab:cyan", lw=1, label="Orbit")
 
    plt.setp(ax, xlim=xlim, ylim=ylim)
    plt.legend()

From ccc5ebdf222a4b99e52e842036e4036f9d6d96c6 Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Thu, 7 May 2026 13:30:07 -0400
Subject: [PATCH 28/34] Move Syndyne.plot example up

---
 docs/sbpy/dynamics.rst | 89 ++++++++++++++++++++++++++++--------------
 1 file changed, 59 insertions(+), 30 deletions(-)

diff --git a/docs/sbpy/dynamics.rst b/docs/sbpy/dynamics.rst
index baa645ce..645f3415 100644
--- a/docs/sbpy/dynamics.rst
+++ b/docs/sbpy/dynamics.rst
@@ -524,8 +524,66 @@ Here is an example that plots the syndynes and synchrones from above as offsets
    plt.legend()
    plt.tight_layout()
 
+For more complex plot logic, e.g., specific line colors and styles, use the plot methods of the individual syndynes/synchrones:
 
-The following example compares syndynes to a Spitzer Space Telesocpe image of comet 48P/Johnson (`Reach et al. 2007 <https://scixplorer.org/abs/2007Icar..191..298R/abstract>`_).  The FITS world coordinate system is used to account for the image orientation and scale.  To precisely align the syndynes with the comet nucleus, we update the world coordinate system to use our calculated comet coordinates.
+.. doctest-requires:: scipy,astroquery,matplotlib
+.. doctest-remote-data::
+
+   >>> ls = ["-", "--", "-."]
+   >>> syndynes = dust.syndynes()
+   >>> for i in range(3):
+   ...     syndynes[i].plot(ax, color="k", ls=ls[i])
+
+.. plot::
+   :show-source-link:
+
+   import numpy as np
+   import matplotlib.pyplot as plt
+
+   import astropy.units as u
+   from astropy.time import Time
+   from sbpy.dynamics import State, SynGenerator
+
+   # define the dust source
+   r = [2, 0, 0] * u.au
+   v = [0, 30, 0] * u.km / u.s
+   t = Time("2023-12-08")
+   frame = "heliocentriceclipticiau76"
+   comet = State(r, v, t, frame=frame)
+
+   # define particle parameters
+   betas = [1, 0.1, 0.01, 0]
+   ages = np.linspace(0, 100, 25) * u.day
+
+   # observe it from a fixed location in the Solar System
+   observer = State(
+       r=[0, 2, 2] * u.au,
+       v=[0, 0, 0] * u.km / u.s,
+       t=comet.t,
+       frame="icrs",
+   )
+
+   # define the dust generator
+   dust = SynGenerator(comet, betas, ages, observer=observer)
+
+   # plot syndynes
+   fig, ax = plt.subplots()
+
+   ls = ["-", "--", "-."]
+   syndynes = dust.syndynes()
+   for i in range(3):
+      syndynes[i].plot(ax, color="k", ls=ls[i])
+      ax.invert_xaxis()
+
+   ax.set(ax,
+          xlim=[100, -10],
+          ylim=[-10, 100],
+          xlabel="$\\Delta$RA (arcsec)",
+          ylabel="$\\Delta$Dec (arcsec)",
+   )
+   fig.tight_layout()
+
+The following complete example compares syndynes to a Spitzer Space Telesocpe image of comet 48P/Johnson (`Reach et al. 2007 <https://scixplorer.org/abs/2007Icar..191..298R/abstract>`_).  The FITS world coordinate system is used to account for the image orientation and scale.  To precisely align the syndynes with the comet nucleus, we update the world coordinate system to use our calculated comet coordinates.
 
 .. note::
    The `sbpy` testing suite shows that arcsecond-level accuracy is possible, but this is generally not enough for direct comparison to typical images of comets, which need sub-arcsecond alignment.  The accuracy of the coordinates object depends on the the comet and observer states, but also on whether or not light travel time is accounted for, and the accuracy of the orbit integrator.
@@ -660,35 +718,6 @@ The following example compares syndynes to a Spitzer Space Telesocpe image of co
    plt.legend()
    plt.tight_layout()
 
-For more complex plot logic, e.g., to use specific line colors and styles, we can use the plot methods of the individual syndynes/synchrones:
-
-.. doctest-requires:: scipy,astroquery,matplotlib
-.. doctest-remote-data::
-
-   >>> ls = ["-", "--", "-."]
-   >>> syndynes = dust.syndynes()
-   >>> for i in range(3):
-   ...     syndynes[i].plot(ax, color="k", ls=ls[i])
-
-.. plot::
-   :context:
-
-   fig, ax = plt.subplots()
-
-   ls = ["-", "--", "-."]
-   syndynes = dust.syndynes()
-   for i in range(3):
-      syndynes[i].plot(ax, color="k", ls=ls[i])
-      ax.invert_xaxis()
-
-   plt.setp(ax,
-            xlim=[100, -10],
-            ylim=[-10, 100],
-            xlabel="$\\Delta$RA (arcsec)",
-            ylabel="$\\Delta$Dec (arcsec)",
-   )
-   plt.tight_layout()
-
 
 Reference/API
 =============

From 6a929e04d071565e859cf93ed23ca78dd0104710 Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Thu, 7 May 2026 13:37:23 -0400
Subject: [PATCH 29/34] Fix ax.set line

---
 docs/sbpy/dynamics.rst | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/docs/sbpy/dynamics.rst b/docs/sbpy/dynamics.rst
index 645f3415..edf4c8bf 100644
--- a/docs/sbpy/dynamics.rst
+++ b/docs/sbpy/dynamics.rst
@@ -575,11 +575,11 @@ For more complex plot logic, e.g., specific line colors and styles, use the plot
       syndynes[i].plot(ax, color="k", ls=ls[i])
       ax.invert_xaxis()
 
-   ax.set(ax,
-          xlim=[100, -10],
-          ylim=[-10, 100],
-          xlabel="$\\Delta$RA (arcsec)",
-          ylabel="$\\Delta$Dec (arcsec)",
+   ax.set(
+      xlim=[100, -10],
+      ylim=[-10, 100],
+      xlabel="$\\Delta$RA (arcsec)",
+      ylabel="$\\Delta$Dec (arcsec)",
    )
    fig.tight_layout()
 

From 5ef0ccde7b3b33f8b317bc76090752d1786d0ad3 Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Thu, 7 May 2026 13:43:54 -0400
Subject: [PATCH 30/34] Give a resolution to the note on syndyne accuracy

---
 docs/sbpy/dynamics.rst | 40 ++++++++++++++++++++--------------------
 1 file changed, 20 insertions(+), 20 deletions(-)

diff --git a/docs/sbpy/dynamics.rst b/docs/sbpy/dynamics.rst
index edf4c8bf..8a0af8f5 100644
--- a/docs/sbpy/dynamics.rst
+++ b/docs/sbpy/dynamics.rst
@@ -463,14 +463,14 @@ Here is an example that plots the syndynes and synchrones from above as offsets
    >>> orbit.plot(ax, color="k", ls=":", label="Orbit")
    >>>
    >>> ax.invert_xaxis()
-   >>> plt.setp(ax,
-   ...          xlim=[100, -10],
-   ...          ylim=[-10, 100],
-   ...          xlabel="$\\Delta$RA (arcsec)",
-   ...          ylabel="$\\Delta$Dec (arcsec)",
+   >>> ax.set(
+   ...     xlim=[100, -10],
+   ...     ylim=[-10, 100],
+   ...     xlabel="$\\Delta$RA (arcsec)",
+   ...     ylabel="$\\Delta$Dec (arcsec)",
    ... ) # doctest: +SKIP
-   >>> plt.legend() # doctest: +SKIP
-   >>> plt.tight_layout()
+   >>> ax.legend() # doctest: +SKIP
+   >>> fig.tight_layout()
 
 .. plot::
    :context:
@@ -515,14 +515,14 @@ Here is an example that plots the syndynes and synchrones from above as offsets
    orbit.plot(ax, color="k", ls=":", label="Orbit")
 
    ax.invert_xaxis()
-   plt.setp(ax,
-            xlim=[100, -10],
-            ylim=[-10, 100],
-            xlabel="$\\Delta$RA (arcsec)",
-            ylabel="$\\Delta$Dec (arcsec)",
+   ax.set(
+      xlim=[100, -10],
+      ylim=[-10, 100],
+      xlabel="$\\Delta$RA (arcsec)",
+      ylabel="$\\Delta$Dec (arcsec)",
    )
-   plt.legend()
-   plt.tight_layout()
+   ax.legend()
+   fig.tight_layout()
 
 For more complex plot logic, e.g., specific line colors and styles, use the plot methods of the individual syndynes/synchrones:
 
@@ -586,7 +586,7 @@ For more complex plot logic, e.g., specific line colors and styles, use the plot
 The following complete example compares syndynes to a Spitzer Space Telesocpe image of comet 48P/Johnson (`Reach et al. 2007 <https://scixplorer.org/abs/2007Icar..191..298R/abstract>`_).  The FITS world coordinate system is used to account for the image orientation and scale.  To precisely align the syndynes with the comet nucleus, we update the world coordinate system to use our calculated comet coordinates.
 
 .. note::
-   The `sbpy` testing suite shows that arcsecond-level accuracy is possible, but this is generally not enough for direct comparison to typical images of comets, which need sub-arcsecond alignment.  The accuracy of the coordinates object depends on the the comet and observer states, but also on whether or not light travel time is accounted for, and the accuracy of the orbit integrator.
+   The `sbpy` testing suite shows that arcsecond-level accuracy is possible, but this is generally not enough for direct comparison to typical images of comets, which need sub-arcsecond alignment.  The accuracy of the coordinates object depends on the the comet and observer states, but also on whether or not light travel time is accounted for, and the accuracy of the orbit integrator.  In the following example, the world coordinate system is edited to manually align the image and syndynes.
 
 .. doctest-requires:: scipy,astroquery,matplotlib
 .. doctest-remote-data::
@@ -642,8 +642,8 @@ The following complete example compares syndynes to a Spitzer Space Telesocpe im
    >>> orbit = dust.source_orbit(dt)
    >>> orbit.plot(ax, wcs=wcs, color="tab:cyan", lw=1, label="Orbit")
    >>> 
-   >>> plt.setp(ax, xlim=xlim, ylim=ylim)  # doctest: +SKIP
-   >>> plt.legend()  # doctest: +SKIP
+   >>> ax.set(xlim=xlim, ylim=ylim)  # doctest: +SKIP
+   >>> ax.legend()  # doctest: +SKIP
 
 .. plot::
    :show-source-link:
@@ -714,9 +714,9 @@ The following complete example compares syndynes to a Spitzer Space Telesocpe im
    orbit = dust.source_orbit(dt)
    orbit.plot(ax, wcs=wcs, color="tab:cyan", lw=1, label="Orbit")
 
-   plt.setp(ax, xlim=xlim, ylim=ylim)
-   plt.legend()
-   plt.tight_layout()
+   ax.set(xlim=xlim, ylim=ylim)
+   ax.legend()
+   fig.tight_layout()
 
 
 Reference/API

From 7ab9960f7b3af4d4fab419f50e78659c45906ad8 Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Thu, 7 May 2026 13:46:54 -0400
Subject: [PATCH 31/34] Fix typo.

---
 docs/sbpy/dynamics.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/sbpy/dynamics.rst b/docs/sbpy/dynamics.rst
index 8a0af8f5..789d3781 100644
--- a/docs/sbpy/dynamics.rst
+++ b/docs/sbpy/dynamics.rst
@@ -583,7 +583,7 @@ For more complex plot logic, e.g., specific line colors and styles, use the plot
    )
    fig.tight_layout()
 
-The following complete example compares syndynes to a Spitzer Space Telesocpe image of comet 48P/Johnson (`Reach et al. 2007 <https://scixplorer.org/abs/2007Icar..191..298R/abstract>`_).  The FITS world coordinate system is used to account for the image orientation and scale.  To precisely align the syndynes with the comet nucleus, we update the world coordinate system to use our calculated comet coordinates.
+The following complete example compares syndynes to a Spitzer Space Telescope image of comet 48P/Johnson (`Reach et al. 2007 <https://scixplorer.org/abs/2007Icar..191..298R/abstract>`_).  The FITS world coordinate system is used to account for the image orientation and scale.  To precisely align the syndynes with the comet nucleus, we update the world coordinate system to use our calculated comet coordinates.
 
 .. note::
    The `sbpy` testing suite shows that arcsecond-level accuracy is possible, but this is generally not enough for direct comparison to typical images of comets, which need sub-arcsecond alignment.  The accuracy of the coordinates object depends on the the comet and observer states, but also on whether or not light travel time is accounted for, and the accuracy of the orbit integrator.  In the following example, the world coordinate system is edited to manually align the image and syndynes.

From 5e48a48244d2e5aac24f8c968ec05c40dca20c8b Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Thu, 7 May 2026 13:54:33 -0400
Subject: [PATCH 32/34] Fix doctests

---
 docs/sbpy/dynamics.rst | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/sbpy/dynamics.rst b/docs/sbpy/dynamics.rst
index 789d3781..c4305656 100644
--- a/docs/sbpy/dynamics.rst
+++ b/docs/sbpy/dynamics.rst
@@ -277,7 +277,7 @@ Zero-ejection velocity syndynes are generated with the `~sbpy.dynamics.syndynes.
    >>> syndynes = dust.syndynes()
    >>> syndynes
    <Syndynes: betas=[1.   0.1  0.01 0.  ]>
-   >>> syndynes[0]
+   >>> syndynes[0]  # doctest: +FLOAT_CMP
    <Syndyne (<ArbitraryFrame Frame>):
     r
      [[ 2.99195741e+08  0.00000000e+00  0.00000000e+00]
@@ -348,7 +348,7 @@ This method produced a `SourceOrbit` object that is functionally the same as the
 .. doctest-requires:: scipy
 
    >>> type(orbit)
-   sbpy.dynamics.syndynes.SourceOrbit
+   <class 'sbpy.dynamics.syndynes.SourceOrbit'>
    >>> orbit.ages
    <Quantity [ 2., -0., -2.] d>
    >>> np.linalg.norm(orbit[0].r - orbit[1].r)

From 0ec12f8c04d608b39eae93735e85b76d78419e97 Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Thu, 7 May 2026 13:58:54 -0400
Subject: [PATCH 33/34] doctest: Another missing float comparison

---
 docs/sbpy/dynamics.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/sbpy/dynamics.rst b/docs/sbpy/dynamics.rst
index c4305656..cef2bef0 100644
--- a/docs/sbpy/dynamics.rst
+++ b/docs/sbpy/dynamics.rst
@@ -351,7 +351,7 @@ This method produced a `SourceOrbit` object that is functionally the same as the
    <class 'sbpy.dynamics.syndynes.SourceOrbit'>
    >>> orbit.ages
    <Quantity [ 2., -0., -2.] d>
-   >>> np.linalg.norm(orbit[0].r - orbit[1].r)
+   >>> np.linalg.norm(orbit[0].r - orbit[1].r)  # doctest: +FLOAT_CMP
    <Quantity 5183919.42549707 km>
 
 

From 408c2589307ac641eb22848405f47193564064ce Mon Sep 17 00:00:00 2001
From: "Michael S. P. Kelley" <msk@astro.umd.edu>
Date: Thu, 7 May 2026 14:06:04 -0400
Subject: [PATCH 34/34] Add a new sub-section for clarity

---
 docs/sbpy/dynamics.rst | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/docs/sbpy/dynamics.rst b/docs/sbpy/dynamics.rst
index cef2bef0..92b5223a 100644
--- a/docs/sbpy/dynamics.rst
+++ b/docs/sbpy/dynamics.rst
@@ -224,9 +224,12 @@ Dust syndynes and synchrones
 
 Syndynes are lines in space connecting particles that are experiencing the same forces.  A syndyne is parameterized by :math:`\beta`, the ratio of the force from solar radiation to the force from solar gravity, :math:`F_r / F_g`, and age (or time of release).  Thus, all particles in a syndyne have a constant :math:`\beta` but variable age.  Similarly, synchrones are lines of constant particle age, but variable :math:`\beta`.
 
+Start with a `SynGenerator`
+---------------------------
+
 sbpy has the `sbpy.dynamics.syndynes.SynGenerator` class to assist the user in the generation of syndynes, synchrones, and points along a source object's orbit.  Use the `~sbpy.dynamics.syndynes.SynGenerator`'s :func:`~sbpy.dynamics.syndynes.SynGenerator.syndynes`, :func:`~sbpy.dynamics.syndynes.SynGenerator.synchrones`, and :func:`~sbpy.dynamices.syndynes.SynGenerator.source_orbit` methods to produce `~sbpy.dynamics.syndynes.Syndynes`, `~sbpy.dynamics.syndynes.Synchrones`, and `~sbpy.dynamics.syndynes.SourceOrbit` objects, described in further detail below.
 
-First, setup a `~sbpy.dynamics.syndynes.SynGenerator` instance.  This requires a dust source described by a `~sbpy.dynamics.state.State` object, :math:`\beta` values, and particle ages from which to generate the syndynes/synchrones.
+Setting up a `~sbpy.dynamics.syndynes.SynGenerator` instance requires a dust source, and dust :math:`\beta` values and ages from which to generate the syndynes/synchrones.
 
 First, define the dust source, a comet at 2 au from the Sun: