/*
 * Euler integration of harmonic oscillator
 *
 *  teori=importdata("oscillator_teori.dat");
 *  sim=importdata("oscillator_sim.dat");
 *  plot(sim(:,1),sim(:,2))
 *  hold on
 *  plot(sim(:,1),teori(:,2))
 *
 *  Algoritm:
 *
 *  Given x(0),v(0)
 *
 *  n=0;
 *  while (t<tmax)
 *    {
 *       x(n)=x(n-1)+v(n-1)*dt
 *       v(n)=v(n-1)+F*dt/m
 *
 *       T+=dt;
 *       n++;
 *    }
 */

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

int main()
{
   double x0=1;  // m
   double v0=0;  // m/s
   double k=2;   // k
   double m=1;   // kg

   double w=sqrt(k/m);            // Cyclic frequency
   double T=2.0*M_PI/w;           // Period
   double dt=1e-4*T;

   ofstream fo("oscillator_sim.dat");
   ofstream fot("oscillator_teori.dat");

   double t=0;   // s
   double x=x0,v=v0;  // Initial position/velocity

   while (t<100*T)
     {
         double xn,vn;     // x(n), v(n)
         double F= -k*x;   //  - k3*x*x*x ;  // HER

         xn=x + dt*v;      // x=x(n-1) and v=v(n-1)
         vn=v + dt*F/m;

         double xteori=x0*cos(w*t); // Analytic solution

         cout << t << " " << xn << " " << xteori << endl;
         fo  << t << " " << xn << endl;
         fot << t << " " << xteori  << endl;

         x=xn;
         v=vn;
         t+=dt;
     }
   fo.close();
   fot.close();

   return 0;
}