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

int main() {
  double t, x, v, a, w2, A, e;
  double dt, dt2;
  double Omega, F; // forces, in units of m
  double FF;
  double gamma; // dumping
  int i, n;
  FILE *dat, *gp;
  
  A = .4;
  w2 = 1;
  dt = 0.02;
  dt2=dt*dt;
  n = 10000;
  F = 1;
  gamma = 0.1;
  Omega = 1.01;
  
 
  // initial conditions
  x = A;
  v = 0;
  a = -w2*x;
  t = 0;
  dat = fopen("ForcedHarmonicOscillator.dat", "w");
  for (i=0; i< n; i++) {
    x += v*dt + 0.5*a*dt2;
    v += 0.5*a*dt;
    FF = F * sin(Omega * t); // forcing
    a = -w2*x - gamma * v + FF;
    v += 0.5*a*dt;
    e = v*v+w2*x*x;    
    t += dt;
    if (i > n/2) {
      fprintf(dat, "%f %f %f %f\n", t, x, v, FF);
    }
  }
  fclose(dat);
  gp = popen("gnuplot", "w");
  fprintf(gp, "set multiplot layout 1,2\n");
  fprintf(gp, "plot 'ForcedHarmonicOscillator.dat'  u 1:2 t 'x' w l, ''  u 1:4 t 'F' w l\n");
  fprintf(gp, "unset key; set size square; plot 'ForcedHarmonicOscillator.dat'  u 2:4 w l\n");
  fflush(gp);
  getchar();
  fclose(gp);
  return(0);
}