I have set up Newtonian Gravity in my Game Engine, allowing me to simulate the gravitational attraction between celestial bodies.
I have the following variables defined:
btScalar whiteHoleMass = 5e+28;
btScalar whiteHoleRadius = 500'000.0;
btScalar sunRadius = 100'000;
btScalar sunOrbitRadius = 50000;
btScalar sunMass = 5e+24;
btScalar sunOrbitTime = 1800; // Orbit time is in seconds
btScalar planetRadius = 50'000;
btScalar planetMass = 5e+21;
btScalar planetOrbitTime = 900;
I add the white hole at position: {0, whiteHoleRadius, 0}
I add the sun at position: {sunOrbitRadius, whiteHoleRadius, 0}
I add the planet at position: {planetOrbitRadius + sunOrbitRadius, whiteHoleRadius, 0}
I add 10 orbiting satellites at position: {thisOrbitalHeight + planetOrbitRadius + sunOrbitRadius, whiteHoleRadius, 0}
sunOrbitRadius, planetOrbitRadius and thisOrbitalHeight are calculated with the following functions:
btScalar sunOrbitRadius = 0;
btScalar sunOrbitVelocity = 0;
calculateOrbit(sunOrbitTime, whiteHoleMass, sunOrbitRadius, sunOrbitVelocity);
btScalar planetOrbitRadius = 0;
btScalar planetOrbitVelocity = 0;
calculateOrbit(planetOrbitTime, sunMass, planetOrbitRadius, planetOrbitVelocity);
btScalar orbitalSatelliteOrbitVelocity = 0;
btScalar thisOrbitalHeight = 0;
calculateOrbit(currentOrbitalLength, planetMass, thisOrbitalHeight, orbitalSatelliteOrbitVelocity);
and the calculateOrbit
function is defined as:
void calculateOrbit(const btScalar &orbitTimeSecs, const btScalar ¢ralMassKg, btScalar &orbitRadius, btScalar &orbitVelocity)
{
orbitRadius = std::cbrt((Common::G * centralMassKg * orbitTimeSecs * orbitTimeSecs) / (4 * Common::Pi * Common::Pi));
orbitVelocity = std::sqrt(Common::G * centralMassKg / orbitRadius);
};
I then set the Linear Velocity of the Sun, Planet and Satellites like so:
sunPointer->rigidBody->setLinearVelocity({0, 0, sunOrbitVelocity});
planetPointer->rigidBody->setLinearVelocity({0, 0, planetOrbitVelocity + sunOrbitVelocity});
orbitalPointer->rigidBody->setLinearVelocity({0, 0, orbitalSatelliteOrbitVelocity + planetOrbitVelocity + sunOrbitVelocity});
The linear velocity is added in the Z dimension.
I would have thought this was sufficient to keep Satellite/Planet/Sun in constant orbit, however 45 seconds into running the simulation the satellites fall out of orbit and crash into the planet.
How can I set up orbits of multiple bodies so they are stable?
I calculate gravity between objects like so:
for (uInteger64 i = 1; i <= entityMap._size; i++)
{
btVector3 accelerations(0, 0, 0);
auto &entityA = *entityMap._data[i]->value.pointer;
auto &bodyA = *entityA.rigidBody;
auto bodyAMass = bodyA.getMass();
auto &positionA = bodyA.getCenterOfMassPosition();
if (bodyA.getActivationState() == ACTIVE_TAG)
{
for (uInteger64 j = 1; j <= entityMap._size; j++)
{
if (i == j)
{
continue; // Skip self
}
auto &entityB = *entityMap._data[j]->value.pointer;
auto &bodyB = *entityB.rigidBody;
auto bodyBMass = bodyB.getMass();
auto &positionB = bodyB.getCenterOfMassPosition();
btVector3 vectorAB = positionB - positionA;
btScalar distanceAB = vectorAB.length();
btVector3 directionAB = vectorAB.normalized();
Floating128 forceMagnitude = (Floating128)Common::G * ((Floating128)bodyAMass * (Floating128)bodyBMass) / ((Floating128)distanceAB * (Floating128)distanceAB);
btVector3 forceAB = directionAB * forceMagnitude;
accelerations += forceAB / bodyAMass;
}
entityA.applyCustomAcceleration(accelerations);
}
}
and applyCustomAcceleration is defined as:
void IEntity::applyCustomAcceleration(const btVector3 &acceleration)
{
if (!rigidBody || rigidBody->isStaticOrKinematicObject())
{
return;
}
if (rigidBody)
{
btVector3 force = acceleration * rigidBody->getMass();
rigidBody->applyCentralForce(force);
}
};