// Equation of motion solver for a particle in 1D using leapfrog method
#include <stdio.h>
#include <math.h>

double f(double k, double x);
double analyticalsol(double x);
void leapfrog(double dt, double t0, double x0, double v0, double m, double k, double tmax);
void plotfunction(double x0, double x1, int N);

// Driver program
int main() 
{ 
    // Initial values
    leapfrog(0.5,0.0,1.0,0.0, 1.0,1.0, 50); 
    plotfunction(0.00000001, 50,1000);
    return 0; 
} 


// Prints root of func(x) with error of EPSILON 
void leapfrog(double dt, double t0, double x0, double v0, double m, double k, double tmax) 
{ 
   double x,v,ene,t;
   FILE *output,*output2;

   output = fopen("emleapfrog.dat", "w");//opening file. 
   output2 = fopen("emleapfrog-energy.dat", "w");//opening file. 

    t=t0;
    x=x0;
    v=v0-dt/2*f(k,x)/m;    
    fprintf(output,"%.15f\t%.15f\n",t0,x0);

    while (t <= tmax) 
    {  
       v=v0+0.5*dt*f(k,x0)/m;
       x=x0+dt*v;
       v=v+0.5*dt*f(k,x)/m;
       v0=v;
       x0=x;
       t=t+dt;
      fprintf(output,"%.15f\t%.15f\n",t,x);
      ene=0.5*k*x*x+0.5*m*v*v;
      fprintf(output2,"%.15f\t%.15f\n",t,ene);

    } 
     fclose(output);//closing file.
     fclose(output2);//closing file.

} 


void plotfunction(double x0, double x1,int N)
{ int i;
  double x,dx;
  FILE *output;

      dx=(x1-x0)/N;


output = fopen("function.dat", "w");//opening file. 
 for(i=0;i<N;i++){
    x=x0+dx*i;
fprintf(output,"%.15f\t%.15f\n",x,analyticalsol(x));
 }
 
 fclose(output);//closing file.

}


// Define function

double f(double k, double x) 
{ 
    double fun;

fun=-k*x;
    return fun;
} 

double analyticalsol(double x) 
{ 
    double y;

y=cos(x);
    return y;
}