#include <stdio.h>
#include <math.h>
typedef struct
{
double x, y;
double vx, vy;
double r, phi;
double rmax, rmin;
} State;
static void update(State *s, double GM, double dt)
{
const double r = sqrt(s->x * s->x + s->y * s->y),
F = -GM / (r * r);
s->vx += F * s->x / r * dt;
s->vy += F * s->y / r * dt;
s->x += s->vx * dt;
s->y += s->vy * dt;
s->r = r;
s->phi = atan2(s->y, s->x);
if (s->rmax == 0 || s->r > s->rmax)
s->rmax = s->r;
if (s->rmin == 0 || s->r < s->rmin)
s->rmin = s->r;
}
int main()
{
const double GM = 1.0, dt = 0.0001;
State planets[] = {
{.x = 10.0, .vx = 0.0, .vy = 0.2},
{.x = 5.0, .vx = 0.1, .vy = 0.1}};
for (int i = 0; i < sizeof(planets) / sizeof(planets[0]); i++)
{
double t = 0, phiprev = 0.0;
while (!(planets[i].phi > 0 && phiprev < 0))
{
phiprev = planets[i].phi;
update(&planets[i], GM, dt);
t += dt;
}
const double a = (planets[i].rmax + planets[i].rmin) / 2.0;
printf("Planet %d: a^3 / T^2 = %f\n", i, (a * a * a) / (t * t));
}
}
