Changeset 3500 in Sophya for trunk/Cosmo/SimLSS/geneutils.cc

Timestamp:

Jun 18, 2008, 6:44:59 PM (17 years ago)

Author:

cmv

Message:

fct sin(x)/x et sin(n*x)/sin(x) cmv 18/06/2008

File:

• trunk/Cosmo/SimLSS/geneutils.cc (modified) (1 diff)

trunk/Cosmo/SimLSS/geneutils.cc

@@ -294 +294 @@
 
 //-------------------------------------------------------------------
+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;