// Utilisation de SOPHYA pour faciliter les tests ...
#include "sopnamsp.h"
#include "machdefs.h"

/* ---------------------------------------------------------- 
   Projet BAORadio/PAON4 - (C) LAL/IRFU  2008-2018

   Programme de lecture des fichiers matrices de visibilites de PAON4 
   produits par mfacq pour les voies Thermometres et RFI (9,10)
   R. Ansari  Univ. Paris Sud, & LAL CNRS/IN2P3 
   Juillet 2018 
   ---------------------------------------------------------- */

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

#include <iostream>
#include <string>

#include "pexceptions.h"
#include "array.h"
#include "histats.h"
#include "fitsioserver.h"

#include "resusage.h"

// include lecteur de fichiers visibilites 
#include "p4autils.h"
#include "visp4winreader.h"

int Usage(void)
{
  cout<<"--- rdthermrfi : Read PPF files for thermometer & RFI signal channnels and produce TF-Maps \n"<<endl;
  cout<<"  Usage: rdthermrfi [-options] \n";
  cout<<"  Specify -inth PathForThermometer or in datacards " <<endl;
  P4AnaParams::UsageOptions();
  cout<<endl;
  return 1;
}

//----------------------------------------------------
int main(int narg, const char* arg[])
{
  // --- Decoding parameters 
  if ((narg<2)||((narg>1)&&(strcmp(arg[1],"-h")==0)))  return Usage();
  P4AnaParams params;
  params.DecodeArgs(narg, arg);
  if (params.inpath_thermRFI_.length()<1) {
    cout<<"rdthermrfi/ERROR no path specified for thermometer/RFI signals ... (rdthermrfi -h for usage)" <<endl;
    return 2;
  }
  string outfile = params.outfile_;
  if (outfile.length()<1)  outfile = "thrfiavg.ppf";

  string fitsoutfile = params.fitsoutfile_;
  if (fitsoutfile.length()>=1) {
    fitsoutfile = "!"+fitsoutfile ; // adds '!' ?
  }

  int deltaIavg = params.TFMtimebin_;
  sa_size_t TFMfbin = params.TFMfreqbin_;
  
  int rc = 0;
  try {
    FitsIOServerInit();
    HiStatsInitiator _inia;

    ResourceUsage resu;

    VisiP4WindowReader wreader(params, true);  // true -> Read Therometer/RFI 

    long Imin = wreader.getReader().getSerialFirst();
    long Imax =  wreader.getReader().getSerialLast();
    long Istep =  wreader.getReader().getSerialStep();
    cout << "vis2ra/Info: processing visibility matrix serial/sequence number range "
	 <<Imin<<" <= seq <= " << Imax << " with step="<<Istep<<endl;
    cout << " WindowSize="<<wreader.getWindowSize()<<"  -> TotalNbWindows="<<wreader.getTotalNbWindows()<<endl;


    //vismtx: visib matrice, meanmtx: mean FFT coef, vismtx_center: visib matrix corrigee de l'offset (normalisee avec NPAQSUM)
    // un vecteur avec les temps 
    TVector< double > timevec(wreader.getTotalNbWindows()/deltaIavg); 
    TVector< double > ravec(wreader.getTotalNbWindows()/deltaIavg); 
    TMatrix< complex<r_4> > meanvismtx;
    TMatrix< complex<r_4> > visum; 

    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; 

    //----- 2 auto-corr (Thermometer/RFI) TimeFrequency maps
    TArray<r_4> TFM_therm;
    TArray<r_4> TFM_RFI;
    //----  1 cross-correlation TFM   (Thermometer <> RFI)
    TArray< complex<r_4> > TFM_cross;

    //---- for the time-freqency map filling    
    sa_size_t TFMtmidx=0;
    sa_size_t tfmSX, tfmSY;

    bool fgok=true;
    while (fgok) {
      //reads next visibility matrix window 
      fgok = wreader.Shift();
      if (!fgok)  break;
      
      TMatrix< complex<r_4> > vismtx = wreader.getAverageVisMtx(cfdate);
      if (cnt==0)  {    //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations
	visum = vismtx;  meanvismtx=vismtx;
	tfmSX=wreader.getTotalNbWindows()/deltaIavg;
	tfmSY=vismtx.NCols()/TFMfbin;
	TFM_therm.SetSize2D(tfmSX, tfmSY);
	TFM_RFI.SetSize2D(tfmSX, tfmSY);
	TFM_cross.SetSize2D(tfmSX, tfmSY);
	dateobs=cfdate;
	cnt++;
	continue;
      }
      if (I==0)  dateobs=cfdate; // start filling a new time bin
      // Accumulating 
      visum += vismtx;
      meanvismtx += vismtx;
      
      I++;    // incrementing DeltaTime counter 
      // 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 << "rdthermrfi/Warning: something wrong in the logic , (TFMtmidx="<<TFMtmidx<<") >= (tfmSX="<<tfmSX<<")"
	     << " for read count="<<cnt<<endl; 
      }
      else if (I==deltaIavg) {   // Filling TFM maps
	TVector<r_4> vac0 = real(visum.Row(0));   // Voie 1 : Thermometre 
	TVector<r_4> vac1 = real(visum.Row(2));   // Voie 2 : RFI
	TVector< complex<r_4> > vacX = visum.Row(1);  //  Cx Thermometere X RFI 
	for(sa_size_t jy=0; jy<tfmSY; jy++) {  // frequency binning 
	  TFM_therm(TFMtmidx, jy) = vac0( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
	  TFM_RFI(TFMtmidx, jy) = vac1( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
	  TFM_cross(TFMtmidx, jy) = vacX( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
	}
	double tdif =  cfdate.TimeDifferenceSeconds(cfdate,dateobs)/2.;
	timevec(TFMtmidx) = dateobs.TimeDifferenceSeconds(dateobs.ShiftSeconds (tdif ),datestart);	// centre du bin 
	ravec(TFMtmidx) = P4Coords::RAFromTimeTU(dateobs.ShiftSeconds (tdif ));
	TFMtmidx++;
	//  ... done 
	I=0;  cntnt++;
	visum = complex<r_4>(0.,0.);
      }
      cnt++;
      if ((cnt>0)&&(cntnt%10==0)&&(cntnt>pcntnt)) {
	cout<<"rdthermrfi/Info: TFM-Map fill cnt="<<cntnt<<" VisMtxCount="<<cnt
	    <<" /Max="<<wreader.getTotalNbWindows()<<" DateObs="<<dateobs<<endl;
	pcntnt=cntnt;
      }
    }  // End of redaing visibility matrices loop
    
    cout<<"rdthermrfi/Info: count="<<cnt*wreader.getWindowSize()<<" Visibility Matrices read "<<endl;
    if (cnt>0)  { 
      meanvismtx /= complex<r_4>((r_4)cnt, (r_4)0.);
      cout<<"  rdthermrfi/Info: Saving vismtx_mean to PPF file "<<outfile<<endl;
      POutPersist po(outfile);
      po<<PPFNameTag("meanvis_Th_RFI")<<meanvismtx;
      cout<<"  rdthermrfi/Info: TFM_therm , TFM_RFI, TFM_cross (TFM_ThXRFI) to PPF file "<<outfile<<endl;
      po<<PPFNameTag("TFM_therm")<<TFM_therm;
      po<<PPFNameTag("TFM_RFI")<<TFM_RFI;
      po<<PPFNameTag("TFM_ThXRFI")<<TFM_cross;
      po << PPFNameTag("TimeVec") << timevec ;
      po << PPFNameTag("RAVec") << ravec ;

      P4FreqBand myp4fre;
      TVector <double> avg_freqs( myp4fre.getP4NbFreqChannels()/TFMfbin);
      double frbase =  myp4fre.freqstart_ + myp4fre.getP4FreqResolution()/2. ;
      for (int kf=0 ; kf< myp4fre.getP4NbFreqChannels()/TFMfbin ; kf++,frbase += myp4fre.getP4FreqResolution()*TFMfbin )
	avg_freqs(kf) = frbase ;
    
      po << PPFNameTag("FreqVec") << avg_freqs;

      FitsInOutFile  * fos = NULL ;
      
      if (fitsoutfile.length()>=1){
	vector <string> ext_names;

	cout << " fitsoutfile :" <<fitsoutfile<<":"<< endl;
	fos = new FitsInOutFile(fitsoutfile, FitsInOutFile::Fits_Create);

	ext_names.push_back("TFM_therm");
	(*fos) << TFM_therm;
	ext_names.push_back("TFM_RFI");
	(*fos) << TFM_RFI;
	ext_names.push_back("TFM_ThXRFI_R");
	(*fos) << real(TFM_cross);
	ext_names.push_back("TFM_ThXRFI_I");
	(*fos) << imag(TFM_cross);

	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);

      }


    } 
    //    resu.Update();
    cout << resu;   // Update est fait lors du print
  }
  catch (PException& exc) {
    cerr << " rdthermrfi.cc catched PException " << exc.Msg() << endl;
    rc = 77;
  }  
  catch (std::exception& sex) {
    cerr << "\n rdthermrfi.cc std::exception :" 
         << (string)typeid(sex).name() << "\n msg= " 
         << sex.what() << endl;
    rc = 78;
  }
  catch (...) {
    cerr << " rdthermrfi.cc catched unknown (...) exception  " << endl; 
    rc = 79; 
  } 

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

}