Index: /trunk/Cosmo/SimLSS/cmvconcherr.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvconcherr.cc	(revision 3192)
+++ /trunk/Cosmo/SimLSS/cmvconcherr.cc	(revision 3193)
@@ -13,5 +13,5 @@
 
 int ObjectType(string nameh,string nameppf);
-int ConcatHistoErr(string nameh,vector<string> ppfname,int wrtyp);
+int ConcatHistoErr(string nameh,vector<string> ppfname,int wrtyp,bool do_cov);
 int ConcatHisto2DErr(string nameh,vector<string> ppfname,int wrtyp);
 void usage(void);
@@ -21,5 +21,6 @@
 <<"cmvconcherr -w wtyp -n name_histoerr file1.ppf file2.ppf ..."<<endl
 <<"  name_histoerr: object name"<<endl
-<<"  wtyp: bit0 write ASCII, bit1 write PPF, bit2 write FITS (all=7)"<<endl;
+<<"  wtyp: bit0 write ASCII, bit1 write PPF, bit2 write FITS (all=7)"<<endl
+<<"  -C : compute covariance for 1D spectrum"<<endl;
 }
 
@@ -32,8 +33,9 @@
  vector<string> ppfname;
  int Write_Type = 3;
+ bool Do_Cov = false;
 
  // --- Decodage arguments
  char c;
- while((c = getopt(narg,arg,"hn:w:")) != -1) {
+ while((c = getopt(narg,arg,"hCn:w:")) != -1) {
   switch (c) {
   case 'n' :
@@ -43,4 +45,7 @@
     sscanf(optarg,"%d",&Write_Type);
     break;
+  case 'C' :
+    Do_Cov = true;
+    break;
   case 'h' :
   default :
@@ -59,5 +64,5 @@
  // --- lecture et moyenne des HistoErr ou Histo2DErr
  int nread=0;
- if(obtyp==1) nread = ConcatHistoErr(nameh,ppfname,Write_Type);
+ if(obtyp==1) nread = ConcatHistoErr(nameh,ppfname,Write_Type,Do_Cov);
  else if(obtyp==2) nread = ConcatHisto2DErr(nameh,ppfname,Write_Type);
  else {
@@ -95,9 +100,10 @@
 
 //---------------------------------------------------------------
-int ConcatHistoErr(string nameh,vector<string> ppfname,int wrtyp)
+int ConcatHistoErr(string nameh,vector<string> ppfname,int wrtyp,bool do_cov)
 {
  if(ppfname.size()<=0) return 0;
 
  HistoErr *herrconc = NULL;
+ TVector<r_8> Tsum; TMatrix<r_8> Tsum2;
  int nherr=0, nread=0, itest=0;
  double sum=0., sum2=0., nsum=0;
@@ -111,4 +117,5 @@
      nherr = herr.NBins();
      itest = int(herrconc->NBins()/5.+0.5);
+     if(do_cov) {Tsum.ReSize(nherr); Tsum=0.; Tsum2.ReSize(nherr,nherr); Tsum2=0.;}
    } else if(nherr!=herr.NBins()) {
      cout<<"BAD NUMBER OF BINS"<<endl;
@@ -120,4 +127,12 @@
    for(int i=0;i<nherr;i++) herrconc->AddBin(i,herr(i));
    sum+=herr(itest); sum2+=herr(itest)*herr(itest); nsum++;
+
+   if(do_cov) {
+     for(int i=0;i<nherr;i++) {
+       Tsum(i) += herr(i);
+       for(int j=0;j<nherr;j++) Tsum2(i,j) += herr(i)*herr(j);
+     }
+   }
+
  }
  herrconc->ToVariance();
@@ -127,4 +142,12 @@
                   <<" sigma^2="<<herrconc->Error2(itest)<<endl;
 
+ // --- Covariance
+ if(do_cov && nread>1) {
+   Tsum /= nread;
+   Tsum2 /= nread;
+   for(int i=0;i<nherr;i++)
+     for(int j=0;j<nherr;j++) Tsum2(i,j) -= Tsum(i)*Tsum(j);
+ }
+
  // --- ecriture
  if(wrtyp&1) {   // ecriture ASCII
@@ -132,8 +155,12 @@
    herrconc->WriteASCII(asname);
  }
- if(wrtyp&2) {   // ecriture PPF
+ if(wrtyp&2 || do_cov) {   // ecriture PPF
    string tagobs = "cmvconcherr.ppf";
    POutPersist posobs(tagobs);
    tagobs = "herrconc"; posobs.PutObject(*herrconc,tagobs);
+   if(do_cov) {
+     tagobs = "mean"; posobs.PutObject(Tsum,tagobs);
+     tagobs = "cov"; posobs.PutObject(Tsum2,tagobs);
+   }
  }
  if(wrtyp&4) {   // ecriture FITS
@@ -220,4 +247,21 @@
 n/plot herrconc.sqrt(err2)/val%log10(x) x>0&&err2>0&&val>0 ! "connectpoints"
 
+disp mean
+imag cov
+
+del cor
+c++exec \
+TMatrix<r_8> cor(cov,false); cor = 0.; \
+for(int i=0;i<cor.NRows();i++) { \
+  for(int j=0;j<cor.NCols();j++) { \
+    double v = cov(i,i)*cov(j,j); \
+    if(v<=0.) continue; \
+    cor(i,j) = cov(i,j)/sqrt(v); \
+} } \
+KeepObj(cor); \
+cout<<"Matrice cor cree "<<endl;
+
+imag cor
+
 #### Histo2DErr 2D
 openppf cmvconcherr2.ppf
Index: /trunk/Cosmo/SimLSS/cmvdefsurv.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvdefsurv.cc	(revision 3192)
+++ /trunk/Cosmo/SimLSS/cmvdefsurv.cc	(revision 3193)
@@ -15,4 +15,8 @@
 #include "planckspectra.h"
 
+/* --- Check Peterson at al. astro-ph/0606104 v1
+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
+ --- */
+
 inline double rad2deg(double trad) {return trad/M_PI*180.;}
 inline double rad2min(double trad) {return trad/M_PI*180.*60.;}
@@ -28,12 +32,15 @@
       <<" -y adty,atylarg : idem selon y, si <=0 meme que x"<<endl
       <<" -z dred,redlarg : resolution en redshift, largeur en redshift"<<endl
-      <<" -P : on donne -x -y -z en Mpc aulieu d\'angles et de redshift"<<endl
+      <<" -P : on donne -x -y -z en Mpc au lieu d\'angles et de redshift"<<endl
+      <<" -L lobewidth : taille du lobe d\'observation (FWHM) en arcmin (def= 1\')"<<endl
       <<" -O surf,tobs : surface effective (m^2) et temps d\'observation (s)"<<endl
-      <<" -A Tbruit : temperature de bruit (K)"<<endl
+      <<" -A Tsys : temperature du system (K)"<<endl
       <<" -S Tsynch,indnu : temperature (K) synch a 408 Mhz, index d\'evolution"<<endl
       <<"                   (indnu==0 no evolution with freq.)"<<endl
-      <<" -M  : masse de HI de reference (MSol), si <=0 mean schechter"<<endl
+      <<" -M  : masse de HI de reference (MSol), si <=0 mean schechter in pixel"<<endl
       <<" -F  : HI flux factor to be applied for our redshift"<<endl
-      <<" -1 : on ne mesure qu\'une polarisation"<<endl
+      <<" -V Vrot : largeur en vitesse (km/s) pour l\'elargissement doppler (def=300km/s)"<<endl
+      <<" -U h100,om0,ol0,w0,or0,flat : cosmology"<<endl
+      <<" -2 : two polarisations measured"<<endl
       <<" redshift : redshift moyen du survey"<<endl
       <<endl;
@@ -45,6 +52,5 @@
  // WMAP
  unsigned short flat = 0;
- double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
- cout<<"\nh100="<<h100<<" Om0="<<om0<<" Or0="<<or0<<" Or0="<<or0<<" Ol0="<<ol0<<" w0="<<w0<<endl;
+ double h100=0.71, om0=0.267804, or0=7.9e-05*0., ol0=0.73,w0=-1.;
  // Schechter
  double h75 = h100 / 0.75;
@@ -61,15 +67,17 @@
  int nx,ny,nz;
  double redshift = 1., dred=0.01, redlarg=0.3;
- double tobs = 3600., surfeff = 400000.;
- double Tbruit=75.;
+ double tobs = 6000., surfeff = 400000.;
+ double lobewidth = 1.;  // taille du lobe d'observation en arcmin
+ double Tsys=75.;
  // a 408 MHz (Haslam) + evol index a -2.6
  double Tsynch408=60., nuhaslam=0.408, indnu = -2.6;
  double mhiref = -1.; // reference Mass en HI (def integ schechter)
  double hifactor = 1.;
- double facpolar = 1.; // si on ne mesure qu'un polarisation 0.5
+ double vrot = 300.; // largeur en vitesse (km/s) pour elargissement doppler
+ double facpolar = 0.5; // si on ne mesure les 2 polars -> 1.0
 
  // --- Decodage arguments
  char c;
-  while((c = getopt(narg,arg,"hP1x:y:z:A:S:O:M:F:")) != -1) {
+  while((c = getopt(narg,arg,"hP2x:y:z:A:S:O:M:F:V:U:L:")) != -1) {
   switch (c) {
   case 'P' :
@@ -88,6 +96,9 @@
     sscanf(optarg,"%lf,%lf",&surfeff,&tobs);
     break;
+  case 'L' :
+    sscanf(optarg,"%lf",&lobewidth);
+    break;
   case 'A' :
-    sscanf(optarg,"%lf",&Tbruit);
+    sscanf(optarg,"%lf",&Tsys);
     break;
   case 'S' :
@@ -100,6 +111,12 @@
     sscanf(optarg,"%lf",&hifactor);
     break;
-  case '1' :
-    facpolar = 0.5;
+  case 'V' :
+    sscanf(optarg,"%lf",&vrot);
+    break;
+  case 'U' :
+    sscanf(optarg,"%lf,%lf,%lf,%lf,%u",&h100,&om0,&ol0,&w0,&flat);
+    break;
+  case '2' :
+    facpolar = 1.0;
     break;
   case 'h' :
@@ -107,5 +124,5 @@
     usage(); return -1;
   }
- }
+ } 
  if(optind>=narg) {usage(); return-1;}
  sscanf(arg[optind],"%lf",&redshift);
@@ -113,4 +130,6 @@
 
  // --- Initialisation de la Cosmologie
+ cout<<"\nh100="<<h100<<" Om0="<<om0<<" Or0="<<or0<<" Or0="
+     <<or0<<" Ol0="<<ol0<<" w0="<<w0<<" flat="<<flat<<endl;
  cout<<"\n--- Cosmology for z = "<<redshift<<endl;
  CosmoCalc univ(flat,true,2.*redshift);
@@ -118,4 +137,5 @@
  univ.SetInteg(perc,dzinc,dzmax,glorder);
  univ.SetDynParam(h100,om0,or0,ol0,w0);
+ univ.Print(0.);
  univ.Print(redshift);
 
@@ -258,4 +278,12 @@
  cout<<"dang lumi = "<<dlum<<" in ["<<dlumlim[0]<<","<<dlumlim[1]<<"] Mpc"<<endl;
 
+ double slobe = lobewidth/2.35482; // sigma du lobe en arcmin
+ double lobecyl = sqrt(8.)*slobe; // diametre du lobe cylindrique equiv en arcmin
+ double lobearea = M_PI*lobecyl*lobecyl/4.;  // en arcmin^2 (hypothese lobe gaussien)
+ cout<<"\nBeam FWHM = "<<lobewidth<<"\' -> sigma = "<<slobe<<"\' -> "
+     <<" Dcyl = "<<lobecyl<<"\' -> area = "<<lobearea<<" arcmin^2"<<endl;
+ double nlobes = rad2min(adtx)*rad2min(adty)/lobearea;
+ cout<<"Number of beams in one transversal pixel = "<<nlobes<<endl;
+
  // --- Power emitted by HI
  cout<<"\n--- Power from HI for M = "<<mhiref<<" Msol at "<<nuhiz<<" GHz"<<endl;
@@ -267,9 +295,11 @@
      <<"      in ["<<hifactor*FluxHI(mhiref,dlumlim[0])
      <<","<<hifactor*FluxHI(mhiref,dlumlim[1])<<"] W/m^2"<<endl;
- double sfhi = fhi / (dnuhiz*1.e+9) / Jansky2Watt_cst;
+ double sfhi = fhi / (dnuhiz*1e9) / Jansky2Watt_cst;
  cout<<"If spread over "<<dnuhiz<<" GHz, flux density = "<<sfhi<<" Jy"<<endl;
 
  // --- Signal analysis
  cout<<"\n--- Signal analysis"<<endl;
+ cout<<"Facteur polar = "<<facpolar<<endl;
+
  PlanckSpectra planck(T_CMB_Par);
  planck.SetApprox(1);  // Rayleigh spectra
@@ -278,50 +308,49 @@
  planck.SetTypSpectra(0); // radiance W/m^2/Sr/Hz
 
- double hnu = h_Planck_Cst*(nuhiz*1.e+9);
- cout<<"Facteur polar = "<<facpolar<<",   h.nu="<<hnu<<" J"<<endl;
-
+ // Signal
  double psig = facpolar * fhi * surfeff;
- double esig = psig * tobs;
- double nsig = esig / hnu;
- double tsig = psig / k_Boltzman_Cst / (dnuhiz*1.e+9);
- double ssig = psig / surfeff / (dnuhiz*1.e+9) / Jansky2Watt_cst;
- cout<<"Signal("<<mhiref<<" Msol): P="<<psig<<" W -> E="<<esig<<" J -> "<<nsig<<" ph"<<endl;
+ double tsig = psig / k_Boltzman_Cst / (dnuhiz*1e9);
+ double ssig = psig / surfeff / (dnuhiz*1e9) / Jansky2Watt_cst;
+ cout<<"Signal("<<mhiref<<" Msol): P="<<psig<<" W"<<endl;
+ cout<<"     flux density = "<<ssig<<" Jy  (for Dnu="<<dnuhiz<<" GHz)"<<endl;
  cout<<"     Antenna temperature: tsig="<<tsig<<" K"<<endl;
- cout<<"     flux density = "<<ssig<<" Jy"<<endl;
-
- double facnoise = facpolar * tobs * surfeff * angsol * (dnuhiz*1.e+9);
-
+
+ // Elargissement doppler de la raie a 21cm: dNu = vrot/C * Nu(21cm) / (1+z)
+ double doplarge = vrot / SpeedOfLight_Cst * nuhiz;
+ cout<<"     Doppler width="<<doplarge*1.e3<<" MHz  for rotation width of "<<vrot<<" km/s"<<endl;
+ if(doplarge>dnuhiz) {
+   cout<<"Warning: doppler width "<<doplarge<<" GHz > "<<dnuhiz<<" GHz redshift bin width"<<endl;
+ }
+
+ // Synchrotron
  double tsynch = Tsynch408;
         if(fabs(indnu)>1.e-50) tsynch *= pow(nuhiz/nuhaslam,indnu);
  planck.SetTemperature(tsynch);
- double psynch = facpolar * planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1.e+9);
- double esynch = psynch * tobs;
- double nsynch = esynch / hnu;
- double ssynch = psynch / surfeff / (dnuhiz*1.e+9) / Jansky2Watt_cst;
+ double psynch = facpolar * planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1e9);
+ double ssynch = psynch / surfeff / (dnuhiz*1e9) / Jansky2Watt_cst;
  cout<<"Synchrotron: T="<<Tsynch408<<" K ("<<nuhaslam<<" GHz), "
      <<tsynch<<" K ("<<nuhiz<<" GHz)"<<endl
-     <<"             P="<<psynch<<" W -> E="<<esynch<<" J -> "<<nsynch<<" ph"<<endl;
- cout<<"     flux density = "<<ssynch<<" Jy"<<endl;
-
+     <<"             P="<<psynch<<" W"<<endl;
+ cout<<"             flux density = "<<ssynch<<" Jy"<<endl;
+
+ // CMB
  double tcmb = T_CMB_Par;
  planck.SetTemperature(tcmb);
- double pcmb = facpolar * planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1.e+9);
- double ecmb = pcmb * tobs;
- double ncmb = ecmb / hnu;
+ double pcmb = facpolar * planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1e9);
  double scmb = pcmb / surfeff / (dnuhiz*1.e+9) / Jansky2Watt_cst;
- cout<<"CMB: T="<<tcmb<<" K -> P="<<pcmb<<" W -> E="<<ecmb<<" J -> "<<ncmb<<" ph"<<endl;
- cout<<"     flux density = "<<scmb<<" Jy"<<endl;
+ cout<<"CMB: T="<<tcmb<<" K -> P="<<pcmb<<" W"<<endl;
+ cout<<"                       flux density = "<<scmb<<" Jy"<<endl;
 
  // --- Noise analysis
  cout<<"\n--- Noise analysis"<<endl;
- double pbruit = k_Boltzman_Cst * Tbruit * (dnuhiz*1.e+9);
- double ebruit = pbruit * tobs;
- cout<<"Bruit: T="<<Tbruit<<" K, P="<<pbruit<<" W -> E="<<ebruit<<" J"<<endl;
-
- double seq = (1./facpolar) * k_Boltzman_Cst * Tbruit / surfeff /Jansky2Watt_cst;
- cout<<"\nSystem equivalent flux density: "<<seq<<" Jy"<<endl;
-
- double slim = seq / sqrt(2.*(dnuhiz*1.e+9)*tobs);
- cout<<"Observation flux density limit: "<<slim<<" Jy"<<endl;
+ double psys = k_Boltzman_Cst * Tsys * (dnuhiz*1.e+9);
+ cout<<"Noise: T="<<Tsys<<" K, P="<<psys<<" W  (for Dnu="<<dnuhiz<<" GHz)"<<endl;
+
+ double slim = 2. * k_Boltzman_Cst * Tsys / surfeff
+               / sqrt(2.*(dnuhiz*1.e+9)*tobs) /Jansky2Watt_cst;
+ cout<<"Observation flux density limit: "<<slim<<" Jy (in 1 lobe)"<<endl;
+
+ slim = slim * sqrt(nlobes);
+ cout<<"Observation flux density limit: "<<slim<<" Jy (in "<<nlobes<<" lobes)"<<endl;
 
  double SsN = ssig/slim;
Index: /trunk/Cosmo/SimLSS/cmvobserv3d.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvobserv3d.cc	(revision 3192)
+++ /trunk/Cosmo/SimLSS/cmvobserv3d.cc	(revision 3193)
@@ -31,5 +31,6 @@
      <<" -s snoise : sigma du bruit en Msol"<<endl
      <<" -2 : compute 2D spectrum"<<endl
-     <<" -F : write cube in FITS format"<<endl
+     <<" -M schmin,schmax,nsch : min,max mass and nb points for schechter HI"<<endl
+     <<" -W : write cube in FITS format"<<endl
      <<" -P : write cube in PPF format"<<endl
      <<" -V : compute variance from real space"<<endl
@@ -65,6 +66,6 @@
 
  // *** Schechter mass function definition
- double h75 = 0.71 / 0.75;
- double nstar = 0.006*pow(h75,3.); 
+ double h75 = h100 / 0.75;
+ double nstar = 0.006*pow(h75,3.);
  double mstar = pow(10.,9.8/(h75*h75));  // MSol
  double alpha = -1.37;
@@ -72,5 +73,4 @@
  double schmin=1e8, schmax=1e12;
  int schnpt = 1000;
- double lschmin=log10(schmin), lschmax=log10(schmax), dlsch=(lschmax-lschmin)/schnpt;
 
  // *** Niveau de bruit
@@ -91,5 +91,5 @@
 
  char c;
- while((c = getopt(narg,arg,"ha0PWV2Gx:y:z:s:Z:")) != -1) {
+ while((c = getopt(narg,arg,"ha0PWV2Gx:y:z:s:Z:M:")) != -1) {
   switch (c) {
   case 'a' :
@@ -119,4 +119,7 @@
   case '2' :
     comp2dspec = true;
+    break;
+  case 'M' :
+    sscanf(optarg,"%lf,%lf,%d",&schmin,&schmax,&schnpt);
     break;
   case 'V' :
@@ -134,4 +137,6 @@
   }
  }
+
+ double lschmin=log10(schmin), lschmax=log10(schmax), dlsch=(lschmax-lschmin)/schnpt;
 
  string tagobs = "cmvobserv3d.ppf";
@@ -157,8 +162,9 @@
  univ.Print();
  double loscomref = univ.Dloscom(zref);
- cout<<"zref = "<<zref<<" -> dloscom = "<<loscomref<<" Mpc"<<endl;
-
- //-----------------------------------------------------------------
-  cout<<endl<<"\n--- Create Spectrum and mass function"<<endl;
+ cout<<"\nzref = "<<zref<<" -> dloscom = "<<loscomref<<" Mpc"<<endl;
+ univ.Print(zref);
+
+ //-----------------------------------------------------------------
+ cout<<endl<<"\n--- Create Spectrum"<<endl;
 
  InitialSpectrum pkini(ns,as);
@@ -168,4 +174,5 @@
 
  GrowthFactor growth(om0,ol0);
+ // GrowthFactor growth(1.,0.); // D(z) = 1/(1+z)
 
  PkSpectrum0 pk0(pkini,tf);
@@ -173,5 +180,9 @@
  PkSpectrumZ pkz(pk0,growth,zref);
  
+ //-----------------------------------------------------------------
+ cout<<endl<<"\n--- Create mass function"<<endl;
+
  Schechter sch(nstar,mstar,alpha);
+ sch.Print();
 
  //-----------------------------------------------------------------
@@ -234,4 +245,6 @@
  double mass_by_mpc3 = IntegrateFuncLog(sch,lschmin,lschmax,0.01,dlsch,10.*dlsch,4);
  cout<<"Galaxy mass density= "<<mass_by_mpc3<<" Msol/Mpc^3 between limits"<<endl;
+ cout<<"Omega_HI at z=0 is "<<mass_by_mpc3/(univ.Rhoc(0.)*GCm3toMsolMpc3_Cst)<<endl
+     <<"         at z="<<zref<<" is "<<mass_by_mpc3/(univ.Rhoc(zref)*GCm3toMsolMpc3_Cst)<<endl;
 
  PrtTim(">>>> End of definition");
@@ -463,10 +476,12 @@
  }
 
- cout<<"\n--- Add noise to HI Flux snoise="<<snoise<<endl;
- fluct3d.AddNoise2Real(snoise);
- nm = fluct3d.MeanSigma2(rm,rs2);
- cout<<"HI mass with noise: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
-     <<rs2<<" -> "<<sqrt(rs2)<<endl;
- PrtTim(">>>> End Add noise");
+ if(snoise>0.) {
+   cout<<"\n--- Add noise to HI Flux snoise="<<snoise<<endl;
+   fluct3d.AddNoise2Real(snoise);
+   nm = fluct3d.MeanSigma2(rm,rs2);
+   cout<<"HI mass with noise: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
+       <<rs2<<" -> "<<sqrt(rs2)<<endl;
+   PrtTim(">>>> End Add noise");
+ }
 
  //-----------------------------------------------------------------
Index: /trunk/Cosmo/SimLSS/cmvtstpk.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvtstpk.cc	(revision 3192)
+++ /trunk/Cosmo/SimLSS/cmvtstpk.cc	(revision 3193)
@@ -22,5 +22,6 @@
      <<" -H h100 -B Ob0 -M Om0 -L Ol0,w0"<<endl
      <<" -k npt,lkmin,lkmax : les valeurs sont en log10"<<endl
-     <<" -s scale : on multiplie pkz par scale"<<endl;
+     <<" -s scale : on multiplie pkz par scale"<<endl
+     <<" -w : write spectra on ASCII file"<<endl;
 }
 
@@ -41,7 +42,8 @@
  double lkmin = -3., lkmax=2.;
  double scale = 1.;
+ bool wrascii = false;
 
  char c;
-  while((c = getopt(narg,arg,"hk:s:H:M:B:L:")) != -1) {
+  while((c = getopt(narg,arg,"hwk:s:H:M:B:L:")) != -1) {
   switch (c) {
   case 'H' :
@@ -64,4 +66,7 @@
   case 's' :
     sscanf(optarg,"%lf",&scale);
+    break;
+  case 'w' :
+    wrascii = true;
     break;
   case 'h' :
@@ -150,6 +155,6 @@
 
  //--------------------------
- cout<<">>>> ASCII"<<endl;
- if(1) {
+ if(wrascii) {
+   cout<<">>>> ASCII"<<endl;
    FILE * fdata = fopen("cmvtstpk.data","w");
    fprintf(fdata,"# z= %g : k(Mpc^-1) pkini, tf pk0 pkz, tfnosc2 pk0nosc2 pkznosc2, ",zval);
Index: /trunk/Cosmo/SimLSS/cmvtuniv.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvtuniv.cc	(revision 3192)
+++ /trunk/Cosmo/SimLSS/cmvtuniv.cc	(revision 3193)
@@ -57,9 +57,9 @@
  univ1.Print();
 
- //tstprint(univ,z1,z2,dz);
+ tstprint(univ,z1,z2,dz);
  //tstprint(univ,univ1,z1,z2,dz);
  //tstspeed(univ,z1,z2,dz);
  //tstntuple(univ,z1,z2,dz);
- tstinterp(univ,z1,z2,dz);
+ //tstinterp(univ,z1,z2,dz);
 
  return 0;
Index: /trunk/Cosmo/SimLSS/constcosmo.h
===================================================================
--- /trunk/Cosmo/SimLSS/constcosmo.h	(revision 3192)
+++ /trunk/Cosmo/SimLSS/constcosmo.h	(revision 3193)
@@ -2,23 +2,25 @@
 #define CONSTCOSMO_SEEN
 
-static const double MpctoMeters_Cst = 3.0856e+22; // Mpc in meters 
-static const double YeartoSec_Cst = 31556925.2; // seconds in a year
 static const double SpeedOfLight_Cst = 299792.458; // speed of light "km/sec"
-static const double LightYear_Cst = 9.46073042e+15; // 1 Light-Year in "m"
 static const double G_Newton_Cst = 6.6742e-11; // G_N in SI units "m^3 / kg /s^2"
 static const double h_Planck_Cst = 6.6260693e-34;   // Planck constant SI units "J s"
 static const double k_Boltzman_Cst = 1.3806503e-23;  // Boltzman constant SI "J / K"
-static const double ProtonMass_Cst = 1.67262171e-24; // Proton mass in "g"
-static const double SolarMass_Cst = 1.98844e+33; // Solar mass in "g"
 //            Sigma = 2*PI^5*K^4/(15.*C^2*H^3) = PI^2*K^4/(60.*Hbar^3*C2)
 static const double StefBoltz_Cst = 5.670400e-8; // constante de Stefan-Boltzman in "W/m^2/K^4"
 static const double Zeta_3_Cst = 1.2020569031595942854; // Zeta(3)
 
-static const double Jansky2Watt_cst = 1.e-26;  // conversion Jansky to "Watt/m^2/Hz == J/m^2"
-static const double Deg2Rad_Cst = (M_PI/180.); // conversion deg to rad
+static const double YeartoSec_Cst = 31556925.2; // 1 year into seconds
+static const double MpctoMeters_Cst = 3.0856e+22; // 1 Mpc into meters
+static const double LightYear_Cst = 9.46073042e+15; // 1 Light-Year into "m"
+static const double Jansky2Watt_cst = 1.e-26;  // 1 Jansky into "Watt/m^2/Hz == J/m^2"
+static const double SolarMass_Cst = 1.98844e+33; // 1 Solar mass into "g"
+//                  GCm3toMsolMpc3_Cst = (MpctoMeters_Cst*100)^3 / SolarMass_Cst
+static const double GCm3toMsolMpc3_Cst = 1.477428208143872e+40; // 1 g/cm^3 into 1 Msol/Mpc^3
+static const double ProtonMass_Cst = 1.67262171e-24; // 1 "proton mass" in "g"
+static const double Deg2Rad_Cst = (M_PI/180.); // 1 deg into rad
 
 static const double T_CMB_Par = 2.725;  // temperature du CMB en "K"
 static const double T_NU_Par = 1.9;  // temperature des neutrinos en "K"
-static const double Fr_HyperFin_Par = 1.4204;  // frequence de la raie hyperfine HI en "GHz"
+static const double Fr_HyperFin_Par = 1.420405751786;  // frequence de la raie hyperfine HI en "GHz"
 
 #endif
Index: /trunk/Cosmo/SimLSS/cosmocalc.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cosmocalc.cc	(revision 3192)
+++ /trunk/Cosmo/SimLSS/cosmocalc.cc	(revision 3193)
@@ -206,6 +206,6 @@
 
   printf("CosmoCalc::Print(spl=%d,zmax=%g,flat=%u)  for z=%g\n",int(_usespline),ZMax(),Flat(),z);
-  printf("h100=%g H0=%g Dhub=%g  H(z)=%g  Rhoc=%g g/cm^3\n"
-        ,h100(),H0(),Dhubble(),H(z),Rhoc(z));
+  printf("h100=%g H0=%g Dhub=%g  H(z)=%g  Rhoc=%g g/cm^3 =%g Msol/Mpc^3\n"
+        ,h100(),H0(),Dhubble(),H(z),Rhoc(z),Rhoc(z)*GCm3toMsolMpc3_Cst);
   printf("Olambda=%g  W0=%g\n",Olambda(z),_W0);
   printf("Omatter=%g  Obaryon=%g\n",Omatter(z),Obaryon(z));
Index: /trunk/Cosmo/SimLSS/pkspectrum.cc
===================================================================
--- /trunk/Cosmo/SimLSS/pkspectrum.cc	(revision 3192)
+++ /trunk/Cosmo/SimLSS/pkspectrum.cc	(revision 3193)
@@ -244,4 +244,5 @@
 
 // From Eisenstein & Hu ApJ 496:605-614 1998 April 1
+// Pour avoir D(z) = 1/(1+z) faire: OmegaMatter0=1 OmegaLambda0=0
 GrowthFactor::GrowthFactor(double OmegaMatter0,double OmegaLambda0)
   : O0_(OmegaMatter0) , Ol_(OmegaLambda0) , Ok_(1.-OmegaMatter0-OmegaLambda0)
Index: /trunk/Cosmo/SimLSS/schechter.cc
===================================================================
--- /trunk/Cosmo/SimLSS/schechter.cc	(revision 3192)
+++ /trunk/Cosmo/SimLSS/schechter.cc	(revision 3193)
@@ -54,4 +54,14 @@
 }
 
+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;
+}
+
 ///////////////////////////////////////////////////////////
 //******************* Le Flux a 21cm ********************//
@@ -64,11 +74,16 @@
 // Return:
 //    le flux total emis a 21 cm en W/m^2
-// Ref: Binney & Merrifield, Galactic Astronomy p474 (ed:1998)
-//      J.Peterson & K.Bandura, astroph-0606104  (eq 7)
-//      F.Abdalla & Rawlings, astroph-0411342 (eq 7 may pb de d'unites?)
+// 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  m / (2.35e+5 * d*d )
-          * Jansky2Watt_cst
-          * (Fr_HyperFin_Par*1e+9)/SpeedOfLight_Cst;
+  return  2.0e-28 * m / (d*d);
 }
-
Index: /trunk/Cosmo/SimLSS/schechter.h
===================================================================
--- /trunk/Cosmo/SimLSS/schechter.h	(revision 3192)
+++ /trunk/Cosmo/SimLSS/schechter.h	(revision 3193)
@@ -16,4 +16,5 @@
   void SetOutValue(unsigned short outvalue=0);
   virtual double operator() (double m);
+  virtual void Print(void);
 protected:
   double nstar_,mstar_,alpha_;
