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

double funcVx(double X,double Y, double A);
double funcVx(double X,double Y, double A)//function for Vx

{
    double a;
    a=A*X/(double)(pow((pow(X,2)+pow(Y,2)),1.5));
    return a;
}

double funcVy(double X,double Y, double A);
double funcVy(double X,double Y, double A)//function for Vy

{
    double a;
    a=A*Y/(double)(pow((pow(X,2)+pow(Y,2)),1.5));
    return a;
}

int main()
{
    FILE *file;
    double G,M,A,deltat,t,tmax,x,y,Vx,Vy,Vxk1,Vxk2,Vxk3,Vxk4,Vyk1,Vyk2,Vyk3,Vyk4,xk1,xk2,xk3,xk4,yk1,yk2,yk3,yk4;
    G=6.673E-11;
    t=0;

    printf("Please enter the mass of the larger body\n");
    scanf("%lf",&M);
    printf("please enter the initial cartesian coordinates of the smaller body\n");
    printf("x=");
    scanf("%lf",&x);
    printf("y=");
    scanf("%lf",&y);
    printf("Please enter the initial velocity in each direction\n");
    printf("Vx=");
    scanf("%lf",&Vx);
    printf("Vy=");
    scanf("%lf",&Vy);
    printf("Please enter the time over which the problem is to be modelled\n");//MORE ELOQUENT!!!!
    scanf("%lf",&tmax);
    printf("please enter the time step\n");
    scanf("%lf",&deltat);
    A=M*G;

    file=fopen("circularmotion.txt","w");

    for (t=0;t<tmax;t=t+deltat)
    {

    Vxk1=funcVx(x,y,A);//Runge-Kutta for Vx
    xk1=Vx;

    Vxk2=funcVx((x+0.5*xk1*deltat),y,A);
    xk2=Vx*Vxk1*0.5*deltat;

    Vxk3=funcVx((x+0.5*xk2*deltat),y,A);
    xk3=Vx*Vxk2*0.5*deltat;

    Vxk4=funcVx((x+xk3*deltat),y,A);
    xk4=Vx*Vxk3*deltat;


    Vyk1=funcVy(x,y,A);//Runge-Kutta for Vx
    yk1=Vy;

    Vyk2=funcVy(x,(y+0.5*yk1*deltat),A);
    yk2=Vy*Vyk1*0.5*deltat;

    Vyk3=funcVy(x,(y+0.5*yk2*deltat),A);
    yk3=Vy*Vyk2*0.5*deltat;

    Vyk4=funcVy(x,(y+yk3*deltat),A);
    yk4=Vy*Vyk3*deltat;


    Vx=Vx+deltat*(Vxk1/(double)6)+(Vxk2/(double)3)+(Vxk3/(double)3)+(Vxk4/(double)6);
    Vy=Vy+deltat*(Vyk1/(double)6)+(Vyk2/(double)3)+(Vyk3/(double)3)+(Vyk4/(double)6);
    x=x+deltat*(xk1/(double)6)+(xk2/(double)3)+(xk3/(double)3)+(xk4/(double)6);
    y=y+deltat*(yk1/(double)6)+(yk2/(double)3)+(yk3/(double)3)+(yk4/(double)6);

    fprintf(file,"%g\t%g\t%g\t%g\n",x,y,Vx,Vy);
    printf("x=%g, y=%g, Vx=%g, Vy=%g\n",x,y,Vx,Vy);

    }

fclose(file);

    return 0;
}