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

/* ---------------------------------------------------------- 
   Programme de lecture multi canaux pour BAORadio  
   R. Ansari, C. Magneville
   V : Mai 2009
   ---------------------------------------------------------- */

// 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 "histinit.h" 
#include "matharr.h" 
#include "timestamp.h"
/*
#include "fftpserver.h"
#include "fftwserver.h"

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


int Usage(bool fgshort=true);
// Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
int ProcSVFiles(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq, int card=1);

//----------------------------------------------------
//----------------------------------------------------
int main(int narg, char* arg[])
{
  if ((narg>1)&&(strcmp(arg[1],"-h")==0))  return Usage(false);
  if (narg<4) return Usage(true);

  HiStatsInitiator _inia;
  //   TArrayInitiator  _inia;

  int rc = 0;
  try {
    string inoutpath = arg[1];    
    int imin=0;
    int imax=0;
    int istep=1;
    sscanf(arg[2],"%d,%d,%d",&imin,&imax,&istep);
    int jf1=0;
    int jf2=0;
    int nfreq=0;
    sscanf(arg[3],"%d,%d,%d",&jf1,&jf2,&nfreq);
    cout << " ----- svv2mtx.cc Start - InOutPath= " << inoutpath << " IMin,Max,Step=" 
         << imin << "," << imax << "," << istep << " JF=" << jf1 << "," << jf2 << "  ------- " << endl;
    int card=1;
    if (narg>4) card=atoi(arg[3]); 
    ResourceUsage resu;
    rc=ProcSVFiles(inoutpath, imin, imax, istep, jf1, jf2, card);
    resu.Update();
    cout << resu;
  }
  catch (PException& exc) {
    cerr << " svv2mtx.cc catched MiniFITSException " << exc.Msg() << endl;
    rc = 77;
  }  
  catch (std::exception& sex) {
    cerr << "\n svv2mtx.cc std::exception :" 
         << (string)typeid(sex).name() << "\n msg= " 
         << sex.what() << endl;
    rc = 78;
  }
  catch (...) {
    cerr << " svv2mtx.cc catched unknown (...) exception  " << endl; 
    rc = 79; 
  } 

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

}



// Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
int ProcSVFiles(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq, int card)
{
  Timer tm("ProcSVFiles");
  char fname[512];
// NTuple 
  const char* nnames[10] = {"fcsm","ttsm","jfreq","s1","s2","s12","s12re","s12im","s12phi","s12mod"};
  NTuple nt(10, nnames);
  double xnt[15];
  uint_4 nmnt = 0;
  double ms1,ms2,ms12,ms12re,ms12im,ms12phi,ms12mod;

  TMatrix<r_4> s1, s2;
  TMatrix<r_4> v12re, v12im, v12phi,v12mod;
  sa_size_t ncols = (imax-imin+1)/istep;
  sa_size_t nrows = 10;
  sa_size_t kc=0;
  for(int ifile=imin; ifile<=imax; ifile+=istep) {
    if (card==2) 
      sprintf(fname, "%s/Ch34_%d.ppf",inoutpath.c_str(),ifile);
    else 
      sprintf(fname, "%s/Ch12_%d.ppf",inoutpath.c_str(),ifile);
    cout << " ProcSVFiles[" << ifile << "] opening file " << fname << endl;
    PInPersist pin(fname);
    string tag1="specV1";
    string tag2="specV2";
    string tag12="visiV12";
    if (card==2) {
      tag1 = "specV3";
      tag2 = "specV4";
      tag12="visiV34";
    }
    TVector<r_4> sv1;
    TVector<r_4> sv2; 
    TVector< complex<r_4> > vv12;
    pin >> PPFNameTag(tag1) >> sv1; 
    pin >> PPFNameTag(tag2) >> sv2; 
    pin >> PPFNameTag(tag12) >> vv12; 
    if (ifile==imin) {
      nrows = sv1.Size();
      cout << " ProcSVFiles/Info: Output s1,s2 matrix size NRows=NFreq=" 
           << nrows << " NCols=NFiles=" << ncols << endl;
      s1.SetSize(nrows, ncols);
      s2.SetSize(nrows, ncols);
      nrows = vv12.Size();
      cout << " ProcSVFiles/Info: Output v12 matrix size NRows=NFreq=" 
           << nrows << " NCols=NFiles=" << ncols << endl;
      v12re.SetSize(nrows, ncols);
      v12im.SetSize(nrows, ncols);
      v12phi.SetSize(nrows, ncols);
      v12mod.SetSize(nrows, ncols);
    }
    s1.Column(kc) = sv1;
    s2.Column(kc) = sv2;
    v12re.Column(kc) = real(vv12);
    v12im.Column(kc) = imag(vv12);
    v12phi.Column(kc) = phase(vv12);
    v12mod.Column(kc) = module(vv12);

// Calcul moyenne dans des bandes en frequence 
    int deltajf=(jf2-jf1)/nfreq;
    if (deltajf<1) deltajf=1; 
    for(int kf=0; kf<nfreq; kf++) {
      sa_size_t jfstart=jf1+kf*deltajf;
      sa_size_t jfend=jfstart+deltajf;
      if (jfend>jf2)  break;
      nmnt=0;  ms1=ms2=ms12=ms12re=ms12im=ms12phi=ms12mod=0.;
      for(sa_size_t jf=jfstart; jf<jfend; jf++) { 
        ms1 += s1(jf,kc);
        ms2 += s2(jf,kc);
        ms12re += v12re(jf,kc);
        ms12im += v12im(jf,kc);
        ms12phi += v12phi(jf,kc);
        ms12mod += v12mod(jf,kc);
      }
      nmnt = (jfend-jfstart);
      if (nmnt>0)  {
        double fnorm = (double)nmnt; 
        xnt[0] = ((int_8)(sv1.Info()["StartFC"])+(int_8)(sv1.Info()["EndFC"]))*0.5;
        xnt[1] = ((int_8)(sv1.Info()["StartTT"])+(int_8)(sv1.Info()["EndTT"]))*0.5;
        xnt[2] = kf;
        xnt[3] = ms1/fnorm;
        xnt[4] = ms2/fnorm;
        xnt[5] = ms12/fnorm;
        xnt[6] = ms12re/fnorm;
        xnt[7] = ms12im/fnorm;
        xnt[8] = ms12phi/fnorm;
        xnt[9] = ms12mod/fnorm;
        nt.Fill(xnt);
      }
    }
    kc++;

  }
  if (card==2) 
    sprintf(fname, "%s/Ch34mtx.ppf",inoutpath.c_str());
  else 
    sprintf(fname, "%s/Ch12mtx.ppf",inoutpath.c_str());

  cout << nt;
  cout << "ProcSVFiles: Opening file " << fname << " for writing" << endl;  
  POutPersist po(fname);
  string tag1="s1";
  string tag2="s2";
  string tag12r="v12re";
  string tag12i="v12im";
  string tag12p="v12phi";
  string tagnt="nt12";
  if (card==2) {
    tag1="s3";
    tag2="s4";
    tag12r="v34re";
    tag12i="v34im";
    tag12p="v34phi";
    tagnt="nt34";
    }
  po <<  PPFNameTag(tag1) << s1;
  po <<  PPFNameTag(tag2) << s2;
  po <<  PPFNameTag(tag12r) << v12re;
  po <<  PPFNameTag(tag12i) << v12im;
  po <<  PPFNameTag(tag12p) << v12phi;
  po <<  PPFNameTag(tagnt) << nt;
  cout << "ProcSVFiles: Matrices s1, s2, v12re, v12im, v12phi, NTuple nt written to file " << fname << endl;  
  return 0;
}




/* --Fonction-- */
int Usage(bool fgshort)
{
  cout << " --- svv2mtx.cc : Read PPF files produced by mcrd to make matrices BAORadio" << endl;
  cout << " Usage:  mcrd InOutPath Imin,Imax,step NumFreq1,NumFreq2,NBinFreq [card=1]" << endl; 
  if (fgshort) {
    cout << " mcrd -h for detailed instructions" << endl;
    return 1;
  }
  cout << " InOutPath : Input/Output directory name " << endl;
  cout << " Imin,Imax,IStep: Input PPF files sequence number \n"
       << "    FileNames=InOutPath/Ch12_II.fits Imin<=II<=Imax II+=IStep \n" 
       << " NumFreq1,NumFreq2,NBinFreq: Freq Zone and number of frequency bins for ntuple\n"
       << " card=1 Ch12 , card=2 Ch34 " << endl;
  return 1;
}