| 1 | /* ------------------------ Projet BAORadio --------------------
|
|---|
| 2 | Calcul de la reponse 2D (plan (u,v) d'un interferometre
|
|---|
| 3 | R. Ansari , C. Magneville - Juin-Dec 2010
|
|---|
| 4 |
|
|---|
| 5 | Usage: repicon [-parname value] configId OutPPFName
|
|---|
| 6 | --------------------------------------------------------------- */
|
|---|
| 7 |
|
|---|
| 8 | #include "machdefs.h"
|
|---|
| 9 | #include "sopnamsp.h"
|
|---|
| 10 | #include <iostream>
|
|---|
| 11 | #include <string>
|
|---|
| 12 | #include <math.h>
|
|---|
| 13 |
|
|---|
| 14 | #include <typeinfo>
|
|---|
| 15 |
|
|---|
| 16 | #include "mdish.h"
|
|---|
| 17 | #include "specpk.h"
|
|---|
| 18 | #include "interfconfigs.h"
|
|---|
| 19 | #include "radutil.h"
|
|---|
| 20 |
|
|---|
| 21 | #include "ntuple.h"
|
|---|
| 22 | #include "histinit.h"
|
|---|
| 23 | // #include "fiosinit.h"
|
|---|
| 24 | // #include "fitsioserver.h"
|
|---|
| 25 |
|
|---|
| 26 | #include "randr48.h"
|
|---|
| 27 |
|
|---|
| 28 | #include "timing.h"
|
|---|
| 29 | #include "ctimer.h"
|
|---|
| 30 |
|
|---|
| 31 | typedef DR48RandGen RandomGenerator ;
|
|---|
| 32 |
|
|---|
| 33 | // pour sauver la reponse mdresp et la config des dishes dans un fichier PPF
|
|---|
| 34 | void SaveDTVecDishesH2Resp(POutPersist& po, vector<Dish>& vdishes, Four2DRespTable& mdresp);
|
|---|
| 35 |
|
|---|
| 36 | // ---------------------------------------------------------------------
|
|---|
| 37 | // main program for computing interferometer response in (u,v) plane
|
|---|
| 38 | // R. Ansari - Avril-Juin 2010
|
|---|
| 39 | // ---------------------------------------------------------------------
|
|---|
| 40 |
|
|---|
| 41 | void Usage()
|
|---|
| 42 | {
|
|---|
| 43 | cout << " Usage: repicon [-parname Value] configId OutPPFName \n"
|
|---|
| 44 | << " configIds: f4x4,f8x8,f11x11,f20x20, confA,confB,confC,confD, hex12,cross11, \n"
|
|---|
| 45 | << " f4cyl,f8cyl,f4cylp,f4cylp, s4x4,s11x11,s20x20 nan12,nan24,nan36,nan40,nan128 \n"
|
|---|
| 46 | << " f4x4 , f8x8 , f11x11 , f20x20 Filled array of nxn dishes \n"
|
|---|
| 47 | << " s4x4 , s11x11 , s20x20 sparse/semi-Filled regular array of nxn dishes \n"
|
|---|
| 48 | << " confA , confB, confC, confD : semi-filled array of dishes \n"
|
|---|
| 49 | << " hex12,cross11 : ASKAP like double hexagonal (12xD=12m), cross config (11xD=12m) \n"
|
|---|
| 50 | << " nan12,nan24,nan36,nan40,nan128: 3,4,6,5,8 cylinder like configurations with 12,24,36,40,128 dishes\n"
|
|---|
| 51 | << " f4cyl, f8cyl , f4cylp, f8cylp : filled array of non perfect/perfect of n cylinders \n"
|
|---|
| 52 | << " f4cylw : filled array of 4 perfect of wide cylinders \n"
|
|---|
| 53 | << " pit2cyl,pit2cylw : 2 cylinders 15mx7m, 20mx10 for z=0.3 at Pittsburg (perfect) \n"
|
|---|
| 54 | << " pit4cyl : 4 cylinders 32mx8m for z=0.3-0.7 at Pittsburg (perfect) \n"
|
|---|
| 55 | << " [ -parname value] : -renmax -z -prt -D -lmax -eta -autocor \n"
|
|---|
| 56 | << " -renmax MaxValue (default : Do NOT renormalize 2D response value \n"
|
|---|
| 57 | << " -z redshift (default=0.7) --> determines Lambda \n"
|
|---|
| 58 | << " -D DishDiameter (default=5 m) \n"
|
|---|
| 59 | << " -lsep Lsep (default=5 m) dish separation for sparse/semi-filled array \n"
|
|---|
| 60 | << " -eta fill_factor (default=0.90) \n"
|
|---|
| 61 | << " -lmax array extension (default=100 m ) for response calculation kmax \n"
|
|---|
| 62 | << " -repnxy Nx,Ny (default=200,200) ResponseTable binning \n"
|
|---|
| 63 | << " -beamnxy Nx,Ny (default=200,200) Beam sampling \n"
|
|---|
| 64 | << " -rot ThetaMaxDeg,NTheta,PhiMaxDeg,NPhi (default NO-rotate/pointing -> 23,10,30,15 ) \n"
|
|---|
| 65 | << " -autocor : keep antenna auto-correlation signals (default NO) \n"
|
|---|
| 66 | << " -prt PrtLev,PrtModulo (default=0,10) \n"
|
|---|
| 67 | << endl;
|
|---|
| 68 | return;
|
|---|
| 69 | }
|
|---|
| 70 |
|
|---|
| 71 | // ------------------ MAIN PROGRAM ------------------------------
|
|---|
| 72 | int main(int narg, const char* arg[])
|
|---|
| 73 | {
|
|---|
| 74 | if (((narg>1)&&(strcmp(arg[1],"-h")==0))||(narg<3)) {
|
|---|
| 75 | Usage();
|
|---|
| 76 | return 1;
|
|---|
| 77 | }
|
|---|
| 78 | // make sure SOPHYA modules are initialized
|
|---|
| 79 | SophyaInit();
|
|---|
| 80 | // FitsIOServerInit();
|
|---|
| 81 | InitTim();
|
|---|
| 82 | //--- decoding command line arguments
|
|---|
| 83 | string config = "f8x8" ;
|
|---|
| 84 | string outfile = "repicon.ppf";
|
|---|
| 85 | bool fgrenorm=false;
|
|---|
| 86 | double rmax=1.;
|
|---|
| 87 | double z_Redshift=0.7 ; // 21 cm at z=0.7 -> 0.357 m
|
|---|
| 88 | int prtlev=0;
|
|---|
| 89 | int prtmod=10;
|
|---|
| 90 |
|
|---|
| 91 | double Ddish=5.;
|
|---|
| 92 | double Lsep=5.;
|
|---|
| 93 | double Eta=0.9;
|
|---|
| 94 | bool fgDfixed=false;
|
|---|
| 95 | double LMAX = 100.; // taille de la zone, pour calcul kmax
|
|---|
| 96 | bool fgLMAXfixed=false;
|
|---|
| 97 | double thetamxdeg=23.; // 23 degres : angle d'inclinaison de l'orbite terrestre
|
|---|
| 98 | double phimxdeg=30.;
|
|---|
| 99 | int nteta=10;
|
|---|
| 100 | int nphi=15;
|
|---|
| 101 | bool fgpoint=false;
|
|---|
| 102 | int NRX=200;
|
|---|
| 103 | int NRY=200;
|
|---|
| 104 | int NBX=200;
|
|---|
| 105 | int NBY=200;
|
|---|
| 106 | bool fgnoauto=true; // do not keep antenna autocorrelation signal
|
|---|
| 107 |
|
|---|
| 108 | int ka=1;
|
|---|
| 109 | while (ka<(narg-1)) {
|
|---|
| 110 | if (strcmp(arg[ka],"-renmax")==0) {
|
|---|
| 111 | rmax=atof(arg[ka+1]); fgrenorm=true; ka+=2;
|
|---|
| 112 | }
|
|---|
| 113 | else if (strcmp(arg[ka],"-z")==0) {
|
|---|
| 114 | z_Redshift=atof(arg[ka+1]); ka+=2;
|
|---|
| 115 | }
|
|---|
| 116 | else if (strcmp(arg[ka],"-D")==0) {
|
|---|
| 117 | Ddish=atof(arg[ka+1]); fgDfixed=true; ka+=2;
|
|---|
| 118 | }
|
|---|
| 119 | else if (strcmp(arg[ka],"-lsep")==0) {
|
|---|
| 120 | Lsep=atof(arg[ka+1]); ka+=2;
|
|---|
| 121 | }
|
|---|
| 122 | else if (strcmp(arg[ka],"-eta")==0) {
|
|---|
| 123 | Eta=atof(arg[ka+1]); ka+=2;
|
|---|
| 124 | }
|
|---|
| 125 | else if (strcmp(arg[ka],"-lmax")==0) {
|
|---|
| 126 | LMAX=atof(arg[ka+1]); fgLMAXfixed=true; ka+=2;
|
|---|
| 127 | }
|
|---|
| 128 | else if (strcmp(arg[ka],"-repnxy")==0) {
|
|---|
| 129 | sscanf(arg[ka+1],"%d,%d",&NRX,&NRY); ka+=2;
|
|---|
| 130 | }
|
|---|
| 131 | else if (strcmp(arg[ka],"-beamnxy")==0) {
|
|---|
| 132 | sscanf(arg[ka+1],"%d,%d",&NBX,&NBY); ka+=2;
|
|---|
| 133 | }
|
|---|
| 134 | else if (strcmp(arg[ka],"-rot")==0) {
|
|---|
| 135 | sscanf(arg[ka+1],"%lg,%d,%lg,%d",&thetamxdeg,&nteta,&phimxdeg,&nphi); fgpoint=true; ka+=2;
|
|---|
| 136 | }
|
|---|
| 137 | else if (strcmp(arg[ka],"-autocor")==0) {
|
|---|
| 138 | fgnoauto=false; ka+=1;
|
|---|
| 139 | }
|
|---|
| 140 | else if (strcmp(arg[ka],"-prt")==0) {
|
|---|
| 141 | sscanf(arg[ka+1],"%d,%d",&prtlev,&prtmod); ka+=2;
|
|---|
| 142 | }
|
|---|
| 143 | else break;
|
|---|
| 144 | }
|
|---|
| 145 |
|
|---|
| 146 | if ((ka+1)>=narg) {
|
|---|
| 147 | cout << " repicon / Argument error " << endl;
|
|---|
| 148 | Usage();
|
|---|
| 149 | return 2;
|
|---|
| 150 | }
|
|---|
| 151 |
|
|---|
| 152 | config = arg[ka];
|
|---|
| 153 | outfile = arg[ka+1];
|
|---|
| 154 | //-- end command line arguments
|
|---|
| 155 |
|
|---|
| 156 | int rc = 1;
|
|---|
| 157 | try { // exception handling try bloc at top level
|
|---|
| 158 | cout << " ==== repicon.cc program , (u,v) plane response ==== " << endl;
|
|---|
| 159 |
|
|---|
| 160 | H21Conversions conv;
|
|---|
| 161 | double LAMBDA=0.357 ; // 21 cm at z=0.7
|
|---|
| 162 | conv.setRedshift(z_Redshift);
|
|---|
| 163 | LAMBDA = conv.getLambda();
|
|---|
| 164 |
|
|---|
| 165 | double cylW=12.; // Largeur des cylindres
|
|---|
| 166 | double cylRL=2*LAMBDA; // Longeur des elements de reception le long du cylindre
|
|---|
| 167 | double etaW=0.9; // Efficacite de couverture en largeur
|
|---|
| 168 | double etaRL=0.8; // Efficacite de couverture le long du cylindre
|
|---|
| 169 |
|
|---|
| 170 | vector<Dish> vdishes;
|
|---|
| 171 | cout << " repicon[1] : creating dish vector/ MultiDish for configuration: " << config
|
|---|
| 172 | << " Ddish=" << Ddish << " m. Eta=" << Eta << endl;
|
|---|
| 173 |
|
|---|
| 174 | if (config=="f4x4") {
|
|---|
| 175 | vdishes=CreateFilledSqConfig(4,Ddish, Eta);
|
|---|
| 176 | }
|
|---|
| 177 | else if (config=="f8x8") {
|
|---|
| 178 | vdishes=CreateFilledSqConfig(8,Ddish, Eta);
|
|---|
| 179 | }
|
|---|
| 180 | else if (config=="f11x11") {
|
|---|
| 181 | vdishes=CreateFilledSqConfig(11,Ddish, Eta);
|
|---|
| 182 | }
|
|---|
| 183 | else if (config=="f20x20") {
|
|---|
| 184 | vdishes=CreateFilledSqConfig(20,Ddish, Eta);
|
|---|
| 185 | }
|
|---|
| 186 |
|
|---|
| 187 | else if (config=="s4x4") {
|
|---|
| 188 | vdishes=CreateSparseSqConfig(4, Lsep, Ddish, Eta);
|
|---|
| 189 | }
|
|---|
| 190 | else if (config=="s11x11") {
|
|---|
| 191 | vdishes=CreateSparseSqConfig(11, Lsep, Ddish, Eta);
|
|---|
| 192 | }
|
|---|
| 193 | else if (config=="s20x20") {
|
|---|
| 194 | vdishes=CreateSparseSqConfig(20, Lsep, Ddish, Eta);
|
|---|
| 195 | }
|
|---|
| 196 |
|
|---|
| 197 | else if (config=="f4cyl") {
|
|---|
| 198 | cylW=12.; cylRL=2*LAMBDA;
|
|---|
| 199 | vdishes=CreateFilledCylConfig(4, 100, cylW, cylRL, etaW, etaRL, true);
|
|---|
| 200 | }
|
|---|
| 201 | else if (config=="f4cylp") {
|
|---|
| 202 | cylW=12.; cylRL=2*LAMBDA;
|
|---|
| 203 | vdishes=CreateFilledCylConfig(4, 100, cylW, cylRL, etaW, etaRL, false);
|
|---|
| 204 | }
|
|---|
| 205 | else if (config=="f4cylw") {
|
|---|
| 206 | cylW=25.; cylRL=3*LAMBDA;
|
|---|
| 207 | vdishes=CreateFilledCylConfig(4, 100, cylW, cylRL, etaW, etaRL, false);
|
|---|
| 208 | }
|
|---|
| 209 | else if (config=="pit2cyl") {
|
|---|
| 210 | etaW=0.9; etaRL=0.9;
|
|---|
| 211 | cylW=7.; cylRL=2*LAMBDA;
|
|---|
| 212 | vdishes=CreateFilledCylConfig(2, 32, cylW, cylRL, etaW, etaRL, false);
|
|---|
| 213 | }
|
|---|
| 214 | else if (config=="pit2cylw") {
|
|---|
| 215 | etaW=0.9; etaRL=0.9;
|
|---|
| 216 | cylW=9.; cylRL=2*LAMBDA;
|
|---|
| 217 | vdishes=CreateFilledCylConfig(2, 40, cylW, cylRL, etaW, etaRL, false);
|
|---|
| 218 | }
|
|---|
| 219 | else if (config=="pit4cyl") {
|
|---|
| 220 | etaW=0.9; etaRL=0.9;
|
|---|
| 221 | cylW=8.; cylRL=2*LAMBDA;
|
|---|
| 222 | vdishes=CreateFilledCylConfig(4, 64, cylW, cylRL, etaW, etaRL, false);
|
|---|
| 223 | }
|
|---|
| 224 | else if (config=="f8cyl") {
|
|---|
| 225 | cylW=12.; cylRL=2*LAMBDA;
|
|---|
| 226 | vdishes=CreateFilledCylConfig(8, 120, cylW, cylRL, etaW, etaRL, true);
|
|---|
| 227 | }
|
|---|
| 228 | else if (config=="f8cylp") {
|
|---|
| 229 | cylW=12.; cylRL=2*LAMBDA;
|
|---|
| 230 | vdishes=CreateFilledCylConfig(8, 120, cylW, cylRL, etaW, etaRL, false);
|
|---|
| 231 | }
|
|---|
| 232 | else if (config=="confA") {
|
|---|
| 233 | vdishes=CreateConfigA(Ddish, Eta);
|
|---|
| 234 | }
|
|---|
| 235 | else if (config=="confB") {
|
|---|
| 236 | vdishes=CreateConfigB(Ddish, Eta);
|
|---|
| 237 | }
|
|---|
| 238 | else if (config=="confC") {
|
|---|
| 239 | vdishes=CreateConfigC(Ddish, Eta);
|
|---|
| 240 | }
|
|---|
| 241 | else if (config=="confD") {
|
|---|
| 242 | vdishes=CreateConfigD(Ddish, Eta);
|
|---|
| 243 | }
|
|---|
| 244 | else if (config=="cross11") {
|
|---|
| 245 | if (!fgDfixed) Ddish = 12.;
|
|---|
| 246 | double base=20.;
|
|---|
| 247 | Eta=0.95;
|
|---|
| 248 | if (!fgLMAXfixed) LMAX = 250.;
|
|---|
| 249 | vdishes=CreateCrossConfig(Ddish,base,Eta);
|
|---|
| 250 | }
|
|---|
| 251 | else if (config=="hex12") {
|
|---|
| 252 | if (!fgDfixed) Ddish = 12.;
|
|---|
| 253 | Eta=0.95;
|
|---|
| 254 | if (!fgLMAXfixed) LMAX = 350.;
|
|---|
| 255 | vdishes=CreateDoubleHexagonConfig(Ddish);
|
|---|
| 256 | }
|
|---|
| 257 | else if (config=="nan12") {
|
|---|
| 258 | vdishes=CreateConfigNancay12(Ddish, Eta);
|
|---|
| 259 | }
|
|---|
| 260 | else if (config=="nan24") {
|
|---|
| 261 | vdishes=CreateConfigNancay24(Ddish, Eta);
|
|---|
| 262 | }
|
|---|
| 263 | else if (config=="nan36") {
|
|---|
| 264 | vdishes=CreateConfigNancay36(Ddish, Eta);
|
|---|
| 265 | }
|
|---|
| 266 | else if (config=="nan40") {
|
|---|
| 267 | vdishes=CreateConfigNancay40(Ddish, Eta);
|
|---|
| 268 | }
|
|---|
| 269 | else if (config=="nan128") {
|
|---|
| 270 | vdishes=CreateConfigNancay128(Ddish, Eta);
|
|---|
| 271 | }
|
|---|
| 272 | else {
|
|---|
| 273 | cout << " NON valid configuration option -> exit" << endl;
|
|---|
| 274 | return 99;
|
|---|
| 275 | }
|
|---|
| 276 |
|
|---|
| 277 | double Dol = LMAX/LAMBDA;
|
|---|
| 278 |
|
|---|
| 279 | cout << " repicon[1] : Lambda=" << LAMBDA << " LMAX= " << LMAX << " NDishes=" << vdishes.size()
|
|---|
| 280 | << " D-Dish=" << Ddish << " m. " << ((fgnoauto)?" NO-":"With-") << "AutoCorrelation" << endl;
|
|---|
| 281 |
|
|---|
| 282 | MultiDish mdish(LAMBDA, LMAX, vdishes, fgnoauto);
|
|---|
| 283 | mdish.SetPrtLevel(prtlev,prtmod);
|
|---|
| 284 | mdish.SetRespHisNBins(NRX,NRY);
|
|---|
| 285 | mdish.SetBeamNSamples(NBX,NBY);
|
|---|
| 286 | if (fgpoint) {
|
|---|
| 287 | cout << " repicon[1.b] : activating pointing , ThetaMaxDeg=" << thetamxdeg << " NTheta=" << nteta
|
|---|
| 288 | << " PhiMaxDeg= " << phimxdeg << " NPhi=" << nphi << endl;
|
|---|
| 289 | mdish.SetThetaPhiRange(Angle(thetamxdeg,Angle::Degree).ToRadian(),nteta,Angle(phimxdeg,Angle::Degree).ToRadian(), nphi);
|
|---|
| 290 | }
|
|---|
| 291 | cout << " repicon[2] : calling mdish.GetResponse() ..."<< endl;
|
|---|
| 292 |
|
|---|
| 293 | Histo2D hrep = mdish.GetResponse();
|
|---|
| 294 | PrtTim("Apres mdish.GetResponse()");
|
|---|
| 295 | HProf h1dnoise = mdish.GetProjNoiseLevel();
|
|---|
| 296 | HProf h1drep = mdish.GetProjResponse();
|
|---|
| 297 |
|
|---|
| 298 | Four2DRespTable mdresp(hrep, Dol, LAMBDA);
|
|---|
| 299 | if (fgrenorm) {
|
|---|
| 300 | cout << " repicon[2.b] call to mdresp.renormalize(" << rmax << ")";
|
|---|
| 301 | double omax=mdresp.renormalize(rmax);
|
|---|
| 302 | cout << " Old Max=" << omax << endl;
|
|---|
| 303 | }
|
|---|
| 304 | cout << " repicon[3] : saving Four2DRespTable for config " << config << " to PPF file " << outfile << endl;
|
|---|
| 305 | mdresp.writeToPPF(outfile);
|
|---|
| 306 |
|
|---|
| 307 | string outfile2 = "hdt_"+outfile;
|
|---|
| 308 | POutPersist po2(outfile2);
|
|---|
| 309 | cout << " repicon[4] : saving H1D/H2D-response, multidish config to PPF file " << outfile2 << endl;
|
|---|
| 310 | po2 << PPFNameTag("h1dnoise") << h1dnoise;
|
|---|
| 311 | po2 << PPFNameTag("h1drep") << h1drep;
|
|---|
| 312 | SaveDTVecDishesH2Resp(po2, vdishes, mdresp);
|
|---|
| 313 |
|
|---|
| 314 | rc = 0;
|
|---|
| 315 | } // End of try bloc
|
|---|
| 316 | catch (PThrowable & exc) { // catching SOPHYA exceptions
|
|---|
| 317 | cerr << " repicon.cc: Catched Exception (PThrowable)" << (string)typeid(exc).name()
|
|---|
| 318 | << "\n...exc.Msg= " << exc.Msg() << endl;
|
|---|
| 319 | rc = 99;
|
|---|
| 320 | }
|
|---|
| 321 | catch (std::exception & e) { // catching standard C++ exceptions
|
|---|
| 322 | cerr << " repicon.cc: Catched std::exception " << " - what()= " << e.what() << endl;
|
|---|
| 323 | rc = 98;
|
|---|
| 324 | }
|
|---|
| 325 | catch (...) { // catching other exceptions
|
|---|
| 326 | cerr << " repicon.cc: some other exception (...) was caught ! " << endl;
|
|---|
| 327 | rc = 97;
|
|---|
| 328 | }
|
|---|
| 329 | PrtTim("End-repicon");
|
|---|
| 330 | cout << " ==== End of repicon.cc program Rc= " << rc << endl;
|
|---|
| 331 | return rc;
|
|---|
| 332 | }
|
|---|
| 333 |
|
|---|
| 334 | /*-- Nouvelle-Fonction --*/
|
|---|
| 335 | void SaveDTVecDishesH2Resp(POutPersist& po, vector<Dish>& vdishes, Four2DRespTable& mdresp)
|
|---|
| 336 | {
|
|---|
| 337 | const char* names[5]={"did","posx","posy","diam","diamy"};
|
|---|
| 338 | NTuple ntvd(5,names,64,false);
|
|---|
| 339 | r_4 xnt[10];
|
|---|
| 340 | for(size_t i=0; i<vdishes.size(); i++) {
|
|---|
| 341 | xnt[0]=vdishes[i].ReflectorId();
|
|---|
| 342 | xnt[1]=vdishes[i].X;
|
|---|
| 343 | xnt[2]=vdishes[i].Y;
|
|---|
| 344 | if (vdishes[i].isCircular()) {
|
|---|
| 345 | xnt[3]=vdishes[i].Diameter();
|
|---|
| 346 | xnt[4]=0.;
|
|---|
| 347 | }
|
|---|
| 348 | else {
|
|---|
| 349 | xnt[3]=vdishes[i].DiameterX();
|
|---|
| 350 | xnt[4]=vdishes[i].DiameterY();
|
|---|
| 351 | }
|
|---|
| 352 | ntvd.Fill(xnt);
|
|---|
| 353 | }
|
|---|
| 354 | Histo2D h2rep=mdresp.GetResponse();
|
|---|
| 355 | po << PPFNameTag("mdish") << ntvd;
|
|---|
| 356 | po << PPFNameTag("h2rep") << h2rep;
|
|---|
| 357 | return;
|
|---|
| 358 | }
|
|---|
| 359 |
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
