Merge pull request #374 from LSSTDESC/tests_for_hiz_curvature_neutrinos

Distance tests to z=1050 for massive neutrinos and curvature

include/ccl_constants.h

@@ -87,6 +87,12 @@ extern "C" {
  */
 #define EV_IN_J  GSL_CONST_MKSA_ELECTRON_VOLT
 
+/**
+ * Temperature of the CMB in K
+ */
+#define TCMB 2.725
+//#define TCMB 2.7255 // CLASS value
+
 /**
  * T_ncdm, as taken from CLASS, explanatory.ini
  */

include/ccl_neutrinos.h

@@ -43,27 +43,27 @@ gsl_spline* calculate_nu_phasespace_spline(int *status);
  * @param a Scale factor
  * @param Neff The effective number of species with neutrino mass mnu.
  * @param mnu Pointer to array containing neutrino mass (can be 0).
- * @param TCMB Temperature of the CMB
+ * @param T_CMB Temperature of the CMB
  * @param accel - Interpolation accelerator to be used with phasespace spline. If not set yet, pass NULL.
  * @param status Status flag. 0 if there are no errors, nonzero otherwise.
  * For specific cases see documentation for ccl_error.c
  * @return OmNuh2 Fractional energy density of neutrions with mass mnu, multiplied by h squared. 
  */
 
-double ccl_Omeganuh2 (double a, int N_nu_mass, double* mnu, double TCMB, gsl_interp_accel* accel, int * status);
+double ccl_Omeganuh2 (double a, int N_nu_mass, double* mnu, double T_CMB, gsl_interp_accel* accel, int * status);
 
 /** 
  * Returns mass of one neutrino species at a scale factor a. 
  * @param a Scale factor
  * @param Neff The effective number of species with neutrino mass mnu.
  * @param OmNuh2 Fractional energy density of neutrions with mass mnu, multiplied by h squared. (can be 0).
- * @param TCMB Temperature of the CMB
+ * @param T_CMB Temperature of the CMB
  * @param accel - Interpolation accelerator to be used with phasespace spline. If not set yet, pass NULL.
  * @param status Status flag. 0 if there are no errors, nonzero otherwise.
  * For specific cases see documentation for ccl_error.c
  * @return Mnu Neutrino mass [eV]. 
  */
-double* ccl_nu_masses (double OmNuh2, ccl_neutrino_mass_splits mass_split, double TCMB, int * status);
+double* ccl_nu_masses (double OmNuh2, ccl_neutrino_mass_splits mass_split, double T_CMB, int * status);
 
 #ifdef __cplusplus
 }

pyccl/ccl_neutrinos.i

@@ -19,7 +19,7 @@
 
 %inline %{
 
-void Omeganuh2_vec(int N_nu_mass, double TCMB,
+void Omeganuh2_vec(int N_nu_mass, double T_CMB,
                    double* a, int na, 
                    double* mnu, int nm, 
                    double* output, int nout,
@@ -28,16 +28,16 @@ void Omeganuh2_vec(int N_nu_mass, double TCMB,
     assert(nout == na);
     assert(nm == 3);
     for(int i=0; i < na; i++){
-        output[i] = ccl_Omeganuh2(a[i], N_nu_mass, mnu, TCMB, NULL, status);
+        output[i] = ccl_Omeganuh2(a[i], N_nu_mass, mnu, T_CMB, NULL, status);
     }   
 }
 
-void nu_masses_vec(double OmNuh2, int label, double TCMB,
+void nu_masses_vec(double OmNuh2, int label, double T_CMB,
                           double* output, int nout,
                           int* status)
 {
     double* mnu;
-    mnu = ccl_nu_masses(OmNuh2, label, TCMB, status);
+    mnu = ccl_nu_masses(OmNuh2, label, T_CMB, status);
     for(int i=0; i < nout; i++){
         output[i] = *(mnu+i);
     }

pyccl/ccl_wrap.c

@@ -4370,7 +4370,7 @@ user_pz_info* specs_create_photoz_info_from_py(PyObject *pyfunc)
 
 
 
-void Omeganuh2_vec(int N_nu_mass, double TCMB,
+void Omeganuh2_vec(int N_nu_mass, double T_CMB,
                    double* a, int na, 
                    double* mnu, int nm, 
                    double* output, int nout,
@@ -4379,16 +4379,16 @@ void Omeganuh2_vec(int N_nu_mass, double TCMB,
     assert(nout == na);
     assert(nm == 3);
     for(int i=0; i < na; i++){
-        output[i] = ccl_Omeganuh2(a[i], N_nu_mass, mnu, TCMB, NULL, status);
+        output[i] = ccl_Omeganuh2(a[i], N_nu_mass, mnu, T_CMB, NULL, status);
     }   
 }
 
-void nu_masses_vec(double OmNuh2, int label, double TCMB,
+void nu_masses_vec(double OmNuh2, int label, double T_CMB,
                           double* output, int nout,
                           int* status)
 {
     double* mnu;
-    mnu = ccl_nu_masses(OmNuh2, label, TCMB, status);
+    mnu = ccl_nu_masses(OmNuh2, label, T_CMB, status);
     for(int i=0; i < nout; i++){
         output[i] = *(mnu+i);
     }
@@ -24008,10 +24008,10 @@ static PyMethodDef SwigMethods[] = {
 	 { (char *)"call_py_photoz_fn", _wrap_call_py_photoz_fn, METH_VARARGS, (char *)"call_py_photoz_fn(double z_ph, double z_s, void * py_func_obj, int * status) -> double"},
 	 { (char *)"specs_create_photoz_info_from_py", _wrap_specs_create_photoz_info_from_py, METH_VARARGS, (char *)"specs_create_photoz_info_from_py(PyObject * pyfunc) -> user_pz_info"},
 	 { (char *)"calculate_nu_phasespace_spline", _wrap_calculate_nu_phasespace_spline, METH_VARARGS, (char *)"calculate_nu_phasespace_spline(int * status) -> gsl_spline *"},
-	 { (char *)"Omeganuh2", _wrap_Omeganuh2, METH_VARARGS, (char *)"Omeganuh2(double a, int N_nu_mass, double * mnu, double TCMB, gsl_interp_accel * accel, int * status) -> double"},
-	 { (char *)"nu_masses", _wrap_nu_masses, METH_VARARGS, (char *)"nu_masses(double OmNuh2, ccl_neutrino_mass_splits mass_split, double TCMB, int * status) -> double *"},
-	 { (char *)"Omeganuh2_vec", _wrap_Omeganuh2_vec, METH_VARARGS, (char *)"Omeganuh2_vec(int N_nu_mass, double TCMB, double * a, double * mnu, double * output, int * status)"},
-	 { (char *)"nu_masses_vec", _wrap_nu_masses_vec, METH_VARARGS, (char *)"nu_masses_vec(double OmNuh2, int label, double TCMB, double * output, int * status)"},
+	 { (char *)"Omeganuh2", _wrap_Omeganuh2, METH_VARARGS, (char *)"Omeganuh2(double a, int N_nu_mass, double * mnu, double T_CMB, gsl_interp_accel * accel, int * status) -> double"},
+	 { (char *)"nu_masses", _wrap_nu_masses, METH_VARARGS, (char *)"nu_masses(double OmNuh2, ccl_neutrino_mass_splits mass_split, double T_CMB, int * status) -> double *"},
+	 { (char *)"Omeganuh2_vec", _wrap_Omeganuh2_vec, METH_VARARGS, (char *)"Omeganuh2_vec(int N_nu_mass, double T_CMB, double * a, double * mnu, double * output, int * status)"},
+	 { (char *)"nu_masses_vec", _wrap_nu_masses_vec, METH_VARARGS, (char *)"nu_masses_vec(double OmNuh2, int label, double T_CMB, double * output, int * status)"},
 	 { (char *)"spline_params_A_SPLINE_NA_set", _wrap_spline_params_A_SPLINE_NA_set, METH_VARARGS, (char *)"spline_params_A_SPLINE_NA_set(spline_params self, int A_SPLINE_NA)"},
 	 { (char *)"spline_params_A_SPLINE_NA_get", _wrap_spline_params_A_SPLINE_NA_get, METH_VARARGS, (char *)"spline_params_A_SPLINE_NA_get(spline_params self) -> int"},
 	 { (char *)"spline_params_A_SPLINE_MIN_set", _wrap_spline_params_A_SPLINE_MIN_set, METH_VARARGS, (char *)"spline_params_A_SPLINE_MIN_set(spline_params self, double A_SPLINE_MIN)"},
@@ -25092,6 +25092,7 @@ SWIG_init(void) {
   SWIG_Python_SetConstant(d, "HPLANCK",SWIG_From_double((double)((6.62606896e-34))));
   SWIG_Python_SetConstant(d, "CLIGHT",SWIG_From_double((double)((2.99792458e8))));
   SWIG_Python_SetConstant(d, "EV_IN_J",SWIG_From_double((double)((1.602176487e-19))));
+  SWIG_Python_SetConstant(d, "TCMB",SWIG_From_double((double)(2.725)));
   SWIG_Python_SetConstant(d, "TNCDM",SWIG_From_double((double)(0.71611)));
   SWIG_Python_SetConstant(d, "DELTAM12_sq",SWIG_From_double((double)(7.62E-5)));
   SWIG_Python_SetConstant(d, "DELTAM13_sq_pos",SWIG_From_double((double)(2.55E-3)));

pyccl/ccllib.py

@@ -1096,6 +1096,7 @@ def clt_fa_vec(cosmo, clt, func_code, aarr, output, status):
 HPLANCK = _ccllib.HPLANCK
 CLIGHT = _ccllib.CLIGHT
 EV_IN_J = _ccllib.EV_IN_J
+TCMB = _ccllib.TCMB
 TNCDM = _ccllib.TNCDM
 DELTAM12_sq = _ccllib.DELTAM12_sq
 DELTAM13_sq_pos = _ccllib.DELTAM13_sq_pos
@@ -1214,21 +1215,21 @@ def calculate_nu_phasespace_spline(status):
     """calculate_nu_phasespace_spline(int * status) -> gsl_spline *"""
     return _ccllib.calculate_nu_phasespace_spline(status)
 
-def Omeganuh2(a, N_nu_mass, mnu, TCMB, accel, status):
-    """Omeganuh2(double a, int N_nu_mass, double * mnu, double TCMB, gsl_interp_accel * accel, int * status) -> double"""
-    return _ccllib.Omeganuh2(a, N_nu_mass, mnu, TCMB, accel, status)
+def Omeganuh2(a, N_nu_mass, mnu, T_CMB, accel, status):
+    """Omeganuh2(double a, int N_nu_mass, double * mnu, double T_CMB, gsl_interp_accel * accel, int * status) -> double"""
+    return _ccllib.Omeganuh2(a, N_nu_mass, mnu, T_CMB, accel, status)
 
-def nu_masses(OmNuh2, mass_split, TCMB, status):
-    """nu_masses(double OmNuh2, ccl_neutrino_mass_splits mass_split, double TCMB, int * status) -> double *"""
-    return _ccllib.nu_masses(OmNuh2, mass_split, TCMB, status)
+def nu_masses(OmNuh2, mass_split, T_CMB, status):
+    """nu_masses(double OmNuh2, ccl_neutrino_mass_splits mass_split, double T_CMB, int * status) -> double *"""
+    return _ccllib.nu_masses(OmNuh2, mass_split, T_CMB, status)
 
-def Omeganuh2_vec(N_nu_mass, TCMB, a, mnu, output, status):
-    """Omeganuh2_vec(int N_nu_mass, double TCMB, double * a, double * mnu, double * output, int * status)"""
-    return _ccllib.Omeganuh2_vec(N_nu_mass, TCMB, a, mnu, output, status)
+def Omeganuh2_vec(N_nu_mass, T_CMB, a, mnu, output, status):
+    """Omeganuh2_vec(int N_nu_mass, double T_CMB, double * a, double * mnu, double * output, int * status)"""
+    return _ccllib.Omeganuh2_vec(N_nu_mass, T_CMB, a, mnu, output, status)
 
-def nu_masses_vec(OmNuh2, label, TCMB, output, status):
-    """nu_masses_vec(double OmNuh2, int label, double TCMB, double * output, int * status)"""
-    return _ccllib.nu_masses_vec(OmNuh2, label, TCMB, output, status)
+def nu_masses_vec(OmNuh2, label, T_CMB, output, status):
+    """nu_masses_vec(double OmNuh2, int label, double T_CMB, double * output, int * status)"""
+    return _ccllib.nu_masses_vec(OmNuh2, label, T_CMB, output, status)
 class spline_params(_object):
     """Proxy of C ccl_spline_params struct."""
 

src/ccl_core.c

@@ -279,16 +279,22 @@ void ccl_parameters_fill_initial(ccl_parameters * params, int *status)
 {
   // Fixed radiation parameters
   // Omega_g * h**2 is known from T_CMB
-  params->T_CMB =  2.725; 
-  params->Omega_g = 4. * STBOLTZ / CLIGHT *pow(params->T_CMB,4.)/(3. * pow(10., 10.) * CLIGHT * CLIGHT *params->h* params->h / (8. * M_PI * GNEWT * MPC_TO_METER * MPC_TO_METER));
+  params->T_CMB =  TCMB;
+  // kg / m^3
+  double rho_g = 4. * STBOLTZ / pow(CLIGHT, 3) * pow(params->T_CMB, 4);
+  // kg / m^3
+  double rho_crit = RHO_CRITICAL * SOLAR_MASS/pow(MPC_TO_METER, 3) * pow(params->h, 2);
+  params->Omega_g = rho_g/rho_crit;
 
   // Get the N_nu_rel from Neff and N_nu_mass
-  params->N_nu_rel = params->Neff - params->N_nu_mass * TNCDM * TNCDM * TNCDM * TNCDM / pow(4./11.,4./3.);
-
-  //Get the relativistic neutrino Omega_nu. It's more efficient to get this in-line, to avoid computing the phase-space integral if not necessary.
-  double Tnu= (params->T_CMB) *pow(4./11.,1./3.); 
-  params-> Omega_n_rel = params->N_nu_rel* 8. * pow(M_PI,5) *pow((KBOLTZ/ HPLANCK),3)* KBOLTZ/(15. *pow( CLIGHT,3))* (8. * M_PI * GNEWT) / (3. * 100.*100.*1000.*1000. /MPC_TO_METER /MPC_TO_METER  * CLIGHT * CLIGHT)  * Tnu * Tnu * Tnu * Tnu *7./8.;
+  params->N_nu_rel = params->Neff - params->N_nu_mass * pow(TNCDM, 4) / pow(4./11.,4./3.);
 
+  // Temperature of the relativistic neutrinos in K
+  double T_nu= (params->T_CMB) * pow(4./11.,1./3.); 
+  // in kg / m^3
+  double rho_nu_rel = params->N_nu_rel* 7.0/8.0 * 4. * STBOLTZ / pow(CLIGHT, 3) * pow(T_nu, 4);
+  params-> Omega_n_rel = rho_nu_rel/rho_crit;
+
   // If non-relativistic neutrinos are present, calculate the phase_space integral.
   if((params->N_nu_mass)>0) {
     // Pass NULL for the accelerator here because we don't have our cosmology object defined yet.

src/ccl_neutrinos.c

@@ -102,12 +102,12 @@ double nu_phasespace_intg(gsl_interp_accel* accel, double mnuOT, int* status)
 }
 
 /* -------- ROUTINE: Omeganuh2 ---------
-INPUTS: a: scale factor, Nnumass: number of massive neutrino species, mnu: total mass in eV of neutrinos, TCMB: CMB temperature, accel: pointer to an accelerator which will evaluate the neutrino phasespace spline if defined, status: pointer to status integer.
+INPUTS: a: scale factor, Nnumass: number of massive neutrino species, mnu: total mass in eV of neutrinos, T_CMB: CMB temperature, accel: pointer to an accelerator which will evaluate the neutrino phasespace spline if defined, status: pointer to status integer.
 TASK: Compute Omeganu * h^2 as a function of time.
 !! To all practical purposes, Neff is simply N_nu_mass !!
 */
 
-double ccl_Omeganuh2 (double a, int N_nu_mass, double* mnu, double TCMB, gsl_interp_accel* accel, int* status)
+double ccl_Omeganuh2 (double a, int N_nu_mass, double* mnu, double T_CMB, gsl_interp_accel* accel, int* status)
 {
   double Tnu, a4, prefix_massless, mnuone, OmNuh2;
   double Tnu_eff, mnuOT, intval, prefix_massive;
@@ -116,7 +116,7 @@ double ccl_Omeganuh2 (double a, int N_nu_mass, double* mnu, double TCMB, gsl_int
   // First check if N_nu_mass is 0
   if (N_nu_mass==0) return 0.0;  
 
-  Tnu=TCMB*pow(4./11.,1./3.);
+  Tnu=T_CMB*pow(4./11.,1./3.);
   a4=a*a*a*a;  
   // Check if mnu=0. We assume that in the massless case mnu is a pointer to a single element and that element is 0. This should in principle never be called.
   if (mnu[0] < 0.00017) {  // Limit taken from Lesgourges et al. 2012
@@ -125,7 +125,7 @@ double ccl_Omeganuh2 (double a, int N_nu_mass, double* mnu, double TCMB, gsl_int
   }
 
   // And the remaining massive case.
-  // Tnu_eff is used in the massive case because CLASS uses an effective temperature of nonLCDM components to match to mnu / Omeganu =93.14eV. Tnu_eff = T_ncdm * TCMB = 0.71611 * TCMB
+  // Tnu_eff is used in the massive case because CLASS uses an effective temperature of nonLCDM components to match to mnu / Omeganu =93.14eV. Tnu_eff = T_ncdm * T_CMB = 0.71611 * T_CMB
   Tnu_eff = Tnu * TNCDM / (pow(4./11.,1./3.));
 
   // Define the prefix using the effective temperature (to get mnu / Omega = 93.14 eV) for the massive case: 
@@ -146,11 +146,11 @@ double ccl_Omeganuh2 (double a, int N_nu_mass, double* mnu, double TCMB, gsl_int
 }
 
 /* -------- ROUTINE: Omeganuh2_to_Mnu ---------
-INPUTS: OmNuh2: neutrino mass density today Omeganu * h^2, label: how you want to split up the masses, see ccl_neutrinos.h for options, TCMB: CMB temperature, accel: pointer to an accelerator which will evaluate the neutrino phasespace spline if defined, status: pointer to status integer.
+INPUTS: OmNuh2: neutrino mass density today Omeganu * h^2, label: how you want to split up the masses, see ccl_neutrinos.h for options, T_CMB: CMB temperature, accel: pointer to an accelerator which will evaluate the neutrino phasespace spline if defined, status: pointer to status integer.
 TASK: Given Omeganuh2 today, the method of splitting into masses, and the temperature of the CMB, output a pointer to the array of neutrino masses (may be length 1 if label asks for sum) 
 */
 
-double* ccl_nu_masses(double OmNuh2, ccl_neutrino_mass_splits mass_split, double TCMB,  int* status){
+double* ccl_nu_masses(double OmNuh2, ccl_neutrino_mass_splits mass_split, double T_CMB,  int* status){
 
   double sumnu;
 

tests/benchmark/chi_hiz_mnu_model6-15.txt

@@ -0,0 +1,11 @@
+# z flat_nonu   pos_curv_nonu   neg_curv_nonu   flat_massnu1   flat_massnu2   flat_massnu3   flat_manynu1   neg_curv_massnu1   pos_curv_manynu1
+1.000000000000000021e-02 4.273088716101704421e+01 4.272875822309976002e+01 4.273301652297963926e+01 4.273019036037864282e+01 4.272977903984298109e+01 4.272963874233249726e+01 4.273087473464769914e+01 4.273210147576675411e+01 4.272749759389967039e+01
+3.613349096754857326e-02 1.534826496361100681e+02 1.534552824843099472e+02 1.535100362974301902e+02 1.534736134787392530e+02 1.534682809803454973e+02 1.534664618499638209e+02 1.534824592225216122e+02 1.534981349729091278e+02 1.534389160092656255e+02
+1.305629169501913711e-01 5.423419259913349606e+02 5.420055639531829001e+02 5.426791199163794772e+02 5.422273521635887619e+02 5.421598072211883164e+02 5.421367535132249031e+02 5.423394049762028999e+02 5.425278758343570189e+02 5.417983884993147967e+02
+4.717693980316532421e-01 1.795382591281910209e+03 1.791986203452796190e+03 1.798804069805638619e+03 1.794099688430913375e+03 1.793346077384032469e+03 1.793088355224907673e+03 1.795350962528290438e+03 1.797099959167988800e+03 1.789677763086143614e+03
+1.704667528254257158e+00 4.721650617302094361e+03 4.704914319325481301e+03 4.738598375794766071e+03 4.713489591559931796e+03 4.708742006344820766e+03 4.707105423165194225e+03 4.721375013717224647e+03 4.727646762947467323e+03 4.690286133803381745e+03
+6.159558873484855646e+00 8.310099119530776079e+03 8.280512195286313727e+03 8.340043152825259313e+03 8.289938274867332439e+03 8.278420661769756407e+03 8.274362151373314191e+03 8.308949941748644051e+03 8.312664834354787672e+03 8.244066805131278670e+03
+2.225663649191484694e+01 1.090450810236021789e+04 1.087123877020451800e+04 1.093815187992903520e+04 1.087558781951848869e+04 1.085955215624802986e+04 1.085368072035978184e+04 1.090145479836753839e+04 1.089781728832543922e+04 1.081790293080846823e+04
+8.042099736486171935e+01 1.240496559920679465e+04 1.237105152135618118e+04 1.243925505318421710e+04 1.237168347053333491e+04 1.235385639566074133e+04 1.234706847566714350e+04 1.239824294137528341e+04 1.239001876020283999e+04 1.230703967922422999e+04
+2.905891381884478619e+02 1.320509096754415441e+04 1.317108158904521588e+04 1.323947576060097344e+04 1.317028730553663263e+04 1.315198598247911832e+04 1.314490590675488784e+04 1.319175249207913657e+04 1.318190256301443333e+04 1.309822126050697807e+04
+1.050000000000000455e+03 1.360851790987604181e+04 1.357449619903580970e+04 1.364291503657256180e+04 1.357336254204349098e+04 1.355490335352935472e+04 1.354777051608289912e+04 1.358483929801837439e+04 1.357460061898953245e+04 1.349077421513099762e+04