source: Sophya/trunk/Cosmo/SimLSS/agnjackson.cc@ 3196

#include "sopnamsp.h"
#include "machdefs.h"
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <vector>

#include "pexceptions.h"

#include "histos.h"
#include "srandgen.h"

#include "geneutils.h"
#include "agnjackson.h"


AGNJackson::AGNJackson(void)
  : nobjang_(0.) , fluxang_(0.) , dndls_(NULL), tirls_(NULL)
{

  //--- Les valeurs lues sur la courbe de Jackson (from JIM)
  const int nval = 10;
  // log10(S en Jy)
  double x[nval] = {-7.,-6.,-5.,-4.,-3.,-2.,-1.,0.,1.,2.};
  // S^2.5*dN/dS (en Jy^1.5 /sr)
  double y[nval] = {0.0035,0.14,0.95,2.5,7.,50.,260.,300.,150.,100.};

  if(sizeof(x)!=sizeof(y) || sizeof(y)!=8*nval) {
    cout<<"AGNJackson::AGNJackson_Error: incompatible x,y,nval sizes"<<endl;
    throw SzMismatchError("AGNJackson::AGNJackson_Error: incompatible x,y,nval sizes");
  }

  xjack_.resize(0); yjack_.resize(0);
  for(int i=0;i<nval;i++)
    {xjack_.push_back(x[i]); yjack_.push_back(y[i]);}

  // On cree l'histo: dndls = dN/dlog10(S) en fct de  log10(S en Jy)
  int nbin = nval*20;
  double xmin = x[0], xmax=x[nval-1];
  dndls_ = new Histo(xmin,xmax,nbin);
  dndls_->Zero();

  // On remplit les histos
  // yjack_: S^2.5 * dN/dS (en Jy^1.5 /sr) en fct de log10(S en Jy)
  // dndls: dN/dlog10(S)  (en 1/log10(Jy)/sr) en fct de log10(Jy en Jy)
  //   avec dN/dlog10(S) = dN/dS * dS/dlog10(S)
  //                     = (S^2.5*dN/dS)*S^(-2.5) * (ln(10)*S)
  // On va interpoler en log(dN/dlog10(S)) versus log10(S)
  vector<double> Y;
  for(int i=0;i<nval;i++) {
    double v = yjack_[i]*pow(10.,-1.5*xjack_[i])*M_LN10;
    Y.push_back(log(v));
  }
  for(int i=0;i<nbin;i++) {
    double lf = dndls_->BinCenter(i);
    // Interpolation lineaire (parabolic not OK!)
    double v = InterpTab(lf,xjack_,Y,1);
    (*dndls_)(i) = exp(v);
  }

  // On remplit la fonction de tirage
  tirls_ = new FunRan(*dndls_,true);

  // On calcul le nombre moyen d'AGN et le flux moyen par sr
  // ATTENTION: dn/dlog10(S) en fct de log10(s)
  nobjang_ = fluxang_ = 0.;
  for(int i=0;i<nbin;i++) {
    double lf = dndls_->BinCenter(i);
    nobjang_ += (*dndls_)(i);
    fluxang_ += (*dndls_)(i)*pow(10.,lf);
  }
  nobjang_ *= dndls_->BinWidth();
  fluxang_ *= dndls_->BinWidth();

}

AGNJackson::~AGNJackson(void)
{
  if(dndls_!=NULL) delete dndls_;
  if(tirls_!=NULL) delete tirls_;
}

void AGNJackson::Print(void)
{
  cout<<"AGNJackson::Print: nobj="<<nobjang_
      <<" /sr, flux="<<fluxang_<<" Jy/sr"<<endl;
  dndls_->Show();
}

void AGNJackson::OrigJack(vector<double>& xjack,vector<double>& yjack)
{
  xjack = xjack_;
  yjack = yjack_;
}