Numerical integration - gravity

## Numerical integration - gravity

```// A program to illustrate numerical integration for solving gravitational
// interactions between two bodies.

#include <stdio.h>
#include <math.h>

// Physical constants.

const double bigG = 6.6726E-11; // Nm^2/kg^2

// Useful typedefs, i.e., our own data types for the problem at hand.

// A 3-D vector.

typedef struct {
double x, y, z;
} vector;

// All the information that we need to know about a test particle.

typedef struct {
double mass;      // The mass of the particle, in kg.
double epoch;     // The time in seconds at which the following
//     two items are specified.
vector pos;       // The particle's position.
vector vel;       // The particle's velocity.
vector force;     // The force acting on the particle.
double potential; // The gravitational potential that the particle is in.
double ke;        // The particle's kinetic energy.
double pe;        // The particle's potential energy.
vector I;         // The particle's angular momentum.
} particle;

// Now we define some functions that we will need.

void gravity(particle *p1, particle *p2) {

// Calculates the force and potential acting on particle p1 due to
// particle p2.

double forceMagnitude;
double distance;
vector d;

// Calculate the position offset between the two particles.

d.x = p1->pos.x - p2->pos.x;
d.y = p1->pos.y - p2->pos.y;
d.z = p1->pos.z - p2->pos.z;

// The magnitude of the position offset vector is the distance between
// the particles.

distance = sqrt(d.x*d.x + d.y*d.y + d.z*d.z);

// Now calculate the magnitude of the force acting on p1 due to p2.

forceMagnitude = bigG * p1->mass * p2->mass / (distance * distance);

// And orient it correctly in 3 dimensions by multiplying by a unit
// vector from p2 to p1.

p1->force.x = - forceMagnitude * d.x / distance;
p1->force.y = - forceMagnitude * d.y / distance;
p1->force.z = - forceMagnitude * d.z / distance;

// Finally, calculate the potential at p1 due to p2.

p1->potential = - bigG * p2->mass / distance;
}

void update(particle *p, double deltaT) {

// Updates the properties of particle p after a time interval
// deltaT [seconds].

vector acc;

// Calculate the acceleration due to the force acting on the particle.

acc.x = p->force.x / p->mass;
acc.y = p->force.y / p->mass;
acc.z = p->force.z / p->mass;

// Update the position vector, assuming a constant acceleration
// over the time interval.

p->pos.x += (p->vel.x + 0.5 * acc.x * deltaT) * deltaT;
p->pos.y += (p->vel.y + 0.5 * acc.y * deltaT) * deltaT;
p->pos.z += (p->vel.z + 0.5 * acc.z * deltaT) * deltaT;

// Update the velocity vector.

p->vel.x += acc.x * deltaT;
p->vel.y += acc.y * deltaT;
p->vel.z += acc.z * deltaT;

// Update the epoch.

p->epoch += deltaT;

// Update the kinetic energy.

p->ke = 0.5 * p->mass * (p->vel.x * p->vel.x +
p->vel.y * p->vel.y +
p->vel.z * p->vel.z);

// Update the potential energy.

p->pe = p->mass * p->potential;

// Update the angular momentum.

p->I.x = p->mass * (p->pos.y * p->vel.z - p->pos.z * p->vel.y);
p->I.y = p->mass * (p->pos.z * p->vel.x - p->pos.x * p->vel.z);
p->I.z = p->mass * (p->pos.x * p->vel.y - p->pos.y * p->vel.x);
}

void displayParticle(particle *p) {

// Displays information about particle p.

printf("m %10.2e; x (%10.2e,%10.2e,%10.2e); v (%10.2e,%10.2e,%10.2e)",
p->mass,
p->pos.x, p->pos.y, p->pos.z,
p->vel.x, p->vel.y, p->vel.z);
}

void displayEnergy(particle *p) {

printf(" E (%10.2e =%10.2e +%10.2e)", p->ke + p->pe, p->ke, p->pe);
}

void displayI(particle *p) {

printf(" I (%10.2e,%10.2e, %10.2e)", p->I.x, p->I.y, p->I.z);
}

void displayCM(particle *p1, particle *p2) {
vector cm;

cm.x = (p1->mass * p1->pos.x  + p2->mass * p2->pos.x) /
(p1->mass + p2->mass);
cm.y = (p1->mass * p1->pos.y  + p2->mass * p2->pos.y) /
(p1->mass + p2->mass);
cm.z = (p1->mass * p1->pos.z  + p2->mass * p2->pos.z) /
(p1->mass + p2->mass);
printf(" cm (%10.2e,%10.2e, %10.2e)", cm.x, cm.y, cm.z);
}

// Here are two representative particles:

particle sun =   {1.989E30,              // Mass in kg
0.0,
{0.0, 0.0, 0.0},       // The initial position
{0.0, 0.0, 0.0}};      // The initial velocity

particle earth = {5.976E24,              // Mass in kg
0.0,
{149.6E9, 0.0, 0.0},   // The initial position
{0.0, 27.786E3, 0.0}}; // The initial velocity

const double deltaT = 3600.0; // The time interval in seconds

int main (void) {

double time;

// Choose the position of the Sun to put the centre of mass at (0,0,0).

sun.pos.x = - earth.pos.x * earth.mass/sun.mass;
sun.pos.y = - earth.pos.y * earth.mass/sun.mass;
sun.pos.z = - earth.pos.z * earth.mass/sun.mass;

// Choose the velocity of the Sun to keep the centre of mass
// at a constant position.

sun.vel.x = 0.0;
sun.vel.y = earth.vel.y * sun.pos.x / earth.pos.x;
sun.vel.z = 0.0;

// Now follow the motion for one year...

for (time = 0.0; time < 365.2422 * 24 * 3600; time += deltaT) {
gravity(&sun, &earth);
gravity(&earth, &sun);

update(&sun,   deltaT);
update(&earth, deltaT);

//		displayParticle(&sun);
//      printf("\n");
#if 0
displayParticle(&sun);
displayEnergy(&sun);
displayI(&sun);
//       	displayCM(&earth, &sun);
printf("\n");
#endif
displayParticle(&earth);
displayEnergy(&earth);
displayI(&earth);

printf("\n");

}
return 0;
}

```