mirror/astronomy
1
0
Fork 0
mirror of https://github.com/cosinekitty/astronomy.git synced 2026-05-19 06:17:03 -04:00
Code Issues Packages Projects Releases Wiki Activity

PY gravsim: initial coding completed

Browse Source 
Finished coding the Python version of the gravity simulator.
No unit tests have been written yet.
Cleaned up documentation in the other languages.
Made some functions static that did not need to be members.
This commit is contained in:
Don Cross
2022-05-22 09:41:44 -04:00
parent 68f118587d
commit 4303137c0a
27 changed files with 681 additions and 269 deletions
Download Patch File Download Diff File Expand all files Collapse all files

35
source/python/README.md
View File

@@ -462,6 +462,13 @@ time steps.
| [`Time`](#Time) | `time` | The initial time at which to start the simulation. |
| [`StateVector[]`](#StateVector[]) | `bodyStates` | An array of zero or more initial state vectors (positions and velocities) of the small bodies to be simulated. The caller must know the positions and velocities of the small bodies at an initial moment in time. Their positions and velocities are expressed with respect to `originBody`, using equatorial J2000 orientation (EQJ). Positions are expressed in astronomical units (AU). Velocities are expressed in AU/day. All the times embedded within the state vectors must exactly match `time`, or this constructor will throw an exception. |
<a name="GravitySimulator.OriginBody"></a>
### GravitySimulator.OriginBody(self)
**The origin of the reference frame. See constructor for more info.**
### Returns: [`Body`](#Body)
<a name="GravitySimulator.SolarSystemBodyState"></a>
### GravitySimulator.SolarSystemBodyState(self, body)
@@ -480,6 +487,9 @@ of the `originBody` parameter that was passed to this object's constructor.
| --- | --- | --- |
| [`Body`](#Body) | `body` | The Sun, Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, or Neptune. |
### Returns: [`StateVector`](#StateVector)
The state vector of the requested Solar System body.
<a name="GravitySimulator.Swap"></a>
### GravitySimulator.Swap(self)
@@ -506,7 +516,30 @@ that has not yet been updated by a call to [`GravitySimulator`](#GravitySimulato
<a name="GravitySimulator.Time"></a>
### GravitySimulator.Time(self)
The time represented by the current step of the gravity simulation.
**The time represented by the current step of the gravity simulation.**
### Returns: [`Time`](#Time)
<a name="GravitySimulator.Update"></a>
### GravitySimulator.Update(self, time)
**Advances the gravity simulation by a small time step.**
Updates the simulation of the user-supplied small bodies
to the time indicated by the `time` parameter.
Returns an array of state vectors for the simulated bodies.
The array is in the same order as the original array that
was used to construct this simulator object.
The positions and velocities in the returned array are
referenced to the `originBody` that was used to construct
this simulator.
| Type | Parameter | Description |
| --- | --- | --- |
| [`Time`](#Time) | `time` | A time that is a small increment away from the current simulation time. It is up to the developer to figure out an appropriate time increment. Depending on the trajectories, a smaller or larger increment may be needed for the desired accuracy. Some experimentation may be needed. Generally, bodies that stay in the outer Solar System and move slowly can use larger time steps. Bodies that pass into the inner Solar System and move faster will need a smaller time step to maintain accuracy. The `time` value may be after or before the current simulation time to move forward or backward in time. |
### Returns: [`StateVector[]`](#StateVector[])
An array of state vectors, one for each small body.
---

128
source/python/astronomy/astronomy.py
View File

@@ -9840,6 +9840,7 @@ def LagrangePointFast(point, major_state, major_mass, minor_state, minor_mass):
scale = (x - r1) / R
return StateVector(scale*dx, scale*dy, scale*dz, scale*vx, scale*vy, scale*vz, major_state.t)
#--------------------------------------------------------------------------------------------------
class GravitySimulator:
"""A simulation of zero or more small bodies moving through the Solar System.
@@ -9884,7 +9885,7 @@ class GravitySimulator:
All the times embedded within the state vectors must exactly match `time`,
or this constructor will throw an exception.
"""
self.originBody = originBody
self._originBody = originBody
# Verify that the state vectors have matching times.
for b in bodyStates:
if b.t.tt != time.tt:
@@ -9915,9 +9916,123 @@ class GravitySimulator:
self.prev = self._Duplicate()
def Time(self):
"""The time represented by the current step of the gravity simulation."""
"""The time represented by the current step of the gravity simulation.
Returns
-------
Time
"""
return self.curr.time
def OriginBody(self):
"""The origin of the reference frame. See constructor for more info.
Returns
-------
Body
"""
return self._originBody
def Update(self, time):
"""Advances the gravity simulation by a small time step.
Updates the simulation of the user-supplied small bodies
to the time indicated by the `time` parameter.
Returns an array of state vectors for the simulated bodies.
The array is in the same order as the original array that
was used to construct this simulator object.
The positions and velocities in the returned array are
referenced to the `originBody` that was used to construct
this simulator.
Parameters
----------
time : Time
A time that is a small increment away from the current simulation time.
It is up to the developer to figure out an appropriate time increment.
Depending on the trajectories, a smaller or larger increment
may be needed for the desired accuracy. Some experimentation may be needed.
Generally, bodies that stay in the outer Solar System and move slowly can
use larger time steps. Bodies that pass into the inner Solar System and
move faster will need a smaller time step to maintain accuracy.
The `time` value may be after or before the current simulation time
to move forward or backward in time.
Returns
-------
StateVector[]
An array of state vectors, one for each small body.
"""
dt = time.tt - self.curr.time.tt
if dt == 0.0:
# Special case: the time has not changed, so skip the usual physics calculations.
# This allows another way for the caller to query the current body states.
# It is also necessary to avoid dividing by `dt` if `dt` is zero.
# To prepare for a possible swap operation, duplicate the current state into the previous state.
self.prev = self._Duplicate()
else:
# Exchange the current state with the previous state. Then calculate the new current state.
self.Swap()
# Update the current time
self.curr.time = time
# Calculate the positions and velocities of the Sun and planets at the given time.
self.curr.gravitators = _CalcSolarSystem(time)
# Estimate the positions of the small bodies as if their existing
# existing accelerations apply across the whole time interval.
nbodies = len(self.curr.bodies)
for i in range(nbodies):
p = self.prev.bodies[i]
c = self.curr.bodies[i]
c.r = _UpdatePosition(dt, p.r, p.v, p.a)
# Calculate the acceleration experienced by the small bodies at
# their respective approximate next locations.
self._CalcBodyAccelerations()
for i in range(nbodies):
# Calculate the average of the acceleration vectors
# experienced by the previous body positions and
# their estimated next positions.
# These become estimates of the mean effective accelerations
# over the whole interval.
p = self.prev.bodies[i]
c = self.curr.bodies[i]
acc = p.a.mean(c.a)
# Refine the estimates of position and velocity at the next time step,
# using the mean acceleration as a better approximation of the
# continuously changing acceleration acting on each body.
c.tt = time.tt
c.r = _UpdatePosition(dt, p.r, p.v, acc)
c.v = _UpdateVelocity(dt, p.v, acc)
# Re-calculate accelerations experienced by each body.
# These will be needed for the next simulation step (if any).
# Also, they will be potentially useful if some day we add
# a function to query the acceleration vectors for the bodies.
self._CalcBodyAccelerations()
# Translate our internal calculations of body positions and velocities
# into state vectors that the caller can understand.
# We have to convert the internal type _body_grav_calc_t to the public type StateVector.
# Also convert from barycentric coordinates to coordinates based on the selected origin body.
bodyStates = []
ostate = self._InternalBodyState(self._originBody)
for bcalc in self.curr.bodies:
bodyStates.append(StateVector(
bcalc.r.x - ostate.r.x,
bcalc.r.y - ostate.r.y,
bcalc.r.z - ostate.r.z,
bcalc.v.x - ostate.v.x,
bcalc.v.y - ostate.v.y,
bcalc.v.z - ostate.v.z,
time
))
return bodyStates
def Swap(self):
"""Exchange the current time step with the previous time step.
@@ -9962,9 +10077,14 @@ class GravitySimulator:
----------
body : Body
The Sun, Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, or Neptune.
Returns
-------
StateVector
The state vector of the requested Solar System body.
"""
bstate = self._InternalBodyState(body)
ostate = self._InternalBodyState(self.originBody)
ostate = self._InternalBodyState(self._originBody)
return _ExportState(bstate - ostate, self.curr.time)
def _InternalBodyState(self, body):
@@ -10039,3 +10159,5 @@ def _AddAcceleration(acc, smallPos, majorPos, gm):
acc.x += dx * pull
acc.y += dy * pull
acc.z += dz * pull
#--------------------------------------------------------------------------------------------------