visiavg.cc 25.9 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;

OP's avatar
OP committed
158 159 160 161 162
    TimeStamp dateobs, cfdate,datestart;
    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 191 192
    
    //---- for the time-freqency map filling    
    sa_size_t TFMtmidx=0;
    sa_size_t tfmSX, tfmSY;

    while (fgok) {
193
      //reads next visibility matrix window 
194 195
      fgok = wreader.Shift();
      
196 197
      if (!fgok)  break;
     
198
      vismtx = wreader.getAverageVisMtx(cfdate);
199
      
200
      if (cnt==0)  {    //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations 
201

OP's avatar
OP committed
202
	acsum.SetSize(8, vismtx.NCols());
203 204
	if (params.doSigma_)   // for computing Sigma, if required 
	  acsum_sq.SetSize(8, vismtx.NCols());
205

OP's avatar
OP committed
206
	cxsum.SetSize(6, vismtx.NCols());
207 208 209 210
	if (params.doSigma_) {   // for computing Sigma, if required 
	  cxsum_sq_rp.SetSize(6, vismtx.NCols());
	  cxsum_sq_ip.SetSize(6, vismtx.NCols());
	}
211 212 213
	// VV if needed 
	if (params.doVV_){
	  cxsum_vv.SetSize(6, vismtx.NCols());
214 215 216 217
	  if (params.doSigma_) {   // for computing Sigma, if required 
	    cxsum_vv_sq_rp.SetSize(6, vismtx.NCols());
	    cxsum_vv_sq_ip.SetSize(6, vismtx.NCols());
	  }
218 219 220 221 222
	}

	// HV if needed 
	if (params.doHV_){
	  cxsum_hv.SetSize(16, vismtx.NCols());
223 224 225 226
	  if (params.doSigma_) {   // for computing Sigma, if required 
	    cxsum_hv_sq_rp.SetSize(16, vismtx.NCols());
	    cxsum_hv_sq_ip.SetSize(16, vismtx.NCols());
	  }
227 228 229
	}


230
	tfmSX=wreader.getTotalNbWindows()/deltaIavg;
OP's avatar
OP committed
231 232 233
	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;
234 235 236 237 238
	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
239 240 241
	//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) );
242 243 244 245 246
	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) );
	}
247 248 249 250 251
	// 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) );
252 253 254 255 256
	  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) );
	  }
257 258 259 260 261 262 263
	}

	// 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) );
264 265 266 267 268
	  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) );
	  }
269 270
	}

OP's avatar
OP committed
271 272 273 274

	// recupere le jour de depart @ 0h
	datestart = TimeStamp(cfdate.DaysPart(),0.);
 
275 276
      } // end if cnt==0

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

	// HV if needed 
	if (params.doHV_){
	  cxsum_hv = complex<r_4>(0.,0.);
300 301 302 303
	  if (params.doSigma_) {   // for computing Sigma, if required 
	    cxsum_hv_sq_rp = 0.;
	    cxsum_hv_sq_ip = 0.;
	  }
304
	}
OP's avatar
OP committed
305 306 307 308
      }
      
      //   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 
309 310 311 312 313
      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
314 315 316 317
      }

      for(size_t k=0; k<KVCXHH.size(); k++){
	cxsum.Row(k) += vismtx.Row(KVCXHH[k]);   // les cross-correlations 
318 319 320 321 322 323
	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
324
      }
325 326 327 328
      
      if (params.doVV_){
	for(size_t k=0; k<KVCXVV.size(); k++){
	  cxsum_vv.Row(k) += vismtx.Row(KVCXVV[k]);   // les cross-correlations 
329 330 331 332 333 334
	  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) ;
	  }
335 336 337 338 339 340
	}
      }

      if (params.doHV_){
	for(size_t k=0; k<KVCXHV.size(); k++){
	  cxsum_hv.Row(k) += vismtx.Row(KVCXHV[k]);   // les cross-correlations HV
341 342 343 344 345 346
	  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) ;
	  }
347 348 349
	}
      }

OP's avatar
OP committed
350 351


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

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

	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);
399 400 401 402 403 404 405
	    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);
406 407 408 409 410
	      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();
411 412
	      }
	    } //-- end of computing sigmas  
413 414 415 416 417 418 419
	  }  //----- 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);
420 421 422 423 424 425 426
	    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);
427 428 429 430 431 432
	      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();
	      } 
433
	    } //-- end of computing sigmas  
434 435 436
	  }  //----- end of loop over the 16 Xcor HV
	} // end HV option 

OP's avatar
OP committed
437 438
	double tdif =  cfdate.TimeDifferenceSeconds(cfdate,dateobs)/2.;
	timevec(TFMtmidx) = dateobs.TimeDifferenceSeconds(dateobs.ShiftSeconds (tdif ),datestart);	// centre du bin 
439
	ravec(TFMtmidx) = P4Coords::RAFromTimeTU(dateobs.ShiftSeconds (tdif ));
OP's avatar
OP committed
440 441 442 443 444 445
	TFMtmidx++;
	//  ... done 
	I=0;  cntnt++;
      }
      cnt++;
      if ((cnt>0)&&(cntnt%10==0)&&(cntnt>pcntnt)) {
446 447
	cout<<"visiavg/Info: TFM-Map fill cnt="<<cntnt<<" VisMtxCount="<<cnt
	    <<" /Max="<<wreader.getTotalNbWindows()<<" DateObs="<<dateobs<<endl;
OP's avatar
OP committed
448 449 450
	pcntnt=cntnt;
      }
    }
451
    cout<<"visiavg/Info: count="<<cnt*wreader.getWindowSize()<<" Visibility Matrices read "<<endl;
OP's avatar
OP committed
452

OP's avatar
OP committed
453
    // --- Sauvegarde cartes temps-frequence en fits 
perdereau's avatar
perdereau committed
454

OP's avatar
OP committed
455 456
    FitsInOutFile  * fos = NULL ;

OP's avatar
OP committed
457
    if (fitsoutfile.length()>=1){
458
      cout << " fitsoutfile :" <<fitsoutfile<<":"<< endl;
OP's avatar
OP committed
459 460 461
      fos = new FitsInOutFile(fitsoutfile, FitsInOutFile::Fits_Create);
    }
    
OP's avatar
OP committed
462

OP's avatar
OP committed
463
    POutPersist potfm(outfile);
OP's avatar
OP committed
464 465
    char bufnam[10];
    int numkey =0;
OP's avatar
OP committed
466 467 468
    // --- 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
469
    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
470 471 472 473
    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
474

OP's avatar
OP committed
475 476 477
      if (fos != NULL) {
	ext_names.push_back(tfm_names[k]);
	(*fos)<<  tfmap;
OP's avatar
OP committed
478
      }
479 480 481 482 483 484 485 486 487 488
      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
489
    }
OP's avatar
OP committed
490

OP's avatar
OP committed
491 492
    // --- 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
493

OP's avatar
OP committed
494
    const char* tfmCC_names[6]={"TFM_1H2H", "TFM_1H3H", "TFM_1H4H", "TFM_2H3H", "TFM_2H4H", "TFM_3H4H"};
OP's avatar
OP committed
495 496 497
    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
498 499
    for(int k=0; k<6; k++)  {   // loop over the 6 Xcor 
      TArray< complex<r_4> > & tfmap = vtfm[k];
OP's avatar
OP committed
500 501 502

      TArray< r_4 > & tfmap_sqpr = vtfm_rp_sq[k];
      TArray< r_4 > & tfmap_sqpi = vtfm_ip_sq[k];
OP's avatar
OP committed
503
      tfmap *= complex<r_4>((r_4)(1./((double)deltaIavg*(double)TFMfbin)), 0.);
OP's avatar
OP committed
504 505 506 507 508 509 510
      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);
      }
511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529
      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
530
    }
531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555
    // ========================== 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);
	}
556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574
	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  
575 576 577 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 604 605 606 607 608 609 610 611 612 613
      }
    } // 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);
	}
614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633

	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 
634 635 636
      }

    } // end option HV 
OP's avatar
OP committed
637 638 639 640 641 642 643 644 645 646
    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_;
    }
647 648

    cout  << " frequences limites "<< lim_freq(0) <<" ; "<< lim_freq(1)<<endl;
OP's avatar
OP committed
649
    potfm << PPFNameTag("FreqLims") << lim_freq ;
OP's avatar
OP committed
650
    potfm << PPFNameTag("TimeVec") << timevec ;
OP's avatar
OP committed
651
    potfm << PPFNameTag("RAVec") << ravec ;
OP's avatar
OP committed
652 653
    // --- FIN sauvegarde cartes temps-frequence 
    //    resu.Update();
OP's avatar
OP committed
654 655 656

    

OP's avatar
OP committed
657
    TVector <double> avg_freqs( myp4fre.getP4NbFreqChannels()/TFMfbin);
658
    
OP's avatar
OP committed
659
    double frbase =  myp4fre.freqstart_ + myp4fre.getP4FreqResolution()/2. ;
OP's avatar
OP committed
660
    for (int kf=0 ; kf< myp4fre.getP4NbFreqChannels()/TFMfbin ; kf++,frbase += myp4fre.getP4FreqResolution()*TFMfbin )
OP's avatar
OP committed
661
      avg_freqs(kf) = frbase ;
662
    
OP's avatar
OP committed
663
    potfm << PPFNameTag("FreqVec") << avg_freqs;
OP's avatar
OP committed
664

OP's avatar
OP committed
665
    if (fos != NULL) {
OP's avatar
OP committed
666 667
      ext_names.push_back("FreqsLims");
      (*fos)<< lim_freq ;
OP's avatar
OP committed
668 669 670 671 672 673 674 675 676 677
      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);
    }
678
    cout << " return code "<<rc<<endl;
OP's avatar
OP committed
679
    cout << resu;   // Update est fait lors du print
680

OP's avatar
OP committed
681 682 683 684 685 686 687 688 689 690 691 692 693 694
  }
  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; 
695
  }
OP's avatar
OP committed
696 697 698 699 700 701

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

}