source: Sophya/trunk/Cosmo/SimLSS/cmvtstsch.cc@ 3315

#include "sopnamsp.h"
#include "machdefs.h"
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <unistd.h>
#include "timing.h"
#include "histos.h"
#include "ntuple.h"
#include "srandgen.h"
#include "perandom.h"

#include "geneutils.h"
#include "cosmocalc.h"
#include "schechter.h"

void usage(void);
void usage(void)
{
 cout<<"cmvtstsch -m xmin,xmax -a -i perc,dlxinc,dlxmax,glorder -n nbin -N nalea"<<endl
     <<"  xmin,xmax : limites en masse"<<endl
     <<"  perc : see IntegrateFunc"<<endl
     <<"  dlnx : see IntegrateFunc"<<endl
     <<"         si <0 nombre de point dlx = (log10(xmax)-log10(xmin))/(-dlx)"<<endl
     <<"  dlxmax : see IntegrateFunc, si <0 -dlxmax*dlx, si ==0 no check"<<endl
     <<"  glorder : odre gauss-legendre"<<endl
     <<"  nbin : nombre de bins"<<endl
     <<"  nalea : nombre de tirages aleatoires"<<endl
     <<"  -a : init auto de l\'aleatoire"<<endl;
}


int main(int narg,char *arg[])
{
 double h75 = 0.71 / 0.75;
 double nstar = 0.006*pow(h75,3.);  //  
 double mstar = pow(10.,9.8/(h75*h75));  // MSol
 double alpha = -1.37;
 cout<<"h75= "<<h75<<" nstar= "<<nstar<<"  mstar="<<mstar<<"  alpha="<<alpha<<endl;

 double xmin=1e6, xmax=1e12;
 double perc=0.01,dlxinc=-1000.,dlxmax=0.; unsigned short glorder=4;
 int npt = 1000;
 int nalea = 100000;

 char c;
 while((c = getopt(narg,arg,"ham:i:n:N:")) != -1) {
   switch (c) {
   case 'm' :
     sscanf(optarg,"%lf,%lf",&xmin,&xmax);
     break;
   case 'i' :
     int ig;
     sscanf(optarg,"%lf,%lf,%lf,%u",&perc,&dlxinc,&dlxmax,&ig);
     if(perc<=0.) perc = 0.01;
     glorder = (ig<=0) ? 4: ig;
     break;
   case 'n' :
     sscanf(optarg,"%d",&npt);
     if(npt<=0) npt = 1000;
     break;
   case 'N' :
     sscanf(optarg,"%d",&nalea);
     break;
   case 'a' :
     Auto_Ini_Ranf(5);
     break;
   case 'h' :
   default :
     usage(); return -1;
   }
 }

 double lnx1=log10(xmin), lnx2=log10(xmax), dlnx = (lnx2-lnx1)/npt;
 cout<<"xmin="<<xmin<<" ("<<lnx1<<") xmax="<<xmax<<" ("<<lnx2<<"), dlnx="<<dlnx<<endl;

 if(dlxinc<0.) dlxinc=(lnx2-lnx1)/(-dlxinc);
 if(dlxmax<0.) dlxmax *= -dlxinc; else if(dlxmax==0.) dlxmax=-1.;
 cout<<"perc="<<perc<<" dlxinc="<<dlxinc<<" dlxmax="<<dlxmax<<" glorder="<<glorder<<endl;

 cout<<"npt="<<npt<<",  nalea="<<nalea<<endl;

 Schechter sch(nstar,mstar,alpha);

 InitTim();

 //-------dn/dm
 cout<<endl;
 cout<<"sch(mstar) = "<<sch(mstar)<<" /Mpc^3/Msol"<<endl;
 sch.SetOutValue(0);
 cout<<"......"<<endl;
 for(double lnx=lnx1-3.;lnx<=lnx2-1.;lnx+=1.) {
   double num = IntegrateFuncLog(sch,lnx,lnx2,perc,dlxinc,dlxmax,glorder);
   cout<<"["<<lnx<<","<<lnx2<<"] integrated number = "<<num<<" Mpc^-3"<<endl;
 }
 cout<<"......"<<endl;
 for(double lnx=lnx1+1.;lnx<=lnx2+3.;lnx+=1.) {
   double num = IntegrateFuncLog(sch,lnx1,lnx,perc,dlxinc,dlxmax,glorder);
   cout<<"["<<lnx1<<","<<lnx<<"] integrated number = "<<num<<" Mpc^-3"<<endl;
 }
 cout<<"......"<<endl;
 for(double lnx=lnx1-3.;lnx<=lnx2+3.;lnx+=1.) {
   double num = IntegrateFuncLog(sch,lnx,lnx+1.,perc,dlxinc,dlxmax,glorder);
   cout<<"["<<lnx<<","<<lnx+1.<<"] integrated number = "<<num<<" Mpc^-3"<<endl;
 }
 Histo hdndm(lnx1,lnx2,npt); hdndm.ReCenterBin();
 FuncToHisto(sch,hdndm,true);

 //-------m*dn/dm
 cout<<endl;
 sch.SetOutValue(1);
 cout<<"mstar*sch(mstar) = "<<sch(mstar)<<" Msol/Mpc^3/Msol"<<endl;
 cout<<"......"<<endl;
 for(double lnx=lnx1-3.;lnx<=lnx2-1.;lnx+=1.) {
   double sum = IntegrateFuncLog(sch,lnx,lnx2,perc,dlxinc,dlxmax,glorder);
   cout<<"["<<lnx<<","<<lnx2<<"] integrated mass = "<<sum<<" Msol.Mpc^-3"<<endl;
 }
 cout<<"......"<<endl;
 for(double lnx=lnx1+1.;lnx<=lnx2+3.;lnx+=1.) {
   double sum = IntegrateFuncLog(sch,lnx1,lnx,perc,dlxinc,dlxmax,glorder);
   cout<<"["<<lnx1<<","<<lnx<<"] integrated mass = "<<sum<<" Msol.Mpc^-3"<<endl;
 }
 cout<<"......"<<endl;
 for(double lnx=lnx1-3.;lnx<=lnx2+3.;lnx+=1.) {
   double sum = IntegrateFuncLog(sch,lnx,lnx+1.,perc,dlxinc,dlxmax,glorder);
   cout<<"["<<lnx<<","<<lnx+1.<<"] integrated mass = "<<sum<<" Msol.Mpc^-3"<<endl;
 }
 Histo hmdndm(lnx1,lnx2,npt); hmdndm.ReCenterBin();
 FuncToHisto(sch,hmdndm,true);

 //------- Le flux HI pour M=mstar a 1 Gpc
 cout<<endl;
 double d = 1000.;
 double f = Msol2FluxHI(mstar,d);
 cout<<"FluxHI: d="<<d<<" Mpc:"<<"  M= "<<mstar<<" Msol  "
     <<"-> Flux= "<<f<<" W/m^2"<<endl;
 cout<<"Check... MassHI: d="<<d<<" Mpc:"<<"  f= "<<f<<" W/m^2  "
     <<"-> Mass= "<<FluxHI2Msol(f,d)<<" Msol"<<endl;

 //------- Random
 Histo hdna(hdndm); hdna.Zero();
 Histo hmdna(hmdndm); hmdna.Zero();
 FunRan tirhdndm = FunRan(hdndm,true);
 FunRan tirhmdndm = FunRan(hmdndm,true);
 if(nalea>0) {
   double summ=0.;
   PrtTim("--- avant tirage aleatoire");
   for(int i=0;i<nalea;i++) {
     double a = tirhdndm.RandomInterp();
     hdna.Add(a);
     a = tirhmdndm.RandomInterp();
     summ += pow(10.,a);
     hmdna.Add(a);
   }
   cout<<"Tirgae aleatoire: <m>="<<summ/(double)nalea<<endl;
   hdna *= hdndm.Sum()/hdna.Sum();
   hmdna *= hmdndm.Sum()/hmdna.Sum();
   PrtTim("--- apres tirage aleatoire");
 }

 //-------
 const int n = 3;
 char *vname[n] = {"m","f","mf"};
 NTuple nt(n,vname);
 double xnt[n];
 for(double lx=lnx1;lx<lnx2+dlnx/2.;lx+=dlnx) {
   double x = pow(10.,lx);
   xnt[0]  = x;
   sch.SetOutValue(0);
   xnt[1]  = sch(x);
   sch.SetOutValue(1);
   xnt[2]  = sch(x);
   nt.Fill(xnt);
 }

 cout<<"Ecriture"<<endl;
 string tag = "cmvtstsch.ppf";
 POutPersist pos(tag);
 tag = "nt"; pos.PutObject(nt,tag);
 tag = "hdndm"; pos.PutObject(hdndm,tag);
 tag = "hmdndm"; pos.PutObject(hmdndm,tag);
 tag = "hdna"; pos.PutObject(hdna,tag);
 tag = "hmdna"; pos.PutObject(hmdna,tag);
 Histo hdum1(tirhdndm);
 tag = "tirhdndm"; pos.PutObject(hdum1,tag);
 Histo hdum2(tirhmdndm);
 tag = "tirhmdndm"; pos.PutObject(hdum2,tag);

 return 0;
}

/*
openppf cmvtstsch.ppf

#--- la fonction de masse numerique
zone 2 2
n/plot nt.f%m ! ! "nsta crossmarker3"
n/plot nt.f%log10(m) ! ! "nsta crossmarker3"
n/plot nt.log10(f)%log10(m) f>0. ! "nsta crossmarker3"

# ce qu'on integre (en log10)
n/plot nt.m*f%log10(m) ! ! "nsta crossmarker3"

#--- la fonction de masse
zone 2 2
n/plot nt.mf%m ! ! "nsta crossmarker3"
n/plot nt.mf%log10(m) f>0. ! "nsta crossmarker3"
n/plot nt.log10(mf)%log10(m) mf>0. ! "nsta crossmarker3"

# ce qu'on integre (en log10)
n/plot nt.m*mf%log10(m) ! ! "nsta crossmarker3"

#--- Check histos
zone 2 1
disp hdndm "red"
n/plot nt.f%log10(m) ! ! "nsta crossmarker3 same"

disp hmdndm "red"
n/plot nt.m*f%log10(m) ! ! "nsta crossmarker3 same"

#--- Check tirage aleatoire
zone 2 1
disp hdndm
disp hdna "same red"

disp hmdndm
disp hmdna "same red"

zone
disp tirhdndm
disp tirhmdndm "same red"

 */