search

NASA-Planetary-Science / sbpy

A Python package for small bodies research
https://sbpy.org/
Other
67 stars 34 forks source link

Add support for dynamical models and dust syndynes/synchrones. #394

Closed mkelley closed 1 month ago

mkelley commented 8 months ago

Adds the capability to generate syndynes and synchrones. This PR is essentially four enhancements:

  1. Add State objects to hold particle position, velocity, and time.
  2. Add SPICE-equivalent "ephemeris time" capability to the astropy.time system.
  3. Add the DynamicalModel abstract base class to define an API for solving the equations of motion, and implement a few simple models.
  4. Add SynGenerator and supporting classes which can generate syndynes/synchrones.

Design choices are given below.

State objects

Dynamical state is 3D position, velocity, and time. The sbpy.dynamics.State object encapsulates these quantities using the astropy.coordinates sub-module. State vectors are internally stored as a coordinate frame object with position and velocity data. Example coordinate frame objects are HeliocentricEclipticIAU76, FK5, and ICRS. Some coordinate frames require observation time for transformation to other coordinate frames, e.g., from a heliocentric to a barycentric frame. Time is added to the frame as needed. In addition, sbpy adds an ArbitraryFrame, which is used by default and cannot be converted to other frames.

The state's time attribute is typically an astropy.time.Time object. However, when the ArbitraryFrame is used the state's time attribute may be an astropy.units.Quantity.

Create a State object for a particle at x=2 au, moving with v_y=30 km/s, on 2023 Dec 08:

>>> 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)
>>> comet
<State (<ArbitraryFrame Frame>):
 r
  [2. 0. 0.] AU
 v
  [ 0. 30.  0.] km / s
 t
  2023-12-08 00:00:00.000>

States can also be arrays of states:

>>> len(comet)
1
>>> comets = State([r, r], [v, v], [t, t])
>>> len(comets)
2
>>> comet.r.shape  # [x, y, z]
(3,)
>>> comets.r.shape  # [[x, y, z], [x, y, z]]
(2, 3)

To specify a coordinate frame, use the frame keyword argument:

>>> State(r, v, t, frame="heliocentriceclipticiau76")
<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>

Alternative implementations

Initially, State stored its data as r, v, and t arrays. This involved extra code for parameter checking and coordinate transformations. Switching to using BaseCoordinateFrame objects simplified both of these aspects. However, t still must be saved internally, as some coordinate frames do not use time.

State could instead use astropy.coordinates.SkyCoord. SkyCoord can store 3D vectors, but is also designed to accommodate 2D representations and operations on the celestial sphere. Dynamical state must only be 3D. BaseCoordinateFrame data could also be 2D, but is more general and lacks many of SkyCoord's 2D-based convenience methods. I also found confusing SkyCoord's tendency to duplicate position and velocity vectors in its coordinate frame objects:

>>> coords = SkyCoord(ra=1 * u.deg, dec=2 * u.deg, distance=1 * u.au)
>>> coords
<SkyCoord (ICRS): (ra, dec, distance) in (deg, deg, AU)
    (1., 2., 1.)>

>>> coords.frame
<ICRS Coordinate: (ra, dec, distance) in (deg, deg, AU)
    (1., 2., 1.)>

Sub-classing SkyCoord or BaseCoordinateFrame would complicate the State API, especially if SkyCoord was sub-classed, whereas the goal is to keep the API simple.

State could have been based on Ephem. Ephem lacks coordinate frame information. Functionally, this would be similar to saving r, v, and t as their own arrays.

Ephemeris time

NAIF SPICE uses TDB seconds from the J2000 epoch for time coordinates. This is not explicitly supported in astropy.time, so a new SpiceEphemerisTime has been created. By importing the sbpy.time sub-module, it becomes available for use as the format "et":

>>> from astropy.time import Time
>>> import sbpy.time
>>>
>>> j2000 = Time(0, format="et")
>>> j2000.iso
'2000-01-01 12:00:00.000'
>>> Time(759307746.954761, format="et")
<Time object: scale='tdb' format='et' value=759307746.954761>

Dynamical integrations (below) will convert Time objects to ephemeris time.

Dynamical models

A lightweight approach to integrating the equations of motion has been implemented in sbpy.dynamics. The new DynamicalModel abstract base class defines the API and implements solving the equations with scipy.integrate.solve_ivp. The implementation is straightforward if we bump sbpy's minimum supported scipy version from 1.3 to 1.4, which allows for arguments to be passed from solve_ivp to the integrated functions.

solve_ivp keyword arguments, such as solution tolerances, or the choice of integrator, are set during DyanmicalModel initialization. The default arguments were designed for high precision: a relative tolerance of 2.3e-14, use of a function for the Jacobian matrix (as opposed to numerical estimates), and use of the LSODA integrator. All of these may be overridden by the user at the time of initialization.

Inheriting classes provide methods that calculate the derivatives of position and velocity, and the Jacobian matrix. Time is given as a float in units of seconds, and position and velocity as a 6-element array ([x, y, z, v_x, v_y, v_z]) in units of kilometers and kilometers per second.

@abc.abstractclassmethod
def dx_dt(cls, t: float, rv: np.ndarray, *args) -> np.ndarray:
    """Derivative of position and velocity.
    ...
    """

@abc.abstractclassmethod
def df_drv(cls, t: float, rv: np.ndarray, *args) -> np.ndarray:
    """Jacobian matrix, :math:`df/drv`.
    ...
    """

The solution is generated using the solve method, whose arguments are the initial dynamical state and the time at which the solution is desired. Additional arguments for dx_dt and df_drv are possible, and necessary in order to pass the beta parameter for integrations considering radiation pressure.

Three dynamical model implementations are defined, FreeExpansion, SolarGravity, and SolarGravityAndRadiationPressure. The following example integrates a pure-gravitational orbit around the Sun:

>>> import astropy.units as u
>>> from sbpy.dynamics import State, 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>

Syndynes and synchrones

The capability to produce syndynes and synchrones has been added to sbpy.dynamics with the SynGenerator class. A State object is used for the dust source, and dust particle orbits are parameterized using beta, the ratio of the force from solar radiation to the force from solar gravity.

The source object state and dust particle states are integrated using the DynamicalModel framework. The SolarGravityAndRadiationPressure model is used by default, but others may be provided at the time the SynGenerator object is initialized. This framework should allow for more complex problems, such as planetary perturbations, or particle fragmentation, provided a dynamical model can be developed to support them.

>>> r = [2, 0, 0] * u.au
>>> v = [0, 30, 0] * u.km / u.s
>>> t = Time("2023-12-08")
>>> comet = State(r, v, t)
>>> betas = [1, 0.1, 0.01, 0]
>>> ages = np.linspace(0, 100, 26) * u.day
>>> dust = SynGenerator(comet, betas, ages)

The integrated particle states are saved in the particles attribute. However, it is more convenient get the syndynes and synchrones packaged as Syndyne and Synchrone objects, which are specialized State objects:

>>> syndyne = dust.syndyne(0)  # get the first syndyne
>>> len(syndyne)
26
>>> syndyne.x
<Quantity [2.99195741e+08, 2.99284224e+08, 2.99549029e+08, 2.99988249e+08,
           3.00598735e+08, 3.01376158e+08, 3.02315083e+08, 3.03409054e+08,
           3.04650704e+08, 3.06031870e+08, 3.07543713e+08, 3.09176847e+08,
           3.10921470e+08, 3.12767487e+08, 3.14704634e+08, 3.16722591e+08,
           3.18811090e+08, 3.20960013e+08, 3.23159471e+08, 3.25399886e+08,
           3.27672045e+08, 3.29967158e+08, 3.32276893e+08, 3.34593405e+08,
           3.36909361e+08, 3.39217944e+08] km>

If an observer is provided to the SynGenerator object at initialization, then the syndyne/synchrone object will have a coords attribute, which is a SkyCoord object containing the observed coordinates of the syndyne/synchrone test particles.

>>> observer = State([0, 1, 0] * u.km, [30, 0, 0] * u.km / u.s, comet.t)
>>> dust = SynGenerator(comet, betas, ages, observer=observer)
>>> syndyne = dust.syndyne(0)
>>> syndyne.coords # equivalent to: observer.observe(syndyne)
<SkyCoord (ArbitraryFrame): (lon, lat, distance) in (deg, deg, km)
    [(359.99999981, 0., 2.99195741e+08),
     (359.99960848, 0., 2.99284224e+08),
     (359.99687763, 0., 2.99549030e+08),
     (359.98950966, 0., 2.99988254e+08),
     (359.97528925, 0., 3.00598763e+08),
     (359.95212067, 0., 3.01376263e+08),
     (359.91806084, 0., 3.02315392e+08),
     (359.87134702, 0., 3.03409819e+08),
...

Syndynes and synchrones may be converted to Ephem objects. The sbpy DataClass fields were updated to include the beta parameter (beta_rad) and relative time (t_relative).

Alternative implementations

Initially, SynGenerator had a two step process:

>>> dust = SynGenerator(comet, betas, ages)
>>> dust.solve()

The initial states were generated at initialization, and final states by the solve() method. Two steps allows for the initial particle states to be modified by the user before being integrated, e.g., to give them a non-zero ejection velocity. This is still possible in the implemented design, but would be instead addressed by sub-classing and overriding initialize_states:

>>> class AlternativeSynGenerator(SynGenerator):
...     def initialize_states(self):
...         super().initialize_states()
...         # then modify self.initial_states as needed

Other changes

There are some improvements to the activity.dust documentation. Some package infrastructure files have been edited with small improvements. The pytest configuration now respects the options in setup.cfg.

Addresses #19.

pep8speaks commented 8 months ago

Hello @mkelley! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

Line 27:80: E501 line too long (88 > 79 characters) Line 38:80: E501 line too long (87 > 79 characters) Line 75:80: E501 line too long (88 > 79 characters) Line 126:80: E501 line too long (84 > 79 characters) Line 189:80: E501 line too long (82 > 79 characters) Line 484:80: E501 line too long (86 > 79 characters) Line 496:80: E501 line too long (85 > 79 characters)

Line 76:80: E501 line too long (80 > 79 characters) Line 234:80: E501 line too long (84 > 79 characters) Line 246:80: E501 line too long (87 > 79 characters) Line 264:80: E501 line too long (87 > 79 characters) Line 334:80: E501 line too long (81 > 79 characters) Line 340:80: E501 line too long (81 > 79 characters) Line 346:80: E501 line too long (83 > 79 characters) Line 894:80: E501 line too long (81 > 79 characters) Line 960:80: E501 line too long (89 > 79 characters)

Line 432:80: E501 line too long (80 > 79 characters) Line 438:80: E501 line too long (80 > 79 characters) Line 442:80: E501 line too long (80 > 79 characters)

Line 53:80: E501 line too long (80 > 79 characters) Line 92:80: E501 line too long (81 > 79 characters) Line 190:80: E501 line too long (80 > 79 characters) Line 225:80: E501 line too long (80 > 79 characters) Line 286:80: E501 line too long (92 > 79 characters) Line 311:80: E501 line too long (80 > 79 characters) Line 344:80: E501 line too long (80 > 79 characters)

Line 32:80: E501 line too long (84 > 79 characters) Line 55:80: E501 line too long (82 > 79 characters) Line 58:80: E501 line too long (88 > 79 characters) Line 79:80: E501 line too long (80 > 79 characters) Line 128:80: E501 line too long (83 > 79 characters) Line 137:80: E501 line too long (80 > 79 characters) Line 140:80: E501 line too long (84 > 79 characters) Line 142:80: E501 line too long (81 > 79 characters) Line 298:80: E501 line too long (88 > 79 characters) Line 415:80: E501 line too long (83 > 79 characters) Line 418:80: E501 line too long (86 > 79 characters)

Line 43:80: E501 line too long (86 > 79 characters) Line 57:80: E501 line too long (84 > 79 characters) Line 135:80: E501 line too long (87 > 79 characters) Line 196:80: E501 line too long (81 > 79 characters) Line 199:80: E501 line too long (81 > 79 characters) Line 220:80: E501 line too long (81 > 79 characters) Line 268:80: E501 line too long (81 > 79 characters) Line 297:24: E721 do not compare types, use 'isinstance()' Line 437:80: E501 line too long (83 > 79 characters) Line 455:80: E501 line too long (80 > 79 characters) Line 517:80: E501 line too long (82 > 79 characters) Line 606:80: E501 line too long (83 > 79 characters)

Line 57:80: E501 line too long (88 > 79 characters) Line 72:80: E501 line too long (88 > 79 characters) Line 163:80: E501 line too long (86 > 79 characters) Line 193:80: E501 line too long (86 > 79 characters)

Line 41:80: E501 line too long (83 > 79 characters) Line 54:80: E501 line too long (82 > 79 characters) Line 114:80: E501 line too long (95 > 79 characters) Line 158:80: E501 line too long (86 > 79 characters) Line 197:80: E501 line too long (83 > 79 characters) Line 220:80: E501 line too long (84 > 79 characters) Line 239:80: E501 line too long (80 > 79 characters) Line 279:80: E501 line too long (82 > 79 characters) Line 285:80: E501 line too long (94 > 79 characters) Line 286:80: E501 line too long (88 > 79 characters) Line 287:80: E501 line too long (110 > 79 characters) Line 290:80: E501 line too long (88 > 79 characters) Line 312:80: E501 line too long (81 > 79 characters) Line 340:80: E501 line too long (88 > 79 characters) Line 372:80: E501 line too long (93 > 79 characters) Line 382:80: E501 line too long (118 > 79 characters) Line 391:80: E501 line too long (95 > 79 characters) Line 396:80: E501 line too long (97 > 79 characters) Line 435:80: E501 line too long (80 > 79 characters) Line 438:80: E501 line too long (80 > 79 characters) Line 446:80: E501 line too long (81 > 79 characters) Line 453:80: E501 line too long (80 > 79 characters) Line 490:80: E501 line too long (82 > 79 characters) Line 491:80: E501 line too long (91 > 79 characters) Line 509:80: E501 line too long (80 > 79 characters) Line 559:80: E501 line too long (86 > 79 characters)

Line 166:80: E501 line too long (83 > 79 characters) Line 190:80: E501 line too long (83 > 79 characters) Line 203:80: E501 line too long (83 > 79 characters) Line 337:80: E501 line too long (87 > 79 characters) Line 373:80: E501 line too long (85 > 79 characters) Line 378:80: E501 line too long (84 > 79 characters) Line 392:80: E501 line too long (83 > 79 characters) Line 396:80: E501 line too long (89 > 79 characters) Line 406:80: E501 line too long (85 > 79 characters) Line 438:80: E501 line too long (85 > 79 characters) Line 440:80: E501 line too long (84 > 79 characters) Line 451:80: E501 line too long (86 > 79 characters) Line 463:80: E501 line too long (81 > 79 characters) Line 478:80: E501 line too long (83 > 79 characters) Line 481:80: E501 line too long (86 > 79 characters) Line 494:80: E501 line too long (80 > 79 characters) Line 495:80: E501 line too long (86 > 79 characters)

Line 69:80: E501 line too long (84 > 79 characters) Line 124:80: E501 line too long (84 > 79 characters)

Comment last updated at 2024-08-20 19:32:27 UTC
codecov[bot] commented 8 months ago

Codecov Report

Attention: Patch coverage is 99.81343% with 2 lines in your changes missing coverage. Please review.

Project coverage is 82.99%. Comparing base (eaff110) to head (057efbc). Report is 1 commits behind head on main.

Files Patch % Lines
sbpy/dynamics/state.py 98.75% 2 Missing :warning:
Additional details and impacted files ```diff @@ Coverage Diff @@ ## main #394 +/- ## ========================================== + Coverage 80.49% 82.99% +2.49% ========================================== Files 79 88 +9 Lines 7070 8108 +1038 ========================================== + Hits 5691 6729 +1038 Misses 1379 1379 ```

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

mkelley commented 8 months ago

@hhsieh00 FYI, readthedocs builds documentation on pull requests (see the "checks" below). For your convenience, the time and dynamics documentation are at:

https://sbpy--394.org.readthedocs.build/en/394/sbpy/time.html https://sbpy--394.org.readthedocs.build/en/394/sbpy/dynamics.html

mkelley commented 1 month ago

Unfortunately, removing the blackslashes causes Sphinx to interpret those lines as parts of an enumerated list, so "(2001) Einstein" becomes "2001. Einstein". However, we can fix that by following the closing parenthesis with a non-breaking space.

https://docutils.sourceforge.io/docs/ref/rst/restructuredtext.html#enumerated-lists

mkelley commented 1 month ago

👍🏻 An example has been added.