| 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 -kzf/-nzf -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 | << " -kzf -nzf : Keep (default) or Not zero space frequency when applying lobes (BeamEffect class) \n"
|
|---|
| 66 | << " -mxr val: Max beam correction factor (default=10.) \n"
|
|---|
| 67 | << " Diameter/Four2DRespTableFile : dish diameter or 2D response PPF file name\n"
|
|---|
| 68 | << " In3DPPFName Out3DPPFName: Input/Output PPF file names (TArray<r_4> 3D) \n"
|
|---|
| 69 | << " TargetBeamDiam: Corrected beam target diameter (D/Lambda at lowest frequency) \n"
|
|---|
| 70 | << " DoL = 100 --> beam ~ 35 arcmin (D=30m @ z~0.5 Lambda~30cm)) \n"
|
|---|
| 71 | << " ResmapleFactor : Resampling (2-> oversampling, 0.5 : 1/2 undersampling) \n" << endl;
|
|---|
| 72 | }
|
|---|
| 73 | return 1;
|
|---|
| 74 | }
|
|---|
| 75 |
|
|---|
| 76 | Timer tm("applobe");
|
|---|
| 77 |
|
|---|
| 78 | // decodage arguments
|
|---|
| 79 | bool fggaussian=true; // true -> gaussian beam
|
|---|
| 80 | bool fixedbeam=false; // true -> apply freq. independent beam
|
|---|
| 81 | double maxratio=10.; // valeur max du rapport des lobes lors de la correction de lobe
|
|---|
| 82 | bool preservefreq0=true; // true -> keep zero frequency when appyling lobe
|
|---|
| 83 |
|
|---|
| 84 | // decodage argument optionnel
|
|---|
| 85 | bool fgoptarg=true;
|
|---|
| 86 | while (fgoptarg) {
|
|---|
| 87 | string fbo = arg[1];
|
|---|
| 88 | if (fbo=="-t") { fggaussian=false; arg++; narg--; }
|
|---|
| 89 | else if (fbo=="-g") { fggaussian=true; arg++; narg--; }
|
|---|
| 90 | else if (fbo=="-fib") { fixedbeam=true; arg++; narg--; }
|
|---|
| 91 | else if (fbo=="-kzf") { preservefreq0=true; arg++; narg--; }
|
|---|
| 92 | else if (fbo=="-nzf") { preservefreq0=false; arg++; narg--; }
|
|---|
| 93 | else if (fbo=="-mxr") { arg++; maxratio=atof(arg[1]); arg++; narg-=2; }
|
|---|
| 94 | else fgoptarg=false;
|
|---|
| 95 | }
|
|---|
| 96 | if (narg < 4) {
|
|---|
| 97 | cout << " applobe/error arguments , applobe -h for help " << endl;
|
|---|
| 98 | return 2;
|
|---|
| 99 | }
|
|---|
| 100 |
|
|---|
| 101 | bool fgresptbl=true;
|
|---|
| 102 | double DIAMETRE=100.;
|
|---|
| 103 | string resptblname;
|
|---|
| 104 | if (isdigit(*arg[1])) {
|
|---|
| 105 | fgresptbl=false;
|
|---|
| 106 | DIAMETRE=atof(arg[1]);
|
|---|
| 107 | }
|
|---|
| 108 | else resptblname=arg[1];
|
|---|
| 109 |
|
|---|
| 110 | string inname = arg[2];
|
|---|
| 111 | string outname = arg[3];
|
|---|
| 112 | bool fgcorrbeam=false;
|
|---|
| 113 | double tbeamDiam=50.;
|
|---|
| 114 | if (narg>4) {
|
|---|
| 115 | tbeamDiam=atof(arg[4]);
|
|---|
| 116 | if (tbeamDiam>1.) fgcorrbeam=true;
|
|---|
| 117 | }
|
|---|
| 118 | bool fgresample=false;
|
|---|
| 119 | double fresamp=1.;
|
|---|
| 120 | if (narg>5) {
|
|---|
| 121 | fresamp=atof(arg[5]);
|
|---|
| 122 | if ((fabs(fresamp)>1.e-2)&&(fabs(fresamp-1.)>1.e-2)) fgresample=true;
|
|---|
| 123 | }
|
|---|
| 124 |
|
|---|
| 125 | int rc = 91;
|
|---|
| 126 |
|
|---|
| 127 | cout << " ====== applobe : Input skycube name= " << inname << " OutName=" << outname << endl;
|
|---|
| 128 | cout << ((fggaussian)?" Gaussian ":" Triangular") << ((fixedbeam)?" Fixed (freq.independent)":"") << " beams" << endl;
|
|---|
| 129 | bool fginmap=true;
|
|---|
| 130 | try {
|
|---|
| 131 | TArray<r_4> incube;
|
|---|
| 132 | cout << "applobe[1]: reading input 3D map (cube) from file " << inname << endl;
|
|---|
| 133 | {
|
|---|
| 134 | PInPersist pin(inname);
|
|---|
| 135 | pin >> incube;
|
|---|
| 136 | }
|
|---|
| 137 | incube.Show();
|
|---|
| 138 |
|
|---|
| 139 | double dxdeg = ThetaSizeDegre/(double)NTheta;
|
|---|
| 140 | double dydeg = PhiSizeDegre/(double)NPhi;
|
|---|
| 141 | double dx = DegreeToRadian(dxdeg);
|
|---|
| 142 | double dy = DegreeToRadian(dydeg);
|
|---|
| 143 | double dfreq = FreqSizeMHz/(double)NFreq;
|
|---|
| 144 |
|
|---|
| 145 | cout << " X,Y map size in degrees , X/Phi=" << PhiSizeDegre << " Y/Theta=" << ThetaSizeDegre
|
|---|
| 146 | << " \n dx=" << dxdeg << " dy=" << dydeg << " degres ( dx_rad=" << dx << " dy_rad=" << dy << ")"
|
|---|
| 147 | << " FreqSize=" << FreqSizeMHz << " dfreq=dz= " << dfreq << " MHz" << endl;
|
|---|
| 148 |
|
|---|
| 149 | double mean, sigma;
|
|---|
| 150 | MeanSigma(incube, mean, sigma);
|
|---|
| 151 | cout << " InCube 3D- : Mean=" << mean << " Sigma=" << sigma << endl;
|
|---|
| 152 |
|
|---|
| 153 | cout << "applobe[2]: creating Four2DResponse and BeamEffect objects..." << endl;
|
|---|
| 154 | H21Conversions conv;
|
|---|
| 155 | conv.setFrequency(Freq0MHz);
|
|---|
| 156 | double lambda = conv.getLambda();
|
|---|
| 157 |
|
|---|
| 158 | int typbeam=(fggaussian)?1:2;
|
|---|
| 159 | Four2DResponse fresp(typbeam, DIAMETRE/lambda, DIAMETRE/lambda, lambda);
|
|---|
| 160 | Four2DResponse* fresp_p=&fresp;
|
|---|
| 161 | Four2DRespTable resptbl;
|
|---|
| 162 | if (fgresptbl) {
|
|---|
| 163 | cout << "applobe[2.b]: initializing Four2DRespTable from file" << resptblname << endl;
|
|---|
| 164 | resptbl.readFromPPF(resptblname);
|
|---|
| 165 | resptbl.renormalize(1.);
|
|---|
| 166 | fresp_p=&resptbl;
|
|---|
| 167 | }
|
|---|
| 168 | else cout << " applobe[2.b]: Four2DResponse ( Diameter=" << DIAMETRE << " Lambda= " << lambda
|
|---|
| 169 | << " DoL=" << DIAMETRE/lambda << " ) " << endl;
|
|---|
| 170 | BeamEffect beam(*fresp_p,preservefreq0);
|
|---|
| 171 |
|
|---|
| 172 | if (fgcorrbeam) {
|
|---|
| 173 | double DoL = tbeamDiam/lambda;
|
|---|
| 174 | double tbeamarcmin = RadianToDegree(1.22/DoL)*60.;
|
|---|
| 175 | int typcb = (fggaussian)?1:2;
|
|---|
| 176 | // if (fgresptbl) typcb=22;
|
|---|
| 177 | Four2DResponse tbeam(typcb, DoL, DoL );
|
|---|
| 178 | cout << "applobe[3]: calling Correct2RefLobe() with target beam Diameter=" << tbeamDiam
|
|---|
| 179 | << " D/Lambda=" << DoL << " -> arcmin " << tbeamarcmin << " TypDishResp=" << typcb
|
|---|
| 180 | << " MaxRatio=" << maxratio << endl;
|
|---|
| 181 | beam.Correct2RefLobe(tbeam, incube, dx, dy, Freq0MHz, dfreq, maxratio);
|
|---|
| 182 | }
|
|---|
| 183 | else if (fixedbeam) {
|
|---|
| 184 | cout << "applobe[3]: calling ApplyLobe() (freq.independent beam) ... " << endl;
|
|---|
| 185 | beam.ApplyLobe(incube, dx, dy, Freq0MHz);
|
|---|
| 186 | }
|
|---|
| 187 | else {
|
|---|
| 188 | cout << "applobe[3]: calling ApplyLobe3D() ... " << endl;
|
|---|
| 189 | beam.ApplyLobe3D(incube, dx, dy, Freq0MHz, dfreq);
|
|---|
| 190 | }
|
|---|
| 191 | TArray< r_4 > outcube;
|
|---|
| 192 | if (fgresample) {
|
|---|
| 193 | cout << "applobe[4]: calling ReSample(incube," << fresamp << "," << ",1.) ... " << endl;
|
|---|
| 194 | outcube.Share(beam.ReSample(incube, fresamp, fresamp, 1.));
|
|---|
| 195 | }
|
|---|
| 196 | else outcube.Share(incube);
|
|---|
| 197 |
|
|---|
| 198 | outcube.Show();
|
|---|
| 199 | MeanSigma(outcube, mean, sigma);
|
|---|
| 200 | cout << " OutCube 3D- : Mean=" << mean << " Sigma=" << sigma << endl;
|
|---|
| 201 |
|
|---|
| 202 | // On sauve le cube de sortie
|
|---|
| 203 | {
|
|---|
| 204 | cout << " applobe[5]: Saving output cube to -> " << outname << endl;
|
|---|
| 205 | POutPersist poc(outname);
|
|---|
| 206 | poc << outcube;
|
|---|
| 207 | }
|
|---|
| 208 |
|
|---|
| 209 | rc = 0;
|
|---|
| 210 | }
|
|---|
| 211 | catch (PThrowable& exc) {
|
|---|
| 212 | cerr << " applobe.cc catched SOPHYA Exception " << exc.Msg() << endl;
|
|---|
| 213 | rc = 77;
|
|---|
| 214 | }
|
|---|
| 215 | catch (std::exception& sex) {
|
|---|
| 216 | cerr << "\n applobe.cc std::exception :"
|
|---|
| 217 | << (string)typeid(sex).name() << "\n msg= "
|
|---|
| 218 | << sex.what() << endl;
|
|---|
| 219 | }
|
|---|
| 220 | catch (...) {
|
|---|
| 221 | cerr << " applobe.cc catched unknown (...) exception " << endl;
|
|---|
| 222 | rc = 78;
|
|---|
| 223 | }
|
|---|
| 224 |
|
|---|
| 225 | cout << ">>>> applobe[9] ------- FIN ----------- Rc=" << rc << endl;
|
|---|
| 226 | return rc;
|
|---|
| 227 | }
|
|---|
| 228 |
|
|---|
| 229 |
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
