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

Last change on this file since 3971 was 3948, checked in by ansari, 15 years ago

modifications de config cylindres pour papier, Reza 16/02/2011

File size: 7.0 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 , -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
33typedef DR48RandGen RandomGenerator ;
34
35// ---------------------------------------------------------------------
36// Test main
37// R. Ansari - Avril-Dec 2010
38// ---------------------------------------------------------------------
39
40void 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 << " -z redshift (default=0.7) \n"
48 << " -scut SCutValue (default= -100.) \n"
49 << " if SCutValue<0. ==> SCut=MinNoisePower*(-SCutValue) \n"
50 << " -prt PrtLev,PrtModulo (default=0,10) " << endl;
51 return;
52}
53
54//-------------------------------------------------------------------------
55// ------------------ MAIN PROGRAM ------------------------------
56//-------------------------------------------------------------------------
57int main(int narg, const char* arg[])
58{
59 if ( (narg<3)||((narg>1)&&(strcmp(arg[1],"-h")==0)) ) {
60 Usage();
61 return 1;
62 }
63 cout << " ==== pknoise.cc : interferometer noise power spectrum computation ==== " << endl;
64 // make sure SOPHYA modules are initialized
65 SophyaInit();
66 // FitsIOServerInit();
67 InitTim();
68 //--- decoding command line arguments
69 string tits="pknoise";
70 char cbuff[64];
71 bool fgresptbl=false;
72 double DIAMETRE=100.;
73 string resptblname;
74 double NoiseLevel=1.;
75
76 bool fgrenorm=false;
77 double rmax=1.;
78 int NMAX = 0;
79 double SCut=-100.;
80 bool fgautoscut=true;
81 double FacSCut=-SCut;
82 double z_Redshift=0.7 ; // 21 cm at z=0.7 -> 0.357 m
83 int prtlev=0;
84 int prtmod=10;
85
86 int ka=1;
87 while (ka<(narg-1)) {
88 if (strcmp(arg[ka],"-noise")==0) {
89 NoiseLevel=atof(arg[ka+1]);
90 ka+=2;
91 }
92 else if (strcmp(arg[ka],"-renmax")==0) {
93 rmax=atof(arg[ka+1]); fgrenorm=true; ka+=2;
94 }
95 else if (strcmp(arg[ka],"-scut")==0) {
96 SCut=atof(arg[ka+1]); ka+=2;
97 if (SCut<0.) { FacSCut=-SCut; fgautoscut=true; }
98 }
99 else if (strcmp(arg[ka],"-ngen")==0) {
100 NMAX=atoi(arg[ka+1]); ka+=2;
101 }
102 else if (strcmp(arg[ka],"-z")==0) {
103 z_Redshift=atof(arg[ka+1]); ka+=2;
104 }
105 else if (strcmp(arg[ka],"-prt")==0) {
106 sscanf(arg[ka+1],"%d,%d",&prtlev,&prtmod); ka+=2;
107 }
108 else break;
109 }
110
111 if ((ka+1)>=narg) {
112 cout << " pknoise / Argument error " << endl;
113 Usage();
114 return 2;
115 }
116 if (isdigit(*arg[ka])) {
117 fgresptbl=false;
118 DIAMETRE=atof(arg[ka]);
119 sprintf(cbuff,"pknoise_Dish(%g m)", DIAMETRE);
120 }
121 else {
122 resptblname=arg[ka];
123 fgresptbl=true;
124 sprintf(cbuff,"pknoise_RespTblName=%s", arg[ka]);
125 }
126 tits=cbuff;
127 string outfile = arg[ka+1];
128 if (outfile==".") outfile = "pknoise.ppf";
129 //-- end command line arguments
130
131 int rc = 1;
132 try { // exception handling try bloc at top level
133 cout << " pknoise[0] : Executing, output file= " << outfile << endl;
134 POutPersist po(outfile);
135
136 H21Conversions conv;
137 conv.setRedshift(z_Redshift);
138 double lambda = conv.getLambda();
139
140 double Dol = DIAMETRE/lambda;
141
142 Four2DResponse arep(2, DIAMETRE/lambda, DIAMETRE/lambda, lambda);
143 Four2DResponse* arep_p=&arep;
144 Four2DRespTable resptbl;
145 if (fgresptbl) {
146 cout << "pknoise[1]: initializing Four2DRespTable from file" << resptblname << endl;
147 resptbl.readFromPPF(resptblname);
148 cout << "pknoise[1.b] Four2DRespTable.LambdaRef=" << resptbl.getLambdaRef() << " setLambda("
149 << lambda << ")" << endl;
150 resptbl.setLambda(lambda);
151 arep_p=&resptbl;
152 if (fgrenorm) {
153 cout << "pknoise[1.c] call to resptbl.renormalize(" << rmax << ")";
154 double omax=resptbl.renormalize(rmax);
155 cout << " ... Old Max=" << omax << endl;
156 }
157 }
158 else cout << " pknoise[1]: Four2DResponse ( Diameter=" << DIAMETRE << " Lambda= " << lambda
159 << " DoL=" << DIAMETRE/lambda << " ) " << endl;
160 Histo2D h2drep = arep_p->GetResponse();
161 double repmax= h2drep.VMax();
162 if (fgautoscut) {
163 SCut = FacSCut/repmax;
164 cout << " pknoise[1.b]: Four2DResponse.RepMax=" << repmax << " --> SCut=" << FacSCut << "/repmax="
165 << SCut << endl;
166 }
167 else cout << " pknoise[1.b]: Four2DResponse.RepMax=" << repmax << " , SCut=" << SCut << endl;
168
169 cout << " pknoise[2]: Instanciating object type Four3DPk " << endl;
170 RandomGenerator rg;
171 Four3DPk m3d(rg);
172 m3d.SetCellSize(2.*DeuxPI, 2.*DeuxPI, 2.*DeuxPI);
173 m3d.SetPrtLevel(prtlev,prtmod);
174
175 cout << " pknoise[3]: Computing Noise P(k) using PkNoiseCalculator ..." << endl;
176 if (NMAX>0) {
177 PkNoiseCalculator pkn(m3d, *(arep_p), SCut, NMAX, tits.c_str());
178 pkn.SetPrtLevel(prtlev,prtmod);
179 HProf hpn = pkn.Compute();
180 cout << " pknoise[3.b]: writing hpn noise profile to " << outfile << endl;
181 po << hpn;
182 }
183 else {
184 Histo fracmodok;
185 DataTable dtnoise;
186 HProf hpn = m3d.ComputeNoisePk(*(arep_p),fracmodok,dtnoise,SCut);
187 HProf h1dnoise=arep_p->GetProjNoiseLevel();
188 HProf h1drep=arep_p->GetProjResponse();
189 cout << " pknoise[3.b]: writing dtnoise,hpn,h2rep... with tags to " << outfile << endl;
190 po << PPFNameTag("dtnoise") << dtnoise;
191 po << PPFNameTag("hpnoise") << hpn;
192 po << PPFNameTag("fracmodok") << fracmodok;
193 po << PPFNameTag("h1dnoise") << h1dnoise;
194 po << PPFNameTag("h1drep") << h1drep;
195 po << PPFNameTag("h2drep") << h2drep;
196 }
197 rc = 0;
198 } // End of try bloc
199 catch (PThrowable & exc) { // catching SOPHYA exceptions
200 cerr << " pknoise.cc: Catched Exception (PThrowable)" << (string)typeid(exc).name()
201 << "\n...exc.Msg= " << exc.Msg() << endl;
202 rc = 99;
203 }
204 catch (std::exception & e) { // catching standard C++ exceptions
205 cerr << " pknoise.cc: Catched std::exception " << " - what()= " << e.what() << endl;
206 rc = 98;
207 }
208 catch (...) { // catching other exceptions
209 cerr << " pknoise.cc: some other exception (...) was caught ! " << endl;
210 rc = 97;
211 }
212 PrtTim("End-pknoise");
213 cout << " ==== End of pknoise.cc program Rc= " << rc << endl;
214 return rc;
215}
216
217
Note: See TracBrowser for help on using the repository browser.