Changeset 3365 in Sophya


Ignore:
Timestamp:
Oct 29, 2007, 7:21:41 PM (18 years ago)
Author:
cmv
Message:

1./ grosses modifs de structure pour cmvdefsurv.cc
2./ PoissRandLimit enleve de geneutils.{h,cc} et mis dans srandgen.{h,cc}

cmv 29/10/2007

Location:
trunk/Cosmo/SimLSS
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/Cosmo/SimLSS/cmvdefsurv.cc

    r3352 r3365  
    2121 --- */
    2222
     23//-----------------------------------------------------------------------------------------------------------
    2324inline double rad2deg(double trad) {return trad/M_PI*180.;}
    2425inline double rad2min(double trad) {return trad/M_PI*180.*60.;}
    … …  
    2829inline double sec2rad(double tsec) {return tsec*M_PI/(180.*3600.);}
    2930
     31double LargeurDoppler(double v /* km/s */, double nu);
     32double DzFrV(double v /* km/s */, double zred);
     33double DNuFrDz(double dzred,double nu_at_0,double zred);
     34double DzFrDNu(double dnu_at_0,double nu_at_0,double zred);
     35double DzFrDNuApprox(double dnu_at_0,double nu_at_0,double zred);
     36double ZFrLos(double loscom /* Mpc com */,CosmoCalc& univ);
     37double AngsolEqTelescope(double nu /* GHz */,double telsurf /* m^2 */);
    3038void usage(void);
     39
    3140void usage(void) {
    32   cout<<"cmvdefsurv [options] -x adtx,atxlarg[,unit_x] -y adty,atylarg[,unit_y] -z dred,redlarg[,unit_z] redshift"<<endl
     41  cout<<"cmvdefsurv [options] -x dx,txlarg[,unit_x] -y dy,tylarg[,unit_y] -z dz,tzlarg[,unit_z] redshift"<<endl
    3342      <<"----------------"<<endl
    34       <<" -x adtx,atxlarg : resolution et largeur dans le plan transverse selon X"<<endl
    35       <<" -y adty,atylarg : idem selon Y, si <=0 meme que X"<<endl
    36       <<" -z dred,redlarg : resolution et largeur sur la ligne de visee"<<endl
     43      <<" -x dx,txlarg[,unit_x] : resolution et largeur dans le plan transverse selon X"<<endl
     44      <<" -y dy,tylarg[,unit_y] : idem selon Y, si <=0 meme que X"<<endl
     45      <<" -z dz,tzlarg[,unit_z] : resolution et largeur sur la ligne de visee"<<endl
    3746      <<"-- Unites pour X-Y:"<<endl
    38       <<"   \'A\' : en angles (pour X-Y)     : resolution=ArcMin, largeur=Degre (defaut)"<<endl
    39       <<"   \'Z\' : en redshift (pour Z)     : resolution et largeur en redshift (defaut)"<<endl
    40       <<"   \'F\' : en frequence (pour Z)    : resolution et largeur MHz"<<endl
    41       <<"   \'M\' : en distance (pour X-Y-Z) : resolution et largeur Mpc"<<endl
     47      <<"   \'A\' : en angles (pour X-Y)  : resolution=ArcMin, largeur=Degre (defaut)"<<endl
     48      <<"   \'Z\' : en redshift (pour Z)  : resolution et largeur en redshift (defaut)"<<endl
     49      <<"   \'F\' : en frequence (pour Z) : resolution et largeur MHz"<<endl
     50      <<"   \'M\' : en longeur comobile (pour X-Y-Z) : resolution et largeur Mpc"<<endl
    4251      <<"----------------"<<endl
    4352      <<" -K k,dk,pk : k(Mpc^-1) dk(Mpc^-1) pk(a z=0 en Mpc^-3) pour estimer la variance cosmique"<<endl
    … …  
    4554      <<" -O surf,tobs,eta : surface effective (m^2), temps d\'observation (s), efficacite d\'ouverture"<<endl
    4655      <<" -N Tsys : temperature du system (K)"<<endl
    47       <<" -L lobewidth,freqlob : taille du lobe d\'observation (FWHM) en arcmin (def= 1\')"<<endl
     56      <<" -L lobewidth,freqlob : taille du lobe d\'observation (FWHM) en arcmin (def= rien)"<<endl
    4857      <<"                        pour la frequence freqlob en MHz"<<endl
    49       <<"            Si lobewidth<=0 : l'angle solide du lobe = celui du pixel"<<endl
    5058      <<"            Si freqlob<=0 : la frequence de reference est celle du redshift etudie"<<endl
    5159      <<"            Si freqlob absent : la frequence de reference 1.4 GHz"<<endl
     60      <<"            lobewidth est pris comme la FWHM du lobe gaussien equivalent"<<endl
    5261      <<" -2 : two polarisations measured"<<endl
    5362      <<" -M  : masse de HI de reference (MSol), si <=0 mean schechter in pixel"<<endl
    … …  
    6574}
    6675
     76//-------------------------------------------------------------------------------------------
    6777int main(int narg,char *arg[])
    6878{
    69  // --- Valeurs fixes
    70  // WMAP
     79 // ---
     80 // --- Valeurs fixes ou par defaut
     81 // ---
     82
     83 //-- WMAP
    7184 unsigned short flat = 0;
    7285 double h100=0.71, om0=0.267804, or0=7.9e-05*0., ol0=0.73,w0=-1.;
    73  // Schechter
     86
     87 //-- Schechter HIMASS
    7488 double h75 = h100 / 0.75;
    7589 double nstar = 0.006*pow(h75,3.);  //
    … …  
    7892 cout<<"nstar= "<<nstar<<"  mstar="<<mstar<<"  alpha="<<alpha<<endl;
    7993
    80  // --- Arguments
    81  double adtx=0., atxlarg=0., dx=0.,txlarg=0.;
    82    int nx=0; char unit_x = 'A';
    83  double adty=-1., atylarg=-1., dy=0.,tylarg=0.;
    84    int ny=0; char unit_y = 'A';
    85  double dred=0., redlarg=0., dz=0.,tzlarg=0.;
    86    int nz=0; char unit_z = 'Z';
    87  double redshift = 0.;
     94 double hifactor = 1.;
     95 double vrot = 300.; // largeur en vitesse (km/s) pour elargissement doppler
     96
     97 //-- CMB , Synchrotron et AGN
     98 double tcmb = T_CMB_Par;
     99 // a 408 MHz (Haslam) + evol index a -2.6
     100 double Tsynch408=60., nuhaslam=0.408, indnu = -2.6;
     101 double lflux_agn = -3.;
     102
     103 //-- Appareillage
     104 bool ya2polar = false;
     105 double facpolar = 0.5; // si on mesure les 2 polars -> 1.0
     106 double Tsys=75.;
    88107 double tobs = 6000., surftot = 400000., eta_eff  = 1.;
    89  // variance cosmique (default = standard SDSSII)
    90  double kcosm = 0.05, dkcosm = -1., pkcosm = 40000.;
    91108 // taille du lobe d'observation en arcmin pour la frequence
    92109 double lobewidth0 = -1., lobefreq0 = Fr_HyperFin_Par*1.e3;
    93  double Tsys=75.;
    94  // a 408 MHz (Haslam) + evol index a -2.6
    95  double Tsynch408=60., nuhaslam=0.408, indnu = -2.6;
     110
     111 //-- Variance cosmique (default = standard SDSSII)
     112 double kcosm = 0.05, dkcosm = -1., pkcosm = 40000.;
     113
     114 // --- Pour taille du survey
     115 double dx=0., dy=-1., dz=0., txlarg=0., tylarg=-1., tzlarg=0.;
     116 int nx=0, ny=0, nz=0;
     117 char unit_x = 'A', unit_y = 'A', unit_z = 'Z';
     118 double redshift0 = 0.;
    96119 double mhiref = -1.; // reference Mass en HI (def integ schechter)
    97  double hifactor = 1.;
    98  double vrot = 300.; // largeur en vitesse (km/s) pour elargissement doppler
    99  bool ya2polar = false;
    100  double facpolar = 0.5; // si on mesure les 2 polars -> 1.0
    101  double lflux_agn = -3.;
    102 
     120
     121 // ---
    103122 // --- Decodage arguments
     123 // ---
    104124 char c;
    105125  while((c = getopt(narg,arg,"h2x:y:z:N:S:O:M:F:V:U:L:A:K:")) != -1) {
    106126  switch (c) {
    107127  case 'x' :
    108     sscanf(optarg,"%lf,%lf,%c",&adtx,&atxlarg,&unit_x);
     128    sscanf(optarg,"%lf,%lf,%c",&dx,&txlarg,&unit_x);
    109129    break;
    110130  case 'y' :
    111     sscanf(optarg,"%lf,%lf,%c",&adty,&atylarg,&unit_y);
     131    sscanf(optarg,"%lf,%lf,%c",&dy,&tylarg,&unit_y);
    112132    break;
    113133  case 'z' :
    114     sscanf(optarg,"%lf,%lf,%c",&dred,&redlarg,&unit_z);
     134    sscanf(optarg,"%lf,%lf,%c",&dz,&tzlarg,&unit_z);
    115135    break;
    116136  case 'O' :
    … …  
    153173  }
    154174 }
    155  if(optind>=narg) {usage(); return-1;}
    156  sscanf(arg[optind],"%lf",&redshift);
    157  if(redshift<=0.) {cout<<"Redshift "<<redshift<<" should be >0"<<endl; return -2;}
    158 
     175 if(optind>=narg) {usage(); return -1;}
     176 sscanf(arg[optind],"%lf",&redshift0);
     177 if(redshift0<=0.) {cout<<"Redshift "<<redshift0<<" should be >0"<<endl; return -2;}
     178
     179 // ---
    159180 // --- Initialisation de la Cosmologie
    160  cout<<"\n>>>>\n>>>> Cosmologie generale\n>>>>"<<endl;
    161  cout<<"h100="<<h100<<" Om0="<<om0<<" Or0="<<or0<<" Or0="
    162      <<or0<<" Ol0="<<ol0<<" w0="<<w0<<" flat="<<flat<<endl;
    163  cout<<"--- Cosmology for z = "<<redshift<<endl;
    164  CosmoCalc univ(flat,true,2.*redshift);
    165  double perc=0.01,dzinc=redshift/100.,dzmax=dzinc*10.; unsigned short glorder=4;
     181 // ---
     182 CosmoCalc univ(flat,true,2.*redshift0);
     183 double perc=0.01,dzinc=redshift0/100.,dzmax=dzinc*10.; unsigned short glorder=4;
    166184 univ.SetInteg(perc,dzinc,dzmax,glorder);
    167185 univ.SetDynParam(h100,om0,or0,ol0,w0);
     186
     187 GrowthFactor growth(om0,ol0);
     188
     189 // ---
     190 // --- Mise en forme dans les unites appropriees
     191 // ---
     192 // ATTENTION: le cube de simulation est en Mpc avec des pixels de taille comobile fixe
     193 cout<<"\n>>>>\n>>>> Geometrie\n>>>>"<<endl;
     194 if(dy<=0. || tylarg<=0.) {dy=dx; tylarg=txlarg; unit_y=unit_x;}
     195 cout<<"X values: resolution="<<dx<<" largeur="<<txlarg<<" unite="<<unit_x<<endl;
     196 if(unit_x == 'A') {
     197   nx = int(txlarg*60./dx+0.5);
     198   dx = min2rad(dx)*univ.Dtrcom(redshift0);
     199   txlarg = dx*nx;
     200 } else if(unit_x == 'M') {
     201   nx = int(txlarg/dx+0.5);
     202   txlarg = dx*nx;
     203 } else {
     204   cout<<"Unknown unit_x = "<<unit_x<<endl; return -2;
     205 }
     206 cout<<"Y values: resolution="<<dy<<" largeur="<<tylarg<<" unite="<<unit_y<<endl;
     207 if(unit_y == 'A') {
     208   ny = int(tylarg*60./dy+0.5);
     209   dy = min2rad(dy)*univ.Dtrcom(redshift0);
     210   tylarg = dy*ny;
     211 } else if(unit_y == 'M') {
     212   ny = int(tylarg/dy+0.5);
     213   tylarg = dy*ny;
     214 } else {
     215   cout<<"Unknown unit_y = "<<unit_y<<endl; return -2;
     216 }
     217 cout<<"Z values: resolution="<<dz<<" largeur="<<tzlarg<<" unite="<<unit_z<<endl;
     218 if(unit_z == 'Z') {
     219   nz = int(tzlarg/dz+0.5);
     220   dz = dz*univ.Dhubble()/univ.E(redshift0);
     221   tzlarg = dz*nz;
     222 } else if(unit_z == 'M') {
     223   nz = int(tzlarg/dz+0.5);
     224   tzlarg = dz*nz;
     225 } else if(unit_z == 'F') { // unite en MHz
     226   nz = int(tzlarg/dz+0.5);
     227   dz = DzFrDNu(dz,Fr_HyperFin_Par*1.e3,redshift0);  //dzred
     228   dz = dz*univ.Dhubble()/univ.E(redshift0);
     229   tzlarg = dz*nz;
     230 } else {
     231   cout<<"Unknown unit_z = "<<unit_z<<endl; return -2;
     232 }
     233
     234 // ---
     235 // --- Calcul des valeurs au centre du cube et aux faces limites av/ar
     236 // ---
     237 cout<<"\n>>>>\n>>>> Cosmologie generale\n>>>>"<<endl;
     238 double adtx[3], adty[3], atxlarg[3], atylarg[3];
     239 double loscom[3], redshift[3], unplusz[3], dred[3], growthfac[3], rhocz[3];
     240 double dang[3], dtrcom[3], dlum[3], dloscom[3], dlosdz[3];
     241 double nuhiz[3], lambdahiz[3], dnuhiz[3];
     242 double angsol_pix[3], angsol_tot[3];
     243
     244 redshift[0] = redshift0;
     245 loscom[0] = univ.Dloscom(redshift0);
     246 loscom[1] = loscom[0]-tzlarg/2.;
     247 loscom[2] = loscom[0]+tzlarg/2.;
     248 //if(loscom[1]<=0.) {cout<<"Lower distance limit "<<loscom[1]<<" should be >0"<<endl; return -3;}
     249 for(int i=0;i<3;i++) {
     250   double l = loscom[i]; if(l<=0.) l = dz/2.;
     251   redshift[i]  = ZFrLos(l,univ);
     252   unplusz[i]   = 1. + redshift[i];
     253   growthfac[i] = growth(redshift[i]);
     254   rhocz[i]     = univ.Rhoc(redshift[i])*GCm3toMsolMpc3_Cst;;
     255   dang[i]      = univ.Dang(redshift[i]);
     256   dtrcom[i]    = univ.Dtrcom(redshift[i]);
     257   dlum[i]      = univ.Dlum(redshift[i]);
     258   dloscom[i]   = univ.Dloscom(redshift[i]);
     259   dlosdz[i]    = univ.Dhubble()/univ.E(redshift[i]);
     260   dred[i]      = dz/dlosdz[i];
     261   adtx[i]      = dx/dtrcom[i];
     262   atxlarg[i]   = adtx[i]*nx;
     263   adty[i]      = dy/dtrcom[i];
     264   atylarg[i]   = adty[i]*ny;
     265   nuhiz[i]     = Fr_HyperFin_Par / unplusz[i];  // GHz
     266   lambdahiz[i] = SpeedOfLight_Cst*1000./(nuhiz[i]*1.e9);  // m
     267   dnuhiz[i]    = DNuFrDz(dred[i],Fr_HyperFin_Par,redshift[i]);  // GHz
     268   angsol_pix[i] = AngSol(adtx[i]/2.,adty[i]/2.,M_PI/2.);
     269   angsol_tot[i] = AngSol(atxlarg[i]/2.,atylarg[i]/2.,M_PI/2.);
     270 }
     271 
     272
     273 cout<<"--- Cosmology for z = "<<0.<<endl;
    168274 univ.Print(0.);
    169  univ.Print(redshift);
    170  GrowthFactor growth(om0,ol0);
    171  double growthfac = growth(redshift);
    172  cout<<"Facteur de croissance lineaire = "<<growthfac
    173      <<" ^2 , ( 1/(1+z)="<<1./(1.+redshift)<<" )"<<endl;
    174275 double rhoc0 = univ.Rhoc(0.)*GCm3toMsolMpc3_Cst;
    175  double rhocz = univ.Rhoc(redshift)*GCm3toMsolMpc3_Cst;
    176  cout<<"Rho_c a z=0: "<<rhoc0<<", a z="<<redshift<<": "<<rhocz<<" Msol/Mpc^3"<<endl;
    177 
    178  double dang    = univ.Dang(redshift);
    179  double dtrcom  = univ.Dtrcom(redshift);
    180  double dlum    = univ.Dlum(redshift);
    181  double dloscom = univ.Dloscom(redshift);
    182  double dlosdz  = univ.Dhubble()/univ.E(redshift);
    183  cout<<"dang="<<dang<<" dlum="<<dlum<<" dtrcom="<<dtrcom
    184      <<" dloscom="<<dloscom<<" dlosdz="<<dlosdz<<" Mpc"<<endl;
    185 
    186  cout<<"\n1\" -> "<<dang*sec2rad(1.)<<" Mpc = "<<dtrcom*sec2rad(1.)<<" Mpc com"<<endl;
    187  cout<<"1\' -> "<<dang*min2rad(1.)<<" Mpc = "<<dtrcom*min2rad(1.)<<" Mpc com"<<endl;
    188  cout<<"1d -> "<<dang*deg2rad(1.)<<" Mpc = "<<dtrcom*deg2rad(1.)<<" Mpc com"<<endl;
    189 
    190  cout<<"dz=1 -> "<<dlosdz<<" Mpc com"<<endl;
    191 
    192  cout<<"1 Mpc los com -> dz = "<<1./dlosdz<<endl;
    193  cout<<"1 Mpc transv com -> "<<rad2sec(1./dtrcom)<<"\" = "
    194      <<rad2min(1./dtrcom)<<" \' = "<<rad2deg(1./dtrcom)<<" d"<<endl;
    195 
    196  // --- Mise en forme dans les unites appropriees
    197  cout<<"\n>>>>\n>>>> Geometrie\n>>>>"<<endl;
    198  if(adty<=0. || atylarg<=0.) {adty=adtx; atylarg=atxlarg; unit_y=unit_x;}
    199  cout<<"X values: resolution="<<adtx<<" largeur="<<atxlarg<<" unite="<<unit_x<<endl;
    200  if(unit_x == 'A') {
    201    nx = int(atxlarg*60./adtx+0.5);
    202    adtx = min2rad(adtx); atxlarg = deg2rad(atxlarg);
    203    dx = adtx*dtrcom; txlarg = dx*nx;
    204  } else if(unit_x == 'M') {
    205    nx = int(atxlarg/adtx+0.5);
    206    dx = adtx; txlarg = atxlarg;
    207    adtx = dx/dtrcom; atxlarg = adtx*nx;
    208  } else {
    209    cout<<"Unknown unit_x = "<<unit_x<<endl;
    210  }
    211  cout<<"Y values: resolution="<<adty<<" largeur="<<atylarg<<" unite="<<unit_y<<endl;
    212  if(unit_y == 'A') {
    213    ny = int(atylarg*60./adty+0.5);
    214    adty = min2rad(adty); atylarg = deg2rad(atylarg);
    215    dy = adty*dtrcom; tylarg = dy*ny;
    216  } else if(unit_y == 'M') {
    217    ny = int(atylarg/adty+0.5);
    218    dy = adty; tylarg = atylarg;
    219    adty = dy/dtrcom; atylarg = adty*ny;
    220  } else {
    221    cout<<"Unknown unit_y = "<<unit_y<<endl;
    222  }
    223  cout<<"Z values: resolution="<<dred<<" largeur="<<redlarg<<" unite="<<unit_z<<endl;
    224  if(unit_z == 'Z') {
    225    nz = int(redlarg/dred+0.5);
    226    dz = dred*dlosdz; tzlarg = dz*nz;
    227  } else if(unit_z == 'M') {
    228    nz = int(redlarg/dred+0.5);
    229    dz = dred; tzlarg = redlarg;
    230    dred = dz/dlosdz; redlarg = dred*nz;
    231  } else if(unit_z == 'F') {
    232    nz = int(redlarg/dred+0.5);
    233    dred = dred/(Fr_HyperFin_Par*1.e3)*pow(1.+redshift,2.); redlarg = dred*nz;
    234    dz = dred*dlosdz; tzlarg = dz*nz;
    235  } else {
    236    cout<<"Unknown unit_z = "<<unit_z<<endl;
    237  }
    238 
     276
     277 for(int i=0;i<3;i++) {
     278   cout<<"\n--- Cosmology for z = "<<redshift[i]<<endl;
     279   univ.Print(redshift[i]);
     280 }
     281
     282 cout<<endl;
     283 cout<<"Facteur de croissance lineaire = "<<growthfac[0]
     284     <<" in ["<<growthfac[1]<<","<<growthfac[2]<<"]"<<endl;
     285 cout<<"                       1/(1+z) = "<<1./(1.+redshift[0])
     286     <<" in ["<<1./(1.+redshift[1])<<","<<1./(1.+redshift[2])<<"]"<<endl;
     287 cout<<"Facteur de croissance lineaire^2 = "<<growthfac[0]*growthfac[0]
     288     <<" in ["<<growthfac[1]*growthfac[1]<<","<<growthfac[2]*growthfac[2]<<"]"<<endl;
     289
     290 cout<<"Rho_c (z=0) =  "<<rhoc0<<", a z="<<0<<": "<<rhoc0<<" Msol/Mpc^3"<<endl;
     291 cout<<"Rho_c a z = "<<rhocz[0]
     292     <<" Msol/Mpc^3 in ["<<rhocz[1]<<","<<rhocz[2]<<"]"<<endl;
     293
     294 cout<<endl;
     295 cout<<"dang= "<<dang[0]<<" Mpc in ["<<dang[1]<<","<<dang[2]<<"]"<<endl;
     296 cout<<"dtrcom= "<<dtrcom[0]<<" Mpc com in ["<<dtrcom[1]<<","<<dtrcom[2]<<"]"<<endl;
     297 cout<<"dlum= "<<dlum[0]<<" Mpc in ["<<dlum[1]<<","<<dlum[2]<<"]"<<endl;
     298 cout<<"dloscom= "<<dloscom[0]<<" Mpc com in ["<<dloscom[1]<<","<<dloscom[2]<<"]"<<endl;
     299 cout<<"dz=1 -> dlosdz= "<<dlosdz[0]<<" Mpc com in ["<<dlosdz[1]<<","<<dlosdz[2]<<"]"<<endl;
     300 cout<<"1 Mpc los com -> dz = "<<1./dlosdz[0]<<" in ["<<1./dlosdz[1]<<","<<1./dlosdz[2]<<"]"<<endl;
     301
     302 cout<<endl;
     303 for(int i=0;i<3;i++) {
     304   cout<<"...Redshift="<<redshift[i]<<endl;
     305   cout<<"1\" -> "<<dang[i]*sec2rad(1.)<<" Mpc = "<<dtrcom[i]*sec2rad(1.)<<" Mpc com"<<endl;
     306   cout<<"1\' -> "<<dang[i]*min2rad(1.)<<" Mpc = "<<dtrcom[i]*min2rad(1.)<<" Mpc com"<<endl;
     307   cout<<"1d -> "<<dang[i]*deg2rad(1.)<<" Mpc = "<<dtrcom[i]*deg2rad(1.)<<" Mpc com"<<endl;
     308   cout<<"1 Mpc transv com -> "<<rad2sec(1./dtrcom[i])<<"\" = "
     309       <<rad2min(1./dtrcom[i])<<" \' = "<<rad2deg(1./dtrcom[i])<<" d"<<endl;
     310   cout<<"dz=1 -> dlosdz= "<<dlosdz[i]<<" Mpc com"<<endl;
     311   cout<<"1 Mpc los com -> dz = "<<1./dlosdz[0]<<endl;
     312 }
     313
     314 // ---
     315 // --- Cosmolographie Transverse
     316 // ---
     317 cout<<"\n>>>>\n>>>> Cosmologie & Geometrie transverse\n>>>>"<<endl;
     318
     319 cout<<"dx = "<<dx<<", txlarg = "<<txlarg<<" Mpc (com), nx="<<nx<<endl;
     320 cout<<"adtx = "<<adtx[0]<<" rd in ["<<adtx[1]<<","<<adtx[2]<<"]"<<endl;
     321 cout<<"     = "<<rad2min(adtx[0])<<"\' in ["<<rad2min(adtx[1])<<","<<rad2min(adtx[2])<<"]"<<endl;
     322 cout<<"atxlarg = "<<atxlarg[0]<<" rd in ["<<atxlarg[1]<<","<<atxlarg[2]<<"]"<<endl;
     323 cout<<"          "<<rad2deg(atxlarg[0])<<" d in ["<<rad2deg(atxlarg[1])<<","<<rad2deg(atxlarg[2])<<"]"<<endl;
     324
     325 if(fabs(dx-dy)>1.e-20 && fabs(txlarg-tylarg)>1.e-20) {
     326   cout<<"\ndy = "<<dy<<", tylarg = "<<tylarg<<" Mpc (com), ny="<<ny<<endl;
     327   cout<<"adty = "<<adty[0]<<" rd in ["<<adty[1]<<","<<adty[2]<<"]"<<endl;
     328   cout<<"     = "<<rad2min(adty[0])<<"\' in ["<<rad2min(adty[1])<<","<<rad2min(adty[2])<<"]"<<endl;
     329   cout<<"atylarg = "<<atylarg[0]<<" rd in ["<<atylarg[1]<<","<<atylarg[2]<<"]"<<endl;
     330   cout<<"          "<<rad2deg(atylarg[0])<<" d in ["<<rad2deg(atylarg[1])<<","<<rad2deg(atylarg[2])<<"]"<<endl;
     331 }
     332
     333 // ---
     334 // --- Cosmolographie Line of sight
     335 // ---
     336 cout<<"\n>>>>\n>>>> Cosmologie & Geometrie ligne de visee\n>>>>"<<endl;
     337 cout<<"dz = "<<dz<<", tzlarg = "<<tzlarg<<" Mpc (com), nz="<<nz<<endl;
     338 cout<<"Redshift = "<<redshift[0]<<" in ["<<redshift[1]<<","<<redshift[2]<<"]"<<endl;
     339 cout<<"dred = "<<dred[0]<<" in ["<<dred[1]<<","<<dred[2]<<"]"<<endl;
     340 cout<<"nu HI = "<<nuhiz[0]<<" GHz in ["<<nuhiz[1]<<","<<nuhiz[2]<<"]"<<endl;
     341 cout<<"lambda HI = "<<lambdahiz[0]<<" m in ["<<lambdahiz[1]<<","<<lambdahiz[2]<<"]"<<endl;
     342 cout<<"dnu HI= "<<dnuhiz[0]*1e3<<" MHz in ["<<dnuhiz[1]*1e3<<","<<dnuhiz[2]*1e3<<"]"<<endl;
     343
     344 // ---
     345 // --- Volume
     346 // ---
     347 cout<<"\n>>>>\n>>>> Volumes\n>>>>"<<endl;
    239348 double Npix = (double)nx*(double)ny*(double)nz;
    240  double redlim[2] = {redshift-redlarg/2.,redshift+redlarg/2.};
    241  if(redlim[0]<=0.)
    242    {cout<<"Lower redshift limit "<<redlim[0]<<" should be >0"<<endl; return -3;}
    243  double dtrlim[2]  = {univ.Dtrcom(redlim[0]),univ.Dtrcom(redlim[1])};
    244  double loslim[2]  = {univ.Dloscom(redlim[0]), univ.Dloscom(redlim[1])};
    245  double dlumlim[2] = {univ.Dlum(redlim[0]),univ.Dlum(redlim[1])};
    246 
    247  cout<<"---- Line of Sight: Redshift = "<<redshift<<endl
    248      <<"dred = "<<dred<<"  redlarg = "<<redlarg<<endl
    249      <<"  dz = "<<dz<<" Mpc  redlarg = "<<tzlarg<<" Mpc com, nz = "<<nz<<" pix"<<endl;
    250  cout<<"---- Transverse X:"<<endl
    251      <<"adtx = "<<rad2min(adtx)<<"\',  atxlarg = "<<rad2deg(atxlarg)<<" d"<<endl
    252      <<"  dx = "<<dx<<" Mpc,  txlarg = "<<txlarg<<" Mpc com, nx = "<<nx<<" pix"<<endl;
    253  cout<<"---- Transverse Y:"<<endl
    254      <<"adty = "<<rad2min(adty)<<"\',  atylarg = "<<rad2deg(atylarg)<<" d"<<endl
    255      <<"  dy = "<<dy<<" Mpc,  tylarg = "<<tylarg<<" Mpc com, ny = "<<ny<<" pix"<<endl;
    256  cout<<"---- Npix total = "<<Npix<<" -> "<<Npix*sizeof(double)/1.e6<<" Mo"<<endl;
    257  cout<<"     Volume pixel = "<<dx*dy*dz<<" Mpc^3"<<endl;
    258  cout<<"     Volume total = "<<Npix*dx*dy*dz<<" Mpc^3"<<endl;
    259 
    260  // --- Cosmolographie Transverse
    261  cout<<"\n>>>>\n>>>> Cosmologie & Geometrie transverse\n>>>>"<<endl;
    262  cout<<"dang comoving = "<<dtrcom<<" Mpc (com) var_in_z ["
    263                          <<dtrlim[0]<<","<<dtrlim[1]<<"]"<<endl;
    264 
    265  cout<<"... dx = "<<dx<<" Mpc (com), with angle "<<adtx*dtrcom<<endl
    266      <<"       with angle var_in_z ["<<adtx*dtrlim[0]<<","<<adtx*dtrlim[1]<<"]"<<endl;
    267  cout<<"... largx = "<<txlarg<<" Mpc (com), with angle "<<atxlarg*dtrcom<<endl
    268      <<"          with angle var_in_z ["<<atxlarg*dtrlim[0]<<","<<atxlarg*dtrlim[1]<<"]"<<endl;
    269 
    270  cout<<"... dy = "<<dy<<" Mpc (com), with angle "<<adty*dtrcom<<endl
    271      <<"       with angle var_in_z ["<<adty*dtrlim[0]<<","<<adty*dtrlim[1]<<"]"<<endl;
    272  cout<<"... largy = "<<tylarg<<" Mpc (com), with angle "<<atylarg*dtrcom<<endl
    273      <<"          with angle var_in_z ["<<atylarg*dtrlim[0]<<","<<atylarg*dtrlim[1]<<"]"<<endl;
    274 
    275  // --- Cosmolographie Line of sight
    276  cout<<"\n>>>>\n>>>> Cosmologie & Geometrie ligne de visee\n>>>>"<<endl;
    277  cout<<"los comoving distance = "<<dloscom<<" Mpc (com) in ["
    278      <<loslim[0]<<","<<loslim[1]<<"]"<<endl
    279      <<"                                    diff = "
    280      <<loslim[1]-loslim[0]<<" Mpc"<<endl;
    281 
    282  cout<<"...dz = "<<dz<<" Mpc (com), with redshift approx "<<dred*dlosdz<<endl;
    283  cout<<"...tzlarg = "<<tzlarg<<" Mpc (com), with redshift approx "<<redlarg*dlosdz<<endl;
    284 
    285  // --- Solid Angle & Volume
    286  cout<<"\n>>>>\n>>>> Angles solides et Volumes\n>>>>"<<endl;
    287  cout<<"--- Solid angle"<<endl;
    288  double angsol = AngSol(adtx/2.,adty/2.,M_PI/2.);
    289  cout<<"Elementary solid angle = "<<angsol<<" sr = "<<angsol/(4.*M_PI)<<" *4Pi sr"<<endl;
    290  double angsoltot = AngSol(atxlarg/2.,atylarg/2.,M_PI/2.);
    291  cout<<"Total solid angle = "<<angsoltot<<" sr = "<<angsoltot/(4.*M_PI)<<" *4Pi sr"<<endl;
    292 
    293  cout<<"\n--- Volume"<<endl;
    294  double dvol = dx*dy*dz;
    295  cout<<"Pixel volume comoving = "<<dvol<<" Mpc^3"<<endl;
    296  double vol = univ.Vol4Pi(redlim[0],redlim[1])/(4.*M_PI)*angsoltot;
    297  cout<<"Volume comoving = "<<vol<<" Mpc^3 = "<<vol/1.e9<<" Gpc^3"<<endl
    298      <<"Pixel volume comoving = vol/Npix = "<<vol/Npix<<" Mpc^3"<<endl;
    299 
    300  // --- Fourier space: k = omega = 2*Pi*Nu
     349 double vol_pixel = dx*dy*dz;
     350 double vol_survey = vol_pixel*Npix;
     351 cout<<"Npix total = "<<Npix<<" -> "<<Npix*sizeof(double)/1.e6<<" Mo"<<endl;
     352 cout<<"Volume pixel = "<<vol_pixel<<" Mpc^3 com"<<endl;
     353 cout<<"Volume total = "<<vol_survey<<" Mpc^3 com = "<<vol_survey/1.e9<<" Gpc^3"<<endl;
     354
     355 double vol = univ.Vol4Pi(redshift[1],redshift[2])/(4.*M_PI)*angsol_tot[0];
     356 cout<<"Calcul avec angsol et redshift: vol = "<<vol<<" Mpc^3 = "<<vol/1.e9<<" Gpc^3"<<endl
     357     <<"       -> pixel volume comoving = vol/Npix = "<<vol/Npix<<" Mpc^3"<<endl;
     358
     359 // ---
     360 // --- Angles solides
     361 // ---
     362 cout<<"\n>>>>\n>>>> Angles solides\n>>>>"<<endl;
     363 cout<<"Pixel solid angle = "<<angsol_pix[0]<<" sr ("<<angsol_pix[0]/(4.*M_PI)<<" *4Pi)"
     364     <<" in ["<<angsol_pix[1]<<","<<angsol_pix[2]<<"]"<<endl;
     365 cout<<"Total solid angle = "<<angsol_tot[0]<<" sr ("<<angsol_tot[0]/(4.*M_PI)<<" *4Pi)"
     366     <<" in ["<<angsol_tot[1]<<","<<angsol_tot[2]<<"]"<<endl;
     367
     368 // ---
     369 // --- Fourier space: k = omega = 2*Pi*Nu = 2*Pi*C/Lambda
     370 // ---
    301371 cout<<"\n>>>>\n>>>> Geometrie dans l'espace de Fourier\n>>>>"<<endl;
    302372 cout<<"Array size: nx = "<<nx<<",  ny = "<<ny<<",  nz = "<<nz<<endl;
    … …  
    304374 double dk_y = 2.*M_PI/(nx*dy), knyq_y = M_PI/dy;
    305375 double dk_z = 2.*M_PI/(nz*dz), knyq_z = M_PI/dz;
    306  cout<<"Resolution: dk_x = "<<dk_x<<" Mpc^-1  (2Pi/dk_x="<<2.*M_PI/dk_x<<" Mpc)"<<endl
    307      <<"            dk_y = "<<dk_y<<" Mpc^-1  (2Pi/dk_y="<<2.*M_PI/dk_y<<" Mpc)"<<endl;
    308  cout<<"Nyquist:    kx = "<<knyq_x<<" Mpc^-1  (2Pi/knyq_x="<<2.*M_PI/knyq_x<<" Mpc)"<<endl
    309      <<"            ky = "<<knyq_y<<" Mpc^-1  (2Pi/knyq_y="<<2.*M_PI/knyq_y<<" Mpc)"<<endl;
    310  cout<<"Resolution: dk_z = "<<dk_z<<" Mpc^-1  (2Pi/dk_z="<<2.*M_PI/dk_z<<" Mpc)"<<endl;
    311  cout<<"Nyquist:    kz = "<<knyq_z<<" Mpc^-1  (2Pi/knyq_z="<<2.*M_PI/knyq_z<<" Mpc)"<<endl;
    312 
     376 cout<<"Transverse:"<<endl
     377     <<"  Resolution: dk_x = "<<dk_x<<" Mpc^-1  (2Pi/dk_x="<<2.*M_PI/dk_x<<" Mpc)"<<endl
     378     <<"              dk_y = "<<dk_y<<" Mpc^-1  (2Pi/dk_y="<<2.*M_PI/dk_y<<" Mpc)"<<endl
     379     <<"  Nyquist:    kx = "<<knyq_x<<" Mpc^-1  (2Pi/knyq_x="<<2.*M_PI/knyq_x<<" Mpc)"<<endl
     380     <<"              ky = "<<knyq_y<<" Mpc^-1  (2Pi/knyq_y="<<2.*M_PI/knyq_y<<" Mpc)"<<endl;
     381 cout<<"Line of sight:"<<endl
     382     <<"  Resolution: dk_z = "<<dk_z<<" Mpc^-1  (2Pi/dk_z="<<2.*M_PI/dk_z<<" Mpc)"<<endl
     383     <<"  Nyquist:    kz = "<<knyq_z<<" Mpc^-1  (2Pi/knyq_z="<<2.*M_PI/knyq_z<<" Mpc)"<<endl;
     384
     385 // ---
    313386 // --- Variance cosmique
     387 // ---
    314388 //                      cosmique                    poisson
    315389 // (sigma/P)^2 = 2*(2Pi)^3 / (4Pi k^2 dk Vsurvey) * [(1+n*P)/(n*P)]^2
    316390 // nombre de mode = 1/2 * Vsurvey/(2Pi)^3 * 4Pi*k^2*dk
    317391 if(kcosm>0. && pkcosm>0.) {
    318    double pk = pkcosm*pow(growthfac,2.);
     392   double pk = pkcosm*pow(growthfac[0],2.);
    319393   cout<<"\n>>>>\n>>>> variance cosmique pour k="<<kcosm<<" Mpc^-1, pk(z=0)="
    320        <<pkcosm<<", pk(z="<<redshift<<")="<<pk<<"\n>>>>"<<endl;
     394       <<pkcosm<<", pk(z="<<redshift[0]<<")="<<pk<<"\n>>>>"<<endl;
    321395   for(int i=0;i<3;i++) {  // la correction de variance du au bruit de poisson
    322396     double v = 1.1; if(i==1) v=1.5; if(i==2) v=2.0;
    … …  
    324398     cout<<"  pour "<<ngal<<" gal/Mpc^3 on multiplie par "<<v<<" sigma/P"<<endl;
    325399   }
     400
     401   cout<<endl;
    326402   vector<double> dk; if(dkcosm>0.) dk.push_back(dkcosm);
    327403   dk.push_back(dk_x); dk.push_back(dk_y); dk.push_back(dk_z);
    328404   for(int i=0;i<(int)dk.size();i++) { // la variance cosmique pure
    329      double vcosm = sqrt( 2.*pow(2.*M_PI,3.)/(4.*M_PI*pow(kcosm,2.)*dk[i]*vol) );
    330      double nmode = 0.5*vol/pow(2.*M_PI,3.) * 4.*M_PI*pow(kcosm,2.)*dk[i];
     405     double vcosm = sqrt( 2.*pow(2.*M_PI,3.)/(4.*M_PI*pow(kcosm,2.)*dk[i]*vol_survey) );
     406     double nmode = 0.5*vol_survey/pow(2.*M_PI,3.) * 4.*M_PI*pow(kcosm,2.)*dk[i];
    331407     cout<<"  pour dk = "<<dk[i]<<" Mpc^-1:  sigma/P = "<<vcosm
    332408         <<" , Nmode = "<<nmode<<endl;
    … …  
    334410 }
    335411
     412 // ---
    336413 // --- Masse de HI
     414 // ---
    337415 cout<<"\n>>>>\n>>>> Mass HI\n>>>>"<<endl;
    338416 Schechter sch(nstar,mstar,alpha);
    … …  
    345423 cout<<"Mass density: "<<masshimpc3<<" Msol/Mpc^3"<<endl;
    346424
    347  double masshipix = masshimpc3*dvol;
    348  double masshitot = masshimpc3*vol;
     425 double masshipix = masshimpc3*vol_pixel;
     426 double masshitot = masshimpc3*vol_survey;
    349427 cout<<"Pixel mass = "<<masshipix<<" Msol"<<endl
    350428     <<"Total mass in survey = "<<masshitot<<" Msol"<<endl;
    351  cout<<"OmegaHI a z=0: "<<masshimpc3/rhoc0
    352      <<", a z="<<redshift<<": "<<masshimpc3/rhocz<<endl;
     429 cout<<"OmegaHI a z=0: "<<masshimpc3/rhoc0<<endl;
     430 for(int i=0;i<3;i++)
     431   cout<<"        a z="<<redshift[i]<<": "<<masshimpc3/rhocz[i]<<endl;
    353432 if(mhiref<=0.)  mhiref = masshipix;
    354433
    … …  
    358437 for(double x=lnx1; x<lnx2-0.5; x+=1.) {
    359438   double n = IntegrateFuncLog(sch,x,lnx2,0.001,dlnx,10.*dlnx,6);
    360    cout<<"  m>"<<x<<" Msol: "<<n<<" /Mpc^3, "<<n*dvol<<" /pixel, "
    361        <<n*vol<<" in survey"<<endl;
     439   cout<<"  m>"<<x<<" Msol: "<<n<<" /Mpc^3, "<<n*vol_pixel<<" /pixel, "
     440       <<n*vol_survey<<" in survey"<<endl;
    362441 }
    363442 sch.SetOutValue(1);
    364443
    365 
     444 // ---
    366445 // --- Survey values
     446 // ---
    367447 cout<<"\n>>>>\n>>>> Observations\n>>>>"<<endl;
    368448 double surfeff = surftot * eta_eff;
    369  double unplusz = 1.+redshift;
    370  double nuhiz = Fr_HyperFin_Par / unplusz;  // GHz
    371  // dnu = NuHi/(1.+z0-dz/2) - NuHi/(1.+z0+dz/2)
    372  //     = NuHi*dz/(1.+z0)^2 * 1/[1-(dz/(2*(1+z0)))^2]
    373  //    ~= NuHi*dz/(1.+z0)^2
    374  double dnuhiz = Fr_HyperFin_Par *dred/(unplusz*unplusz)
    375                / (1.- pow(dred/.2/unplusz,2.));
    376  cout<<"      surf="<<surftot<<" m^2, eta="<<eta_eff<<"  surf_eff="<<surfeff<<" m^2, tobs="<<tobs<<" s"<<endl
    377      <<"      nu="<<nuhiz<<" GHz, dnu="<<dnuhiz*1.e3<<" Mhz"<<endl;
    378  cout<<"dang lumi = "<<dlum<<" in ["<<dlumlim[0]<<","<<dlumlim[1]<<"] Mpc"<<endl;
    379 
    380  double nlobes = 1.;
     449 cout<<"surf_tot="<<surftot<<" m^2, eta="<<eta_eff<<"  surf_eff="<<surfeff<<" m^2, tobs="<<tobs<<" s"<<endl;
     450
     451 // Angles solides pour un telescope plein
     452 double angSEq[4], angEq[4];
     453 for(int i=0;i<4;i++) {
     454   double nu =Fr_HyperFin_Par;
     455   if(i<3) nu = nuhiz[i];
     456   angSEq[i] = AngsolEqTelescope(nu,surftot);
     457   angEq[i] = 2.*FrAngSol(angSEq[i]);
     458 }
     459 cout<<"\nFor a "<<surftot<<" m^2 telescope:"<<endl
     460     <<"    equivalent Omega = "<<angSEq[0]<<" sr in ["<<angSEq[1]<<","<<angSEq[2]<<"]"<<endl
     461     <<"    angular diameter = "<<angEq[0]<<" rd = "<<rad2min(angEq[0])
     462     <<"\' in ["<<angEq[2]<<","<<angEq[2]<<"]"<<endl;
     463 cout<<"At z=0: equivalent Omega = "<<angSEq[3]<<" sr"<<endl
     464     <<"        angular diameter = "<<angEq[3]<<" rd = "<<rad2min(angEq[3])<<"\'"<<endl;
     465
     466 // Pour une observation ou le lobe est + petit que le pixel de simulation
     467 // La taille angulaire du lobe change avec la frequence donc avec le redshift
     468 // La taille du lobe "lobewidth" est donnee pour une frequence de reference "lobefreq0"
     469 double nlobes[3] = {1.,1.,1.};
    381470 if(lobewidth0>0.) {
    382    double lobewidth = lobewidth0;   // ArcMin
    383    if(lobefreq0<=0.) lobefreq0 = nuhiz*1.e3; // MHz
    384    // La taille angulaire du lobe change avec la frequence donc avec le redshift
    385    lobewidth *= lobefreq0/(nuhiz*1.e3);
    386    cout<<"\n--- Lobe: width="<<lobewidth0<<" pour "<<lobefreq0<<" MHz"<<endl
    387        <<"          changed to "<<lobewidth<<" pour "<<nuhiz*1.e3<<" MHz"<<endl;
    388    double slobe = lobewidth/2.35482; // sigma du lobe en arcmin
    389    double lobecyl = sqrt(8.)*slobe; // diametre du lobe cylindrique equiv en arcmin
    390    double lobearea = M_PI*lobecyl*lobecyl/4.;  // en arcmin^2 (hypothese lobe gaussien)
    391    nlobes = rad2min(adtx)*rad2min(adty)/lobearea;
    392    cout<<"Beam FWHM = "<<lobewidth<<"\' -> sigma = "<<slobe<<"\' -> "
    393        <<" Dcyl = "<<lobecyl<<"\' -> area = "<<lobearea<<" arcmin^2"<<endl;
    394    cout<<"Number of beams in one transversal pixel = "<<nlobes<<endl;
    395  }
    396 
    397  // --- Power emitted by HI
    398  cout<<"\n--- Power from HI for M = "<<mhiref<<" Msol at "<<nuhiz<<" GHz"<<endl;
    399  cout<<"flux factor = "<<hifactor<<" at redshift = "<<redshift<<endl;
    400 
    401  double fhi = hifactor*Msol2FluxHI(mhiref,dlum);
    402  cout<<"FluxHI("<<dlum<<" Mpc) all polar:"<<endl
    403      <<"  Flux= "<<fhi<<" W/m^2 = "<<fhi/Jansky2Watt_cst<<" Jy.Hz"<<endl
    404      <<"      in ["<<hifactor*Msol2FluxHI(mhiref,dlumlim[0])
    405      <<","<<hifactor*Msol2FluxHI(mhiref,dlumlim[1])<<"] W/m^2"<<endl;
    406  double sfhi = fhi / (dnuhiz*1e9) / Jansky2Watt_cst;
    407  cout<<"If spread over pixel depth ("<<dnuhiz<<" GHz), flux density = "<<sfhi*1.e6<<" mu_Jy"<<endl;
    408 
     471   if(lobefreq0<=0.) lobefreq0 = nuhiz[0]*1.e3; // MHz
     472   cout<<"\n--- Lobe: width="<<lobewidth0<<"\' pour "<<lobefreq0<<" MHz"<<endl;
     473   double lobewidth[3], slobe[3], lobecyl[3], lobearea[3];
     474   for(int i=0;i<3;i++) {
     475     lobewidth[i] = lobewidth0*lobefreq0/(nuhiz[i]*1.e3); // arcmin
     476     slobe[i] = lobewidth[i]/2.35482; // sigma du lobe en arcmin
     477     lobecyl[i] = sqrt(8.)*slobe[i]; // diametre du lobe cylindrique equiv en arcmin
     478     lobearea[i] = M_PI*lobecyl[i]*lobecyl[i]/4.; // en arcmin^2 (hypothese lobe gaussien)
     479     nlobes[i] = rad2min(adtx[i])*rad2min(adty[i])/lobearea[i];
     480   }
     481   cout<<"Lobe width changed to "<<lobewidth[0]<<"\' pour "<<nuhiz[0]*1.e3<<" MHz"
     482       <<" in ["<<lobewidth[1]<<","<<lobewidth[2]<<"]"<<endl;
     483   cout<<"Lobe sigma "<<slobe[0]<<"\' in ["<<slobe[1]<<","<<slobe[2]<<"]"<<endl;
     484   cout<<"Lobe cylindrique diametre "<<lobecyl[0]<<"\' in ["<<lobecyl[1]<<","<<lobecyl[2]<<"]"<<endl;
     485   cout<<"Lobe cylindrique area "<<lobearea[0]<<"\'^2 in ["<<lobearea[1]<<","<<lobearea[2]<<"]"<<endl;
     486   cout<<"Number of beams in one transversal pixel "<<nlobes[0]<<" in ["<<nlobes[1]<<","<<nlobes[2]<<"]"<<endl;
     487 }
     488
     489 // ---
    409490 // --- Signal analysis
    410  cout<<"\n--- Signal analysis"<<endl;
     491 // ---
     492 // --- Temperature d'antenne Ta et temperature de brillance Tb
     493 // La puissance d'une source de brillance I non polarisee est:
     494 // P = I * S * Omega * dNu
     495 // La puissance recue pour une seule polarisation est:
     496 // P = 1/2 * I * S * Omega * dNu
     497 // La puissance recue pour un telescope (plein) est
     498 //    en remplacant Omega = lambda^2/S
     499 // P = 1/2 * I * lambda^2 * dNu
     500 // En appliquant la loi de Rayleigh: I = 2*k*Tb/lambda^2
     501 //   on obtient
     502 // P = 1/2 * 2*k*Tb * dNu = k * Tb * dNu
     503 // La temperature d'antenne est definie comme
     504 // P = k * Ta * dNu
     505 // Donc pour un ciel de temperature de brillance Tb
     506 //   et pour une antenne qui mesure une seule polarisation
     507 // on a: Ta = Tb
     508
     509 cout<<"\n>>>>\n>>>> Signal Analysis\n>>>>"<<endl;
    411510 cout<<"Facteur polar = "<<facpolar<<endl;
    412511
    … …  
    417516 planck.SetSpectraUnit(PlanckSpectra::ANGSFLUX); // radiance W/m^2/Sr/Hz
    418517
    419  // Signal
    420  double psig_2polar = fhi * surfeff;
    421  double tsig_2polar = psig_2polar / k_Boltzman_Cst / (dnuhiz*1e9);
    422  double ssig_2polar = psig_2polar / surfeff / (dnuhiz*1e9) / Jansky2Watt_cst;
    423  double psig = facpolar * psig_2polar;
    424  double tsig = facpolar * tsig_2polar;
    425  double ssig = facpolar * ssig_2polar;
    426  cout<<"\nSignal("<<mhiref<<" Msol):"<<endl
    427      <<"     P="<<psig<<" W"<<endl
    428      <<"     flux density = "<<ssig*1.e6<<" mu_Jy  (for Dnu="<<dnuhiz<<" GHz)"<<endl
    429      <<"                   ("<<ssig_2polar<<" mu_Jy for 2 polars)"<<endl
    430      <<"     Antenna temperature: tsig="<<tsig<<" K"<<endl;
     518 // ---
     519 // --- Signal HI
     520 // ---
     521 // Power emitted by HI
     522 cout<<"--- Power from HI for M = "<<mhiref<<" Msol"<<endl;
     523 cout<<"Flux factor = "<<hifactor<<" at redshift = "<<redshift[0]<<endl;
     524 cout<<"Luminosite = "<<hifactor*Msol2LumiHI(mhiref)<<" W"<<endl;
     525
     526 double fhi_ref[3], sfhi_ref[3];
     527 for(int i=0;i<3;i++) {
     528   fhi_ref[i] = hifactor*Msol2FluxHI(mhiref,dlum[i]);
     529   sfhi_ref[i] = fhi_ref[i] / (dnuhiz[i]*1e9) / Jansky2Watt_cst;
     530 }
     531 cout<<"FluxHI("<<dlum[0]<<" Mpc) all polar: "<<fhi_ref[0]<<" W/m^2 = "
     532     <<fhi_ref[0]/Jansky2Watt_cst<<" Jy.Hz"
     533     <<" in ["<<fhi_ref[1]/Jansky2Watt_cst<<","<<fhi_ref[2]/Jansky2Watt_cst<<"]"<<endl;
     534 cout<<"If spread over pixel depth ("<<dnuhiz[0]<<" GHz):"<<endl
     535     <<"  flux density = "<<sfhi_ref[0]*1.e6<<" mu_Jy"
     536     <<" in ["<<sfhi_ref[1]<<","<<sfhi_ref[2]<<"]"<<endl;
     537
     538 // Signal HI
     539 double psig_2polar[3], tasig_2polar[3], ssig_2polar[3], isig_2polar[3], tsig_2polar[3];
     540 double psig[3], tasig[3], ssig[3], isig[3], tsig[3];
     541 double doplarge[3], dzvrot[3];
     542 
     543 for(int i=0;i<3;i++) {
     544   psig_2polar[i] = fhi_ref[i] * surfeff;
     545   tasig_2polar[i] = psig_2polar[i] / k_Boltzman_Cst / (dnuhiz[i]*1e9);
     546   ssig_2polar[i] = psig_2polar[i] / surfeff / (dnuhiz[i]*1e9) / Jansky2Watt_cst;
     547   isig_2polar[i] = ssig_2polar[i] / angsol_pix[i];
     548   tsig_2polar[i] = isig_2polar[i]*Jansky2Watt_cst
     549                    / (2.*pow(nuhiz[i]*1e9/(SpeedOfLight_Cst*1000.),2.)*k_Boltzman_Cst);
     550   psig[i] = facpolar * psig_2polar[i];
     551   tasig[i] = facpolar * tasig_2polar[i];
     552   ssig[i] = facpolar * ssig_2polar[i];
     553   isig[i] = facpolar * isig_2polar[i];
     554   tsig[i] = facpolar * tsig_2polar[i];
     555   doplarge[i] = LargeurDoppler(vrot,nuhiz[i]);
     556   dzvrot[i] = DzFrV(vrot,redshift[i]);
     557 }
     558 cout<<"\n--- Signal HI("<<mhiref<<" Msol) for Dnu="<<dnuhiz[0]<<" GHz :"<<endl
     559     <<" Observation:"<<endl
     560     <<"    Power="<<psig[0]<<" W in ["<<psig[1]<<","<<psig[2]<<"]"<<endl
     561     <<"    Flux density = "<<ssig[0]*1.e6<<" mu_Jy in ["<<ssig[1]*1.e6<<","<<ssig[2]*1.e6<<"]"<<endl
     562     <<"    Intensity = "<<isig[0]<<" Jy/sr in ["<<isig[1]<<","<<isig[2]<<"]"<<endl
     563     <<"    Antenna temperature = "<<tasig[0]<<" K in ["<<tasig[1]<<","<<tasig[2]<<"]"<<endl
     564     <<"    Brightness temperature = "<<tsig[0]<<" K in ["<<tsig[1]<<","<<tsig[2]<<"]"<<endl
     565     <<" 2 polars:"<<endl
     566     <<"    Power="<<psig_2polar[0]<<" W in ["<<psig_2polar[1]<<","<<psig_2polar[2]<<"]"<<endl
     567     <<"    Flux density = "<<ssig_2polar[0]*1.e6<<" mu_Jy in ["<<ssig_2polar[1]*1.e6<<","<<ssig_2polar[2]*1.e6<<"]"<<endl
     568     <<"    Intensity = "<<isig_2polar[0]<<" Jy/sr in ["<<isig_2polar[1]<<","<<isig_2polar[2]<<"]"<<endl
     569     <<"    Antenna temperature = "<<tasig_2polar[0]<<" K in ["<<tasig_2polar[1]<<","<<tasig_2polar[2]<<"]"<<endl
     570     <<"    Brightness temperature = "<<tsig_2polar[0]<<" K in ["<<tsig_2polar[1]<<","<<tsig_2polar[2]<<"]"<<endl;
    431571
    432572 // Elargissement doppler de la raie a 21cm: dNu = vrot/C * Nu(21cm) / (1+z)
    433  double doplarge = vrot / SpeedOfLight_Cst * nuhiz;
    434  double dzvrot = vrot / SpeedOfLight_Cst * unplusz;
    435  cout<<"     Doppler width="<<doplarge*1.e3<<" MHz  for rotation width of "<<vrot<<" km/s"<<endl
    436      <<"                 dx= "<<dzvrot<<" a z="<<redshift<<endl;
    437  if(doplarge>dnuhiz)
    438    cout<<"Warning: doppler width "<<doplarge<<" GHz > "<<dnuhiz<<" GHz redshift bin width"<<endl;
    439 
    440  // Synchrotron  (T en -2.7  -> Flux en -0.7  dans l'approximation Rayleigh)
    441  double tsynch = Tsynch408;
    442         if(fabs(indnu)>1.e-50) tsynch *= pow(nuhiz/nuhaslam,indnu);
    443  planck.SetTemperature(tsynch);
    444  double psynch_2polar = planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1e9);
    445  double ssynch_2polar = psynch_2polar / surfeff / (dnuhiz*1e9) / Jansky2Watt_cst;
    446  double psynch = facpolar * psynch_2polar;
    447  double ssynch = facpolar * ssynch_2polar;
    448  cout<<"\nSynchrotron: T="<<Tsynch408<<" K ("<<nuhaslam<<" GHz), "
    449      <<tsynch<<" K ("<<nuhiz<<" GHz)"<<endl
    450      <<"             P="<<psynch<<" W for pixel"<<endl
    451      <<"             flux density = "<<ssynch<<" Jy for pixel solid angle"<<endl;
    452 
    453  // CMB
    454  double tcmb = T_CMB_Par;
     573 cout<<"    Doppler for rotation width of "<<vrot<<" km/s"<<endl
     574     <<"            width="<<doplarge[0]*1.e3<<" MHz  in ["<<doplarge[1]*1.e3<<","<<doplarge[2]*1.e3<<"]"<<endl
     575     <<"            dzvrot= "<<dzvrot[0]<<" in ["<<dzvrot[1]<<","<<dzvrot[2]<<"]"<<endl;
     576 if(doplarge[0]>dnuhiz[0])
     577   cout<<"Warning: doppler width "<<doplarge[0]<<" GHz > "<<dnuhiz[0]<<" GHz redshift bin width"<<endl;
     578
     579 // ---
     580 // --- Synchrotron  (T en -2.7  -> Flux en -0.7  dans l'approximation Rayleigh)
     581 // ---
     582 double tsynch[3];
     583 double psynch_2polar[3], tasynch_2polar[3], ssynch_2polar[3], isynch_2polar[3];
     584 double psynch[3], tasynch[3], ssynch[3], isynch[3];
     585
     586 for(int i=0;i<3;i++) {
     587   tsynch[i] = Tsynch408;
     588   if(fabs(indnu)>1.e-50) tsynch[i] *= pow(nuhiz[i]/nuhaslam,indnu);
     589   planck.SetTemperature(tsynch[i]);
     590   psynch_2polar[i] = planck(nuhiz[i]*1.e+9) * surfeff * angsol_pix[i] * (dnuhiz[i]*1e9);
     591   tasynch_2polar[i] = psynch_2polar[i] / k_Boltzman_Cst / (dnuhiz[i]*1e9);
     592   ssynch_2polar[i] = psynch_2polar[i] / surfeff / (dnuhiz[i]*1e9) / Jansky2Watt_cst;
     593   isynch_2polar[i] = ssynch_2polar[i] / angsol_pix[i];
     594   psynch[i] = facpolar * psynch_2polar[i];
     595   tasynch[i] = facpolar * tasynch_2polar[i];
     596   ssynch[i] = facpolar * ssynch_2polar[i];
     597   isynch[i] = ssynch[i] / angsol_pix[i];
     598 }
     599 cout<<"\n--- Synchrotron: T="<<Tsynch408<<" K ("<<nuhaslam<<" GHz), "<<endl
     600     <<"    tsynch="<<tsynch[0]<<" K ("<<nuhiz[0]<<" GHz) in ["<<tsynch[1]<<","<<tsynch[2]<<"]"<<endl
     601     <<" Observation:"<<endl
     602     <<"    Power="<<psynch[0]<<" W for pixel in ["<<psynch[1]<<","<<psynch[2]<<"]"<<endl
     603     <<"    Flux density = "<<ssynch[0]<<" Jy for pixel solid angle in ["<<ssynch[1]<<","<<ssynch[2]<<"]"<<endl
     604     <<"    Intensity = "<<isynch[0]<<" Jy/sr in ["<<isynch[1]<<","<<isynch[2]<<"]"<<endl
     605     <<"    Antenna temperature = "<<tasynch[0]<<" K in ["<<tasynch[1]<<","<<tasynch[2]<<"]"<<endl
     606     <<" 2 polars:"<<endl
     607     <<"    Power="<<psynch_2polar[0]<<" W for pixel in ["<<psynch_2polar[1]<<","<<psynch_2polar[2]<<"]"<<endl
     608     <<"    Flux density = "<<ssynch_2polar[0]<<" Jy for pixel solid angle in ["<<ssynch_2polar[1]<<","<<ssynch_2polar[2]<<"]"<<endl
     609     <<"    Intensity = "<<isynch_2polar[0]<<" Jy/sr in ["<<isynch_2polar[1]<<","<<isynch_2polar[2]<<"]"<<endl
     610     <<"    Antenna temperature = "<<tasynch_2polar[0]<<" K in ["<<tasynch_2polar[1]<<","<<tasynch_2polar[2]<<"]"<<endl;
     611
     612 // ---
     613 // --- CMB
     614 // ---
    455615 planck.SetTemperature(tcmb);
    456  double pcmb_2polar = planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1e9);
    457  double scmb_2polar = pcmb_2polar / surfeff / (dnuhiz*1.e+9) / Jansky2Watt_cst;
    458  double pcmb = facpolar * pcmb_2polar;
    459  double scmb = facpolar * scmb_2polar;
    460  cout<<"\nCMB: T="<<tcmb<<" K"<<endl
    461      <<"     P="<<pcmb<<" W for pixel"<<endl
    462      <<"     flux density = "<<scmb<<" Jy for pixel solid angle"<<endl;
    463 
    464  // AGN
     616 double pcmb_2polar[3], tacmb_2polar[3], scmb_2polar[3], icmb_2polar[3];
     617 double pcmb[3], tacmb[3], scmb[3], icmb[3];
     618
     619 for(int i=0;i<3;i++) {
     620   pcmb_2polar[i] = planck(nuhiz[i]*1.e+9) * surfeff * angsol_pix[i] * (dnuhiz[i]*1e9);
     621   tacmb_2polar[i] = pcmb_2polar[i] / k_Boltzman_Cst / (dnuhiz[i]*1e9);
     622   scmb_2polar[i] = pcmb_2polar[i] / surfeff / (dnuhiz[i]*1e9) / Jansky2Watt_cst;
     623   icmb_2polar[i] = scmb_2polar[i] / angsol_pix[i];
     624   pcmb[i] = facpolar * pcmb_2polar[i];
     625   tacmb[i] = facpolar * tacmb_2polar[i];
     626   scmb[i] = facpolar * scmb_2polar[i];
     627   icmb[i] = scmb[i] / angsol_pix[i];
     628 }
     629 cout<<"\n--- CMB: T="<<tcmb<<" K"<<endl
     630     <<" Observation:"<<endl
     631     <<"    Power="<<pcmb[0]<<" W for pixel in ["<<pcmb[1]<<","<<pcmb[2]<<"]"<<endl
     632     <<"    Flux density = "<<scmb[0]<<" Jy for pixel solid angle in ["<<scmb[1]<<","<<scmb[2]<<"]"<<endl
     633     <<"    Intensity = "<<icmb[0]<<" Jy/sr in ["<<icmb[1]<<","<<icmb[2]<<"]"<<endl
     634     <<"    Antenna temperature = "<<tacmb[0]<<" K in ["<<tacmb[1]<<","<<tacmb[2]<<"]"<<endl
     635     <<" 2 polars:"<<endl
     636     <<"    Power="<<pcmb_2polar[0]<<" W for pixel in ["<<pcmb_2polar[1]<<","<<pcmb_2polar[2]<<"]"<<endl
     637     <<"    Flux density = "<<scmb_2polar[0]<<" Jy for pixel solid angle in ["<<scmb_2polar[1]<<","<<scmb_2polar[2]<<"]"<<endl
     638     <<"    Intensity = "<<icmb_2polar[0]<<" Jy/sr in ["<<icmb_2polar[1]<<","<<icmb_2polar[2]<<"]"<<endl
     639     <<"    Antenna temperature = "<<tacmb_2polar[0]<<" K in ["<<tacmb_2polar[1]<<","<<tacmb_2polar[2]<<"]"<<endl;
     640
     641 // ---
     642 // --- AGN
     643 // ---
    465644 double flux_agn = pow(10.,lflux_agn);
    466  double mass_agn = FluxHI2Msol(flux_agn*Jansky2Watt_cst,dlum);
    467  cout<<"\nAGN: log10(S_agn)="<<lflux_agn<<" -> S_agn="
    468      <<flux_agn<<" Jy -> "<<mass_agn<<" equiv. Msol/Hz"<<endl;
    469  double flux_agn_pix = flux_agn*(dnuhiz*1e9);
    470  double mass_agn_pix = FluxHI2Msol(flux_agn_pix*Jansky2Watt_cst,dlum);
    471  double lmass_agn_pix = log10(mass_agn_pix);
    472  cout<<"...pixel: f="<<flux_agn_pix<<" 10^-26 W/m^2"
    473      <<" -> "<<mass_agn_pix<<" Msol -> log10 = "<<lmass_agn_pix<<endl;
     645 double mass_agn[3], flux_agn_pix[3], mass_agn_pix[3], lmass_agn_pix[3];
     646 for(int i=0;i<3;i++) {
     647   mass_agn[i] = FluxHI2Msol(flux_agn*Jansky2Watt_cst,dlum[i]);
     648   flux_agn_pix[i] = flux_agn*(dnuhiz[i]*1e9);
     649   mass_agn_pix[i] = FluxHI2Msol(flux_agn_pix[i]*Jansky2Watt_cst,dlum[i]);
     650   lmass_agn_pix[i] = log10(mass_agn_pix[i]);
     651 }
     652 cout<<"\n--- AGN: log10(S_agn)="<<lflux_agn<<"   S_agn="<<flux_agn<<" Jy :"<<endl
     653     <<"  mass_agn = "<<mass_agn[0]<<" equiv. Msol/Hz in ["<<mass_agn[1]<<","<<mass_agn[2]<<"]"<<endl
     654     <<"  flux_agn_pix = "<<flux_agn_pix[0]<<" 10^-26 W/m^2 in ["<<flux_agn_pix[1]<<","<<flux_agn_pix[2]<<"]"<<endl
     655     <<"  mass_agn_pix = "<<mass_agn_pix[0]<<" Msol in ["<<mass_agn_pix[1]<<","<<mass_agn_pix[2]<<"]"<<endl
     656     <<"  log10(mass_agn_pix) = "<<lmass_agn_pix[0]<<"  in ["<<lmass_agn_pix[1]<<","<<lmass_agn_pix[2]<<"]"<<endl;
    474657
    475658 // =====================================================================
    … …  
    541724 // =====================================================================
    542725
    543  cout<<"\n---\n--- Noise analysis \n---"<<endl;
    544  double psys = k_Boltzman_Cst * Tsys * (dnuhiz*1.e+9);
    545  cout<<"Noise: T="<<Tsys<<" K, P="<<psys<<" W  (for Dnu="<<dnuhiz<<" GHz)"<<endl;
     726 cout<<"\n>>>>\n>>>> Noise analysis\n>>>>"<<endl;
     727 double psys[3];
     728 for(int i=0;i<3;i++) psys[i] = k_Boltzman_Cst * Tsys * (dnuhiz[i]*1.e+9);
     729 cout<<"Noise: T="<<Tsys<<" K"<<endl
     730     <<"       P="<<psys[0]<<" W in ["<<psys[1]<<","<<psys[2]<<"]"<<endl;
    546731
    547732 cout<<"...Computation assume that noise dominate the signal."<<endl;
    … …  
    549734   cout<<"...Assuming 2 polarisations measurements with 2 different electronics."<<endl;
    550735
    551  double slim,slim_nl,SsN,SsN_nl,smass,smass_nl;
    552736
    553737 //---
    … …  
    566750   } else continue;
    567751
    568    slim = 2. * k_Boltzman_Cst * Tsys / surfeff
    569                / sqrt(fac*(dnuhiz*1.e+9)*tobs) /Jansky2Watt_cst;
    570    if(ya2polar) slim /= sqrt(2.);
    571    SsN = ssig_2polar / slim;
    572    smass = mhiref / ssig_2polar * slim;
     752   double slim[3], SsN[3], smass[3], pkbruit[3];
     753   for(int i=0;i<3;i++) {
     754     slim[i] = 2. * k_Boltzman_Cst * Tsys / surfeff
     755               / sqrt(fac*(dnuhiz[i]*1.e+9)*tobs) /Jansky2Watt_cst;
     756     if(ya2polar) slim[i] /= sqrt(2.);
     757     SsN[i] = ssig_2polar[i] / slim[i];  // independant de  angsol_pix*surfeff
     758     smass[i] = mhiref / ssig_2polar[i] * slim[i];
     759     pkbruit[i] = pow(smass[i]/mhiref,2.)*vol_pixel;
     760   }
    573761   cout<<"for 1 lobe:"<<endl
    574        <<"   Slim    = "<<slim*1.e6<<" mu_Jy"<<endl
    575        <<"   S/N     = "<<SsN<<endl
    576        <<"   Mass HI = "<<smass<<" Msol"<<endl;
    577 
    578    slim_nl = slim * sqrt(nlobes);
    579    SsN_nl = ssig_2polar / slim_nl;
    580    smass_nl = mhiref / ssig_2polar * slim_nl;
    581    cout<<"for "<<nlobes<<" lobes:"<<endl
    582        <<"   Flux    = "<<slim_nl*1.e6<<" mu_Jy"<<endl
    583        <<"   S/N     = "<<SsN_nl<<endl
    584        <<"   Mass HI = "<<smass_nl<<" Msol"<<endl;
     762       <<"   Slim    = "<<slim[0]*1.e6<<" mu_Jy in ["<<slim[1]*1.e6<<","<<slim[2]*1.e6<<"]"<<endl
     763       <<"   S/N     = "<<SsN[0]<<" in ["<<SsN[1]<<","<<SsN[2]<<"]"<<endl
     764       <<"   Mass HI = "<<smass[0]<<" Msol in ["<<smass[1]<<","<<smass[2]<<"]"<<endl
     765       <<"   Pk      = "<<pkbruit[0]<<" Mpc^3 in ["<<pkbruit[1]<<","<<pkbruit[2]<<"]"<<endl;
     766
     767   if(fabs(nlobes[0]-1.)>1e-6) {
     768     double slim_nl[3], SsN_nl[3], smass_nl[3], pkbruit_nl[3];
     769     for(int i=0;i<3;i++) {
     770       slim_nl[i] = slim[i] * sqrt(nlobes[i]);
     771       SsN_nl[i] = ssig_2polar[i] / slim_nl[i];
     772       smass_nl[i] = mhiref / ssig_2polar[i] * slim_nl[i];
     773       pkbruit_nl[i] = pow(smass_nl[i]/mhiref,2.)*vol_pixel;
     774     }
     775     cout<<"for "<<nlobes[0]<<" lobes[i] in ["<<nlobes[1]<<","<<nlobes[2]<<"] :"<<endl
     776       <<"   Slim    = "<<slim_nl[0]*1.e6<<" mu_Jy in ["<<slim_nl[1]*1.e6<<","<<slim_nl[2]*1.e6<<"]"<<endl
     777       <<"   S/N     = "<<SsN_nl[0]<<" in ["<<SsN_nl[1]<<","<<SsN_nl[2]<<"]"<<endl
     778       <<"   Mass HI = "<<smass_nl[0]<<" Msol in ["<<smass_nl[1]<<","<<smass_nl[2]<<"]"<<endl
     779       <<"   Pk      = "<<pkbruit_nl[0]<<" Mpc^3 in ["<<pkbruit_nl[1]<<","<<pkbruit_nl[2]<<"]"<<endl;
     780   }
    585781
    586782 }
    587783
    588784 return 0;
     785 }
     786
     787
     788//-------------------------------------------------------------------------------------------
     789double LargeurDoppler(double v, double nu)
     790// largeur doppler pour une vitesse v en km/s et une frequence nu
     791{
     792 return v / SpeedOfLight_Cst * nu;
    589793}
     794
     795double DzFrV(double v, double zred)
     796// largeur en redshift pour une vitesse v en km/s au redshift zred
     797{
     798  return v / SpeedOfLight_Cst * (1. + zred);
     799}
     800
     801double DNuFrDz(double dzred,double nu_at_0,double zred)
     802// Largeur DNu pour une largeur en redshift "dzred" au redshift "zred"
     803//    pour la frequence "nu_at_0" a z=0
     804// nu =  NuHi(z=0)/(1.+z0)
     805// dnu = NuHi(z=0)/(1.+z0-dz/2) - NuHi/(1.+z0+dz/2)
     806//     = NuHi(z=0)*dz/[ (1+z0)^2 - (dz/2)^2 ]
     807//     = NuHi(z=0)*dz/(1.+z0)^2 / [ 1 - [dz/(1+z0)/2)]^2 ]
     808//     = NuHi(z=0)*dz/(1.+z0)^2 / [1 - dz/(1+z0)/2] / [1 + dz/(1+z0)/2]
     809//    ~= NuHi(z=0)*dz/(1.+z0)^2   (approx. pour dz<<z0 a l'ordre (dz/z0)^2)
     810{
     811  double zp1 = 1.+zred;
     812  return nu_at_0*dzred/(zp1*zp1)/(1.-dzred/zp1/2.)/(1.+dzred/zp1/2.);
     813}
     814
     815double DzFrDNu(double dnu_at_0,double nu_at_0,double zred)
     816// Largeur en redshift au redshift "zred" pour une largeur
     817// en frequence "dnu_at_0" a la frequence "nu_at_0" a z=0
     818{
     819  if(dnu_at_0<=0.) return 0.;
     820  double zp1 = 1.+zred;
     821  double dnusnu0 = dnu_at_0/nu_at_0;
     822  return 2./dnusnu0 * (sqrt(1.+(dnusnu0*zp1)*(dnusnu0*zp1)) - 1.);
     823}
     824double DzFrDNuApprox(double dnu_at_0,double nu_at_0,double zred)
     825// idem DzFrDNu mais on utilise l'approximation: dnu=NuHi(z=0)*dz/(1.+z0)^2
     826{
     827  double zp1 = 1.+zred;
     828  return dnu_at_0/nu_at_0 *(zp1*zp1);
     829}
     830
     831double ZFrLos(double loscom,CosmoCalc& univ)
     832// Recherche du redshift correspondant a une distance comobile
     833// le long de la ligne de visee egale a "loscom" Mpc
     834// et pour un univers "univ"
     835{
     836  double dz = univ.ZMax()/10.; if(dz<=0.) dz = 0.1;
     837  double zmin=0., zmax=0.;
     838  while(univ.Dloscom(zmax)<loscom) zmax += dz;
     839  if(zmax==0.) return 0.;
     840  for(int i=0; i<6; i++) {
     841    zmin=zmax-dz; if(zmin<0.) zmin=0.;
     842    dz /= 10.;
     843    for(double z=zmin; z<zmax+dz; z+=dz) {
     844      double d = univ.Dloscom(z);
     845      if(d<loscom) continue;
     846      zmax = z;
     847      //cout<<"ZFrLos: z="<<zmax<<" d="<<d<<" / "<<loscom<<endl;
     848      break;
     849    }
     850  }
     851  return zmax;
     852}
     853
     854double AngsolEqTelescope(double nu /* GHz */,double telsurf /* m^2 */)
     855/*
     856Calcul de l'angle solide (sr) equivalent pour un telescope
     857de surface totale "telsurf" (m^2) a la frequence "nu" (GHz)
     858- Soit D(t) la figure de diffraction du telescope telle que D(t=0)=1
     859       (t = angle depuis l'axe optique)
     860  telescope circulaire de diametre D:
     861      D(t) = [2J1(Pi*D*t/l)/(Pi*D*t/l)]
     862  telescope rectangulaire axb:
     863      D(t) = [sin(Pi*a*t/l)/(Pi*a*t/l)]*[sin(Pi*b*t/l)/(Pi*b*t/l)]
     864- On cherche l'angle solide equivalent (par ex d'un cylindre de hauteur 1)
     865  Int[ D(t)^2 dOmega ] =  Int[ D(t)^2 2Pi t dt ] = Lambda^2/S
     866- En conclusion, pour un ciel d'intensite uniforme I, on a:
     867  P = I * S * Omega * dNu = I * S * (Lambda^2/S) * dNu
     868*/
     869{
     870  double lambda = SpeedOfLight_Cst*1000./(nu*1.e9);
     871  return lambda*lambda / telsurf;
     872}

  • trunk/Cosmo/SimLSS/genefluct3d.cc

    r3358 r3365  
    12081208
    12091209double GeneFluct3D::TurnFluct2MeanNumber(double val_by_mpc3)
     1210// ATTENTION: la gestion des pixels<0 proposee ici induit une perte de variance
     1211// sur la carte, le spectre Pk reconstruit sera plus faible!
     1212// L'effet sera d'autant plus grand que le nombre de pixels<0 sera grand.
    12101213{
    12111214 if(lp_>0) cout<<"--- TurnFluct2MeanNumber : "<<val_by_mpc3<<" quantity (gal or mass)/Mpc^3"<<endl;

  • trunk/Cosmo/SimLSS/geneutils.cc

    r3330 r3365  
    278278//               dtheta in [0,Pi]
    279279// Return: l'angle solide en steradian
    280 // Approx pour theta petit: PI theta^2
     280// Approx pour theta petit: PI * theta^2
    281281{
    282282  return   2.*M_PI * (1.-cos(dtheta));
    283283}
    284284
    285 //-------------------------------------------------------------------
    286 unsigned long PoissRandLimit(double mu,double mumax)
    287 {
    288   double pp,ppi,x;
    289   unsigned long n;
    290 
    291   if(mu>=mumax) {
    292     pp = sqrt(mu);
    293     while( (x=pp*NorRand()) < -mu );
    294     return (unsigned long)(mu+x+0.5);
    295   }
    296 
    297   ppi = pp = exp(-mu);
    298   x = drand01();
    299   n = 0;
    300   while (x > ppi) {
    301     n++;
    302     pp = mu*pp/(double)n;
    303     ppi += pp;
    304   }
    305   return n;
     285double FrAngSol(double angsol)
     286// Retourne la demi-ouverture "dtheta" d'une calotte spherique d'angle solide "angsol"
     287// Input: angle solide de la calotte spherique en steradians
     288// Return: demi-ouverture de la calotte spherique en radians
     289{
     290  angsol = 1. - angsol/(2.*M_PI);
     291  if(angsol<-1. || angsol>1.) return -1.;
     292  return acos(angsol);
    306293}
    307294

  • trunk/Cosmo/SimLSS/geneutils.h

    r3330 r3365  
    8686double AngSol(double dtheta,double dphi,double theta0=M_PI/2.);
    8787double AngSol(double dtheta);
    88 
    89 //----------------------------------------------------
    90 unsigned long PoissRandLimit(double mu,double mumax=10.);
     88double FrAngSol(double angsol);
    9189
    9290//----------------------------------------------------

  • trunk/Cosmo/SimLSS/schechter.cc

    r3345 r3365  
    429429//******************* Le Flux a 21cm ********************//
    430430///////////////////////////////////////////////////////////
     431double Msol2LumiHI(double m)
     432// Input:
     433//    m : masse de HI en "Msol"
     434// Return:
     435//    la luminosite en W
     436// L(m) =  2.393e+18* m(en masse solaire) en W
     437{
     438  return 2.393e+18 * m;
     439}
    431440
    432441double Msol2FluxHI(double m,double d)

  • trunk/Cosmo/SimLSS/schechter.h

    r3345 r3365  
    9393
    9494//-----------------------------------------------------------------------------------
     95double Msol2LumiHI(double m);
    9596double Msol2FluxHI(double m,double d);
    9697double FluxHI2Msol(double f,double d);
Note: See TracChangeset for help on using the changeset viewer.