copy code from Tristram

src/camel/Astro/1halo_2halo.cc

@@ -13,6 +13,9 @@
 
 namespace cib {
   
+  /********************************************/
+  /* COSMOLOGY                                */
+  /********************************************/
   double cosmology::dchi_dz(double z)
   //Units in Mpc
   {
@@ -29,6 +32,8 @@ namespace cib {
     Om0 = klass.Omega_m();
     Ode0 = klass.Omega_L();
     Tcmb = klass.Tcmb();
+    rho0 = klass.critical_density(0);
+    redshift = z_;
     
     for( size_t i=0; i<z_.size(); i++) {
       comoving_distances.push_back( klass.com_distance(z_[i]));
@@ -48,6 +53,9 @@ namespace cib {
   }
   
 
+  /********************************************/
+  /* CL_CIB_TSZ_CROSS                         */
+  /********************************************/
   cl_cib_tsz_cross::cl_cib_tsz_cross(const cib::array2d & power_,
                                      const cib::array1d & z_,
                                      const cib::array1d & sig_z_arr_,
@@ -72,7 +80,6 @@ namespace cib {
                                      const cib::array1d & x_
                                      )
   {
-//     this->_kk_ = kk_;
     this->_power_ = power_;
     this->_z_ = z_;
     this->_sig_z_arr_ = sig_z_arr_;
@@ -603,4 +610,289 @@ namespace cib {
     }
   }
 
-} // end of namespace cib
+  void cl_cib_tsz_cross::hmf(const cib::array1d mass, 
+			     const cib::array2d k, const cib::array2d Pk, 
+			     cib::cosmology *cosmo, double delta_halo,
+			     cib::array2d & hmf_) const
+  {
+    cib::array1d zvalues = cosmo->redshift;
+
+    cib::zeros(hmf_, mass.size(), zvalues.size());
+    for( size_t iz=0; iz<zvalues.size(); iz++) {
+      double rho_z = cosmo->critical_densities[iz];
+      double Omega_m_z = cosmo->Om[iz];
+      cib::Filter filt = Filter( cosmo, k[iz], Pk[iz], zvalues[iz], rho_z, Omega_m_z, delta_halo);
+      for( size_t im=0; im<mass.size(); im++) {
+	hmf_[im][iz] = filt.dn_dlogm(mass[im]);
+      }
+    }
+  }
+  
+
+
+  void cl_cib_tsz_cross::biasdm( const cib::array1d mass, 
+				 const cib::array2d k, const cib::array2d Pk, 
+				 cib::cosmology *cosmo, double delta_halo,
+				 cib::array2d & bmz_) const
+  {
+    double y, A, aa, C, nuu;
+    cib::array1d zvalues = cosmo->redshift;
+    double B = 0.183, b = 1.5, c = 2.4;
+    double dc = 1.686;  /* neglecting the redshift evolution */
+
+    cib::zeros(bmz_, mass.size(), zvalues.size());
+    for( size_t iz=0; iz<zvalues.size(); iz++) {
+      double rho_z = cosmo->critical_densities[iz];
+      double Omega_m_z = cosmo->Om[iz];
+      cib::Filter filt = Filter( cosmo, k[iz], Pk[iz], zvalues[iz], rho_z, Omega_m_z, delta_halo);
+      for( size_t im=0; im<mass.size(); im++) {
+	y = std::log10(filt.delta_halo());
+	A = 1.0 + 0.24*y*std::exp(-std::pow(4./y,4));
+	aa = 0.44*y - 0.88;
+	C = 0.019 + 0.107*y + 0.19*std::exp(-std::pow(4./y,4));
+	nuu = filt.nu(mass[im]);
+	bmz_[im][iz] = 1. - (A*std::pow(nuu,aa)/(std::pow(nuu,aa) + std::pow(dc,aa))) + B*std::pow(nuu,b) + C*std::pow(nuu,c);
+      }
+    }
+  }
+  
+
+
+
+
+
+
+  /********************************************/
+  /* ONE HALO FILTER                          */
+  /********************************************/
+  /* for a redshift */
+  Filter::Filter( cib::cosmology *cosmo, 
+		  cib::array1d k,
+		  cib::array1d Pk,
+		  double z,
+		  double rho_z,      /* critical density at z */
+		  double Omega_m_z,
+		  double delta_h)
+  {
+    _k_ = k;
+    _Pk_ = Pk;
+    _z_ = z;
+    _Omz_ = Omega_m_z;
+    _rhoz_ = rho_z;
+    _delta_h_ = delta_h;
+    _mean_density0_ = cosmo->Om0 * cosmo->rho0;
+  }  
+
+  double Filter::r_delta( double m) const
+  /*
+    radius of the halo containing amount of matter corresponding to delta times
+    the critical density of the universe inside that halo
+  */
+  {
+    double r3 = 3.*m/(4*CLHEP::pi*_delta_h_*_rhoz_);
+    return( std::pow( r3, 1./3.));
+  }
+  
+  double Filter::mass_to_radius( double m)
+  /* Lagrangian radius of a dark matter halo */
+  {
+    double r3 = 3.*m/(4*CLHEP::pi*_mean_density0_);
+    return( std::pow(r3,1./3.) );
+  }
+
+  void Filter::W( cib::array1d rk, cib::array1d & wrk)
+  /* Fourier transform of top hat window function
+     limit of 1.4e-6 put as it doesn't add much to the final answer and
+     helps for faster convergence
+  */
+  {
+    double temp;
+
+    for( int i=0; i<rk.size(); i++) {
+      if( rk[i] > 1.4e-6) temp = 3.*(std::sin(rk[i])-rk[i]*std::cos(rk[i])) / std::pow(rk[i],3);
+      else temp = 1.;
+      wrk.push_back( temp);      
+    }
+  }
+
+  double Filter::sigma( double rad)
+  /*
+    matter variance for given power spectrum, wavenumbers k and radius of
+    the halo. radius is calculated below from mass of the halo. It has to
+    be noted that there's no factor of \Delta = 200 to calculate the radius.
+  */
+  {
+    cib::array1d rk;
+    cib::array1d rest;
+    for( int i=0; i<_k_.size(); i++) {
+      rk.push_back( rad*_k_[i]);
+      rest.push_back( _Pk_[i] * std::pow(_k_[i],3));
+    }
+    
+    // dk = 1e-7
+    const double dlnk = std::log(_k_[1] / _k_[0]);
+    cib::array1d integ;
+    cib::array1d wrk;
+    W(rk,wrk);
+    for( int i=0; i<_k_.size(); i++)
+      integ.push_back( rest[i] * std::pow(wrk[i],2));
+    
+    double sigm = 0.5 / std::pow(CLHEP::pi,2) * integration_simps(integ, dlnk);
+    return( std::sqrt(sigm));
+  }
+
+  double Filter::nu_delta( double m)
+  /*
+    to calculate the peak heights: nu_delta,  we use r_delta rather than the
+    simple Lagrangian radius calculated using mass_to_radius function.
+    This will be used in c-M relation to calculate the NFW profile.
+  */
+  {
+    double rad = r_delta(m);
+    double delta_c = 1.686;  /* critical density of the universe. Redshift evolution is small and neglected */
+    double sig = sigma(rad);
+    return( delta_c / sig); /* length of mass array */
+  }
+  
+  double Filter::nu( double m)
+  /*
+    to calculate the peak heights: nu, we use the
+    simple Lagrangian radius calculated using mass_to_radius function.
+  */
+  {
+    double rad = mass_to_radius(m);
+    double delta_c = 1.686;  /* critical density of the universe. Redshift evolution is small and neglected */
+    double sig = sigma(rad);
+    return( delta_c / sig); /* length of mass array */
+  }
+  
+  void Filter::dw_dlnkr( cib::array1d rk, cib::array1d &dw)
+  {
+    double dwdlnkr;
+
+    for( int i=0; i<rk.size(); i++) {
+      if( rk[i] > 1e-3) dwdlnkr = ( 9.*rk[i]*std::cos(rk[i]) + 3*std::sin(rk[i])*(std::pow(rk[i],2)-3.) ) / std::pow(rk[i],3);
+      else dwdlnkr = 0.;
+      dw.push_back( dwdlnkr);      
+    }
+  } 
+
+
+  double Filter::dlns2_dlnr( double rad)
+  {
+    cib::array1d rk, w, dw;
+    cib::array1d rest;
+    for( int i=0; i<_k_.size(); i++) {
+      rk.push_back( rad*_k_[i]);
+      rest.push_back( _Pk_[i] * std::pow(_k_[i],3));
+    }
+
+    W(rk,w);
+    dw_dlnkr(rk,dw);
+
+    cib::array1d inte;
+    for( int i=0; i<_k_.size(); i++)
+      inte.push_back( w[i]*dw[i]*rest[i]);
+
+    // dk = 1e-7
+    const double dlnk = std::log(_k_[1] / _k_[0]);
+    double s = sigma(rad);
+    return( integration_simps( inte, dlnk) / std::pow( CLHEP::pi*s, 2));
+  }
+  
+
+  /* The derivative of log radius with log mass. */
+  double Filter::dlnr_dlnm() 
+  {
+    return(1./3.);
+  }
+  
+  double Filter::dlns2_dlnm( double rad)
+  {
+    return( dlns2_dlnr(rad) * dlnr_dlnm());
+  }
+  
+  double Filter::dlns_dlnm( double rad)
+  {
+    return(0.5 * dlns2_dlnm(rad));
+  }
+  
+  double Filter::delta_halo() 
+  {
+    /* Overdensity of a halo w.r.t mean density */
+    //     return( _delta_h_);
+    /* Overdensity of a halo w.r.t critical density */
+    return( _delta_h_ / _Omz_);
+  }
+  
+  double Filter::fsigma( double rad)
+  /* https://arxiv.org/pdf/0803.2706.pdf
+     this is the function giving f(sigma)
+     All the paprameters mentioned here are taken for \Delta = 200
+     z is the redshift and redshift evolution has to be considered here
+     It has to be noted that values provided in Tinker 2008 paper are
+     wrt to mean background density. If we want to calculate the best fit
+     parameter values wrt to critical background density, we have to change
+     the corresponding delta_halo value and interpolate the values for the
+     corresponding delta_halo values.
+  */
+  {
+    cib::array1d delta_virs = {200, 300, 400, 600, 800, 1200, 1600, 2400, 3200};
+    cib::array1d A_array = { 1.858659e-01, 1.995973e-01, 2.115659e-01, 2.184113e-01, 2.480968e-01, 2.546053e-01, 2.600000e-01, 2.600000e-01, 2.600000e-01};
+    cib::array1d a_array = { 1.466904, 1.521782, 1.559186, 1.614585, 1.869936, 2.128056, 2.301275, 2.529241, 2.661983};
+    cib::array1d b_array = { 2.571104, 2.254217, 2.048674, 1.869559, 1.588649, 1.507134, 1.464374, 1.436827, 1.405210};
+    cib::array1d c_array = { 1.193958, 1.270316, 1.335191, 1.446266, 1.581345, 1.795050, 1.965613, 2.237466, 2.439729};
+    
+    double A_exp = -0.14;
+    double a_exp = 0.06;
+    double dhalo = delta_halo();
+    
+    //interp linear
+    int idx = 0;
+    while( delta_virs[idx+1] < dhalo) idx++;
+    if( idx >= delta_virs.size()) idx = delta_virs.size()-1;
+    double dx = double(dhalo - delta_virs[idx])/double(delta_virs[idx+1]-delta_virs[idx]);
+    double A0 = A_array[idx] + (A_array[idx+1]-A_array[idx])*dx;
+    double a0 = a_array[idx] + (a_array[idx+1]-a_array[idx])*dx;
+    double b0 = b_array[idx] + (b_array[idx+1]-b_array[idx])*dx;
+    double c0 = c_array[idx] + (c_array[idx+1]-c_array[idx])*dx;
+
+    double s = sigma(rad);
+    double A = A0 * std::pow(1.+_z_,A_exp);
+    double a = a0 * std::pow(1.+_z_,a_exp);
+    double alpha = std::pow(10, -std::pow(0.75/std::log10(dhalo/75.),1.2));
+    double b = b0 * std::pow(1.+_z_,alpha);
+    return( A * (std::pow(s/b,-a) + 1)*std::exp(-c0/std::pow(s,2)));
+  }
+  
+
+  double Filter::dn_dm( double m) {
+    double rad = mass_to_radius(m);
+    return( fsigma(rad) * _mean_density0_ * std::fabs( dlns_dlnm(rad)) / std::pow(m,2));
+  }
+
+  double Filter::dn_dlnm( double m)
+  {
+    return( m * dn_dm(m));
+  }
+  
+  double Filter::dn_dlogm( double m)
+  {
+    return( m * dn_dm(m) * std::log(10));
+  }
+  
+  double Filter::n_eff( double rad)
+  /*  Effective spectral index at scale of halo radius, ``len=len(m)``
+      Notes
+      -----
+      This function, and any derived quantities, can show small non-physical
+      'wiggles' at the 0.1% level, if too coarse a grid in ln(k) is used. If
+      applications are sensitive at this level, please use a very fine
+      k-space grid. Uses eq. 42 in Lukic et. al 2007. */
+  {
+    return( -3.0 * (2.0 * dlns_dlnm(rad) + 1.0));
+  }
+  
+
+
+}

src/camel/Astro/1halo_2halo.h

@@ -32,14 +32,16 @@ namespace cib {
     double H0;
     double Om0;
     double Ode0;
+    double rho0;
 
-    array1d comoving_distances;
-    array1d critical_densities;
-    array1d Ob;
-    array1d Om;
-    array1d Ez;
-    array1d dchidz;
-    array1d dVcdz;
+    cib::array1d redshift;
+    cib::array1d comoving_distances;
+    cib::array1d critical_densities;
+    cib::array1d Ob;
+    cib::array1d Om;
+    cib::array1d Ez;
+    cib::array1d dchidz;
+    cib::array1d dVcdz;
 
     double Tcmb;
     
@@ -136,6 +138,16 @@ namespace cib {
 
     void subhmf(const double mhalo_, const cib::array1d & ms_, cib::array1d & subhmf_) const;
 
+    void hmf(const cib::array1d mass, 
+	     const cib::array2d k, const cib::array2d Pk, 
+	     cib::cosmology * cosmo, double delta_halo,
+	     cib::array2d & hmf_) const;
+
+    void biasdm( const cib::array1d mass, 
+		 const cib::array2d k, const cib::array2d Pk, 
+		 cib::cosmology *cosmo, double delta_halo,
+		 cib::array2d & bmz_) const;
+
   private:
     // Default constructor can not instantiated
     cl_cib_tsz_cross();
@@ -176,6 +188,52 @@ namespace cib {
     double Kcmb_MJy[7] = {244.1, 371.74, 483.69, 287.45, 58.04, 2.27, 1.};
 
   };
-}// end of namespace cib
+
+
+  double mean_density0( ClassEngine & klass);
+  double mass_to_radius( double m, ClassEngine & klass);
+  double r_delta( double m, double z, ClassEngine & klass, double delta_h);
+
+  class Filter {
+  public:
+
+    //constuctor
+    Filter( cib::cosmology * cosmo, 
+	    cib::array1d k,
+	    cib::array1d Pk,
+	    double z,
+	    double rho_z,      /* critical density at z */
+	    double Omega_m_z,
+	    double delta_h);
+
+    double _delta_h_;
+    double _z_, _Omz_, _rhoz_;
+    double _mean_density0_;
+    cib::array1d _k_;
+    cib::array1d _Pk_;
+ 
+    double r_delta( double m) const;
+    double mass_to_radius( double m);
+    void W( cib::array1d rk, cib::array1d & wrk);
+    double sigma( double rad);
+    double nu_delta( double m);
+    double nu( double m);
+    void dw_dlnkr( cib::array1d rk, cib::array1d &dw);
+    double dlns2_dlnr( double rad);
+    double dlnr_dlnm();
+    double dlns2_dlnm( double rad);
+    double dlns_dlnm( double rad);
+    double delta_halo();
+    double fsigma( double rad);
+    double dn_dm( double m);
+    double dn_dlnm( double m);
+    double dn_dlogm( double m);
+    double n_eff( double rad);
+
+
+  };
+
+}
+
 
 #endif