| 1 | /* ------------------------ Projet BAORadio --------------------
|
|---|
| 2 | Programme de calcul du spectre de puissance de bruit pour
|
|---|
| 3 | un interferometre (spectre moyenne sur 3D -> P_noise(k) )
|
|---|
| 4 |
|
|---|
| 5 | R. Ansari , C. Magneville - Juin 2010
|
|---|
| 6 |
|
|---|
| 7 | Usage: pknoise pknoise [-parname value] Diameter/Four2DRespTableFile OutPPFName
|
|---|
| 8 | -parname : -noise , -renmax, -scut, -ngen, -z, -bsize, -cszmpc, -nbin -prt
|
|---|
| 9 | --------------------------------------------------------------- */
|
|---|
| 10 |
|
|---|
| 11 | #include "machdefs.h"
|
|---|
| 12 | #include "sopnamsp.h"
|
|---|
| 13 | #include <iostream>
|
|---|
| 14 | #include <string>
|
|---|
| 15 | #include <math.h>
|
|---|
| 16 |
|
|---|
| 17 | #include <typeinfo>
|
|---|
| 18 |
|
|---|
| 19 | #include "mdish.h"
|
|---|
| 20 | #include "specpk.h"
|
|---|
| 21 | #include "interfconfigs.h"
|
|---|
| 22 | #include "radutil.h"
|
|---|
| 23 |
|
|---|
| 24 | #include "histinit.h"
|
|---|
| 25 | // #include "fiosinit.h"
|
|---|
| 26 | // #include "fitsioserver.h"
|
|---|
| 27 |
|
|---|
| 28 | #include "randr48.h"
|
|---|
| 29 |
|
|---|
| 30 | #include "timing.h"
|
|---|
| 31 | #include "ctimer.h"
|
|---|
| 32 | #include "slininterp.h"
|
|---|
| 33 |
|
|---|
| 34 | typedef DR48RandGen RandomGenerator ;
|
|---|
| 35 |
|
|---|
| 36 | // ---------------------------------------------------------------------
|
|---|
| 37 | // Test main
|
|---|
| 38 | // R. Ansari - Avril-Dec 2010
|
|---|
| 39 | // ---------------------------------------------------------------------
|
|---|
| 40 |
|
|---|
| 41 | void Usage()
|
|---|
| 42 | {
|
|---|
| 43 | cout << " Usage: pknoise [-parname value] Diameter/Four2DRespTableFile OutPPFName \n"
|
|---|
| 44 | << " -noise NoiseLevel (default=1.) \n"
|
|---|
| 45 | << " -renmax MaxValue (default : Do NOT renormalize 2D response value \n"
|
|---|
| 46 | << " -ngen NGen (default=0) number of noise fourier amp generations \n"
|
|---|
| 47 | << " NGen=0 -> Call ComputeNoisePk(), else generate Fourier Amplitudes (random) \n"
|
|---|
| 48 | << " NGen=-1 -> Call ComputeNoisePk(), dont change Four3DPk cell size \n"
|
|---|
| 49 | << " -bsize sx,sy,sz : 3D real space box size (default=512x512x256) \n"
|
|---|
| 50 | << " -cszmpc sx,sy,sz : 3D real space box cell size in Mpc (default=3x3x1.5) \n"
|
|---|
| 51 | << " -z redshift : redshift (default=1.0) \n"
|
|---|
| 52 | << " -pnf : Apply frequency dependent noise factor=ComovDa^2*(1+z)^2/Hz (default= NON) \n"
|
|---|
| 53 | << " -scut SCutValue (default= -100.) \n"
|
|---|
| 54 | << " if SCutValue<0. ==> SCut=MinNoisePower*(-SCutValue) \n"
|
|---|
| 55 | << " -prt PrtLev,PrtModulo (default=0,10) \n"
|
|---|
| 56 | << " -nbin PkNBin : nomber of histogram bins for Pnoise(k) (default= 256) \n"
|
|---|
| 57 | << " NOTE : this program needs redshift_da.txt andredshift_Hz.txt files" << endl;
|
|---|
| 58 | return;
|
|---|
| 59 | }
|
|---|
| 60 |
|
|---|
| 61 | //-------------------------------------------------------------------------
|
|---|
| 62 | // ------------------ MAIN PROGRAM ------------------------------
|
|---|
| 63 | //-------------------------------------------------------------------------
|
|---|
| 64 | int main(int narg, const char* arg[])
|
|---|
| 65 | {
|
|---|
| 66 | if ( (narg<3)||((narg>1)&&(strcmp(arg[1],"-h")==0)) ) {
|
|---|
| 67 | Usage();
|
|---|
| 68 | return 1;
|
|---|
| 69 | }
|
|---|
| 70 | cout << " ==== pknoise.cc : interferometer noise power spectrum computation ==== " << endl;
|
|---|
| 71 | // make sure SOPHYA modules are initialized
|
|---|
| 72 | SophyaInit();
|
|---|
| 73 | // FitsIOServerInit();
|
|---|
| 74 | InitTim();
|
|---|
| 75 | //--- decoding command line arguments
|
|---|
| 76 | string tits="pknoise";
|
|---|
| 77 | char cbuff[64];
|
|---|
| 78 | bool fgresptbl=false;
|
|---|
| 79 | double DIAMETRE=100.;
|
|---|
| 80 | string resptblname;
|
|---|
| 81 | double NoiseLevel=1.;
|
|---|
| 82 |
|
|---|
| 83 | bool fgrenorm=false;
|
|---|
| 84 | double rmax=1.;
|
|---|
| 85 | int NMAX = 0;
|
|---|
| 86 | double SCut=-100.;
|
|---|
| 87 | bool fgautoscut=true;
|
|---|
| 88 | double FacSCut=-SCut;
|
|---|
| 89 | double z_Redshift=1.0 ; // z=0.7 : 21 cm at z=0.7 -> 0.357 m
|
|---|
| 90 | double comov_dA_z=3330.; // Comoving radial distance = (1+z)dA
|
|---|
| 91 | double H_z=120.5; // Hubble_param(z)
|
|---|
| 92 | int box3dsz[3]={512,512,256};
|
|---|
| 93 | double cellsz[3]={3.,3.,3.};
|
|---|
| 94 | bool fgnoisefreq=false;
|
|---|
| 95 | int NBINPK=256;
|
|---|
| 96 | int prtlev=0;
|
|---|
| 97 | int prtmod=10;
|
|---|
| 98 |
|
|---|
| 99 | int ka=1;
|
|---|
| 100 | while (ka<(narg-1)) {
|
|---|
| 101 | if (strcmp(arg[ka],"-noise")==0) {
|
|---|
| 102 | NoiseLevel=atof(arg[ka+1]);
|
|---|
| 103 | ka+=2;
|
|---|
| 104 | }
|
|---|
| 105 | else if (strcmp(arg[ka],"-renmax")==0) {
|
|---|
| 106 | rmax=atof(arg[ka+1]); fgrenorm=true; ka+=2;
|
|---|
| 107 | }
|
|---|
| 108 | else if (strcmp(arg[ka],"-scut")==0) {
|
|---|
| 109 | SCut=atof(arg[ka+1]); ka+=2;
|
|---|
| 110 | if (SCut<0.) { FacSCut=-SCut; fgautoscut=true; }
|
|---|
| 111 | }
|
|---|
| 112 | else if (strcmp(arg[ka],"-ngen")==0) {
|
|---|
| 113 | NMAX=atoi(arg[ka+1]); ka+=2;
|
|---|
| 114 | }
|
|---|
| 115 | else if (strcmp(arg[ka],"-z")==0) {
|
|---|
| 116 | z_Redshift=atof(arg[ka+1]); ka+=2;
|
|---|
| 117 | }
|
|---|
| 118 | else if (strcmp(arg[ka],"-bsize")==0) {
|
|---|
| 119 | sscanf(arg[ka+1],"%d,%d,%d",box3dsz,box3dsz+1,box3dsz+2); ka+=2;
|
|---|
| 120 | }
|
|---|
| 121 | else if (strcmp(arg[ka],"-cszmpc")==0) {
|
|---|
| 122 | sscanf(arg[ka+1],"%lg,%lg,%lg",cellsz,cellsz+1,cellsz+2); ka+=2;
|
|---|
| 123 | }
|
|---|
| 124 | else if (strcmp(arg[ka],"-nbin")==0) {
|
|---|
| 125 | NBINPK=atoi(arg[ka+1]); ka+=2;
|
|---|
| 126 | }
|
|---|
| 127 | else if (strcmp(arg[ka],"-pnf")==0) {
|
|---|
| 128 | fgnoisefreq=true; ka++;
|
|---|
| 129 | }
|
|---|
| 130 | else if (strcmp(arg[ka],"-prt")==0) {
|
|---|
| 131 | sscanf(arg[ka+1],"%d,%d",&prtlev,&prtmod); ka+=2;
|
|---|
| 132 | }
|
|---|
| 133 | else break;
|
|---|
| 134 | }
|
|---|
| 135 |
|
|---|
| 136 | if ((ka+1)>=narg) {
|
|---|
| 137 | cout << " pknoise / Argument error " << endl;
|
|---|
| 138 | Usage();
|
|---|
| 139 | return 2;
|
|---|
| 140 | }
|
|---|
| 141 | if (isdigit(*arg[ka])) {
|
|---|
| 142 | fgresptbl=false;
|
|---|
| 143 | DIAMETRE=atof(arg[ka]);
|
|---|
| 144 | sprintf(cbuff,"pknoise_Dish(%g m)", DIAMETRE);
|
|---|
| 145 | }
|
|---|
| 146 | else {
|
|---|
| 147 | resptblname=arg[ka];
|
|---|
| 148 | fgresptbl=true;
|
|---|
| 149 | sprintf(cbuff,"pknoise_RespTblName=%s", arg[ka]);
|
|---|
| 150 | }
|
|---|
| 151 | tits=cbuff;
|
|---|
| 152 | string outfile = arg[ka+1];
|
|---|
| 153 | if (outfile==".") outfile = "pknoise.ppf";
|
|---|
| 154 | //-- end command line arguments
|
|---|
| 155 |
|
|---|
| 156 | int rc = 1;
|
|---|
| 157 | try { // exception handling try bloc at top level
|
|---|
| 158 | cout << " pknoise[0] : Executing, output file= " << outfile << endl;
|
|---|
| 159 | POutPersist po(outfile);
|
|---|
| 160 |
|
|---|
| 161 | cout << " pknoise[0.b] : Creating interpolating objects from redshift_da.txt and redshift_Hz.txt " << endl;
|
|---|
| 162 | SLinInterp1D redshift2da;
|
|---|
| 163 | string interpfilename="redshift_da.txt";
|
|---|
| 164 | redshift2da.ReadXYFromFile(interpfilename);
|
|---|
| 165 | SLinInterp1D redshift2hz;
|
|---|
| 166 | interpfilename="redshift_Hz.txt";
|
|---|
| 167 | redshift2hz.ReadXYFromFile(interpfilename);
|
|---|
| 168 | comov_dA_z=redshift2da(z_Redshift);
|
|---|
| 169 | H_z=redshift2hz(z_Redshift);
|
|---|
| 170 | cout << " Redshift= " << z_Redshift << " -> ComovdA= " << comov_dA_z
|
|---|
| 171 | << " Mpc H(z)= " << H_z << " km/s/Mpc " << endl;
|
|---|
| 172 |
|
|---|
| 173 | H21Conversions conv;
|
|---|
| 174 | conv.setRedshift(z_Redshift);
|
|---|
| 175 | double lambda = conv.getLambda();
|
|---|
| 176 |
|
|---|
| 177 | double Dol = DIAMETRE/lambda;
|
|---|
| 178 |
|
|---|
| 179 | Four2DResponse arep(2, DIAMETRE/lambda, DIAMETRE/lambda, lambda);
|
|---|
| 180 | Four2DResponse* arep_p=&arep;
|
|---|
| 181 | Four2DRespTable resptbl;
|
|---|
| 182 | if (fgresptbl) {
|
|---|
| 183 | cout << "pknoise[1]: initializing Four2DRespTable from file" << resptblname << endl;
|
|---|
| 184 | resptbl.readFromPPF(resptblname);
|
|---|
| 185 | cout << "pknoise[1.b] Four2DRespTable.LambdaRef=" << resptbl.getLambdaRef() << " setLambda("
|
|---|
| 186 | << lambda << ")" << endl;
|
|---|
| 187 | resptbl.setLambda(lambda);
|
|---|
| 188 | arep_p=&resptbl;
|
|---|
| 189 | if (fgrenorm) {
|
|---|
| 190 | cout << "pknoise[1.c] call to resptbl.renormalize(" << rmax << ")";
|
|---|
| 191 | double omax=resptbl.renormalize(rmax);
|
|---|
| 192 | cout << " ... Old Max=" << omax << endl;
|
|---|
| 193 | }
|
|---|
| 194 | }
|
|---|
| 195 | else cout << " pknoise[1]: Four2DResponse ( Diameter=" << DIAMETRE << " Lambda= " << lambda
|
|---|
| 196 | << " DoL=" << DIAMETRE/lambda << " ) " << endl;
|
|---|
| 197 | Histo2D h2drep = arep_p->GetResponse();
|
|---|
| 198 | Histo2D h2drepuv = arep_p->GetResponseUV();
|
|---|
| 199 | double repmax= h2drep.VMax();
|
|---|
| 200 | if (fgautoscut) {
|
|---|
| 201 | SCut = FacSCut/repmax;
|
|---|
| 202 | cout << " pknoise[1.b]: Four2DResponse.RepMax=" << repmax << " --> SCut=" << FacSCut << "/repmax="
|
|---|
| 203 | << SCut << endl;
|
|---|
| 204 | }
|
|---|
| 205 | else cout << " pknoise[1.b]: Four2DResponse.RepMax=" << repmax << " , SCut=" << SCut << endl;
|
|---|
| 206 |
|
|---|
| 207 | cout << " pknoise[2]: Instanciating object type Four3DPk " << endl;
|
|---|
| 208 | RandomGenerator rg;
|
|---|
| 209 |
|
|---|
| 210 | double dfreq=H_z*cellsz[2]/(1+z_Redshift)/lambda/1000.;
|
|---|
| 211 | double freq0=conv.getFrequency()-dfreq*box3dsz[2]/2;
|
|---|
| 212 |
|
|---|
| 213 | double dkxmpc=2.*DeuxPI;
|
|---|
| 214 | double dkympc=2.*DeuxPI;
|
|---|
| 215 | double dkzmpc=2.*DeuxPI;
|
|---|
| 216 | double angscale=1.;
|
|---|
| 217 | Vector angscales(box3dsz[2]);
|
|---|
| 218 | double pnoise0=comov_dA_z*comov_dA_z*(1.+z_Redshift)*(1.+z_Redshift)/H_z;
|
|---|
| 219 | Vector pnoisevec(box3dsz[2]);
|
|---|
| 220 | pnoisevec=1.;
|
|---|
| 221 | angscales=angscale;
|
|---|
| 222 | if (NMAX<0) { box3dsz[0]=256; box3dsz[1]=256; box3dsz[2]=128; }
|
|---|
| 223 | else {
|
|---|
| 224 | angscale=comov_dA_z;
|
|---|
| 225 | for(int kz=0; kz<box3dsz[2]; kz++) angscales(kz)=redshift2da(1420.4/(freq0+(double)kz*dfreq)-1.);
|
|---|
| 226 | dkxmpc = DeuxPI/(double)box3dsz[0]/cellsz[0];
|
|---|
| 227 | dkympc = DeuxPI/(double)box3dsz[1]/cellsz[1];
|
|---|
| 228 | dkzmpc = DeuxPI/(double)box3dsz[2]/cellsz[2];
|
|---|
| 229 | }
|
|---|
| 230 | cout << " pknoise[2.b]: ComovDist=" << comov_dA_z << " Mpc H(z)=" << H_z << " Mpc/km/s -> angscale=" << angscale << endl;
|
|---|
| 231 | cout << " pknoise[2.c]: Freq0=" << freq0 << " dFreq=" << dfreq << " freq(z=" << z_Redshift << ")="
|
|---|
| 232 | << conv.getFrequency() << " zmin=" << 1420.4/freq0-1. << " zmax=" << 1420.4/(freq0+box3dsz[2]*dfreq)-1. << endl;
|
|---|
| 233 | if(fgnoisefreq) {
|
|---|
| 234 | for(int kz=0; kz<box3dsz[2]; kz++) {
|
|---|
| 235 | double zreds=1420.4/(freq0+(double)kz*dfreq)-1.;
|
|---|
| 236 | double comda=redshift2da(zreds);
|
|---|
| 237 | pnoisevec(kz)=comda*comda*(1.+zreds)*(1.+zreds)/redshift2hz(zreds)/pnoise0;
|
|---|
| 238 | }
|
|---|
| 239 | cout << " pknoise[2.d]: PnoiseFactor:" << pnoisevec(0) << " ... " << pnoisevec(pnoisevec.Size()-1) << endl;
|
|---|
| 240 | }
|
|---|
| 241 | cout << " pknoise[2.e]: Four3DPk m3d(rg," << box3dsz[0]/2 << "," << box3dsz[1] << ","
|
|---|
| 242 | << box3dsz[2] << ")" << endl;
|
|---|
| 243 | Four3DPk m3d(rg,box3dsz[0]/2,box3dsz[1],box3dsz[2]);
|
|---|
| 244 | cout << " pknoise[2.f]: m3d.SetCellSize(" << dkxmpc << "," << dkympc << "," << dkzmpc
|
|---|
| 245 | << ") cell size (Mpc) : " << cellsz[0] << "x" << cellsz[1] << "x" << cellsz[2] << endl;
|
|---|
| 246 | m3d.SetCellSize(dkxmpc, dkympc, dkzmpc);
|
|---|
| 247 | m3d.SetPrtLevel(prtlev,prtmod);
|
|---|
| 248 |
|
|---|
| 249 | if (NMAX>0) {
|
|---|
| 250 | cout << " pknoise[3]: Computing Noise P(k) using PkNoiseCalculator ..." << endl;
|
|---|
| 251 | PkNoiseCalculator pkn(m3d, *(arep_p), SCut, NMAX, tits.c_str());
|
|---|
| 252 | pkn.SetFreqRange(freq0, dfreq);
|
|---|
| 253 | pkn.SetAngScaleConversion(angscales);
|
|---|
| 254 | pkn.SetPNoiseFactor(pnoisevec);
|
|---|
| 255 | pkn.SetPrtLevel(prtlev,prtmod);
|
|---|
| 256 | HProf hpn = pkn.Compute(NBINPK);
|
|---|
| 257 | }
|
|---|
| 258 | else {
|
|---|
| 259 | cout << " pknoise[3]: Computing Noise P(k) m3d.ComputeNoisePk(...) " << endl;
|
|---|
| 260 | HProf hpn = m3d.ComputeNoisePk(*(arep_p),angscale,SCut,NBINPK);
|
|---|
| 261 | }
|
|---|
| 262 |
|
|---|
| 263 | HProf hpnoise=m3d.GetPk();
|
|---|
| 264 | DataTable dtnoise;
|
|---|
| 265 | Histo fracmodok=m3d.FillPkDataTable(dtnoise);
|
|---|
| 266 | HProf h1dnoise=arep_p->GetProjNoiseLevel();
|
|---|
| 267 | HProf h1drep=arep_p->GetProjResponse();
|
|---|
| 268 | cout << " pknoise[3.b]: writing dtnoise,hpn,h2rep... with tags to " << outfile << endl;
|
|---|
| 269 | po << PPFNameTag("dtnoise") << dtnoise;
|
|---|
| 270 | po << PPFNameTag("hpnoise") << hpnoise;
|
|---|
| 271 | po << PPFNameTag("fracmodok") << fracmodok;
|
|---|
| 272 | po << PPFNameTag("h1dnoise") << h1dnoise;
|
|---|
| 273 | po << PPFNameTag("h1drep") << h1drep;
|
|---|
| 274 | po << PPFNameTag("h2drep") << h2drep;
|
|---|
| 275 | po << PPFNameTag("h2drepuv") << h2drepuv;
|
|---|
| 276 |
|
|---|
| 277 | rc = 0;
|
|---|
| 278 | } // End of try bloc
|
|---|
| 279 | catch (PThrowable & exc) { // catching SOPHYA exceptions
|
|---|
| 280 | cerr << " pknoise.cc: Catched Exception (PThrowable)" << (string)typeid(exc).name()
|
|---|
| 281 | << "\n...exc.Msg= " << exc.Msg() << endl;
|
|---|
| 282 | rc = 99;
|
|---|
| 283 | }
|
|---|
| 284 | catch (std::exception & e) { // catching standard C++ exceptions
|
|---|
| 285 | cerr << " pknoise.cc: Catched std::exception " << " - what()= " << e.what() << endl;
|
|---|
| 286 | rc = 98;
|
|---|
| 287 | }
|
|---|
| 288 | catch (...) { // catching other exceptions
|
|---|
| 289 | cerr << " pknoise.cc: some other exception (...) was caught ! " << endl;
|
|---|
| 290 | rc = 97;
|
|---|
| 291 | }
|
|---|
| 292 | PrtTim("End-pknoise");
|
|---|
| 293 | cout << " ==== End of pknoise.cc program Rc= " << rc << endl;
|
|---|
| 294 | return rc;
|
|---|
| 295 | }
|
|---|
| 296 |
|
|---|
| 297 |
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
