changes in lensings quadrature and in transfer k list

Makefile

@@ -77,6 +77,8 @@ TEST_LOOPS = test_loops.o
 
 TEST_SPECTRA = test_spectra.o
 
+TEST_LENSING = test_lensing.o
+
 TEST_TRANSFER = test_transfer.o
 
 TEST_BESSEL = test_bessel.o
@@ -94,13 +96,16 @@ TEST_KARIM = test_karim.o
 TEST_PBC = test_pbc.o
 
 all: class libclass.a
-	
+
 libclass.a: $(TOOLS) $(INPUT) $(BACKGROUND) $(THERMO) $(PERTURBATIONS) $(BESSEL) $(TRANSFER) $(PRIMORDIAL) $(SPECTRA) $(NONLINEAR) $(LENSING) $(OUTPUT)
 	$(AR)  $@ $(addprefix build/,$(notdir $^))
-	
+
 class: $(TOOLS) $(INPUT) $(BACKGROUND) $(THERMO) $(PERTURBATIONS) $(BESSEL) $(TRANSFER) $(PRIMORDIAL) $(SPECTRA) $(NONLINEAR) $(LENSING) $(OUTPUT) $(CLASS)
 	$(CC) $(LDFLAG) -o  $@ $(addprefix build/,$(notdir $^)) -lm
 
+test_lensing: $(TOOLS) $(INPUT) $(BACKGROUND) $(THERMO) $(PERTURBATIONS) $(BESSEL) $(TRANSFER) $(PRIMORDIAL) $(SPECTRA) $(NONLINEAR) $(LENSING) $(OUTPUT) $(TEST_LENSING)
+	$(CC) $(LDFLAG) -o  $@ $(addprefix build/,$(notdir $^)) -lm
+
 test_pbc: $(TOOLS) $(INPUT) $(BACKGROUND) $(THERMO) $(PERTURBATIONS) $(BESSEL) $(TRANSFER) $(PRIMORDIAL) $(SPECTRA) $(LENSING) $(OUTPUT) $(TEST_PBC)
 	$(CC) $(LDFLAG) -o  $@ $(addprefix build/,$(notdir $^)) -lm
 

REFCLASS.pre

@@ -54,11 +54,11 @@ thermo_rate_smoothing_radius=50
 
 gauge=1 #synchronous
 
-k_scalar_min_eta0=0.003
-k_scalar_max_eta0_over_l_max=5.
-k_scalar_step_sub=0.01
-k_scalar_step_super=0.0005
-k_scalar_step_transition=0.2
+k_scalar_min_eta0=0.002
+k_scalar_max_eta0_over_l_max=3.
+k_scalar_step_sub=0.015
+k_scalar_step_super=0.0001 # was 0.0005
+k_scalar_step_transition=0.2 # was 0.2
 
 #start_small_k_at_eta_g_over_eta_h = 0.006
 #start_large_k_at_eta_g_over_eta_k = 1.e-5
@@ -79,7 +79,7 @@ l_max_ncdm1=28
 
 tol_eta_approx=1.e-5
 tol_perturb_integration=1.e-6
-perturb_sampling_stepsize=0.005
+perturb_sampling_stepsize=0.01
 
 free_streaming_approximation = 2
 free_streaming_trigger_eta_h_over_eta_k = 120.
@@ -88,21 +88,21 @@ free_streaming_trigger_Omega_r = 0.07
 l_logstep=1.026
 l_linstep=25
 
-bessel_x_step=0.01
-bessel_j_cut=5.e-10   # WATCH IT!
+bessel_x_step=0.005
+bessel_j_cut=5.e-10  # WATCH IT!
 bessel_delta_x_min =1.e-4
 bessel_file_name=bessel_large.dat
 
 k_per_decade_primordial = 10.
 
-k_step_trans_scalars=0.1 # 0.04 to smooth oscillations at l<200, but 0.2 reasonnable; smaller -> job killed on superb
+k_step_trans_scalars=0.002
 
-transfer_cut=2 #0=none,1=osc,2=cl #segfault with zero, watch it!
+transfer_cut=1 #0=none,1=osc,2=cl #segfault with zero, watch it!
 transfer_cut_threshold_osc=0.005 # with 0.005, more smooth than when smaller for l>2000 WATCH IT
-transfer_cut_threshold_cl=1.e-8
+transfer_cut_threshold_cl=1.e-8 # not robust if tilted or if one Delta_l(k) oscillation tangents zero
 
 evolver=0
 
 l_switch_limber = 40.
 num_mu_minus_lmax = 1000.
-delta_l_max = 2000.
+delta_l_max = 1000.

explanatory.ini

@@ -125,7 +125,7 @@ non linear = no
 
 3) if you want the spectrum of lensed Cls, enter a word containing the letter 'y' or 'Y' (default: no lensed Cls)
 
-lensing = no
+lensing = yes
 
 4) list of modes ('s' for scalars, 'v' for vectors, 't' for tensors). More than one letter allowed, can be attached or separated by arbitrary characters; letters can be small or capital. (default: set to 's')
 

include/common.h

@@ -532,11 +532,8 @@ struct precision
 
   //@{
 
-  double k_step_trans_scalars; /**< sampling step in k space, in units of \f$ 2\pi/(\eta_0-\eta_{rec}) \f$, which is the typical period of oscillations of \f$ \Delta_l(k) \f$ */
-  double k_step_trans_tensors; /**< sampling step in k space, in units of \f$ 2\pi/(\eta_0-\eta_{rec}) \f$, which is the typical period of oscillations of \f$ \Delta_l(k) \f$ */
-
-/*   int k_oversampling_scalars; */
-/*   int k_oversampling_tensors; */
+  double k_step_trans_scalars; /**< upper bound on linear sampling step in k space, in units of one period of acoustic oscillation at decoupling (usually chosen to be between k_scalar_step_sub and k_scalar_step_super) */
+  double k_step_trans_tensors; /**< upper bound on linear sampling step in k space, in units of one period of acoustic oscillation at decoupling (usually chosen to be between k_tensor_step_sub and k_tensor_step_super) */
 
   enum transfer_cutting transfer_cut; /**< flag telling how to cut automatically the transfer function computation at a given \f$ k_{max} \f$ value */
 
@@ -546,8 +543,6 @@ struct precision
 
   int l_switch_limber;
 
-  double smallest_allowed_variation; /**< machine-dependent, assigned automatically by the code */
-
   //@}
 
   /** @name - parameters related to lensing */
@@ -556,7 +551,7 @@ struct precision
 
   int num_mu_minus_lmax; /**< difference between num_mu and l_max, increase for more precision */
   int delta_l_max; /**<difference between l_max in unlensed and lensed spectra */
-
+  double tol_gauss_legendre; /**< tolerance with which quadrature points are found: must be very small for an accurate integration (if not entered manually, set automatically to match machine precision) */
   //@}
 
 
@@ -582,6 +577,14 @@ struct precision
 
   //@}
 
+  /** @name - general precision parameters */
+
+  //@{
+
+  double smallest_allowed_variation; /**< machine-dependent, assigned automatically by the code */
+
+  //@}
+
   /** @name - zone for writing error messages */
 
   //@{

include/lensing.h

@@ -284,14 +284,7 @@ extern "C" {
                    int lmax,
                    double ** d4m4
                    );
-
-  int lensing_gauss_legendre(
-                             double *mu,
-                             double *w8,
-                             int nmu
-                             );
-
-
+
 #ifdef __cplusplus
 }
 #endif

include/quadrature.h

@@ -70,7 +70,12 @@ extern "C" {
 	      double b);
   double testfun(double x);
 
-
+int quadrature_gauss_legendre(
+                           double *mu,
+                           double *w8,
+                           int n,
+			   double tol,
+			   ErrorMsg error_message);
 
 
 #ifdef __cplusplus

include/transfer.h

@@ -121,8 +121,7 @@ extern "C" {
 				    struct perturbs * ppt,
 				    struct bessels * pbs,
 				    struct transfers * ptr,
-				    double eta0,
-				    double eta_rec
+				    double rs_rec
 				    );
 
   int transfer_get_l_list(
@@ -137,8 +136,7 @@ extern "C" {
 			  struct precision * ppr,
 			  struct perturbs * ppt,
 			  struct transfers * ptr,
-			  double eta0,
-			  double eta_rec,
+			  double rs_rec,
 			  int index_mode
 			  );
 

source/input.c

@@ -1059,9 +1059,6 @@ int input_init(
 
   class_read_double("k_step_trans_scalars",ppr->k_step_trans_scalars);
   class_read_double("k_step_trans_tensors",ppr->k_step_trans_tensors);
-  /*   class_read_int("k_oversampling_scalars",ppr->k_oversampling_scalars); */
-  /*   class_read_int("k_oversampling_tensors",ppr->k_oversampling_tensors); */
-
   class_read_int("transfer_cut",ppr->transfer_cut);
   class_read_double("transfer_cut_threshold_osc",ppr->transfer_cut_threshold_osc);
   class_read_double("transfer_cut_threshold_cl",ppr->transfer_cut_threshold_cl);
@@ -1089,6 +1086,7 @@ int input_init(
 
   class_read_int("num_mu_minus_lmax",ppr->num_mu_minus_lmax);
   class_read_int("delta_l_max",ppr->delta_l_max);
+  class_read_int("tol_gauss_legendre",ppr->tol_gauss_legendre);
 
 
   /* check various l_max */
@@ -1484,12 +1482,8 @@ int input_default_precision ( struct precision * ppr ) {
    * - parameter related to the transfer functions
    */
 
-  ppr->k_step_trans_scalars=0.2; /* 03.12.10 for chi2plT0.01: difficult to optimize, numerical instability below 0.1, need to study this better */
-  ppr->k_step_trans_tensors=0.15;
-
-  /*   ppr->k_oversampling_scalars=3; */
-  /*   ppr->k_oversampling_tensors=1; */
-
+  ppr->k_step_trans_scalars=0.002;
+  ppr->k_step_trans_tensors=0.0015;
   ppr->transfer_cut=tc_cl; /* 03.12.10 for chi2plT0.01: tc_cl slightly faster (by 20%) for equal precision, but also slightly less robust (better to switch to tc_osc if primordial tilt can depart significantly from one) */
   ppr->transfer_cut_threshold_osc=0.015; /* 03.12.10 for chi2plT0.01 */
   ppr->transfer_cut_threshold_cl=2.e-6; /* 14.12.10 for chi2plT0.01 */
@@ -1533,6 +1527,8 @@ int input_default_precision ( struct precision * ppr ) {
 	     ppr->error_message,
 	     "smallest_allowed_variation = %e < 0",ppr->smallest_allowed_variation);
 
+  ppr->tol_gauss_legendre = ppr->smallest_allowed_variation;
+
   return _SUCCESS_;
 
 }

source/lensing.c

@@ -1,11 +1,13 @@
 /** @file lensing.c Documented lensing module
  *
- * Julien Lesgourgues and Simon Prunet, 6.12.2010    
+ * Simon Prunet and Julien Lesgourgues, 6.12.2010    
  *
  * This module computes the lensed temperature and polarization
  * anisotropy power spectra \f$ C_l^{X}, P(k), ... \f$'s given the
  * unlensed temperature, polarization and lensing potential spectra.
  *
+ * Follows Challinor and Lewis full-sky method, astro-ph/0502425
+ *
  * The following functions can be called from other modules:
  *
  * -# lensing_init() at the beginning (but after spectra_init())
@@ -17,7 +19,7 @@
 
 /** 
  * Anisotropy power spectra C_l's for all types, modes and initial conditions.
- * SO FAR: ONLY SCALAR ADIABATIC 
+ * SO FAR: ONLY SCALAR
  *
  * This routine evaluates all the lensed C_l's at a given value of l by
  * picking it in the pre-computed table. When relevant, it also
@@ -148,6 +150,8 @@ int lensing_init(
   num_mu += num_mu%2; /* Force it to be even */ 
   /** - allocate array of mu values, as well as quadrature weights */
 
+  printf("num_mu=%d\n",num_mu);
+
   class_alloc(mu,
               num_mu*sizeof(double),
               ple->error_message);
@@ -158,7 +162,11 @@ int lensing_init(
               (num_mu-1)*sizeof(double),
               ple->error_message);
 
-  class_call(lensing_gauss_legendre(mu,w8,num_mu-1),
+  class_call(quadrature_gauss_legendre(mu,
+				    w8,
+				    num_mu-1,
+				    ppr->tol_gauss_legendre,
+				    ple->error_message),
              ple->error_message,
              ple->error_message); 
 
@@ -1934,52 +1942,3 @@ int lensing_d4m4(
   free(fac1); free(fac2); free(fac3); free(fac4);
   return _SUCCESS_;
 }
-
-/**
- * This routine computes the weights and abscissas of a Gauss-Legendre quadrature
- *
- * @param mu     Input/output: Vector of cos(beta) values
- * @param w8     Input/output: Vector of quadrature weights
- * @param n      Input       : Number of quadrature points
- *
- * From Numerical recipes 
- **/
-
-int lensing_gauss_legendre(
-                           double *mu,
-                           double *w8,
-                           int n) {
-
-  int m,j,i;
-  double z1,z,xm,xl,pp,p3,p2,p1;
-  double x1,x2,EPS;
-  x1 = -1.0;
-  x2 = 1.0;
-  EPS = 1.0e-11;
-
-  m=(n+1)/2;
-  xm=0.5*(x2+x1);
-  xl=0.5*(x2-x1);
-  for (i=1;i<=m;i++) {
-    z=cos(_PI_*(i-0.25)/(n+0.5));
-    do {
-      p1=1.0;
-      p2=0.0;
-      for (j=1;j<=n;j++) {
-        p3=p2;
-        p2=p1;
-        p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
-      }
-      pp=n*(z*p1-p2)/(z*z-1.0);
-      z1=z;
-      z=z1-p1/pp;
-    } while (fabs(z-z1) > EPS);
-    mu[i-1]=xm-xl*z;
-    mu[n-i]=xm+xl*z;
-    w8[i-1]=2.0*xl/((1.0-z*z)*pp*pp);
-    w8[n-i]=w8[i-1];
-
-  }
-  return _SUCCESS_;
-}
-