Changeset 3287 in Sophya for trunk/Cosmo/SimLSS/cmvdefsurv.cc
- Timestamp:
- Jul 31, 2007, 3:54:47 PM (18 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Cosmo/SimLSS/cmvdefsurv.cc
r3271 r3287 27 27 void usage(void); 28 28 void 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 34 40 <<" -L lobewidth : taille du lobe d\'observation (FWHM) en arcmin (def= 1\')"<<endl 35 41 <<" Si <=0 alors angle solide du lobe = celui du pixel"<<endl 36 42 <<" -O surf,tobs : surface effective (m^2) et temps d\'observation (s)"<<endl 37 43 <<" -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 40 45 <<" -M : masse de HI de reference (MSol), si <=0 mean schechter in pixel"<<endl 41 46 <<" -F : HI flux factor to be applied for our redshift"<<endl 42 47 <<" -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 43 52 <<" -U h100,om0,ol0,w0,or0,flat : cosmology"<<endl 44 <<" -2 : two polarisations measured"<<endl53 <<"----------------"<<endl 45 54 <<" -A <log10(S_agn)> : moyenne du flux AGN en Jy dans le pixel"<<endl 46 55 <<" redshift : redshift moyen du survey"<<endl … … 62 71 63 72 // --- 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.; 69 80 double tobs = 6000., surfeff = 400000.; 70 81 double lobewidth = 1.; // taille du lobe d'observation en arcmin … … 80 91 // --- Decodage arguments 81 92 char c; 82 while((c = getopt(narg,arg,"h P2x: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) { 83 94 switch (c) { 84 case 'P' :85 inmpc = true;86 break;87 95 case 'x' : 88 sscanf(optarg,"%lf,%lf ",&adtx,&atxlarg);96 sscanf(optarg,"%lf,%lf,%c",&adtx,&atxlarg,&unit_x); 89 97 break; 90 98 case 'y' : 91 sscanf(optarg,"%lf,%lf ",&adty,&atylarg);99 sscanf(optarg,"%lf,%lf,%c",&adty,&atylarg,&unit_y); 92 100 break; 93 101 case 'z' : 94 sscanf(optarg,"%lf,%lf ",&dred,&redlarg);102 sscanf(optarg,"%lf,%lf,%c",&dred,&redlarg,&unit_z); 95 103 break; 96 104 case 'O' : … … 134 142 135 143 // --- 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=" 137 146 <<or0<<" Ol0="<<ol0<<" w0="<<w0<<" flat="<<flat<<endl; 138 cout<<" \n--- Cosmology for z = "<<redshift<<endl;147 cout<<"--- Cosmology for z = "<<redshift<<endl; 139 148 CosmoCalc univ(flat,true,2.*redshift); 140 149 double perc=0.01,dzinc=redshift/100.,dzmax=dzinc*10.; unsigned short glorder=4; … … 163 172 164 173 // --- 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; 171 184 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; 172 196 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; 173 207 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; 174 212 } 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; 181 214 } 215 182 216 double Npix = (double)nx*(double)ny*(double)nz; 183 184 217 double redlim[2] = {redshift-redlarg/2.,redshift+redlarg/2.}; 185 218 if(redlim[0]<=0.) … … 189 222 double dlumlim[2] = {univ.Dlum(redlim[0]),univ.Dlum(redlim[1])}; 190 223 191 cout<<"\n---- Type de donnees: inmpc = "<<inmpc<<endl;192 224 cout<<"---- Line of Sight: Redshift = "<<redshift<<endl 193 225 <<"dred = "<<dred<<" redlarg = "<<redlarg<<endl … … 202 234 203 235 // --- Cosmolographie Transverse 204 cout<<"\n --- Transverse"<<endl;236 cout<<"\n>>>>\n>>>> Cosmologie & Geometrie transverse\n>>>>"<<endl; 205 237 cout<<"dang comoving = "<<dtrcom<<" Mpc (com) var_in_z [" 206 238 <<dtrlim[0]<<","<<dtrlim[1]<<"]"<<endl; … … 217 249 218 250 // --- Cosmolographie Line of sight 219 cout<<"\n --- Line of Sight"<<endl;251 cout<<"\n>>>>\n>>>> Cosmologie & Geometrie ligne de visee\n>>>>"<<endl; 220 252 cout<<"los comoving distance = "<<dloscom<<" Mpc (com) in [" 221 253 <<loslim[0]<<","<<loslim[1]<<"]"<<endl … … 227 259 228 260 // --- Solid Angle & Volume 229 cout<<"\n--- Solid angle"<<endl; 261 cout<<"\n>>>>\n>>>> Angles solides et Volumes\n>>>>"<<endl; 262 cout<<"--- Solid angle"<<endl; 230 263 double angsol = AngSol(adtx/2.,adty/2.,M_PI/2.); 231 264 cout<<"Elementary solid angle = "<<angsol<<" sr = "<<angsol/(4.*M_PI)<<" *4Pi sr"<<endl; … … 241 274 242 275 // --- 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; 244 277 cout<<"Array size: nx = "<<nx<<", ny = "<<ny<<", nz = "<<nz<<endl; 245 278 double dk_x = 2.*M_PI/(nx*dx), knyq_x = M_PI/dx; … … 254 287 255 288 // --- Masse de HI 256 cout<<"\n --- Mass HI"<<endl;289 cout<<"\n>>>>\n>>>> Mass HI\n>>>>"<<endl; 257 290 Schechter sch(nstar,mstar,alpha); 258 291 sch.SetOutValue(1); … … 271 304 272 305 // --- Survey values 306 cout<<"\n>>>>\n>>>> Observations\n>>>>"<<endl; 273 307 double unplusz = 1.+redshift; 274 308 double nuhiz = Fr_HyperFin_Par / unplusz; // GHz … … 278 312 double dnuhiz = Fr_HyperFin_Par *dred/(unplusz*unplusz) 279 313 / (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 282 315 <<" nu="<<nuhiz<<" GHz, dnu="<<dnuhiz*1.e3<<" Mhz"<<endl; 283 316 cout<<"dang lumi = "<<dlum<<" in ["<<dlumlim[0]<<","<<dlumlim[1]<<"] Mpc"<<endl; … … 330 363 cout<<"Warning: doppler width "<<doplarge<<" GHz > "<<dnuhiz<<" GHz redshift bin width"<<endl; 331 364 332 // Synchrotron 365 // Synchrotron (T en -2.7 -> Flux en -0.7 dans l'approximation Rayleigh) 333 366 double tsynch = Tsynch408; 334 367 if(fabs(indnu)>1.e-50) tsynch *= pow(nuhiz/nuhaslam,indnu);
Note:
See TracChangeset
for help on using the changeset viewer.