Changeset 3287 in Sophya for trunk/Cosmo/SimLSS/cmvdefsurv.cc


Ignore:
Timestamp:
Jul 31, 2007, 3:54:47 PM (18 years ago)
Author:
cmv
Message:

introduction utnies pour arguments cmv 31/07/2007

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Cosmo/SimLSS/cmvdefsurv.cc

    r3271 r3287  
    2727void usage(void);
    2828void usage(void) {
    29   cout<<"cmvdefsurv [-r] -x adtx,atxlarg [-y adty,atylarg] -z dred,redlarg redshift"<<endl
    30       <<" -x adtx,atxlarg : resolution en Theta_x (arcmin), largeur (degre)"<<endl
    31       <<" -y adty,atylarg : idem selon y, si <=0 meme que x"<<endl
    32       <<" -z dred,redlarg : resolution en redshift, largeur en redshift"<<endl
    33       <<" -P : on donne -x -y -z en Mpc au lieu d\'angles et de redshift"<<endl
     29  cout<<"cmvdefsurv [-r] -x adtx,atxlarg[,unit_x] -y adty,atylarg[,unit_y] -z dred,redlarg[,unit_z] redshift"<<endl
     30      <<"----------------"<<endl
     31      <<" -x adtx,atxlarg : resolution et largeur dans le plan transverse selon X"<<endl
     32      <<" -y adty,atylarg : idem selon Y, si <=0 meme que X"<<endl
     33      <<" -z dred,redlarg : resolution et largeursur la ligne de visee"<<endl
     34      <<"-- Unites pour X-Y:"<<endl
     35      <<"   \'A\' : en angles (pour X-Y)     : resolution=ArcMin, largeur=Degre (defaut)"<<endl
     36      <<"   \'Z\' : en redshift (pour Z)     : resolution et largeur en redshift (defaut)"<<endl
     37      <<"   \'F\' : en frequence (pour Z)    : resolution et largeur MHz"<<endl
     38      <<"   \'M\' : en distance (pour X-Y-Z) : resolution et largeur Mpc"<<endl
     39      <<"----------------"<<endl
    3440      <<" -L lobewidth : taille du lobe d\'observation (FWHM) en arcmin (def= 1\')"<<endl
    3541      <<"                Si <=0 alors angle solide du lobe = celui du pixel"<<endl
    3642      <<" -O surf,tobs : surface effective (m^2) et temps d\'observation (s)"<<endl
    3743      <<" -N Tsys : temperature du system (K)"<<endl
    38       <<" -S Tsynch,indnu : temperature (K) synch a 408 Mhz, index d\'evolution"<<endl
    39       <<"                   (indnu==0 no evolution with freq.)"<<endl
     44      <<" -2 : two polarisations measured"<<endl
    4045      <<" -M  : masse de HI de reference (MSol), si <=0 mean schechter in pixel"<<endl
    4146      <<" -F  : HI flux factor to be applied for our redshift"<<endl
    4247      <<" -V Vrot : largeur en vitesse (km/s) pour l\'elargissement doppler (def=300km/s)"<<endl
     48      <<"----------------"<<endl
     49      <<" -S Tsynch,indnu : temperature (K) synch a 408 Mhz, index d\'evolution"<<endl
     50      <<"                   (indnu==0 no evolution with freq.)"<<endl
     51      <<"----------------"<<endl
    4352      <<" -U h100,om0,ol0,w0,or0,flat : cosmology"<<endl
    44       <<" -2 : two polarisations measured"<<endl
     53      <<"----------------"<<endl
    4554      <<" -A <log10(S_agn)> : moyenne du flux AGN en Jy dans le pixel"<<endl
    4655      <<" redshift : redshift moyen du survey"<<endl
     
    6271
    6372 // --- Arguments
    64  bool inmpc = false;
    65  double adtx=1., atxlarg=90., adty=-1., atylarg=-1.;
    66  double dx=1.,txlarg=1000., dy=-1.,tylarg=1000., dz=1.,tzlarg=100.;
    67  int nx,ny,nz;
    68  double redshift = 1., dred=0.01, redlarg=0.3;
     73 double adtx=0., atxlarg=0., dx=0.,txlarg=0.;
     74   int nx=0; char unit_x = 'A';
     75 double adty=-1., atylarg=-1., dy=0.,tylarg=0.;
     76   int ny=0; char unit_y = 'A';
     77 double dred=0., redlarg=0., dz=0.,tzlarg=0.;
     78   int nz=0; char unit_z = 'Z';
     79 double redshift = 0.;
    6980 double tobs = 6000., surfeff = 400000.;
    7081 double lobewidth = 1.;  // taille du lobe d'observation en arcmin
     
    8091 // --- Decodage arguments
    8192 char c;
    82   while((c = getopt(narg,arg,"hP2x:y:z:N:S:O:M:F:V:U:L:A:")) != -1) {
     93  while((c = getopt(narg,arg,"h2x:y:z:N:S:O:M:F:V:U:L:A:")) != -1) {
    8394  switch (c) {
    84   case 'P' :
    85     inmpc = true;
    86     break;
    8795  case 'x' :
    88     sscanf(optarg,"%lf,%lf",&adtx,&atxlarg);
     96    sscanf(optarg,"%lf,%lf,%c",&adtx,&atxlarg,&unit_x);
    8997    break;
    9098  case 'y' :
    91     sscanf(optarg,"%lf,%lf",&adty,&atylarg);
     99    sscanf(optarg,"%lf,%lf,%c",&adty,&atylarg,&unit_y);
    92100    break;
    93101  case 'z' :
    94     sscanf(optarg,"%lf,%lf",&dred,&redlarg);
     102    sscanf(optarg,"%lf,%lf,%c",&dred,&redlarg,&unit_z);
    95103    break;
    96104  case 'O' :
     
    134142
    135143 // --- Initialisation de la Cosmologie
    136  cout<<"\nh100="<<h100<<" Om0="<<om0<<" Or0="<<or0<<" Or0="
     144 cout<<"\n>>>>\n>>>> Cosmologie generale\n>>>>"<<endl;
     145 cout<<"h100="<<h100<<" Om0="<<om0<<" Or0="<<or0<<" Or0="
    137146     <<or0<<" Ol0="<<ol0<<" w0="<<w0<<" flat="<<flat<<endl;
    138  cout<<"\n--- Cosmology for z = "<<redshift<<endl;
     147 cout<<"--- Cosmology for z = "<<redshift<<endl;
    139148 CosmoCalc univ(flat,true,2.*redshift);
    140149 double perc=0.01,dzinc=redshift/100.,dzmax=dzinc*10.; unsigned short glorder=4;
     
    163172
    164173 // --- Mise en forme dans les unites appropriees
    165  if(adty<=0.) adty=adtx;
    166  if(atylarg<=0.) atylarg=atxlarg;
    167  if(inmpc) {
    168    dx = adtx; txlarg = atxlarg; nx = int(txlarg/dx+0.5);
    169    dy = adty; txlarg = atxlarg; ny = int(tylarg/dy+0.5);
    170    dz = dred; tzlarg = redlarg; nz = int(tzlarg/dz+0.5);
     174 cout<<"\n>>>>\n>>>> Geometrie\n>>>>"<<endl;
     175 if(adty<=0. || atylarg<=0.) {adty=adtx; atylarg=atxlarg; unit_y=unit_x;}
     176 cout<<"X values: resolution="<<adtx<<" largeur="<<atxlarg<<" unite="<<unit_x<<endl;
     177 if(unit_x == 'A') {
     178   nx = int(atxlarg*60./adtx+0.5);
     179   adtx = min2rad(adtx); atxlarg = deg2rad(atxlarg);
     180   dx = adtx*dtrcom; txlarg = dx*nx;
     181 } else if(unit_x == 'M') {
     182   nx = int(atxlarg/adtx+0.5);
     183   dx = adtx; txlarg = atxlarg;
    171184   adtx = dx/dtrcom; atxlarg = adtx*nx;
     185 } else {
     186   cout<<"Unknown unit_x = "<<unit_x<<endl;
     187 }
     188 cout<<"Y values: resolution="<<adty<<" largeur="<<atylarg<<" unite="<<unit_y<<endl;
     189 if(unit_y == 'A') {
     190   ny = int(atylarg*60./adty+0.5);
     191   adty = min2rad(adty); atylarg = deg2rad(atylarg);
     192   dy = adty*dtrcom; tylarg = dy*ny;
     193 } else if(unit_y == 'M') {
     194   ny = int(atylarg/adty+0.5);
     195   dy = adty; tylarg = atylarg;
    172196   adty = dy/dtrcom; atylarg = adty*ny;
     197 } else {
     198   cout<<"Unknown unit_y = "<<unit_y<<endl;
     199 }
     200 cout<<"Z values: resolution="<<dred<<" largeur="<<redlarg<<" unite="<<unit_z<<endl;
     201 if(unit_z == 'Z') {
     202   nz = int(redlarg/dred+0.5);
     203   dz = dred*dlosdz; tzlarg = dz*nz;
     204 } else if(unit_z == 'M') {
     205   nz = int(redlarg/dred+0.5);
     206   dz = dred; tzlarg = redlarg;
    173207   dred = dz/dlosdz; redlarg = dred*nz;
     208 } else if(unit_z == 'F') {
     209   nz = int(redlarg/dred+0.5);
     210   dred = dred/(Fr_HyperFin_Par*1.e3)*pow(1.+redshift,2.); redlarg = dred*nz;
     211   dz = dred*dlosdz; tzlarg = dz*nz;
    174212 } else {
    175    adtx = min2rad(adtx); atxlarg = deg2rad(atxlarg); nx = int(atxlarg/adtx+0.5);
    176    adty = min2rad(adty); atylarg = deg2rad(atylarg); ny = int(atylarg/adty+0.5);
    177    nz = int(redlarg/dred+0.5);
    178    dx = adtx*dtrcom; txlarg = dx*nx;
    179    dy = adty*dtrcom; tylarg = dy*ny;
    180    dz = dred*dlosdz; tzlarg = dz*nz;
     213   cout<<"Unknown unit_z = "<<unit_z<<endl;
    181214 }
     215
    182216 double Npix = (double)nx*(double)ny*(double)nz;
    183 
    184217 double redlim[2] = {redshift-redlarg/2.,redshift+redlarg/2.};
    185218 if(redlim[0]<=0.)
     
    189222 double dlumlim[2] = {univ.Dlum(redlim[0]),univ.Dlum(redlim[1])};
    190223
    191  cout<<"\n---- Type de donnees: inmpc = "<<inmpc<<endl;
    192224 cout<<"---- Line of Sight: Redshift = "<<redshift<<endl
    193225     <<"dred = "<<dred<<"  redlarg = "<<redlarg<<endl
     
    202234
    203235 // --- Cosmolographie Transverse
    204  cout<<"\n--- Transverse"<<endl;
     236 cout<<"\n>>>>\n>>>> Cosmologie & Geometrie transverse\n>>>>"<<endl;
    205237 cout<<"dang comoving = "<<dtrcom<<" Mpc (com) var_in_z ["
    206238                         <<dtrlim[0]<<","<<dtrlim[1]<<"]"<<endl;
     
    217249
    218250 // --- Cosmolographie Line of sight
    219  cout<<"\n--- Line of Sight"<<endl;
     251 cout<<"\n>>>>\n>>>> Cosmologie & Geometrie ligne de visee\n>>>>"<<endl;
    220252 cout<<"los comoving distance = "<<dloscom<<" Mpc (com) in ["
    221253     <<loslim[0]<<","<<loslim[1]<<"]"<<endl
     
    227259
    228260 // --- Solid Angle & Volume
    229  cout<<"\n--- Solid angle"<<endl;
     261 cout<<"\n>>>>\n>>>> Angles solides et Volumes\n>>>>"<<endl;
     262 cout<<"--- Solid angle"<<endl;
    230263 double angsol = AngSol(adtx/2.,adty/2.,M_PI/2.);
    231264 cout<<"Elementary solid angle = "<<angsol<<" sr = "<<angsol/(4.*M_PI)<<" *4Pi sr"<<endl;
     
    241274
    242275 // --- Fourier space: k = omega = 2*Pi*Nu
    243  cout<<"\n--- Fourier space"<<endl;
     276 cout<<"\n>>>>\n>>>> Geometrie dans l'espace de Fourier\n>>>>"<<endl;
    244277 cout<<"Array size: nx = "<<nx<<",  ny = "<<ny<<",  nz = "<<nz<<endl;
    245278 double dk_x = 2.*M_PI/(nx*dx), knyq_x = M_PI/dx;
     
    254287
    255288 // --- Masse de HI
    256  cout<<"\n--- Mass HI"<<endl;
     289 cout<<"\n>>>>\n>>>> Mass HI\n>>>>"<<endl;
    257290 Schechter sch(nstar,mstar,alpha);
    258291 sch.SetOutValue(1);
     
    271304
    272305 // --- Survey values
     306 cout<<"\n>>>>\n>>>> Observations\n>>>>"<<endl;
    273307 double unplusz = 1.+redshift;
    274308 double nuhiz = Fr_HyperFin_Par / unplusz;  // GHz
     
    278312 double dnuhiz = Fr_HyperFin_Par *dred/(unplusz*unplusz)
    279313               / (1.- pow(dred/.2/unplusz,2.));
    280  cout<<"\n--- Observation:"<<endl
    281      <<"      surf_eff="<<surfeff<<" m^2, tobs="<<tobs<<" s"<<endl
     314 cout<<"      surf_eff="<<surfeff<<" m^2, tobs="<<tobs<<" s"<<endl
    282315     <<"      nu="<<nuhiz<<" GHz, dnu="<<dnuhiz*1.e3<<" Mhz"<<endl;
    283316 cout<<"dang lumi = "<<dlum<<" in ["<<dlumlim[0]<<","<<dlumlim[1]<<"] Mpc"<<endl;
     
    330363   cout<<"Warning: doppler width "<<doplarge<<" GHz > "<<dnuhiz<<" GHz redshift bin width"<<endl;
    331364
    332  // Synchrotron
     365 // Synchrotron  (T en -2.7  -> Flux en -0.7  dans l'approximation Rayleigh)
    333366 double tsynch = Tsynch408;
    334367        if(fabs(indnu)>1.e-50) tsynch *= pow(nuhiz/nuhaslam,indnu);
Note: See TracChangeset for help on using the changeset viewer.