added support for velocity transfer functions added with the latest (early 2015) CAMB release

plugins/transfer_camb.cc

@@ -1,8 +1,7 @@
 /*
  
  transfer_camb.cc - This file is part of MUSIC -
- a code to generate multi-scale initial conditions 
+ a code to generate multi-scale initial conditions for cosmological simulations 
- for cosmological simulations 
  
  Copyright (C) 2010  Oliver Hahn
  
@@ -10,25 +9,32 @@
 
 #include "transfer_function.hh"
 
+const double tiny = 1e-30;
+
 class transfer_CAMB_plugin : public transfer_function_plugin
 {
   
 private:
   std::string m_filename_Pk, m_filename_Tk;
   std::vector<double> m_tab_k, m_tab_Tk_tot, m_tab_Tk_cdm, m_tab_Tk_baryon;
+  std::vector<double> m_tab_Tvk_cdm, m_tab_Tvk_baryon; //>[150609SH: add]
   gsl_interp_accel *acc_tot, *acc_cdm, *acc_baryon;
+  gsl_interp_accel *acc_vcdm, *acc_vbaryon; //>[150609SH: add]
   gsl_spline *spline_tot, *spline_cdm, *spline_baryon;
+  gsl_spline *spline_vcdm, *spline_vbaryon; //>[150609SH: add]
   
   double m_kmin, m_kmax;
   unsigned m_nlines;
 
+  bool m_linbaryoninterp;
+
+  double m_velunits;
+  
   void read_table( void ){
 #ifdef WITH_MPI
     if( MPI::COMM_WORLD.Get_rank() == 0 ){
 #endif 
-			std::cerr 
+      LOGINFO("Reading tabulated transfer function data from file \n    \'%s\'",m_filename_Tk.c_str() );
-			<< " - reading tabulated transfer function data from file \n"
-			<< "    \'" << m_filename_Tk << "\'\n";
       
       std::string line;
       std::ifstream ifs( m_filename_Tk.c_str() );
@@ -40,10 +46,13 @@ private:
       m_tab_Tk_tot.clear();
       m_tab_Tk_cdm.clear();
       m_tab_Tk_baryon.clear();
+      m_tab_Tvk_cdm.clear(); //>[150609SH: add]
+      m_tab_Tvk_baryon.clear(); //>[150609SH: add]
   
+
+      //>[150609SH: add]
       m_kmin = 1e30;
       m_kmax = -1e30;
-            m_nlines = 0;
 
       while( !ifs.eof() ){
 	getline(ifs,line);
@@ -52,30 +61,56 @@ private:
 	
 	std::stringstream ss(line);
 	
-				double k, Tkc, Tkb, Tkg, Tkr, Tknu, Tktot;
+	double k, Tkc, Tkb, Tktot, Tkvc, Tkvb, dummy;
+
 	ss >> k;
 	ss >> Tkc;
 	ss >> Tkb;
-				ss >> Tkg;
+	ss >> dummy;
-				ss >> Tkr;
+	ss >> dummy;
-				ss >> Tknu;
+	ss >> dummy;
 	ss >> Tktot;
+	ss >> dummy;
+	ss >> dummy;
+	ss >> dummy;
+	ss >> Tkvc; //>[150609SH: add]
+	ss >> Tkvb; //>[150609SH: add]
+	ss >> dummy;
 
-                if( k < m_kmin ) m_kmin = k;
+	m_linbaryoninterp |= Tkb < 0.0 || Tkvb < 0.0;
-                if( k > m_kmax ) m_kmax = k;
 	
 	m_tab_k.push_back( log10(k) );
 	
-				m_tab_Tk_tot.push_back( log10(Tktot) );
+	m_tab_Tk_tot.push_back( Tktot );
-				m_tab_Tk_baryon.push_back( log10(Tkb) );
+	m_tab_Tk_baryon.push_back( Tkb );
-				m_tab_Tk_cdm.push_back( log10(Tkc) );
+	m_tab_Tk_cdm.push_back( Tkc );
+	m_tab_Tvk_baryon.push_back( Tkvb );
+	m_tab_Tvk_cdm.push_back( Tkvc );
+
 	++m_nlines;
+
+	if( k < m_kmin ) m_kmin = k;
+	if( k > m_kmax ) m_kmax = k;
+      }
+      
+      for( size_t i=0; i<m_tab_k.size(); ++i ){
+	m_tab_Tk_tot[i]     = log10( m_tab_Tk_tot[i] );
+	m_tab_Tk_cdm[i]     = log10( m_tab_Tk_cdm[i] );
+	m_tab_Tvk_cdm[i]    = log10( m_tab_Tvk_cdm[i] );
+
+	if( !m_linbaryoninterp )
+	  {
+	    m_tab_Tk_baryon[i]  = log10( m_tab_Tk_baryon[i] );
+	    m_tab_Tvk_baryon[i] = log10( m_tab_Tvk_baryon[i] );
+	  }
 	}
       
       ifs.close();
 
+      LOGINFO("Read CAMB transfer function table with %d rows", m_nlines);
 
-			
+      if( m_linbaryoninterp )
+	LOGINFO("Using log-lin interpolation for baryons\n    (TF is not positive definite)");
       
 #ifdef WITH_MPI
     }
@@ -88,13 +123,16 @@ private:
       m_tab_Tk_tot.assign(n,0);
       m_tab_Tk_cdm.assign(n,0);
       m_tab_Tk_baryon.assign(n,0);
-
+      m_tab_Tvk_cdm.assign(n,0); //>[150609SH: add]
+      m_tab_Tvk_baryon.assign(n,0); //>[150609SH: add]
     }
     
     MPI::COMM_WORLD.Bcast( &m_tab_k[0],  n, MPI_DOUBLE, 0 );
     MPI::COMM_WORLD.Bcast( &m_tab_Tk_tot[0], n, MPI_DOUBLE, 0 );
     MPI::COMM_WORLD.Bcast( &m_tab_Tk_cdm[0], n, MPI_DOUBLE, 0 );
     MPI::COMM_WORLD.Bcast( &m_tab_Tk_baryon[0], n, MPI_DOUBLE, 0 );
+    MPI::COMM_WORLD.Bcast( &m_tab_Tvk_cdm[0], n, MPI_DOUBLE, 0 ); //>[150609SH: add]
+    MPI::COMM_WORLD.Bcast( &m_tab_Tvk_baryon[0], n, MPI_DOUBLE, 0 ); //>[150609SH: add]
     
 #endif
 		
@@ -111,18 +149,23 @@ public:
     acc_tot     = gsl_interp_accel_alloc();
     acc_cdm     = gsl_interp_accel_alloc();
     acc_baryon  = gsl_interp_accel_alloc();
-		
+    acc_vcdm    = gsl_interp_accel_alloc(); //>[150609SH: add]
+    acc_vbaryon = gsl_interp_accel_alloc(); //>[150609SH: add]
         
     spline_tot     = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() );
     spline_cdm     = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() );
     spline_baryon  = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() );
+    spline_vcdm    = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() ); //>[150609SH: add]
+    spline_vbaryon = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() ); //>[150609SH: add]
     
     gsl_spline_init (spline_tot,     &m_tab_k[0], &m_tab_Tk_tot[0],     m_tab_k.size() );
     gsl_spline_init (spline_cdm,     &m_tab_k[0], &m_tab_Tk_cdm[0],     m_tab_k.size() );
     gsl_spline_init (spline_baryon,  &m_tab_k[0], &m_tab_Tk_baryon[0],  m_tab_k.size() );
+    gsl_spline_init (spline_vcdm,    &m_tab_k[0], &m_tab_Tvk_cdm[0],    m_tab_k.size() ); //>[150609SH: add]
+    gsl_spline_init (spline_vbaryon, &m_tab_k[0], &m_tab_Tvk_baryon[0], m_tab_k.size() ); //>[150609SH: add]
     
-		tf_distinct_ = true;
+    tf_distinct_ = true;   // [150612SH: different density between CDM v.s. Baryon]
-		tf_withvel_  = false;
+    tf_withvel_  = true;   // [150612SH: using velocity transfer function]
   }
   
   ~transfer_CAMB_plugin()
@@ -130,59 +173,82 @@ public:
     gsl_spline_free (spline_tot);
     gsl_spline_free (spline_cdm);
     gsl_spline_free (spline_baryon);
+    gsl_spline_free (spline_vcdm); //>[150609SH: add]
+    gsl_spline_free (spline_vbaryon); //>[150609SH: add]
     
     gsl_interp_accel_free (acc_tot);
     gsl_interp_accel_free (acc_cdm);
     gsl_interp_accel_free (acc_baryon);
+    gsl_interp_accel_free (acc_vcdm); //>[150609SH: add]
+    gsl_interp_accel_free (acc_vbaryon); //>[150609SH: add]
   }
   
   // linear interpolation in log-log
   inline double extrap_right( double k, const tf_type& type )
   {
-      double v1(1.0), v2(1.0);
-      
     int n=m_tab_k.size()-1, n1=n-1;
-      switch( type )
-      {
-        case cdm:
-          v1 = m_tab_Tk_cdm[n1];
-          v2 = m_tab_Tk_cdm[n];
-          break;
-        case baryon:
-          v1 = m_tab_Tk_baryon[n1];
-          v2 = m_tab_Tk_baryon[n];
-          break;
-        case vcdm:
-        case vbaryon:
-        case total:
-          v1 = m_tab_Tk_tot[n1];
-          v2 = m_tab_Tk_tot[n];
-          break;
 
-        default:
+    double v1(1.0), v2(1.0);
-          throw std::runtime_error("Invalid type requested in transfer function evaluation");
-      }
 
     double lk   = log10(k);
     double dk   = m_tab_k[n]-m_tab_k[n1];
     double delk = lk-m_tab_k[n];
 
+    switch( type )
+      {
+      case cdm:
+	v1 = m_tab_Tk_cdm[n1];
+	v2 = m_tab_Tk_cdm[n];
 	return pow(10.0,(v2-v1)/dk*(delk)+v2);
+      case baryon:
+	v1 = m_tab_Tk_baryon[n1];
+	v2 = m_tab_Tk_baryon[n];
+	if( m_linbaryoninterp)
+	  return std::max((v2-v1)/dk*(delk)+v2,tiny);
+	return pow(10.0,(v2-v1)/dk*(delk)+v2);
+      case vcdm: //>[150609SH: add]
+	v1 = m_tab_Tvk_baryon[n1];
+	v2 = m_tab_Tvk_baryon[n];
+	return pow(10.0,(v2-v1)/dk*(delk)+v2);
+      case vbaryon: //>[150609SH: add]
+	v1 = m_tab_Tvk_baryon[n1];
+	v2 = m_tab_Tvk_baryon[n];
+	if( m_linbaryoninterp)
+	  return std::max((v2-v1)/dk*(delk)+v2,tiny);
+	return pow(10.0,(v2-v1)/dk*(delk)+v2);
+      case total:
+	v1 = m_tab_Tk_tot[n1];
+	v2 = m_tab_Tk_tot[n];
+	return pow(10.0,(v2-v1)/dk*(delk)+v2);
+      default:
+	throw std::runtime_error("Invalid type requested in transfer function evaluation");
+      }
+
+    return 0.0;
   }
   
   inline double compute( double k, tf_type type )
   {
     // use constant interpolation on the left side of the tabulated values
-        if( k < m_kmin )
+    if( k < m_kmin ){
-        {
+      switch( type ){
-          if( type == cdm )
+      case cdm:
 	return pow(10.0,m_tab_Tk_cdm[0]);
-          
+      case baryon:
-          else if( type == baryon )
+	if( m_linbaryoninterp )
+	  return m_tab_Tk_baryon[0];
 	return pow(10.0,m_tab_Tk_baryon[0]);
-          
+      case vcdm:
+	return pow(10.0,m_tab_Tvk_cdm[0]);
+      case vbaryon:
+	if( m_linbaryoninterp )
+	  return m_tab_Tvk_baryon[0];
+	return pow(10.0,m_tab_Tvk_baryon[0]);
+      case total:
 	return pow(10.0,m_tab_Tk_tot[0]);
-          
+      default:
+	throw std::runtime_error("Invalid type requested in transfer function evaluation");
+      }
     }
     // use linear interpolation on the right side of the tabulated values
     else if( k>m_kmax )
@@ -190,13 +256,24 @@ public:
     
     
     double lk = log10(k);
-		if( type == cdm )
+    switch( type ){
+    case cdm:
       return pow(10.0, gsl_spline_eval (spline_cdm,    lk, acc_cdm) );
-
+    case baryon:
-		if( type == baryon )
+      if( m_linbaryoninterp )
+	return gsl_spline_eval (spline_baryon, lk, acc_baryon);
       return pow(10.0, gsl_spline_eval (spline_baryon, lk, acc_baryon) );
-		
+    case vcdm:
+      return pow(10.0, gsl_spline_eval (spline_vcdm,    lk, acc_vcdm) );
+    case vbaryon:
+       if( m_linbaryoninterp )
+	 return gsl_spline_eval (spline_vbaryon, lk, acc_vbaryon);
+      return pow(10.0, gsl_spline_eval (spline_vbaryon, lk, acc_vbaryon) );
+    case total:
       return pow(10.0, gsl_spline_eval (spline_tot,    lk, acc_tot) );  
+    default:
+      throw std::runtime_error("Invalid type requested in transfer function evaluation");
+    }
   }
   
   inline double get_kmin( void ){