// 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