I’m within the midst of creating a gravitational physique demo for my physics engine. For the demo, I’m making an attempt to simulate each all of the interior planets’ (Mercury-Venus-Earth-Mars) orbits across the solar, in addition to the orbit of the moon across the Earth. The planets are all working as anticipated, nevertheless the Moon doesn’t appear to need to orbit the Earth. I’m at an entire loss for why that is. I am fairly certain it is not that the gravity is calculated incorrectly, in any other case the planets wouldn’t all accurately orbit the solar. The calculations to find out the drive appearing on the moon start slightly below the whereas(!WindowShouldClose()) loop. Under is a (moderately prolonged) MCVE to see what I imply. I am hoping that it’s a conceptual difficulty about how I’ve arrange my preliminary situations greater than a coding mistake. Please word that it requires an set up of raylib
At present, the moon begins with the proper preliminary velocity relative to the solar (the mix of the Earth’s velocity relative to the Solar and the moon’s velocity relative to the Earth), nevertheless that velocity by no means considerably modifications, so the moon finally ends up getting into mainly a straight line endlessly. The “mainly” is essential. The Moon is being affected by the Earth, simply on a a lot a lot smaller scale than anticipated
To make use of the simulation, there are the usual WASD keys, and the mouse to go searching.
#embrace "raylib.h"
#embrace "rlgl.h"
#embrace <math.h>
#embrace <vector>
#outline real_sqrt sqrt
#outline real_pow powf
// Set some bodily values for the simulation
#outline G 6.67408e-11 // [m^3 kg^-1 s^-2]
#outline SOLARMASS 1.989e30 // [kg]
#outline SOLARRADIUS 6.957e8 // [m]
#outline EARTHMASS 5.97219e24 // [kg]
#outline EARTHORBIT 1.496e11 // [m]
#outline EARTHVELOCITY 2.978e4 // [m/s] (Relative to Solar)
#outline EARTHRADIUS 6.371e6 // [m]
#outline MOONMASS 7.34767309e22 // [kg]
#outline MOONORBIT 3.844e8 // [m]
#outline MOONVELOCITY 1.023e3 // [m/s] (Relative to Earth)
#outline MOONRADIUS 1.737e6 // [m]
utilizing namespace std;
namespace engine {
typedef float actual;
class Vec3 {
public:
// Spatial coordinates
actual x;
actual y;
actual z;
personal:
// Padding to make sure 4 phrase alignment
actual pad;
public:
// Default constructor creates a 0 vector
Vec3() : x(0), y(0), z(0) {}
// Constructor for when values are handed
Vec3(const actual x, const actual y, const actual z) : x(x), y(y), z(z) {}
// Get magnitude of vector
actual magnitude() const {
return real_sqrt(x*x + y*y + z*z);
}
// Generally it's helpful and quicker to only have the sq. of the magnitude
actual squareMagnitude() const {
return x * x + y * y + z * z;
}
// Normalize a non-zero vector
void normalize() {
actual l = magnitude();
if (l > 0) {
(*this) *= ((actual)1)/l;
}
}
// Multiplies vector by given scalar
void operator *= (const actual scalar) {
x *= scalar;
y *= scalar;
z *= scalar;
}
// Returns vector scaled by worth
Vec3 operator * (const actual worth) const {
return Vec3(x * worth, y * worth, z * worth);
}
// Provides given vector
Vec3 operator + (const Vec3& v) const {
return Vec3(x + v.x, y + v.y, z + v.z);
}
void operator += (const Vec3& v) {
x += v.x;
y += v.y;
z += v.z;
}
// Subtracts given vector
void operator -= (const Vec3& v) {
x -= v.x;
y -= v.y;
z -= v.z;
}
Vec3 operator - (const Vec3& v) const {
return Vec3(x - v.x, y - v.y, z - v.z);
}
// Provides a given scaled vector
void addScaledVector(const Vec3& v, actual scale) {
x += (v.x * scale);
y += (v.y * scale);
z += (v.z * scale);
}
actual operator * (const Vec3& v) const {
return x * v.x + y * v.y + z * v.z;
}
void operator = (const Vec3& v) {
x = v.x;
y = v.y;
z = v.z;
}
void clear() {
x = y = z = 0;
}
}; // class Vec3
class Particle {
protected:
// Maintain observe of place and its time derivatives
Vec3 pos;
Vec3 vel;
Vec3 acc;
// Maintain observe of damping of linear movement
actual damping;
/**
* Holds the inverse of the mass of the particle. It
* is extra helpful to carry the inverse mass as a result of
* integration is less complicated, and since in real-time
* simulation it's extra helpful to have objects with
* infinite mass (immovable) than zero mass
* (utterly unstable in numerical simulation).
*/
actual inverseMass;
Vec3 forceAccum;
public:
Particle() : pos(0, 0, 0), vel(0, 0, 0), acc(0, 0, 0), damping((actual)1.0), inverseMass(0) {};
Particle(const Vec3 pos,
const Vec3 vel,
const Vec3 acc,
const actual damping,
const actual inverseMass) : pos(pos), vel(vel), acc(acc), damping(damping), inverseMass(inverseMass) {};
void combine(actual period) {
// We can't combine particles with infinite or adverse mass
if (inverseMass <= 0.0f) return;
assert(period > 0.0);
// Replace linear place.
pos.addScaledVector(vel, period);
pos.addScaledVector(acc, period * period * 0.5);
//Vec3 resultingAcc = acc;
//resultingAcc.addScaledVector(forceAccum, inverseMass);
acc = forceAccum * inverseMass;
// Impose drag.
vel *= real_pow(damping, period);
// Replace linear velocity from the acceleration.
//vel.addScaledVector(resultingAcc, period);
vel.addScaledVector(acc, period);
// Clear forces
clearAccumulator();
} // combine()
// Return place of particle
Vec3 getPosition() const {
return pos;
}
void clearAccumulator() {
forceAccum.clear();
}
void addForce(const Vec3& f) {
forceAccum += f;
}
bool hasFiniteMass() const {
return inverseMass > 0.0f;
}
actual getMass() const {
if (inverseMass == 0) {
return std::numeric_limits<float>::max();
} else {
return ((actual)1.0)/inverseMass;
}
}
}; // class Particle
class ParticleForceGenerator {
public:
digital void updateForce(Particle *particle, actual period) = 0;
};
class ParticleForceRegistry {
protected:
// Maintain observe of drive generator and the particle that referred to as for it
struct ParticleForceRegistration {
Particle *particle;
ParticleForceGenerator *fg;
};
// Holds record of registrations
typedef std::vector<ParticleForceRegistration> Registry;
Registry registrations;
public:
// Registers the given drive generator to use to the given particle
void add(Particle* particle, ParticleForceGenerator* fg) {
ParticleForceRegistration registration;
registration.particle = particle;
registration.fg = fg;
registrations.push_back(registration);
}
// Clear all registrations from the registry
void clear() {
registrations.clear();
}
// Calls all of the drive turbines to replace the forces of their corresponding particles
void updateForces(actual period) {
Registry::iterator i = registrations.start();
for (; i != registrations.finish(); i++) {
i->fg->updateForce(i->particle, period);
}
}
}; // class ParticleForceRegistry
class ParticlePointGravity : public ParticleForceGenerator {
Vec3 origin;
actual mass;
public:
ParticlePointGravity(const Vec3& origin, const actual mass)
: origin(origin), mass(mass)
{
}
void setOrigin(const Vec3& origin) {
this->origin = origin;
}
void updateForce(Particle* particle, actual period) {
// Examine if particle has a finite mass
if (!particle->hasFiniteMass()) return;
Vec3 drive;
// Calculate the course of the drive
Vec3 course = particle->getPosition() - origin;
actual distance = course.magnitude();
course.normalize();
// Calculate the magnitude of the drive
actual forceMagnitude = -G * particle->getMass() * mass / (distance * distance);
drive = course * forceMagnitude;
// Apply the drive
particle->addForce(drive);
}
}; // class ParticlePointGravity
}; // namespace engine
utilizing namespace engine;
int major() {
// Preliminary situations are on the rightmost level of the orbit (wanting down on the system),
// the place the place is barely in x and the rate is barely in z
Vec3 moon_xi = Vec3(MOONORBIT + EARTHORBIT, 0, 0); // [m]
Vec3 moon_vi = Vec3(0, 0, MOONVELOCITY + EARTHVELOCITY); // [m/s]
Vec3 moon_ai = Vec3(0, 0, 0); // [m/s^2]
Vec3 earth_xi = Vec3(EARTHORBIT, 0, 0); // [m]
Vec3 earth_vi = Vec3(0, 0, EARTHVELOCITY); // [m/s]
Vec3 earth_ai = Vec3(0, 0, 0); // [m/s^2]
Vec3 solGravityOrigin = Vec3(0, 0, 0);
actual solMass = SOLARMASS; // [kg]
Particle *earth = new Particle(earth_xi, earth_vi, earth_ai, 1, 1/EARTHMASS);
Particle *moon = new Particle(moon_xi, moon_vi, moon_ai, 1, 1/MOONMASS);
// Raylib Initialization
//--------------------------------------------------------------------------------------
const int screenWidth = 720;
const int screenHeight = 450;
InitWindow(screenWidth, screenHeight, "Djinn - Photo voltaic Orbit Demo");
// Outline the digital camera to look into our 3d world (place, goal, up vector)
Digicam digital camera = { 0 };
digital camera.place = (Vector3){ 15.0f, 5.0f, 0.0f };
digital camera.goal = (Vector3){ 15.0f, 0.0f, 0.0f };
digital camera.up = (Vector3){ 0.0f, 1.0f, 0.0f };
digital camera.fovy = 60.0f;
digital camera.projection = CAMERA_PERSPECTIVE;
SetCameraMode(digital camera, CAMERA_FIRST_PERSON); // Set a primary particular person digital camera mode
SetCameraMoveControls(KEY_W, KEY_S, KEY_D, KEY_A, KEY_SPACE, KEY_LEFT_SHIFT);
SetTargetFPS(60); // Set our recreation to run at 60 frames-per-second
//--------------------------------------------------------------------------------------
// Time decision
actual dt = 1e4; // [s]
// Outline drive registry
ParticleForceRegistry registry;
// Outline drive turbines to be utilized to the particle
ParticlePointGravity* solGravity = new ParticlePointGravity(solGravityOrigin, solMass); // Args: (Vec3) origin, (actual) mass
ParticlePointGravity* earthGravity = new ParticlePointGravity(earth_xi, EARTHMASS);
// Register drive turbines with the particle
registry.add(earth, solGravity);
registry.add(moon, earthGravity);
whereas(!WindowShouldClose()) {
// Replace the forces per loop
registry.updateForces(dt);
// Calculate the brand new pos, vel, and acc of the moon and earth
earth->combine(dt);
moon->combine(dt);
// Set the earth gravity's new origin to the brand new place of the earth
earthGravity->setOrigin(earth->getPosition());
// Get positions and scale down to suit on display
Vec3 earth_x = earth->getPosition() * 1e-10; // Scale all the way down to tens of meters
Vec3 moon_x = moon->getPosition() * 1e-10; // Scale all the way down to tens of meters
// Convert my Vec3 object to a Raylib Vector3 object
Vector3 *rl_earth_x = reinterpret_cast<Vector3*>(&earth_x); // Raylib vector for drawing
Vector3 *rl_moon_x = reinterpret_cast<Vector3*>(&moon_x);
UpdateCamera(&digital camera);
BeginDrawing();
ClearBackground(BLACK);
BeginMode3D(digital camera);
DrawGrid(1000, 1.0f); // Draw a grid
DrawSphere(*rl_earth_x, 0.25, BLUE); // Draw earth
DrawSphere(*rl_moon_x, 0.1, GRAY); // Draw moon
DrawSphere((Vector3){0, 0, 0}, 1, YELLOW); // Draw Sol
EndMode3D();
DrawRectangle(10, 30, 150, 75, Fade(SKYBLUE, 0.5f));
DrawRectangleLines(10, 30, 150, 75, BLUE);
DrawText("Earth:", 20, 40, 10, WHITE);
DrawText(TextFormat("X: %02.02f", rl_earth_x->x), 20, 55, 10, WHITE);
DrawText(TextFormat("Y: %02.02f", rl_earth_x->y), 20, 70, 10, WHITE);
DrawText(TextFormat("Z: %02.02f", rl_earth_x->z), 20, 85, 10, WHITE);
DrawText("Moon:", 80, 40, 10, WHITE);
DrawText(TextFormat("X: %02.02f", rl_moon_x->x), 80, 55, 10, WHITE);
DrawText(TextFormat("Y: %02.02f", rl_moon_x->y), 80, 70, 10, WHITE);
DrawText(TextFormat("Z: %02.02f", rl_moon_x->z), 80, 85, 10, WHITE);
DrawFPS(10, 10);
EndDrawing();
}
CloseWindow();
return 0;
}
NOTE: The flags to compile this code seems to be like : -framework IOKit -framework Cocoa -framework OpenGL <backtick>pkg-config --libs --cflags raylib<backtick>
