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

int main() {
  double t, x, v, a, w2, A, e;
  double dt, dt2;
  int i, n;
  double gamma;// dumping in units of m
  double F, Omega;
  double K;
  int trans;
  double maxx;
  FILE *dat, *gp;
  
  A = 1;
  w2 = 1;
  dt = 0.1;
  dt2=dt*dt;
  n = 10000;
  gamma = 0.1;
  trans = 9000;
  F = 1;
  dat = fopen("ScanForcedDumpedHarmonicOscillator.dat", "w");

  for (Omega=0.5; Omega < 1.5; Omega += 0.01) {
    // initial conditions
    x = A;
    v = 0;
    a = -w2*x;
    t = 0;
    maxx = 0;
    for (i=0; i< n; i++) {
      x += v*dt + 0.5*a*dt2;
      v += 0.5*a*dt;
      K = F * sin(Omega * t);
      a = -w2*x - gamma * v + K;
      v += 0.5*a*dt;
      e = v*v+w2*x*x;    
      t += dt;
      if (i>trans) {
        if (x > maxx) maxx = x;
      }
    }
    fprintf(dat, "%f %f\n", Omega, maxx );
 
  }
  fclose(dat);
  gp = popen("gnuplot", "w");
  fprintf(gp, "unset key\n");
  fprintf(gp, "plot 'ScanForcedDumpedHarmonicOscillator.dat' u 1:2 w l \n");
  fflush(gp);
  getchar();
  fclose(gp);
  return(0);
}