source: Sophya/trunk/Cosmo/RadioBeam/pknoise.cc@ 4044

Last change on this file since 4044 was 4030, checked in by ansari, 14 years ago

Corrections papiers avec les commentaires du referee (1) - Reza 26/10/2011

File size: 10.7 KB
Line 
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
34typedef DR48RandGen RandomGenerator ;
35
36// ---------------------------------------------------------------------
37// Test main
38// R. Ansari - Avril-Dec 2010
39// ---------------------------------------------------------------------
40
41void 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//-------------------------------------------------------------------------
64int 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
Note: See TracBrowser for help on using the repository browser.