Index: /trunk/Cosmo/SimLSS/Makefile
===================================================================
--- /trunk/Cosmo/SimLSS/Makefile	(revision 3499)
+++ /trunk/Cosmo/SimLSS/Makefile	(revision 3500)
@@ -38,10 +38,11 @@
 # ---- Les programmes de test
 PROGSTEST = \
-        $(EXE)cmvchkwhu $(EXE)hu_sigma8
+        $(EXE)cmvchkwhu $(EXE)hu_sigma8 $(EXE)cmvsinxsx
 #$(EXE)cmvtluc
 
 PROGSTESTOBJ = \
         $(OBJ)hu_tf_fit.o $(OBJ)hu_power.o \
-        $(OBJ)cmvtluc.o $(OBJ)cmvchkwhu.o $(OBJ)hu_sigma8.o
+        $(OBJ)cmvtluc.o $(OBJ)cmvchkwhu.o $(OBJ)hu_sigma8.o \
+        $(OBJ)cmvsinxsx.o
 
 #--------------------------------------------------------------------------
@@ -240,2 +241,12 @@
 $(OBJ)cmvfitpk.o: cmvfitpk.cc
 	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvfitpk.cc
+
+##############################################################################
+cmvsinxsx: $(EXE)cmvsinxsx
+	echo $@ " done"
+$(EXE)cmvsinxsx: $(OBJ)cmvsinxsx.o $(LIBR)
+	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvsinxsx.o $(MYLIB) 
+$(OBJ)cmvsinxsx.o: cmvsinxsx.cc
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvsinxsx.cc
+
+##############################################################################
Index: /trunk/Cosmo/SimLSS/cmvsinxsx.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvsinxsx.cc	(revision 3500)
+++ /trunk/Cosmo/SimLSS/cmvsinxsx.cc	(revision 3500)
@@ -0,0 +1,163 @@
+#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"
+
+*/
Index: /trunk/Cosmo/SimLSS/geneutils.cc
===================================================================
--- /trunk/Cosmo/SimLSS/geneutils.cc	(revision 3499)
+++ /trunk/Cosmo/SimLSS/geneutils.cc	(revision 3500)
@@ -294,4 +294,66 @@
 
 //-------------------------------------------------------------------
+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
+{
+  double x2 = x*x;
+  if(app || x2<1.7e-4) return 1. - x2/6.*(1. - x2/20.*(1. - x2/42.));
+  return sin(x)/x;
+}
+
+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
+{
+  double x2 = x*x;
+  if(app || x2<6.8e-5) return 1. - x2/3.*(1. - 2.*x2/15.*(1. - x2/14.));
+  x2 = sin(x)/x;
+  return x2*x2;
+}
+
+double SinNXsX(double x,unsigned long N,bool app)
+// Calcul de sin(N*x)/sin(x)  avec N entier positif
+// ATTENTION: N est entier 
+//  1. pour que sin(N*x) et sin(x) soient nuls en meme temps
+//     (on peut faire le DL popur sin(N*x) et sin(x))
+//  2. pour traiter correctement le DL en x = p*Pi+e  avec e<<1 et p entier
+//     sin(N*x)/sin(x) = sin(N*p*Pi+N*e)/sin(p*Pi+e)
+//                     = [sin(N*p*Pi)*cos(N*e)+cos(N*p*Pi)*sin(N*e)]
+//                         / [sin(p*Pi)*cos(e)+cos(p*Pi)*sin(e)]
+//                     comme sin(N*p*Pi)=0
+//                     = [cos(N*p*Pi)*sin(N*e)] / [cos(p*Pi)*sin(e)]
+//                     = [sin(N*e)/sin(e)] * [cos(N*p*Pi)/cos(p*Pi)]
+//                     = [DL autour de x=0] * (+1 ou -1)
+//     Le seul cas ou on a "-1" est quand "p=impair" (cos(p*Pi)=-1) et "N=pair"
+// Approx: dernier terme x^4*N^4/120 ~ 1e-13 -> x^2 ~ 3.5e-6/N^2
+{
+  if(N==0) return 0;
+  double sx = sin(x), N2 = N*N;
+  if(app || sx*sx<3.5e-6/N2) {
+    double x2 = asin(sx); x2 *= x2;
+    double s = 1.;
+    if(N%2==0 && cos(x)<0.) s = -1.; // cas x ~ (2p+1)*Pi et N pair
+    return s*N*(1.-(N-1.)*(N+1.)/6.*x2*(1.-(3.*N2-7.)/60.*x2));
+  }
+  return sin((double)N*x)/sx;
+}
+
+double SinNXsX_Sqr(double x,unsigned long N,bool app)
+// Calcul de [sin(N*x)/sin(x)]^2  avec N entier positif
+// ATTENTION: cf remarque pour N entier dans SinNXsX
+// Approx: dernier terme x^4*2*N^4/45 ~ 1e-13 -> x^2 ~ 1.5e-6/N^2
+{
+  if(N==0) return 0;
+  double sx = sin(x), N2 = N*N;
+  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));
+  }
+  sx = sin((double)N*x)/sx;
+  return sx*sx;
+}
+
+//-------------------------------------------------------------------
 
 static unsigned short IntegrateFunc_GlOrder = 0;
Index: /trunk/Cosmo/SimLSS/geneutils.h
===================================================================
--- /trunk/Cosmo/SimLSS/geneutils.h	(revision 3499)
+++ /trunk/Cosmo/SimLSS/geneutils.h	(revision 3500)
@@ -88,4 +88,10 @@
 double FrAngSol(double angsol);
 
+double SinXsX(double x,bool app=false);
+double SinXsX_Sqr(double x,bool app=false);
+
+double SinNXsX(double x,unsigned long N,bool app=false);
+double SinNXsX_Sqr(double x,unsigned long N,bool app=false);
+
 //----------------------------------------------------
 double IntegrateFunc(GenericFunc& func,double xmin,double xmax
