| 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 , -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 |
|
|---|
| 33 | typedef DR48RandGen RandomGenerator ;
|
|---|
| 34 |
|
|---|
| 35 | // ---------------------------------------------------------------------
|
|---|
| 36 | // Test main
|
|---|
| 37 | // R. Ansari - Avril-Dec 2010
|
|---|
| 38 | // ---------------------------------------------------------------------
|
|---|
| 39 |
|
|---|
| 40 | void Usage()
|
|---|
| 41 | {
|
|---|
| 42 | cout << " Usage: pknoise [-parname value] Diameter/Four2DRespTableFile OutPPFName \n"
|
|---|
| 43 | << " -noise NoiseLevel (default=1.) \n"
|
|---|
| 44 | << " -renmax MaxValue (default : Do NOT renormalize 2D response value \n"
|
|---|
| 45 | << " -ngen NGen (default=0) number of noise fourier amp generations \n"
|
|---|
| 46 | << " NGen=0 -> Call ComputeNoisePk(), else generate Fourier Amplitudes (random) \n"
|
|---|
| 47 | << " NGen=-1 -> Call ComputeNoisePk(), dont change Four3DPk cell size \n"
|
|---|
| 48 | << " -bsize sx,sy,sz : 3D real space box size (default=512x512x256) \n"
|
|---|
| 49 | << " -cszmpc sx,sy,sz : 3D real space box cell size in Mpc (default=3x3x3) \n"
|
|---|
| 50 | << " -z z,dA,H_z : redshift,ComovDist(z),H(z) (default=1.0,3330,120.5) \n"
|
|---|
| 51 | << " -scut SCutValue (default= -100.) \n"
|
|---|
| 52 | << " if SCutValue<0. ==> SCut=MinNoisePower*(-SCutValue) \n"
|
|---|
| 53 | << " -prt PrtLev,PrtModulo (default=0,10) " << endl;
|
|---|
| 54 | return;
|
|---|
| 55 | }
|
|---|
| 56 |
|
|---|
| 57 | //-------------------------------------------------------------------------
|
|---|
| 58 | // ------------------ MAIN PROGRAM ------------------------------
|
|---|
| 59 | //-------------------------------------------------------------------------
|
|---|
| 60 | int main(int narg, const char* arg[])
|
|---|
| 61 | {
|
|---|
| 62 | if ( (narg<3)||((narg>1)&&(strcmp(arg[1],"-h")==0)) ) {
|
|---|
| 63 | Usage();
|
|---|
| 64 | return 1;
|
|---|
| 65 | }
|
|---|
| 66 | cout << " ==== pknoise.cc : interferometer noise power spectrum computation ==== " << endl;
|
|---|
| 67 | // make sure SOPHYA modules are initialized
|
|---|
| 68 | SophyaInit();
|
|---|
| 69 | // FitsIOServerInit();
|
|---|
| 70 | InitTim();
|
|---|
| 71 | //--- decoding command line arguments
|
|---|
| 72 | string tits="pknoise";
|
|---|
| 73 | char cbuff[64];
|
|---|
| 74 | bool fgresptbl=false;
|
|---|
| 75 | double DIAMETRE=100.;
|
|---|
| 76 | string resptblname;
|
|---|
| 77 | double NoiseLevel=1.;
|
|---|
| 78 |
|
|---|
| 79 | bool fgrenorm=false;
|
|---|
| 80 | double rmax=1.;
|
|---|
| 81 | int NMAX = 0;
|
|---|
| 82 | double SCut=-100.;
|
|---|
| 83 | bool fgautoscut=true;
|
|---|
| 84 | double FacSCut=-SCut;
|
|---|
| 85 | double z_Redshift=1.0 ; // z=0.7 : 21 cm at z=0.7 -> 0.357 m
|
|---|
| 86 | double comov_dA_z=3330.; // Comoving radial distance = (1+z)dA
|
|---|
| 87 | double H_z=120.5; // Hubble_param(z)
|
|---|
| 88 | int box3dsz[3]={512,512,256};
|
|---|
| 89 | double cellsz[3]={3.5,3.5,3.5};
|
|---|
| 90 |
|
|---|
| 91 | int prtlev=0;
|
|---|
| 92 | int prtmod=10;
|
|---|
| 93 |
|
|---|
| 94 | int ka=1;
|
|---|
| 95 | while (ka<(narg-1)) {
|
|---|
| 96 | if (strcmp(arg[ka],"-noise")==0) {
|
|---|
| 97 | NoiseLevel=atof(arg[ka+1]);
|
|---|
| 98 | ka+=2;
|
|---|
| 99 | }
|
|---|
| 100 | else if (strcmp(arg[ka],"-renmax")==0) {
|
|---|
| 101 | rmax=atof(arg[ka+1]); fgrenorm=true; ka+=2;
|
|---|
| 102 | }
|
|---|
| 103 | else if (strcmp(arg[ka],"-scut")==0) {
|
|---|
| 104 | SCut=atof(arg[ka+1]); ka+=2;
|
|---|
| 105 | if (SCut<0.) { FacSCut=-SCut; fgautoscut=true; }
|
|---|
| 106 | }
|
|---|
| 107 | else if (strcmp(arg[ka],"-ngen")==0) {
|
|---|
| 108 | NMAX=atoi(arg[ka+1]); ka+=2;
|
|---|
| 109 | }
|
|---|
| 110 | else if (strcmp(arg[ka],"-z")==0) {
|
|---|
| 111 | sscanf(arg[ka+1],"%lg,%lg,%lg",&z_Redshift,&comov_dA_z,&H_z); ka+=2;
|
|---|
| 112 | }
|
|---|
| 113 | else if (strcmp(arg[ka],"-bsize")==0) {
|
|---|
| 114 | sscanf(arg[ka+1],"%d,%d,%d",box3dsz,box3dsz+1,box3dsz+2); ka+=2;
|
|---|
| 115 | }
|
|---|
| 116 | else if (strcmp(arg[ka],"-cszmpc")==0) {
|
|---|
| 117 | sscanf(arg[ka+1],"%lg,%lg,%lg",cellsz,cellsz+1,cellsz+2); ka+=2;
|
|---|
| 118 | }
|
|---|
| 119 | else if (strcmp(arg[ka],"-prt")==0) {
|
|---|
| 120 | sscanf(arg[ka+1],"%d,%d",&prtlev,&prtmod); ka+=2;
|
|---|
| 121 | }
|
|---|
| 122 | else break;
|
|---|
| 123 | }
|
|---|
| 124 |
|
|---|
| 125 | if ((ka+1)>=narg) {
|
|---|
| 126 | cout << " pknoise / Argument error " << endl;
|
|---|
| 127 | Usage();
|
|---|
| 128 | return 2;
|
|---|
| 129 | }
|
|---|
| 130 | if (isdigit(*arg[ka])) {
|
|---|
| 131 | fgresptbl=false;
|
|---|
| 132 | DIAMETRE=atof(arg[ka]);
|
|---|
| 133 | sprintf(cbuff,"pknoise_Dish(%g m)", DIAMETRE);
|
|---|
| 134 | }
|
|---|
| 135 | else {
|
|---|
| 136 | resptblname=arg[ka];
|
|---|
| 137 | fgresptbl=true;
|
|---|
| 138 | sprintf(cbuff,"pknoise_RespTblName=%s", arg[ka]);
|
|---|
| 139 | }
|
|---|
| 140 | tits=cbuff;
|
|---|
| 141 | string outfile = arg[ka+1];
|
|---|
| 142 | if (outfile==".") outfile = "pknoise.ppf";
|
|---|
| 143 | //-- end command line arguments
|
|---|
| 144 |
|
|---|
| 145 | int rc = 1;
|
|---|
| 146 | try { // exception handling try bloc at top level
|
|---|
| 147 | cout << " pknoise[0] : Executing, output file= " << outfile << endl;
|
|---|
| 148 | POutPersist po(outfile);
|
|---|
| 149 |
|
|---|
| 150 | H21Conversions conv;
|
|---|
| 151 | conv.setRedshift(z_Redshift);
|
|---|
| 152 | double lambda = conv.getLambda();
|
|---|
| 153 |
|
|---|
| 154 | double Dol = DIAMETRE/lambda;
|
|---|
| 155 |
|
|---|
| 156 | Four2DResponse arep(2, DIAMETRE/lambda, DIAMETRE/lambda, lambda);
|
|---|
| 157 | Four2DResponse* arep_p=&arep;
|
|---|
| 158 | Four2DRespTable resptbl;
|
|---|
| 159 | if (fgresptbl) {
|
|---|
| 160 | cout << "pknoise[1]: initializing Four2DRespTable from file" << resptblname << endl;
|
|---|
| 161 | resptbl.readFromPPF(resptblname);
|
|---|
| 162 | cout << "pknoise[1.b] Four2DRespTable.LambdaRef=" << resptbl.getLambdaRef() << " setLambda("
|
|---|
| 163 | << lambda << ")" << endl;
|
|---|
| 164 | resptbl.setLambda(lambda);
|
|---|
| 165 | arep_p=&resptbl;
|
|---|
| 166 | if (fgrenorm) {
|
|---|
| 167 | cout << "pknoise[1.c] call to resptbl.renormalize(" << rmax << ")";
|
|---|
| 168 | double omax=resptbl.renormalize(rmax);
|
|---|
| 169 | cout << " ... Old Max=" << omax << endl;
|
|---|
| 170 | }
|
|---|
| 171 | }
|
|---|
| 172 | else cout << " pknoise[1]: Four2DResponse ( Diameter=" << DIAMETRE << " Lambda= " << lambda
|
|---|
| 173 | << " DoL=" << DIAMETRE/lambda << " ) " << endl;
|
|---|
| 174 | Histo2D h2drep = arep_p->GetResponse();
|
|---|
| 175 | double repmax= h2drep.VMax();
|
|---|
| 176 | if (fgautoscut) {
|
|---|
| 177 | SCut = FacSCut/repmax;
|
|---|
| 178 | cout << " pknoise[1.b]: Four2DResponse.RepMax=" << repmax << " --> SCut=" << FacSCut << "/repmax="
|
|---|
| 179 | << SCut << endl;
|
|---|
| 180 | }
|
|---|
| 181 | else cout << " pknoise[1.b]: Four2DResponse.RepMax=" << repmax << " , SCut=" << SCut << endl;
|
|---|
| 182 |
|
|---|
| 183 | cout << " pknoise[2]: Instanciating object type Four3DPk " << endl;
|
|---|
| 184 | RandomGenerator rg;
|
|---|
| 185 |
|
|---|
| 186 | double dkxmpc=2.*DeuxPI;
|
|---|
| 187 | double dkympc=2.*DeuxPI;
|
|---|
| 188 | double dkzmpc=2.*DeuxPI;
|
|---|
| 189 | double angscale=1.;
|
|---|
| 190 | if (NMAX<0) { box3dsz[0]=256; box3dsz[1]=256; box3dsz[3]=128; }
|
|---|
| 191 | else {
|
|---|
| 192 | angscale=comov_dA_z;
|
|---|
| 193 | dkxmpc = DeuxPI/(double)box3dsz[0]/cellsz[0];
|
|---|
| 194 | dkympc = DeuxPI/(double)box3dsz[1]/cellsz[1];
|
|---|
| 195 | dkzmpc = DeuxPI/(double)box3dsz[2]/cellsz[2];
|
|---|
| 196 | }
|
|---|
| 197 | cout << " pknoise[2.b]: ComovDist" << comov_dA_z << " Mpc H(z)=" << H_z << " Mpc/km/s -> angscale=" << angscale << endl;
|
|---|
| 198 | cout << " pknoise[2.c]: Four3DPk m3d(rg," << box3dsz[0]/2 << "," << box3dsz[1] << ","
|
|---|
| 199 | << box3dsz[2] << ")" << endl;
|
|---|
| 200 | Four3DPk m3d(rg,box3dsz[0]/2,box3dsz[1],box3dsz[2]);
|
|---|
| 201 | cout << " pknoise[2.d]: m3d.SetCellSize(" << dkxmpc << "," << dkympc << "," << dkzmpc << endl;
|
|---|
| 202 | m3d.SetCellSize(dkxmpc, dkympc, dkzmpc);
|
|---|
| 203 | m3d.SetPrtLevel(prtlev,prtmod);
|
|---|
| 204 |
|
|---|
| 205 | cout << " pknoise[3]: Computing Noise P(k) using PkNoiseCalculator ..." << endl;
|
|---|
| 206 | if (NMAX>0) {
|
|---|
| 207 | PkNoiseCalculator pkn(m3d, *(arep_p), SCut, NMAX, tits.c_str());
|
|---|
| 208 | double dfreq=H_z*cellsz[2]/(1+z_Redshift)/lambda/1000.;
|
|---|
| 209 | double freq0=conv.getFrequency()-dfreq*box3dsz[2]/2;
|
|---|
| 210 | cout << " pknoise[3.b]: Freq0=" << freq0 << " dFreq=" << dfreq << " freq(z=" << z_Redshift << ")="
|
|---|
| 211 | << conv.getFrequency() << " AngScale=" << angscale << endl;
|
|---|
| 212 | pkn.SetFreqRange(freq0, dfreq);
|
|---|
| 213 | pkn.SetAngScaleConversion(angscale);
|
|---|
| 214 | pkn.SetPrtLevel(prtlev,prtmod);
|
|---|
| 215 | HProf hpn = pkn.Compute();
|
|---|
| 216 | cout << " pknoise[3.c]: writing hpn noise profile to " << outfile << endl;
|
|---|
| 217 | po << hpn;
|
|---|
| 218 | }
|
|---|
| 219 | else {
|
|---|
| 220 | Histo fracmodok;
|
|---|
| 221 | DataTable dtnoise;
|
|---|
| 222 | HProf hpn = m3d.ComputeNoisePk(*(arep_p),fracmodok,dtnoise,angscale,SCut);
|
|---|
| 223 | HProf h1dnoise=arep_p->GetProjNoiseLevel();
|
|---|
| 224 | HProf h1drep=arep_p->GetProjResponse();
|
|---|
| 225 | cout << " pknoise[3.b]: writing dtnoise,hpn,h2rep... with tags to " << outfile << endl;
|
|---|
| 226 | po << PPFNameTag("dtnoise") << dtnoise;
|
|---|
| 227 | po << PPFNameTag("hpnoise") << hpn;
|
|---|
| 228 | po << PPFNameTag("fracmodok") << fracmodok;
|
|---|
| 229 | po << PPFNameTag("h1dnoise") << h1dnoise;
|
|---|
| 230 | po << PPFNameTag("h1drep") << h1drep;
|
|---|
| 231 | po << PPFNameTag("h2drep") << h2drep;
|
|---|
| 232 | }
|
|---|
| 233 | rc = 0;
|
|---|
| 234 | } // End of try bloc
|
|---|
| 235 | catch (PThrowable & exc) { // catching SOPHYA exceptions
|
|---|
| 236 | cerr << " pknoise.cc: Catched Exception (PThrowable)" << (string)typeid(exc).name()
|
|---|
| 237 | << "\n...exc.Msg= " << exc.Msg() << endl;
|
|---|
| 238 | rc = 99;
|
|---|
| 239 | }
|
|---|
| 240 | catch (std::exception & e) { // catching standard C++ exceptions
|
|---|
| 241 | cerr << " pknoise.cc: Catched std::exception " << " - what()= " << e.what() << endl;
|
|---|
| 242 | rc = 98;
|
|---|
| 243 | }
|
|---|
| 244 | catch (...) { // catching other exceptions
|
|---|
| 245 | cerr << " pknoise.cc: some other exception (...) was caught ! " << endl;
|
|---|
| 246 | rc = 97;
|
|---|
| 247 | }
|
|---|
| 248 | PrtTim("End-pknoise");
|
|---|
| 249 | cout << " ==== End of pknoise.cc program Rc= " << rc << endl;
|
|---|
| 250 | return rc;
|
|---|
| 251 | }
|
|---|
| 252 |
|
|---|
| 253 |
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
