corrected bug for computing lCl alone, changed management of xmax in …

…bessel

REFCLASS.pre

@@ -55,7 +55,7 @@ thermo_rate_smoothing_radius=50
 gauge=1 #synchronous
 
 k_scalar_min_eta0=0.003
-k_scalar_max_eta0_over_l_max=3.
+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
@@ -79,7 +79,7 @@ l_max_ncdm1=28
 
 tol_eta_approx=1.e-5
 tol_perturb_integration=1.e-6
-perturb_sampling_stepsize=0.01
+perturb_sampling_stepsize=0.005
 
 free_streaming_approximation = 2
 free_streaming_trigger_eta_h_over_eta_k = 120.
@@ -91,19 +91,18 @@ l_linstep=25
 bessel_x_step=0.01
 bessel_j_cut=5.e-10   # WATCH IT!
 bessel_delta_x_min =1.e-4
-bessel_x_max_over_l_max=1.8
 bessel_file_name=bessel_large.dat
 
 k_per_decade_primordial = 10.
 
-k_step_trans_scalars=0.04 # 0.04 to smooth oscillations at l<200, but 0.2 reasonnable; smaller -> job killed on superb
+k_step_trans_scalars=0.1 # 0.04 to smooth oscillations at l<200, but 0.2 reasonnable; smaller -> job killed on superb
 
 transfer_cut=2 #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
 
 evolver=0
 
-l_switch_limber = 10.
-num_mu_minus_lmax = 400.
-delta_l_max = 100.
+l_switch_limber = 40.
+num_mu_minus_lmax = 1000.
+delta_l_max = 2000.

REFCLASS_tClpCl.pre

@@ -0,0 +1,108 @@
+# this precision file obtained as follows:
+# - computing adiabatic scalar unlensed ClT only (ndf15 integrator with AMD)
+# - for each precision parameter varied individually, Deltachi2(Planck) between high-precision limit and selected value = 1.e-3  
+
+a_ini_over_a_today_default = 1.e-14
+back_integration_stepsize = 7.e-3
+tol_background_integration = 1.e-2
+
+recfast_z_initial=1.e4
+
+recfast_Nz0=100000
+tol_thermo_integration=1.e-5
+
+recfast_Heswitch=6
+recfast_fudge_He=0.86
+
+recfast_Hswitch = 1 #_TRUE_
+recfast_fudge_H = 1.14
+recfast_delta_fudge_H = -0.035
+recfast_AGauss1 = -0.14
+recfast_AGauss2 =  0.05
+recfast_zGauss1 =  7.28
+recfast_zGauss2 =  6.75
+recfast_wGauss1 =  0.18
+recfast_wGauss2 =  0.33
+
+recfast_z_He_1 = 8000.
+recfast_delta_z_He_1 = 50.
+recfast_z_He_2 = 5000.
+recfast_delta_z_He_2 = 100.
+recfast_z_He_3 = 3500.
+recfast_delta_z_He_3 = 50.
+recfast_x_He0_trigger = 0.995
+recfast_x_He0_trigger2 = 0.995
+recfast_x_He0_trigger_delta = 0.01
+recfast_x_H0_trigger = 0.995
+recfast_x_H0_trigger2 = 0.995
+recfast_x_H0_trigger_delta = 0.01
+
+recfast_H_frac=1.e-3
+
+reionization_z_start_max = 50.
+reionization_sampling=1.e-2
+reionization_optical_depth_tol=1.e-2
+reionization_exponent=1.5
+
+reionization_width=0.5
+
+reionization_start_factor=8.
+helium_fullreio_redshift=3.5
+helium_fullreio_width=0.5
+
+thermo_rate_smoothing_radius=50
+
+gauge=1 #synchronous
+
+k_scalar_min_eta0=0.003
+k_scalar_max_eta0_over_l_max=3.
+k_scalar_step_sub=0.01
+k_scalar_step_super=0.0005
+k_scalar_step_transition=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
+#tight_coupling_trigger_eta_g_over_eta_h=0.008
+#tight_coupling_trigger_eta_g_over_eta_k=0.05
+#start_sources_at_eta_g_over_eta_h = 0.01
+start_small_k_at_eta_g_over_eta_h = 0.0004
+start_large_k_at_eta_h_over_eta_k = 0.15
+tight_coupling_trigger_eta_g_over_eta_h=0.005
+tight_coupling_trigger_eta_g_over_eta_k=0.008
+start_sources_at_eta_g_over_eta_h = 0.006
+tight_coupling_approximation=5 #(int)second_order_CRS;
+
+l_max_g=25
+l_max_pol_g=25
+l_max_nur=35
+l_max_ncdm1=28
+
+tol_eta_approx=1.e-5
+tol_perturb_integration=1.e-6
+perturb_sampling_stepsize=0.01
+
+free_streaming_approximation = 2
+free_streaming_trigger_eta_h_over_eta_k = 120.
+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_delta_x_min =1.e-4
+bessel_file_name=bessel_large.dat
+
+k_per_decade_primordial = 10.
+
+k_step_trans_scalars=0.04 # 0.04 to smooth oscillations at l<200, but 0.2 reasonnable; smaller -> job killed on superb
+
+transfer_cut=2 #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
+
+evolver=0
+
+l_switch_limber = 10.
+num_mu_minus_lmax = 1000.
+delta_l_max = 2000.

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 = yes
+lensing = no
 
 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

@@ -516,7 +516,6 @@ struct precision
   double bessel_x_step; /**< step dx for sampling Bessel functions \f$ j_l(x) \f$ */
   double bessel_j_cut; /**< value of \f$ j_l \f$ below which it is approximated by zero (in the region \f$ x \ll l \f$) */
   double bessel_tol_x_min;  /**< precision with which x_min such that j_l(x_min)=j_cut is found (order of magnitude set by k_min) */
-  double bessel_x_max_over_l_max; /**< x_max is defined as this parameter times largest l in Cls computation */ 
   FileName bessel_file_name; /**< name of file where Bessel functions will evnetually be written or read */
 
   //@}

source/bessel.c

@@ -163,7 +163,7 @@ int bessel_init(
 	     pbs->error_message,
 	     "x_step=%e, stop to avoid segmentation fault",pbs->x_step);
 
-  pbs->x_max = ((int)(pbs->l_max * ppr->bessel_x_max_over_l_max / pbs->x_step)+1)*pbs->x_step;
+  pbs->x_max = ((int)(pbs->l_max * ppr->k_scalar_max_eta0_over_l_max / pbs->x_step)+1)*pbs->x_step;
 
   pbs->j_cut = ppr->bessel_j_cut;
 

source/input.c

@@ -1049,7 +1049,6 @@ int input_init(
   class_read_double("bessel_x_step",ppr->bessel_x_step);
   class_read_double("bessel_j_cut",ppr->bessel_j_cut);
   class_read_double("bessel_tol_x_min",ppr->bessel_tol_x_min);
-  class_read_double("bessel_x_max_over_l_max",ppr->bessel_x_max_over_l_max);
   class_read_string("bessel_file_name",ppr->bessel_file_name);
 
   /** h.5. parameter related to the primordial spectra */
@@ -1473,7 +1472,6 @@ int input_default_precision ( struct precision * ppr ) {
   ppr->bessel_x_step=0.045; /* 14.12.10 for chi2plT0.1 */
   ppr->bessel_j_cut=1.5e-4; /* 14.12.10 for chi2plT0.1 */
   ppr->bessel_tol_x_min =1.e-4; /* 03.12.10 for chi2plT0.01 */
-  ppr->bessel_x_max_over_l_max=1.8; /* 03.12.10 for chi2plT0.01 */
   sprintf(ppr->bessel_file_name,"bessels.dat");
 
   /**

source/perturbations.c

@@ -796,8 +796,32 @@ int perturb_timesampling_for_sources(
     eta_ini = eta_mid;
 
   }
-  else 
+  else {
+
+    /* case when CMB not requested: start at recombination time */
     eta_ini = pth->eta_rec;
+
+    /* set values of first_index_back/thermo */
+    class_call(background_at_eta(pba,
+				 eta_ini, 
+				 short_info, 
+				 normal, 
+				 &first_index_back, 
+				 pvecback),
+	       pba->error_message,
+	       ppt->error_message);
+
+    class_call(thermodynamics_at_z(pba,
+				   pth,
+				   1./pvecback[pba->index_bg_a]-1.,  /* redshift z=1/a-1 */
+				   normal,
+				   &first_index_thermo,
+				   pvecback,
+				   pvecthermo),
+	       pth->error_message,
+	       ppt->error_message);
+  }    
+
 
   counter = 1;
 

source/transfer.c

@@ -415,6 +415,18 @@ int transfer_init(
 	       in the workspace) */
 	    ddj_l = j_l + x_size_l;
 
+	    /* check that the computation will never need values of
+	       j_l(x) with x > x_max (should never happen, since x_max
+	       is chosen to be greater than eta0*k_max in bessel
+	       module) */
+
+	    class_test_parallel((int)((eta0_minus_eta[0] * ptr->k[index_mode][ptr->k_size[index_mode]-1] - (*x_min_l))/pbs->x_step)+1 >= x_size_l,
+				ptr->error_message,
+				"Increase x_max in bessel functions! The computation needs index_x up to %d while x_size[%d]=%d\n",
+				(int)((eta0_minus_eta[0] * ptr->k[index_mode][ptr->k_size[index_mode]-1] - (*x_min_l))/pbs->x_step)+1,
+				index_l,
+				x_size_l);
+
 	    /* compute the transfer function for this l */
 	    class_call_parallel(transfer_compute_for_each_l(ppr,
 							    ppt,