mirror of
https://github.com/cosinekitty/astronomy.git
synced 2026-05-19 14:27:52 -04:00
Kotlin gravsim: more work in progress.
Added the function GravitySimulator.update(). Not yet tested.
This commit is contained in:
Don Cross
@@ -2731,12 +2731,12 @@ static void CalcSolarSystem(astro_grav_sim_t *sim)
|
||||
static void CalcBodyAccelerations(astro_grav_sim_t *sim)
|
||||
{
|
||||
int i;
|
||||
const body_state_t *grav = sim->curr->gravitators;
|
||||
|
||||
/* Calculate the gravitational acceleration experienced by the simulated bodies. */
|
||||
for (i = 0; i < sim->numBodies; ++i)
|
||||
{
|
||||
body_grav_calc_t *calc = &sim->curr->bodies[i];
|
||||
const body_state_t *grav = sim->curr->gravitators;
|
||||
|
||||
calc->a = VecZero;
|
||||
|
||||
@@ -2978,7 +2978,7 @@ fail:
|
||||
|
||||
|
||||
/**
|
||||
* @brief Advances a gravity simulation of a small body by a small time step.
|
||||
* @brief Advances a gravity simulation by a small time step.
|
||||
*
|
||||
* @param sim
|
||||
* A simulation object that was created by a prior call to #Astronomy_GravSimInit.
|
||||
|
||||
@@ -705,6 +705,12 @@ internal data class TerseVector(var x: Double, var y: Double, var z: Double) {
|
||||
z = -z
|
||||
}
|
||||
|
||||
fun copyFrom(other: TerseVector) {
|
||||
x = other.x
|
||||
y = other.y
|
||||
z = other.z
|
||||
}
|
||||
|
||||
companion object {
|
||||
@JvmStatic
|
||||
fun zero() = TerseVector(0.0, 0.0, 0.0)
|
||||
@@ -2944,10 +2950,21 @@ internal class BodyState(
|
||||
val r: TerseVector,
|
||||
val v: TerseVector
|
||||
) {
|
||||
fun increment(other: BodyState) {
|
||||
r.increment(other.r)
|
||||
v.increment(other.v)
|
||||
}
|
||||
|
||||
fun decrement(other: BodyState) {
|
||||
r.decrement(other.r)
|
||||
v.decrement(other.v)
|
||||
}
|
||||
|
||||
fun copyFrom(other: BodyState) {
|
||||
tt = other.tt
|
||||
r.copyFrom(other.r)
|
||||
v.copyFrom(other.v)
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@@ -2959,6 +2976,14 @@ private fun exportState(bodyState: BodyState, time: Time) =
|
||||
)
|
||||
|
||||
|
||||
private fun exportGravCalc(calc: BodyGravCalc, time: Time) =
|
||||
StateVector(
|
||||
calc.r.x, calc.r.y, calc.r.z,
|
||||
calc.v.x, calc.v.y, calc.v.z,
|
||||
time
|
||||
)
|
||||
|
||||
|
||||
private fun calcVsopPosVel(model: VsopModel, tt: Double): BodyState {
|
||||
val t = tt / DAYS_PER_MILLENNIUM
|
||||
|
||||
@@ -3231,11 +3256,18 @@ private class MoonContext(time: Time) {
|
||||
// Pluto/SSB gravitation simulator
|
||||
|
||||
private class BodyGravCalc(
|
||||
val tt: Double, // J2000 terrestrial time [days]
|
||||
val r: TerseVector, // position [au]
|
||||
val v: TerseVector, // velocity [au/day]
|
||||
val a: TerseVector // acceleration [au/day]
|
||||
)
|
||||
var tt: Double, // J2000 terrestrial time [days]
|
||||
var r: TerseVector, // position [au]
|
||||
var v: TerseVector, // velocity [au/day]
|
||||
var a: TerseVector // acceleration [au/day]
|
||||
) {
|
||||
fun copyFrom(other: BodyGravCalc) {
|
||||
tt = other.tt
|
||||
r.copyFrom(other.r)
|
||||
v.copyFrom(other.v)
|
||||
a.copyFrom(other.a)
|
||||
}
|
||||
}
|
||||
|
||||
private class MajorBodies(
|
||||
val sun: BodyState,
|
||||
@@ -7872,17 +7904,189 @@ class GravitySimulator {
|
||||
|
||||
// Verify that all the state vectors are valid and have matching times.
|
||||
for (b in bodyStates) {
|
||||
if (b.t.tt != time.tt)
|
||||
if (b.t.tt != time.tt) {
|
||||
throw IllegalArgumentException("Inconsistent time(s) in bodyStates array.")
|
||||
}
|
||||
}
|
||||
|
||||
curr = initialEndpoint(time, bodyStates)
|
||||
prev = initialEndpoint(time, bodyStates)
|
||||
|
||||
// Calculate the states of the Sun and planets.
|
||||
calcSolarSystem()
|
||||
|
||||
// We need to do all the physics calculations in barycentric coordinates.
|
||||
// But the caller provides the input vectors with respect to `originBody`.
|
||||
// Correct the input body state vectors for the specified origin.
|
||||
if (originBody != Body.SSB) {
|
||||
// Determine the barycentric state of the origin body.
|
||||
val ostate = bodyState(originBody)
|
||||
|
||||
// Add barycentric origin to origin-centric bodies to obtain barycentric bodies.
|
||||
for (bstate in curr.bodies) {
|
||||
bstate.r.x += ostate.r.x
|
||||
bstate.r.y += ostate.r.y
|
||||
bstate.r.z += ostate.r.z
|
||||
bstate.v.x += ostate.v.x
|
||||
bstate.v.y += ostate.v.y
|
||||
bstate.v.z += ostate.v.z
|
||||
}
|
||||
}
|
||||
|
||||
// Calculate the net acceleration experienced by the small bodies.
|
||||
calcBodyAccelerations()
|
||||
|
||||
// To prepare for a possible swap operation, duplicate the current state into the previous state.
|
||||
duplicate()
|
||||
}
|
||||
|
||||
/**
|
||||
* Advances a gravity simulation by a small time step.
|
||||
*
|
||||
* Updates the simulation of the user-supplied small bodies
|
||||
* to the time indicated by the `time` parameter.
|
||||
* Retuns an updated array of state vectors for the small bodies.
|
||||
* The positions and velocities in the returned array are referenced
|
||||
* to the `originBody` that was used to construct this simulator.
|
||||
*
|
||||
* @param 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.
|
||||
* Some experimentation may be necessary to find a good value.
|
||||
* The `time` value may be after or before the current simulation time
|
||||
* to move forward or backward in time.
|
||||
*/
|
||||
fun update(time: Time): Array<StateVector> {
|
||||
val dt = time.tt - curr.time.tt
|
||||
if (dt == 0.0) {
|
||||
// Special case: the time has not changed, so skip the usual physics calculations.
|
||||
// This allows a way for the caller to query the current state if desired.
|
||||
// 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.
|
||||
duplicate()
|
||||
} else {
|
||||
// Swap the current state and the previous state. Then calculate the new state.
|
||||
swap()
|
||||
|
||||
// Update the current time.
|
||||
curr.time = time
|
||||
|
||||
// Now that the time is set, it is safe to update the solar system.
|
||||
calcSolarSystem()
|
||||
|
||||
// Estimate the position of each small body as if the current
|
||||
// acceleration applies across the whole time interval.
|
||||
prev.bodies.forEachIndexed { i, p ->
|
||||
curr.bodies[i].r = updatePosition(dt, p.r, p.v, p.a)
|
||||
}
|
||||
|
||||
// Calculate the acceleration experienced by the small bodies
|
||||
// at their respective approximate next locations.
|
||||
calcBodyAccelerations()
|
||||
|
||||
prev.bodies.forEachIndexed { i, p ->
|
||||
val c = curr.bodies[i]
|
||||
|
||||
// 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.
|
||||
val 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.
|
||||
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 BodyGravCalc to the public
|
||||
// type StateVector.
|
||||
val bodyStateArray = curr.bodies.map { exportGravCalc(it, time) }.toTypedArray()
|
||||
|
||||
if (originBody != Body.SSB) {
|
||||
// Now we have to convert the coordinate system to the caller's chosen origin body.
|
||||
// Determine the barycentric state of the origin body.
|
||||
val ostate = bodyState(originBody)
|
||||
|
||||
// Subtract vectors to convert barycentric states to origin-centric states.
|
||||
for (i in bodyStateArray.indices) {
|
||||
bodyStateArray[i] = StateVector(
|
||||
bodyStateArray[i].x - ostate.r.x,
|
||||
bodyStateArray[i].y - ostate.r.y,
|
||||
bodyStateArray[i].z - ostate.r.z,
|
||||
bodyStateArray[i].vx - ostate.v.x,
|
||||
bodyStateArray[i].vy - ostate.v.y,
|
||||
bodyStateArray[i].vz - ostate.v.z,
|
||||
time
|
||||
)
|
||||
}
|
||||
}
|
||||
|
||||
return bodyStateArray
|
||||
}
|
||||
|
||||
/**
|
||||
* 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].
|
||||
*/
|
||||
fun swap() {
|
||||
val s = curr
|
||||
curr = prev
|
||||
prev = s
|
||||
}
|
||||
|
||||
private fun bodyState(body: Body): BodyState =
|
||||
// Return a reference to the object where we cache the barycentric
|
||||
// state of a given body, or fail if this body is not valid.
|
||||
if (body == Body.Sun || (body.ordinal >= Body.Mercury.ordinal && body.ordinal <= Body.Neptune.ordinal))
|
||||
curr.gravitators[body.ordinal]
|
||||
else
|
||||
throw InvalidBodyException(body)
|
||||
|
||||
private fun initialEndpoint(time: Time, bodyStates: List<StateVector>): GravSimEndpoint {
|
||||
// Create an initial "endpoint" data structure, but without correcting
|
||||
// the coordinate system yet for the caller's chosen originBody.
|
||||
// That has to wait until after we have called `calcSolarSystem`
|
||||
// to know the positions of the planets.
|
||||
// We also make placeholders for the Sun and planet body states,
|
||||
// which are calculated later by `calcSolarSystem`.
|
||||
|
||||
val gravitators = Array<BodyState>(Body.Sun.ordinal + 1) {
|
||||
BodyState(
|
||||
time.tt,
|
||||
@@ -7932,6 +8136,47 @@ class GravitySimulator {
|
||||
sun.r.negate()
|
||||
sun.v.negate()
|
||||
}
|
||||
|
||||
private fun calcBodyAccelerations() {
|
||||
// Calculate the gravitational acceleration experienced by the simulated bodies.
|
||||
for (calc in curr.bodies) {
|
||||
calc.a.setToZero()
|
||||
addAcceleration(calc.a, calc.r, SUN_GM, curr.gravitators[Body.Sun.ordinal ].r)
|
||||
addAcceleration(calc.a, calc.r, MERCURY_GM, curr.gravitators[Body.Mercury.ordinal].r)
|
||||
addAcceleration(calc.a, calc.r, VENUS_GM, curr.gravitators[Body.Venus.ordinal ].r)
|
||||
addAcceleration(calc.a, calc.r, EARTH_GM + MOON_GM, curr.gravitators[Body.Earth.ordinal ].r)
|
||||
addAcceleration(calc.a, calc.r, MARS_GM, curr.gravitators[Body.Mars.ordinal ].r)
|
||||
addAcceleration(calc.a, calc.r, JUPITER_GM, curr.gravitators[Body.Jupiter.ordinal].r)
|
||||
addAcceleration(calc.a, calc.r, SATURN_GM, curr.gravitators[Body.Saturn.ordinal ].r)
|
||||
addAcceleration(calc.a, calc.r, URANUS_GM, curr.gravitators[Body.Uranus.ordinal ].r)
|
||||
addAcceleration(calc.a, calc.r, NEPTUNE_GM, curr.gravitators[Body.Neptune.ordinal].r)
|
||||
}
|
||||
}
|
||||
|
||||
private fun addAcceleration(acc: TerseVector, smallPos: TerseVector, gm: Double, majorPos: TerseVector) {
|
||||
val dx = majorPos.x - smallPos.x
|
||||
val dy = majorPos.y - smallPos.y
|
||||
val dz = majorPos.z - smallPos.z
|
||||
val r2 = dx*dx + dy*dy + dz*dz
|
||||
val pull = gm / (r2 * sqrt(r2))
|
||||
acc.x += dx * pull
|
||||
acc.y += dy * pull
|
||||
acc.z += dz * pull
|
||||
}
|
||||
|
||||
private fun duplicate() {
|
||||
// Copy the current state into the previous state, so that both become the same moment in time.
|
||||
|
||||
prev.time = curr.time
|
||||
|
||||
for (i in curr.gravitators.indices) {
|
||||
prev.gravitators[i].copyFrom(curr.gravitators[i])
|
||||
}
|
||||
|
||||
for (i in curr.bodies.indices) {
|
||||
prev.bodies[i].copyFrom(curr.bodies[i])
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
@@ -3743,12 +3743,12 @@ static void CalcSolarSystem(astro_grav_sim_t *sim)
|
||||
static void CalcBodyAccelerations(astro_grav_sim_t *sim)
|
||||
{
|
||||
int i;
|
||||
const body_state_t *grav = sim->curr->gravitators;
|
||||
|
||||
/* Calculate the gravitational acceleration experienced by the simulated bodies. */
|
||||
for (i = 0; i < sim->numBodies; ++i)
|
||||
{
|
||||
body_grav_calc_t *calc = &sim->curr->bodies[i];
|
||||
const body_state_t *grav = sim->curr->gravitators;
|
||||
|
||||
calc->a = VecZero;
|
||||
|
||||
@@ -3990,7 +3990,7 @@ fail:
|
||||
|
||||
|
||||
/**
|
||||
* @brief Advances a gravity simulation of a small body by a small time step.
|
||||
* @brief Advances a gravity simulation by a small time step.
|
||||
*
|
||||
* @param sim
|
||||
* A simulation object that was created by a prior call to #Astronomy_GravSimInit.
|
||||
|
||||
@@ -14,6 +14,13 @@ A gravity simulator object simulates a series of incremental time steps, calcula
|
||||
|---|---|
|
||||
| [GravitySimulator](-gravity-simulator.md)<br>fun [GravitySimulator](-gravity-simulator.md)(originBody: [Body](../-body/index.md), time: [Time](../-time/index.md), bodyStates: [List](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin.collections/-list/index.html)<[StateVector](../-state-vector/index.md)>)<br>Creates a gravity simulation object. |
|
||||
|
||||
## Functions
|
||||
|
||||
| Name | Summary |
|
||||
|---|---|
|
||||
| [swap](swap.md)<br>fun [swap](swap.md)()<br>Exchange the current time step with the previous time step. |
|
||||
| [update](update.md)<br>fun [update](update.md)(time: [Time](../-time/index.md)): [Array](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-array/index.html)<[StateVector](../-state-vector/index.md)><br>Advances a gravity simulation by a small time step. |
|
||||
|
||||
## Properties
|
||||
|
||||
| Name | Summary |
|
||||
|
||||
15
source/kotlin/doc/-gravity-simulator/swap.md
Normal file
15
source/kotlin/doc/-gravity-simulator/swap.md
Normal file
@@ -0,0 +1,15 @@
|
||||
//[astronomy](../../../index.md)/[io.github.cosinekitty.astronomy](../index.md)/[GravitySimulator](index.md)/[swap](swap.md)
|
||||
|
||||
# swap
|
||||
|
||||
fun [swap](swap.md)()
|
||||
|
||||
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](update.md). 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](update.md).
|
||||
15
source/kotlin/doc/-gravity-simulator/update.md
Normal file
15
source/kotlin/doc/-gravity-simulator/update.md
Normal file
@@ -0,0 +1,15 @@
|
||||
//[astronomy](../../../index.md)/[io.github.cosinekitty.astronomy](../index.md)/[GravitySimulator](index.md)/[update](update.md)
|
||||
|
||||
# update
|
||||
|
||||
fun [update](update.md)(time: [Time](../-time/index.md)): [Array](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-array/index.html)<[StateVector](../-state-vector/index.md)>
|
||||
|
||||
Advances a gravity simulation by a small time step.
|
||||
|
||||
Updates the simulation of the user-supplied small bodies to the time indicated by the time parameter. Retuns an updated array of state vectors for the small bodies. The positions and velocities in the returned array are referenced to the originBody that was used to construct this simulator.
|
||||
|
||||
## Parameters
|
||||
|
||||
| | |
|
||||
|---|---|
|
||||
| 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. Some experimentation may be necessary to find a good value. The time value may be after or before the current simulation time to move forward or backward in time. |
|
||||
@@ -705,6 +705,12 @@ internal data class TerseVector(var x: Double, var y: Double, var z: Double) {
|
||||
z = -z
|
||||
}
|
||||
|
||||
fun copyFrom(other: TerseVector) {
|
||||
x = other.x
|
||||
y = other.y
|
||||
z = other.z
|
||||
}
|
||||
|
||||
companion object {
|
||||
@JvmStatic
|
||||
fun zero() = TerseVector(0.0, 0.0, 0.0)
|
||||
@@ -2944,10 +2950,21 @@ internal class BodyState(
|
||||
val r: TerseVector,
|
||||
val v: TerseVector
|
||||
) {
|
||||
fun increment(other: BodyState) {
|
||||
r.increment(other.r)
|
||||
v.increment(other.v)
|
||||
}
|
||||
|
||||
fun decrement(other: BodyState) {
|
||||
r.decrement(other.r)
|
||||
v.decrement(other.v)
|
||||
}
|
||||
|
||||
fun copyFrom(other: BodyState) {
|
||||
tt = other.tt
|
||||
r.copyFrom(other.r)
|
||||
v.copyFrom(other.v)
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@@ -2959,6 +2976,14 @@ private fun exportState(bodyState: BodyState, time: Time) =
|
||||
)
|
||||
|
||||
|
||||
private fun exportGravCalc(calc: BodyGravCalc, time: Time) =
|
||||
StateVector(
|
||||
calc.r.x, calc.r.y, calc.r.z,
|
||||
calc.v.x, calc.v.y, calc.v.z,
|
||||
time
|
||||
)
|
||||
|
||||
|
||||
private fun calcVsopPosVel(model: VsopModel, tt: Double): BodyState {
|
||||
val t = tt / DAYS_PER_MILLENNIUM
|
||||
|
||||
@@ -3231,11 +3256,18 @@ private class MoonContext(time: Time) {
|
||||
// Pluto/SSB gravitation simulator
|
||||
|
||||
private class BodyGravCalc(
|
||||
val tt: Double, // J2000 terrestrial time [days]
|
||||
val r: TerseVector, // position [au]
|
||||
val v: TerseVector, // velocity [au/day]
|
||||
val a: TerseVector // acceleration [au/day]
|
||||
)
|
||||
var tt: Double, // J2000 terrestrial time [days]
|
||||
var r: TerseVector, // position [au]
|
||||
var v: TerseVector, // velocity [au/day]
|
||||
var a: TerseVector // acceleration [au/day]
|
||||
) {
|
||||
fun copyFrom(other: BodyGravCalc) {
|
||||
tt = other.tt
|
||||
r.copyFrom(other.r)
|
||||
v.copyFrom(other.v)
|
||||
a.copyFrom(other.a)
|
||||
}
|
||||
}
|
||||
|
||||
private class MajorBodies(
|
||||
val sun: BodyState,
|
||||
@@ -7872,17 +7904,189 @@ class GravitySimulator {
|
||||
|
||||
// Verify that all the state vectors are valid and have matching times.
|
||||
for (b in bodyStates) {
|
||||
if (b.t.tt != time.tt)
|
||||
if (b.t.tt != time.tt) {
|
||||
throw IllegalArgumentException("Inconsistent time(s) in bodyStates array.")
|
||||
}
|
||||
}
|
||||
|
||||
curr = initialEndpoint(time, bodyStates)
|
||||
prev = initialEndpoint(time, bodyStates)
|
||||
|
||||
// Calculate the states of the Sun and planets.
|
||||
calcSolarSystem()
|
||||
|
||||
// We need to do all the physics calculations in barycentric coordinates.
|
||||
// But the caller provides the input vectors with respect to `originBody`.
|
||||
// Correct the input body state vectors for the specified origin.
|
||||
if (originBody != Body.SSB) {
|
||||
// Determine the barycentric state of the origin body.
|
||||
val ostate = bodyState(originBody)
|
||||
|
||||
// Add barycentric origin to origin-centric bodies to obtain barycentric bodies.
|
||||
for (bstate in curr.bodies) {
|
||||
bstate.r.x += ostate.r.x
|
||||
bstate.r.y += ostate.r.y
|
||||
bstate.r.z += ostate.r.z
|
||||
bstate.v.x += ostate.v.x
|
||||
bstate.v.y += ostate.v.y
|
||||
bstate.v.z += ostate.v.z
|
||||
}
|
||||
}
|
||||
|
||||
// Calculate the net acceleration experienced by the small bodies.
|
||||
calcBodyAccelerations()
|
||||
|
||||
// To prepare for a possible swap operation, duplicate the current state into the previous state.
|
||||
duplicate()
|
||||
}
|
||||
|
||||
/**
|
||||
* Advances a gravity simulation by a small time step.
|
||||
*
|
||||
* Updates the simulation of the user-supplied small bodies
|
||||
* to the time indicated by the `time` parameter.
|
||||
* Retuns an updated array of state vectors for the small bodies.
|
||||
* The positions and velocities in the returned array are referenced
|
||||
* to the `originBody` that was used to construct this simulator.
|
||||
*
|
||||
* @param 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.
|
||||
* Some experimentation may be necessary to find a good value.
|
||||
* The `time` value may be after or before the current simulation time
|
||||
* to move forward or backward in time.
|
||||
*/
|
||||
fun update(time: Time): Array<StateVector> {
|
||||
val dt = time.tt - curr.time.tt
|
||||
if (dt == 0.0) {
|
||||
// Special case: the time has not changed, so skip the usual physics calculations.
|
||||
// This allows a way for the caller to query the current state if desired.
|
||||
// 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.
|
||||
duplicate()
|
||||
} else {
|
||||
// Swap the current state and the previous state. Then calculate the new state.
|
||||
swap()
|
||||
|
||||
// Update the current time.
|
||||
curr.time = time
|
||||
|
||||
// Now that the time is set, it is safe to update the solar system.
|
||||
calcSolarSystem()
|
||||
|
||||
// Estimate the position of each small body as if the current
|
||||
// acceleration applies across the whole time interval.
|
||||
prev.bodies.forEachIndexed { i, p ->
|
||||
curr.bodies[i].r = updatePosition(dt, p.r, p.v, p.a)
|
||||
}
|
||||
|
||||
// Calculate the acceleration experienced by the small bodies
|
||||
// at their respective approximate next locations.
|
||||
calcBodyAccelerations()
|
||||
|
||||
prev.bodies.forEachIndexed { i, p ->
|
||||
val c = curr.bodies[i]
|
||||
|
||||
// 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.
|
||||
val 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.
|
||||
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 BodyGravCalc to the public
|
||||
// type StateVector.
|
||||
val bodyStateArray = curr.bodies.map { exportGravCalc(it, time) }.toTypedArray()
|
||||
|
||||
if (originBody != Body.SSB) {
|
||||
// Now we have to convert the coordinate system to the caller's chosen origin body.
|
||||
// Determine the barycentric state of the origin body.
|
||||
val ostate = bodyState(originBody)
|
||||
|
||||
// Subtract vectors to convert barycentric states to origin-centric states.
|
||||
for (i in bodyStateArray.indices) {
|
||||
bodyStateArray[i] = StateVector(
|
||||
bodyStateArray[i].x - ostate.r.x,
|
||||
bodyStateArray[i].y - ostate.r.y,
|
||||
bodyStateArray[i].z - ostate.r.z,
|
||||
bodyStateArray[i].vx - ostate.v.x,
|
||||
bodyStateArray[i].vy - ostate.v.y,
|
||||
bodyStateArray[i].vz - ostate.v.z,
|
||||
time
|
||||
)
|
||||
}
|
||||
}
|
||||
|
||||
return bodyStateArray
|
||||
}
|
||||
|
||||
/**
|
||||
* 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].
|
||||
*/
|
||||
fun swap() {
|
||||
val s = curr
|
||||
curr = prev
|
||||
prev = s
|
||||
}
|
||||
|
||||
private fun bodyState(body: Body): BodyState =
|
||||
// Return a reference to the object where we cache the barycentric
|
||||
// state of a given body, or fail if this body is not valid.
|
||||
if (body == Body.Sun || (body.ordinal >= Body.Mercury.ordinal && body.ordinal <= Body.Neptune.ordinal))
|
||||
curr.gravitators[body.ordinal]
|
||||
else
|
||||
throw InvalidBodyException(body)
|
||||
|
||||
private fun initialEndpoint(time: Time, bodyStates: List<StateVector>): GravSimEndpoint {
|
||||
// Create an initial "endpoint" data structure, but without correcting
|
||||
// the coordinate system yet for the caller's chosen originBody.
|
||||
// That has to wait until after we have called `calcSolarSystem`
|
||||
// to know the positions of the planets.
|
||||
// We also make placeholders for the Sun and planet body states,
|
||||
// which are calculated later by `calcSolarSystem`.
|
||||
|
||||
val gravitators = Array<BodyState>(Body.Sun.ordinal + 1) {
|
||||
BodyState(
|
||||
time.tt,
|
||||
@@ -7932,6 +8136,47 @@ class GravitySimulator {
|
||||
sun.r.negate()
|
||||
sun.v.negate()
|
||||
}
|
||||
|
||||
private fun calcBodyAccelerations() {
|
||||
// Calculate the gravitational acceleration experienced by the simulated bodies.
|
||||
for (calc in curr.bodies) {
|
||||
calc.a.setToZero()
|
||||
addAcceleration(calc.a, calc.r, SUN_GM, curr.gravitators[Body.Sun.ordinal ].r)
|
||||
addAcceleration(calc.a, calc.r, MERCURY_GM, curr.gravitators[Body.Mercury.ordinal].r)
|
||||
addAcceleration(calc.a, calc.r, VENUS_GM, curr.gravitators[Body.Venus.ordinal ].r)
|
||||
addAcceleration(calc.a, calc.r, EARTH_GM + MOON_GM, curr.gravitators[Body.Earth.ordinal ].r)
|
||||
addAcceleration(calc.a, calc.r, MARS_GM, curr.gravitators[Body.Mars.ordinal ].r)
|
||||
addAcceleration(calc.a, calc.r, JUPITER_GM, curr.gravitators[Body.Jupiter.ordinal].r)
|
||||
addAcceleration(calc.a, calc.r, SATURN_GM, curr.gravitators[Body.Saturn.ordinal ].r)
|
||||
addAcceleration(calc.a, calc.r, URANUS_GM, curr.gravitators[Body.Uranus.ordinal ].r)
|
||||
addAcceleration(calc.a, calc.r, NEPTUNE_GM, curr.gravitators[Body.Neptune.ordinal].r)
|
||||
}
|
||||
}
|
||||
|
||||
private fun addAcceleration(acc: TerseVector, smallPos: TerseVector, gm: Double, majorPos: TerseVector) {
|
||||
val dx = majorPos.x - smallPos.x
|
||||
val dy = majorPos.y - smallPos.y
|
||||
val dz = majorPos.z - smallPos.z
|
||||
val r2 = dx*dx + dy*dy + dz*dz
|
||||
val pull = gm / (r2 * sqrt(r2))
|
||||
acc.x += dx * pull
|
||||
acc.y += dy * pull
|
||||
acc.z += dz * pull
|
||||
}
|
||||
|
||||
private fun duplicate() {
|
||||
// Copy the current state into the previous state, so that both become the same moment in time.
|
||||
|
||||
prev.time = curr.time
|
||||
|
||||
for (i in curr.gravitators.indices) {
|
||||
prev.gravitators[i].copyFrom(curr.gravitators[i])
|
||||
}
|
||||
|
||||
for (i in curr.bodies.indices) {
|
||||
prev.bodies[i].copyFrom(curr.bodies[i])
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user