#include "sopnamsp.h"
#include "machdefs.h"
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <unistd.h>

#include "ntuple.h"
#include "geneutils.h"
/*
...Check autour de zero
> cmvsinxsx -n 4 1e-6,-0.1,0.1
> cmvsinxsx -n 5 1e-6,-0.1,0.1
...Check autour de Pi
> cmvsinxsx -n 4 1e-6,3.041592659,3.241592659
> cmvsinxsx -n 5 1e-6,3.041592659,3.241592659
...Check autour de -Pi
> cmvsinxsx -n 4 1e-6,-3.241592659,-3.041592659
> cmvsinxsx -n 5 1e-6,-3.241592659,-3.041592659

...Check longue echelle
> cmvsinxsx -n 10 0.0001,-6.5,6.5
> cmvsinxsx -n 11 0.0001,-6.5,6.5

...Check sans plot
> cmvsinxsx -n 100 0.00001,-6.5,6.5
> cmvsinxsx -n 1000 0.00001,-6.5,6.5
 */


int main(int narg,char *arg[])
{
 double dt = 0.000001;
 //double tmin = -0.1, tmax = 0.1;
 double tmin = M_PI-0.1, tmax = M_PI+0.1;
 unsigned long N = 4;
 bool only_test = false;
 double pr_test = 1.e99;

 char c;
 while((c = getopt(narg,arg,"htp:n:")) != -1) {
  switch (c) {
  case 'n' :
    sscanf(optarg,"%lu",&N);
    break;
  case 't' :
    only_test = true;
    break;
  case 'p' :
    sscanf(optarg,"%lf",&pr_test);
    break;
  case 'h' :
  default :
    cout<<"cmvsinxsx [-t] [-p eps_print] [-n N] [dt,tmin,tmax (rad)]"<<endl;
    return -1;
  }
 }
 if(optind<narg) sscanf(arg[optind],"%lf,%lf,%lf",&dt,&tmin,&tmax);

 cout<<"N="<<N<<" dt="<<dt<<" tmin="<<tmin<<" tmax="<<tmax<<" rad"<<endl;
 cout<<"only_test="<<only_test<<" pr_test="<<pr_test<<endl;
 if(tmax<=tmin) return -1;

 const int nvar = 11;
 char *vname[nvar] = {"t","s","sn","sx","sx2","snx","snx2","asx","asx2","asnx","asnx2"};
 double xnt[nvar];
 NTuple nt(nvar,vname);

 long n = 0;
 double diffsx=0., diffsx2=0., diffsnx=0., diffsnx2=0.;
 for(double t=tmin; t<tmax+dt/3.; t+=dt) {
   xnt[0] = t;
   xnt[1] = sin(t);
   xnt[2] = sin(N*t);
   xnt[3] = SinXsX(t);
   xnt[4] = SinXsX_Sqr(t);
   xnt[5] = SinNXsX(t,N);
   xnt[6] = SinNXsX_Sqr(t,N);
   xnt[7] = SinXsX(t,true);
   xnt[8] = SinXsX_Sqr(t,true);
   xnt[9] = SinNXsX(t,N,true);
   xnt[10] = SinNXsX_Sqr(t,N,true);
   if(!only_test) nt.Fill(xnt);

   // les tests de difference maxi
   double av;
   if(xnt[0]>1e-15) {
     double ref = xnt[1]/xnt[0];
     av = fabs(xnt[3]-ref);
     if(av>pr_test) printf("t=%.15e sx/x=%.15e, %.15e d=%g\n",t,xnt[3],ref,av);
     if(av>diffsx) diffsx = av;
     av = fabs(xnt[4]-ref*ref);
     if(av>pr_test) printf("t=%.15e (sx/x)^2=%.15e, %.15e d=%g\n",t,xnt[4],ref*ref,av);
     if(av>diffsx2) diffsx2 = av;
   }
   if(fabs(xnt[1])>1e-15) {
     double ref = xnt[2]/xnt[1];
     av = fabs(xnt[5]-ref);
     if(av>pr_test) printf("t=%.15e snx/sx=%.15e, %.15e d=%g\n",t,xnt[5],ref,av);
     if(fabs(av)>diffsnx) diffsnx = av;
     av = fabs(xnt[6]-ref*ref);
     if(av>pr_test) printf("t=%.15e (snx/sx)^2=%.15e, %.15e d=%g\n",t,xnt[6],ref*ref,av);
     if(fabs(av)>diffsnx2) diffsnx2 = av;
   }

   n++;
 }
 cout<<"Number of entries = "<<n<<endl;
 cout<<"Diff  sx/x       = "<<diffsx<<endl;
 cout<<"     (sx/x)^2    = "<<diffsx2<<endl;
 //Attention, vers Pi c'est la precision machine sur sin(N*x)/sin(x) qui donne l'erreur
 cout<<"      sNx/x    = "<<diffsnx<<" = "<<diffsnx/N<<" * "<<N<<endl;
 cout<<"     (sNx/x)^2 = "<<diffsnx2<<" = "<<diffsnx2/(N*N)<<" * "<<N<<" * "<<N<<endl;

 if(!only_test) {
   cout<<"writing ppf"<<endl;
   POutPersist pos("cmvsinxsx.ppf");
   pos.PutObject(nt,"nt");
 }

 return 0;
}

/*
openppf cmvsinxsx.ppf

zone
n/plot nt.s%t ! ! "connectpoints"
n/plot nt.sn%t ! ! "connectpoints same red"

#-------
zone 1 3
n/plot nt.sx%t ! ! "connectpoints"
n/plot nt.asx%t ! ! "connectpoints same green"
n/plot nt.sx-s/t%t fabs(t)>1e-15 ! "connectpoints"
n/plot nt.sx-asx%t ! ! "connectpoints"

#-------
zone 1 3
n/plot nt.sx2%t ! ! "connectpoints"
n/plot nt.sx*sx%t ! ! "connectpoints same red"
n/plot nt.asx2%t ! ! "connectpoints same green"
n/plot nt.sx2-pow(s/t,2.)%t fabs(t)>1e-15 ! "connectpoints"
n/plot nt.sx2-asx2%t ! ! "connectpoints"

#-------
zone 1 3
n/plot nt.snx%t ! ! "connectpoints"
n/plot nt.asnx%t ! ! "connectpoints same green"
n/plot nt.snx-sn/s%t fabs(s)>1e-15 ! "connectpoints"
n/plot nt.snx-asnx%t ! ! "connectpoints"

#-------
zone 1 3
n/plot nt.snx2%t ! ! "connectpoints"
n/plot nt.snx*snx%t ! ! "connectpoints same red"
n/plot nt.asnx2%t ! ! "connectpoints same green"
n/plot nt.snx2-pow(sn/s,2.)%t fabs(s)>1e-15 ! "connectpoints"
n/plot nt.snx2-asnx2%t ! ! "connectpoints"

*/