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

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

   visi2ntac: programme de lecture des fichiers matrices de 
   visibilites de PAON4, remplissage et ecriture d'un 
   DataTable avec la puissance des 8 voies en fonction du temps 
   R. Ansari, J.E. Campagne, C. Magneville   -  LAL/Irfu
   ---------------------------------------------------------- */

// 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 
#include "visip4reader.h"

//--------------------------- Fonctions de ce fichier   ------------------- 
int Usage(bool fgshort=true);
// int DecodeArgs(int narg, char* arg[]);

/* --Fonction-- */
int Usage(bool fgea)
{
  cout << " --- visi2ntac.cc : Read PPF files produced by mfacq time-frequency\n" << endl;
  if (fgea) cout << " visi2ntac: Missing or wrong argument ! " << endl;
  cout << " Usage: visi2ntac InPathBAO5 InPathBAO6 Imin,Imax GainPPFFile OutPPFFile Jf1A,Jf2A JF1B,Jf2B DeltaIAvg [PrtLev=0] \n"
       << "  InPathBAO5/6: path for input vismtx_J_iii.ppf on bao5/6 \n"
       << "  Imin,Imax : read files for Imin<=iii<=Imax  iii+=step \n"
       << "  GainPPFFile: Input gains PPF file name \n"
       << "  OutPPFFile: Output PPF file name \n"
       << "  Jf1A/B,Jf2A/B : frequency index (Jf1<=f<=Jf2) range for power integration, A:1300,1460 B:2130,2290 \n"
       << "  DeltaIAvg: compute average power every DeltaI input vismtx files, def=1 \n"<<endl;
  /*
  if (fgshort) {
    cout << " rdvisip4 -h for detailed instructions" << endl;
    return 1;
  }
  */
  return 1;
}

//----------------------------------------------------
//----------------------------------------------------
int main(int narg, char* arg[])
{
  if ((narg>1)&&(strcmp(arg[1],"-h")==0))  return Usage(false);
  if (narg<9) return Usage(true);
  string path5 = arg[1];
  string path6 = arg[2];
  int Imin,Imax,Istep=1; 
  sscanf(arg[3],"%d,%d",&Imin,&Imax);
  Istep=1;
  string gainfile=arg[4];
  string outfile=arg[5];
  // frequency range definition 
  sa_size_t JFmin=1300,JFmax=1460;
  int jf1,jf2;
  if ((narg>6)&&(strcmp(arg[6],"-")!=0))  sscanf(arg[6],"%d,%d",&jf1,&jf2);
  JFmin=jf1;  JFmax=jf2;
  sa_size_t JFminB=2130,JFmaxB=2290;  
  if ((narg>7)&&(strcmp(arg[7],"-")!=0))  sscanf(arg[7],"%d,%d",&jf1,&jf2);
  JFminB=jf1;  JFmaxB=jf2;
  // time averaging interval definition, in term of visibility files
  int deltaIavg=10;
  if (narg>8)  deltaIavg=atoi(arg[8]);
  if (deltaIavg<1)  deltaIavg=1;
  int prtlev=1;
  if (narg>9) prtlev=atoi(arg[9]);

  cout << " visi2ntac/Info: Path BAO5:"<<path5<<" BAO6:"<<path6<<"\n"
       << " Imin,max,step="<<Imin<<","<<Imax<<","<<Istep<<"  OutFile:"<<outfile<<"\n"
       << JFmin<<" <=NumFreq<= "<<JFmax<<"  "<<JFminB<<" <=NumFreqB<= "<<JFmaxB
       <<"  DeltaIAvg="<<deltaIavg<<" PrtLev="<<prtlev<<endl;

  cout << " visi2ntac/Info: reading input gain file"<<gainfile<<endl;
  Matrix gains;
  Vector gn;
  {
    PInPersist pin(gainfile);
    pin>>PPFNameTag("gains")>>gains;
    pin>>PPFNameTag("gn")>>gn;
    for (int p=0; p<8; p++) cout<<"gn["<<p<<"]="<<gn(p)<<" ";  
    cout<<endl;
  }
  // Bande de frequence de 2 MHz autour du 21 cm (1420 MHz -> JF=2792) 
  sa_size_t JFmin21=2776,JFmax21=2808;   // +/- 1 MHz autour de 1420 MHz 
  sa_size_t JFmin21_5=2710,JFmax21_5=2874;  //  +/- 5 MHz autour de 1420 MHz 
  cout << " pjHI +/-1 MHz @1420 MHz " << JFmin21<<" <=NumFreq<= "<<JFmax21
       << " pjHI +/-5 MHz @1420 MHz " << JFmin21_5<<" <=NumFreq<= "<<JFmax21_5 << endl;
 
  HiStatsInitiator _inia;
  //   TArrayInitiator  _inia;

  int rc = 0;
  try {
    ResourceUsage resu;
    vector<string> paths; 
    paths.push_back(path5); 
    paths.push_back(path6);

    const char* ntnames[34]={"k","time",
			     "p1","p2","p3","p4","p5","p6","p7","p8",
			     "p1HI","p2HI","p3HI","p4HI","p5HI","p6HI","p7HI","p8HI",
			     "p1HI5","p2HI5","p3HI5","p4HI5","p5HI5","p6HI5","p7HI5","p8HI5",
			     "p1B","p2B","p3B","p4B","p5B","p6B","p7B","p8B" };

    NTuple nt(34,ntnames);
    r_8 xnt[40];
    //  Numero de ligne des auto-correlations dans la matrice des visibilites 
    sa_size_t KVAC[8]={0,8,15,21,26,30,33,35};
    Range freqs(JFmin,JFmax);
    Range freqs21(JFmin21,JFmax21);
    Range freqs21_5(JFmin21_5,JFmax21_5);
    Range freqsB(JFminB,JFmaxB);

    VisiP4Reader vreader(paths, Imin,Imax,Istep,true);
    vreader.setPrintLevel(prtlev);
    bool fgok=true;
    TMatrix< complex<r_4> > vismtx;
    TimeStamp dateobs;
    double mttag;
    int cnt=0, cntnt=0, pcntnt=0;
    int I=0; 
    while (fgok) {
      fgok=vreader.ReadNext(vismtx, dateobs, mttag);
      if (fgok)  {
	if (I==0) {
	  xnt[0]=cnt;   xnt[1]=dateobs.SecondsPart();
	  if (prtlev>0) 
	    cout << "visi2ntac/Info:  dateobs="<<dateobs<<" SecondsPart()="<<dateobs.SecondsPart()<<endl;
	  for(int k=0; k<24; k++)  xnt[2+k]=0.;
	}
	for(int k=0; k<8; k++)  {    
	  TVector< complex<r_4> > vacz = vismtx.Row(KVAC[k]);
	  Vector vac(vismtx.NCols());
	  // correct for gain(nu)
	  for(sa_size_t c=0; c<vac.Size(); c++) vac(c)=vacz(c).real()/gains(k,c);
	  xnt[2+k] += vac(freqs).Sum();    // integration en bande large 
	  xnt[2+8+k] += vac(freqs21).Sum();   // integration a +/- 1 MHz @ 1420 MHz 
	  xnt[2+16+k] += vac(freqs21_5).Sum();   // integration a +/- 5 MHz @ 1420 MHz 
	  xnt[2+24+k] += vac(freqsB).Sum();    // integration en bande large B
	}
	I++; 
	if (I==deltaIavg) {
	  // puissance moyenne dans la bande large
	  double fnorm=(1./(double)deltaIavg)/(double)(JFmax-JFmin+1);
	  for(int k=0; k<8; k++)  xnt[2+k]*=fnorm;
	  // puissance moyenne ds +/5 MHz @1420 , exclu +/- 1 MHz centre a 1420 
	  fnorm=(1./(double)deltaIavg)/(double)((JFmax21_5-JFmin21_5)-(JFmax21-JFmin21));
	  for(int k=0; k<8; k++)  xnt[2+16+k] = (xnt[2+16+k]-xnt[2+8+k])*fnorm;
	  // puissance moyenne ds +/1 MHz @1420 
	  fnorm=(1./(double)deltaIavg)/(double)(JFmax21-JFmin21+1);
	  for(int k=0; k<8; k++)  xnt[2+8+k]*=fnorm;
	  // puissance moyenne en bande large B
	  fnorm=(1./(double)deltaIavg)/(double)(JFmaxB-JFminB+1);
	  for(int k=0; k<8; k++)  xnt[2+24+k]*=fnorm;
	  nt.Fill(xnt);
	  I=0;  cntnt++;
	}
	cnt++;
	if ((cntnt>0)&&(cntnt%10==0)&&(cntnt>pcntnt)) {
	  cout << " visi2ntac/Info: fill-count="<<cntnt<<" DateObs="<<dateobs<<endl;
	  pcntnt=cntnt;
	}
      } 
    }
    cout << " visi2ntac/Info: count="<<cnt<<" visimtx read "<<endl;
    cout << nt;
    cout<<"  rdvisip4/Info: Saving auto-corr=f(time) ntuple to PPF file "<<outfile<<endl;
    POutPersist po(outfile);
    po<<nt; 
    //    resu.Update();
    cout << resu;   // Update est fait lors du print
  }
  catch (PException& exc) {
    cerr << " visi2ntac.cc catched PException " << exc.Msg() << endl;
    rc = 77;
  }  
  catch (std::exception& sex) {
    cerr << "\n visi2ntac.cc std::exception :" 
         << (string)typeid(sex).name() << "\n msg= " 
         << sex.what() << endl;
    rc = 78;
  }
  catch (...) {
    cerr << " visi2ntac.cc catched unknown (...) exception  " << endl; 
    rc = 79; 
  } 

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

}