visiavg.cc 26.5 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"
OP's avatar
OP committed
90 91 92 93 94 95 96 97 98 99
       <<"outfile="<<outfile<<" PrtLev="<<prtlev<<endl;

  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();
100 101
  vector<sa_size_t> KVCXVV = visiencod.getAllVCrossCor();
  vector<sa_size_t> KVCXHV = visiencod.getAllHVCrossCor();
102 103

  
OP's avatar
OP committed
104 105 106 107
  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;
  }
108

OP's avatar
OP committed
109 110 111 112
  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
113

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

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

122
    VisiP4WindowReader wreader(params);
123

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

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

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

146 147 148 149 150 151 152 153 154 155 156 157
    // 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;

158
    TimeStamp dateobs, cfdate,datestart,datused;
OP's avatar
OP committed
159 160 161 162
    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
163
    // for
164

OP's avatar
OP committed
165 166
    //----- 6 H-H cross-cor TimeFrequency maps 
    vector< TArray< complex<r_4> > > vtfm;
OP's avatar
OP committed
167 168 169 170
    //----- 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;

171 172 173 174 175 176 177 178 179 180 181 182
    //----- 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
183 184
    //----- 8 auto-corr TimeFrequency maps 
    vector< TArray< r_4 > > vtfmac;
OP's avatar
OP committed
185 186
    //----- 8 auto-corr TimeFrequency variance maps 
    vector< TArray< r_4 > > vtfmac_sq;
OP's avatar
OP committed
187 188 189 190
    
    //---- for the time-freqency map filling    
    sa_size_t TFMtmidx=0;
    sa_size_t tfmSX, tfmSY;
191 192
    long serialnum; 
    size_t rdidx, lastrdidx=0;
OP's avatar
OP committed
193
    while (fgok) {
194
      //reads next visibility matrix window 
195
      fgok = wreader.Shift(rdidx, serialnum);
196
      
197
      if (!fgok)  break;
198 199 200 201 202 203 204
      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;

205
      vismtx = wreader.getAverageVisMtx(cfdate);
206
      
207
      if (cnt==0)  {    //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations 
208

OP's avatar
OP committed
209
	acsum.SetSize(8, vismtx.NCols());
210 211
	if (params.doSigma_)   // for computing Sigma, if required 
	  acsum_sq.SetSize(8, vismtx.NCols());
212

OP's avatar
OP committed
213
	cxsum.SetSize(6, vismtx.NCols());
214 215 216 217
	if (params.doSigma_) {   // for computing Sigma, if required 
	  cxsum_sq_rp.SetSize(6, vismtx.NCols());
	  cxsum_sq_ip.SetSize(6, vismtx.NCols());
	}
218 219 220
	// VV if needed 
	if (params.doVV_){
	  cxsum_vv.SetSize(6, vismtx.NCols());
221 222 223 224
	  if (params.doSigma_) {   // for computing Sigma, if required 
	    cxsum_vv_sq_rp.SetSize(6, vismtx.NCols());
	    cxsum_vv_sq_ip.SetSize(6, vismtx.NCols());
	  }
225 226 227 228 229
	}

	// HV if needed 
	if (params.doHV_){
	  cxsum_hv.SetSize(16, vismtx.NCols());
230 231 232 233
	  if (params.doSigma_) {   // for computing Sigma, if required 
	    cxsum_hv_sq_rp.SetSize(16, vismtx.NCols());
	    cxsum_hv_sq_ip.SetSize(16, vismtx.NCols());
	  }
234 235 236
	}


237
	tfmSX=wreader.getTotalNbWindows()/deltaIavg;
OP's avatar
OP committed
238 239 240
	tfmSY=vismtx.NCols()/TFMfbin;
	//allocating 8 Auto-Corr time-frequency maps 
	cout<<"visiavg/Info: allocating 8 AutoCor Time-Frequency maps : Time->NX="<<tfmSX<<" x Freq->NY="<<tfmSY<<endl;
241 242 243 244 245
	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
246 247 248
	//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) );
249 250 251 252 253
	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) );
	}
254 255 256 257 258
	// 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) );
259 260 261 262 263
	  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) );
	  }
264 265 266 267 268 269 270
	}

	// 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) );
271 272 273 274 275
	  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) );
	  }
276 277
	}

OP's avatar
OP committed
278 279 280

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

282 283
      } // end if cnt==0

OP's avatar
OP committed
284 285 286 287 288 289
      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.);
290 291 292 293 294
	if (params.doSigma_) {   // for computing Sigma, if required 
	  acsum_sq = 0.;
	  cxsum_sq_rp = 0.;
	  cxsum_sq_ip = 0.;
	}
295 296 297
	// VV if needed 
	if (params.doVV_){
	  cxsum_vv = complex<r_4>(0.,0.);
298 299 300 301
	  if (params.doSigma_) {   // for computing Sigma, if required 
	    cxsum_vv_sq_rp = 0.;
	    cxsum_vv_sq_ip = 0.;
	  }
302 303 304 305 306
	}

	// HV if needed 
	if (params.doHV_){
	  cxsum_hv = complex<r_4>(0.,0.);
307 308 309 310
	  if (params.doSigma_) {   // for computing Sigma, if required 
	    cxsum_hv_sq_rp = 0.;
	    cxsum_hv_sq_ip = 0.;
	  }
311
	}
312
      }// end if(I==0)
OP's avatar
OP committed
313 314 315
      
      //   sum (integration) along the time axis 
      for(size_t k=0; k<KVAC.size(); k++)      acsum.Row(k) += vismtx.Row(KVAC[k]);     // Les auto-correlations 
316 317 318 319 320
      if (params.doSigma_) {   // for computing Sigma, if required 
	for(size_t k=0; k<KVAC.size(); k++)      {
	  TVector<r_4> tmp = real(vismtx.Row(KVAC[k]));
	  acsum_sq.Row(k) += tmp.MulElt(tmp,tmp) ;
	}
OP's avatar
OP committed
321 322 323 324
      }

      for(size_t k=0; k<KVCXHH.size(); k++){
	cxsum.Row(k) += vismtx.Row(KVCXHH[k]);   // les cross-correlations 
325 326 327 328 329 330
	if (params.doSigma_) {   // for computing Sigma, if required 
	  TVector<r_4> tmp = real(vismtx.Row(KVCXHH[k]));
	  cxsum_sq_rp.Row(k) +=  tmp.MulElt(tmp,tmp) ;
	  tmp = imag(vismtx.Row(KVCXHH[k]));
	  cxsum_sq_ip.Row(k) +=  tmp.MulElt(tmp,tmp) ;
	}
OP's avatar
OP committed
331
      }
332 333 334 335
      
      if (params.doVV_){
	for(size_t k=0; k<KVCXVV.size(); k++){
	  cxsum_vv.Row(k) += vismtx.Row(KVCXVV[k]);   // les cross-correlations 
336 337 338 339 340 341
	  if (params.doSigma_) {   // for computing Sigma, if required 
	    TVector<r_4> tmp = real(vismtx.Row(KVCXVV[k]));
	    cxsum_vv_sq_rp.Row(k) +=  tmp.MulElt(tmp,tmp) ;
	    tmp = imag(vismtx.Row(KVCXVV[k]));
	    cxsum_vv_sq_ip.Row(k) +=  tmp.MulElt(tmp,tmp) ;
	  }
342 343 344 345 346 347
	}
      }

      if (params.doHV_){
	for(size_t k=0; k<KVCXHV.size(); k++){
	  cxsum_hv.Row(k) += vismtx.Row(KVCXHV[k]);   // les cross-correlations HV
348 349 350 351 352 353
	  if (params.doSigma_) {   // for computing Sigma, if required 
	    TVector<r_4> tmp = real(vismtx.Row(KVCXHV[k]));
	    cxsum_hv_sq_rp.Row(k) +=  tmp.MulElt(tmp,tmp) ; 
	    tmp = imag(vismtx.Row(KVCXHV[k]));
	    cxsum_hv_sq_ip.Row(k) +=  tmp.MulElt(tmp,tmp) ;
	  }
354 355 356
	}
      }

OP's avatar
OP committed
357 358


OP's avatar
OP committed
359
      I++;    // incrementing DeltaTime counter 
360 361 362 363 364 365 366
      // 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
367 368 369
	//---- 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));
370 371 372 373 374 375
	  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
376
	    TArray< r_4 > & tfmap_sq = vtfmac_sq[k];
OP's avatar
OP committed
377
	    for(sa_size_t jy=0; jy<tfmSY; jy++) {  // frequency binning 
OP's avatar
OP committed
378
	      tfmap_sq(TFMtmidx, jy) = vacsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
379 380
	    }
	  }  //-- end of computing sigmas   
OP's avatar
OP committed
381
	}  //----- end of loop over the 8 AutoCor
OP's avatar
OP committed
382

OP's avatar
OP committed
383
	//---- On s'occupe des 6 cross-correlations  1H-2H ... 3H-4H 
OP's avatar
OP committed
384 385
 	for(size_t k=0; k<KVCXHH.size(); k++)   {   // loop over the 6 Xcor 	  
	  TVector< complex<r_4> > vcx = cxsum.Row(k);
386 387 388 389 390 391 392
	  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
393 394
	    TArray< r_4 > & tfmapsqpr = vtfm_rp_sq[k];
	    TArray< r_4 > & tfmapsqpi = vtfm_ip_sq[k];
OP's avatar
OP committed
395
	    for(sa_size_t jy=0; jy<tfmSY; jy++) {
OP's avatar
OP committed
396 397
	      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
398
	    } 
399
	  }    //-- end of computing sigmas  
OP's avatar
OP committed
400
	}  //----- end of loop over the 6 Xcor 
401 402 403 404 405

	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);
406 407 408 409 410 411 412
	    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);
413 414 415 416 417
	      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();
418 419
	      }
	    } //-- end of computing sigmas  
420 421 422 423 424 425 426
	  }  //----- 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);
427 428 429 430 431 432 433
	    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);
434 435 436 437 438 439
	      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();
	      } 
440
	    } //-- end of computing sigmas  
441 442 443
	  }  //----- end of loop over the 16 Xcor HV
	} // end HV option 

444 445 446 447 448 449 450
	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 ;
451
	ravec(TFMtmidx) = P4Coords::RAFromTimeTU(dateobs.ShiftSeconds (tdif ));
OP's avatar
OP committed
452 453 454 455 456
	TFMtmidx++;
	//  ... done 
	I=0;  cntnt++;
      }
      cnt++;
457 458 459
      if ((cnt>0)&&(cntnt%10==0)&&(cntnt>pcntnt)) 
	{
	cout<<"visiavg/Info: TFM-Map fill cnt="<<cntnt<<" VisMtxCount="<<cnt<<" tfmidx "<< TFMtmidx
460
	    <<" /Max="<<wreader.getTotalNbWindows()<<" DateObs="<<dateobs<<endl;
OP's avatar
OP committed
461
	pcntnt=cntnt;
462 463 464 465
	}

    }// end while(fgok)

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

OP's avatar
OP committed
468
    // --- Sauvegarde cartes temps-frequence en fits 
perdereau's avatar
perdereau committed
469

OP's avatar
OP committed
470 471
    FitsInOutFile  * fos = NULL ;

OP's avatar
OP committed
472
    if (fitsoutfile.length()>=1){
473
      cout << " fitsoutfile :" <<fitsoutfile<<":"<< endl;
OP's avatar
OP committed
474 475 476
      fos = new FitsInOutFile(fitsoutfile, FitsInOutFile::Fits_Create);
    }
    
OP's avatar
OP committed
477

OP's avatar
OP committed
478
    POutPersist potfm(outfile);
OP's avatar
OP committed
479 480
    char bufnam[10];
    int numkey =0;
OP's avatar
OP committed
481 482 483
    // --- 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
484
    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
485 486 487 488
    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
489

OP's avatar
OP committed
490 491 492
      if (fos != NULL) {
	ext_names.push_back(tfm_names[k]);
	(*fos)<<  tfmap;
OP's avatar
OP committed
493
      }
494 495 496 497 498 499 500 501 502 503
      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
504
    }
OP's avatar
OP committed
505

OP's avatar
OP committed
506 507
    // --- 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
508

OP's avatar
OP committed
509
    const char* tfmCC_names[6]={"TFM_1H2H", "TFM_1H3H", "TFM_1H4H", "TFM_2H3H", "TFM_2H4H", "TFM_3H4H"};
OP's avatar
OP committed
510 511 512
    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
513 514
    for(int k=0; k<6; k++)  {   // loop over the 6 Xcor 
      TArray< complex<r_4> > & tfmap = vtfm[k];
OP's avatar
OP committed
515 516 517

      TArray< r_4 > & tfmap_sqpr = vtfm_rp_sq[k];
      TArray< r_4 > & tfmap_sqpi = vtfm_ip_sq[k];
OP's avatar
OP committed
518
      tfmap *= complex<r_4>((r_4)(1./((double)deltaIavg*(double)TFMfbin)), 0.);
OP's avatar
OP committed
519 520 521 522 523 524 525
      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);
      }
526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544
      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
545
    }
546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570
    // ========================== 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);
	}
571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589
	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  
590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628
      }
    } // 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);
	}
629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648

	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 
649 650 651
      }

    } // end option HV 
OP's avatar
OP committed
652 653 654 655 656 657 658 659 660 661
    P4FreqBand myp4fre;
    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{
      lim_freq(0) = myp4fre.freqstart_;
      lim_freq(1) = myp4fre.freqend_;
    }
662 663

    cout  << " frequences limites "<< lim_freq(0) <<" ; "<< lim_freq(1)<<endl;
OP's avatar
OP committed
664
    potfm << PPFNameTag("FreqLims") << lim_freq ;
OP's avatar
OP committed
665
    potfm << PPFNameTag("TimeVec") << timevec ;
OP's avatar
OP committed
666
    potfm << PPFNameTag("RAVec") << ravec ;
OP's avatar
OP committed
667 668
    // --- FIN sauvegarde cartes temps-frequence 
    //    resu.Update();
OP's avatar
OP committed
669 670 671

    

OP's avatar
OP committed
672
    TVector <double> avg_freqs( myp4fre.getP4NbFreqChannels()/TFMfbin);
673
    
OP's avatar
OP committed
674
    double frbase =  myp4fre.freqstart_ + myp4fre.getP4FreqResolution()/2. ;
OP's avatar
OP committed
675
    for (int kf=0 ; kf< myp4fre.getP4NbFreqChannels()/TFMfbin ; kf++,frbase += myp4fre.getP4FreqResolution()*TFMfbin )
OP's avatar
OP committed
676
      avg_freqs(kf) = frbase ;
677
    
OP's avatar
OP committed
678
    potfm << PPFNameTag("FreqVec") << avg_freqs;
OP's avatar
OP committed
679

OP's avatar
OP committed
680
    if (fos != NULL) {
OP's avatar
OP committed
681 682
      ext_names.push_back("FreqsLims");
      (*fos)<< lim_freq ;
OP's avatar
OP committed
683 684 685 686 687 688 689 690 691 692
      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);
    }
693
    cout << " return code "<<rc<<endl;
OP's avatar
OP committed
694
    cout << resu;   // Update est fait lors du print
695

OP's avatar
OP committed
696 697 698 699 700 701 702 703 704 705 706 707 708 709
  }
  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; 
710
  }
OP's avatar
OP committed
711 712 713 714 715 716

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

}