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

// versione con vettori

int main() {
  double t, *x, *v, *a, w2;
  double dt, dt2;
  int i, j, n;
  int N; // number of initial conditions
  FILE *dat, *gp;
  
  w2 = 1;
  dt = 0.1;
  dt2=dt*dt;
  n = 100;
  N=8;
  dat = fopen("Pendulum.dat", "w");
  x = calloc(N, sizeof(double));
  v = calloc(N, sizeof(double));
  a = calloc(N, sizeof(double));
  
  // initial conditions
  for (j=0; j<N; j++) {
    x[j] = 0;
    v[j] = (double)(j+25)/16;
    a[j] = -w2*x[j];
  }
  
  t = 0;
  for (i=0; i<n; i++) {
    fprintf(dat, "%f ", t);
    for (j=0; j<N; j++) {
      x[j] += v[j]*dt + 0.5*a[j]*dt2;
      v[j] += 0.5*a[j]*dt;
      a[j] = -sin(w2*x[j]);
      v[j] += 0.5*a[j]*dt;
      fprintf(dat, "%f %f ", x[j], v[j]);
    }
    fprintf(dat, "\n");
    t += dt;
  }
  fclose(dat);
  gp = popen("gnuplot", "w");
  fprintf(gp, "plot 'Pendulum.dat' u 2:3 w l");
  for (j=1; j<N; j++) {
    fprintf(gp, ", '' u %d:%d w l", 2*j+2, 2*j+3);
  }
  fprintf(gp, "\n");
  fflush(gp);
  getchar();
  fclose(gp);
  return(0);
}