mirror of
https://github.com/cosinekitty/astronomy.git
synced 2026-03-29 03:43:34 -04:00
PY gravim: work in progress.
Started implementation of the Python gravity simulator. Updated the `pydown` markdown generator to include class constructors `__init__` when they contain docstrings.
This commit is contained in:
Don Cross
@@ -436,6 +436,80 @@ not be used.
|
||||
|
||||
---
|
||||
|
||||
<a name="GravitySimulator"></a>
|
||||
### class GravitySimulator
|
||||
|
||||
**A simulation of zero or more small bodies moving through the Solar System.**
|
||||
|
||||
This class calculates the movement of arbitrary small bodies,
|
||||
such as asteroids or comets, that move through the Solar System.
|
||||
It does so by calculating the gravitational forces on the bodies
|
||||
from the Sun and planets. The user of this class supplies a
|
||||
list of initial positions and velocities for the bodies.
|
||||
Then the class can update the positions and velocities over small
|
||||
time steps.
|
||||
|
||||
#### member functions
|
||||
|
||||
<a name="GravitySimulator.__init__"></a>
|
||||
### GravitySimulator.__init__(self, originBody, time, bodyStates)
|
||||
|
||||
**Creates a gravity simulation object.**
|
||||
|
||||
| Type | Parameter | Description |
|
||||
| --- | --- | --- |
|
||||
| [`Body`](#Body) | `originBody` | Specifies the origin of the reference frame. All position vectors and velocity vectors will use `originBody` as the origin of the coordinate system. This origin applies to all the input vectors provided in the `bodyStates` parameter of this function, along with all output vectors returned by [`GravitySimulator`](#GravitySimulator).Update. Most callers will want to provide one of the following: `Body.Sun` for heliocentric coordinates, `Body.SSB` for solar system barycentric coordinates, or `Body.Earth` for geocentric coordinates. Note that the gravity simulator does not correct for light travel time; all state vectors are tied to a Newtonian "instantaneous" time. |
|
||||
| [`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.SolarSystemBodyState"></a>
|
||||
### GravitySimulator.SolarSystemBodyState(self, body)
|
||||
|
||||
**Get the position and velocity of a Solar System body included in the simulation.**
|
||||
|
||||
In order to simulate the movement of small bodies through the Solar System,
|
||||
the simulator needs to calculate the state vectors for the Sun and planets.
|
||||
If an application wants to know the positions of one or more of the planets
|
||||
in addition to the small bodies, this function provides a way to obtain
|
||||
their state vectors. This is provided for the sake of efficiency, to avoid
|
||||
redundant calculations.
|
||||
The state vector is returned relative to the position and velocity
|
||||
of the `originBody` parameter that was passed to this object's constructor.
|
||||
|
||||
| Type | Parameter | Description |
|
||||
| --- | --- | --- |
|
||||
| [`Body`](#Body) | `body` | The Sun, Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, or Neptune. |
|
||||
|
||||
<a name="GravitySimulator.Swap"></a>
|
||||
### GravitySimulator.Swap(self)
|
||||
|
||||
**Exchange the current time step with the previous time step.**
|
||||
|
||||
Sometimes it is helpful to "explore" various times near a given
|
||||
simulation time step, while repeatedly returning to the original
|
||||
time step. For example, when backdating a position for light travel
|
||||
time, the caller may wish to repeatedly try different amounts of
|
||||
backdating. When the backdating solver has converged, the caller
|
||||
wants to leave the simulation in its original state.
|
||||
This function allows a single "undo" of a simulation, and does so
|
||||
very efficiently.
|
||||
Usually this function will be called immediately after a matching
|
||||
call to [`GravitySimulator`](#GravitySimulator).Update. It has the effect of rolling
|
||||
back the most recent update. If called twice in a row, it reverts
|
||||
the swap and thus has no net effect.
|
||||
The constructor initializes the current state and previous
|
||||
state to be identical. Both states represent the `time` parameter that was
|
||||
passed into the constructor. Therefore, `Swap` will
|
||||
have no effect from the caller's point of view when passed a simulator
|
||||
that has not yet been updated by a call to [`GravitySimulator`](#GravitySimulator).Update.
|
||||
|
||||
<a name="GravitySimulator.Time"></a>
|
||||
### GravitySimulator.Time(self)
|
||||
|
||||
The time represented by the current step of the gravity simulation.
|
||||
|
||||
---
|
||||
|
||||
<a name="HorizontalCoordinates"></a>
|
||||
### class HorizontalCoordinates
|
||||
|
||||
|
||||
@@ -3602,6 +3602,11 @@ class _TerseVector:
|
||||
self.y = y
|
||||
self.z = z
|
||||
|
||||
@staticmethod
|
||||
def zero():
|
||||
'''Return a zero vector.'''
|
||||
return _TerseVector(0.0, 0.0, 0.0)
|
||||
|
||||
def ToAstroVector(self, time):
|
||||
'''Convert _TerseVector object to Vector object.'''
|
||||
return Vector(self.x, self.y, self.z, time)
|
||||
@@ -9834,3 +9839,203 @@ def LagrangePointFast(point, major_state, major_mass, minor_state, minor_mass):
|
||||
break
|
||||
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.
|
||||
|
||||
This class calculates the movement of arbitrary small bodies,
|
||||
such as asteroids or comets, that move through the Solar System.
|
||||
It does so by calculating the gravitational forces on the bodies
|
||||
from the Sun and planets. The user of this class supplies a
|
||||
list of initial positions and velocities for the bodies.
|
||||
Then the class can update the positions and velocities over small
|
||||
time steps.
|
||||
"""
|
||||
|
||||
def __init__(self, originBody, time, bodyStates):
|
||||
"""Creates a gravity simulation object.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
originBody: Body
|
||||
Specifies the origin of the reference frame.
|
||||
All position vectors and velocity vectors will use `originBody`
|
||||
as the origin of the coordinate system.
|
||||
This origin applies to all the input vectors provided in the
|
||||
`bodyStates` parameter of this function, along with all
|
||||
output vectors returned by #GravitySimulator.Update.
|
||||
Most callers will want to provide one of the following:
|
||||
`Body.Sun` for heliocentric coordinates,
|
||||
`Body.SSB` for solar system barycentric coordinates,
|
||||
or `Body.Earth` for geocentric coordinates. Note that the
|
||||
gravity simulator does not correct for light travel time;
|
||||
all state vectors are tied to a Newtonian "instantaneous" time.
|
||||
time: Time
|
||||
The initial time at which to start the simulation.
|
||||
bodyStates: StateVector[]
|
||||
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.
|
||||
"""
|
||||
self.originBody = originBody
|
||||
# Verify that the state vectors have matching times.
|
||||
for b in bodyStates:
|
||||
if b.t.tt != time.tt:
|
||||
raise Error('Inconsistent times in bodyStates')
|
||||
|
||||
# Create a stub list of small body states that we will append to later.
|
||||
# We just need the stub to put into `self.curr`
|
||||
smallBodyList = []
|
||||
|
||||
# Calculate the states of the Sun and planets at the initial time.
|
||||
largeBodyDict = _CalcSolarSystem(time)
|
||||
|
||||
# Create a simulation endpoint for the initial time.
|
||||
self.curr = _GravSimEndpoint(time, largeBodyDict, smallBodyList)
|
||||
|
||||
# Convert origin-centric bodyStates vectors into a barycentric `_body_grav_calc_t` array.
|
||||
o = self._InternalBodyState(originBody)
|
||||
for b in bodyStates:
|
||||
r = _TerseVector(b.x + o.r.x, b.y + o.r.y, b.z + o.r.z)
|
||||
v = _TerseVector(b.vx + o.v.x, b.vy + o.v.y, b.vz + o.v.z)
|
||||
a = _TerseVector.zero()
|
||||
smallBodyList.append(_body_grav_calc_t(time.tt, r, v, a))
|
||||
|
||||
# Calculate the net acceleration experienced by the small bodies.
|
||||
self._CalcBodyAccelerations()
|
||||
|
||||
# To prepare for a possible swap operation, duplicate the current state into the previous state.
|
||||
self.prev = self._Duplicate()
|
||||
|
||||
def Time(self):
|
||||
"""The time represented by the current step of the gravity simulation."""
|
||||
return self.curr.time
|
||||
|
||||
def Swap(self):
|
||||
"""Exchange the current time step with the previous time step.
|
||||
|
||||
Sometimes it is helpful to "explore" various times near a given
|
||||
simulation time step, while repeatedly returning to the original
|
||||
time step. For example, when backdating a position for light travel
|
||||
time, the caller may wish to repeatedly try different amounts of
|
||||
backdating. When the backdating solver has converged, the caller
|
||||
wants to leave the simulation in its original state.
|
||||
|
||||
This function allows a single "undo" of a simulation, and does so
|
||||
very efficiently.
|
||||
|
||||
Usually this function will be called immediately after a matching
|
||||
call to #GravitySimulator.Update. It has the effect of rolling
|
||||
back the most recent update. If called twice in a row, it reverts
|
||||
the swap and thus has no net effect.
|
||||
|
||||
The constructor initializes the current state and previous
|
||||
state to be identical. Both states represent the `time` parameter that was
|
||||
passed into the constructor. Therefore, `Swap` will
|
||||
have no effect from the caller's point of view when passed a simulator
|
||||
that has not yet been updated by a call to #GravitySimulator.Update.
|
||||
"""
|
||||
(self.curr, self.prev) = (self.prev, self.curr)
|
||||
|
||||
def SolarSystemBodyState(self, body):
|
||||
"""Get the position and velocity of a Solar System body included in the simulation.
|
||||
|
||||
In order to simulate the movement of small bodies through the Solar System,
|
||||
the simulator needs to calculate the state vectors for the Sun and planets.
|
||||
|
||||
If an application wants to know the positions of one or more of the planets
|
||||
in addition to the small bodies, this function provides a way to obtain
|
||||
their state vectors. This is provided for the sake of efficiency, to avoid
|
||||
redundant calculations.
|
||||
|
||||
The state vector is returned relative to the position and velocity
|
||||
of the `originBody` parameter that was passed to this object's constructor.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
body : Body
|
||||
The Sun, Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, or Neptune.
|
||||
"""
|
||||
bstate = self._InternalBodyState(body)
|
||||
ostate = self._InternalBodyState(self.originBody)
|
||||
return _ExportState(bstate - ostate, self.curr.time)
|
||||
|
||||
def _InternalBodyState(self, body):
|
||||
if body == Body.SSB:
|
||||
return _body_state_t(self.curr.time.tt, _TerseVector.zero(), _TerseVector.zero())
|
||||
if body in self.curr.gravitators:
|
||||
return self.curr.gravitators[body]
|
||||
raise Error('Invalid body: {}'.format(body))
|
||||
|
||||
def _CalcBodyAccelerations(self):
|
||||
for b in self.curr.bodies.values():
|
||||
b.a = _TerseVector.zero()
|
||||
_AddAcceleration(b.a, b.r, self.curr.gravitators[Body.Sun ].r, _SUN_GM)
|
||||
_AddAcceleration(b.a, b.r, self.curr.gravitators[Body.Mercury].r, _MERCURY_GM)
|
||||
_AddAcceleration(b.a, b.r, self.curr.gravitators[Body.Venus ].r, _VENUS_GM)
|
||||
_AddAcceleration(b.a, b.r, self.curr.gravitators[Body.Earth ].r, _EARTH_GM + _MOON_GM)
|
||||
_AddAcceleration(b.a, b.r, self.curr.gravitators[Body.Mars ].r, _MARS_GM)
|
||||
_AddAcceleration(b.a, b.r, self.curr.gravitators[Body.Jupiter].r, _JUPITER_GM)
|
||||
_AddAcceleration(b.a, b.r, self.curr.gravitators[Body.Saturn ].r, _SATURN_GM)
|
||||
_AddAcceleration(b.a, b.r, self.curr.gravitators[Body.Uranus ].r, _URANUS_GM)
|
||||
_AddAcceleration(b.a, b.r, self.curr.gravitators[Body.Neptune].r, _NEPTUNE_GM)
|
||||
|
||||
def _Duplicate(self):
|
||||
# Copy the current stateinto the previous state, so that both become the same moment in time.
|
||||
gravitators = {}
|
||||
for body in self.curr.gravitators.values():
|
||||
gravitators[body] = self.curr.gravitators[body].clone()
|
||||
|
||||
bodies = []
|
||||
for b in self.curr.bodies:
|
||||
bodies.append(b.clone())
|
||||
|
||||
return _GravSimEndpoint(self.curr.time, gravitators, bodies)
|
||||
|
||||
class _GravSimEndpoint:
|
||||
def __init__(self, time, gravitators, bodies):
|
||||
self.time = time
|
||||
self.gravitators = gravitators
|
||||
self.bodies = bodies
|
||||
|
||||
def _CalcSolarSystem(time):
|
||||
d = {}
|
||||
# Start with the SSB at zero position and velocity.
|
||||
ssb = _body_state_t(time.tt, _TerseVector.zero(), _TerseVector.zero())
|
||||
|
||||
# Calculate the heliocentric position of each planet, and adjust the SSB
|
||||
# based on each planet's pull on the Sun.
|
||||
d[Body.Mercury] = _AdjustBarycenterPosVel(ssb, time.tt, Body.Mercury, _MERCURY_GM)
|
||||
d[Body.Venus ] = _AdjustBarycenterPosVel(ssb, time.tt, Body.Venus, _VENUS_GM)
|
||||
d[Body.Earth ] = _AdjustBarycenterPosVel(ssb, time.tt, Body.Earth, _EARTH_GM + _MOON_GM)
|
||||
d[Body.Mars ] = _AdjustBarycenterPosVel(ssb, time.tt, Body.Mars, _MARS_GM)
|
||||
d[Body.Jupiter] = _AdjustBarycenterPosVel(ssb, time.tt, Body.Jupiter, _JUPITER_GM)
|
||||
d[Body.Saturn ] = _AdjustBarycenterPosVel(ssb, time.tt, Body.Saturn, _SATURN_GM)
|
||||
d[Body.Uranus ] = _AdjustBarycenterPosVel(ssb, time.tt, Body.Uranus, _URANUS_GM)
|
||||
d[Body.Neptune] = _AdjustBarycenterPosVel(ssb, time.tt, Body.Neptune, _NEPTUNE_GM)
|
||||
|
||||
# Convert planet states from heliocentric to barycentric.
|
||||
for b in d.values():
|
||||
b.r -= ssb.r
|
||||
b.v -= ssb.v
|
||||
|
||||
# Convert heliocentric SSB to barycentric Sun.
|
||||
d[Body.Sun] = _body_state_t(time.tt, -1.0 * ssb.r, -1.0 * ssb.v)
|
||||
return d
|
||||
|
||||
def _AddAcceleration(acc, smallPos, majorPos, gm):
|
||||
dx = majorPos.x - smallPos.x
|
||||
dy = majorPos.y - smallPos.y
|
||||
dz = majorPos.z - smallPos.z
|
||||
r2 = dx*dx + dy*dy + dz*dz
|
||||
pull = gm / (r2 * math.sqrt(r2))
|
||||
acc.x += dx * pull
|
||||
acc.y += dy * pull
|
||||
acc.z += dz * pull
|
||||
|
||||
Reference in New Issue
Block a user