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