Changeset 3193 in Sophya


Ignore:
Timestamp:
Mar 21, 2007, 7:18:25 PM (19 years ago)
Author:
cmv
Message:

petis changements cmv 21/03/2007

Location:
trunk/Cosmo/SimLSS
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/Cosmo/SimLSS/cmvconcherr.cc

    r3184 r3193  
    1313
    1414int ObjectType(string nameh,string nameppf);
    15 int ConcatHistoErr(string nameh,vector<string> ppfname,int wrtyp);
     15int ConcatHistoErr(string nameh,vector<string> ppfname,int wrtyp,bool do_cov);
    1616int ConcatHisto2DErr(string nameh,vector<string> ppfname,int wrtyp);
    1717void usage(void);
    … …  
    2121<<"cmvconcherr -w wtyp -n name_histoerr file1.ppf file2.ppf ..."<<endl
    2222<<"  name_histoerr: object name"<<endl
    23 <<"  wtyp: bit0 write ASCII, bit1 write PPF, bit2 write FITS (all=7)"<<endl;
     23<<"  wtyp: bit0 write ASCII, bit1 write PPF, bit2 write FITS (all=7)"<<endl
     24<<"  -C : compute covariance for 1D spectrum"<<endl;
    2425}
    2526
    … …  
    3233 vector<string> ppfname;
    3334 int Write_Type = 3;
     35 bool Do_Cov = false;
    3436
    3537 // --- Decodage arguments
    3638 char c;
    37  while((c = getopt(narg,arg,"hn:w:")) != -1) {
     39 while((c = getopt(narg,arg,"hCn:w:")) != -1) {
    3840  switch (c) {
    3941  case 'n' :
    … …  
    4345    sscanf(optarg,"%d",&Write_Type);
    4446    break;
     47  case 'C' :
     48    Do_Cov = true;
     49    break;
    4550  case 'h' :
    4651  default :
    … …  
    5964 // --- lecture et moyenne des HistoErr ou Histo2DErr
    6065 int nread=0;
    61  if(obtyp==1) nread = ConcatHistoErr(nameh,ppfname,Write_Type);
     66 if(obtyp==1) nread = ConcatHistoErr(nameh,ppfname,Write_Type,Do_Cov);
    6267 else if(obtyp==2) nread = ConcatHisto2DErr(nameh,ppfname,Write_Type);
    6368 else {
    … …  
    95100
    96101//---------------------------------------------------------------
    97 int ConcatHistoErr(string nameh,vector<string> ppfname,int wrtyp)
     102int ConcatHistoErr(string nameh,vector<string> ppfname,int wrtyp,bool do_cov)
    98103{
    99104 if(ppfname.size()<=0) return 0;
    100105
    101106 HistoErr *herrconc = NULL;
     107 TVector<r_8> Tsum; TMatrix<r_8> Tsum2;
    102108 int nherr=0, nread=0, itest=0;
    103109 double sum=0., sum2=0., nsum=0;
    … …  
    111117     nherr = herr.NBins();
    112118     itest = int(herrconc->NBins()/5.+0.5);
     119     if(do_cov) {Tsum.ReSize(nherr); Tsum=0.; Tsum2.ReSize(nherr,nherr); Tsum2=0.;}
    113120   } else if(nherr!=herr.NBins()) {
    114121     cout<<"BAD NUMBER OF BINS"<<endl;
    … …  
    120127   for(int i=0;i<nherr;i++) herrconc->AddBin(i,herr(i));
    121128   sum+=herr(itest); sum2+=herr(itest)*herr(itest); nsum++;
     129
     130   if(do_cov) {
     131     for(int i=0;i<nherr;i++) {
     132       Tsum(i) += herr(i);
     133       for(int j=0;j<nherr;j++) Tsum2(i,j) += herr(i)*herr(j);
     134     }
     135   }
     136
    122137 }
    123138 herrconc->ToVariance();
    … …  
    127142                  <<" sigma^2="<<herrconc->Error2(itest)<<endl;
    128143
     144 // --- Covariance
     145 if(do_cov && nread>1) {
     146   Tsum /= nread;
     147   Tsum2 /= nread;
     148   for(int i=0;i<nherr;i++)
     149     for(int j=0;j<nherr;j++) Tsum2(i,j) -= Tsum(i)*Tsum(j);
     150 }
     151
    129152 // --- ecriture
    130153 if(wrtyp&1) {   // ecriture ASCII
    … …  
    132155   herrconc->WriteASCII(asname);
    133156 }
    134  if(wrtyp&2) {   // ecriture PPF
     157 if(wrtyp&2 || do_cov) {   // ecriture PPF
    135158   string tagobs = "cmvconcherr.ppf";
    136159   POutPersist posobs(tagobs);
    137160   tagobs = "herrconc"; posobs.PutObject(*herrconc,tagobs);
     161   if(do_cov) {
     162     tagobs = "mean"; posobs.PutObject(Tsum,tagobs);
     163     tagobs = "cov"; posobs.PutObject(Tsum2,tagobs);
     164   }
    138165 }
    139166 if(wrtyp&4) {   // ecriture FITS
    … …  
    220247n/plot herrconc.sqrt(err2)/val%log10(x) x>0&&err2>0&&val>0 ! "connectpoints"
    221248
     249disp mean
     250imag cov
     251
     252del cor
     253c++exec \
     254TMatrix<r_8> cor(cov,false); cor = 0.; \
     255for(int i=0;i<cor.NRows();i++) { \
     256  for(int j=0;j<cor.NCols();j++) { \
     257    double v = cov(i,i)*cov(j,j); \
     258    if(v<=0.) continue; \
     259    cor(i,j) = cov(i,j)/sqrt(v); \
     260} } \
     261KeepObj(cor); \
     262cout<<"Matrice cor cree "<<endl;
     263
     264imag cor
     265
    222266#### Histo2DErr 2D
    223267openppf cmvconcherr2.ppf

  • trunk/Cosmo/SimLSS/cmvdefsurv.cc

    r3115 r3193  
    1515#include "planckspectra.h"
    1616
     17/* --- Check Peterson at al. astro-ph/0606104 v1
     18cmvdefsurv -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
     19 --- */
     20
    1721inline double rad2deg(double trad) {return trad/M_PI*180.;}
    1822inline double rad2min(double trad) {return trad/M_PI*180.*60.;}
    … …  
    2832      <<" -y adty,atylarg : idem selon y, si <=0 meme que x"<<endl
    2933      <<" -z dred,redlarg : resolution en redshift, largeur en redshift"<<endl
    30       <<" -P : on donne -x -y -z en Mpc aulieu d\'angles et de redshift"<<endl
     34      <<" -P : on donne -x -y -z en Mpc au lieu d\'angles et de redshift"<<endl
     35      <<" -L lobewidth : taille du lobe d\'observation (FWHM) en arcmin (def= 1\')"<<endl
    3136      <<" -O surf,tobs : surface effective (m^2) et temps d\'observation (s)"<<endl
    32       <<" -A Tbruit : temperature de bruit (K)"<<endl
     37      <<" -A Tsys : temperature du system (K)"<<endl
    3338      <<" -S Tsynch,indnu : temperature (K) synch a 408 Mhz, index d\'evolution"<<endl
    3439      <<"                   (indnu==0 no evolution with freq.)"<<endl
    35       <<" -M  : masse de HI de reference (MSol), si <=0 mean schechter"<<endl
     40      <<" -M  : masse de HI de reference (MSol), si <=0 mean schechter in pixel"<<endl
    3641      <<" -F  : HI flux factor to be applied for our redshift"<<endl
    37       <<" -1 : on ne mesure qu\'une polarisation"<<endl
     42      <<" -V Vrot : largeur en vitesse (km/s) pour l\'elargissement doppler (def=300km/s)"<<endl
     43      <<" -U h100,om0,ol0,w0,or0,flat : cosmology"<<endl
     44      <<" -2 : two polarisations measured"<<endl
    3845      <<" redshift : redshift moyen du survey"<<endl
    3946      <<endl;
    … …  
    4552 // WMAP
    4653 unsigned short flat = 0;
    47  double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
    48  cout<<"\nh100="<<h100<<" Om0="<<om0<<" Or0="<<or0<<" Or0="<<or0<<" Ol0="<<ol0<<" w0="<<w0<<endl;
     54 double h100=0.71, om0=0.267804, or0=7.9e-05*0., ol0=0.73,w0=-1.;
    4955 // Schechter
    5056 double h75 = h100 / 0.75;
    … …  
    6167 int nx,ny,nz;
    6268 double redshift = 1., dred=0.01, redlarg=0.3;
    63  double tobs = 3600., surfeff = 400000.;
    64  double Tbruit=75.;
     69 double tobs = 6000., surfeff = 400000.;
     70 double lobewidth = 1.;  // taille du lobe d'observation en arcmin
     71 double Tsys=75.;
    6572 // a 408 MHz (Haslam) + evol index a -2.6
    6673 double Tsynch408=60., nuhaslam=0.408, indnu = -2.6;
    6774 double mhiref = -1.; // reference Mass en HI (def integ schechter)
    6875 double hifactor = 1.;
    69  double facpolar = 1.; // si on ne mesure qu'un polarisation 0.5
     76 double vrot = 300.; // largeur en vitesse (km/s) pour elargissement doppler
     77 double facpolar = 0.5; // si on ne mesure les 2 polars -> 1.0
    7078
    7179 // --- Decodage arguments
    7280 char c;
    73   while((c = getopt(narg,arg,"hP1x:y:z:A:S:O:M:F:")) != -1) {
     81  while((c = getopt(narg,arg,"hP2x:y:z:A:S:O:M:F:V:U:L:")) != -1) {
    7482  switch (c) {
    7583  case 'P' :
    … …  
    8896    sscanf(optarg,"%lf,%lf",&surfeff,&tobs);
    8997    break;
     98  case 'L' :
     99    sscanf(optarg,"%lf",&lobewidth);
     100    break;
    90101  case 'A' :
    91     sscanf(optarg,"%lf",&Tbruit);
     102    sscanf(optarg,"%lf",&Tsys);
    92103    break;
    93104  case 'S' :
    … …  
    100111    sscanf(optarg,"%lf",&hifactor);
    101112    break;
    102   case '1' :
    103     facpolar = 0.5;
     113  case 'V' :
     114    sscanf(optarg,"%lf",&vrot);
     115    break;
     116  case 'U' :
     117    sscanf(optarg,"%lf,%lf,%lf,%lf,%u",&h100,&om0,&ol0,&w0,&flat);
     118    break;
     119  case '2' :
     120    facpolar = 1.0;
    104121    break;
    105122  case 'h' :
    … …  
    107124    usage(); return -1;
    108125  }
    109  }
     126 } 
    110127 if(optind>=narg) {usage(); return-1;}
    111128 sscanf(arg[optind],"%lf",&redshift);
    … …  
    113130
    114131 // --- Initialisation de la Cosmologie
     132 cout<<"\nh100="<<h100<<" Om0="<<om0<<" Or0="<<or0<<" Or0="
     133     <<or0<<" Ol0="<<ol0<<" w0="<<w0<<" flat="<<flat<<endl;
    115134 cout<<"\n--- Cosmology for z = "<<redshift<<endl;
    116135 CosmoCalc univ(flat,true,2.*redshift);
    … …  
    118137 univ.SetInteg(perc,dzinc,dzmax,glorder);
    119138 univ.SetDynParam(h100,om0,or0,ol0,w0);
     139 univ.Print(0.);
    120140 univ.Print(redshift);
    121141
    … …  
    258278 cout<<"dang lumi = "<<dlum<<" in ["<<dlumlim[0]<<","<<dlumlim[1]<<"] Mpc"<<endl;
    259279
     280 double slobe = lobewidth/2.35482; // sigma du lobe en arcmin
     281 double lobecyl = sqrt(8.)*slobe; // diametre du lobe cylindrique equiv en arcmin
     282 double lobearea = M_PI*lobecyl*lobecyl/4.;  // en arcmin^2 (hypothese lobe gaussien)
     283 cout<<"\nBeam FWHM = "<<lobewidth<<"\' -> sigma = "<<slobe<<"\' -> "
     284     <<" Dcyl = "<<lobecyl<<"\' -> area = "<<lobearea<<" arcmin^2"<<endl;
     285 double nlobes = rad2min(adtx)*rad2min(adty)/lobearea;
     286 cout<<"Number of beams in one transversal pixel = "<<nlobes<<endl;
     287
    260288 // --- Power emitted by HI
    261289 cout<<"\n--- Power from HI for M = "<<mhiref<<" Msol at "<<nuhiz<<" GHz"<<endl;
    … …  
    267295     <<"      in ["<<hifactor*FluxHI(mhiref,dlumlim[0])
    268296     <<","<<hifactor*FluxHI(mhiref,dlumlim[1])<<"] W/m^2"<<endl;
    269  double sfhi = fhi / (dnuhiz*1.e+9) / Jansky2Watt_cst;
     297 double sfhi = fhi / (dnuhiz*1e9) / Jansky2Watt_cst;
    270298 cout<<"If spread over "<<dnuhiz<<" GHz, flux density = "<<sfhi<<" Jy"<<endl;
    271299
    272300 // --- Signal analysis
    273301 cout<<"\n--- Signal analysis"<<endl;
     302 cout<<"Facteur polar = "<<facpolar<<endl;
     303
    274304 PlanckSpectra planck(T_CMB_Par);
    275305 planck.SetApprox(1);  // Rayleigh spectra
    … …  
    278308 planck.SetTypSpectra(0); // radiance W/m^2/Sr/Hz
    279309
    280  double hnu = h_Planck_Cst*(nuhiz*1.e+9);
    281  cout<<"Facteur polar = "<<facpolar<<",   h.nu="<<hnu<<" J"<<endl;
    282 
     310 // Signal
    283311 double psig = facpolar * fhi * surfeff;
    284  double esig = psig * tobs;
    285  double nsig = esig / hnu;
    286  double tsig = psig / k_Boltzman_Cst / (dnuhiz*1.e+9);
    287  double ssig = psig / surfeff / (dnuhiz*1.e+9) / Jansky2Watt_cst;
    288  cout<<"Signal("<<mhiref<<" Msol): P="<<psig<<" W -> E="<<esig<<" J -> "<<nsig<<" ph"<<endl;
     312 double tsig = psig / k_Boltzman_Cst / (dnuhiz*1e9);
     313 double ssig = psig / surfeff / (dnuhiz*1e9) / Jansky2Watt_cst;
     314 cout<<"Signal("<<mhiref<<" Msol): P="<<psig<<" W"<<endl;
     315 cout<<"     flux density = "<<ssig<<" Jy  (for Dnu="<<dnuhiz<<" GHz)"<<endl;
    289316 cout<<"     Antenna temperature: tsig="<<tsig<<" K"<<endl;
    290  cout<<"     flux density = "<<ssig<<" Jy"<<endl;
    291 
    292  double facnoise = facpolar * tobs * surfeff * angsol * (dnuhiz*1.e+9);
    293 
     317
     318 // Elargissement doppler de la raie a 21cm: dNu = vrot/C * Nu(21cm) / (1+z)
     319 double doplarge = vrot / SpeedOfLight_Cst * nuhiz;
     320 cout<<"     Doppler width="<<doplarge*1.e3<<" MHz  for rotation width of "<<vrot<<" km/s"<<endl;
     321 if(doplarge>dnuhiz) {
     322   cout<<"Warning: doppler width "<<doplarge<<" GHz > "<<dnuhiz<<" GHz redshift bin width"<<endl;
     323 }
     324
     325 // Synchrotron
    294326 double tsynch = Tsynch408;
    295327        if(fabs(indnu)>1.e-50) tsynch *= pow(nuhiz/nuhaslam,indnu);
    296328 planck.SetTemperature(tsynch);
    297  double psynch = facpolar * planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1.e+9);
    298  double esynch = psynch * tobs;
    299  double nsynch = esynch / hnu;
    300  double ssynch = psynch / surfeff / (dnuhiz*1.e+9) / Jansky2Watt_cst;
     329 double psynch = facpolar * planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1e9);
     330 double ssynch = psynch / surfeff / (dnuhiz*1e9) / Jansky2Watt_cst;
    301331 cout<<"Synchrotron: T="<<Tsynch408<<" K ("<<nuhaslam<<" GHz), "
    302332     <<tsynch<<" K ("<<nuhiz<<" GHz)"<<endl
    303      <<"             P="<<psynch<<" W -> E="<<esynch<<" J -> "<<nsynch<<" ph"<<endl;
    304  cout<<"     flux density = "<<ssynch<<" Jy"<<endl;
    305 
     333     <<"             P="<<psynch<<" W"<<endl;
     334 cout<<"             flux density = "<<ssynch<<" Jy"<<endl;
     335
     336 // CMB
    306337 double tcmb = T_CMB_Par;
    307338 planck.SetTemperature(tcmb);
    308  double pcmb = facpolar * planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1.e+9);
    309  double ecmb = pcmb * tobs;
    310  double ncmb = ecmb / hnu;
     339 double pcmb = facpolar * planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1e9);
    311340 double scmb = pcmb / surfeff / (dnuhiz*1.e+9) / Jansky2Watt_cst;
    312  cout<<"CMB: T="<<tcmb<<" K -> P="<<pcmb<<" W -> E="<<ecmb<<" J -> "<<ncmb<<" ph"<<endl;
    313  cout<<"     flux density = "<<scmb<<" Jy"<<endl;
     341 cout<<"CMB: T="<<tcmb<<" K -> P="<<pcmb<<" W"<<endl;
     342 cout<<"                       flux density = "<<scmb<<" Jy"<<endl;
    314343
    315344 // --- Noise analysis
    316345 cout<<"\n--- Noise analysis"<<endl;
    317  double pbruit = k_Boltzman_Cst * Tbruit * (dnuhiz*1.e+9);
    318  double ebruit = pbruit * tobs;
    319  cout<<"Bruit: T="<<Tbruit<<" K, P="<<pbruit<<" W -> E="<<ebruit<<" J"<<endl;
    320 
    321  double seq = (1./facpolar) * k_Boltzman_Cst * Tbruit / surfeff /Jansky2Watt_cst;
    322  cout<<"\nSystem equivalent flux density: "<<seq<<" Jy"<<endl;
    323 
    324  double slim = seq / sqrt(2.*(dnuhiz*1.e+9)*tobs);
    325  cout<<"Observation flux density limit: "<<slim<<" Jy"<<endl;
     346 double psys = k_Boltzman_Cst * Tsys * (dnuhiz*1.e+9);
     347 cout<<"Noise: T="<<Tsys<<" K, P="<<psys<<" W  (for Dnu="<<dnuhiz<<" GHz)"<<endl;
     348
     349 double slim = 2. * k_Boltzman_Cst * Tsys / surfeff
     350               / sqrt(2.*(dnuhiz*1.e+9)*tobs) /Jansky2Watt_cst;
     351 cout<<"Observation flux density limit: "<<slim<<" Jy (in 1 lobe)"<<endl;
     352
     353 slim = slim * sqrt(nlobes);
     354 cout<<"Observation flux density limit: "<<slim<<" Jy (in "<<nlobes<<" lobes)"<<endl;
    326355
    327356 double SsN = ssig/slim;

  • trunk/Cosmo/SimLSS/cmvobserv3d.cc

    r3182 r3193  
    3131     <<" -s snoise : sigma du bruit en Msol"<<endl
    3232     <<" -2 : compute 2D spectrum"<<endl
    33      <<" -F : write cube in FITS format"<<endl
     33     <<" -M schmin,schmax,nsch : min,max mass and nb points for schechter HI"<<endl
     34     <<" -W : write cube in FITS format"<<endl
    3435     <<" -P : write cube in PPF format"<<endl
    3536     <<" -V : compute variance from real space"<<endl
    … …  
    6566
    6667 // *** Schechter mass function definition
    67  double h75 = 0.71 / 0.75;
    68  double nstar = 0.006*pow(h75,3.); 
     68 double h75 = h100 / 0.75;
     69 double nstar = 0.006*pow(h75,3.);
    6970 double mstar = pow(10.,9.8/(h75*h75));  // MSol
    7071 double alpha = -1.37;
    … …  
    7273 double schmin=1e8, schmax=1e12;
    7374 int schnpt = 1000;
    74  double lschmin=log10(schmin), lschmax=log10(schmax), dlsch=(lschmax-lschmin)/schnpt;
    7575
    7676 // *** Niveau de bruit
    … …  
    9191
    9292 char c;
    93  while((c = getopt(narg,arg,"ha0PWV2Gx:y:z:s:Z:")) != -1) {
     93 while((c = getopt(narg,arg,"ha0PWV2Gx:y:z:s:Z:M:")) != -1) {
    9494  switch (c) {
    9595  case 'a' :
    … …  
    119119  case '2' :
    120120    comp2dspec = true;
     121    break;
     122  case 'M' :
     123    sscanf(optarg,"%lf,%lf,%d",&schmin,&schmax,&schnpt);
    121124    break;
    122125  case 'V' :
    … …  
    134137  }
    135138 }
     139
     140 double lschmin=log10(schmin), lschmax=log10(schmax), dlsch=(lschmax-lschmin)/schnpt;
    136141
    137142 string tagobs = "cmvobserv3d.ppf";
    … …  
    157162 univ.Print();
    158163 double loscomref = univ.Dloscom(zref);
    159  cout<<"zref = "<<zref<<" -> dloscom = "<<loscomref<<" Mpc"<<endl;
    160 
    161  //-----------------------------------------------------------------
    162   cout<<endl<<"\n--- Create Spectrum and mass function"<<endl;
     164 cout<<"\nzref = "<<zref<<" -> dloscom = "<<loscomref<<" Mpc"<<endl;
     165 univ.Print(zref);
     166
     167 //-----------------------------------------------------------------
     168 cout<<endl<<"\n--- Create Spectrum"<<endl;
    163169
    164170 InitialSpectrum pkini(ns,as);
    … …  
    168174
    169175 GrowthFactor growth(om0,ol0);
     176 // GrowthFactor growth(1.,0.); // D(z) = 1/(1+z)
    170177
    171178 PkSpectrum0 pk0(pkini,tf);
    … …  
    173180 PkSpectrumZ pkz(pk0,growth,zref);
    174181 
     182 //-----------------------------------------------------------------
     183 cout<<endl<<"\n--- Create mass function"<<endl;
     184
    175185 Schechter sch(nstar,mstar,alpha);
     186 sch.Print();
    176187
    177188 //-----------------------------------------------------------------
    … …  
    234245 double mass_by_mpc3 = IntegrateFuncLog(sch,lschmin,lschmax,0.01,dlsch,10.*dlsch,4);
    235246 cout<<"Galaxy mass density= "<<mass_by_mpc3<<" Msol/Mpc^3 between limits"<<endl;
     247 cout<<"Omega_HI at z=0 is "<<mass_by_mpc3/(univ.Rhoc(0.)*GCm3toMsolMpc3_Cst)<<endl
     248     <<"         at z="<<zref<<" is "<<mass_by_mpc3/(univ.Rhoc(zref)*GCm3toMsolMpc3_Cst)<<endl;
    236249
    237250 PrtTim(">>>> End of definition");
    … …  
    463476 }
    464477
    465  cout<<"\n--- Add noise to HI Flux snoise="<<snoise<<endl;
    466  fluct3d.AddNoise2Real(snoise);
    467  nm = fluct3d.MeanSigma2(rm,rs2);
    468  cout<<"HI mass with noise: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
    469      <<rs2<<" -> "<<sqrt(rs2)<<endl;
    470  PrtTim(">>>> End Add noise");
     478 if(snoise>0.) {
     479   cout<<"\n--- Add noise to HI Flux snoise="<<snoise<<endl;
     480   fluct3d.AddNoise2Real(snoise);
     481   nm = fluct3d.MeanSigma2(rm,rs2);
     482   cout<<"HI mass with noise: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
     483       <<rs2<<" -> "<<sqrt(rs2)<<endl;
     484   PrtTim(">>>> End Add noise");
     485 }
    471486
    472487 //-----------------------------------------------------------------

  • trunk/Cosmo/SimLSS/cmvtstpk.cc

    r3120 r3193  
    2222     <<" -H h100 -B Ob0 -M Om0 -L Ol0,w0"<<endl
    2323     <<" -k npt,lkmin,lkmax : les valeurs sont en log10"<<endl
    24      <<" -s scale : on multiplie pkz par scale"<<endl;
     24     <<" -s scale : on multiplie pkz par scale"<<endl
     25     <<" -w : write spectra on ASCII file"<<endl;
    2526}
    2627
    … …  
    4142 double lkmin = -3., lkmax=2.;
    4243 double scale = 1.;
     44 bool wrascii = false;
    4345
    4446 char c;
    45   while((c = getopt(narg,arg,"hk:s:H:M:B:L:")) != -1) {
     47  while((c = getopt(narg,arg,"hwk:s:H:M:B:L:")) != -1) {
    4648  switch (c) {
    4749  case 'H' :
    … …  
    6466  case 's' :
    6567    sscanf(optarg,"%lf",&scale);
     68    break;
     69  case 'w' :
     70    wrascii = true;
    6671    break;
    6772  case 'h' :
    … …  
    150155
    151156 //--------------------------
    152  cout<<">>>> ASCII"<<endl;
    153  if(1) {
     157 if(wrascii) {
     158   cout<<">>>> ASCII"<<endl;
    154159   FILE * fdata = fopen("cmvtstpk.data","w");
    155160   fprintf(fdata,"# z= %g : k(Mpc^-1) pkini, tf pk0 pkz, tfnosc2 pk0nosc2 pkznosc2, ",zval);

  • trunk/Cosmo/SimLSS/cmvtuniv.cc

    r3115 r3193  
    5757 univ1.Print();
    5858
    59  //tstprint(univ,z1,z2,dz);
     59 tstprint(univ,z1,z2,dz);
    6060 //tstprint(univ,univ1,z1,z2,dz);
    6161 //tstspeed(univ,z1,z2,dz);
    6262 //tstntuple(univ,z1,z2,dz);
    63  tstinterp(univ,z1,z2,dz);
     63 //tstinterp(univ,z1,z2,dz);
    6464
    6565 return 0;

  • trunk/Cosmo/SimLSS/constcosmo.h

    r3120 r3193  
    22#define CONSTCOSMO_SEEN
    33
    4 static const double MpctoMeters_Cst = 3.0856e+22; // Mpc in meters
    5 static const double YeartoSec_Cst = 31556925.2; // seconds in a year
    64static const double SpeedOfLight_Cst = 299792.458; // speed of light "km/sec"
    7 static const double LightYear_Cst = 9.46073042e+15; // 1 Light-Year in "m"
    85static const double G_Newton_Cst = 6.6742e-11; // G_N in SI units "m^3 / kg /s^2"
    96static const double h_Planck_Cst = 6.6260693e-34;   // Planck constant SI units "J s"
    107static const double k_Boltzman_Cst = 1.3806503e-23;  // Boltzman constant SI "J / K"
    11 static const double ProtonMass_Cst = 1.67262171e-24; // Proton mass in "g"
    12 static const double SolarMass_Cst = 1.98844e+33; // Solar mass in "g"
    138//            Sigma = 2*PI^5*K^4/(15.*C^2*H^3) = PI^2*K^4/(60.*Hbar^3*C2)
    149static const double StefBoltz_Cst = 5.670400e-8; // constante de Stefan-Boltzman in "W/m^2/K^4"
    1510static const double Zeta_3_Cst = 1.2020569031595942854; // Zeta(3)
    1611
    17 static const double Jansky2Watt_cst = 1.e-26;  // conversion Jansky to "Watt/m^2/Hz == J/m^2"
    18 static const double Deg2Rad_Cst = (M_PI/180.); // conversion deg to rad
     12static const double YeartoSec_Cst = 31556925.2; // 1 year into seconds
     13static const double MpctoMeters_Cst = 3.0856e+22; // 1 Mpc into meters
     14static const double LightYear_Cst = 9.46073042e+15; // 1 Light-Year into "m"
     15static const double Jansky2Watt_cst = 1.e-26;  // 1 Jansky into "Watt/m^2/Hz == J/m^2"
     16static const double SolarMass_Cst = 1.98844e+33; // 1 Solar mass into "g"
     17//                  GCm3toMsolMpc3_Cst = (MpctoMeters_Cst*100)^3 / SolarMass_Cst
     18static const double GCm3toMsolMpc3_Cst = 1.477428208143872e+40; // 1 g/cm^3 into 1 Msol/Mpc^3
     19static const double ProtonMass_Cst = 1.67262171e-24; // 1 "proton mass" in "g"
     20static const double Deg2Rad_Cst = (M_PI/180.); // 1 deg into rad
    1921
    2022static const double T_CMB_Par = 2.725;  // temperature du CMB en "K"
    2123static const double T_NU_Par = 1.9;  // temperature des neutrinos en "K"
    22 static const double Fr_HyperFin_Par = 1.4204;  // frequence de la raie hyperfine HI en "GHz"
     24static const double Fr_HyperFin_Par = 1.420405751786;  // frequence de la raie hyperfine HI en "GHz"
    2325
    2426#endif

  • trunk/Cosmo/SimLSS/cosmocalc.cc

    r3115 r3193  
    206206
    207207  printf("CosmoCalc::Print(spl=%d,zmax=%g,flat=%u)  for z=%g\n",int(_usespline),ZMax(),Flat(),z);
    208   printf("h100=%g H0=%g Dhub=%g  H(z)=%g  Rhoc=%g g/cm^3\n"
    209         ,h100(),H0(),Dhubble(),H(z),Rhoc(z));
     208  printf("h100=%g H0=%g Dhub=%g  H(z)=%g  Rhoc=%g g/cm^3 =%g Msol/Mpc^3\n"
     209        ,h100(),H0(),Dhubble(),H(z),Rhoc(z),Rhoc(z)*GCm3toMsolMpc3_Cst);
    210210  printf("Olambda=%g  W0=%g\n",Olambda(z),_W0);
    211211  printf("Omatter=%g  Obaryon=%g\n",Omatter(z),Obaryon(z));

  • trunk/Cosmo/SimLSS/pkspectrum.cc

    r3115 r3193  
    244244
    245245// From Eisenstein & Hu ApJ 496:605-614 1998 April 1
     246// Pour avoir D(z) = 1/(1+z) faire: OmegaMatter0=1 OmegaLambda0=0
    246247GrowthFactor::GrowthFactor(double OmegaMatter0,double OmegaLambda0)
    247248  : O0_(OmegaMatter0) , Ol_(OmegaLambda0) , Ok_(1.-OmegaMatter0-OmegaLambda0)

  • trunk/Cosmo/SimLSS/schechter.cc

    r3115 r3193  
    5454}
    5555
     56void Schechter::Print(void)
     57{
     58  cout<<"Schechter::Print: nstar="<<nstar_<<" Mpc^-3"
     59      <<"  mstar="<<mstar_<<" MSol"
     60      <<"  alpha="<<alpha_
     61      <<"  (outvalue="<<outvalue_<<" -> return ";
     62  if(outvalue_) cout<<"m*dn/dm)"; else cout<<"dn/dm)";
     63  cout<<endl;
     64}
     65
    5666///////////////////////////////////////////////////////////
    5767//******************* Le Flux a 21cm ********************//
    … …  
    6474// Return:
    6575//    le flux total emis a 21 cm en W/m^2
    66 // Ref: Binney & Merrifield, Galactic Astronomy p474 (ed:1998)
    67 //      J.Peterson & K.Bandura, astroph-0606104  (eq 7)
    68 //      F.Abdalla & Rawlings, astroph-0411342 (eq 7 may pb de d'unites?)
     76// Ref:
     77// -- Binney & Merrifield, Galactic Astronomy p474 (ed:1998)
     78// S(W/m^2) = 1e-26 * Nu_21cm(Hz) * m(en masse solaire) /(2.35e5 * C(km/s) * Dlum^2)
     79//          = 2.0e-28 * m / Dlum^2
     80// -- J.Peterson & K.Bandura, astroph-0606104  (eq 7)
     81//    F.Abdalla & Rawlings, astroph-0411342 (eq 7 mais pb de d'unites?)
     82//          (A_21cm = 2.86888e-15 s^-1)
     83// S(W/m^2) = 3/4 * h(J.s) * Nu_21cm(Hz) * A_21cm(s^-1) * Msol(kg)/m_proton(kg)
     84//                         / (4 Pi) / (Dlum(Mpc) * 3.0857e22(m/Mpc))^2
     85//          = 2.0e-28 * m / Dlum^2
     86//-------------
    6987{
    70   return  m / (2.35e+5 * d*d )
    71           * Jansky2Watt_cst
    72           * (Fr_HyperFin_Par*1e+9)/SpeedOfLight_Cst;
     88  return  2.0e-28 * m / (d*d);
    7389}
    74 

  • trunk/Cosmo/SimLSS/schechter.h

    r3115 r3193  
    1616  void SetOutValue(unsigned short outvalue=0);
    1717  virtual double operator() (double m);
     18  virtual void Print(void);
    1819protected:
    1920  double nstar_,mstar_,alpha_;
Note: See TracChangeset for help on using the changeset viewer.