| 1 | // Utilisation de SOPHYA pour faciliter les tests ...
|
|---|
| 2 | #include "sopnamsp.h"
|
|---|
| 3 | #include "machdefs.h"
|
|---|
| 4 |
|
|---|
| 5 | /* ----------------------------------------------------------
|
|---|
| 6 | Projet BAORadio - (C) LAL/IRFU 2011
|
|---|
| 7 |
|
|---|
| 8 | Programme de lecture des fichiers matrices de visibilites
|
|---|
| 9 | au format fits pour ecrire une DataTable de visib_moyenne(temps)
|
|---|
| 10 | Fichier de visibilites produits par le prog vismf
|
|---|
| 11 | R. Ansari, C. Magneville - LAL/Irfu
|
|---|
| 12 | V : Mai 2009
|
|---|
| 13 | ---------------------------------------------------------- */
|
|---|
| 14 |
|
|---|
| 15 | // include standard c/c++
|
|---|
| 16 | #include <math.h>
|
|---|
| 17 | #include <stdio.h>
|
|---|
| 18 | #include <stdlib.h>
|
|---|
| 19 | #include <string.h>
|
|---|
| 20 |
|
|---|
| 21 | #include <iostream>
|
|---|
| 22 | #include <string>
|
|---|
| 23 |
|
|---|
| 24 | #include "pexceptions.h"
|
|---|
| 25 | #include "tvector.h"
|
|---|
| 26 | #include "fioarr.h"
|
|---|
| 27 | // #include "tarrinit.h"
|
|---|
| 28 | #include "ntuple.h"
|
|---|
| 29 | #include "datatable.h"
|
|---|
| 30 | #include "histinit.h"
|
|---|
| 31 | #include "matharr.h"
|
|---|
| 32 | #include "timestamp.h"
|
|---|
| 33 | #include <utilarr.h>
|
|---|
| 34 | #include "fitsarrhand.h"
|
|---|
| 35 | #include "fitshdtable.h"
|
|---|
| 36 |
|
|---|
| 37 | // include sophya mesure ressource CPU/memoire ...
|
|---|
| 38 | #include "resusage.h"
|
|---|
| 39 | #include "ctimer.h"
|
|---|
| 40 | #include "timing.h"
|
|---|
| 41 |
|
|---|
| 42 |
|
|---|
| 43 | //--------------------------- Fonctions de ce fichier -------------------
|
|---|
| 44 | int Usage(bool fgshort=true);
|
|---|
| 45 |
|
|---|
| 46 | int DecodeProc(int narg, char* arg[]);
|
|---|
| 47 | int ProcVisMtxFiles(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq,
|
|---|
| 48 | vector<sa_size_t> tfrlist);
|
|---|
| 49 | int FillVisDTable(BaseDataTable& visdt_, TMatrix< r_4 > vismtx_, TVector< int_4> chanum_,
|
|---|
| 50 | sa_size_t jf1_, sa_size_t jf2_, sa_size_t djf_) ;
|
|---|
| 51 |
|
|---|
| 52 | //------------------------------------------------------------------------------------------------------------
|
|---|
| 53 |
|
|---|
| 54 | /* --Fonction-- */
|
|---|
| 55 | int Usage(bool fgshort)
|
|---|
| 56 | {
|
|---|
| 57 | cout << " --- visfits2dt : Read fits visibility matrices produced by visfmib to make \n"
|
|---|
| 58 | << " time-frequency matrices and Visibilites=f(time) datatable " << endl;
|
|---|
| 59 | cout << " Usage: visfits2dt InOutPath Imin,Imax,step NumFreq1,NumFreq2,NBinFreq [VisMtxRowList] \n" << endl;
|
|---|
| 60 | if (fgshort) {
|
|---|
| 61 | cout << " svv2mtx -h for detailed instructions" << endl;
|
|---|
| 62 | return 1;
|
|---|
| 63 | }
|
|---|
| 64 | cout << " InOutPath : Input/Output directory name \n"
|
|---|
| 65 | << " Imin,Imax,IStep: Input FITS files sequence number (vismtxII.fits) \n"
|
|---|
| 66 | << " FileNames=InOutPath/vismtxII.fits Imin<=II<=Imax II+=IStep \n"
|
|---|
| 67 | << " NumFreq1,NumFreq2,NBinFreq: Freq Zone and number of frequency bins for DataTable\n"
|
|---|
| 68 | << " VisMtxRowList : List of visibiliy matrix rows (0,2,...) -> time-freq map \n" << endl;
|
|---|
| 69 | return 1;
|
|---|
| 70 | }
|
|---|
| 71 |
|
|---|
| 72 | //----------------------------------------------------
|
|---|
| 73 | //----------------------------------------------------
|
|---|
| 74 | int main(int narg, char* arg[])
|
|---|
| 75 | {
|
|---|
| 76 | if ((narg>1)&&(strcmp(arg[1],"-h")==0)) return Usage(false);
|
|---|
| 77 | if (narg<4) return Usage(true);
|
|---|
| 78 | HiStatsInitiator _inia;
|
|---|
| 79 | // TArrayInitiator _inia;
|
|---|
| 80 |
|
|---|
| 81 | int rc = 0;
|
|---|
| 82 | try {
|
|---|
| 83 | bool fgvjun09=false;
|
|---|
| 84 | ResourceUsage resu;
|
|---|
| 85 | DecodeProc(narg-1, arg+1);
|
|---|
| 86 | resu.Update();
|
|---|
| 87 | cout << resu;
|
|---|
| 88 | }
|
|---|
| 89 | catch (PException& exc) {
|
|---|
| 90 | cerr << " visfits2dt.cc catched PException " << exc.Msg() << endl;
|
|---|
| 91 | rc = 77;
|
|---|
| 92 | }
|
|---|
| 93 | catch (std::exception& sex) {
|
|---|
| 94 | cerr << "\n visfits2dt.cc std::exception :"
|
|---|
| 95 | << (string)typeid(sex).name() << "\n msg= "
|
|---|
| 96 | << sex.what() << endl;
|
|---|
| 97 | rc = 78;
|
|---|
| 98 | }
|
|---|
| 99 | catch (...) {
|
|---|
| 100 | cerr << " visfits2dt.cc catched unknown (...) exception " << endl;
|
|---|
| 101 | rc = 79;
|
|---|
| 102 | }
|
|---|
| 103 |
|
|---|
| 104 | cout << ">>>> visfits2dt.cc ------- END ----------- RC=" << rc << endl;
|
|---|
| 105 | return rc;
|
|---|
| 106 |
|
|---|
| 107 | }
|
|---|
| 108 |
|
|---|
| 109 | //--------------------------------------------------------------------
|
|---|
| 110 | // Traitement fichiers produits par vismfib (V Nov09)
|
|---|
| 111 | /* --Fonction-- */
|
|---|
| 112 | int DecodeProc(int narg, char* arg[])
|
|---|
| 113 | {
|
|---|
| 114 | // Decodage des arguments et traitement
|
|---|
| 115 | string inoutpath = arg[0];
|
|---|
| 116 | int imin=0;
|
|---|
| 117 | int imax=0;
|
|---|
| 118 | int istep=1;
|
|---|
| 119 | sscanf(arg[1],"%d,%d,%d",&imin,&imax,&istep);
|
|---|
| 120 | int jf1=0;
|
|---|
| 121 | int jf2=0;
|
|---|
| 122 | int nfreq=0;
|
|---|
| 123 | sscanf(arg[2],"%d,%d,%d",&jf1,&jf2,&nfreq);
|
|---|
| 124 | int card=1;
|
|---|
| 125 | vector<sa_size_t> rowlist;
|
|---|
| 126 | if (narg>3) {
|
|---|
| 127 | EnumeratedSequence es;
|
|---|
| 128 | int nbad;
|
|---|
| 129 | es.Append(arg[3], nbad, ",");
|
|---|
| 130 | for(int j=0; j<es.Size(); j++) rowlist.push_back((int)es.Value(j));
|
|---|
| 131 | }
|
|---|
| 132 | if (rowlist.size()<1) rowlist.push_back(0);
|
|---|
| 133 |
|
|---|
| 134 | cout << " ----- visfits2dt/DecodeProc - Start - InOutPath= " << inoutpath << " IMin,Max,Step="
|
|---|
| 135 | << imin << "," << imax << "," << istep << " Card=" << card << endl;
|
|---|
| 136 | cout << "Frequency num range JF=" << jf1 << "," << jf2 << "," << nfreq << " ------- " << endl;
|
|---|
| 137 | cout << " RowList: " ;
|
|---|
| 138 | for(int j=0; j<rowlist.size(); j++) cout << rowlist[j] << " , " ;
|
|---|
| 139 | cout << endl;
|
|---|
| 140 | int rc=ProcVisMtxFiles(inoutpath, imin, imax, istep, jf1, jf2, nfreq, rowlist);
|
|---|
| 141 | return rc;
|
|---|
| 142 | }
|
|---|
| 143 |
|
|---|
| 144 | // Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
|
|---|
| 145 | /* --Fonction-- */
|
|---|
| 146 | int ProcVisMtxFiles(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq,
|
|---|
| 147 | vector<sa_size_t> rowlist)
|
|---|
| 148 | {
|
|---|
| 149 | Timer tm("ProcVisMtxFiles");
|
|---|
| 150 |
|
|---|
| 151 | DataTable visdt;
|
|---|
| 152 | visdt.AddDoubleColumn("mfc");
|
|---|
| 153 | visdt.AddDoubleColumn("mtt");
|
|---|
| 154 | visdt.AddIntegerColumn("jfreq");
|
|---|
| 155 | visdt.AddIntegerColumn("numch");
|
|---|
| 156 | visdt.AddFloatColumn("vre");
|
|---|
| 157 | visdt.AddFloatColumn("vim");
|
|---|
| 158 |
|
|---|
| 159 | vector< TMatrix< r_4 > > vmtf;
|
|---|
| 160 |
|
|---|
| 161 | if (jf1<1) jf1=1;
|
|---|
| 162 | if ((jf2<1)||(jf2<jf1)) jf2=jf1;
|
|---|
| 163 | if (nfreq<1) nfreq=1;
|
|---|
| 164 | int djf=(jf2-jf1)/nfreq;
|
|---|
| 165 | if (djf<1) djf=1;
|
|---|
| 166 |
|
|---|
| 167 | cout << " --- ProcVisMtxFiles jf1=" << jf1 << " jf2= " << jf2 << " djf=" << djf << endl;
|
|---|
| 168 | char fname[1024];
|
|---|
| 169 |
|
|---|
| 170 | cout << " ---- ProcVisMtxFiles()-Begin - Reading chanum vector" << endl;
|
|---|
| 171 | TVector< int_4 > chanum;
|
|---|
| 172 | {
|
|---|
| 173 | sprintf(fname, "%s/chanum.fits",inoutpath.c_str());
|
|---|
| 174 | FitsInOutFile fic(fname, FitsInOutFile::Fits_RO);
|
|---|
| 175 | fic.MoveAbsToHDU(3);
|
|---|
| 176 | fic >> chanum;
|
|---|
| 177 | }
|
|---|
| 178 |
|
|---|
| 179 | cout << " ---- ProcVisMtxFiles()-Read chanum done " << endl;
|
|---|
| 180 | sa_size_t nrows = (imax-imin+1)/istep;
|
|---|
| 181 | sa_size_t kr=0;
|
|---|
| 182 |
|
|---|
| 183 | sa_size_t mtf_binfreq=25;
|
|---|
| 184 | sa_size_t mtf_bintime=5;
|
|---|
| 185 |
|
|---|
| 186 | sa_size_t ncols=1;
|
|---|
| 187 |
|
|---|
| 188 | for(int ifile=imin; ifile<=imax; ifile+=istep) {
|
|---|
| 189 | sprintf(fname, "%s/vismtx%d.fits",inoutpath.c_str(),ifile);
|
|---|
| 190 | cout << " ProcVisMtxFiles[" << ifile << "] opening file " << fname << endl;
|
|---|
| 191 | FitsInOutFile fin(fname, FitsInOutFile::Fits_RO);
|
|---|
| 192 | TMatrix< r_4 > vismtx;
|
|---|
| 193 | fin >> vismtx;
|
|---|
| 194 | // cout << vismtx.Info();
|
|---|
| 195 |
|
|---|
| 196 | if (ifile==imin) {
|
|---|
| 197 | ncols = vismtx.NCols();
|
|---|
| 198 | cout << " ProcVisMtxFiles/Info: Output Time-Frequency matrices NRows(time) "
|
|---|
| 199 | << nrows/mtf_bintime+1 << " NCols(freq) =" << ncols/mtf_binfreq+1 << endl;
|
|---|
| 200 | for(size_t j=0; j<rowlist.size(); j++)
|
|---|
| 201 | vmtf.push_back(TMatrix< r_4 >(nrows/mtf_bintime+1, ncols/mtf_binfreq+1));
|
|---|
| 202 | }
|
|---|
| 203 |
|
|---|
| 204 | for(size_t j=0; j<rowlist.size(); j++)
|
|---|
| 205 | for(sa_size_t icf=0; icf<ncols; icf++) {
|
|---|
| 206 | vmtf[j](icf/mtf_binfreq,kr/mtf_bintime)+=vismtx(rowlist[j],icf);
|
|---|
| 207 | }
|
|---|
| 208 | // cout << " DBG* kr=" << kr << " kr/mtf_bintime=" << kr/mtf_bintime
|
|---|
| 209 | // << " ncols/2/mtf_binfreq=" << ncols/2/mtf_binfreq << endl;
|
|---|
| 210 | kr++;
|
|---|
| 211 |
|
|---|
| 212 |
|
|---|
| 213 | // Calcul moyenne dans des bandes en frequence
|
|---|
| 214 | FillVisDTable(visdt, vismtx, chanum, jf1, jf2, djf);
|
|---|
| 215 | }
|
|---|
| 216 |
|
|---|
| 217 | if (rowlist.size()>0) {
|
|---|
| 218 | sprintf(fname, "!%s/vistfreqmtx.fits",inoutpath.c_str());
|
|---|
| 219 | cout << "ProcVisMtxFiles: Opening file " << fname << " for writing Visi(Freq,Time) matrices" << endl;
|
|---|
| 220 | FitsInOutFile fo(fname, FitsInOutFile::Fits_Create);
|
|---|
| 221 | for(size_t j=0; j<rowlist.size(); j++)
|
|---|
| 222 | fo << vmtf[j];
|
|---|
| 223 | }
|
|---|
| 224 | cout << visdt;
|
|---|
| 225 | sprintf(fname, "!%s/visidt.ppf",inoutpath.c_str());
|
|---|
| 226 | cout << "ProcVisMtxFiles: writing visibility datatable to file " << fname << endl;
|
|---|
| 227 | FitsInOutFile fod(fname, FitsInOutFile::Fits_Create);
|
|---|
| 228 | fod << visdt;
|
|---|
| 229 |
|
|---|
| 230 | return 0;
|
|---|
| 231 | }
|
|---|
| 232 |
|
|---|
| 233 | /* --Fonction-- */
|
|---|
| 234 | int FillVisDTable(BaseDataTable& visdt_, TMatrix< r_4 > vismtx_, TVector< int_4> chanum_,
|
|---|
| 235 | sa_size_t jf1_, sa_size_t jf2_, sa_size_t djf_)
|
|---|
| 236 | {
|
|---|
| 237 | double xnt_[20];
|
|---|
| 238 | xnt_[0]=(double)vismtx_.Info()["MEANFC"];
|
|---|
| 239 | xnt_[1]=(double)vismtx_.Info()["MEANTT"]/1.25e8;
|
|---|
| 240 |
|
|---|
| 241 | uint_4 nmean_=vismtx_.Info()["NPAQSUM"];
|
|---|
| 242 | if (djf_<2) {
|
|---|
| 243 | for(sa_size_t rv=0; rv<vismtx_.NRows(); rv++) {
|
|---|
| 244 | for(sa_size_t jf=jf1_; jf<jf2_; jf++) {
|
|---|
| 245 | xnt_[2]=jf;
|
|---|
| 246 | xnt_[3]=chanum_(rv);
|
|---|
| 247 | xnt_[4]=vismtx_(rv,jf*2)/(r_4)(nmean_);
|
|---|
| 248 | xnt_[5]=vismtx_(rv,jf*2+1)/(r_4)(nmean_);
|
|---|
| 249 | visdt_.AddRow(xnt_);
|
|---|
| 250 | }
|
|---|
| 251 | }
|
|---|
| 252 | }
|
|---|
| 253 | else {
|
|---|
| 254 | for(sa_size_t rv=0; rv<vismtx_.NRows(); rv++) {
|
|---|
| 255 | for(sa_size_t jf=jf1_; jf<jf2_; jf+=djf_) {
|
|---|
| 256 | r_4 moyreal=0.;
|
|---|
| 257 | r_4 moyimag=0.;
|
|---|
| 258 | sa_size_t jjfmx=jf+djf_;
|
|---|
| 259 | if (jjfmx > vismtx_.NCols()) jjfmx=vismtx_.NCols();
|
|---|
| 260 | for(sa_size_t jjf=jf; jjf<jjfmx; jjf++) {
|
|---|
| 261 | moyreal+=vismtx_(rv,jjf*2);
|
|---|
| 262 | moyimag+=vismtx_(rv,jjf*2+1);
|
|---|
| 263 | }
|
|---|
| 264 | xnt_[2]=jf+djf_/2;
|
|---|
| 265 | xnt_[3]=chanum_(rv);
|
|---|
| 266 | xnt_[4]=moyreal/(r_4)(nmean_*djf_);
|
|---|
| 267 | xnt_[5]=moyimag/(r_4)(nmean_*djf_);
|
|---|
| 268 | visdt_.AddRow(xnt_);
|
|---|
| 269 | }
|
|---|
| 270 | }
|
|---|
| 271 | }
|
|---|
| 272 | return 0;
|
|---|
| 273 | }
|
|---|
| 274 |
|
|---|
| 275 |
|
|---|
| 276 |
|
|---|
| 277 |
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
