/* Questo esempio non e' incluso nel testo, ma e' utile per capire in
dettaglio come funziona Runge-Kutta del quarto ordine nel caso 
dell'integrazione delle equazioni di Newton. Il caso specifico che e'
considerato qui e' un oscillatore anarmonico, con un contributo quartico, 
in presenza di un termine di attrito */

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

struct xv{
  double x;
  double v;
};

struct xv  rk2(struct xv xV,  double gamma, double dt, long int steps);
struct xv  rk4(struct xv xV,  double gamma, double dt, long int steps);
double forzaItinere(struct xv xV,  double gamma);

/***********************************************************************/
int main(){
  struct xv xV;
  double x0, v0;
  double gamma, dt, tMax = 100.;
  long int steps;
  printf("Please input gamma\n"); fflush(stdout);
  scanf("%lg",&gamma);
  printf("Please input dt\n"); fflush(stdout);
  scanf("%lg",&dt);
  printf("Please input x0\n"); fflush(stdout);
  scanf("%lg",&x0);
  printf("Please input v0\n"); fflush(stdout);
  scanf("%lg",&v0);
  steps = tMax / dt;
  printf("# gamma %lg dt %lg tMax %lg steps %ld x0 %lg v0 %lg\n",
	 gamma, dt, tMax, steps, x0, v0);
  xV.x = x0;
  xV.v = v0;
  xV = rk2(xV, gamma, dt, steps);
  xV.x = x0;
  xV.v = v0;
  xV = rk4(xV, gamma, dt, steps);
  return 0;
}

/***********************************************************************/
struct xv  rk2(struct xv xVOld,  double gamma, double dt, long int steps){
  struct xv xVScratch, xVNew = {0.0, 0.0};
  double xP, vP;
  int t;
  for(t=0; t< steps; t++){
    xP = xVOld.v * dt;
    vP = forzaItinere(xVOld, gamma) * dt;
    xVScratch.x = xVOld.x + 0.5 * xP;
    xVScratch.v = xVOld.v + 0.5 * vP;
    xVNew.x = xVOld.x + xVScratch.v * dt;
    xVNew.v = xVOld.v + forzaItinere(xVScratch, gamma) * dt;
    printf(" t %lg x %lg v %lg gamma %lg dt %lg RK2\n",
	   dt * t, xVNew.x, xVNew.v, gamma, dt);
    xVOld.x = xVNew.x;
    xVOld.v = xVNew.v;
  }
  return xVNew;
}

/***********************************************************************/
struct xv  rk4(struct xv xVOld,  double gamma, double dt, long int steps){
  struct xv xVSum, x1, x2, x3, x4, xVNew = {0.0, 0.0};
  int t;
  for(t=0; t< steps; t++){
    x1.x = xVOld.v * dt;
    x1.v = forzaItinere(xVOld, gamma) * dt;
    xVSum.x = xVOld.x + 0.5 * x1.x;
    xVSum.v = xVOld.v + 0.5 * x1.v;
    x2.x = xVSum.v * dt;
    x2.v = forzaItinere(xVSum, gamma) * dt;
    xVSum.x = xVOld.x + 0.5 * x2.x;
    xVSum.v = xVOld.v + 0.5 * x2.v;
    x3.x = xVSum.v * dt;
    x3.v = forzaItinere(xVSum, gamma) * dt;
    xVSum.x = xVOld.x + x3.x;
    xVSum.v = xVOld.v + x3.v;
    x4.x = xVSum.v * dt;
    x4.v = forzaItinere(xVSum, gamma) * dt;
    xVNew.v = xVOld.v + 1./6. * (x1.v + 2. * x2.v + 2. * x3.v + x4.v);
    xVNew.x = xVOld.x + 1./6. * (x1.x + 2. * x2.x + 2. * x3.x + x4.x);
    printf(" t %lg x %lg v %lg gamma %lg dt %lg RK4\n",
	   dt * t, xVNew.x, xVNew.v, gamma, dt);
    xVOld.x = xVNew.x;
    xVOld.v = xVNew.v;
  }
  return xVNew;
}

/***********************************************************************/
double forzaItinere(struct xv xV,  double gamma){
  double f, x1, x2, x3;
  x1 = xV.x;
  x2 = x1 * x1;
  x3 = x1 * x2;
  f = - x1 + 3. * x2 - 2. * x3 - gamma * xV.v;
  return f;
}