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

#define Max 1000000

double x[Max],v[Max];
double xerr[Max],verr[Max];
double xth[Max],vth[Max];
double conv[Max][2];
double dt=0.01,g=1,k=1,m=1,x0=0,v0=0;
double T=1;
int ItNo=20,NConv=10;
int pointno=0;

double F(double x){
  return -k*x-g;
}

void Verlet(double *x,double *v,double x1,double v1,double dt){
  double Fx=F(x1);
  *x = x1+(v1+0.5*Fx*dt)*dt;
  *v = v1+0.5*(Fx+F(*x))*dt;
}

void EulerIt(double *x,double *v,double x1,double v1,double dt){
  *x = x1+v1*dt;
  *v = v1+F(x1)/m*dt;
}

void (*It)(double *,double *,double ,double ,double);

void setEuler(){
  It=&EulerIt;
}

void setVerlet(){
  It=&Verlet;
}

void Newton(){
  x[0]=x0;
  v[0]=v0;
  for (int i=1; i<ItNo;i++){
    It(&x[i],&v[i],x[i-1],v[i-1],dt);
  }
}

double Xth(double t){
  double km=sqrt(k/m);
  return sqrt(1/km)*v0*sin(km*t)+x0*cos(km*t);
}

double Vth(double t){
  double km=sqrt(k/m);
  return v0*cos(km*t)-x0*km*sin(km*t);
}

void Newton_th(){
  for (int i=0; i<ItNo;i++){
    xth[i]=Xth(i*dt);
    vth[i]=Vth(i*dt);
  }
}

void Newton_err(){
  for (int i=0; i<ItNo;i++){
    xerr[i]=x[i]-xth[i];
    verr[i]=v[i]-vth[i];
  }
}

void Do_conv(){
  for (int i=5;i<NConv;i++){
    dt=T/i;
    double x=x0, v=v0;
    for (int j=0;j<i;j++){
      It(&x,&v,x,v,dt);
    }
    conv[i][0]=dt;
    conv[i][1]=x-Xth(i*dt);
  }
}

int main(){
  int done=0;
  setVerlet();
  DefineGraphN_R("x",&x[0],&ItNo,NULL);
  DefineGraphN_R("v",&v[0],&ItNo,NULL);
  DefineGraphN_R("xth",&xth[0],&ItNo,NULL);
  DefineGraphN_R("vth",&vth[0],&ItNo,NULL);
  DefineGraphN_R("xerr",&xerr[0],&ItNo,NULL);
  DefineGraphN_R("verr",&verr[0],&ItNo,NULL);
  DefineGraphN_RxR("Conv",&conv[0][0],&NConv,NULL);
  StartMenu("Start",1);
  DefineGraph(curve2d_,"Graphics");
  DefineInt("ItNo",&ItNo);
  DefineDouble("dt",&dt);
  DefineDouble("x0",&x0);
  DefineDouble("v0",&v0);
  DefineDouble("g",&g);
  DefineDouble("k",&k);
  DefineDouble("m",&m);
  DefineFunction("Newton",&Newton);
  DefineFunction("Newton Th",&Newton_th);
  DefineFunction("Newton Err",&Newton_err);
  DefineInt("NConv",&NConv);
  DefineDouble("T",&T);
  DefineFunction("Conv",&Do_conv);
  DefineFunction("Euler",&setEuler);
  DefineFunction("Verlet",&setVerlet);
  DefineBool("done",&done);
  EndMenu();
  while (!done){
    Events(1);
    DrawGraphs();
    sleep(1);
  }
}