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

#define Max 1000000
#define N 10
double x[N],v[N]; // state of the system
//double xg[N][Max][2],vg[N][Max][2]; //graphical representation of time-evolution
double dt=0.01,g=1,k=1,m=1,x0[N],v0[N];
int ItNo=1;
int r=10;
double fact=2;
void FF(double F[N],double x[N]){
  F[0]= -g -k*(x[0])+k*(x[0+1]-x[0]);
  for (int i=1; i<N-1; i++){
    F[i]= -g -k*(x[i]-x[i-1])+k*(x[i+1]-x[i]);
  }
  F[N-1]= -g -k*(x[N-1]-x[N-1-1]);
}

void Verlet(double x[N],double v[N],double dt){
  double F[N];
  FF(F,x);
  for (int i=0; i<N; i++){
    v[i]+= 0.5*F[i]/m*dt;
    x[i]+= v[i]*dt;
  }
  FF(F,x);
  for (int i=0;i<N; i++){
    v[i]+= 0.5*F[i]/m*dt;
  }
}

void Init(){
  for (int i=0;i<N;i++){
    x[i]=x0[i];
    v[i]=v0[i];
  }
}
void Newton(){
  for (int i=0; i<ItNo;i++){
    Verlet(x,v,dt);
  }
}

void Draw(int xdim,int ydim){
  double hexp=fact*(N*(N-1)/2.*m*g/k);
  for (int i=0; i<N; i++){
    mycircle(1,xdim/2,ydim/2+x[i]/hexp*ydim,r);
  }
}

int main(){
  int sstep=0,cont=0,done=0;
  Init();
  AddFreedraw("Oscillators",&Draw);
  //  DefineGraphN_R("x",&x[0],&ItNo,NULL);
  //  DefineGraphN_R("v",&v[0],&ItNo,NULL);
  StartMenu("Start",1);
  //  DefineGraph(curve2d_,"Graphics");
  DefineGraph(freedraw_,"Oscillators");
  DefineInt("ItNo",&ItNo);
  DefineDouble("dt",&dt);
  StartMenu("Init",0);
  for (int i=0; i<N; i++)
    DefineDouble("x0",&x0[i]);
  for (int i=0; i<N; i++)
    DefineDouble("v0",&v0[i]);
  DefineFunction("Init",&Init);
  EndMenu();
  DefineDouble("g",&g);
  DefineDouble("k",&k);
  DefineDouble("m",&m);
  DefineDouble("scale",&fact);
  DefineInt("radius",&r);
  DefineBool("sstep",&sstep);
  DefineBool("cont",&cont);
  DefineBool("done",&done);
  EndMenu();
  while (!done){
    Events(1);
    DrawGraphs();
    if (cont||sstep){
      sstep=0;
      Newton();
    }
    else sleep(1);
  }
}