diff --git a/CHANGES.rst b/CHANGES.rst
index f5bc53d44..aa7209769 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -8,12 +8,55 @@ 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]
+
+- 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]
 
 
+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. [#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. [#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
+---------
+
+sbpy.dynamics
+^^^^^^^^^^^^^
+
+- Added missing `sbpy.dynamics.state.State.to_ephem`. [#434]
+
+
 v0.6.0 (2025-12-02)
 ===================
 
diff --git a/docs/sbpy/dynamics.rst b/docs/sbpy/dynamics.rst
index 60549fc2b..92b5223a8 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:
 
@@ -168,6 +170,25 @@ 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
+.. doctest-remote-data::
+
+   >>> eph = Ephem.from_horizons(
+   ...     "48P",
+   ...     id_type="designation",
+   ...     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")
+
+
 Dynamical integrators
 =====================
 
@@ -203,13 +224,14 @@ 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`
+---------------------------
 
-Syndynes
---------
+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.
 
-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.
+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 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::
 
@@ -222,7 +244,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
 
@@ -247,7 +269,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]  # doctest: +FLOAT_CMP
+   <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
 
@@ -259,7 +302,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
 
@@ -293,6 +336,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)
+   <class 'sbpy.dynamics.syndynes.SourceOrbit'>
+   >>> orbit.ages
+   <Quantity [ 2., -0., -2.] d>
+   >>> np.linalg.norm(orbit[0].r - orbit[1].r)  # doctest: +FLOAT_CMP
+   <Quantity 5183919.42549707 km>
+
+
 Adding ejection velocity
 ------------------------
 
@@ -355,17 +420,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) * u.d
-   >>> orbit, coords = dust.source_orbit(dt)
-
-
 Using other dynamical models
 ----------------------------
 
@@ -385,52 +439,44 @@ 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 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 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 an 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
-   >>> 
-   >>> 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()
-   >>> 
+   >>>
    >>> # plot all but the last (beta=0) syndyne
-   >>> for syndyne in dust.syndynes()[:-1]:
-   ...     plot(ax, syndyne.coords, label=f"$\\beta={syndyne.beta:.2g}$")
-   >>> 
+   >>> syndynes = dust.syndynes()
+   >>> 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
+   >>> synchrones = dust.synchrones()
+   >>> synchrones[4::5].plot(ax, ls="--")
+   >>>
+   >>> # plot the orbit
    >>> dt = np.linspace(-2, 2) * u.d
-   >>> states, coords = dust.source_orbit(dt)
-   >>> plot(ax, coords, color="k", ls=":", label="Orbit")
+   >>> orbit = dust.source_orbit(dt)
+   >>> 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:
 
    import numpy as np
    import matplotlib.pyplot as plt
@@ -455,43 +501,225 @@ Generally, we are interested in plotting syndynes and synchrones on an image of
    )
    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
+   syndynes = dust.syndynes()
+   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
+   synchrones = dust.synchrones()
+   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")
+   orbit = dust.source_orbit(dt)
+   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:
+
+.. 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(
+      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 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.
+
+.. 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"])
+   >>>
+   >>> # 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
+   >>>
+   >>> # 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])
+   >>> 
+   >>> # 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")  # docte
+   >>> 
+   >>> # save xlim and ylim for later
+   >>> xlim = ax.get_xlim()
+   >>> ylim = ax.get_ylim()
+   >>> 
+   >>> # plot syndynes
+   >>> syndynes = dust.syndynes()
+   >>> syndynes.plot(ax, wcs=wcs)
+   >>> 
+   >>> # plot the orbit
+   >>> dt = np.linspace(-1, 1) * u.d
+   >>> orbit = dust.source_orbit(dt)
+   >>> orbit.plot(ax, wcs=wcs, color="tab:cyan", lw=1, label="Orbit")
+   >>> 
+   >>> ax.set(xlim=xlim, ylim=ylim)  # doctest: +SKIP
+   >>> ax.legend()  # doctest: +SKIP
+
+.. 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=49.1, vmax=49.5, cmap="gray_r")
+
+   # save xlim and ylim for later
+   xlim = ax.get_xlim()
+   ylim = ax.get_ylim()
+
+   # plot syndynes
+   syndynes = dust.syndynes()
+   syndynes.plot(ax, wcs=wcs)
+
+   # plot the orbit
+   dt = np.linspace(-1, 1) * u.d
+   orbit = dust.source_orbit(dt)
+   orbit.plot(ax, wcs=wcs, color="tab:cyan", lw=1, label="Orbit")
+
+   ax.set(xlim=xlim, ylim=ylim)
+   ax.legend()
+   fig.tight_layout()
 
 
 Reference/API
diff --git a/sbpy/data/__init__.py b/sbpy/data/__init__.py
index 0a5a471b2..310a84035 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)",
diff --git a/sbpy/dynamics/state.py b/sbpy/dynamics/state.py
index a41aad8ac..a1c3fdae1 100644
--- a/sbpy/dynamics/state.py
+++ b/sbpy/dynamics/state.py
@@ -294,6 +294,112 @@ def to_skycoord(self) -> SkyCoord:
         else:
             return SkyCoord(self._data, obstime=self.t, representation_type="cartesian")
 
+    def to_ephem(
+        self, observer: StateType | None = None, coords: SkyCoord | None = None
+    ) -> Ephem:
+        r"""Convert to an sbpy ephemeris object.
+
+
+        Parameters
+        ----------
+        observer : `State`, optional
+            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
+        -------
+        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": observer.r,
+                "v": observer.v,
+                "t": observer.t,
+                "frame": observer.frame,
+            }
+
+            if coords is None:
+                coords = observer.observe(self)
+
+            # use SkyCoord's to_table() method, which will account for when
+            # RA/Dec vs lon/lat are used.
+            tab = coords.to_table()
+
+            # 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",
+                "pm_lat": "lat_rate",
+                "distance": "delta",
+                "radial_velocity": "deltadot",
+            }
+            for col in tab.colnames:
+                field = skycoord_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.
 
diff --git a/sbpy/dynamics/syndynes.py b/sbpy/dynamics/syndynes.py
index a198fbcbf..52ffefcfe 100644
--- a/sbpy/dynamics/syndynes.py
+++ b/sbpy/dynamics/syndynes.py
@@ -7,7 +7,10 @@
 
 """
 
+from __future__ import annotations
+
 __all__ = [
+    "SourceOrbit",
     "SynGenerator",
     "SynStates",
     "SynCollection",
@@ -22,17 +25,26 @@
 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.core 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):
@@ -94,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
@@ -118,59 +133,93 @@ def to_ephem(self) -> Ephem:
 
         """
 
-        data: dict = {}
-        data["beta_rad"] = self.betas
-        data["age"] = self.ages
+        eph = super().to_ephem(observer=self.observer, coords=self.coords)
 
-        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(
+        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)
+
+            # 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)
+            x, y = wcs.world_to_pixel(_unmasked(self.coords))
+
+        ax.plot(x, y, **kwargs)
 
 
 class Syndyne(SynStates):
@@ -179,7 +228,7 @@ class Syndyne(SynStates):
 
     Parameters
     ----------
-    source : State
+    source : `State`
         The source of the syndyne dust.
 
     beta : float
@@ -197,6 +246,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 +297,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 +329,52 @@ 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.
+
+    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``.
+
+    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,
+        dt: u.Quantity,
+        r: u.Quantity,
+        v: u.Quantity,
+        t: Time,
+        observer: Optional[State] = None,
+    ) -> None:
+        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:
+        """Epoch of each orbit point."""
+        return self.source.t + self.dt
+
+
 class SynCollection:
     """Immutable collection of syndynes or synchrones.
 
@@ -310,8 +411,16 @@ 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:
+        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])
 
     def to_ephem(self) -> Ephem:
         """Convert to an sbpy ephemeris object.
@@ -323,12 +432,65 @@ 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(
+        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):
@@ -614,13 +776,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] = []
@@ -629,8 +787,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)
diff --git a/sbpy/dynamics/tests/test_state.py b/sbpy/dynamics/tests/test_state.py
index 085c86477..d21069c79 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,37 @@ def test_from_ephem(self):
             del incomplete.table[k]
             with pytest.raises(ValueError):
                 State.from_ephem(incomplete)
+
+    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
+        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)
+
+        # 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)
diff --git a/sbpy/dynamics/tests/test_syndynes.py b/sbpy/dynamics/tests/test_syndynes.py
index 8e239d484..c19de24a4 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):
+    """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,
+        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))
+        initial = fake_state(comet, np.zeros(10))
+        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)
+        initial = fake_state(comet, np.zeros(10))
+        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):
@@ -212,6 +304,30 @@ 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 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):
@@ -281,96 +397,170 @@ 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)
 
+    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(
-        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()
 
-def test_syndynes():
-    comet = random_state()
-    betas = np.linspace(0, 1, 3)
-    ages = np.arange(10) * u.day
+        # create from list
+        syndynes = Syndynes(items)
 
-    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 syndynes object
+        syndynes2 = Syndynes(syndynes)
+        assert np.all(syndynes2[0].r == syndynes[0].r)
+
+        # 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])
+        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])
+        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:
@@ -485,7 +675,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()
@@ -493,8 +683,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)
diff --git a/sbpy/utils/core.py b/sbpy/utils/core.py
index caec9205d..458d5d2f5 100644
--- a/sbpy/utils/core.py
+++ b/sbpy/utils/core.py
@@ -10,6 +10,16 @@
 
 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 272b0cada..d316aa83d 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()