source: Sophya/trunk/Cosmo/SimLSS/cmvsinxsx.cc@ 4049

#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 -x 1e-6,-0.1,0.1
> cmvsinxsx -n 5 -x 1e-6,-0.1,0.1
...Check autour de Pi
> cmvsinxsx -n 4 -x 1e-6,3.041592659,3.241592659
> cmvsinxsx -n 5 -x 1e-6,3.041592659,3.241592659
...Check autour de -Pi
> cmvsinxsx -n 4 -x 1e-6,-3.241592659,-3.041592659
> cmvsinxsx -n 5 -x 1e-6,-3.241592659,-3.041592659

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

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

...Check general
> cmvsinxsx -n 20 -p 1e-8 -v 65536,-2.,0. -x -1,-20.,20.
 */


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

 //--- Read arguments
 char c;
 while((c = getopt(narg,arg,"htp:n:v:x:")) != -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 'x' :
    sscanf(optarg,"%lf,%lf,%lf",&dt,&tmin,&tmax);
    break;
  case 'v' :
    sscanf(optarg,"%d,%lf,%lf",&nvec,&vecmax,&veccent);
    break;
  case 'h' :
  default :
    cout<<"cmvsinxsx [-t] [-p eps_print] [-n N] [-v nv,(+-)vmax,vcent (rad)]"<<endl;
    cout<<"          [dt,tmin,tmax (rad)]"<<endl;
    return -1;
     break;
 }
 }

 if(N<=0) N=1;
 // on calcule pour avoir 20 points par periode de sin(n*x)
 if(dt<=0.) dt = 2.*M_PI/(20.*N);

 cout<<"N="<<N<<" dt="<<dt<<" tmin="<<tmin<<" tmax="<<tmax<<" rad"<<endl;
 cout<<"only_test="<<only_test<<" pr_test="<<pr_test<<endl;
 if(vecmax<0.) vecmax *= -M_PI;
 cout<<"nv="<<nvec<<" vmax="<<vecmax<<" vcent="<<veccent
     <<" -> ["<<veccent-vecmax/2<<","<<veccent+vecmax/2<<"]"<<endl;
 if(tmax<=tmin) return -1;

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

 //--- Fill NTuple and do precision test
 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;

 //--- Fill vector for FFT
 TVector<r_8> vsx,vsx2,vsnx,vsnx2;
 if(nvec>1) {
   vsx.ReSize(nvec); vsx2.ReSize(nvec); vsnx.ReSize(nvec); vsnx2.ReSize(nvec);
   for(int i=0;i<nvec;i++) {
     double t = veccent -vecmax/2. + i*vecmax/(double)nvec;
     vsx(i)   = SinXsX(t);
     vsx2(i)  = SinXsX_Sqr(t);
     vsnx(i)  = SinNXsX(t,N);
     vsnx2(i) = SinNXsX_Sqr(t,N);
   }
 }

 //--- Write PPF
 POutPersist pos("cmvsinxsx.ppf");
 if(!only_test) pos.PutObject(nt,"nt");
 if(vsx.Size()>0)   pos.PutObject(vsx,"vsx");
 if(vsx2.Size()>0)  pos.PutObject(vsx2,"vsx2");
 if(vsnx.Size()>0)  pos.PutObject(vsnx,"vsnx");
 if(vsnx2.Size()>0) pos.PutObject(vsnx2,"vsnx2");

 return 0;
}

/*
openppf cmvsinxsx.ppf

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

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

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

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

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

#-------
openppf cmvsinxsx.ppf

zone 2 2
n/plot vsx.val%n ! ! 
n/plot vsx2.val%n ! ! 
n/plot vsnx.val%n ! ! 
n/plot vsnx2.val%n ! !

zone
fftforw vsx fvsx
fftforw vsx2 fvsx2
fftforw vsnx fvsnx
fftforw vsnx2 fvsnx2

zone
n/plot fvsx.val%n ! ! "nsta connectpoints"
n/plot fvsx2.val%n ! ! "nsta connectpoints"
n/plot fvsnx.val%n ! ! "nsta connectpoints"
n/plot fvsnx2.val%n ! ! "nsta connectpoints"

zone
n/plot fvsx.log10(mod*mod)%n mod>0. ! "nsta connectpoints"
n/plot fvsx2.log10(mod*mod)%n mod>0. ! "nsta connectpoints"
n/plot fvsnx.log10(mod*mod)%n mod>0. ! "nsta connectpoints"
n/plot fvsnx2.log10(mod*mod)%n mod>0. ! "nsta connectpoints"

*/