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

// versione con strutture

typedef struct {
  double x;
  double v;
  double a;
} particle;

int main() {
  double t, w2;
  particle *p;
  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");
  p = calloc(N, sizeof(particle));
  
  // initial conditions
  for (j=0; j<N; j++) {
    p[j].x = 0;
    p[j].v = (double)(j+25)/16;
    p[j].a = -w2*sin(p[j].x);
  }
  
  t = 0;
  for (i=0; i<n; i++) {
    fprintf(dat, "%f ", t);
    for (j=0; j<N; j++) {
      p[j].x += p[j].v*dt + 0.5*p[j].a*dt2;
      p[j].v += 0.5*p[j].a*dt;
      p[j].a = -w2*sin(p[j].x);
      p[j].v += 0.5*p[j].a*dt;
      fprintf(dat, "%f %f ", p[j].x, p[j].v);
    }
    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);
}