Changeset 3595 in Sophya

Timestamp:

Apr 19, 2009, 2:47:10 PM (16 years ago)

Author:

cmv

Message:

improve DL precison and tests main routine, cmv 19/04/2009

Location:

trunk/Cosmo/SimLSS

Files:

• cmvsinxsx.cc (modified) (8 diffs)
• geneutils.cc (modified) (4 diffs)

trunk/Cosmo/SimLSS/cmvsinxsx.cc

@@ -12 +12 @@
 /*
 ...Check autour de zero
-> cmvsinxsx -n 4 1e-6,-0.1,0.1
-> cmvsinxsx -n 5 1e-6,-0.1,0.1
+> 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 1e-6,3.041592659,3.241592659
-> cmvsinxsx -n 5 1e-6,3.041592659,3.241592659
+> 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 1e-6,-3.241592659,-3.041592659
-> cmvsinxsx -n 5 1e-6,-3.241592659,-3.041592659
+> 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 0.0001,-6.5,6.5
-> cmvsinxsx -n 11 0.0001,-6.5,6.5
+> 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 0.00001,-6.5,6.5
-> cmvsinxsx -n 1000 0.00001,-6.5,6.5
+> 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.
  */
 
@@ -34 +37 @@
 {
  double dt = 0.000001;
- //double tmin = -0.1, tmax = 0.1;
- double tmin = M_PI-0.1, tmax = M_PI+0.1;
+ 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;
@@ -44 +47 @@
  //--- Read arguments
  char c;
- while((c = getopt(narg,arg,"htp:n:v:")) != -1) {
+ while((c = getopt(narg,arg,"htp:n:v:x:")) != -1) {
   switch (c) {
   case 'n' :
@@ -55 +58 @@
     sscanf(optarg,"%lf",&pr_test);
     break;
+  case 'x' :
+    sscanf(optarg,"%lf,%lf,%lf",&dt,&tmin,&tmax);
+    break;
   case 'v' :
-    sscanf(optarg,"%d,%lf",&nvec,&vecmax);
+    sscanf(optarg,"%d,%lf,%lf",&nvec,&vecmax,&veccent);
     break;
   case 'h' :
   default :
-    cout<<"cmvsinxsx [-t] [-p eps_print] [-n N] [-v nv,(+-)vmax (rad)]"<<endl;
+    cout<<"cmvsinxsx [-t] [-p eps_print] [-n N] [-v nv,(+-)vmax,vcent (rad)]"<<endl;
     cout<<"          [dt,tmin,tmax (rad)]"<<endl;
     return -1;
-  }
- }
- if(optind<narg) sscanf(arg[optind],"%lf,%lf,%lf",&dt,&tmin,&tmax);
+     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<<" -> ["<<-vecmax/2<<","<<vecmax/2<<"]"<<endl;
+ cout<<"nv="<<nvec<<" vmax="<<vecmax<<" vcent="<<veccent
+     <<" -> ["<<veccent-vecmax/2<<","<<veccent+vecmax/2<<"]"<<endl;
  if(tmax<=tmin) return -1;
 
@@ -83 +94 @@
  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[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);
@@ -142 +153 @@
  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(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");
 
@@ -154 +165 @@
 
 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"
+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"
 
 #-------
@@ -203 +214 @@
 
 zone
-n/plot fvsx.val%n ! ! "connectpoints"
-n/plot fvsx2.val%n ! ! "connectpoints"
-n/plot fvsnx.val%n ! ! "connectpoints"
-n/plot fvsnx2.val%n ! ! "connectpoints"
+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"
 
 */

trunk/Cosmo/SimLSS/geneutils.cc

@@ -296 +296 @@
 double SinXsX(double x,bool app)
 // Calcul de sin(x)/x
-// Approx: dernier terme x^6/(6*20*42) ~ 1e-13 -> x^2 ~ 1.7e-4
+// Le Dl est en O[x]^10  (cad on va jusqu'au terme en x^8 compris)
+// Approx: terme en x^6/(6*20*42) ~ 1e-13 -> x^2 ~ 1.7e-4
 {
   double x2 = x*x;
-  if(app || x2<1.7e-4) return 1. - x2/6.*(1. - x2/20.*(1. - x2/42.));
+  if(app || x2<1.7e-4) return 1.-x2/6.*(1.-x2/20.*(1.-x2/42.*(1.-x2/72.)));
   return sin(x)/x;
 }
@@ -305 +306 @@
 double SinXsX_Sqr(double x,bool app)
 // Calcul de (sin(x)/x)^2
-// Approx: dernier terme 2*x^6/(3*15*14) ~ 1e-13 -> x^2 ~ 6.8e-5
+// Le Dl est en O[x]^10  (cad on va jusqu'au terme en x^8 compris)
+// Approx: terme 2*x^6/(3*15*14) ~ 1e-13 -> x^2 ~ 6.8e-5
 {
   double x2 = x*x;
-  if(app || x2<6.8e-5) return 1. - x2/3.*(1. - 2.*x2/15.*(1. - x2/14.));
+  if(app || x2<6.8e-5) return 1.-x2/3.*(1.-2.*x2/15.*(1.-x2/14.*(1.-2.*x2/45.)));
   x2 = sin(x)/x;
   return x2*x2;
@@ -343 +345 @@
 // Calcul de [sin(N*x)/sin(x)]^2  avec N entier positif
 // ATTENTION: cf remarque pour N entier dans SinNXsX
+//            maximum de la fonction: N^2
 // Approx: dernier terme x^4*2*N^4/45 ~ 1e-13 -> x^2 ~ 1.5e-6/N^2
 {
@@ -349 +352 @@
   if(app || sx*sx<1.5e-6/N2) {
     double x2 = asin(sx); x2 *= x2;
-    return N2*(1. - (N-1.)*(N+1.)/3.*x2*(1. - (2.*N2-3.)/15.*x2));
+    //return N2*(1.-(N*N-1.)/3.*x2);
+    return N2*(1.-(N-1.)*(N+1.)/3.*x2*(1.-(2.*N2-3.)/15.*x2));
   }
   sx = sin((double)N*x)/sx;