visiavg.cc 28.6 KB
Newer Older
OP's avatar
OP committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
// Utilisation de SOPHYA pour faciliter les tests ...
#include "sopnamsp.h"
#include "machdefs.h"

/* ---------------------------------------------------------- 
   Projet BAORadio/PAON4 - (C) LAL/IRFU  2017

   visiavg: programme de lecture des fichiers matrices de 
   visibilites de PAON4, calcul de visibilities moyennes 
    en bin de temps et de frequence
   O. Perdereau, R.Ansari   -  LAL
   ---------------------------------------------------------- */

// include standard c/c++
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <iostream>
#include <string>

#include "pexceptions.h"
#include "tvector.h"
#include "fioarr.h"
// #include "tarrinit.h"
#include "ntuple.h" 
#include "datatable.h" 
#include "histinit.h" 
#include "matharr.h" 
#include "timestamp.h"
#include <utilarr.h>

// include sophya mesure ressource CPU/memoire ...
#include "resusage.h"
#include "ctimer.h"
#include "timing.h"

// include lecteur de fichiers visibilites 
40
#include "visp4winreader.h"
OP's avatar
OP committed
41

OP's avatar
OP committed
42 43 44
#include "fitsioserver.h"
#include "fiosinit.h"

45
int Usage(bool fglong)
OP's avatar
OP committed
46
{
47 48 49
  cout << " --- visiavg.cc : Computes time-frequency visibility maps from time stream visibility data\n" 
       << " Usage: visiavg [-arguments] \n" << endl;
  P4AnaParams::UsageOptions(fglong);
OP's avatar
OP committed
50 51 52 53 54 55 56 57
  cout<< endl;
  return 1;
}

//----------------------------------------------------
int main(int narg, const char* arg[])
{
  // --- Decoding parameters 
58 59 60 61 62 63 64 65 66
  if (narg<2) {
    cout << " --- visiavg : computes TFM (time-frequency) maps from visibilities \n" 
	 << " visiavg/ Missing argument - To get short/long usage information type \n"
	 << " visiavg -h (list of arguments) OR visiavg -help  (list of arguments + description)"<<endl; 
    return 1;
  }
  else if ((narg>1)&&(strcmp(arg[1],"-h")==0) )   return Usage(false);
  else if ((narg>1)&&(strcmp(arg[1],"-help")==0) )   return Usage(true);

OP's avatar
OP committed
67 68
  P4AnaParams params;
  params.DecodeArgs(narg, arg);
OP's avatar
OP committed
69

OP's avatar
OP committed
70 71
  string outfile = params.outfile_;
  if (outfile.length()<1)  outfile = "visavg.ppf";
OP's avatar
OP committed
72

OP's avatar
OP committed
73 74 75 76
  string fitsoutfile = params.fitsoutfile_;
  if (fitsoutfile.length()>=1) {
    fitsoutfile = "!"+fitsoutfile ; // adds '!' ?
  }
OP's avatar
modifs  
OP committed
77

OP's avatar
OP committed
78 79 80 81 82 83 84 85 86 87 88
  int deltaIavg = params.TFMtimebin_;
  sa_size_t TFMfbin = params.TFMfreqbin_;
  int prtlev = params.prtlev_;
  bool FgTFMAC = true;
  bool FgTFMCX = true;
  string desctfmap;
  bool FgTFM = params.fgTFM_;   // true -> create time-frequency maps

  params.Print(cout);
  cout <<"visiavg/Info: Path BAO5:"<<params.inpath5_<<" BAO6:"<<params.inpath6_<<"\n"
       <<"fgreorderfreq="<<params.fgreorderfreq_<<"\n"
89
       <<" DeltaIAvg="<<deltaIavg<<"\n"
perdereau's avatar
perdereau committed
90 91
       << " fband # "<<params.fbands_.size() 
       <<" outfile="<<outfile<<" PrtLev="<<prtlev<<endl;
OP's avatar
OP committed
92 93 94 95 96 97 98 99 100

  if (!FgTFM) {
    cout<<" visiavg/parameter error : specify Time-Frequency map parameter with -tfm "<<endl;
    return 5;
  }

  P4AVisiNumEncoder  visiencod;
  vector<sa_size_t> KVAC = visiencod.getAllAutoCor();
  vector<sa_size_t> KVCXHH = visiencod.getAllHCrossCor();
101 102
  vector<sa_size_t> KVCXVV = visiencod.getAllVCrossCor();
  vector<sa_size_t> KVCXHV = visiencod.getAllHVCrossCor();
103 104

  
OP's avatar
OP committed
105 106 107 108
  cout << " List of AutoCorrelation rows:"<<endl;
  for(size_t k=0; k<KVAC.size(); k++) {
    cout << "KVAC["<<k<<"]="<<KVAC[k]<<"  ->"<<visiencod.Convert2VisiName(KVAC[k])<<endl;
  }
109

OP's avatar
OP committed
110 111 112 113
  cout << " List of HH X-cor rows:"<<endl;
  for(size_t k=0; k<KVCXHH.size(); k++) {
    cout << "KVCXHH["<<k<<"]="<<KVCXHH[k]<<"  ->"<<visiencod.Convert2VisiName(KVCXHH[k])<<endl;
  }
OP's avatar
modifs  
OP committed
114

OP's avatar
OP committed
115
  int rc = 0;
OP's avatar
modifs  
OP committed
116

OP's avatar
OP committed
117
  try {
118 119 120
    // ---  Init SOPHYA 
    HiStatsInitiator _inia;
    FitsIOServerInit();
OP's avatar
OP committed
121 122
    ResourceUsage resu;

123
    VisiP4WindowReader wreader(params);
124

125 126 127
    long Imin = wreader.getReader().getSerialFirst();
    long Imax =  wreader.getReader().getSerialLast();
    long Istep =  wreader.getReader().getSerialStep();
OP's avatar
OP committed
128 129
    cout << "visiavg/Info: processing visibility matrix serial/sequence number range "
	 <<Imin<<" <= seq <= " << Imax << " with step="<<Istep<<endl;
130
    cout << " WindowSize="<<wreader.getWindowSize()<<"  -> TotalNbWindows="<<wreader.getTotalNbWindows()<<endl;
OP's avatar
OP committed
131 132

    bool fgok=true;
OP's avatar
OP committed
133 134
    // vecteur de noms 
    vector <string> ext_names;
perdereau's avatar
perdereau committed
135
    // un vecteur avec les temps & un avec RA
136 137
    TVector< double > timevec(wreader.getTotalNbWindows()/deltaIavg); 
    TVector< double > ravec(wreader.getTotalNbWindows()/deltaIavg); 
OP's avatar
OP committed
138 139
    TMatrix< complex<r_4> > vismtx;
    TMatrix< complex<r_4> > acsum;
OP's avatar
OP committed
140
    TMatrix< r_4 > acsum_sq; // wil sum only the real part ^2
141

OP's avatar
OP committed
142
    TMatrix< complex<r_4> > cxsum;
OP's avatar
OP committed
143 144 145
    // for sums of real and imag parts 
    TMatrix< r_4 > cxsum_sq_rp;
    TMatrix< r_4 > cxsum_sq_ip;
OP's avatar
OP committed
146

147 148 149 150 151 152 153 154 155 156 157 158
    // for VV if needed
    TMatrix< complex<r_4> > cxsum_vv;
    // for sums of real and imag parts 
    TMatrix< r_4 > cxsum_vv_sq_rp;
    TMatrix< r_4 > cxsum_vv_sq_ip;

    // for HV if needed
    TMatrix< complex<r_4> > cxsum_hv;
    // for sums of real and imag parts 
    TMatrix< r_4 > cxsum_hv_sq_rp;
    TMatrix< r_4 > cxsum_hv_sq_ip;

159
    TimeStamp dateobs, cfdate,datestart,datused;
OP's avatar
OP committed
160 161 162 163
    TimeStamp dateorg(2015,1,1,12,0,0.);  // Date origine 1 jan 2015
    double mttag;
    int cnt=0, cntnt=0, pcntnt=0;
    int I=0; 
OP's avatar
OP committed
164
    // for
165

OP's avatar
OP committed
166 167
    //----- 6 H-H cross-cor TimeFrequency maps 
    vector< TArray< complex<r_4> > > vtfm;
OP's avatar
OP committed
168 169 170 171
    //----- 6 H-H cross-cor TimeFrequency maps for the variances of real and imag parts 
    vector< TArray< r_4 > > vtfm_rp_sq;
    vector< TArray< r_4 > > vtfm_ip_sq;

172 173 174 175 176 177 178 179 180 181 182 183
    //----- 6 V-V cross-cor TimeFrequency maps 
    vector< TArray< complex<r_4> > > vtfm_vv;
    //----- 6 V-V cross-cor TimeFrequency maps for the variances of real and imag parts 
    vector< TArray< r_4 > > vtfm_vv_rp_sq;
    vector< TArray< r_4 > > vtfm_vv_ip_sq;

    //----- 16 H-V cross-cor TimeFrequency maps 
    vector< TArray< complex<r_4> > > vtfm_hv;
    //----- 16 H-V cross-cor TimeFrequency maps for the variances of real and imag parts 
    vector< TArray< r_4 > > vtfm_hv_rp_sq;
    vector< TArray< r_4 > > vtfm_hv_ip_sq;

OP's avatar
OP committed
184 185
    //----- 8 auto-corr TimeFrequency maps 
    vector< TArray< r_4 > > vtfmac;
OP's avatar
OP committed
186 187
    //----- 8 auto-corr TimeFrequency variance maps 
    vector< TArray< r_4 > > vtfmac_sq;
OP's avatar
OP committed
188 189 190 191
    
    //---- for the time-freqency map filling    
    sa_size_t TFMtmidx=0;
    sa_size_t tfmSX, tfmSY;
192 193
    long serialnum; 
    size_t rdidx, lastrdidx=0;
perdereau's avatar
perdereau committed
194 195 196 197 198 199

    double freqmin =  P4FreqBand::freqstart_;
    double freqmax = P4FreqBand::freqend_;
    sa_size_t ifrmin =  0 ;
    sa_size_t ifrmax = P4FreqBand::nbnufft_ ;

OP's avatar
OP committed
200
    while (fgok) {
201
      //reads next visibility matrix window 
202
      fgok = wreader.Shift(rdidx, serialnum);
203
      
204
      if (!fgok)  break;
205 206 207 208 209 210 211
      if ( (rdidx <= lastrdidx)&&(lastrdidx>0)) throw ForbiddenError("BAD wreader.Shift(rdidx) rdidx <= lastrdidx");  // BUG si ca arrive !
      if ((rdidx-lastrdidx)>1  ){
	TFMtmidx+=(rdidx-lastrdidx-1)/deltaIavg;
	cout<<" visiavg : skipped some data "<<  (rdidx-lastrdidx-1) << endl;
      }
      lastrdidx=rdidx;

212
      vismtx = wreader.getAverageVisMtx(cfdate);
perdereau's avatar
perdereau committed
213 214
      ///cout << vismtx << endl;
      sa_size_t NTFCols = vismtx.NCols(); // bin en freq a ecrire 
215
      
216
      if (cnt==0)  {    //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations 
217

perdereau's avatar
perdereau committed
218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236
	tfmSX=wreader.getTotalNbWindows()/deltaIavg;

	tfmSY=vismtx.NCols()/TFMfbin;
	// default 

	if ( params.fbands_.size()>0 ){
	  sa_size_t XtfmSY = (sa_size_t) floor(params.fbands_[0].df_/params.fbands_[0].deltanufft_ /TFMfbin ) ;
	  cout << " taille finale en freq "<< XtfmSY << " bin "<< TFMfbin <<endl;
	  tfmSY=XtfmSY;
	  ifrmin = (sa_size_t) params.fbands_[0].jfmin_ ;
	  ifrmax = (sa_size_t) params.fbands_[0].jfmax_ ;
	  freqmin = params.fbands_[0].getP4Frequency(ifrmin);
	  freqmax = params.fbands_[0].getP4Frequency(ifrmax);
	  NTFCols = ifrmax - ifrmin +1 ; //## CHECK THIS 
	  cout << " fmin,max "<< freqmin<<","<< freqmax<<" ident "<<ifrmin<<" "<< ifrmax << endl;
	}


	acsum.SetSize(8, NTFCols);
237
	if (params.doSigma_)   // for computing Sigma, if required 
perdereau's avatar
perdereau committed
238
	  acsum_sq.SetSize(8, NTFCols);
239

perdereau's avatar
perdereau committed
240
	cxsum.SetSize(6, NTFCols);
241
	if (params.doSigma_) {   // for computing Sigma, if required 
perdereau's avatar
perdereau committed
242 243
	  cxsum_sq_rp.SetSize(6, NTFCols);
	  cxsum_sq_ip.SetSize(6, NTFCols);
244
	}
245 246
	// VV if needed 
	if (params.doVV_){
perdereau's avatar
perdereau committed
247
	  cxsum_vv.SetSize(6, NTFCols);
248
	  if (params.doSigma_) {   // for computing Sigma, if required 
perdereau's avatar
perdereau committed
249 250
	    cxsum_vv_sq_rp.SetSize(6, NTFCols);
	    cxsum_vv_sq_ip.SetSize(6, NTFCols);
251
	  }
252 253 254 255
	}

	// HV if needed 
	if (params.doHV_){
perdereau's avatar
perdereau committed
256
	  cxsum_hv.SetSize(16, NTFCols);
257
	  if (params.doSigma_) {   // for computing Sigma, if required 
perdereau's avatar
perdereau committed
258 259
	    cxsum_hv_sq_rp.SetSize(16, NTFCols);
	    cxsum_hv_sq_ip.SetSize(16, NTFCols);
260
	  }
261 262 263
	}


perdereau's avatar
perdereau committed
264

OP's avatar
OP committed
265 266
	//allocating 8 Auto-Corr time-frequency maps 
	cout<<"visiavg/Info: allocating 8 AutoCor Time-Frequency maps : Time->NX="<<tfmSX<<" x Freq->NY="<<tfmSY<<endl;
267 268 269 270 271
	for(int k=0; k<8; k++) vtfmac.push_back( TArray< r_4 >(tfmSX, tfmSY) );
	if (params.doSigma_) {   // for computing Sigma, if required 
	  cout <<" and 8 for the the variance maps "<<endl;
	  for(int k=0; k<8; k++) vtfmac_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) ); 
	}
OP's avatar
OP committed
272 273 274
	//allocating 6 Cross-Corr H-H time-frequency maps 
	cout<<"visiavg/Info: allocating H-H cross-cor Time-Frequency maps : Time->NX="<<tfmSX<<" x Freq->NY="<<tfmSY<<endl;
	for(int k=0; k<6; k++) vtfm.push_back( TArray< complex<r_4> >(tfmSX, tfmSY) );
275 276 277 278 279
	if (params.doSigma_) {   // for computing Sigma, if required 
	  cout << "and the 2x6 for squares of real & imaginary parts " << endl;
	  for(int k=0; k<6; k++) vtfm_rp_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
	  for(int k=0; k<6; k++) vtfm_ip_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
	}
280 281 282 283 284
	// VV if needed 
	if (params.doVV_){
	  //allocating 6 Cross-Corr V-V time-frequency maps 
	  cout<<"visiavg/Info: allocating V-V cross-cor Time-Frequency maps : Time->NX="<<tfmSX<<" x Freq->NY="<<tfmSY<<endl;
	  for(int k=0; k<6; k++) vtfm_vv.push_back( TArray< complex<r_4> >(tfmSX, tfmSY) );
285 286 287 288 289
	  if (params.doSigma_) {   // for computing Sigma, if required 
	    cout << "and the 2x6 for squares of real & imaginary parts " << endl;
	    for(int k=0; k<6; k++) vtfm_vv_rp_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
	    for(int k=0; k<6; k++) vtfm_vv_ip_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
	  }
290 291 292 293 294 295 296
	}

	// HV if needed 
	if (params.doHV_){
	  //allocating 16 Cross-Corr H-V time-frequency maps 
	  cout<<"visiavg/Info: allocating 16 H-V cross-cor Time-Frequency maps : Time->NX="<<tfmSX<<" x Freq->NY="<<tfmSY<<endl;
	  for(int k=0; k<16; k++) vtfm_hv.push_back( TArray< complex<r_4> >(tfmSX, tfmSY) );
297 298 299 300 301
	  if (params.doSigma_) {   // for computing Sigma, if required 
	    cout << "and the 2x16 for squares of real & imaginary parts " << endl;
	    for(int k=0; k<16; k++) vtfm_hv_rp_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
	    for(int k=0; k<16; k++) vtfm_hv_ip_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
	  }
302 303
	}

OP's avatar
OP committed
304 305 306

	// recupere le jour de depart @ 0h
	datestart = TimeStamp(cfdate.DaysPart(),0.);
307

308 309
      } // end if cnt==0

OP's avatar
OP committed
310 311 312 313 314 315
      if (I==0) {   // start filling a new time bin 
	dateobs=cfdate;
	if (prtlev>0) 
	  cout<<"visiavg/Info:  dateobs="<<dateobs<<" SecondsPart()="<<dateobs.SecondsPart()<<endl;
	acsum = complex<r_4>(0.,0.);
	cxsum = complex<r_4>(0.,0.);
316 317 318 319 320
	if (params.doSigma_) {   // for computing Sigma, if required 
	  acsum_sq = 0.;
	  cxsum_sq_rp = 0.;
	  cxsum_sq_ip = 0.;
	}
321 322 323
	// VV if needed 
	if (params.doVV_){
	  cxsum_vv = complex<r_4>(0.,0.);
324 325 326 327
	  if (params.doSigma_) {   // for computing Sigma, if required 
	    cxsum_vv_sq_rp = 0.;
	    cxsum_vv_sq_ip = 0.;
	  }
328 329 330 331 332
	}

	// HV if needed 
	if (params.doHV_){
	  cxsum_hv = complex<r_4>(0.,0.);
333 334 335 336
	  if (params.doSigma_) {   // for computing Sigma, if required 
	    cxsum_hv_sq_rp = 0.;
	    cxsum_hv_sq_ip = 0.;
	  }
337
	}
338
      }// end if(I==0)
OP's avatar
OP committed
339 340
      
      //   sum (integration) along the time axis 
perdereau's avatar
perdereau committed
341 342 343 344 345 346 347 348
      for(size_t k=0; k<KVAC.size(); k++){
	//cout << vismtx.Row(KVAC[k]) << endl;
	//cout <<	(vismtx.Row(KVAC[k])).SubMatrix(Range::all(),Range(ifrmin,ifrmax))<< endl;
       	//cout <<acsum.Row(k) << endl;
	acsum.Row(k) += (vismtx.Row(KVAC[k])).SubMatrix( Range::all(),Range(ifrmin,ifrmax)); 

      }
      //cout << " end ac "<<endl;
349 350
      if (params.doSigma_) {   // for computing Sigma, if required 
	for(size_t k=0; k<KVAC.size(); k++)      {
perdereau's avatar
perdereau committed
351
	  TVector<r_4> tmp = real(vismtx.Row(KVAC[k]).SubMatrix(Range::all(),Range(ifrmin,ifrmax))) ; //(Range(ifrmin,ifrmax)));
352 353
	  acsum_sq.Row(k) += tmp.MulElt(tmp,tmp) ;
	}
OP's avatar
OP committed
354
      }
perdereau's avatar
perdereau committed
355
      //cout << " end ac "<<endl;
OP's avatar
OP committed
356
      for(size_t k=0; k<KVCXHH.size(); k++){
perdereau's avatar
perdereau committed
357
	cxsum.Row(k) += vismtx.Row(KVCXHH[k]).SubMatrix(Range::all(),Range(ifrmin,ifrmax)); // (Range(ifrmin,ifrmax));   // les cross-correlations 
358
	if (params.doSigma_) {   // for computing Sigma, if required 
perdereau's avatar
perdereau committed
359
	  TVector<r_4> tmp = real(vismtx.Row(KVCXHH[k]).SubMatrix(Range::all(),Range(ifrmin,ifrmax)));//(Range(ifrmin,ifrmax)));
360
	  cxsum_sq_rp.Row(k) +=  tmp.MulElt(tmp,tmp) ;
perdereau's avatar
perdereau committed
361
	  tmp = imag(vismtx.Row(KVCXHH[k]).SubMatrix(Range::all(),Range(ifrmin,ifrmax)));//(Range(ifrmin,ifrmax)));
362 363
	  cxsum_sq_ip.Row(k) +=  tmp.MulElt(tmp,tmp) ;
	}
OP's avatar
OP committed
364
      }
365 366 367
      
      if (params.doVV_){
	for(size_t k=0; k<KVCXVV.size(); k++){
perdereau's avatar
perdereau committed
368
	  cxsum_vv.Row(k) += vismtx.Row(KVCXVV[k]).SubMatrix(Range::all(),Range(ifrmin,ifrmax)) ;//(Range(ifrmin,ifrmax));   // les cross-correlations 
369
	  if (params.doSigma_) {   // for computing Sigma, if required 
perdereau's avatar
perdereau committed
370
	    TVector<r_4> tmp = real(vismtx.Row(KVCXVV[k]).SubMatrix(Range::all(),Range(ifrmin,ifrmax)));//(Range(ifrmin,ifrmax)));
371
	    cxsum_vv_sq_rp.Row(k) +=  tmp.MulElt(tmp,tmp) ;
perdereau's avatar
perdereau committed
372
	    tmp = imag(vismtx.Row(KVCXVV[k]).SubMatrix(Range::all(),Range(ifrmin,ifrmax)));//(Range(ifrmin,ifrmax)));
373 374
	    cxsum_vv_sq_ip.Row(k) +=  tmp.MulElt(tmp,tmp) ;
	  }
375 376 377 378 379
	}
      }

      if (params.doHV_){
	for(size_t k=0; k<KVCXHV.size(); k++){
perdereau's avatar
perdereau committed
380
	  cxsum_hv.Row(k) += vismtx.Row(KVCXHV[k]).SubMatrix(Range::all(),Range(ifrmin,ifrmax));//(Range(ifrmin,ifrmax));   // les cross-correlations HV
381
	  if (params.doSigma_) {   // for computing Sigma, if required 
perdereau's avatar
perdereau committed
382
	    TVector<r_4> tmp = real(vismtx.Row(KVCXHV[k]).SubMatrix(Range::all(),Range(ifrmin,ifrmax)));//(Range(ifrmin,ifrmax)));
383
	    cxsum_hv_sq_rp.Row(k) +=  tmp.MulElt(tmp,tmp) ; 
perdereau's avatar
perdereau committed
384
	    tmp = imag(vismtx.Row(KVCXHV[k]).SubMatrix(Range::all(),Range(ifrmin,ifrmax)));//(Range(ifrmin,ifrmax)));
385 386
	    cxsum_hv_sq_ip.Row(k) +=  tmp.MulElt(tmp,tmp) ;
	  }
387 388 389
	}
      }

perdereau's avatar
perdereau committed
390
      //cout << " at i++ "<<endl;
OP's avatar
OP committed
391

OP's avatar
OP committed
392
      I++;    // incrementing DeltaTime counter 
393 394 395 396 397 398 399
      // we check that our time index did not go beyond the allocated array size (might not be necessary)
      if ((I==deltaIavg)&&(TFMtmidx>=tfmSX)) {  // Cela ne devrait pas arriver en principe 
	TFMtmidx++;	I=0;  
	cout << "visiavg/Warning: something wrong in the logic , (TFMtmidx="<<TFMtmidx<<") >= (tfmSX="<<tfmSX<<")"
	     << " for read count="<<cnt<<endl; 
      }
      else if (I==deltaIavg) {   // Filling TFM maps 
OP's avatar
OP committed
400 401 402
	//---- On s'occupe d'abord des autocorrelations P1 ... P8 
	for(size_t k=0; k<KVAC.size(); k++) {  // Loop over the 8 auto-correlations 
	  TVector<r_4> vac = real(acsum.Row(k));
403 404 405 406 407 408
	  TArray< r_4 > & tfmap = vtfmac[k];
	  for(sa_size_t jy=0; jy<tfmSY; jy++) {  // frequency binning 
	    tfmap(TFMtmidx, jy) = vac( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
	  }
	  if (params.doSigma_) {   // for computing Sigma, if required 
	    TVector<r_4> vacsq = acsum_sq.Row(k);
OP's avatar
OP committed
409
	    TArray< r_4 > & tfmap_sq = vtfmac_sq[k];
OP's avatar
OP committed
410
	    for(sa_size_t jy=0; jy<tfmSY; jy++) {  // frequency binning 
OP's avatar
OP committed
411
	      tfmap_sq(TFMtmidx, jy) = vacsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
412 413
	    }
	  }  //-- end of computing sigmas   
OP's avatar
OP committed
414
	}  //----- end of loop over the 8 AutoCor
OP's avatar
OP committed
415

OP's avatar
OP committed
416
	//---- On s'occupe des 6 cross-correlations  1H-2H ... 3H-4H 
OP's avatar
OP committed
417 418
 	for(size_t k=0; k<KVCXHH.size(); k++)   {   // loop over the 6 Xcor 	  
	  TVector< complex<r_4> > vcx = cxsum.Row(k);
419 420 421 422 423 424 425
	  TArray< complex<r_4> > & tfmap = vtfm[k];
	  for(sa_size_t jy=0; jy<tfmSY; jy++) {
	    tfmap(TFMtmidx, jy) = vcx( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
	  } 
	  if (params.doSigma_) {   // for computing Sigma, if required 
	    TVector<r_4>  vcxprsq = cxsum_sq_rp.Row(k);
	    TVector<r_4>  vcxpisq = cxsum_sq_ip.Row(k);
OP's avatar
OP committed
426 427
	    TArray< r_4 > & tfmapsqpr = vtfm_rp_sq[k];
	    TArray< r_4 > & tfmapsqpi = vtfm_ip_sq[k];
OP's avatar
OP committed
428
	    for(sa_size_t jy=0; jy<tfmSY; jy++) {
OP's avatar
OP committed
429 430
	      tfmapsqpr(TFMtmidx, jy) = vcxprsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
	      tfmapsqpi(TFMtmidx, jy) = vcxpisq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
OP's avatar
OP committed
431
	    } 
432
	  }    //-- end of computing sigmas  
OP's avatar
OP committed
433
	}  //----- end of loop over the 6 Xcor 
434 435 436 437 438

	if (params.doVV_){ // Option VV 
	  //---- On s'occupe des 6 cross-correlations  1V-2V ... 3V-4V 
	  for(size_t k=0; k<KVCXVV.size(); k++)   {   // loop over the 6 Xcor 	  
	    TVector< complex<r_4> > vcx = cxsum_vv.Row(k);
439 440 441 442 443 444 445
	    TArray< complex<r_4> > & tfmap = vtfm_vv[k];
	    for(sa_size_t jy=0; jy<tfmSY; jy++) {
	      tfmap(TFMtmidx, jy) = vcx( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
	    }
	    if (params.doSigma_) {   // for computing Sigma, if required 
	      TVector<r_4>  vcxprsq = cxsum_vv_sq_rp.Row(k);
	      TVector<r_4>  vcxpisq = cxsum_vv_sq_ip.Row(k);
446 447 448 449 450
	      TArray< r_4 > & tfmapsqpr = vtfm_vv_rp_sq[k];
	      TArray< r_4 > & tfmapsqpi = vtfm_vv_ip_sq[k];
	      for(sa_size_t jy=0; jy<tfmSY; jy++) {
		tfmapsqpr(TFMtmidx, jy) = vcxprsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
		tfmapsqpi(TFMtmidx, jy) = vcxpisq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
451 452
	      }
	    } //-- end of computing sigmas  
453 454 455 456 457 458 459
	  }  //----- end of loop over the 6 Xcor VV
	} // end VV option 

	if (params.doHV_){ // Option HV 
	  //---- On s'occupe des 16 cross-correlations  1H-1V ... 4H-4V 
	  for(size_t k=0; k<KVCXHV.size(); k++)   {   // loop over the 16 Xcor 	  
	    TVector< complex<r_4> > vcx = cxsum_hv.Row(k);
460 461 462 463 464 465 466
	    TArray< complex<r_4> > & tfmap = vtfm_hv[k];
	    for(sa_size_t jy=0; jy<tfmSY; jy++) {
	      tfmap(TFMtmidx, jy) = vcx( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
	    } 
	    if (params.doSigma_) {   // for computing Sigma, if required 
	      TVector<r_4>  vcxprsq = cxsum_hv_sq_rp.Row(k);
	      TVector<r_4>  vcxpisq = cxsum_hv_sq_ip.Row(k);
467 468 469 470 471 472
	      TArray< r_4 > & tfmapsqpr = vtfm_hv_rp_sq[k];
	      TArray< r_4 > & tfmapsqpi = vtfm_hv_ip_sq[k];
	      for(sa_size_t jy=0; jy<tfmSY; jy++) {
		tfmapsqpr(TFMtmidx, jy) = vcxprsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
		tfmapsqpi(TFMtmidx, jy) = vcxpisq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
	      } 
473
	    } //-- end of computing sigmas  
474 475 476
	  }  //----- end of loop over the 16 Xcor HV
	} // end HV option 

477 478 479 480 481 482 483
	double tdif =  cfdate.TimeDifferenceSeconds(cfdate,dateobs)/2.; // 1/2 largeur du bin en temps 
	double mytime = dateobs.TimeDifferenceSeconds(dateobs.ShiftSeconds (tdif ),datestart);	// centre du bin 

	if ( mytime == 0.) {
	  cout<<"visiavg/PROBLEME : time=0 ?!? " <<cntnt<<endl; 
	}
	timevec(TFMtmidx) = mytime ;
484
	ravec(TFMtmidx) = P4Coords::RAFromTimeTU(dateobs.ShiftSeconds (tdif ));
OP's avatar
OP committed
485 486 487 488 489
	TFMtmidx++;
	//  ... done 
	I=0;  cntnt++;
      }
      cnt++;
490 491 492
      if ((cnt>0)&&(cntnt%10==0)&&(cntnt>pcntnt)) 
	{
	cout<<"visiavg/Info: TFM-Map fill cnt="<<cntnt<<" VisMtxCount="<<cnt<<" tfmidx "<< TFMtmidx
493
	    <<" /Max="<<wreader.getTotalNbWindows()<<" DateObs="<<dateobs<<endl;
OP's avatar
OP committed
494
	pcntnt=cntnt;
495 496 497 498
	}

    }// end while(fgok)

499
    cout<<"visiavg/Info: count="<<cnt*wreader.getWindowSize()<<" Visibility Matrices read "<<endl;
OP's avatar
OP committed
500

OP's avatar
OP committed
501
    // --- Sauvegarde cartes temps-frequence en fits 
perdereau's avatar
perdereau committed
502

OP's avatar
OP committed
503 504
    FitsInOutFile  * fos = NULL ;

OP's avatar
OP committed
505
    if (fitsoutfile.length()>=1){
506
      cout << " fitsoutfile :" <<fitsoutfile<<":"<< endl;
OP's avatar
OP committed
507 508 509
      fos = new FitsInOutFile(fitsoutfile, FitsInOutFile::Fits_Create);
    }
    
OP's avatar
OP committed
510

OP's avatar
OP committed
511
    POutPersist potfm(outfile);
OP's avatar
OP committed
512 513
    char bufnam[10];
    int numkey =0;
OP's avatar
OP committed
514 515 516
    // --- renormalizing and saving AutoCorr time-frequency maps 
    cout<<"  visiavg/Info: Saving 8 AutoCorr time-frequency maps to PPF file "<<outfile<<endl;
    const char* tfm_names[8]={"TFM_1H", "TFM_2H", "TFM_3H", "TFM_4H", "TFM_1V", "TFM_2V", "TFM_3V", "TFM_4V"};
OP's avatar
OP committed
517
    const char* tfmsq_names[8]={"VARTFM_1H", "VARTFM_2H", "VARTFM_3H", "VARTFM_4H", "VARTFM_1V", "VARTFM_2V", "VARTFM_3V", "VARTFM_4V"};
OP's avatar
OP committed
518 519 520 521
    for(int k=0; k<8; k++)  {   // loop over the 8 AutoCorr 
      TArray< r_4 > & tfmap = vtfmac[k];
      tfmap *= (r_4)(1./((double)deltaIavg*(double)TFMfbin));
      potfm << PPFNameTag(tfm_names[k]) << tfmap;
OP's avatar
OP committed
522

OP's avatar
OP committed
523 524 525
      if (fos != NULL) {
	ext_names.push_back(tfm_names[k]);
	(*fos)<<  tfmap;
OP's avatar
OP committed
526
      }
527 528 529 530 531 532 533 534 535 536
      if (params.doSigma_) {   // for saving TFM-Sigma, if required 
	TArray< r_4 > & tfmapsq = vtfmac_sq[k];
	tfmapsq *= (r_4)(1./((double)deltaIavg*(double)TFMfbin));
	tfmapsq = tfmapsq - tfmap.MulElt(tfmap,tfmap) ; // tfmap ->tfmap*tfmap !! 
	potfm << PPFNameTag(tfmsq_names[k]) << tfmapsq;
	if (fos != NULL) {
	  ext_names.push_back(tfmsq_names[k]);
	  (*fos)<<  tfmapsq;
	}
      } // end of saving TFM-of-sigma 
OP's avatar
OP committed
537
    }
OP's avatar
OP committed
538

OP's avatar
OP committed
539 540
    // --- renormalizing and saving H-H Cross-Corr time-frequency maps 
    cout<<"  visiavg/Info: Saving 6 H-H cross-corr time-frequency maps to PPF file "<<outfile<<endl;
OP's avatar
OP committed
541

OP's avatar
OP committed
542
    const char* tfmCC_names[6]={"TFM_1H2H", "TFM_1H3H", "TFM_1H4H", "TFM_2H3H", "TFM_2H4H", "TFM_3H4H"};
OP's avatar
OP committed
543 544 545
    const char* vrtfmCC_names[6]={"RVARTFM_1H2H", "RVARTFM_1H3H", "RVARTFM_1H4H", "RVARTFM_2H3H", "RVARTFM_2H4H", "RVARTFM_3H4H"};
    const char* vitfmCC_names[6]={"IVARTFM_1H2H", "IVARTFM_1H3H", "IVARTFM_1H4H", "IVARTFM_2H3H", "IVARTFM_2H4H", "IVARTFM_3H4H"};

OP's avatar
OP committed
546 547
    for(int k=0; k<6; k++)  {   // loop over the 6 Xcor 
      TArray< complex<r_4> > & tfmap = vtfm[k];
OP's avatar
OP committed
548 549 550

      TArray< r_4 > & tfmap_sqpr = vtfm_rp_sq[k];
      TArray< r_4 > & tfmap_sqpi = vtfm_ip_sq[k];
OP's avatar
OP committed
551
      tfmap *= complex<r_4>((r_4)(1./((double)deltaIavg*(double)TFMfbin)), 0.);
OP's avatar
OP committed
552 553 554 555 556 557 558
      potfm << PPFNameTag(tfmCC_names[k]) << tfmap;
      if (fos != NULL) {
	ext_names.push_back(string(tfmCC_names[k])+"_real");
	(*fos)<<  real(tfmap);
	ext_names.push_back(string(tfmCC_names[k])+"_imag");
	(*fos)<<  imag(tfmap);
      }
559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577
      if (params.doSigma_) {   // for saving TFM-Sigma, if required 
	TArray< r_4 >  tfmapr = real(tfmap);
	TArray< r_4 >  tfmapi = imag(tfmap);
	
	tfmap_sqpr *= ((r_4)(1./((double)deltaIavg*(double)TFMfbin)));
	tfmap_sqpi *= ((r_4)(1./((double)deltaIavg*(double)TFMfbin)));
	tfmapr = tfmapr.MulElt(tfmapr,tfmapr) ;
	tfmap_sqpr -= tfmapr ;
	tfmapi = tfmapi.MulElt(tfmapi,tfmapi) ;
	tfmap_sqpi -= tfmapi;
	potfm << PPFNameTag(vrtfmCC_names[k]) << tfmap_sqpr;
	potfm << PPFNameTag(vitfmCC_names[k]) << tfmap_sqpi;
	if (fos != NULL) {
	  ext_names.push_back(vrtfmCC_names[k]);
	  (*fos)<<  tfmap_sqpr;
	  ext_names.push_back(vitfmCC_names[k]);
	  (*fos)<<  tfmap_sqpi;
	}
      }  // end of saving TFM-of-sigma  
OP's avatar
OP committed
578
    }
579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603
    // ========================== end saving HH XCorr

    //===========================
    if (params.doVV_){ // Option VV 
      // --- renormalizing and saving V-V Cross-Corr time-frequency maps 
      cout<<"  visiavg/Info: Saving 6 V-V cross-corr time-frequency maps to PPF file "<<outfile<<endl;

      const char* tfmVV_names[6]={"TFM_1V2V", "TFM_1V3V", "TFM_1V4V", "TFM_2V3V", "TFM_2V4V", "TFM_3V4V"};
      const char* vrtfmVV_names[6]={"RVARTFM_1V2V", "RVARTFM_1V3V", "RVARTFM_1V4V", "RVARTFM_2V3V", "RVARTFM_2V4V", "RVARTFM_3V4V"};
      const char* vitfmVV_names[6]={"IVARTFM_1V2V", "IVARTFM_1V3V", "IVARTFM_1V4V", "IVARTFM_2V3V", "IVARTFM_2V4V", "IVARTFM_3V4V"};
      
      for(int k=0; k<6; k++)  {   // loop over the 6 Xcor 
	TArray< complex<r_4> > & tfmap = vtfm_vv[k];
	
	TArray< r_4 > & tfmap_sqpr = vtfm_vv_rp_sq[k];
	TArray< r_4 > & tfmap_sqpi = vtfm_vv_ip_sq[k];
	tfmap *= complex<r_4>((r_4)(1./((double)deltaIavg*(double)TFMfbin)), 0.);
	potfm << PPFNameTag(tfmVV_names[k]) << tfmap;

	if (fos != NULL) {
	  ext_names.push_back(string(tfmVV_names[k])+"_real");
	  (*fos)<<  real(tfmap);
	  ext_names.push_back(string(tfmVV_names[k])+"_imag");
	  (*fos)<<  imag(tfmap);
	}
604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622
	if (params.doSigma_) {   // for saving TFM-Sigma, if required 
	  TArray< r_4 >  tfmapr = real(tfmap);
	  TArray< r_4 >  tfmapi = imag(tfmap);
	  
	  tfmap_sqpr *= ((r_4)(1./((double)deltaIavg*(double)TFMfbin)));
	  tfmap_sqpi *= ((r_4)(1./((double)deltaIavg*(double)TFMfbin)));
	  tfmapr = tfmapr.MulElt(tfmapr,tfmapr) ;
	  tfmap_sqpr -= tfmapr ;
	  tfmapi = tfmapi.MulElt(tfmapi,tfmapi) ;
	  tfmap_sqpi -= tfmapi;
	  potfm << PPFNameTag(vrtfmVV_names[k]) << tfmap_sqpr;
	  potfm << PPFNameTag(vitfmVV_names[k]) << tfmap_sqpi;
	  if (fos != NULL) {
	    ext_names.push_back(vrtfmVV_names[k]);
	    (*fos)<<  tfmap_sqpr;
	    ext_names.push_back(vitfmVV_names[k]);
	    (*fos)<<  tfmap_sqpi;
	  }
	}   // end of saving TFM-of-sigma  
623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661
      }
    } // end option VV 


    //===========================
    if (params.doHV_){ // Option HV 
      // --- renormalizing and saving H-V Cross-Corr time-frequency maps 
      cout<<"  visiavg/Info: Saving 16 H-V cross-corr time-frequency maps to PPF file "<<outfile<<endl;

      const char* tfmHV_names[16]={"TFM_1H1V", "TFM_1H2V", "TFM_1H3V", "TFM_1H4V", 
				   "TFM_2H1V", "TFM_2H2V", "TFM_2H3V", "TFM_2H4V", 
				   "TFM_3H1V", "TFM_3H2V", "TFM_3H3V", "TFM_3H4V",  
				   "TFM_4H1V", "TFM_4H2V", "TFM_4H3V", "TFM_4H4V"
      };
      const char* vrtfmHV_names[16]={"RVARTFM_1H1V", "RVARTFM_1H2V", "RVARTFM_1H3V", "RVARTFM_1H4V",
				     "RVARTFM_2H1V", "RVARTFM_2H2V", "RVARTFM_2H3V", "RVARTFM_2H4V",
				     "RVARTFM_3H1V", "RVARTFM_3H2V", "RVARTFM_3H3V", "RVARTFM_3H4V",
				     "RVARTFM_4H1V", "RVARTFM_4H2V", "RVARTFM_4H3V", "RVARTFM_4H4V"
      };
      const char* vitfmHV_names[16]={"IVARTFM_1H1V", "IVARTFM_1H2V", "IVARTFM_1H3V", "IVARTFM_1H4V", 
				     "IVARTFM_2H1V", "IVARTFM_2H2V", "IVARTFM_2H3V", "IVARTFM_2H4V", 
				     "IVARTFM_3H1V", "IVARTFM_3H2V", "IVARTFM_3H3V", "IVARTFM_3H4V", 
				     "IVARTFM_4H1V", "IVARTFM_4H2V", "IVARTFM_4H3V", "IVARTFM_4H4V" 
      };
      
      for(int k=0; k<16; k++)  {   // loop over the 6 Xcor 
	TArray< complex<r_4> > & tfmap = vtfm_hv[k];
	
	TArray< r_4 > & tfmap_sqpr = vtfm_hv_rp_sq[k];
	TArray< r_4 > & tfmap_sqpi = vtfm_hv_ip_sq[k];
	tfmap *= complex<r_4>((r_4)(1./((double)deltaIavg*(double)TFMfbin)), 0.);
	potfm << PPFNameTag(tfmHV_names[k]) << tfmap;

	if (fos != NULL) {
	  ext_names.push_back(string(tfmHV_names[k])+"_real");
	  (*fos)<<  real(tfmap);
	  ext_names.push_back(string(tfmHV_names[k])+"_imag");
	  (*fos)<<  imag(tfmap);
	}
662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681

	if (params.doSigma_) {   // for saving TFM-Sigma, if required 
	  TArray< r_4 >  tfmapr = real(tfmap);
	  TArray< r_4 >  tfmapi = imag(tfmap);
	  
	  tfmap_sqpr *= ((r_4)(1./((double)deltaIavg*(double)TFMfbin)));
	  tfmap_sqpi *= ((r_4)(1./((double)deltaIavg*(double)TFMfbin)));
	  tfmapr = tfmapr.MulElt(tfmapr,tfmapr) ;
	  tfmap_sqpr -= tfmapr ;
	  tfmapi = tfmapi.MulElt(tfmapi,tfmapi) ;
	  tfmap_sqpi -= tfmapi;
	  potfm << PPFNameTag(vrtfmHV_names[k]) << tfmap_sqpr;
	  potfm << PPFNameTag(vitfmHV_names[k]) << tfmap_sqpi;
	  if (fos != NULL) {
	    ext_names.push_back(vrtfmHV_names[k]);
	    (*fos)<<  tfmap_sqpr;
	    ext_names.push_back(vitfmHV_names[k]);
	    (*fos)<<  tfmap_sqpi;
	  }
	}  // end of saving TFM-of-sigma 
682 683 684
      }

    } // end option HV 
perdereau's avatar
perdereau committed
685

OP's avatar
OP committed
686
    P4FreqBand myp4fre;
perdereau's avatar
perdereau committed
687 688
    if(params.fbands_.size()>0 )
      myp4fre=params.fbands_[0];
OP's avatar
OP committed
689 690 691 692 693 694
    TVector <double> lim_freq(2);
    if (params.gain_gnu_file_.length()>0) {
      P4gnuGain p4gnu( params.gain_gnu_file_ );
      lim_freq(0) =  p4gnu. minGoodF();
      lim_freq(1) =  p4gnu. maxGoodF();
    }else{
perdereau's avatar
perdereau committed
695 696 697
      lim_freq(0) = myp4fre.getP4Frequency(params.fbands_[0].jfmin_);
      lim_freq(1) = myp4fre.getP4Frequency(params.fbands_[0].jfmax_);
      
OP's avatar
OP committed
698
    }
699 700

    cout  << " frequences limites "<< lim_freq(0) <<" ; "<< lim_freq(1)<<endl;
OP's avatar
OP committed
701
    potfm << PPFNameTag("FreqLims") << lim_freq ;
OP's avatar
OP committed
702
    potfm << PPFNameTag("TimeVec") << timevec ;
OP's avatar
OP committed
703
    potfm << PPFNameTag("RAVec") << ravec ;
OP's avatar
OP committed
704 705
    // --- FIN sauvegarde cartes temps-frequence 
    //    resu.Update();
OP's avatar
OP committed
706 707

    
perdereau's avatar
perdereau committed
708 709 710 711 712
    size_t szFreqs = myp4fre.getP4NbFreqChannels()/TFMfbin;
    if (params.fbands_.size()>0 )
      szFreqs = (params.fbands_[0].jfmax_ - params.fbands_[0].jfmin_)/TFMfbin;

    TVector <double> avg_freqs( szFreqs );
OP's avatar
OP committed
713

714
    
perdereau's avatar
perdereau committed
715 716
    double frbase =  myp4fre.getP4Frequency(params.fbands_[0].jfmin_) + myp4fre.getP4FreqResolution()/2. ;
    for (int kf=0 ; kf< szFreqs ; kf++,frbase += myp4fre.getP4FreqResolution()*TFMfbin )
OP's avatar
OP committed
717
      avg_freqs(kf) = frbase ;
718
    
OP's avatar
OP committed
719
    potfm << PPFNameTag("FreqVec") << avg_freqs;
OP's avatar
OP committed
720

OP's avatar
OP committed
721
    if (fos != NULL) {
OP's avatar
OP committed
722 723
      ext_names.push_back("FreqsLims");
      (*fos)<< lim_freq ;
OP's avatar
OP committed
724 725 726 727 728 729 730 731 732 733
      ext_names.push_back("Frequences");
      (*fos)<< avg_freqs ;
      ext_names.push_back("RAs");
      (*fos)<<  ravec;
      ext_names.push_back("Times");
      (*fos)<<  timevec;
      cout << " number of objs in fits "<< ext_names.size() << endl;
      cout << ext_names << endl;
      delete(fos);
    }
734
    cout << " return code "<<rc<<endl;
OP's avatar
OP committed
735
    cout << resu;   // Update est fait lors du print
736

OP's avatar
OP committed
737 738 739 740 741 742 743 744 745 746 747 748 749 750
  }
  catch (PException& exc) {
    cerr << " visiavg.cc catched PException " << exc.Msg() << endl;
    rc = 77;
  }  
  catch (std::exception& sex) {
    cerr << "\n visiavg.cc std::exception :" 
         << (string)typeid(sex).name() << "\n msg= " 
         << sex.what() << endl;
    rc = 78;
  }
  catch (...) {
    cerr << " visiavg.cc catched unknown (...) exception  " << endl; 
    rc = 79; 
751
  }
OP's avatar
OP committed
752 753 754 755 756 757

  cout << ">>>> visiavg.cc ------- END ----------- RC=" << rc << endl;
  return rc;

}