Dynamics (sbpy.dynamics)¶
sbpy has the capability to describe dynamical states (position and velocity vectors), and to integrate test particle orbits. These are primarily intended to support the generation of Dust syndynes and synchrones.
Dynamical state¶
State objects encapsulate the position and velocity of an object at a given time. Create a State for a comet at \(x=2\) au, moving along the y-axis at a speed of 30 km/s:
>>> from astropy.time import Time
>>> import astropy.units as u
>>> from sbpy.dynamics import State
>>>
>>> r = [2, 0, 0] * u.au
>>> v = [0, 30, 0] * u.km / u.s
>>> t = Time("2023-12-08")
>>> comet = State(r, v, t)
State objects may also represent an array of objects:
>>> r = ([2, 0, 0], [0, 2, 0]) * u.au
>>> v = ([0, 30, 0], [30, 0, 0]) * u.km / u.s
>>> t = Time(["2023-12-08", "2023-12-09"])
>>> comets = State(r, v, t)
>>> len(comets)
2
The r, v, and t attributes hold the position, velocity, and time for the object(s). The first index iterates over the object, the second iterates over the x-, y-, and z-axes.
>>> comets.r.shape
(2, 3)
>>> comets.r[0]
<Quantity [2., 0., 0.] AU>
Or, index the State object itself:
>>> comets[0].r # equivalent to comets.r[0]
<Quantity [2., 0., 0.] AU>
Convert to/from Ephem and SkyCoord¶
State objects may be initialized from ephemeris objects (Ephem), provided they contain time, and 3D position and velocity:
>>> from sbpy.data import Ephem
>>> eph = Ephem.from_horizons("9P",
... epochs=Time("2005-07-04"),
... id_type="designation",
... closest_apparition=True)
>>> tempel1 = State.from_ephem(eph)
And State may be converted to an Ephem object:
>>> eph = tempel1.to_ephem()
astropy’s SkyCoord objects may also be used, assuming the time and 3D vectors are fully defined:
>>> from astropy.coordinates import SkyCoord
>>> coords = SkyCoord(ra="01:02:03h",
... dec="+04:05:06d",
... pm_ra_cosdec=100 * u.arcsec / u.hr,
... pm_dec=-10 * u.arcsec / u.hr,
... distance=1 * u.au,
... radial_velocity=15 * u.km / u.s,
... obstime=Time("2024-01-19"))
>>> state = State.from_skycoord(coords)
And back to a SkyCoord object:
>>> state.to_skycoord()
<SkyCoord (ICRS): (x, y, z) in AU
(0.96112414, 0.26676914, 0.07123631)
(v_x, v_y, v_z) in km / s
(9.16701924, 23.45244422, -0.94097856)>
Reference frames¶
Coordinate reference frames can be specified with the frame keyword argument during initialization. Most of astropy reference frames are supported (see astropy’s Built-in Frame Classes):
Note
When working with heliocentric ecliptic coordinates from JPL Horizons or NAIF SPICE, you may want to use the HeliocentricEclipticIAU76 reference frame.
>>> r = [2, 0, 0] * u.au
>>> v = [0, 30, 0] * u.km / u.s
>>> t = Time("2023-12-08")
>>> state = State(r, v, t, frame="heliocentriceclipticiau76")
>>> state
<State (<HeliocentricEclipticIAU76 Frame (obstime=2023-12-08 00:00:00.000)>):
r
[2. 0. 0.] AU
v
[ 0. 30. 0.] km / s
t
2023-12-08 00:00:00.000>
Use transform_to() to transform between reference frames:
>>> new_state = state.transform_to("icrs")
>>> new_state
<State (<ICRS Frame>):
r
[ 1.99191790e+00 -2.59236051e-03 -8.93633307e-04] AU
v
[8.11760173e-03 2.75129298e+01 1.19282350e+01] km / s
t
2023-12-08 00:00:00.000>
State lengths, subtraction, and observations¶
A few mathematical operations are possible. Get the magnitude of the heliocentric distance and velocity with abs:
>>> print(abs(comet))
(<Quantity 2. AU>, <Quantity 30. km / s>)
Get the Earth-comet state vector by subtraction:
>>> earth = State([0, 1, 0] * u.km, [30, 0, 0] * u.km / u.s, comet.t)
>>> print(comet - earth)
<State (<ArbitraryFrame Frame>):
r
[ 2.00000000e+00 -6.68458712e-09 0.00000000e+00] AU
v
[-30. 30. 0.] km / s
t
2023-12-08 00:00:00.000>
Or, get the sky coordinates of the comet as seen from the Earth:
>>> earth.observe(comet)
<SkyCoord (ArbitraryFrame): (lon, lat, distance) in (deg, deg, AU)
(359.99999981, 0., 2.)
(pm_lon, pm_lat, radial_velocity) in (mas / yr, mas / yr, km / s)
(6.52671946e+08, 0., -30.0000001)>
The result, a SkyCoord object, is expressed in the reference frame of the observer.
Dynamical integrators¶
A state object may be propagated to a new time using a dynamical integrator. Three integrators are defined. Use FreeExpansion for motion in free space, SolarGravity for orbits around the Sun, and SolarGravityAndRadiationPressure for orbits around the Sun considering radiation pressure.
>>> from sbpy.dynamics import SolarGravity
>>> state = State([1, 0, 0] * u.au, [0, 30, 0] * u.km / u.s, 0 * u.s)
>>> integrator = SolarGravity()
>>> t_final = 1 * u.year
>>> integrator.solve(state, t_final)
<State (<ArbitraryFrame Frame>):
r
[ 1.48146925e+08 -2.09358771e+07 0.00000000e+00] km
v
[ 4.137801 29.70907176 0. ] km / s
t
1.0 yr>
Other integrators may be defined. Use the above classes as templates or as base classes. For example, to simulate orbits around the Earth, update the _GM private attribute from SolarGravity:
>>> class EarthGravity(SolarGravity):
... _GM = 3.9860043543609598e5 # km3/s2
Dust syndynes and synchrones¶
Syndynes are lines in space connecting particles that are experiencing the same forces. A syndyne is parameterized by \(\beta\), the ratio of the force from solar radiation to the force from solar gravity, \(F_r / F_g\), and age (or time of release). Thus, all particles in a syndyne have a constant \(\beta\) but variable age. Similarly, synchrones are lines of constant particle age, but variable \(\beta\).
Syndynes¶
Zero-ejection velocity syndynes are generated with the SynGenerator class. The class requires a dust source described by a State object, \(\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:
>>> from astropy.time import Time
>>> import astropy.units as u
>>> from sbpy.dynamics import State
>>>
>>> r = [2, 0, 0] * u.au
>>> v = [0, 30, 0] * u.km / u.s
>>> t = Time("2023-12-08")
>>> comet = State(r, v, t)
Next, initialize the syndyne object:
>>> import numpy as np
>>> from sbpy.dynamics import SynGenerator
>>>
>>> betas = [1, 0.1, 0.01, 0]
>>> ages = np.linspace(0, 100, 26) * u.day
>>> dust = SynGenerator(comet, betas, ages)
>>> dust
<SynGenerator:
betas
[1. 0.1 0.01 0. ]
ages
[ 0. 4. 8. 12. 16. 20. 24. 28. 32. 36. 40. 44. 48. 52.
56. 60. 64. 68. 72. 76. 80. 84. 88. 92. 96. 100.] d>
The computed particle positions are saved in the particles attribute. For our example, the 4 \(\beta\)-values and the 26 ages produce 104 particles:
>>> print(len(dust.particles))
104
Get the results with the syndynes() method, which returns a list-like collection of 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:
>>> for syndyne in dust.syndynes():
... r, v = abs(syndyne - comet)
... print(syndyne.beta, "{:.3f}".format(r.max().to("au")), sep=", ")
1.0, 0.309 AU
0.1, 0.032 AU
0.01, 0.003 AU
0.0, 0.000 AU
Individual syndynes may be retrieved with the syndyne() method and a syndyne index. The index for the syndyne matches the index of the betas array, i.e., to get the \(\beta=0.1\) syndyne from our example:
>>> print(dust.betas)
[1. 0.1 0.01 0. ]
>>> syndyne = dust.syndyne(1)
>>> print(syndyne.beta)
0.1
Synchrones¶
Synchrones are also simulated with the SynGenerator class, but instead retrieved with the synchrone() and synchrones() methods, which return Synchrone objects, or collections thereof:
>>> synchrone = dust.synchrone(24)
>>> r, v = abs(synchrone)
>>> print(synchrone.age, "{:.3g}".format(r.max().to("au")), sep=", ")
96.0 d, 2.26 AU
To support generating synchrones at specific times, the at_epochs() method may be used to initialize a SynGenerator:
>>> dates = Time(["2023-10-01", "2023-11-01", "2023-12-01"])
>>> dust = SynGenerator.at_epochs(comet, betas, dates)
>>> synchrone = dust.synchrone(0)
>>> synchrone.epoch
<Time object: scale='utc' format='iso' value=2023-10-01 00:00:00.000>
Adding ejection velocity¶
An ejection velocity can be added to the dust particles. This is not implemented in sbpy, but is possible by sub-classing SynGenerator and overriding the initialize_states method. In the following example we add 10 m/s to the z-component of the velocity:
>>> class AlternativeSynGenerator(SynGenerator):
... def initialize_states(self):
... super().initialize_states()
... delta_v = State([0, 0, 0] * u.km, [0, 0, 0.01] * u.km / u.s, 0)
... self.initial_states = self.initial_states + delta_v
>>>
>>> dust_with_delta_v = AlternativeSynGenerator.at_epochs(comet, betas, dates)
>>> dust_with_delta_v.initial_states - dust.initial_states
<State (<ArbitraryFrame Frame>):
r
[[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]] km
v
[[0. 0. 0.01]
[0. 0. 0.01]
[0. 0. 0.01]] km / s
t
['2023-10-01 00:00:00.000' '2023-11-01 00:00:00.000'
'2023-12-01 00:00:00.000']>
Projecting onto the sky¶
Syndynes and synchrones may be projected onto the sky as seen by a telescope. This requires an observer. For precision work, the SynGenerator source object’s State should be defined in a heliocentric reference frame. Typically, the observer will observe in an equatorial reference frame, but your needs may vary. Here, we generate syndynes in the HeliocentricEclipticIAU76 coordinate frame, which is the same J2000 heliocentric ecliptic coordinate frame that JPL Horizons and the NAIF SPICE toolkit use. We then observe the results in the ICRS reference frame:
>>> r = [2, 0, 0] * u.au
>>> v = [0, 30, 0] * u.km / u.s
>>> t = Time("2023-12-08")
>>> comet = State(r, v, t, frame="heliocentriceclipticiau76")
>>> 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)
With the observer and coordinate frames defined, the generated syndynes and synchrones will have a coords attribute, which is an astropy.coordinates.SkyCoord object representing the sky positions of the test particles. Here, we print a sample of the coordinates:
>>> syndyne = dust.syndyne(0)
>>> print("\n".join(syndyne.coords[::5].to_string("hmsdms", precision=0)))
20h59m23s -35d18m48s
21h00m08s -35d12m50s
21h01m53s -34d55m44s
21h03m49s -34d30m22s
21h05m19s -34d00m41s
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 source_orbit() method:
>>> dt = np.linspace(-2, 2) * u.d
>>> orbit, coords = dust.source_orbit(dt)
Using other dynamical models¶
sbpy’s built-in models solve the equations of motion for dust grains given two-body dynamics. Users may provide their own models in order to, e.g., improve code performance, add planetary perturbations, model grain fragmentation, etc.. Use the SolarGravityAndRadiationPressure class as a template.
In this example, we compute the syndynes of a comet orbiting β Pic (1.8 solar masses) by sub-classing SolarGravityAndRadiationPressure and updating \(GM\), the mass of the star times the gravitational constant:
>>> import astropy.constants as const
>>> from sbpy.dynamics import SolarGravityAndRadiationPressure
>>>
>>> class BetaPicGravityAndRadiationPressure(SolarGravityAndRadiationPressure):
... _GM = 1.8 * SolarGravityAndRadiationPressure._GM
>>>
>>> solver = BetaPicGravityAndRadiationPressure()
>>> betapic_dust = SynGenerator(comet, [1, 0], [0, 100] * u.d, solver=solver)
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.
>>> 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}$")
>>>
>>> # 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()
(Source code, png, hires.png, pdf)
Reference/API¶
sbpy.dynamics.state Module¶
Classes¶
|
Coordinate frame with arbitrary space and time references. |
|
Abstract base class for dynamical state. |
|
Dynamical state. |
DynamicalModel solver failed. |
Class Inheritance Diagram¶
sbpy.dynamics.models Module¶
Classes¶
|
Super-class for dynamical models. |
|
Equation of motion solver for particle motion in free space. |
|
Equation of motion solver for a particle orbiting the Sun. |
|
Equation of motion solver for a particle orbiting the Sun, including radiation force. |
DynamicalModel solver failed. |
Class Inheritance Diagram¶
sbpy.dynamics.syndynes Module¶
Classes¶
|
Syndyne / synchrone generator for cometary dust. |
|
Abstract base class for particle states that make up a syndyne or synchrone. |
|
Immutable collection of syndynes or synchrones. |
|
Collection of particle states that make up a syndyne. |
|
Collection of syndynes. |
|
Collection of particle states that make up a synchrone. |
|
Collection of synchrones. |