Orbital Aero Model: Indirect Forces

The orbital aero model calculates two forces I'm calling "Indirect" - Gravity and Atmospheric Drag.

Gravity

The force of gravity is dependent on the mass of the two bodies and distance between them. The simulation calculates the influence of gravity from every object to every other object every frame. Not only will satellites be pulled to the Earth, the Earth will be pulled to the satellites!

Fg is the force of gravity, in newtons.
G is the gravitational constant, 6.67408 x 10^-11 in m^3/(kg * s^2).
m1 is one body, in kilograms.
m2 is another body, in kilograms.
r is the distance between the bodies, in meters.

https://en.wikipedia.org/wiki/Gravity

Atmospheric Drag

Some bodies have an atmosphere modeled around them. For any body within the atmosphere of another body, the drag force is calculated and applied.

Fd is the force of drag, in newtons.

rho is the density of the fluid, in kg / m^3.
A is the drag reference area, in m^2.
Cd is the coefficient of drag, which is unitless.
v is the velocity relative to the object, in m/s.

https://en.wikipedia.org/wiki/Drag_(physics)

Orbital Aero Code - update()

To recap from a previous post, there is an update() function in the OrbitalAeroModel class:

// Primary function to process the physics and kinematics all bodies.
bool OrbitalAeroModel::update(double timeStep_seconds)
{
 // It's important all forces are calculated before all positions are updated before all collisions
 // are checked, etc. That way the order in which bodies are looped through does not affect how they
 // affect other bodies in space.

 // First calculate indirect forces, such as those from gravity and atmospheric drag.
 // Currently includes the thrust of engines, but that may be separated eventually.
 processIndirectForces();

 // Next add frictional forces.
 // This must be done after gravity/atmospheric forces to know the normal force for friction.
 processDirectForces();

 // Can now translate and rotate each body.
 processKinematics(timeStep_seconds);

 // After moving the body, check for collisions and adjust positions/velocities as appropriate.
 processCollisions(timeStep_seconds);

 return true;
}

Orbital Aero Code - processIndirectForces()

// Function to loop through all bodies and calculate indirect forces from the gravity
// and atmosphere of other bodies.
bool OrbitalAeroModel::processIndirectForces()
{
// Loop through all bodies and calculate their forces upon each other.

// The outer loop is for the primary body that we're calculating forces for.
map <uint64, BodyModel>::iterator iteratorPrimaryBody = tableBodies.begin();
while (iteratorPrimaryBody != tableBodies.end())
{
BodyModel &primaryBody = iteratorPrimaryBody->second;

// Only process the primary body if it's active.
if (primaryBody.active)
{
// Create a vector to contain the sum of all forces.
Vector sumForces_newtons;

// Start with the thrust force of any onboard engines.
sumForces_newtons = primaryBody.forceThrust_newtons;

// The inner loop is for all other bodies that influence the primary.
map <uint64, BodyModel>::iterator iteratorSecondaryBody = tableBodies.begin();
while (iteratorSecondaryBody != tableBodies.end())
{
BodyModel &secondaryBody = iteratorSecondaryBody->second;

// Only contribute the secondary body's influence if it's active and not the primary body.
if (secondaryBody.active && (&primaryBody != &secondaryBody))
{
// Get the relative distance between the two bodies.
Vector relativeDistance_meters = primaryBody.location_meters - secondaryBody.location_meters;
double relativeDistanceMagnitude_meters = relativeDistance_meters.magnitude();

if (relativeDistanceMagnitude_meters > 0.0)
{
// Calculate the force of gravity.
// Fg = GMm/r^2
// G = Gravitational Constant
// M = Mass of first body
// m = Mass of second body
// r = Distance between bodies
double secondayGravitationalForce_newtons = Constants::gravitationalConstant_meters3PerKilogramSecond2 * (primaryBody.massTotal_kilograms() * secondaryBody.massTotal_kilograms()) / 
relativeDistanceMagnitude_meters / relativeDistanceMagnitude_meters;

// Add passive body gravitational influence to the total force.
sumForces_newtons += -relativeDistance_meters.unit() * secondayGravitationalForce_newtons;

// Are we within the atmosphere of the secondary body?
double altitudeMSL_meters = relativeDistanceMagnitude_meters - secondaryBody.radius_meters;
if (altitudeMSL_meters < secondaryBody.atmosphereMaximumAltitudeMSL_meters())
{
// Get the air density at our current altitude.
double atmosphereDensityMSL_kilogramsPerMeter3 = secondaryBody.getAtmosphereDensityMSL_kilogramsPerMeter3(altitudeMSL_meters);

// Calculate and add passive body atmospheric forces. (Should implement 3-D drag & area.)
// Fd = 0.5 * rho * A * Cd * v * v
// rho = Density of Fluid
// Cd = Coefficient of Drag
// A = Reference Area
// v = Relative velocity through fluid
Vector relativeVelocity_metersPerSecond = primaryBody.velocity_metersPerSecond - secondaryBody.velocity_metersPerSecond;
sumForces_newtons += -relativeVelocity_metersPerSecond.sign() * relativeVelocity_metersPerSecond * relativeVelocity_metersPerSecond * 0.5 * atmosphereDensityMSL_kilogramsPerMeter3 * primaryBody.dragCoefficient * primaryBody.dragReferenceArea_meters2;
}
} // relativeDistanceMagnitude_meters > 0.0

} // ((secondaryBody.active) && (primaryBody.index != secondaryBody.index))
iteratorSecondaryBody++;
} // (iteratorSecondaryBody != tableBodies.end())
primaryBody.forceIndirect_newtons = sumForces_newtons;
} // (primaryBody.active)
iteratorPrimaryBody++;
} // (iteratorPrimaryBody != tableBodies.end())
return true;
}

Results

A couple screenshots from the simulation...

The moon satellite is orbiting the moon while the moon is orbiting the Earth. The smaller satellites are contributing gravitational effects on the planet and moon, but their influence is so small they do practically nothing.

The two satellites closest to the Earth are within the atmosphere and are influenced by atmospheric drag. The lower the altitude, the higher the air density, and thus the higher drag force.

The code as-is supports the concept of a terminal velocity. A body dropped to another body with an atmosphere will stop accelerating when its force of gravity equals the force of drag.

Copyright (c) 2016 Clinton Kam
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.