#include <iostream>
#include <cmath>
using namespace std;

struct astral_body{
    double position[3];
    double velocity[3];
    double mass;
};

//function declaration
void compute_leapfrog(astral_body &planet, astral_body &star_1, astral_body &star_2, double dt);

int main()
{
    double dt=1;
    int N=10;

    astral_body planet;
    astral_body star_1;
    astral_body star_2;

    //initial parameters
    star_1.position[0]=1;
    star_1.position[1]=0;
    star_1.position[2]=0;
    star_2.position[0]=-1;
    star_2.position[1]=0;
    star_2.position[2]=0;
    planet.position[0]=0;
    planet.position[1]=0;
    planet.position[2]=0;

    planet.mass=1;
    star_1.mass=100;
    star_2.mass=100;

    cout<<"enter the x component of the planet's velocity"<<endl;
    cin>>planet.velocity[0];
    cout<<"enter the y component of the planet's velocity"<<endl;
    cin>>planet.velocity[1];
    cout<<"enter the z component of the planet's velocity"<<endl;
    cin>>planet.velocity[2];


    for(int count=1; count<N; count++)
    {
        cout<<planet.position[0]<<" "<<planet.position[1]<<"    "<<planet.position[2]<<endl;
        compute_leapfrog(planet, star_1, star_2, dt);
    }

}

void compute_leapfrog(astral_body &planet, astral_body &star_1, astral_body &star_2, double dt)
{
    double G=100;
    double f0, f1, f2;

    planet.position[0]+=0.5*dt*planet.velocity[0];
    planet.position[1]+=0.5*dt*planet.velocity[1];
    planet.position[2]+=0.5*dt*planet.velocity[2];

    f0=(G*planet.mass*star_1.mass)*((1.0/pow((star_1.position[0]-planet.position[0]),2.0)))+((G*planet.mass*star_2.mass)/(1.0/pow((star_2.position[0]-planet.position[0]),2.0)));
    f1=(G*planet.mass*star_1.mass)*((1.0/pow((star_1.position[1]-planet.position[1]),2.0)))+((G*planet.mass*star_2.mass)/(1.0/pow((star_2.position[1]-planet.position[1]),2.0)));
    f2=(G*planet.mass*star_1.mass)*((1.0/pow((star_1.position[2]-planet.position[2]),2.0)))+((G*planet.mass*star_2.mass)/(1.0/pow((star_2.position[2]-planet.position[2]),2.0)));

    cout<<f0<<f1<<f2<<endl;

    planet.velocity[0]+=dt*(f0/planet.mass);
    planet.velocity[1]+=dt*(f1/planet.mass);
    planet.velocity[2]+=dt*(f2/planet.mass);

    planet.position[0]+=0.5*dt*planet.velocity[0];
    planet.position[1]+=0.5*dt*planet.velocity[1];
    planet.position[2]+=0.5*dt*planet.velocity[2];

    //any steps about incrementing t have been left out as I am aiming to plot the planet's path, not anything with time.
    }