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


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


typedef struct {
    double mass;
    Point position, velocity, acceleration;
} Particle;


const double G = 6.6742367e-11;

void integrator_leapfrog_part1(int n_particles, Particle* particles, double half_time_step){
    for (int i = 0; i < n_particles; i++) {
        particles[i].position.x += half_time_step * particles[i].velocity.x;
        particles[i].position.y += half_time_step * particles[i].velocity.y;
        particles[i].position.z += half_time_step * particles[i].velocity.z;
    }
}


void integrator_leapfrog_part2(int n_particles, Particle* particles, double time_step, double half_time_step){
    for (int i = 0; i < n_particles; i++) {
        particles[i].velocity.x += time_step * particles[i].acceleration.x;
        particles[i].velocity.y += time_step * particles[i].acceleration.y;
        particles[i].velocity.z += time_step * particles[i].acceleration.z;

        particles[i].position.x += half_time_step * particles[i].velocity.x;
        particles[i].position.y += half_time_step * particles[i].velocity.y;
        particles[i].position.z += half_time_step * particles[i].velocity.z;
    }
}


void gravity_calculate_acceleration(int n_particles, Particle* particles) {
    for (int i = 0; i < n_particles; i++) {
        particles[i].acceleration.x = 0.0;
        particles[i].acceleration.y = 0.0;
        particles[i].acceleration.z = 0.0;
        for (int j = 0; j < n_particles; j++) {
            if (j == i) {
                continue;
            }
            double dx = particles[i].position.x - particles[j].position.x;
            double dy = particles[i].position.y - particles[j].position.y;
            double dz = particles[i].position.z - particles[j].position.z;
            double r = sqrt(dx*dx + dy*dy + dz*dz);
            double prefact = -G/(r*r*r) * particles[j].mass;
            particles[i].acceleration.x += prefact * dx;
            particles[i].acceleration.y += prefact * dy;
            particles[i].acceleration.z += prefact * dz;
        }
    }
}


int main() {

    double time = 0.0;
    double time_step = 0.08;
    double half_time_step = time_step/2;
    double time_limit = 365.25 * 1e4;

    Particle star = {
            .mass = 0.08,
            .position = {.x = 0.0, .y = 0.0, .z = 0.0},
            .velocity = {.x =  0.0, .y = 0.0, .z = 0.0},
            .acceleration = {.x = 0.0, .y = 0.0, .z = 0.0}
    };
    Particle planet = {
            .mass = 3.0e-6,
            .position = {.x = 0.0162, .y = 6.57192058353e-15, .z = 5.74968548652e-16},
            .velocity = {.x = -1.48427302304e-14, .y = 0.0399408809121, .z = 0.00349437429104},
            .acceleration = {.x = 0.0, .y = 0.0, .z = 0.0},
    };

    int num_particles = 2;
    Particle* particles = malloc(2*sizeof(Particle));
    particles[0] = star;
    particles[1] = planet;

    while (time < time_limit) {
        integrator_leapfrog_part1(num_particles, particles, half_time_step);
        time += half_time_step;
        gravity_calculate_acceleration(num_particles, particles);
        integrator_leapfrog_part2(num_particles, particles, time_step, half_time_step);
        time += half_time_step;
    }

    printf("%f", particles[1].position.x);
    return 0;
}