| 1 | /* ------------------------ Projet BAORadio --------------------
|
|---|
| 2 | Programme de convolution avec le lobe d'un cube 3D (angles,freq)
|
|---|
| 3 | R. Ansari , C. Magneville - Juin 2010
|
|---|
| 4 |
|
|---|
| 5 | Usage: applobe [-g -t -fib] Diameter/Four2DRespTableFile In3DPPFName Out3DPPFName
|
|---|
| 6 | [TargetBeamDiam] [ResmapleFactor=0.5,0.333...]
|
|---|
| 7 | o -t : triangle beam shape in k-space
|
|---|
| 8 | o -g : gaussian beam shape in k-space
|
|---|
| 9 | o -fib: application d'un lobe fixe (independant de la frequence)
|
|---|
| 10 | o Diameter/Four2DRespTableFile : diametre dish ou nom de fichier PPF reponse 2D
|
|---|
| 11 | o In3DPPFName Out3DPPFName : Noms des fichiers (TArray<r_4> 3D) entree-sortie
|
|---|
| 12 | o TargetBeamDiam : Correction à un beam fiduciel, independant de la frequence
|
|---|
| 13 | avec specification de la valeur du diametre pour la frequence la plus basse
|
|---|
| 14 | o ResmapleFactor : reechantillonnage du cube (2 -> surechantillonnage, 0.5 : 1 sur 2)
|
|---|
| 15 | --------------------------------------------------------------- */
|
|---|
| 16 |
|
|---|
| 17 | #include "sopnamsp.h"
|
|---|
| 18 | #include "machdefs.h"
|
|---|
| 19 | #include <math.h>
|
|---|
| 20 | #include <iostream>
|
|---|
| 21 | #include <typeinfo>
|
|---|
| 22 |
|
|---|
| 23 | #include "array.h"
|
|---|
| 24 | #include "histats.h"
|
|---|
| 25 |
|
|---|
| 26 | #include "swfitsdtable.h"
|
|---|
| 27 | #include "fitshdtable.h"
|
|---|
| 28 |
|
|---|
| 29 | #include "randr48.h"
|
|---|
| 30 | #include "vector3d.h"
|
|---|
| 31 |
|
|---|
| 32 | // #include "xastropack.h" -- Pour faire les conversions de coordonnees celestes
|
|---|
| 33 |
|
|---|
| 34 | #include "radutil.h"
|
|---|
| 35 | #include "lobe.h"
|
|---|
| 36 |
|
|---|
| 37 | // Pour l'initialisation des modules
|
|---|
| 38 | #include "tarrinit.h"
|
|---|
| 39 | #include "histinit.h"
|
|---|
| 40 | #include "fiosinit.h"
|
|---|
| 41 |
|
|---|
| 42 | #include "timing.h"
|
|---|
| 43 | #include "ctimer.h"
|
|---|
| 44 |
|
|---|
| 45 | #include "cubedef.h"
|
|---|
| 46 |
|
|---|
| 47 | //----------------------------------------------------------------------------
|
|---|
| 48 | //----------------------------------------------------------------------------
|
|---|
| 49 | int main(int narg, char* arg[])
|
|---|
| 50 | {
|
|---|
| 51 | // Sophya modules initialization
|
|---|
| 52 | TArrayInitiator _inia;
|
|---|
| 53 | HiStatsInitiator _inih;
|
|---|
| 54 | FitsIOServerInitiator _inif;
|
|---|
| 55 | //------- AU LIEU DE ------> SophyaInit();
|
|---|
| 56 |
|
|---|
| 57 | InitTim(); // Initializing the CPU timer
|
|---|
| 58 |
|
|---|
| 59 | if ((narg < 3)||(strcmp(arg[1],"-h")==0)) {
|
|---|
| 60 | cout << "Usage: applobe [-t -g -fib -mxr val] Diameter/Four2DRespTableFile In3DPPFName Out3DPPFName \n"
|
|---|
| 61 | << " [TargetBeamDiam] [ResmapleFactor=0.5,0.333...] \n" << endl;
|
|---|
| 62 | if ((narg>1)&&(strcmp(arg[1],"-h")==0)) {
|
|---|
| 63 | cout << " -t -g : Triangular / gaussian beam shape (def=gaussian) \n"
|
|---|
| 64 | << " -fib : Application of a fixed (freq.independent) lobe dish-triangle or gaussian \n"
|
|---|
| 65 | << " -mxr val: Max beam correction factor (default=10.) \n"
|
|---|
| 66 | << " Diameter/Four2DRespTableFile : dish diameter or 2D response PPF file name\n"
|
|---|
| 67 | << " In3DPPFName Out3DPPFName: Input/Output PPF file names (TArray<r_4> 3D) \n"
|
|---|
| 68 | << " TargetBeamDiam: Corrected beam target diameter (D/Lambda at lowest frequency) \n"
|
|---|
| 69 | << " DoL = 100 --> beam ~ 35 arcmin (D=30m @ z~0.5 Lambda~30cm)) \n"
|
|---|
| 70 | << " ResmapleFactor : Resampling (2-> oversampling, 0.5 : 1/2 undersampling) \n" << endl;
|
|---|
| 71 | }
|
|---|
| 72 | return 1;
|
|---|
| 73 | }
|
|---|
| 74 |
|
|---|
| 75 | Timer tm("applobe");
|
|---|
| 76 |
|
|---|
| 77 | // decodage arguments
|
|---|
| 78 | bool fggaussian=true; // true -> gaussian beam
|
|---|
| 79 | bool fixedbeam=false; // true -> apply freq. independent beam
|
|---|
| 80 | double maxratio=10.; // valeur max du rapport des lobes lors de la correction de lobe
|
|---|
| 81 |
|
|---|
| 82 | // decodage argument optionnel
|
|---|
| 83 | bool fgoptarg=true;
|
|---|
| 84 | while (fgoptarg) {
|
|---|
| 85 | string fbo = arg[1];
|
|---|
| 86 | if (fbo=="-t") { fggaussian=false; arg++; narg--; }
|
|---|
| 87 | else if (fbo=="-g") { fggaussian=true; arg++; narg--; }
|
|---|
| 88 | else if (fbo=="-fib") { fixedbeam=true; arg++; narg--; }
|
|---|
| 89 | else if (fbo=="-mxr") { arg++; maxratio=atof(arg[1]); arg++; narg-=2; }
|
|---|
| 90 | else fgoptarg=false;
|
|---|
| 91 | }
|
|---|
| 92 | if (narg < 4) {
|
|---|
| 93 | cout << " applobe/error arguments , applobe -h for help " << endl;
|
|---|
| 94 | return 2;
|
|---|
| 95 | }
|
|---|
| 96 |
|
|---|
| 97 | bool fgresptbl=true;
|
|---|
| 98 | double DIAMETRE=100.;
|
|---|
| 99 | string resptblname;
|
|---|
| 100 | if (isdigit(*arg[1])) {
|
|---|
| 101 | fgresptbl=false;
|
|---|
| 102 | DIAMETRE=atof(arg[1]);
|
|---|
| 103 | }
|
|---|
| 104 | else resptblname=arg[1];
|
|---|
| 105 |
|
|---|
| 106 | string inname = arg[2];
|
|---|
| 107 | string outname = arg[3];
|
|---|
| 108 | bool fgcorrbeam=false;
|
|---|
| 109 | double tbeamDiam=50.;
|
|---|
| 110 | if (narg>4) {
|
|---|
| 111 | tbeamDiam=atof(arg[4]);
|
|---|
| 112 | if (tbeamDiam>1.) fgcorrbeam=true;
|
|---|
| 113 | }
|
|---|
| 114 | bool fgresample=false;
|
|---|
| 115 | double fresamp=1.;
|
|---|
| 116 | if (narg>5) {
|
|---|
| 117 | fresamp=atof(arg[5]);
|
|---|
| 118 | if ((fabs(fresamp)>1.e-2)&&(fabs(fresamp-1.)>1.e-2)) fgresample=true;
|
|---|
| 119 | }
|
|---|
| 120 |
|
|---|
| 121 | int rc = 91;
|
|---|
| 122 |
|
|---|
| 123 | cout << " ====== applobe : Input skycube name= " << inname << " OutName=" << outname << endl;
|
|---|
| 124 | cout << ((fggaussian)?" Gaussian ":" Triangular") << ((fixedbeam)?" Fixed (freq.independent)":"") << " beams" << endl;
|
|---|
| 125 | bool fginmap=true;
|
|---|
| 126 | try {
|
|---|
| 127 | TArray<r_4> incube;
|
|---|
| 128 | cout << "applobe[1]: reading input 3D map (cube) from file " << inname << endl;
|
|---|
| 129 | {
|
|---|
| 130 | PInPersist pin(inname);
|
|---|
| 131 | pin >> incube;
|
|---|
| 132 | }
|
|---|
| 133 | incube.Show();
|
|---|
| 134 |
|
|---|
| 135 | double dxdeg = ThetaSizeDegre/(double)NTheta;
|
|---|
| 136 | double dydeg = PhiSizeDegre/(double)NPhi;
|
|---|
| 137 | double dx = DegreeToRadian(dxdeg);
|
|---|
| 138 | double dy = DegreeToRadian(dydeg);
|
|---|
| 139 | double dfreq = FreqSizeMHz/(double)NFreq;
|
|---|
| 140 |
|
|---|
| 141 | cout << " X,Y map size in degrees , X/Phi=" << PhiSizeDegre << " Y/Theta=" << ThetaSizeDegre
|
|---|
| 142 | << " \n dx=" << dxdeg << " dy=" << dydeg << " degres ( dx_rad=" << dx << " dy_rad=" << dy << ")"
|
|---|
| 143 | << " FreqSize=" << FreqSizeMHz << " dfreq=dz= " << dfreq << " MHz" << endl;
|
|---|
| 144 |
|
|---|
| 145 | double mean, sigma;
|
|---|
| 146 | MeanSigma(incube, mean, sigma);
|
|---|
| 147 | cout << " InCube 3D- : Mean=" << mean << " Sigma=" << sigma << endl;
|
|---|
| 148 |
|
|---|
| 149 | cout << "applobe[2]: creating Four2DResponse and BeamEffect objects..." << endl;
|
|---|
| 150 | H21Conversions conv;
|
|---|
| 151 | conv.setFrequency(Freq0MHz);
|
|---|
| 152 | double lambda = conv.getLambda();
|
|---|
| 153 |
|
|---|
| 154 | int typbeam=(fggaussian)?1:2;
|
|---|
| 155 | Four2DResponse fresp(typbeam, DIAMETRE/lambda, DIAMETRE/lambda, lambda);
|
|---|
| 156 | Four2DResponse* fresp_p=&fresp;
|
|---|
| 157 | Four2DRespTable resptbl;
|
|---|
| 158 | if (fgresptbl) {
|
|---|
| 159 | cout << "applobe[2.b]: initializing Four2DRespTable from file" << resptblname << endl;
|
|---|
| 160 | resptbl.readFromPPF(resptblname);
|
|---|
| 161 | resptbl.renormalize(1.);
|
|---|
| 162 | fresp_p=&resptbl;
|
|---|
| 163 | }
|
|---|
| 164 | else cout << " applobe[2.b]: Four2DResponse ( Diameter=" << DIAMETRE << " Lambda= " << lambda
|
|---|
| 165 | << " DoL=" << DIAMETRE/lambda << " ) " << endl;
|
|---|
| 166 | BeamEffect beam(*fresp_p);
|
|---|
| 167 |
|
|---|
| 168 | if (fgcorrbeam) {
|
|---|
| 169 | double DoL = tbeamDiam/lambda;
|
|---|
| 170 | double tbeamarcmin = RadianToDegree(1.22/DoL)*60.;
|
|---|
| 171 | int typcb = (fggaussian)?1:2;
|
|---|
| 172 | // if (fgresptbl) typcb=22;
|
|---|
| 173 | Four2DResponse tbeam(typcb, DoL, DoL );
|
|---|
| 174 | cout << "applobe[3]: calling Correct2RefLobe() with target beam Diameter=" << tbeamDiam
|
|---|
| 175 | << " D/Lambda=" << DoL << " -> arcmin " << tbeamarcmin << " TypDishResp=" << typcb
|
|---|
| 176 | << " MaxRatio=" << maxratio << endl;
|
|---|
| 177 | beam.Correct2RefLobe(tbeam, incube, dx, dy, Freq0MHz, dfreq, maxratio);
|
|---|
| 178 | }
|
|---|
| 179 | else if (fixedbeam) {
|
|---|
| 180 | cout << "applobe[3]: calling ApplyLobe() (freq.independent beam) ... " << endl;
|
|---|
| 181 | beam.ApplyLobe(incube, dx, dy, Freq0MHz);
|
|---|
| 182 | }
|
|---|
| 183 | else {
|
|---|
| 184 | cout << "applobe[3]: calling ApplyLobe3D() ... " << endl;
|
|---|
| 185 | beam.ApplyLobe3D(incube, dx, dy, Freq0MHz, dfreq);
|
|---|
| 186 | }
|
|---|
| 187 | TArray< r_4 > outcube;
|
|---|
| 188 | if (fgresample) {
|
|---|
| 189 | cout << "applobe[4]: calling ReSample(incube," << fresamp << "," << ",1.) ... " << endl;
|
|---|
| 190 | outcube.Share(beam.ReSample(incube, fresamp, fresamp, 1.));
|
|---|
| 191 | }
|
|---|
| 192 | else outcube.Share(incube);
|
|---|
| 193 |
|
|---|
| 194 | outcube.Show();
|
|---|
| 195 | MeanSigma(outcube, mean, sigma);
|
|---|
| 196 | cout << " OutCube 3D- : Mean=" << mean << " Sigma=" << sigma << endl;
|
|---|
| 197 |
|
|---|
| 198 | // On sauve le cube de sortie
|
|---|
| 199 | {
|
|---|
| 200 | cout << " applobe[5]: Saving output cube to -> " << outname << endl;
|
|---|
| 201 | POutPersist poc(outname);
|
|---|
| 202 | poc << outcube;
|
|---|
| 203 | }
|
|---|
| 204 |
|
|---|
| 205 | rc = 0;
|
|---|
| 206 | }
|
|---|
| 207 | catch (PThrowable& exc) {
|
|---|
| 208 | cerr << " applobe.cc catched SOPHYA Exception " << exc.Msg() << endl;
|
|---|
| 209 | rc = 77;
|
|---|
| 210 | }
|
|---|
| 211 | catch (std::exception& sex) {
|
|---|
| 212 | cerr << "\n applobe.cc std::exception :"
|
|---|
| 213 | << (string)typeid(sex).name() << "\n msg= "
|
|---|
| 214 | << sex.what() << endl;
|
|---|
| 215 | }
|
|---|
| 216 | catch (...) {
|
|---|
| 217 | cerr << " applobe.cc catched unknown (...) exception " << endl;
|
|---|
| 218 | rc = 78;
|
|---|
| 219 | }
|
|---|
| 220 |
|
|---|
| 221 | cout << ">>>> applobe[9] ------- FIN ----------- Rc=" << rc << endl;
|
|---|
| 222 | return rc;
|
|---|
| 223 | }
|
|---|
| 224 |
|
|---|
| 225 |
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
