source: Sophya/trunk/Cosmo/SimLSS/schechter.cc@ 3892

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

#include "pexceptions.h"
#include "histos.h"
#include "perandom.h"
#include "tvector.h"
#include "fioarr.h"

#include "constcosmo.h"
#include "geneutils.h"
#include "schechter.h"

namespace SOPHYA {

///////////////////////////////////////////////////////////
//***************** Schechter Functions *****************//
///////////////////////////////////////////////////////////

// HI mass function:
//    see M.Zwaan astroph-0502257 (MNRAS Letters, Volume 359, Issue 1, pp. L30-L34.)
//    see M.Zwaan astroph-9707109 (ApJ, Volume 509, Issue 2, pp. 687-702.)

Schechter::Schechter(double nstar,double mstar,double alpha)
  : nstar_(nstar) , mstar_(mstar) , alpha_(alpha) , outvalue_(0)
{
}

Schechter::Schechter(Schechter& f)
  : nstar_(f.nstar_) , mstar_(f.mstar_) , alpha_(f.alpha_) , outvalue_(f.outvalue_)
{
}

Schechter::Schechter(void)
  : nstar_(0.) , mstar_(0.) , alpha_(0.) , outvalue_(0)
{
}

Schechter::~Schechter(void)
{
}

void Schechter::SetOutValue(unsigned short outvalue)
// outvalue = 0 : give dn/dm
//          = 1 : give m*dn/dm
{
  if(outvalue>1) {
    cout<<"Schechter::SetOutValue: Error bad outvalue: "<<outvalue<<endl;
    throw ParmError("Schechter::SetOutValue: Error bad outvalue");
  }
  outvalue_ = outvalue;
}

unsigned short Schechter::GetOutValue(void)
{
  return outvalue_;
}

void Schechter::SetParam(double nstar,double mstar,double alpha)
{
  nstar_ = nstar; mstar_ = mstar; alpha_ = alpha;
}

void Schechter::GetParam(double& nstar,double& mstar,double& alpha)
{
  nstar = nstar_; mstar = mstar_; alpha = alpha_;
}

double Schechter::operator() (double m)
// Return : "dn/dm = f(m)" or "m*dn/dm = f(m)"
{
  double x = m/mstar_;
  double dndm = nstar_/mstar_ * pow(x,alpha_) * exp(-x);
  if(outvalue_) dndm *= m;  // on veut m*dn/dm
  return dndm;
}


double Schechter::Integrate(double massmin,double massmax,int npt)
// Integrate from massmin to massmax with at least npt points log10 spaced
{
 if(massmin<=0. || massmax<=0. || massmax<=massmin) return 0.;
 if(npt<1) npt = 100;
 double lmassmin = log10(massmin), lmassmax = log10(massmax);
 double perc=0.01, dlxinc=(lmassmax-lmassmin)/npt, dlxmax=10.*dlxinc; unsigned short glorder=4;
 double sum = IntegrateFuncLog(*this,lmassmin,lmassmax,perc,dlxinc,dlxmax,glorder);
 return sum;
}

void Schechter::Print(void)
{
  cout<<"Schechter::Print: nstar="<<nstar_<<" Mpc^-3"
      <<"  mstar="<<mstar_<<" MSol"
      <<"  alpha="<<alpha_
      <<"  (outvalue="<<outvalue_<<" -> return ";
  if(outvalue_) cout<<"m*dn/dm)"; else cout<<"dn/dm)";
  cout<<endl;
}

///////////////////////////////////////////////////////////
//******* Les facilites pour tirer sur Schechter ********//
///////////////////////////////////////////////////////////

SchechterMassDist::SchechterMassDist(Schechter sch,double massmin,double massmax,int nbinmass)
// On veut une fonction de tirage aleatoire de la masse de N galaxies sur une Schechter
// Pour optimiser on cree (eventuellement) des histos de tirage
// ATTENTION: on donne les limites en masses pour les histos mais leurs abscisses
//   sont en log10(masse): histo [log10(massmin),log10(massmax)] avec nbinmass bins
// Si nbinmass<0 alors c'est le nombre de points par decade
  : sch_(sch) , sch_outvalue_(sch.GetOutValue())
  , massmin_(massmin) , massmax_(massmax) , nbinmass_(nbinmass)
  , ngalmin_(0) , ngalmax_(0) , nvalngal_(0)
  , ntrial_dir(0) , ntrial_tab(0)
  , hmdndm_(NULL) , tirhmdndm_(NULL)
{
  if(massmin_>massmax_  || massmin_<=0. || massmax_<=0.|| nbinmass_==0) {
    cout<<"SchechterMassDist::SchechterMassDist: error in input values"<<endl;
    throw ParmError("SchechterMassDist::SchechterMassDist: error in input values");
  }

  // Creation du tirage aleatoire sur la Schechter
  sch_outvalue_ = sch.GetOutValue();  // on sauve la configuration de la schechter
  sch.SetOutValue(1);  // on veut m*dN/dm
  double lnx1 = log10(massmin_), lnx2 = log10(massmax_); // Binning en log10 de la masse
  if(nbinmass_<0) nbinmass_ = int((-nbinmass_)*(lnx2-lnx1+1.));
  hmdndm_ = new Histo(lnx1,lnx2,nbinmass_); hmdndm_->ReCenterBin();
  FuncToHisto(sch,*hmdndm_,true);  // true -> bin en log10(x)
  tirhmdndm_ = new FunRan(*hmdndm_,true);  // true -> histo est pdf
  sch.SetOutValue(sch_outvalue_);  // on remet a la valeur initiale
}

SchechterMassDist::SchechterMassDist(void)
  : sch_outvalue_(0)
  , massmin_(0.) , massmax_(0.) , nbinmass_(0)
  , ngalmin_(0) , ngalmax_(0) , nvalngal_(0)
  , ntrial_dir(0) , ntrial_tab(0)
  , hmdndm_(NULL) , tirhmdndm_(NULL)
{
}

SchechterMassDist::~SchechterMassDist(void)
{
 Delete();
}

void SchechterMassDist::Delete(void)
{
 if(hmdndm_) delete hmdndm_;
 if(tirhmdndm_) delete tirhmdndm_;
 hmass_.resize(0);
 tmass_.resize(0);
 ntrial_dir = ntrial_tab = 0;
}

int SchechterMassDist::GetMassLim(double& massmin,double& massmax)
{
 massmin = massmin_;
 massmax = massmax_;
 return nbinmass_;
}

int SchechterMassDist::SetNgalLim(int ngalmax,int ngalmin,unsigned long nalea)
// Creation des Histos de tirage pour ngalmin a ngalmax galaxies
{
 int lp=2;
 ngalmin_=ngalmax_=nvalngal_=0;
 if(ngalmin<=0) ngalmin=1;
 if(ngalmax<ngalmin || ngalmax==1) return 0;
 ngalmin_ = ngalmin;
 ngalmax_ = ngalmax;
 nvalngal_ = ngalmax-ngalmin+1;
  
 if(nalea<1) nalea = 100000;

 if(lp>0) cout<<"SchechterMassDist::SetNgalLim: ngal=["
              <<ngalmin_<<","<<ngalmax_<<"] n="<<nvalngal_
              <<" filling with "<<nalea<<" trials"<<endl;

 //------- Construct histo
 double lnx1 = log10(massmin_), lnx2 = log10(massmax_);
 if(lp>0) cout<<"> Creating "<<nvalngal_<<" histos ["<<lnx1<<","<<lnx2<<"] n="<<nbinmass_<<endl;
 for(int i=ngalmin_;i<=ngalmax_;i++) {
   Histo h(*hmdndm_); h.Zero();
   hmass_.push_back(h);
 }
 if(lp>1) cout<<"...number of histos is "<<hmass_.size()<<endl;

 //------- Remplissage random
 sch_.SetOutValue(1);  // on veut m*dN/dm
 int lpmod = nalea/20; if(lpmod<=0) lpmod=1;
 double s1=0., sc1=0.; unsigned long ns1=0;
 double sax=0., scax=0.; unsigned long nsax=0;
 for(unsigned long ia=0;ia<nalea;ia++) {
   if(lp>1 && ia%lpmod==0) cout<<"...tirage "<<ia<<endl;
   double sum = 0.;
   for(int i=1;i<=ngalmax_;i++) {
     //double l10m = tirhmdndm_->Random();
     double l10m = tirhmdndm_->RandomInterp();
     double m = pow(10.,l10m);
     sum += m;
     s1 += m; sc1 += m*m; ns1++;
     int ipo = i-ngalmin_;
     if(ipo<0) continue;
     // ATTENTION: l'histo de tirage stoque le log10(moyenne=somme_masses/ngal).
     //            Ca permet d'avoir le meme binning quelque soit ngal
     double v = log10(sum/(double)i);
     hmass_[ipo].Add(v);
     if(i==ngalmax) {sax += sum/(double)i; scax += sum/(double)i*sum/(double)i; nsax++;}
   }
 }
 sch_.SetOutValue(sch_outvalue_);  // on remet a la valeur initiale

 if(ns1>1) {
   s1 /= ns1; sc1 = sc1/ns1 - s1*s1;
   cout<<"...Mean mass for ngal=1: "<<s1<<" ("<<log10(fabs(s1))<<")"
       <<" s="<<sqrt(fabs(sc1))<<" (ntrials "<<ns1<<")"<<endl;
 }
 if(nsax>1) {
   sax /= nsax; scax = scax/nsax - sax*sax;
   cout<<"...Mean mass for ngal="<<ngalmax_<<": "<<sax<<" ("<<log10(fabs(sax))<<")"
       <<" s="<<sqrt(fabs(scax))<<" (ntrials "<<nsax<<")"<<endl;
 }

 //------- Generation des classes de tirage aleatoire et des histos de verif
 if(lp>0) cout<<"> Creating "<<nvalngal_<<" FunRan"<<endl;
 for(unsigned int i=0;i<hmass_.size();i++) {
   FunRan t(hmass_[i],true);
   tmass_.push_back(t);
 }
 if(lp>1) cout<<"...number of funran is "<<tmass_.size()<<endl;

 return nvalngal_;
}

int SchechterMassDist::GetNgalLim(int& ngalmax,int& ngalmin)
{
 ngalmax = ngalmax_;
 ngalmin = ngalmin_;
 return nvalngal_;
}

Histo SchechterMassDist::GetHmDnDm(void) const
{
  return *hmdndm_;
}

FunRan SchechterMassDist::GetTmDnDm(void) const
{
  return *tirhmdndm_;
}

Histo SchechterMassDist::GetHisto(int i) const
{
  if(i<0 || i>=nvalngal_) {
    cout<<"SchechterMassDist::GetHisto: error in input values"<<endl;
    throw ParmError("SchechterMassDist::GetHisto: error in input values");
  }
  return hmass_[i];
}

FunRan SchechterMassDist::GetFunRan(int i) const
{
  if(i<0 || i>=nvalngal_) {
    cout<<"SchechterMassDist::GetFunRan: error in input values"<<endl;
    throw ParmError("SchechterMassDist::GetFunRan: error in input values");
  }
  return tmass_[i];
}

double SchechterMassDist::TirMass(int ngal)
{
 if(ngal<1) return 0.;

 int ipo = IndexFrNGal(ngal);
 double masse_des_ngal = 0.;
 if(ipo<0) {  // Pas d'histos de tirage pour ce nombre de galaxies
   for(long i=0;i<ngal;i++) {  // On tire donc ngal fois sur la Schechter
     // double lm = tirhmdndm_->Random();
     double lm = tirhmdndm_->RandomInterp();
     masse_des_ngal += pow(10.,lm);  // ATTENTION abscisse en log10(masse)
   }
   ntrial_dir++;
 } else {
   // ATTENTION l'histo de tirage stoque le log10(moyenne=somme_masses/ngal)
   //double lmngal = tmass_[ipo].Random();
   double lmngal = tmass_[ipo].RandomInterp();
   masse_des_ngal = pow(10.,lmngal) * ngal;
   ntrial_tab++;
 }

 return masse_des_ngal;
}

void SchechterMassDist::Print(void)
{
 cout<<"SchechterMassDist::Print: mass=["<<massmin_<<","<<massmax_<<"] n="<<nbinmass_<<endl
     <<"                 ngal=["<<ngalmin_<<","<<ngalmax_<<"] n="<<nvalngal_<<endl;
 sch_.Print();
}

void SchechterMassDist::PrintStatus(void)
{
 cout<<"SchechterMassDist::PrintStatus: number of trials: direct="<<ntrial_dir
     <<" tabulated="<<ntrial_tab<<endl;
}

void SchechterMassDist::WritePPF(string ppfname)
{
 char str[64];
 cout<<"SchechterMassDist::WritePPF into "<<ppfname<<endl;
 POutPersist pos(ppfname.c_str());

 double nstar,mstar,alpha;
 sch_.GetParam(nstar,mstar,alpha);
 TVector<r_8> tdum(20); tdum = 0.;
 tdum(0) = nstar;
 tdum(1) = mstar;
 tdum(2) = alpha;
 tdum(3) = sch_outvalue_;
 tdum(4) = massmin_;
 tdum(5) = massmax_;
 tdum(6) = nbinmass_;
 tdum(7) = ngalmin_;
 tdum(8) = ngalmax_;
 tdum(9) = nvalngal_;
 tdum(10) = hmass_.size();
 tdum(11) = tmass_.size();
 pos << PPFNameTag("SMDparam") << tdum;

 pos << PPFNameTag("SMDhmdndm") << *hmdndm_;
 pos << PPFNameTag("SMDtirhmdndm") << *tirhmdndm_;

 if(hmass_.size()>0) {
   for(unsigned int i=0;i<hmass_.size();i++) {
     sprintf(str,"SMDh%d",NGalFrIndex(i));
     pos << PPFNameTag(str) << hmass_[i];
   }
 }

 if(tmass_.size()>0) {
   for(unsigned int i=0;i<tmass_.size();i++) {
     sprintf(str,"SMDt%d",NGalFrIndex(i));
     Histo hdum(tmass_[i]);
     pos << PPFNameTag(str) << hdum;
   }
 }

}

void SchechterMassDist::ReadPPF(string ppfname)
{
 Delete(); // On des-alloue si deja remplit!

 char str[64];
 cout<<"SchechterMassDist::ReadPPF from "<<ppfname<<endl;
 PInPersist pis(ppfname.c_str());

 TVector<r_8> tdum;
 pis >> PPFNameTag("SMDparam") >> tdum;
 sch_.SetParam(tdum(0),tdum(1),tdum(2));
 sch_.SetOutValue((unsigned short)(tdum(3)+0.1));
 massmin_ = tdum(4);
 massmax_ = tdum(5);
 nbinmass_ = int(tdum(6)+0.1);
 ngalmin_ = int(tdum(7)+0.1);
 ngalmax_ = int(tdum(8)+0.1);
 nvalngal_ = int(tdum(9)+0.1);
 unsigned int nhmass = (unsigned int)(tdum(10)+0.1);
 unsigned int ntmass = (unsigned int)(tdum(11)+0.1);

 {
 Histo hdum;
 pis >> PPFNameTag("SMDhmdndm") >> hdum;
 hmdndm_ = new Histo(hdum);
 pis >> PPFNameTag("SMDtirhmdndm") >> hdum;
 tirhmdndm_ = new FunRan(hdum,false);
 }

 if(nhmass>0) {
   for(unsigned int i=0;i<nhmass;i++) {
     sprintf(str,"SMDh%d",NGalFrIndex(i));
     Histo hdum;
     pis >> PPFNameTag(str) >> hdum;
     hmass_.push_back(hdum);
   }
 }

 if(ntmass>0) {
   for(unsigned int i=0;i<ntmass;i++) {
     sprintf(str,"SMDt%d",NGalFrIndex(i));
     Histo hdum;
     pis >> PPFNameTag(str) >> hdum;
     FunRan fdum(hdum,false);
     tmass_.push_back(fdum);
   }
 }

}

///////////////////////////////////////////////////////////
//***************** Fonctions de Check ******************//
///////////////////////////////////////////////////////////

bool IsCompatible(Schechter& sch1,Schechter& sch2,double eps)
// on compare les differences a eps pres
{
  if(eps<=0.) eps=1.e-4;
  double nstar1,mstar1,alpha1;
  sch1.GetParam(nstar1,mstar1,alpha1);
  double nstar2,mstar2,alpha2;
  sch2.GetParam(nstar2,mstar2,alpha2);

  // nstar et mstar ne sont jamais nuls
  if(fabs(nstar1-nstar2)>fabs(nstar1+nstar2)/2.*eps) return false;
  if(fabs(mstar1-mstar2)>fabs(mstar1+mstar2)/2.*eps) return false;

  // alpha peut etre eventuellement nul
  if(fabs(alpha1)<1.e-100 && fabs(alpha2)<1.e-100 && fabs(alpha1-alpha2)>eps) return false;
  if(fabs(alpha1-alpha2)>fabs(alpha1+alpha2)/2.*eps) return false;
  return true;
}


///////////////////////////////////////////////////////////
//******************* Le Flux a 21cm ********************//
///////////////////////////////////////////////////////////
double Msol2LumiHI(double m)
// Input:
//    m : masse de HI en "Msol"
// Return:
//    la luminosite en W
// L(m) =  2.393e+18* m(en masse solaire) en W
{
  return 2.393e+18 * m;
}

double Msol2FluxHI(double m,double d)
// Input:
//    m : masse de HI en "Msol"
//    d : distance en "Mpc"  (si cosmo d=DLum(z))
// Return:
//    le flux total emis a 21 cm en W/m^2
// Ref:
// -- Binney & Merrifield, Galactic Astronomy p474 (ed:1998)
// S(W/m^2) = 1e-26 * Nu_21cm(Hz) * m(en masse solaire) /(2.35e5 * C(km/s) * Dlum^2)
//          = 2.0e-28 * m / Dlum^2
// -- J.Peterson & K.Bandura, astroph-0606104  (eq 7)
//    F.Abdalla & Rawlings, astroph-0411342 (eq 7 mais pb de d'unites?)
//          (A_21cm = 2.86888e-15 s^-1)
// S(W/m^2) = 3/4 * h(J.s) * Nu_21cm(Hz) * A_21cm(s^-1) * Msol(kg)/m_proton(kg)
//                         / (4 Pi) / (Dlum(Mpc) * 3.0857e22(m/Mpc))^2
//          = 2.0e-28 * m / Dlum^2
//-------------
{
  return  2.0e-28 * m / (d*d);
}

double FluxHI2Msol(double f,double d)
// Input:
//    f : flux total emis a 21 cm en W/m^2
//    d : distance en "Mpc"  (si cosmo d=DLum(z))
// Return:
//    m : masse de HI en "Msol"
{
  return f *d*d / 2.0e-28;
}


}  // Fin namespace SOPHYA