// Appendix D.1 (p121) of Klein and Godunov Introductory Computational Physics
// program to calculate the mass on a spring problem.
//m=k=1; therefore a=-x
//with the simple Euler method
//ak 9/7/00
//all in double precision
//small adjustments isv Jan. 2017
#include <iostream>
#include <math.h>
#include <stdio.h>
using namespace std;

void deriv(double y[], double f[]);

void deriv(double y[], double f[])
{
   f[0]=y[1];  // x' = v,   x is y0
   f[1]=-y[0]; // v' = -x   v is y1
   return;
}


int main (int argc, char **argv)
{
   double y[2]; // y[0] is 'x', y[1] is 'v'
   double vinit= 5.;// initial velocity
   double xinit= 0.;// initial position
   double energy; // the energy from mv**2/2+kx**2/2
   double time=0.; //the time
   double ymid[2], f[2];

   if(argc != 3) {
      printf("Need two arguments: n and step\n");
      exit(0);
   }

   int n=atoi(argv[1]);   // number of steps
   double step = atof(argv[2]);  //step size


   int i = 0;//the loop counter
// initial conditions
   y[0]=xinit;
   y[1]=vinit;

   energy=0.5*xinit*xinit + 0.5*vinit*vinit;
   cout << time << " " << y[0] << " " << y[1] << " " << energy << endl;

   while(i < n)  //make n steps
      {
      deriv(y, f);
      ymid[0] = y[0] + 0.5*step*f[0];
      ymid[1] = y[1] + 0.5*step*f[1];

      deriv(ymid, f);
      y[0]=y[0]+step*f[0]; //calculate the forward position 
      y[1]=y[1]+step*f[1]; //calculate the forward velocity

      energy=pow(y[0],2)/2.+pow(y[1],2)/2.; //energy conservation
      time=time+step;
      i=i+1;
      cout << time << " " << y[0] << " " << y[1] << " " << energy << endl;
      }
}