| 1 | #include "machdefs.h"
|
|---|
| 2 | #include "sopnamsp.h"
|
|---|
| 3 | #include <iostream>
|
|---|
| 4 | #include <string>
|
|---|
| 5 | #include <math.h>
|
|---|
| 6 |
|
|---|
| 7 | #include <typeinfo>
|
|---|
| 8 |
|
|---|
| 9 | #include "mdish.h"
|
|---|
| 10 | #include "specpk.h"
|
|---|
| 11 | #include "interfconfigs.h"
|
|---|
| 12 |
|
|---|
| 13 | #include "histinit.h"
|
|---|
| 14 | // #include "fiosinit.h"
|
|---|
| 15 | // #include "fitsioserver.h"
|
|---|
| 16 |
|
|---|
| 17 | #include "randr48.h"
|
|---|
| 18 |
|
|---|
| 19 | #include "timing.h"
|
|---|
| 20 | #include "ctimer.h"
|
|---|
| 21 |
|
|---|
| 22 | typedef DR48RandGen RandomGenerator ;
|
|---|
| 23 |
|
|---|
| 24 | // ---------------------------------------------------------------------
|
|---|
| 25 | // Test main program for computing interferometer noise power spectrum
|
|---|
| 26 | // R. Ansari - Avril 2010
|
|---|
| 27 | // ---------------------------------------------------------------------
|
|---|
| 28 |
|
|---|
| 29 | class PkNoiseCalculator
|
|---|
| 30 | {
|
|---|
| 31 | public:
|
|---|
| 32 | PkNoiseCalculator(Four3DPk& pk3, Four2DResponse& rep,
|
|---|
| 33 | double s2cut=100., int ngen=1, const char* tit="PkNoise")
|
|---|
| 34 | : pkn3d(pk3), frep(rep), S2CUT(s2cut), NGEN(ngen), title(tit)
|
|---|
| 35 | { }
|
|---|
| 36 | inline void SetS2Cut(double s2cut=100.)
|
|---|
| 37 | { S2CUT=s2cut; }
|
|---|
| 38 | HProf Compute()
|
|---|
| 39 | {
|
|---|
| 40 | Timer tm(title.c_str());
|
|---|
| 41 | tm.Nop();
|
|---|
| 42 | HProf hnd;
|
|---|
| 43 | cout << "PkNoiseCalculator::Compute() " << title << " NGEN=" << NGEN << " S2CUT=" << S2CUT << endl;
|
|---|
| 44 | for(int igen=0; igen<NGEN; igen++) {
|
|---|
| 45 | pkn3d.ComputeNoiseFourierAmp(frep);
|
|---|
| 46 | if (igen==0) hnd = pkn3d.ComputePk(S2CUT);
|
|---|
| 47 | else pkn3d.ComputePkCumul(hnd,S2CUT);
|
|---|
| 48 | }
|
|---|
| 49 | return hnd;
|
|---|
| 50 | }
|
|---|
| 51 |
|
|---|
| 52 | Four3DPk& pkn3d;
|
|---|
| 53 | Four2DResponse& frep;
|
|---|
| 54 | double S2CUT;
|
|---|
| 55 | int NGEN;
|
|---|
| 56 | string title;
|
|---|
| 57 | };
|
|---|
| 58 |
|
|---|
| 59 |
|
|---|
| 60 |
|
|---|
| 61 | //-------------------------------------------------------------------------
|
|---|
| 62 | // ------------------ MAIN PROGRAM ------------------------------
|
|---|
| 63 | //-------------------------------------------------------------------------
|
|---|
| 64 | int main(int narg, const char* arg[])
|
|---|
| 65 | {
|
|---|
| 66 | if ((narg>1)&&(strcmp(arg[1],"-h")==0)) {
|
|---|
| 67 | cout<< " Usage: pknoise [OutPPFName NGen S2Cut Lambda] " << endl;
|
|---|
| 68 | cout<< " Default: OutPPFName=pknoise.ppf, NGen=1 " << endl;
|
|---|
| 69 | cout<< " S2CUT=0. , Lambda=0.357 \n" << endl;
|
|---|
| 70 |
|
|---|
| 71 | return 1;
|
|---|
| 72 | }
|
|---|
| 73 | cout << " ==== pknoise.cc program , test of SpectralShape and MassDist2D classes ==== " << endl;
|
|---|
| 74 | // make sure SOPHYA modules are initialized
|
|---|
| 75 | SophyaInit();
|
|---|
| 76 | // FitsIOServerInit();
|
|---|
| 77 | InitTim();
|
|---|
| 78 | //--- decoding command line arguments
|
|---|
| 79 | string outfile = "pknoise.ppf";
|
|---|
| 80 | if (narg>1) outfile = arg[1];
|
|---|
| 81 | if (outfile==".") outfile = "pknoise.ppf";
|
|---|
| 82 | int NMAX = 1;
|
|---|
| 83 | if (narg>2) NMAX = atoi(arg[2]);
|
|---|
| 84 | double SCut=0.;
|
|---|
| 85 | if (narg>3) SCut = atof(arg[3]);
|
|---|
| 86 | double LAMBDA=0.357 ; // 21 cm at z=0.7
|
|---|
| 87 | if (narg>4) LAMBDA = atof(arg[4]);
|
|---|
| 88 |
|
|---|
| 89 | //-- end command line arguments
|
|---|
| 90 |
|
|---|
| 91 | int rc = 1;
|
|---|
| 92 | try { // exception handling try bloc at top level
|
|---|
| 93 | cout << "0/ pknoise.cc: Executing, output file= " << outfile << endl;
|
|---|
| 94 | POutPersist po(outfile);
|
|---|
| 95 | cout << " 1.a/ Instanciating object type SpectralShape " << endl;
|
|---|
| 96 | SpectralShape spec(2);
|
|---|
| 97 | cout << " 1.b/ Wrinting spectral shape vector (name= Pk) to output PPF " << endl;
|
|---|
| 98 | Histo hpk = spec.GetPk(1024);
|
|---|
| 99 | po << PPFNameTag("Pk") << hpk;
|
|---|
| 100 |
|
|---|
| 101 | double D = 100.;
|
|---|
| 102 | double lambda = LAMBDA;
|
|---|
| 103 | double Dol = D/lambda;
|
|---|
| 104 | cout << " 2.a/ Instanciating Four2DResponse(1/2/3...) " << endl;
|
|---|
| 105 | Four2DResponse dishg(1,Dol,Dol);
|
|---|
| 106 | Four2DResponse dish(2,Dol,Dol);
|
|---|
| 107 | Four2DResponse dish2(2,Dol*2.,Dol*2.);
|
|---|
| 108 | Four2DResponse dishsq(3,Dol,Dol/5.);
|
|---|
| 109 | cout << " 2.b/ Writing Four2DResponse Histo2D to output ppf " << endl;
|
|---|
| 110 | Histo2D hdg = dishg.GetResponse();
|
|---|
| 111 | Histo2D hd = dish.GetResponse();
|
|---|
| 112 | Histo2D hd2 = dish2.GetResponse();
|
|---|
| 113 | Histo2D hdsq = dishsq.GetResponse();
|
|---|
| 114 | po << PPFNameTag("dishg") << hdg;
|
|---|
| 115 | po << PPFNameTag("dish") << hd;
|
|---|
| 116 | po << PPFNameTag("dish2") << hd2;
|
|---|
| 117 | po << PPFNameTag("dishsq") << hdsq;
|
|---|
| 118 |
|
|---|
| 119 | cout << " 2.c/ Creating MultiDish Filled Array " << endl;
|
|---|
| 120 | double Ddish = 5.;
|
|---|
| 121 | double Ddish2 = 7.5;
|
|---|
| 122 | double Eta=0.95;
|
|---|
| 123 | int cnt=0;
|
|---|
| 124 | vector<Dish> vdplein = CreateFilledSqConfig(20, Ddish, Eta);
|
|---|
| 125 | vector<Dish> vdpl64 = CreateFilledSqConfig(8, Ddish, Eta);
|
|---|
| 126 |
|
|---|
| 127 | vector<Dish> vdsparse = CreateConfigA(Ddish, Eta);
|
|---|
| 128 | vector<Dish> vdsparseD7 = CreateConfigA(Ddish2, Eta);
|
|---|
| 129 | // vector<Dish> vdsparseB = CreateConfigB(Ddish, Eta);
|
|---|
| 130 | vector<Dish> vdsparseB = CreateConfigB(Ddish, Eta);
|
|---|
| 131 | vector<Dish> vdsparseC = CreateConfigC(Ddish, Eta);
|
|---|
| 132 |
|
|---|
| 133 |
|
|---|
| 134 | double cylW=12.; // Largeur des cylindres
|
|---|
| 135 | double cylRL=0.5; // Longeur des elements de reception le long du cylindre
|
|---|
| 136 | double etaW=0.95; // Efficacite de couverture en largeur
|
|---|
| 137 | double etaRL=0.9; // Efficacite de couverture le long du cylindre
|
|---|
| 138 | vector<Dish> vcylplein = CreateFilledCylConfig(8, 192, cylW, cylRL, etaW, etaRL, true);
|
|---|
| 139 | vector<Dish> vcylplP = CreateFilledCylConfig(8, 192, cylW, cylRL, etaW, etaRL, false);
|
|---|
| 140 |
|
|---|
| 141 | cylW=10.;
|
|---|
| 142 | cylRL=0.5;
|
|---|
| 143 | vector<Dish> v3cyl = CreateFilledCylConfig(3, 128, cylW, cylRL, etaW, etaRL, true);
|
|---|
| 144 | vector<Dish> v3cylP = CreateFilledCylConfig(3, 128, cylW, cylRL, etaW, etaRL, false);
|
|---|
| 145 | cylW=25.;
|
|---|
| 146 | cylRL=0.5;
|
|---|
| 147 | etaW=0.3;
|
|---|
| 148 | etaRL=0.9;
|
|---|
| 149 | vector<Dish> v2cyl = CreateFilledCylConfig(2, 32, cylW, cylRL, etaW, etaRL, true);
|
|---|
| 150 | vector<Dish> v2cylP = CreateFilledCylConfig(2, 32, cylW, cylRL, etaW, etaRL, false);
|
|---|
| 151 |
|
|---|
| 152 | double LMAX = D;
|
|---|
| 153 | bool fgnoauto = true;
|
|---|
| 154 | int NRX=100;
|
|---|
| 155 | int NRY=100;
|
|---|
| 156 |
|
|---|
| 157 | MultiDish mdfill(lambda, LMAX, vdplein, fgnoauto);
|
|---|
| 158 | mdfill.SetRespHisNBins(NRX,NRY);
|
|---|
| 159 | Histo2D hrfill = mdfill.GetResponse();
|
|---|
| 160 | PrtTim("Apres mdfill.GetResponse()");
|
|---|
| 161 |
|
|---|
| 162 | MultiDish mdfill64(lambda, LMAX, vdpl64, fgnoauto);
|
|---|
| 163 | mdfill64.SetRespHisNBins(NRX,NRY);
|
|---|
| 164 | {
|
|---|
| 165 | Histo2D hpos=mdfill64.PosDist(10,10,10.*Ddish);
|
|---|
| 166 | po << PPFNameTag("posf64") << hpos;
|
|---|
| 167 | }
|
|---|
| 168 | Histo2D hrf64 = mdfill64.GetResponse();
|
|---|
| 169 | PrtTim("Apres mdfill64.GetResponse()");
|
|---|
| 170 |
|
|---|
| 171 | MultiDish mdsparse(lambda, LMAX, vdsparse, fgnoauto);
|
|---|
| 172 | mdsparse.SetThetaPhiRange(M_PI/6.,12, M_PI/6., 12);
|
|---|
| 173 | mdsparse.SetRespHisNBins(NRX,NRY);
|
|---|
| 174 | {
|
|---|
| 175 | Histo2D hpos=mdsparse.PosDist(22,22,22.*Ddish);;
|
|---|
| 176 | po << PPFNameTag("posspA") << hpos;
|
|---|
| 177 | }
|
|---|
| 178 | Histo2D hrsp = mdsparse.GetResponse();
|
|---|
| 179 | PrtTim("Apres mdsparse.GetResponse()");
|
|---|
| 180 |
|
|---|
| 181 | /*
|
|---|
| 182 | MultiDish mdsparseD7(lambda, LMAX, vdsparseD7, fgnoauto);
|
|---|
| 183 | mdsparseD7.SetThetaPhiRange(M_PI/4.,16, M_PI/4., 16);
|
|---|
| 184 | mdsparseD7.SetRespHisNBins(NRX,NRY);
|
|---|
| 185 | Histo2D hrspd7 = mdsparseD7.GetResponse();
|
|---|
| 186 | PrtTim("Apres mdsparseD7.GetResponse()");
|
|---|
| 187 | */
|
|---|
| 188 | MultiDish mdsparseB(lambda, LMAX, vdsparseB, fgnoauto);
|
|---|
| 189 | mdsparseB.SetThetaPhiRange(M_PI/6.,12, M_PI/6., 12);
|
|---|
| 190 | mdsparseB.SetRespHisNBins(NRX,NRY);
|
|---|
| 191 | {
|
|---|
| 192 | Histo2D hpos=mdsparseB.PosDist(15,15,15.*Ddish);
|
|---|
| 193 | po << PPFNameTag("posspB") << hpos;
|
|---|
| 194 | }
|
|---|
| 195 | Histo2D hrspB = mdsparseB.GetResponse();
|
|---|
| 196 | PrtTim("Apres mdsparseB.GetResponse()");
|
|---|
| 197 |
|
|---|
| 198 | MultiDish mdsparseC(lambda, LMAX, vdsparseC, fgnoauto);
|
|---|
| 199 | mdsparseC.SetThetaPhiRange(M_PI/6.,12, M_PI/6., 12);
|
|---|
| 200 | mdsparseC.SetRespHisNBins(NRX,NRY);
|
|---|
| 201 | {
|
|---|
| 202 | Histo2D hpos=mdsparseC.PosDist(20,20,20.*Ddish);
|
|---|
| 203 | po << PPFNameTag("posspC") << hpos;
|
|---|
| 204 | }
|
|---|
| 205 | Histo2D hrspC = mdsparseC.GetResponse();
|
|---|
| 206 | PrtTim("Apres mdsparseC.GetResponse()");
|
|---|
| 207 |
|
|---|
| 208 | MultiDish mdsparseBfp(lambda, LMAX, vdsparseB, fgnoauto);
|
|---|
| 209 | mdsparseBfp.SetRespHisNBins(NRX,NRY);
|
|---|
| 210 | Histo2D hrspBfp = mdsparseBfp.GetResponse();
|
|---|
| 211 | PrtTim("Apres mdsparseBfp.GetResponse()");
|
|---|
| 212 |
|
|---|
| 213 |
|
|---|
| 214 | MultiDish mcylfill(lambda, LMAX, vcylplein, fgnoauto);
|
|---|
| 215 | mcylfill.SetRespHisNBins(NRX,NRY);
|
|---|
| 216 | Histo2D hfcyl = mcylfill.GetResponse();
|
|---|
| 217 | PrtTim("Apres mcylfill.GetResponse()");
|
|---|
| 218 | MultiDish mcylfillP(lambda, LMAX, vcylplP, fgnoauto);
|
|---|
| 219 | mcylfillP.SetRespHisNBins(NRX,NRY);
|
|---|
| 220 | Histo2D hfcylP = mcylfillP.GetResponse();
|
|---|
| 221 | PrtTim("Apres mcylfillP.GetResponse()");
|
|---|
| 222 |
|
|---|
| 223 | MultiDish md3cyl(lambda, LMAX, v3cyl, fgnoauto);
|
|---|
| 224 | md3cyl.SetRespHisNBins(NRX,NRY);
|
|---|
| 225 | Histo2D h3cyl = md3cyl.GetResponse();
|
|---|
| 226 | PrtTim("Apres md3cyl.GetResponse()");
|
|---|
| 227 | MultiDish md3cylP(lambda, LMAX, v3cylP, fgnoauto);
|
|---|
| 228 | md3cylP.SetRespHisNBins(NRX,NRY);
|
|---|
| 229 | Histo2D h3cylP = md3cylP.GetResponse();
|
|---|
| 230 | PrtTim("Apres md3cylP.GetResponse()");
|
|---|
| 231 |
|
|---|
| 232 | MultiDish md2cyl(lambda, LMAX, v2cyl, fgnoauto);
|
|---|
| 233 | md2cyl.SetRespHisNBins(NRX,NRY);
|
|---|
| 234 | Histo2D h2cyl = md2cyl.GetResponse();
|
|---|
| 235 | PrtTim("Apres md2cyl.GetResponse()");
|
|---|
| 236 | MultiDish md2cylP(lambda, LMAX, v2cylP, fgnoauto);
|
|---|
| 237 | md2cylP.SetRespHisNBins(NRX,NRY);
|
|---|
| 238 | Histo2D h2cylP = md2cylP.GetResponse();
|
|---|
| 239 | PrtTim("Apres md2cylP.GetResponse()");
|
|---|
| 240 |
|
|---|
| 241 | po << PPFNameTag("mfill") << hrfill;
|
|---|
| 242 | po << PPFNameTag("mfill64") << hrf64;
|
|---|
| 243 | po << PPFNameTag("mspars") << hrsp;
|
|---|
| 244 | // po << PPFNameTag("msparsd7") << hrspd7;
|
|---|
| 245 | po << PPFNameTag("msparsB") << hrspB;
|
|---|
| 246 | po << PPFNameTag("msparsC") << hrspC;
|
|---|
| 247 | po << PPFNameTag("msparsBfp") << hrspBfp;
|
|---|
| 248 | po << PPFNameTag("mcylf") << hfcyl;
|
|---|
| 249 | po << PPFNameTag("m3cyl") << h3cyl;
|
|---|
| 250 | po << PPFNameTag("m2cyl") << h2cyl;
|
|---|
| 251 | po << PPFNameTag("mcylfP") << hfcylP;
|
|---|
| 252 | po << PPFNameTag("m3cylP") << h3cylP;
|
|---|
| 253 | po << PPFNameTag("m2cylP") << h2cylP;
|
|---|
| 254 |
|
|---|
| 255 | PrtTim("Done computing multi-dish response");
|
|---|
| 256 |
|
|---|
| 257 |
|
|---|
| 258 | Four2DRespTable mdf(hrfill, Dol);
|
|---|
| 259 | Four2DRespTable mdf64(hrf64, Dol);
|
|---|
| 260 | Four2DRespTable mds(hrsp, Dol);
|
|---|
| 261 | // Four2DRespTable mdsfp(hrspfp, Dol);
|
|---|
| 262 | // Four2DRespTable mdsd7(hrspd7, Dol);
|
|---|
| 263 | Four2DRespTable mdsB(hrspB, Dol);
|
|---|
| 264 | Four2DRespTable mdsC(hrspC, Dol);
|
|---|
| 265 | Four2DRespTable mdsBfp(hrspBfp, Dol);
|
|---|
| 266 |
|
|---|
| 267 | Four2DRespTable mcylf(hfcyl, Dol);
|
|---|
| 268 | Four2DRespTable m3cyl(h3cyl, Dol);
|
|---|
| 269 | Four2DRespTable m2cyl(h2cyl, Dol);
|
|---|
| 270 | Four2DRespTable mcylfP(hfcylP, Dol);
|
|---|
| 271 | Four2DRespTable m3cylP(h3cylP, Dol);
|
|---|
| 272 | Four2DRespTable m2cylP(h2cylP, Dol);
|
|---|
| 273 |
|
|---|
| 274 |
|
|---|
| 275 | cout << " 3.a/ Instanciating object type Four3DPk " << endl;
|
|---|
| 276 | RandomGenerator rg;
|
|---|
| 277 | Four3DPk m3d(rg);
|
|---|
| 278 | m3d.SetCellSize(2.*DeuxPI, 2.*DeuxPI, 2.*DeuxPI);
|
|---|
| 279 | cout << " 3.b/ call ComputeFourierAmp() NGEN=" << NMAX << endl;
|
|---|
| 280 | HProf hrpk;
|
|---|
| 281 | for(int igen=0; igen<NMAX; igen++) {
|
|---|
| 282 | m3d.ComputeFourierAmp(spec);
|
|---|
| 283 | if (igen==0) hrpk = m3d.ComputePk();
|
|---|
| 284 | else m3d.ComputePkCumul(hrpk);
|
|---|
| 285 | }
|
|---|
| 286 | PrtTim("md.ComputeFourierAmp() done");
|
|---|
| 287 | po << PPFNameTag("recPk") << hrpk;
|
|---|
| 288 |
|
|---|
| 289 | cout << " 4/ Computing Noise P(k) using PkNoiseCalculator ..." << endl;
|
|---|
| 290 | #define NCONFIG 14
|
|---|
| 291 | Four2DResponse* f2rep[NCONFIG]={&dish, &dish2, &mdf, &mdf64, &mds, &mdsB, &mdsC, &mdsBfp,
|
|---|
| 292 | &mcylf, &mcylfP, &m3cyl, &m3cylP, &m2cyl, &m2cylP};
|
|---|
| 293 | const char* tits[NCONFIG]={"Dish100m", "Dish200m","F20x20Dish5m","F8x8Dish5m",
|
|---|
| 294 | "S68Dish5m","S72Dish5m","S129CDish5m","S72BDish5mFP",
|
|---|
| 295 | "F8Cyl","F8CylP","F3Cyl","F3CylP","BiCyl","BiCylP"};
|
|---|
| 296 | const char* tags[NCONFIG]={"noiseD", "noiseD2","noisemdf","noisemdf64","noisemds",
|
|---|
| 297 | "noisemdsB","noisemdsC","noisemdsBfp",
|
|---|
| 298 | "noisefcyl","noisefcylP","noise3cyl","noise3cylP", "noise2cyl","noise2cylP"};
|
|---|
| 299 | vector<int> nbdishes;
|
|---|
| 300 | nbdishes.push_back(1);
|
|---|
| 301 | nbdishes.push_back(1);
|
|---|
| 302 | nbdishes.push_back(vdplein.size());
|
|---|
| 303 | nbdishes.push_back(vdpl64.size());
|
|---|
| 304 | nbdishes.push_back(vdsparse.size());
|
|---|
| 305 | // nbdishes.push_back(vdsparse.size());
|
|---|
| 306 | nbdishes.push_back(vdsparseB.size());
|
|---|
| 307 | nbdishes.push_back(vdsparseC.size());
|
|---|
| 308 | nbdishes.push_back(vdsparseB.size());
|
|---|
| 309 | nbdishes.push_back(vcylplein.size());
|
|---|
| 310 | nbdishes.push_back(vcylplP.size());
|
|---|
| 311 | nbdishes.push_back(v3cyl.size());
|
|---|
| 312 | nbdishes.push_back(v3cylP.size());
|
|---|
| 313 | nbdishes.push_back(v2cyl.size());
|
|---|
| 314 | nbdishes.push_back(v2cylP.size());
|
|---|
| 315 |
|
|---|
| 316 | for(int lc=0; lc<NCONFIG; lc++) {
|
|---|
| 317 | PkNoiseCalculator pkn(m3d, *(f2rep[lc]), SCut/(double)nbdishes[lc], NMAX, tits[lc]);
|
|---|
| 318 | HProf hpn = pkn.Compute();
|
|---|
| 319 | po << PPFNameTag(tags[lc]) << hpn;
|
|---|
| 320 | }
|
|---|
| 321 | rc = 0;
|
|---|
| 322 | } // End of try bloc
|
|---|
| 323 | catch (PThrowable & exc) { // catching SOPHYA exceptions
|
|---|
| 324 | cerr << " pknoise.cc: Catched Exception (PThrowable)" << (string)typeid(exc).name()
|
|---|
| 325 | << "\n...exc.Msg= " << exc.Msg() << endl;
|
|---|
| 326 | rc = 99;
|
|---|
| 327 | }
|
|---|
| 328 | catch (std::exception & e) { // catching standard C++ exceptions
|
|---|
| 329 | cerr << " pknoise.cc: Catched std::exception " << " - what()= " << e.what() << endl;
|
|---|
| 330 | rc = 98;
|
|---|
| 331 | }
|
|---|
| 332 | catch (...) { // catching other exceptions
|
|---|
| 333 | cerr << " pknoise.cc: some other exception (...) was caught ! " << endl;
|
|---|
| 334 | rc = 97;
|
|---|
| 335 | }
|
|---|
| 336 | PrtTim("End-pknoise");
|
|---|
| 337 | cout << " ==== End of pknoise.cc program Rc= " << rc << endl;
|
|---|
| 338 | return rc;
|
|---|
| 339 | }
|
|---|
| 340 |
|
|---|
| 341 |
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
