cib_interface.cc 4.9 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
// Ourselves
#include "cib_interface.h"

// Standard libraries
#include <iostream>
#include <cmath>

namespace cib {

  cib_interface::cib_interface()
  {
    _initialized_ = false;
  }

  cib_interface::~cib_interface()
  {
  }

  void cib_interface::initialize(// File to be parsed)
)
  {
    // law == "lognormal_sigevol"
    _parameters_.Meffmax = 8753289339381.791;
    _parameters_.etamax = 0.4028353504978569;
    _parameters_.sigmaMh = 1.807080723258688;
    //_parameters_.alpha = _parameters_.delta = _parameters_.beta = 0.0;
    _parameters_.tau = 1.2040244128818796;

    _cc_ = {1.076, 1.017, 1.119, 1.097, 1.068, 0.995, 0.960};
    _fc_ = {1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000};
    _nu_ = {100e9, 143e9, 217e9, 353e9, 545e9, 857e9, 3000e9}; // Hz
    _nfreq_ = _nu_.size();
    _ell_ = {187.,  320.,  502.,  684.,   890., 1158., 1505., 1956};
    _nell_ = _ell_.size();

    for (double m = 6; m <= 15; m+=0.1) _mh_.push_back(std::pow(10, m));

    cib::read_from_file("data/z.txt", _z_);

    _initialized_ = true;
  }

  void cib_interface::update(const cib::cib_parameters & parameters_)
  {
    if (! _initialized_) {
      throw std::logic_error("cib_interface::update: CIB model not initialized !");
    }
    _parameters_ = parameters_;
  }

  void cib_interface::compute_cls(cib::array3d & cls_) const
  {
    if (! _initialized_) {
      throw std::logic_error("cib_interface::update: CIB model not initialized !");
    }
    std::cout << "cib_interface::compute_cls" << std::endl;
    cib::zeros(cls_, _nfreq_, _nfreq_, _nell_);

    cib::array3d dj_cen; // [nfreq][mhalo][z]
    _compute_djc_dlnMh_(dj_cen);


    for (size_t l = 0; l < _nell_; l++) {
      for (size_t f1 = 0; f1 < _nfreq_; f1++) {
        for (size_t f2 = 0; f2 < _nfreq_; f2++) {

          // CIB

          // tSZ

          // CIBxtSZ


        } // end of f2 loop
      } // end of f1 loop
    } // end of ell loop
  }

  void cib_interface::_compute_djc_dlnMh_(cib::array3d & dj_cen_) const
  {
    const double fsub = 0.134;
    cib::array1d mhalo;
    for (size_t i = 0; i < _mh_.size(); i++) mhalo.push_back(_mh_[i]*(1-fsub));
    cib::array2d sfr;
    _compute_sfr_(mhalo, sfr);

    // cib::zeros(dj_cen_, this->_snu_eff_.size(), this->_mh_.size(), this->_z_.size());
    // for (size_t i = 0; i < this->_snu_eff_.size(); i++) {
    //   for (size_t j = 0; j < this->_mh_.size(); j++) {
    //     for (size_t k = 0; k < this->_z_.size(); k++) {
    //       const double z = this->_z_[k];
    //       const double cd_Mpc = this->_cosmo_->comoving_distances[k];
    //       const double rest = this->_hmfmz_[j][k]*sfr[j][k]*(1+z)*std::pow(cd_Mpc, 2)/KC;
    //       dj_cen_[i][j][k] = rest*this->_snu_eff_[i][k];
    //     }
    //   }
    // }
  }


  void cib_interface::_compute_sfr_(const cib::array1d & mhalo_, cib::array2d & sfr_) const
  {
    cib::array2d sfr_mhdot;
    _compute_sfr_mhdot_(mhalo_, sfr_mhdot);
    cib::dump(sfr_mhdot);
    // cib::array2d mhdot;
    // _compute_Mdot_(mhalo_, mhdot);
    // cib::array1d ob = this->_cosmo_->Ob;
    // cib::array1d om = this->_cosmo_->Om;
    // cib::zeros(sfr_, mhalo_.size(), this->_z_.size());
    // for (size_t i = 0; i < mhalo_.size(); i++) {
    //   for (size_t j = 0; j < this->_z_.size(); j++) {
    //     sfr_[i][j] = mhdot[i][j] * ob[i]/om[i] * sfr_mhdot[i][j];
    //   }
    // }
  }

  void cib_interface::_compute_sfr_mhdot_(const cib::array1d & mhalo_, cib::array2d & sfr_mhdot_) const
  {
    const double Meffmax = _parameters_.Meffmax;
    const double etamax = _parameters_.etamax;
    const double sigmaMh = _parameters_.sigmaMh;
    const double tau = _parameters_.tau;

    sfr_mhdot_.resize(mhalo_.size());
    for (size_t i = 0; i < mhalo_.size(); i++) {
      if (mhalo_[i] < Meffmax) {
        const double val = etamax * std::exp(-std::pow(std::log(mhalo_[i]/Meffmax)/sigmaMh, 2)/2);
        sfr_mhdot_[i].resize(_z_.size(), val);
      } else {
        sfr_mhdot_[i].resize(_z_.size());
        for (size_t j = 0; j < _z_.size(); j++) {
          const double z_c = 1.5;
          const double sigpow = sigmaMh - std::max(0., z_c-_z_[j])*tau;
          const double val = etamax * std::exp(-std::pow(std::log(mhalo_[i]/Meffmax)/sigpow, 2)/2);
          sfr_mhdot_[i][j] = val;
        } // end of z loop
      }
    } // end of halo mass loop
  }

  // void cl_cib_tsz_cross::Mdot(const cib::array1d & mhalo_, cib::array2d & mdot_) const
  // {
  //   cib::zeros(mdot_, mhalo_.size(), this->_z_.size());
  //   const bool use_mean = true;
  //   if (use_mean) {
  //     for (size_t i = 0; i < mhalo_.size(); i++) {
  //       const double b = std::pow(mhalo_[i]/1e12, 1.1);
  //       for (size_t j = 0; j < this->_z_.size(); j++) {
  //         const double a = 46.1*(1 + 1.11*this->_z_[j]) *       \
  //           std::sqrt(this->_cosmo_->Om0 * std::pow(1 + this->_z_[j], 3) + this->_cosmo_->Ode0);
  //         mdot_[i][j] = b*a;
  //       }
  //     }
  //   }
  // }



} // end of namespace cib