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

const float G=-6.673E-11;//stores the gravitational constant

double funcx(double M, double x, double y);//the functions for the velocities
double funcy(double M, double x, double y);

int main()
{
    FILE *path;
    double x,y,t,T,dt,vx,vy,M,k1x,k2x,k3x,k4x,k1y,k2y,k3y,k4y,k1vx,k2vx,k3vx,k4vx,k1vy,k2vy,k3vy,k4vy;

    printf("Calculating a basic orbit of a moving body\n");

    printf("\nEnter the mass of the large object: ");
    scanf("%lf",&M);
    printf("\nEnter the initial coordinates of the moving object\n");
    printf("x = ");
    scanf("%lf",&x);
    printf("y = ");
    scanf("%lf",&y);
    printf("\nEnter the initial velocities of the moving object\n");
    printf("x velocity = ");
    scanf("%lf",&vx);
    printf("y velocity = ");
    scanf("%lf",&vy);
    printf("\nEnter the time over which to map the orbit: ");
    scanf("%lf",&T);
    printf("Enter the size of the time interval to calculate over: ");
    scanf("%lf",&dt);

    path=fopen("path.txt","w");

    for(t=0;t<T;t+=dt)
    {
        k1vx=funcx(M,x,y);
        k1x=vx;
        k2vx=funcx(M,(x+(k1x*dt/2)),y);
        k2x=vx*k1vx*(dt/2);
        k3vx=funcx(M,(x+(k2x*dt/2)),y);
        k3x=vx*k2vx*(dt/2);
        k4vx=funcx(M,(x+(k3x*dt)),y);
        k4x=vx*k3x*dt;

        k1vy=funcy(M,x,y);
        k1y=vy;
        k2vy=funcy(M,x,(y+(k1y*dt/2)));
        k2y=vy*k1vy*(dt/2);
        k3vy=funcy(M,x,(y+(k2y*dt/2)));
        k3y=vy*k2vy*(dt/2);
        k4vy=funcy(M,x,(y+(k3y*dt)));
        k4y=vy*k3vy*dt;

        vx+=(dt/6)*(k1vx+(2*k2vx)+(2*k3vx)+k4vx);
        vy+=(dt/6)*(k1vy+(2*k2vy)+(2*k3vy)+k4vy);

        x+=dt*vx;
        y+=dt*vy;

        fprintf(path,"%f%c%f%c%f%c%f\n",x,9,y,9,vx,9,vy);
    }

    return 0;
}
double funcx(double M, double x, double y)
{
    return (M*G*x)/pow((x*x)+(y*y),3/(double)2);
}
double funcy(double M, double x, double y)
{
    return (M*G*y)/pow((x*x)+(y*y),3/(double)2);
}