now, when CMB is requested, z_pk or z_max_pk can be arbitrarily large

abacusorg · Jun 14, 2015 · 4c8c09a · 4c8c09a
1 parent eb2264b
commit 4c8c09a
Show file tree

Hide file tree

Showing 4 changed files with 26 additions and 11 deletions.
diff --git a/include/perturbations.h b/include/perturbations.h
@@ -181,6 +181,8 @@ struct perturbs
   double three_ceff2_ur;/**< 3 x effective squared sound speed for the ultrarelativistic perturbations */
   double three_cvis2_ur;/**< 3 x effective viscosity parameter for the ultrarelativistic perturbations */
 
+  double z_max_pk; /**< when we compute only the matter spectrum / transfer functions, but not the CMB, we are sometimes interested to sample source functions at very high redshift, way before recombination. This z_max_pk will then fix the initial sampling time of the sources. */
+
   //@}
 
   /** @name - useful flags infered from the ones above */

diff --git a/source/input.c b/source/input.c
@@ -2098,13 +2098,15 @@ int input_read_parameters(
                errmsg);
 
     if (flag1==_TRUE_) {
-      psp->z_max_pk = param1;
+      ppt->z_max_pk = param1;
     }
     else {
-      psp->z_max_pk = 0.;
-      for (i=0; i<pop->z_pk_num; i++)
-        psp->z_max_pk = MAX(psp->z_max_pk,pop->z_pk[i]);
+      ppt->z_max_pk = 0.;
+      for (i=0; i<pop->z_pk_num; i++) {
+        ppt->z_max_pk = MAX(ppt->z_max_pk,pop->z_pk[i]);
+      }
     }
+    psp->z_max_pk = ppt->z_max_pk;
   }
 
   /* deal with selection functions */
@@ -2889,6 +2891,8 @@ int input_default_params(
   ppt->three_ceff2_ur=1.;
   ppt->three_cvis2_ur=1.;
 
+  ppt->z_max_pk=0.;
+
   ppt->selection_num=1;
   ppt->selection=gaussian;
   ppt->selection_mean[0]=1.;

diff --git a/source/perturbations.c b/source/perturbations.c
@@ -960,8 +960,15 @@ int perturb_timesampling_for_sources(
   }
   else {
 
-    /* case when CMB not requested: start at recombination time */
-    tau_ini = pth->tau_rec;
+    /* check the time corresponding to the highest redshift requested in output plus one */
+    class_call(background_tau_of_z(pba,
+                                   ppt->z_max_pk+1,
+                                   &tau_ini),
+               pba->error_message,
+               ppt->error_message);
+
+    /* obsolete: previous choice was to start always at recombination time */
+    /* tau_ini = pth->tau_rec; */
 
     /* set values of first_index_back/thermo */
     class_call(background_at_tau(pba,
@@ -984,9 +991,6 @@ int perturb_timesampling_for_sources(
                ppt->error_message);
   }
 
-
-  counter = 1;
-
   /** (b) next sampling point = previous + ppr->perturb_sampling_stepsize * timescale_source, where:
       - if CMB requested:
       timescale_source1 = \f$ |g/\dot{g}| = |\dot{\kappa}-\ddot{\kappa}/\dot{\kappa}|^{-1} \f$;
@@ -996,6 +1000,7 @@ int perturb_timesampling_for_sources(
       timescale_source = 1/aH; repeat till today.
   */
 
+  counter = 1;
   last_index_back = first_index_back;
   last_index_thermo = first_index_thermo;
   tau = tau_ini;

diff --git a/source/spectra.c b/source/spectra.c
@@ -2471,14 +2471,18 @@ int spectra_k_and_tau(
                psp->error_message);
 
     index_tau=0;
-    class_test((tau_min < ppt->tau_sampling[index_tau]),
+    class_test((tau_min <= ppt->tau_sampling[index_tau]),
                psp->error_message,
-               "you asked for zmax=%e, i.e. taumin=%e, smaller than first possible value =%e",psp->z_max_pk,tau_min,ppt->tau_sampling[0]);
+               "you asked for zmax=%e, i.e. taumin=%e, smaller than or equal to the first possible value =%e; it should be strictly bigger for a successfull interpolation",psp->z_max_pk,tau_min,ppt->tau_sampling[0]);
 
     while (ppt->tau_sampling[index_tau] < tau_min){
       index_tau++;
     }
     index_tau --;
+    class_test(index_tau<0,
+               psp->error_message,
+               "by construction, this should never happen, a bug must have been introduced somewhere");
+
     /* whenever possible, take a few more values in to avoid boundary effects in the interpolation */
     if (index_tau>0) index_tau--;
     if (index_tau>0) index_tau--;