Adding getTkm method to ClassEngine

cpp/ClassEngine.cc

@@ -465,6 +465,84 @@ ClassEngine::getTk( double z,
   }
 }
 
+void 
+ClassEngine::getTkm( double z, 
+   std::vector<double>& k,
+   std::vector<double>& d_cdm,
+   std::vector<double>& d_b,
+   std::vector<double>& d_ncdm,
+   std::vector<double>& d_m,
+   std::vector<double>& t_cdm,
+   std::vector<double>& t_b,
+   std::vector<double>& t_ncdm,
+   std::vector<double>& t_m )
+{
+
+  if (!dofree) throw out_of_range("no Tk available because CLASS failed");
+
+  double tau;
+  int index;
+  //transform redshift in conformal time
+  background_tau_of_z(&ba,z,&tau);
+
+  if(log(tau) < pt.ln_tau[0]){
+    cerr << "Asking sources at a z bigger than z_max_pk, something probably went wrong\n";
+    throw out_of_range(pt.error_message);
+  }
+
+  double *pvecback=new double[ba.bg_size];
+  background_at_tau(&ba,tau,ba.long_info,ba.inter_normal, &index, pvecback);
+  double fHa = pvecback[ba.index_bg_f] * (pvecback[ba.index_bg_a]*pvecback[ba.index_bg_H]);
+  delete[] pvecback;
+
+  // copy transfer func data to temporary
+  const size_t index_md = pt.index_md_scalars;
+  d_cdm.assign( pt.k_size[index_md], 0.0 );
+  d_b.assign( pt.k_size[index_md], 0.0 );
+  d_ncdm.assign( pt.k_size[index_md], 0.0 );
+  d_m.assign( pt.k_size[index_md], 0.0 );
+  t_cdm.assign( pt.k_size[index_md], 0.0 );
+  t_b.assign( pt.k_size[index_md], 0.0 );
+  t_ncdm.assign( pt.k_size[index_md], 0.0 );
+  t_m.assign( pt.k_size[index_md], 0.0 );
+
+  if( pt.ic_size[index_md] > 1 ){
+    cerr << ">>>have more than 1 ICs, will use first and ignore others" << endl;
+  }
+
+  call_perturb_sources_at_tau(index_md, 0, pt.index_tp_delta_cdm, tau, &d_cdm[0]);
+  call_perturb_sources_at_tau(index_md, 0, pt.index_tp_delta_b, tau, &d_b[0]);
+  call_perturb_sources_at_tau(index_md, 0, pt.index_tp_delta_ncdm1, tau, &d_ncdm[0]);
+  call_perturb_sources_at_tau(index_md, 0, pt.index_tp_delta_m, tau, &d_m[0]);
+  call_perturb_sources_at_tau(index_md, 0, pt.index_tp_theta_b, tau, &t_b[0]);
+  call_perturb_sources_at_tau(index_md, 0, pt.index_tp_theta_ncdm1, tau, &t_ncdm[0]);
+  call_perturb_sources_at_tau(index_md, 0, pt.index_tp_theta_m, tau, &t_m[0]);
+
+  //
+  std::vector<double> h_prime(pt.k_size[index_md],0.0), eta_prime(pt.k_size[index_md],0.0);
+  call_perturb_sources_at_tau(index_md, 0, pt.index_tp_eta_prime, tau, &eta_prime[0]);
+  call_perturb_sources_at_tau(index_md, 0, pt.index_tp_h_prime, tau, &h_prime[0]);
+
+  // gauge trafo velocities, store k-vector
+  for (int index_k=0; index_k<pt.k_size[index_md]; index_k++) 
+  {
+    auto ak = pt.k[index_md][index_k];
+
+    // write data to vectors
+    k.push_back( ak );
+
+    // use the conformal Newtonian gauge for velocities
+    // not correct, but N-body gauge currently not implemented
+    double alphak2 = (h_prime[index_k]+6*eta_prime[index_k])/2; 
+
+    t_cdm[index_k]  = (-alphak2) / fHa;
+    t_b[index_k]    = (-alphak2 + t_b[index_k]) / fHa;
+    t_ncdm[index_k] = (-alphak2 + t_ncdm[index_k]) / fHa;
+    t_m[index_k]  = (-alphak2 + t_m[index_k]) / fHa;
+  }
+}
+
+
 double
 ClassEngine::getCl(Engine::cltype t,const long &l){
 

cpp/ClassEngine.hh

@@ -114,6 +114,16 @@ public:
         std::vector<double>& t_b,
         std::vector<double>& t_ncdm,
         std::vector<double>& t_tot );
+  void getTkm( double z, 
+        std::vector<double>& k,
+        std::vector<double>& d_cdm,
+        std::vector<double>& d_b,
+        std::vector<double>& d_ncdm,
+        std::vector<double>& d_m,
+        std::vector<double>& t_cdm,
+        std::vector<double>& t_b,
+        std::vector<double>& t_ncdm,
+        std::vector<double>& t_m );
 
  //for BAO
   inline double z_drag() const {return th.z_d;}