| 1 | #include "sopnamsp.h"
|
|---|
| 2 | #include "machdefs.h"
|
|---|
| 3 | #include <iostream>
|
|---|
| 4 | #include <stdlib.h>
|
|---|
| 5 | #include <stdio.h>
|
|---|
| 6 | #include <string.h>
|
|---|
| 7 | #include <math.h>
|
|---|
| 8 | #include <unistd.h>
|
|---|
| 9 |
|
|---|
| 10 | #include <vector>
|
|---|
| 11 |
|
|---|
| 12 | #include "constcosmo.h"
|
|---|
| 13 | #include "cosmocalc.h"
|
|---|
| 14 | #include "geneutils.h"
|
|---|
| 15 | #include "schechter.h"
|
|---|
| 16 | #include "pkspectrum.h"
|
|---|
| 17 | #include "planckspectra.h"
|
|---|
| 18 |
|
|---|
| 19 | /* --- Check Peterson at al. astro-ph/0606104 v1 (pb facteur sqrt(2) sur S/N !)
|
|---|
| 20 | cmvdefsurv -U 0.75,0.3,0.7,-1,1 -V 300 -z 0.0025,0.2,Z -x 1,90,A -O 400000,6000 -N 75 -M 6.156e9 -F 3 -2 1.5
|
|---|
| 21 | --- */
|
|---|
| 22 |
|
|---|
| 23 | //-----------------------------------------------------------------------------------------------------------
|
|---|
| 24 | inline double rad2deg(double trad) {return trad/M_PI*180.;}
|
|---|
| 25 | inline double rad2min(double trad) {return trad/M_PI*180.*60.;}
|
|---|
| 26 | inline double rad2sec(double trad) {return trad/M_PI*180.*3600.;}
|
|---|
| 27 | inline double deg2rad(double tdeg) {return tdeg*M_PI/180.;}
|
|---|
| 28 | inline double min2rad(double tmin) {return tmin*M_PI/(180.*60.);}
|
|---|
| 29 | inline double sec2rad(double tsec) {return tsec*M_PI/(180.*3600.);}
|
|---|
| 30 |
|
|---|
| 31 | double LargeurDoppler(double v /* km/s */, double nu);
|
|---|
| 32 | double DzFrV(double v /* km/s */, double zred);
|
|---|
| 33 | double DNuFrDz(double dzred,double nu_at_0,double zred);
|
|---|
| 34 | double DzFrDNu(double dnu_at_0,double nu_at_0,double zred);
|
|---|
| 35 | double DzFrDNuApprox(double dnu_at_0,double nu_at_0,double zred);
|
|---|
| 36 | double ZFrLos(double loscom /* Mpc com */,CosmoCalc& univ);
|
|---|
| 37 | double AngsolEqTelescope(double nu /* GHz */,double telsurf /* m^2 */);
|
|---|
| 38 | void usage(void);
|
|---|
| 39 |
|
|---|
| 40 | void usage(void) {
|
|---|
| 41 | cout<<"cmvdefsurv [options] -x dx,txlarg[,unit_x] -y dy,tylarg[,unit_y] -z dz,tzlarg[,unit_z] redshift"<<endl
|
|---|
| 42 | <<"----------------"<<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
|
|---|
| 46 | <<"-- Unites pour X-Y:"<<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
|
|---|
| 51 | <<"----------------"<<endl
|
|---|
| 52 | <<" -K k,dk,pk : k(Mpc^-1) dk(Mpc^-1) pk(a z=0 en Mpc^-3) pour estimer la variance cosmique"<<endl
|
|---|
| 53 | <<"----------------"<<endl
|
|---|
| 54 | <<" -O surf,tobs,eta : surface effective (m^2), temps d\'observation (s), efficacite d\'ouverture"<<endl
|
|---|
| 55 | <<" -N Tsys : temperature du system (K)"<<endl
|
|---|
| 56 | <<" -L lobewidth,freqlob : taille du lobe d\'observation (FWHM) en arcmin (def= rien)"<<endl
|
|---|
| 57 | <<" pour la frequence freqlob en MHz"<<endl
|
|---|
| 58 | <<" Si freqlob<=0 : la frequence de reference est celle du redshift etudie"<<endl
|
|---|
| 59 | <<" Si freqlob absent : la frequence de reference 1.4 GHz"<<endl
|
|---|
| 60 | <<" lobewidth est pris comme la FWHM du lobe gaussien equivalent"<<endl
|
|---|
| 61 | <<" -2 : two polarisations measured"<<endl
|
|---|
| 62 | <<" -M : masse de HI de reference (MSol), si <=0 mean schechter in pixel"<<endl
|
|---|
| 63 | <<" -F : HI flux factor to be applied for our redshift"<<endl
|
|---|
| 64 | <<" -V Vrot : largeur en vitesse (km/s) pour l\'elargissement doppler (def=300km/s)"<<endl
|
|---|
| 65 | <<"----------------"<<endl
|
|---|
| 66 | <<" -S Tsynch,indnu : temperature (K) synch a 408 Mhz, index d\'evolution"<<endl
|
|---|
| 67 | <<" (indnu==0 no evolution with freq.)"<<endl
|
|---|
| 68 | <<"----------------"<<endl
|
|---|
| 69 | <<" -U h100,om0,ol0,w0,or0,flat : cosmology"<<endl
|
|---|
| 70 | <<"----------------"<<endl
|
|---|
| 71 | <<" -A <log10(S_agn)> : moyenne du flux AGN en Jy dans le pixel"<<endl
|
|---|
| 72 | <<" redshift : redshift moyen du survey"<<endl
|
|---|
| 73 | <<endl;
|
|---|
| 74 | }
|
|---|
| 75 |
|
|---|
| 76 | //-------------------------------------------------------------------------------------------
|
|---|
| 77 | int main(int narg,char *arg[])
|
|---|
| 78 | {
|
|---|
| 79 | // ---
|
|---|
| 80 | // --- Valeurs fixes ou par defaut
|
|---|
| 81 | // ---
|
|---|
| 82 |
|
|---|
| 83 | //-- WMAP
|
|---|
| 84 | unsigned short flat = 0;
|
|---|
| 85 | double h100=0.71, om0=0.267804, or0=7.9e-05*0., ol0=0.73,w0=-1.;
|
|---|
| 86 |
|
|---|
| 87 | //-- Schechter HIMASS
|
|---|
| 88 | double h75 = h100 / 0.75;
|
|---|
| 89 | double nstar = 0.006*pow(h75,3.); //
|
|---|
| 90 | double mstar = pow(10.,9.8); // MSol
|
|---|
| 91 | double alpha = -1.37;
|
|---|
| 92 | cout<<"nstar= "<<nstar<<" mstar="<<mstar<<" alpha="<<alpha<<endl;
|
|---|
| 93 |
|
|---|
| 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.;
|
|---|
| 107 | double tobs = 6000., surftot = 400000., eta_eff = 1.;
|
|---|
| 108 | // taille du lobe d'observation en arcmin pour la frequence
|
|---|
| 109 | double lobewidth0 = -1., lobefreq0 = Fr_HyperFin_Par*1.e3;
|
|---|
| 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.;
|
|---|
| 119 | double mhiref = -1.; // reference Mass en HI (def integ schechter)
|
|---|
| 120 |
|
|---|
| 121 | // ---
|
|---|
| 122 | // --- Decodage arguments
|
|---|
| 123 | // ---
|
|---|
| 124 | char c;
|
|---|
| 125 | while((c = getopt(narg,arg,"h2x:y:z:N:S:O:M:F:V:U:L:A:K:")) != -1) {
|
|---|
| 126 | switch (c) {
|
|---|
| 127 | case 'x' :
|
|---|
| 128 | sscanf(optarg,"%lf,%lf,%c",&dx,&txlarg,&unit_x);
|
|---|
| 129 | break;
|
|---|
| 130 | case 'y' :
|
|---|
| 131 | sscanf(optarg,"%lf,%lf,%c",&dy,&tylarg,&unit_y);
|
|---|
| 132 | break;
|
|---|
| 133 | case 'z' :
|
|---|
| 134 | sscanf(optarg,"%lf,%lf,%c",&dz,&tzlarg,&unit_z);
|
|---|
| 135 | break;
|
|---|
| 136 | case 'O' :
|
|---|
| 137 | sscanf(optarg,"%lf,%lf,%lf",&surftot,&tobs,&eta_eff);
|
|---|
| 138 | break;
|
|---|
| 139 | case 'L' :
|
|---|
| 140 | sscanf(optarg,"%lf,%lf",&lobewidth0,&lobefreq0);
|
|---|
| 141 | break;
|
|---|
| 142 | case 'N' :
|
|---|
| 143 | sscanf(optarg,"%lf",&Tsys);
|
|---|
| 144 | break;
|
|---|
| 145 | case 'S' :
|
|---|
| 146 | sscanf(optarg,"%lf,%lf",&Tsynch408,&indnu);
|
|---|
| 147 | break;
|
|---|
| 148 | case 'M' :
|
|---|
| 149 | sscanf(optarg,"%lf",&mhiref);
|
|---|
| 150 | break;
|
|---|
| 151 | case 'F' :
|
|---|
| 152 | sscanf(optarg,"%lf",&hifactor);
|
|---|
| 153 | break;
|
|---|
| 154 | case 'V' :
|
|---|
| 155 | sscanf(optarg,"%lf",&vrot);
|
|---|
| 156 | break;
|
|---|
| 157 | case 'U' :
|
|---|
| 158 | sscanf(optarg,"%lf,%lf,%lf,%lf,%hu",&h100,&om0,&ol0,&w0,&flat);
|
|---|
| 159 | break;
|
|---|
| 160 | case '2' :
|
|---|
| 161 | ya2polar = true;
|
|---|
| 162 | facpolar = 1.0;
|
|---|
| 163 | break;
|
|---|
| 164 | case 'A' :
|
|---|
| 165 | sscanf(optarg,"%lf",&lflux_agn);
|
|---|
| 166 | break;
|
|---|
| 167 | case 'K' :
|
|---|
| 168 | sscanf(optarg,"%lf,%lf,%lf",&kcosm,&dkcosm,&pkcosm);
|
|---|
| 169 | break;
|
|---|
| 170 | case 'h' :
|
|---|
| 171 | default :
|
|---|
| 172 | usage(); return -1;
|
|---|
| 173 | }
|
|---|
| 174 | }
|
|---|
| 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 | // ---
|
|---|
| 180 | // --- Initialisation de la Cosmologie
|
|---|
| 181 | // ---
|
|---|
| 182 | CosmoCalc univ(flat,true,2.*redshift0);
|
|---|
| 183 | double perc=0.01,dzinc=redshift0/100.,dzmax=dzinc*10.; unsigned short glorder=4;
|
|---|
| 184 | univ.SetInteg(perc,dzinc,dzmax,glorder);
|
|---|
| 185 | 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;
|
|---|
| 274 | univ.Print(0.);
|
|---|
| 275 | double rhoc0 = univ.Rhoc(0.)*GCm3toMsolMpc3_Cst;
|
|---|
| 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;
|
|---|
| 348 | double Npix = (double)nx*(double)ny*(double)nz;
|
|---|
| 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 | // ---
|
|---|
| 371 | cout<<"\n>>>>\n>>>> Geometrie dans l'espace de Fourier\n>>>>"<<endl;
|
|---|
| 372 | cout<<"Array size: nx = "<<nx<<", ny = "<<ny<<", nz = "<<nz<<endl;
|
|---|
| 373 | double dk_x = 2.*M_PI/(nx*dx), knyq_x = M_PI/dx;
|
|---|
| 374 | double dk_y = 2.*M_PI/(nx*dy), knyq_y = M_PI/dy;
|
|---|
| 375 | double dk_z = 2.*M_PI/(nz*dz), knyq_z = M_PI/dz;
|
|---|
| 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 | // ---
|
|---|
| 386 | // --- Variance cosmique
|
|---|
| 387 | // ---
|
|---|
| 388 | // cosmique poisson
|
|---|
| 389 | // (sigma/P)^2 = 2*(2Pi)^3 / (4Pi k^2 dk Vsurvey) * [(1+n*P)/(n*P)]^2
|
|---|
| 390 | // nombre de mode = 1/2 * Vsurvey/(2Pi)^3 * 4Pi*k^2*dk
|
|---|
| 391 | if(kcosm>0. && pkcosm>0.) {
|
|---|
| 392 | double pk = pkcosm*pow(growthfac[0],2.);
|
|---|
| 393 | cout<<"\n>>>>\n>>>> variance cosmique pour k="<<kcosm<<" Mpc^-1, pk(z=0)="
|
|---|
| 394 | <<pkcosm<<", pk(z="<<redshift[0]<<")="<<pk<<"\n>>>>"<<endl;
|
|---|
| 395 | for(int i=0;i<3;i++) { // la correction de variance du au bruit de poisson
|
|---|
| 396 | double v = 1.1; if(i==1) v=1.5; if(i==2) v=2.0;
|
|---|
| 397 | double ngal = 1./(v-1.)/pk;
|
|---|
| 398 | cout<<" pour "<<ngal<<" gal/Mpc^3 on multiplie par "<<v<<" sigma/P"<<endl;
|
|---|
| 399 | }
|
|---|
| 400 |
|
|---|
| 401 | cout<<endl;
|
|---|
| 402 | vector<double> dk; if(dkcosm>0.) dk.push_back(dkcosm);
|
|---|
| 403 | dk.push_back(dk_x); dk.push_back(dk_y); dk.push_back(dk_z);
|
|---|
| 404 | for(int i=0;i<(int)dk.size();i++) { // la variance cosmique pure
|
|---|
| 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];
|
|---|
| 407 | cout<<" pour dk = "<<dk[i]<<" Mpc^-1: sigma/P = "<<vcosm
|
|---|
| 408 | <<" , Nmode = "<<nmode<<endl;
|
|---|
| 409 | }
|
|---|
| 410 | }
|
|---|
| 411 |
|
|---|
| 412 | // ---
|
|---|
| 413 | // --- Masse de HI
|
|---|
| 414 | // ---
|
|---|
| 415 | cout<<"\n>>>>\n>>>> Mass HI\n>>>>"<<endl;
|
|---|
| 416 | Schechter sch(nstar,mstar,alpha);
|
|---|
| 417 | sch.SetOutValue(1);
|
|---|
| 418 | cout<<"nstar= "<<nstar<<" mstar="<<mstar<<" alpha="<<alpha<<endl;
|
|---|
| 419 | cout<<"mstar*sch(mstar) = "<<sch(mstar)<<" Msol/Mpc^3/Msol"<<endl;
|
|---|
| 420 | int npt = 10000;
|
|---|
| 421 | double lnx1=log10(1.e+6), lnx2=log10(1.e+13), dlnx=(lnx2-lnx1)/npt;
|
|---|
| 422 | double masshimpc3 = IntegrateFuncLog(sch,lnx1,lnx2,0.001,dlnx,10.*dlnx,6);
|
|---|
| 423 | cout<<"Mass density: "<<masshimpc3<<" Msol/Mpc^3"<<endl;
|
|---|
| 424 |
|
|---|
| 425 | double masshipix = masshimpc3*vol_pixel;
|
|---|
| 426 | double masshitot = masshimpc3*vol_survey;
|
|---|
| 427 | cout<<"Pixel mass = "<<masshipix<<" Msol"<<endl
|
|---|
| 428 | <<"Total mass in survey = "<<masshitot<<" Msol"<<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;
|
|---|
| 432 | if(mhiref<=0.) mhiref = masshipix;
|
|---|
| 433 |
|
|---|
| 434 | sch.SetOutValue(0);
|
|---|
| 435 | cout<<"\nsch(mstar) = "<<sch(mstar)<<" /Mpc^3/Msol"<<endl;
|
|---|
| 436 | cout<<"Galaxy number density:"<<endl;
|
|---|
| 437 | for(double x=lnx1; x<lnx2-0.5; x+=1.) {
|
|---|
| 438 | double n = IntegrateFuncLog(sch,x,lnx2,0.001,dlnx,10.*dlnx,6);
|
|---|
| 439 | cout<<" m>"<<x<<" Msol: "<<n<<" /Mpc^3, "<<n*vol_pixel<<" /pixel, "
|
|---|
| 440 | <<n*vol_survey<<" in survey"<<endl;
|
|---|
| 441 | }
|
|---|
| 442 | sch.SetOutValue(1);
|
|---|
| 443 |
|
|---|
| 444 | // ---
|
|---|
| 445 | // --- Survey values
|
|---|
| 446 | // ---
|
|---|
| 447 | cout<<"\n>>>>\n>>>> Observations\n>>>>"<<endl;
|
|---|
| 448 | double surfeff = surftot * eta_eff;
|
|---|
| 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.};
|
|---|
| 470 | if(lobewidth0>0.) {
|
|---|
| 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 | // ---
|
|---|
| 490 | // --- Signal analysis
|
|---|
| 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;
|
|---|
| 510 | cout<<"Facteur polar = "<<facpolar<<endl;
|
|---|
| 511 |
|
|---|
| 512 | PlanckSpectra planck(T_CMB_Par);
|
|---|
| 513 | planck.SetSpectraApprox(PlanckSpectra::RAYLEIGH); // Rayleigh spectra
|
|---|
| 514 | planck.SetSpectraVar(PlanckSpectra::NU); // frequency
|
|---|
| 515 | planck.SetSpectraPower(PlanckSpectra::POWER); // output en W/....
|
|---|
| 516 | planck.SetSpectraUnit(PlanckSpectra::ANGSFLUX); // radiance W/m^2/Sr/Hz
|
|---|
| 517 |
|
|---|
| 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;
|
|---|
| 571 |
|
|---|
| 572 | // Elargissement doppler de la raie a 21cm: dNu = vrot/C * Nu(21cm) / (1+z)
|
|---|
| 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 | // ---
|
|---|
| 615 | planck.SetTemperature(tcmb);
|
|---|
| 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 | // ---
|
|---|
| 644 | double flux_agn = pow(10.,lflux_agn);
|
|---|
| 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;
|
|---|
| 657 |
|
|---|
| 658 | // =====================================================================
|
|---|
| 659 | // ---
|
|---|
| 660 | // --- Noise analysis
|
|---|
| 661 | // ---
|
|---|
| 662 | // --- Puissance du bruit pour un telescope de surface Ae et de BW dNu
|
|---|
| 663 | // Par definition la puissance du bruit est:
|
|---|
| 664 | // Pb = k * Tsys * dNu (W)
|
|---|
| 665 | // Pour une source (non-polarisee) de densite de flux (totale 2 polars)
|
|---|
| 666 | // St (exprimee en Jy = 10^-26 W/m^2/Hz)
|
|---|
| 667 | // Pt = St * Ae * dNu (puissance totale emise en W pour 2 polars)
|
|---|
| 668 | // P1 = 1/2 * St * Ae * dNu (puissance emise en W pour une polar)
|
|---|
| 669 | // la SEFD (system equivalent flux density en Jy) est definie comme
|
|---|
| 670 | // la densite de flux total (2 polars) "St" d'une source (non-polarisee)
|
|---|
| 671 | // dont la puissance P1 mesuree pour une seule polarisation
|
|---|
| 672 | // serait egale a la puissance du bruit. De P1 = Pb on deduit:
|
|---|
| 673 | // SEFD = 2 * k * Tsys / Ae (en Jy)
|
|---|
| 674 | // la puissance du bruit est: Pb = 1/2 * SEFD * Ae * dNu (en W)
|
|---|
| 675 | // la sensibilite Slim tient compte du temps d'integration et de la BW:
|
|---|
| 676 | // le nombre de mesures independantes est "2*dNu*Tobs" donc
|
|---|
| 677 | // Slim = SEFD / sqrt(2*dNu*Tobs) = 2*k*Tsys/[Ae*sqrt(2*dNu*Tobs) (en Jy)
|
|---|
| 678 | // --- Puissance du bruit pour un interferometre
|
|---|
| 679 | // Ae = surface d'un telescope elementaire
|
|---|
| 680 | // N = nombre de telescopes dans l'interferometre (Atot = N*Ae)
|
|---|
| 681 | // La sensibilite Slim en Jy est:
|
|---|
| 682 | // Slim = 2 * k * Tsys / [ Ae * Sqrt(2*N(N-1)/2 *dnu*Tobs) ]
|
|---|
| 683 | // = 2 * k * Tsys / [ Atot/N * Sqrt(2*N(N-1)/2*dnu*Tobs) ]
|
|---|
| 684 | // = 2 * k * Tsys / [ Atot * Sqrt((N-1)/N *dnu*Tobs) ]
|
|---|
| 685 | // - Interferometre a deux antennes:
|
|---|
| 686 | // Slim = 2 * k * Tsys / [ Atot * Sqrt(1/2 *dnu*Tobs) ]
|
|---|
| 687 | // - Interferometre a N antennes (N grand):
|
|---|
| 688 | // Slim -> 2 * k * Tsys / [ Atot * Sqrt(dnu*Tobs) ]
|
|---|
| 689 | // C'est aussi la formule pour un telescope unique de surface Atot
|
|---|
| 690 | // --- On ne mesure qu'une seule polarisation
|
|---|
| 691 | // Ces formules sont valables si on mesure 1 polarisation:
|
|---|
| 692 | // Slim est la densite de flux total "St" (2 polars) d'une source (non-polarisee)
|
|---|
| 693 | // qui donne la meme puissance que le bruit dans un detecteur qui ne
|
|---|
| 694 | // mesure qu'une seule polarisation:
|
|---|
| 695 | // Le rapport S/N pour une source de densite de flux St (totale 2 polars):
|
|---|
| 696 | // S/N = St / Slim
|
|---|
| 697 | // La puissance de bruit est, par definition:
|
|---|
| 698 | // Pb = 1/2 *Slim*Atot*dNu
|
|---|
| 699 | // = k*Tsys*sqrt(2*dNu/Tobs) pour N=2
|
|---|
| 700 | // = k*Tsys*sqrt(dNu/Tobs) pour N>>grand
|
|---|
| 701 | // La densite de flux d'une source a S/N=1 est:
|
|---|
| 702 | // St = Slim
|
|---|
| 703 | // La puissance d'une source a S/N=1 mesuree par un detecteur
|
|---|
| 704 | // qui ne mesure qu'une polar est:
|
|---|
| 705 | // P1_lim = 1/2 *Slim*Atot*dNu
|
|---|
| 706 | // --- On mesure les 2 polarisations avec deux voies d'electronique distinctes
|
|---|
| 707 | // la puissance du signal mesure est multipliee par 2
|
|---|
| 708 | // la puissance du bruit est multipliee par sqrt(2)
|
|---|
| 709 | // on a donc un gain d'un facteur sqrt(2) sur le rapport S/N
|
|---|
| 710 | // (cela revient d'ailleur a doubler le temps de pose: Tobs -> 2*Tobs)
|
|---|
| 711 | // En notant arbitrairement: Slim' = Slim / sqrt(2)
|
|---|
| 712 | // ou Slim est defini par les formules ci-dessus
|
|---|
| 713 | // Le rapport S/N pour une source de densite de flux St (totale 2 polars):
|
|---|
| 714 | // (S/N)_2 = (S/N)_1 * sqrt(2) = (St / Slim) * sqrt(2) = St / Slim'
|
|---|
| 715 | // La densite de flux d'une source a S/N=1 est:
|
|---|
| 716 | // St = Slim' = Slim / sqrt(2)
|
|---|
| 717 | // La puissance d'une source a S/N=1 cumulee par les 2 detecteurs est:
|
|---|
| 718 | // P_lim = St*Atot*dNu = Slim'*Atot*dNu = 1/sqrt(2) *Slim*Atot*dNu
|
|---|
| 719 | // = P1_lim * sqrt(2)
|
|---|
| 720 | // La puissance de bruit cumulee par les 2 detecteurs est, par definition:
|
|---|
| 721 | // Pb = P_lim = Slim'*Atot*dNu = P1_lim * sqrt(2)
|
|---|
| 722 | // = 2*k*Tsys*sqrt(dNu/Tobs) pour N=2
|
|---|
| 723 | // = k*Tsys*sqrt(2*dNu/Tobs) pour N>>grand
|
|---|
| 724 | // =====================================================================
|
|---|
| 725 |
|
|---|
| 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;
|
|---|
| 731 |
|
|---|
| 732 | cout<<"...Computation assume that noise dominate the signal."<<endl;
|
|---|
| 733 | if(ya2polar)
|
|---|
| 734 | cout<<"...Assuming 2 polarisations measurements with 2 different electronics."<<endl;
|
|---|
| 735 |
|
|---|
| 736 |
|
|---|
| 737 | //---
|
|---|
| 738 | for(unsigned short it=0;it<2;it++) {
|
|---|
| 739 |
|
|---|
| 740 | double fac = 1.;
|
|---|
| 741 | if(it==0) { // Interferometre a 2 telescopes
|
|---|
| 742 | fac = 0.5;
|
|---|
| 743 | cout<<"\n...Observation limits for a 2 telescope interferometer (with complex correlator)"<<endl
|
|---|
| 744 | <<" (sensitivity is given for real or complex correlator output)"<<endl;
|
|---|
| 745 | } else if (it==1) { // Interferometre a N>> telescopes
|
|---|
| 746 | fac = 1.;
|
|---|
| 747 | cout<<"\n...Observation limits for a N (large) telescope interferometer (with complex correlator)"<<endl
|
|---|
| 748 | <<" (weak source limit sensitivity in a synthetised image)"<<endl
|
|---|
| 749 | <<" Also valid for a single dish telescope."<<endl;
|
|---|
| 750 | } else continue;
|
|---|
| 751 |
|
|---|
| 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 | }
|
|---|
| 761 | cout<<"for 1 lobe:"<<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 | }
|
|---|
| 781 |
|
|---|
| 782 | }
|
|---|
| 783 |
|
|---|
| 784 | return 0;
|
|---|
| 785 | }
|
|---|
| 786 |
|
|---|
| 787 |
|
|---|
| 788 | //-------------------------------------------------------------------------------------------
|
|---|
| 789 | double 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;
|
|---|
| 793 | }
|
|---|
| 794 |
|
|---|
| 795 | double 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 |
|
|---|
| 801 | double 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 |
|
|---|
| 815 | double 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 | }
|
|---|
| 824 | double 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 |
|
|---|
| 831 | double 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 |
|
|---|
| 854 | double AngsolEqTelescope(double nu /* GHz */,double telsurf /* m^2 */)
|
|---|
| 855 | /*
|
|---|
| 856 | Calcul de l'angle solide (sr) equivalent pour un telescope
|
|---|
| 857 | de 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 | }
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
