Index: /trunk/Cosmo/SimLSS/cmvdefsurv.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvdefsurv.cc	(revision 3286)
+++ /trunk/Cosmo/SimLSS/cmvdefsurv.cc	(revision 3287)
@@ -27,20 +27,29 @@
 void usage(void);
 void usage(void) {
-  cout<<"cmvdefsurv [-r] -x adtx,atxlarg [-y adty,atylarg] -z dred,redlarg redshift"<<endl
-      <<" -x adtx,atxlarg : resolution en Theta_x (arcmin), largeur (degre)"<<endl
-      <<" -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 au lieu d\'angles et de redshift"<<endl
+  cout<<"cmvdefsurv [-r] -x adtx,atxlarg[,unit_x] -y adty,atylarg[,unit_y] -z dred,redlarg[,unit_z] redshift"<<endl
+      <<"----------------"<<endl
+      <<" -x adtx,atxlarg : resolution et largeur dans le plan transverse selon X"<<endl
+      <<" -y adty,atylarg : idem selon Y, si <=0 meme que X"<<endl
+      <<" -z dred,redlarg : resolution et largeursur la ligne de visee"<<endl
+      <<"-- Unites pour X-Y:"<<endl
+      <<"   \'A\' : en angles (pour X-Y)     : resolution=ArcMin, largeur=Degre (defaut)"<<endl
+      <<"   \'Z\' : en redshift (pour Z)     : resolution et largeur en redshift (defaut)"<<endl
+      <<"   \'F\' : en frequence (pour Z)    : resolution et largeur MHz"<<endl
+      <<"   \'M\' : en distance (pour X-Y-Z) : resolution et largeur Mpc"<<endl
+      <<"----------------"<<endl
       <<" -L lobewidth : taille du lobe d\'observation (FWHM) en arcmin (def= 1\')"<<endl
       <<"                Si <=0 alors angle solide du lobe = celui du pixel"<<endl
       <<" -O surf,tobs : surface effective (m^2) et temps d\'observation (s)"<<endl
       <<" -N 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
+      <<" -2 : two polarisations measured"<<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
       <<" -V Vrot : largeur en vitesse (km/s) pour l\'elargissement doppler (def=300km/s)"<<endl
+      <<"----------------"<<endl
+      <<" -S Tsynch,indnu : temperature (K) synch a 408 Mhz, index d\'evolution"<<endl
+      <<"                   (indnu==0 no evolution with freq.)"<<endl
+      <<"----------------"<<endl
       <<" -U h100,om0,ol0,w0,or0,flat : cosmology"<<endl
-      <<" -2 : two polarisations measured"<<endl
+      <<"----------------"<<endl
       <<" -A <log10(S_agn)> : moyenne du flux AGN en Jy dans le pixel"<<endl
       <<" redshift : redshift moyen du survey"<<endl
@@ -62,9 +71,11 @@
 
  // --- Arguments
- bool inmpc = false;
- double adtx=1., atxlarg=90., adty=-1., atylarg=-1.;
- double dx=1.,txlarg=1000., dy=-1.,tylarg=1000., dz=1.,tzlarg=100.;
- int nx,ny,nz;
- double redshift = 1., dred=0.01, redlarg=0.3;
+ double adtx=0., atxlarg=0., dx=0.,txlarg=0.;
+   int nx=0; char unit_x = 'A';
+ double adty=-1., atylarg=-1., dy=0.,tylarg=0.;
+   int ny=0; char unit_y = 'A';
+ double dred=0., redlarg=0., dz=0.,tzlarg=0.;
+   int nz=0; char unit_z = 'Z';
+ double redshift = 0.;
  double tobs = 6000., surfeff = 400000.;
  double lobewidth = 1.;  // taille du lobe d'observation en arcmin
@@ -80,17 +91,14 @@
  // --- Decodage arguments
  char c;
-  while((c = getopt(narg,arg,"hP2x:y:z:N:S:O:M:F:V:U:L:A:")) != -1) {
+  while((c = getopt(narg,arg,"h2x:y:z:N:S:O:M:F:V:U:L:A:")) != -1) {
   switch (c) {
-  case 'P' :
-    inmpc = true;
-    break;
   case 'x' :
-    sscanf(optarg,"%lf,%lf",&adtx,&atxlarg);
+    sscanf(optarg,"%lf,%lf,%c",&adtx,&atxlarg,&unit_x);
     break;
   case 'y' :
-    sscanf(optarg,"%lf,%lf",&adty,&atylarg);
+    sscanf(optarg,"%lf,%lf,%c",&adty,&atylarg,&unit_y);
     break;
   case 'z' :
-    sscanf(optarg,"%lf,%lf",&dred,&redlarg);
+    sscanf(optarg,"%lf,%lf,%c",&dred,&redlarg,&unit_z);
     break;
   case 'O' :
@@ -134,7 +142,8 @@
 
  // --- Initialisation de la Cosmologie
- cout<<"\nh100="<<h100<<" Om0="<<om0<<" Or0="<<or0<<" Or0="
+ cout<<"\n>>>>\n>>>> Cosmologie generale\n>>>>"<<endl;
+ cout<<"h100="<<h100<<" Om0="<<om0<<" Or0="<<or0<<" Or0="
      <<or0<<" Ol0="<<ol0<<" w0="<<w0<<" flat="<<flat<<endl;
- cout<<"\n--- Cosmology for z = "<<redshift<<endl;
+ cout<<"--- Cosmology for z = "<<redshift<<endl;
  CosmoCalc univ(flat,true,2.*redshift);
  double perc=0.01,dzinc=redshift/100.,dzmax=dzinc*10.; unsigned short glorder=4;
@@ -163,23 +172,47 @@
 
  // --- Mise en forme dans les unites appropriees
- if(adty<=0.) adty=adtx;
- if(atylarg<=0.) atylarg=atxlarg;
- if(inmpc) {
-   dx = adtx; txlarg = atxlarg; nx = int(txlarg/dx+0.5);
-   dy = adty; txlarg = atxlarg; ny = int(tylarg/dy+0.5);
-   dz = dred; tzlarg = redlarg; nz = int(tzlarg/dz+0.5);
+ cout<<"\n>>>>\n>>>> Geometrie\n>>>>"<<endl;
+ if(adty<=0. || atylarg<=0.) {adty=adtx; atylarg=atxlarg; unit_y=unit_x;}
+ cout<<"X values: resolution="<<adtx<<" largeur="<<atxlarg<<" unite="<<unit_x<<endl;
+ if(unit_x == 'A') {
+   nx = int(atxlarg*60./adtx+0.5);
+   adtx = min2rad(adtx); atxlarg = deg2rad(atxlarg);
+   dx = adtx*dtrcom; txlarg = dx*nx;
+ } else if(unit_x == 'M') {
+   nx = int(atxlarg/adtx+0.5);
+   dx = adtx; txlarg = atxlarg;
    adtx = dx/dtrcom; atxlarg = adtx*nx;
+ } else {
+   cout<<"Unknown unit_x = "<<unit_x<<endl;
+ }
+ cout<<"Y values: resolution="<<adty<<" largeur="<<atylarg<<" unite="<<unit_y<<endl;
+ if(unit_y == 'A') {
+   ny = int(atylarg*60./adty+0.5);
+   adty = min2rad(adty); atylarg = deg2rad(atylarg);
+   dy = adty*dtrcom; tylarg = dy*ny;
+ } else if(unit_y == 'M') {
+   ny = int(atylarg/adty+0.5);
+   dy = adty; tylarg = atylarg;
    adty = dy/dtrcom; atylarg = adty*ny;
+ } else {
+   cout<<"Unknown unit_y = "<<unit_y<<endl;
+ }
+ cout<<"Z values: resolution="<<dred<<" largeur="<<redlarg<<" unite="<<unit_z<<endl;
+ if(unit_z == 'Z') {
+   nz = int(redlarg/dred+0.5);
+   dz = dred*dlosdz; tzlarg = dz*nz;
+ } else if(unit_z == 'M') {
+   nz = int(redlarg/dred+0.5);
+   dz = dred; tzlarg = redlarg;
    dred = dz/dlosdz; redlarg = dred*nz;
+ } else if(unit_z == 'F') {
+   nz = int(redlarg/dred+0.5);
+   dred = dred/(Fr_HyperFin_Par*1.e3)*pow(1.+redshift,2.); redlarg = dred*nz;
+   dz = dred*dlosdz; tzlarg = dz*nz;
  } else {
-   adtx = min2rad(adtx); atxlarg = deg2rad(atxlarg); nx = int(atxlarg/adtx+0.5);
-   adty = min2rad(adty); atylarg = deg2rad(atylarg); ny = int(atylarg/adty+0.5);
-   nz = int(redlarg/dred+0.5);
-   dx = adtx*dtrcom; txlarg = dx*nx;
-   dy = adty*dtrcom; tylarg = dy*ny;
-   dz = dred*dlosdz; tzlarg = dz*nz;
+   cout<<"Unknown unit_z = "<<unit_z<<endl;
  }
+
  double Npix = (double)nx*(double)ny*(double)nz;
-
  double redlim[2] = {redshift-redlarg/2.,redshift+redlarg/2.};
  if(redlim[0]<=0.)
@@ -189,5 +222,4 @@
  double dlumlim[2] = {univ.Dlum(redlim[0]),univ.Dlum(redlim[1])};
 
- cout<<"\n---- Type de donnees: inmpc = "<<inmpc<<endl;
  cout<<"---- Line of Sight: Redshift = "<<redshift<<endl
      <<"dred = "<<dred<<"  redlarg = "<<redlarg<<endl
@@ -202,5 +234,5 @@
 
  // --- Cosmolographie Transverse
- cout<<"\n--- Transverse"<<endl;
+ cout<<"\n>>>>\n>>>> Cosmologie & Geometrie transverse\n>>>>"<<endl;
  cout<<"dang comoving = "<<dtrcom<<" Mpc (com) var_in_z ["
                          <<dtrlim[0]<<","<<dtrlim[1]<<"]"<<endl;
@@ -217,5 +249,5 @@
 
  // --- Cosmolographie Line of sight
- cout<<"\n--- Line of Sight"<<endl;
+ cout<<"\n>>>>\n>>>> Cosmologie & Geometrie ligne de visee\n>>>>"<<endl;
  cout<<"los comoving distance = "<<dloscom<<" Mpc (com) in ["
      <<loslim[0]<<","<<loslim[1]<<"]"<<endl
@@ -227,5 +259,6 @@
 
  // --- Solid Angle & Volume
- cout<<"\n--- Solid angle"<<endl;
+ cout<<"\n>>>>\n>>>> Angles solides et Volumes\n>>>>"<<endl;
+ cout<<"--- Solid angle"<<endl;
  double angsol = AngSol(adtx/2.,adty/2.,M_PI/2.);
  cout<<"Elementary solid angle = "<<angsol<<" sr = "<<angsol/(4.*M_PI)<<" *4Pi sr"<<endl;
@@ -241,5 +274,5 @@
 
  // --- Fourier space: k = omega = 2*Pi*Nu
- cout<<"\n--- Fourier space"<<endl;
+ cout<<"\n>>>>\n>>>> Geometrie dans l'espace de Fourier\n>>>>"<<endl;
  cout<<"Array size: nx = "<<nx<<",  ny = "<<ny<<",  nz = "<<nz<<endl;
  double dk_x = 2.*M_PI/(nx*dx), knyq_x = M_PI/dx;
@@ -254,5 +287,5 @@
 
  // --- Masse de HI
- cout<<"\n--- Mass HI"<<endl;
+ cout<<"\n>>>>\n>>>> Mass HI\n>>>>"<<endl;
  Schechter sch(nstar,mstar,alpha);
  sch.SetOutValue(1);
@@ -271,4 +304,5 @@
 
  // --- Survey values
+ cout<<"\n>>>>\n>>>> Observations\n>>>>"<<endl;
  double unplusz = 1.+redshift;
  double nuhiz = Fr_HyperFin_Par / unplusz;  // GHz
@@ -278,6 +312,5 @@
  double dnuhiz = Fr_HyperFin_Par *dred/(unplusz*unplusz)
                / (1.- pow(dred/.2/unplusz,2.));
- cout<<"\n--- Observation:"<<endl
-     <<"      surf_eff="<<surfeff<<" m^2, tobs="<<tobs<<" s"<<endl
+ cout<<"      surf_eff="<<surfeff<<" m^2, tobs="<<tobs<<" s"<<endl
      <<"      nu="<<nuhiz<<" GHz, dnu="<<dnuhiz*1.e3<<" Mhz"<<endl;
  cout<<"dang lumi = "<<dlum<<" in ["<<dlumlim[0]<<","<<dlumlim[1]<<"] Mpc"<<endl;
@@ -330,5 +363,5 @@
    cout<<"Warning: doppler width "<<doplarge<<" GHz > "<<dnuhiz<<" GHz redshift bin width"<<endl;
 
- // Synchrotron
+ // Synchrotron  (T en -2.7  -> Flux en -0.7  dans l'approximation Rayleigh)
  double tsynch = Tsynch408;
         if(fabs(indnu)>1.e-50) tsynch *= pow(nuhiz/nuhaslam,indnu);
