Changeset 3196 in Sophya for trunk/Cosmo/SimLSS


Ignore:
Timestamp:
Apr 3, 2007, 4:06:50 PM (18 years ago)
Author:
cmv
Message:

les AGN selon C.Jackson, une premiere approche simplifiee, recodage from Jim Rich. cmv 03/04/2007

Location:
trunk/Cosmo/SimLSS
Files:
3 added
14 edited

Legend:

Unmodified
Added
Removed

  • trunk/Cosmo/SimLSS/Makefile

    r3157 r3196  
    1919       $(EXE)cmvtstsch $(EXE)cmvtstblack $(EXE)cmvtvarspec \
    2020       $(EXE)cmvdefsurv $(EXE)cmvobserv3d $(EXE)cmvtintfun \
    21        $(EXE)cmvtpoisson $(EXE)cmvconcherr $(EXE)cmvtinterp
     21       $(EXE)cmvtpoisson $(EXE)cmvconcherr $(EXE)cmvtinterp \
     22       $(EXE)cmvtstagn
    2223#$(EXE)cmvtluc
    2324 
    … …  
    2627          $(OBJ)cmvtstsch.o $(OBJ)cmvtstblack.o $(OBJ)cmvtvarspec.o $(OBJ)cmvdefsurv.o \
    2728          $(OBJ)cmvobserv3d.o $(OBJ)cmvtintfun.o $(OBJ)cmvtinterp.o \
    28           $(OBJ)cmvtpoisson.o $(OBJ)cmvconcherr.o $(OBJ)cmvtluc.o
     29          $(OBJ)cmvtpoisson.o $(OBJ)cmvconcherr.o $(OBJ)cmvtluc.o \
     30          $(OBJ)cmvtstagn.o
    2931 
    3032LIBROBJ = \
    3133          $(OBJ)cosmocalc.o $(OBJ)pkspectrum.o $(OBJ)schechter.o \
    32           $(OBJ)planckspectra.o $(OBJ)integfunc.o $(OBJ)geneutils.o \
     34          $(OBJ)planckspectra.o $(OBJ)geneutils.o $(OBJ)agnjackson.o \
    3335          $(OBJ)genefluct3d.o
    3436
    … …  
    5658$(OBJ)schechter.o: schechter.cc schechter.h constcosmo.h
    5759        $(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ schechter.cc
    58 $(OBJ)integfunc.o: integfunc.cc integfunc.h
    59         $(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ integfunc.cc
    6060$(OBJ)planckspectra.o: planckspectra.cc planckspectra.h constcosmo.h
    6161        $(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ planckspectra.cc
    62 $(OBJ)pkspectrum.o: pkspectrum.cc pkspectrum.h constcosmo.h integfunc.h
     62$(OBJ)pkspectrum.o: pkspectrum.cc pkspectrum.h constcosmo.h
    6363        $(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ pkspectrum.cc
    6464$(OBJ)geneutils.o: geneutils.cc geneutils.h
    6565        $(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ geneutils.cc
     66$(OBJ)agnjackson.o: agnjackson.cc agnjackson.h
     67        $(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ agnjackson.cc
    6668$(OBJ)genefluct3d.o: genefluct3d.cc genefluct3d.h
    6769        $(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ genefluct3d.cc
    … …  
    173175
    174176##############################################################################
     177cmvtstagn: $(EXE)cmvtstagn
     178        echo $@ " done"
     179$(EXE)cmvtstagn: $(OBJ)cmvtstagn.o $(LIB)libcmvsimbao.a
     180        $(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvtstagn.o $(MYLIB)
     181$(OBJ)cmvtstagn.o: cmvtstagn.cc
     182        $(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvtstagn.cc
     183
     184##############################################################################
    175185cmvtluc: $(EXE)cmvtluc
    176186        echo $@ " done"

  • trunk/Cosmo/SimLSS/cmvdefsurv.cc

    r3193 r3196  
    1111#include "cosmocalc.h"
    1212#include "geneutils.h"
    13 #include "integfunc.h"
    1413#include "schechter.h"
    1514#include "planckspectra.h"
    1615
    1716/* --- Check Peterson at al. astro-ph/0606104 v1
    18 cmvdefsurv -z 0.0025 -x 1 -U 0.75,0.3,0.7,-1,1 -V 300 -O 400000,6000 -A 75 -M 6.156e9 -F 3 -2 1.5
     17cmvdefsurv -z 0.0025 -x 1 -U 0.75,0.3,0.7,-1,1 -V 300 -O 400000,6000 -N 75 -M 6.156e9 -F 3 -2 1.5
    1918 --- */
    2019
    … …  
    3534      <<" -L lobewidth : taille du lobe d\'observation (FWHM) en arcmin (def= 1\')"<<endl
    3635      <<" -O surf,tobs : surface effective (m^2) et temps d\'observation (s)"<<endl
    37       <<" -A Tsys : temperature du system (K)"<<endl
     36      <<" -N Tsys : temperature du system (K)"<<endl
    3837      <<" -S Tsynch,indnu : temperature (K) synch a 408 Mhz, index d\'evolution"<<endl
    3938      <<"                   (indnu==0 no evolution with freq.)"<<endl
    … …  
    4342      <<" -U h100,om0,ol0,w0,or0,flat : cosmology"<<endl
    4443      <<" -2 : two polarisations measured"<<endl
     44      <<" -A <log10(S_agn)> : moyenne du flux AGN en Jy dans le pixel"<<endl
    4545      <<" redshift : redshift moyen du survey"<<endl
    4646      <<endl;
    … …  
    5959 double alpha = -1.37;
    6060 cout<<"nstar= "<<nstar<<"  mstar="<<mstar<<"  alpha="<<alpha<<endl;
    61 
    6261
    6362 // --- Arguments
    … …  
    7675 double vrot = 300.; // largeur en vitesse (km/s) pour elargissement doppler
    7776 double facpolar = 0.5; // si on ne mesure les 2 polars -> 1.0
     77 double lflux_agn = -3.;
    7878
    7979 // --- Decodage arguments
    8080 char c;
    81   while((c = getopt(narg,arg,"hP2x:y:z:A:S:O:M:F:V:U:L:")) != -1) {
     81  while((c = getopt(narg,arg,"hP2x:y:z:N:S:O:M:F:V:U:L:A:")) != -1) {
    8282  switch (c) {
    8383  case 'P' :
    … …  
    9999    sscanf(optarg,"%lf",&lobewidth);
    100100    break;
    101   case 'A' :
     101  case 'N' :
    102102    sscanf(optarg,"%lf",&Tsys);
    103103    break;
    … …  
    119119  case '2' :
    120120    facpolar = 1.0;
     121    break;
     122  case 'A' :
     123    sscanf(optarg,"%lf",&lflux_agn);
    121124    break;
    122125  case 'h' :
    … …  
    281284 double lobecyl = sqrt(8.)*slobe; // diametre du lobe cylindrique equiv en arcmin
    282285 double lobearea = M_PI*lobecyl*lobecyl/4.;  // en arcmin^2 (hypothese lobe gaussien)
     286 double nlobes = rad2min(adtx)*rad2min(adty)/lobearea;
     287 if(lobewidth<=0.) nlobes = 1.;
    283288 cout<<"\nBeam FWHM = "<<lobewidth<<"\' -> sigma = "<<slobe<<"\' -> "
    284289     <<" Dcyl = "<<lobecyl<<"\' -> area = "<<lobearea<<" arcmin^2"<<endl;
    285  double nlobes = rad2min(adtx)*rad2min(adty)/lobearea;
    286290 cout<<"Number of beams in one transversal pixel = "<<nlobes<<endl;
    287291
    … …  
    290294 cout<<"flux factor = "<<hifactor<<" at redshift = "<<redshift<<endl;
    291295
    292  double fhi = hifactor*FluxHI(mhiref,dlum);
     296 double fhi = hifactor*Msol2FluxHI(mhiref,dlum);
    293297 cout<<"FluxHI("<<dlum<<" Mpc) all polar:"<<endl
    294298     <<"  Flux= "<<fhi<<" W/m^2 = "<<fhi/Jansky2Watt_cst<<" Jy.Hz"<<endl
    295      <<"      in ["<<hifactor*FluxHI(mhiref,dlumlim[0])
    296      <<","<<hifactor*FluxHI(mhiref,dlumlim[1])<<"] W/m^2"<<endl;
     299     <<"      in ["<<hifactor*Msol2FluxHI(mhiref,dlumlim[0])
     300     <<","<<hifactor*Msol2FluxHI(mhiref,dlumlim[1])<<"] W/m^2"<<endl;
    297301 double sfhi = fhi / (dnuhiz*1e9) / Jansky2Watt_cst;
    298302 cout<<"If spread over "<<dnuhiz<<" GHz, flux density = "<<sfhi<<" Jy"<<endl;
    … …  
    342346 cout<<"                       flux density = "<<scmb<<" Jy"<<endl;
    343347
     348 // AGN
     349 double flux_agn = pow(10.,lflux_agn);
     350 cout<<"AGN: log10(S_agn)="<<lflux_agn<<" -> S_agn="<<flux_agn<<" Jy"<<endl;
     351 double flux_agn_pix = flux_agn*(dnuhiz*1e9);
     352 double mass_agn_pix = FluxHI2Msol(flux_agn_pix*Jansky2Watt_cst,dlum);
     353 double lmass_agn_pix = log10(mass_agn_pix);
     354 cout<<"...pixel: f="<<flux_agn_pix<<" 10^-26 W/m^2"
     355     <<" -> "<<mass_agn_pix<<" Msol -> log10 = "<<lmass_agn_pix<<endl;
     356
    344357 // --- Noise analysis
    345358 cout<<"\n--- Noise analysis"<<endl;

  • trunk/Cosmo/SimLSS/cmvobserv3d.cc

    r3193 r3196  
    1515#include "schechter.h"
    1616#include "geneutils.h"
    17 #include "integfunc.h"
    1817#include "genefluct3d.h"
    1918
    … …  
    3231     <<" -2 : compute 2D spectrum"<<endl
    3332     <<" -M schmin,schmax,nsch : min,max mass and nb points for schechter HI"<<endl
     33     <<" -A <log10(S_agn en Msol)>,sigma : distribution des AGN par pixel"<<endl
    3434     <<" -W : write cube in FITS format"<<endl
    3535     <<" -P : write cube in PPF format"<<endl
    … …  
    7777 double snoise= 0.;   // en equivalent MSol
    7878
     79 // *** AGN
     80 bool do_agn = false;
     81 double lmsol_agn=-99., lsigma_agn=0.;   // en equivalent MSol
     82
    7983 // *** type de generation
    8084 bool computefourier0=false;
    … …  
    9195
    9296 char c;
    93  while((c = getopt(narg,arg,"ha0PWV2Gx:y:z:s:Z:M:")) != -1) {
     97 while((c = getopt(narg,arg,"ha0PWV2Gx:y:z:s:Z:M:A:")) != -1) {
    9498  switch (c) {
    9599  case 'a' :
    … …  
    122126  case 'M' :
    123127    sscanf(optarg,"%lf,%lf,%d",&schmin,&schmax,&schnpt);
     128    break;
     129  case 'A' :
     130    do_agn = true;
     131    sscanf(optarg,"%lf,%lf",&lmsol_agn,&lsigma_agn);
    124132    break;
    125133  case 'V' :
    … …  
    153161     <<", schnpt="<<schnpt<<endl;
    154162 cout<<"snoise="<<snoise<<" equivalent Msol"<<endl;
     163 if(do_agn) cout<<"AGN: <log10(Msol)>="<<lmsol_agn<<" , sigma="<<lsigma_agn<<endl;
    155164
    156165 //-----------------------------------------------------------------
    … …  
    476485 }
    477486
     487 if(do_agn) {
     488   cout<<"\n--- Add AGN: <mass>="<<lmsol_agn<<" , sigma="<<lsigma_agn<<endl;
     489   fluct3d.AddAGN(lmsol_agn,lsigma_agn);
     490   nm = fluct3d.MeanSigma2(rm,rs2);
     491   cout<<"HI mass with AGN: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
     492       <<rs2<<" -> "<<sqrt(rs2)<<endl;
     493   PrtTim(">>>> End Add AGN");
     494 }
     495
    478496 if(snoise>0.) {
    479497   cout<<"\n--- Add noise to HI Flux snoise="<<snoise<<endl;

  • trunk/Cosmo/SimLSS/cmvtintfun.cc

    r3115 r3196  
    1010#include "ntuple.h"
    1111
    12 #include "integfunc.h"
     12#include "geneutils.h"
    1313
    1414// cmvtintfun [n] [xmin,xmax,npt]

  • trunk/Cosmo/SimLSS/cmvtstblack.cc

    r3115 r3196  
    1010#include "ntuple.h"
    1111
    12 #include "integfunc.h"
    1312#include "constcosmo.h"
     13#include "geneutils.h"
    1414#include "planckspectra.h"
    1515

  • trunk/Cosmo/SimLSS/cmvtstsch.cc

    r3120 r3196  
    1414
    1515#include "geneutils.h"
    16 #include "integfunc.h"
    1716#include "cosmocalc.h"
    1817#include "schechter.h"
    … …  
    113112 cout<<endl;
    114113 double d = 1000.;
     114 double f = Msol2FluxHI(mstar,d);
    115115 cout<<"FluxHI: d="<<d<<" Mpc:"<<"  M= "<<mstar<<" Msol  "
    116      <<"-> Flux= "<<FluxHI(mstar,d)<<" W/m^2"<<endl;
     116     <<"-> Flux= "<<f<<" W/m^2"<<endl;
     117 cout<<"Check... MassHI: d="<<d<<" Mpc:"<<"  f= "<<f<<" W/m^2  "
     118     <<"-> Mass= "<<FluxHI2Msol(f,d)<<" Msol"<<endl;
    117119
    118120 //------- Random

  • trunk/Cosmo/SimLSS/cosmocalc.cc

    r3193 r3196  
    1313#include "constcosmo.h"
    1414#include "cosmocalc.h"
    15 #include "integfunc.h"
     15#include "geneutils.h"
    1616
    1717

  • trunk/Cosmo/SimLSS/genefluct3d.cc

    r3157 r3196  
    2020
    2121#include "constcosmo.h"
    22 #include "integfunc.h"
    2322#include "geneutils.h"
    2423
    … …  
    970969}
    971970
     971void GeneFluct3D::AddAGN(double lmsol,double lsigma)
     972// Add AGN flux into simulation:
     973// --- Procedure:
     974// 1. lancer "cmvdefsurv" avec les parametres du survey
     975//    et recuperer l'angle solide "angsol sr" du pixel elementaire
     976//    au centre du cube.
     977// 2. lancer "cmvtstagn" pour cet angle solide -> cmvtstagn.ppf
     978// 3. regarder l'histo "hlfang" et en deduire un equivalent gaussienne
     979//    cad une moyenne <log10(S)> et un sigma "sig"
     980//    Attention: la distribution n'etant pas gaussienne
     981// 4. re-lancer "cmvdefsurv" en ajoutant l'info sur les AGN
     982//    "cmvdefsurv ... -G <log10(S)> ..."
     983//    et recuperer le log10(masse solaire equivalente)
     984// --- Limitations actuelle du code:
     985// a. l'angle solide du pixel est pris au centre du cube
     986//    et ne varie pas avec la distance a l'interieur du cube
     987// b. la taille en dNu des pixels (selon z) est supposee constante
     988//    et egale a celle calculee pour le centre du cube
     989// c. les AGN sont supposes plats en flux
     990// d. le flux des AGN est mis dans une colonne Oz (indice k)
     991//    et pas sur la ligne de visee
     992// e. la distribution est approximee a une gaussienne
     993// .. C'est une approximation pour un observateur loin du centre du cube
     994//    et pour un cube peu epais / distance observateur
     995// --- Parametres de la routine:
     996// lmsol : c'est le <log10(Msol)> qui est la conversion en masse solaire
     997//         du flux depose dans un pixel par les AGN
     998// lsigma : c'est le sigma de la distribution
     999// - Comme on est en echelle log10():
     1000//   on tire log10(Msol) + X
     1001//   ou X est une realisation sur une gaussienne de variance "sig^2"
     1002//   La masse realisee est donc: Msol*10^X
     1003// - Pas de probleme de pixel negatif car on a une multiplication!
     1004{
     1005  if(lp_>0) cout<<"--- AddAGN: lmsol = "<<lmsol<<" , sigma = "<<lsigma<<endl;
     1006  check_array_alloc();
     1007
     1008  double m = pow(10.,lmsol);
     1009
     1010  double sum=0., sum2=0., nsum=0.;
     1011  for(long l=0;l<Nz_;l++) {
     1012    double a = lsigma*NorRand();
     1013    a = m*pow(10.,a);
     1014    // On met le meme tirage le long de Oz (indice k)
     1015    for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) {
     1016      int_8 ip = IndexR(i,j,l);
     1017      data_[ip] += a;
     1018    }
     1019    sum += a; sum2 += a*a; nsum += 1.;
     1020  }
     1021
     1022 if(nsum>1.) {
     1023   sum /= nsum;
     1024   sum2 = sum2/nsum - sum*sum;
     1025   cout<<"...Mean mass="<<sum<<" Msol , s^2="<<sum2<<" s="<<sqrt(fabs(sum2))<<endl;
     1026 }
     1027 
     1028}
     1029
    9721030void GeneFluct3D::AddNoise2Real(double snoise)
    9731031// add noise to every pixels (meme les <=0 !)

  • trunk/Cosmo/SimLSS/genefluct3d.h

    r3157 r3196  
    8383  double TurnNGal2Mass(FunRan& massdist,bool axeslog=false);
    8484  double TurnMass2Flux(void);
     85  void AddAGN(double lmsol,double lsigma);
    8586  void AddNoise2Real(double snoise);
    8687

  • trunk/Cosmo/SimLSS/geneutils.cc

    r3157 r3196  
    150150}
    151151
     152//----------------------------------------------------
     153double InterpTab(double x0,vector<double>& X,vector<double>& Y,unsigned short typint)
     154// Interpole in x0 the table Y = f(X)
     155//           X doit etre ordonne par ordre croissant (strictement)
     156// typint = 0 : nearest value
     157//          1 : linear interpolation
     158//          2 : parabolique interpolation
     159{
     160 long n = X.size();
     161 if(Y.size()<n || n<2)
     162   throw ParmError("InterpTab_Error :  incompatible size between X and Y tables!");
     163
     164 if(x0<X[0] || x0>X[n-1]) return 0.;
     165 if(typint>2) typint = 0;
     166
     167 // Recherche des indices encadrants par dichotomie
     168 long k, klo=0, khi=n-1;
     169 while (khi-klo > 1) {
     170   k = (khi+klo) >> 1;
     171   if (X[k] > x0) khi=k; else klo=k;
     172 }
     173
     174 // Quel est le plus proche?
     175 k = (x0-X[klo]<X[khi]-x0) ? klo: khi;
     176
     177 // On retourne le plus proche
     178 if(typint==0) return Y[k];
     179
     180 // On retourne l'extrapolation lineaire
     181 if(typint==1 || n<3)
     182   return Y[klo] + (Y[khi]-Y[klo])/(X[khi]-X[klo])*(x0-X[klo]);
     183
     184 // On retourne l'extrapolation parabolique
     185 if(k==0) k++; else if(k==n-1) k--;
     186 klo = k-1; khi = k+1;
     187 double x1 = X[klo]-X[k], x2 = X[khi]-X[k];
     188 double y1 = Y[klo]-Y[k], y2 = Y[khi]-Y[k];
     189 double den = x1*x2*(x1-x2);
     190 double a = (y1*x2-y2*x1)/den;
     191 double b = (y2*x1*x1-y1*x2*x2)/den;
     192 x0 -= X[k];
     193 return Y[k] + (a*x0+b)*x0;;
     194
     195}
     196
    152197//-------------------------------------------------------------------
    153198int FuncToHisto(GenericFunc& func,Histo& h,bool logaxex)
    … …  
    252297  return n;
    253298}
     299
     300//-------------------------------------------------------------------
     301
     302static unsigned short IntegrateFunc_GlOrder = 0;
     303static vector<double> IntegrateFunc_x;
     304static vector<double> IntegrateFunc_w;
     305
     306///////////////////////////////////////////////////////////
     307//************** Integration of Functions ***************//
     308///////////////////////////////////////////////////////////
     309
     310
     311////////////////////////////////////////////////////////////////////////////////////
     312double IntegrateFunc(GenericFunc& func,double xmin,double xmax
     313       ,double perc,double dxinc,double dxmax,unsigned short glorder)
     314// --- Integration adaptative ---
     315//     Integrate[func[x], {x,xmin,xmax}]
     316// ..xmin,xmax are the integration limits
     317// ..dxinc is the searching increment x = xmin+i*dxinc
     318//         if <0  npt = int(|dxinc|)   (si<2 alors npt=100)
     319//                et dxinc = (xmax-xmin)/npt
     320// ..dxmax is the maximum possible increment (if dxmax<=0 no test)
     321// ---
     322// Partition de [xmin,xmax] en intervalles [x(i),x(i+1)]:
     323// on parcourt [xmin,xmax] par pas de "dxinc" : x(i) = xmin + i*dxinc
     324// on cree un intervalle [x(i),x(i+1)]
     325//     - si |f[x(i+1)] - f[x(i)]| > perc*|f[[x(i)]|
     326//     - si |x(i+1)-x(i)| >= dxmax  (si dxmax>0.)
     327// Dans un intervalle [x(i),x(i+1)] la fonction est integree
     328//   avec une methode de gauss-legendre d'ordre "glorder"
     329{
     330  double signe = 1.;
     331  if(xmin>xmax) {double tmp=xmax; xmax=xmin; xmin=tmp; signe=-1.;}
     332
     333  if(glorder==0) glorder = 4;
     334  if(perc<=0.) perc=0.1;
     335  if(dxinc<=0.) {int n=int(-dxinc); if(n<2) n=100; dxinc=(xmax-xmin)/n;}
     336  if(glorder != IntegrateFunc_GlOrder) {
     337    IntegrateFunc_GlOrder = glorder;
     338    Compute_GaussLeg(glorder,IntegrateFunc_x,IntegrateFunc_w,0.,1.);
     339  }
     340
     341  // Recherche des intervalles [x(i),x(i+1)]
     342  int_4 ninter = 0;
     343  double sum = 0., xbas=xmin, fbas=func(xbas), fabsfbas=fabs(fbas);
     344  for(double x=xmin+dxinc; x<xmax+dxinc/2.; x += dxinc) {
     345    double f = func(x);
     346    double dx = x-xbas;
     347    bool doit = false;
     348    if( x>xmax ) {doit = true; x=xmax;}
     349    else if( dxmax>0. && dx>dxmax ) doit = true;
     350    else if( fabs(f-fbas)>perc*fabsfbas ) doit = true;
     351    if( ! doit ) continue;
     352    double s = 0.;
     353    for(unsigned short i=0;i<IntegrateFunc_GlOrder;i++)
     354      s += IntegrateFunc_w[i]*func(xbas+IntegrateFunc_x[i]*dx);
     355    sum += s*dx;
     356    xbas = x; fbas = f; fabsfbas=fabs(fbas); ninter++;
     357  }
     358  //cout<<"Ninter="<<ninter<<endl;
     359
     360  return sum*signe;
     361}
     362
     363////////////////////////////////////////////////////////////////////////////////////
     364double IntegrateFuncLog(GenericFunc& func,double lxmin,double lxmax
     365         ,double perc,double dlxinc,double dlxmax,unsigned short  glorder)
     366// --- Integration adaptative ---
     367// Idem IntegrateFunc but integrate on logarithmic (base 10) intervals:
     368//    Int[ f(x) dx ] = Int[ x*f(x) dlog10(x) ] * log(10)
     369// ..lxmin,lxmax are the log10() of the integration limits
     370// ..dlxinc is the searching logarithmic (base 10) increment lx = lxmin+i*dlxinc
     371// ..dlxmax is the maximum possible logarithmic (base 10) increment (if dlxmax<=0 no test)
     372// Remarque: to be use if "x*f(x) versus log10(x)" looks like a polynomial
     373//            better than "f(x) versus x"
     374// ATTENTION: la fonction func que l'on passe en argument
     375//            est "func(x)" et non pas "func(log10(x))"
     376{
     377  double signe = 1.;
     378  if(lxmin>lxmax) {double tmp=lxmax; lxmax=lxmin; lxmin=tmp; signe=-1.;}
     379
     380  if(glorder==0) glorder = 4;
     381  if(perc<=0.) perc=0.1;
     382  if(dlxinc<=0.) {int n=int(-dlxinc); if(n<2) n=100; dlxinc=(lxmax-lxmin)/n;}
     383  if(glorder != IntegrateFunc_GlOrder) {
     384    IntegrateFunc_GlOrder = glorder;
     385    Compute_GaussLeg(glorder,IntegrateFunc_x,IntegrateFunc_w,0.,1.);
     386  }
     387
     388  // Recherche des intervalles [lx(i),lx(i+1)]
     389  int_4 ninter = 0;
     390  double sum = 0., lxbas=lxmin, fbas=func(pow(10.,lxbas)), fabsfbas=fabs(fbas);
     391  for(double lx=lxmin+dlxinc; lx<lxmax+dlxinc/2.; lx += dlxinc) {
     392    double f = func(pow(10.,lx));
     393    double dlx = lx-lxbas;
     394    bool doit = false;
     395    if( lx>lxmax ) {doit = true; lx=lxmax;}
     396    else if( dlxmax>0. && dlx>dlxmax ) doit = true;
     397    else if( fabs(f-fbas)>perc*fabsfbas ) doit = true;
     398    if( ! doit ) continue;
     399    double s = 0.;
     400    for(unsigned short i=0;i<IntegrateFunc_GlOrder;i++) {
     401      double y = pow(10.,lxbas+IntegrateFunc_x[i]*dlx);
     402      s += IntegrateFunc_w[i]*y*func(y);
     403    }
     404    sum += s*dlx;
     405    lxbas = lx; fbas = f; fabsfbas=fabs(fbas); ninter++;
     406  }
     407  //cout<<"Ninter="<<ninter<<endl;
     408
     409  return M_LN10*sum*signe;
     410}
     411
     412////////////////////////////////////////////////////////////////////////////////////
     413/*
     414Integration des fonctions a une dimension y=f(x) par la Methode de Gauss-Legendre.
     415  --> Calcul des coefficients du developpement pour [x1,x2]
     416| INPUT:
     417|  x1,x2 : bornes de l'intervalle (dans nbinteg.h -> x1=-0.5 x2=0.5)
     418|  glorder = degre n du developpement de Gauss-Legendre
     419| OUTPUT:
     420|  x[] = valeur des abscisses ou l'on calcule (dim=n)
     421|  w[] = valeur des coefficients associes
     422| REMARQUES:
     423|  - x et w seront dimensionnes a n.
     424|  - l'integration est rigoureuse si sur l'intervalle d'integration
     425|    la fonction f(x) peut etre approximee par un polynome
     426|    de degre 2*m (monome le + haut x**(2*n-1) )
     427*/
     428void Compute_GaussLeg(unsigned short glorder,vector<double>& x,vector<double>& w,double x1,double x2)
     429{
     430  if(glorder==0) return;
     431  int n = (int)glorder;
     432  x.resize(n,0.); w.resize(n,0.);
     433
     434   int m,j,i;
     435   double z1,z,xm,xl,pp,p3,p2,p1;
     436
     437   m=(n+1)/2;
     438   xm=0.5*(x2+x1);
     439   xl=0.5*(x2-x1);
     440   for (i=1;i<=m;i++)  {
     441      z=cos(M_PI*(i-0.25)/(n+0.5));
     442      do {
     443         p1=1.0;
     444         p2=0.0;
     445         for (j=1;j<=n;j++) {
     446            p3=p2;
     447            p2=p1;
     448            p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
     449         }
     450         pp=n*(z*p1-p2)/(z*z-1.0);
     451         z1=z;
     452         z=z1-p1/pp;
     453      } while (fabs(z-z1) > 3.0e-11);  // epsilon
     454      x[i-1]=xm-xl*z;
     455      x[n-i]=xm+xl*z;
     456      w[i-1]=2.0*xl/((1.0-z*z)*pp*pp);
     457      w[n-i]=w[i-1];
     458   }
     459}

  • trunk/Cosmo/SimLSS/geneutils.h

    r3157 r3196  
    7979
    8080//----------------------------------------------------
     81double InterpTab(double x0,vector<double>& X,vector<double>& Y,unsigned short typint=0);
     82
     83//----------------------------------------------------
    8184int FuncToHisto(GenericFunc& func,Histo& h,bool logaxex=false);
    8285int FuncToVec(GenericFunc& func,TVector<r_8>& h,double xmin,double xmax,bool logaxex=false);
    … …  
    8992unsigned long PoissRandLimit(double mu,double mumax=10.);
    9093
     94//----------------------------------------------------
     95double IntegrateFunc(GenericFunc& func,double xmin,double xmax
     96         ,double perc=0.1,double dxinc=-1.,double dxmax=-1.,unsigned short glorder=4);
     97
     98double IntegrateFuncLog(GenericFunc& func,double lxmin,double lxmax
     99         ,double perc=0.1,double dlxinc=-1.,double dlxmax=-1.,unsigned short glorder=4);
     100
     101void Compute_GaussLeg(unsigned short glorder,vector<double>& x,vector<double>& w,double x1=0.,double x2=1.);
     102
    91103#endif

  • trunk/Cosmo/SimLSS/pkspectrum.cc

    r3193 r3196  
    1111
    1212#include "constcosmo.h"
    13 #include "integfunc.h"
     13#include "geneutils.h"
    1414#include "pkspectrum.h"
    1515
    … …  
    390390void VarianceSpectrum::SetInteg(double dperc,double dlogkinc,double dlogkmax,unsigned short glorder)
    391391// ATTENTION: on n'integre pas f(k)*dk mais k*f(k)*d(log10(k))
    392 // see argument details in function IntegrateFuncLog (integfunc.cc)
     392// see argument details in function IntegrateFuncLog (geneutils.cc)
    393393{
    394394  dperc_ = dperc;  if(dperc_<=0.) dperc_ = 0.1;

  • trunk/Cosmo/SimLSS/schechter.cc

    r3193 r3196  
    6868///////////////////////////////////////////////////////////
    6969
    70 double FluxHI(double m,double d)
     70double Msol2FluxHI(double m,double d)
    7171// Input:
    7272//    m : masse de HI en "Msol"
    … …  
    8888  return  2.0e-28 * m / (d*d);
    8989}
     90
     91double FluxHI2Msol(double f,double d)
     92// Input:
     93//    f : flux total emis a 21 cm en W/m^2
     94//    d : distance en "Mpc"  (si cosmo d=DLum(z))
     95// Return:
     96//    m : masse de HI en "Msol"
     97{
     98  return f *d*d / 2.0e-28;
     99}

  • trunk/Cosmo/SimLSS/schechter.h

    r3193 r3196  
    2525
    2626//-----------------------------------------------------------------------------------
    27 double FluxHI(double m,double d);
     27double Msol2FluxHI(double m,double d);
     28double FluxHI2Msol(double f,double d);
    2829
    2930#endif
Note: See TracChangeset for help on using the changeset viewer.