/*---------------------------- 
   PAON4 analysis : 
   Determining antenna pointing and dish size fitting parameters ... 
   (source track corresponds to a satellite track (or planet or sun ...)
   Does also fit of phases for the 6 cross-correlations
   Uses gacfit.h and gcxfit.h acbeam.h cxbeam.h  
   R. Ansari  , February 2019 

   Evolution of trkacfit.cc , intended to replace it 

   Assuming that the file trk_1.d contains :
@zenang 10 
@trk  dt_A21  1358,1410  1278.5  trk_A21_43057
@trk  dt_A21  1512,1555  1278.5  trk_A21_43055    

   Example command to run : 
csh> ../AnaPaon4/Objs/trkacxfit -D 4.5 -out A21_fac.txt  trk_1.d  


  visi-track-set datacard format 
----------------------------------
# input path (directory for files) listed here (optional)
@inpath  path 
## Zenith angle or declination shift 
# Zenith angle_in_degree (positive toward north)  initial value for fit 
@zenang ZenithAngle 
### Input data and track files - multiple @trk cards 
# filenames with .ppf extension , Tstart,Tend in minutes 
@trk AutoCrossVisiDataTableFile  Tstart,Tend  Freq  TrackFile [NOAC NOCX  NOACX]
### Optional last argument [ NOAC NOCX  NOACX ] not used currently 
 ----------------------------------------------------------------------------*/

#include <stdlib.h>
#include <string.h>

#include <iostream>
#include <fstream>

#include <cmath>
#include <complex>

#include "array.h"
#include "randr48.h"
#include "skymap.h"
#include "samba.h"
#include "fftwserver.h"
#include "strutilxx.h"

#include "ctimer.h"
#include "tarrinit.h"
#include "skymapinit.h"
#include "fitsioserver.h"
#include "slininterp.h"

#include "p4autils.h"

#include "acbeam.h"
#include "gacfit.h"
#include "gcxfit.h"
#include "gcxfitbaseline.h"
#include "trkfit.h"

using namespace std; 
using namespace SOPHYA; 

//---------------------------------------------------------------------------
//---  Grouping here beam & computation parameters -------------
//--- Global parameters 

// static double latitude=45.;   // latitude in degree 
// static double longitude=0.;   // longitude in degree 


static double D_dish = 4.5;       // effective dish size = effD*D 
static bool fggaussbeam=true;

// --- if true, perform Beam normalisation 
// static bool fgdobeamnormalise=false;

static string outfilename;            // output filename for auto-correlation fit parameters  
static string checkfilename;       // output filename for PAON4 and computed auto-cor using fitted parameters  

static bool fguseAac4Cx=true;       //if true, use fitted Amplitudes on autocor for initial Cross-cor fit parameter value

// static bool fgoutfits = false;   // true -> write output in fits format 

static int prtlevel=0;  // print level 

static string default_input_path;
static vector<string> trksetfiles;   // datacard files defining the track/data sets 

static bool do_cxfit = true ;       // if false, perform AutoCor fit only
// fg_cxB0: if true, fix offset B to zero (complex(0,0)) when fitting cross-correlations , false : fit B0 
static bool fg_cxB0 = true ;       // if false, fit the  offset B  when fitting cross-correlations 
static bool fg_phi0only = true;     // if false, linear varying phase with frequency, in AutoCor and Cross-Cor fits

// ---- fit simultane sur les 6 cross-correlations 
static bool do_baselinefit = false;   // if true , perform baseline fit 
static bool fg_fixbaseline = false;   // if true , do phases fit with fixed baselines 
static bool do_baselinesimplex = false;   // if true , perform baseline determination using Simplex minimisation  

//---   End of list of global parameters 



//-------- Declaration of functions in this file - code after the main, below 
int decode_args(int narg, char* arg[]);


//--------------------------------------------------------------
//---------  Main program -------------
//--------------------------------------------------------------
int main (int narg, char* arg[])
{
  cout << " ------ trkacxfit.cc : Fitting auto-correlations using track (satellites...) files  ------- " << endl;
  int rc = 0;
  try {
    SophyaInit();
    int rcda=decode_args(narg, arg);
    if (rcda)  return rcda;

    TrkFit_SetPrintLevel(prtlevel);
    TrkFit_FitLibInfo(); 
    vector<AcxDataSet> v_acxd;
    vector<TrackSet> v_trk;
    for(size_t i=0; i<trksetfiles.size(); i++) {
      TrkInputDataSet tids(trksetfiles[i], default_input_path);
      cout<<" --- Processing data-track set # "<<i+1<<" /max="<<trksetfiles.size()<<endl;
      cout<<tids;
      AcxDataSet  acxd(tids);
      TrackSet  tks(tids);
      ACxSetFitter acxfit(acxd, tks);
      char fnbuff[32];
      sprintf(fnbuff,"ac_%d_",(int)(i+1));
      string acofile = fnbuff+outfilename;
      acxfit.doACfit(acofile);  
      string acckfile = fnbuff+checkfilename;
      acxfit.saveExpectedAC(acckfile);
      if (do_cxfit)  {  //  Fit des cross-cor aussi 
	sprintf(fnbuff,"cx_%d_",(int)(i+1));
	string cxofile = fnbuff+outfilename;
	string cxckfile = fnbuff+checkfilename;
	acxfit.doCxfit(cxofile, fguseAac4Cx, fg_cxB0, fg_phi0only); 
	acxfit.saveExpectedCx(cxckfile);
      }
      v_acxd.push_back(acxd);
      v_trk.push_back(tks);
    }

    if (do_cxfit)  {  //  Resume resultat fit cross-cor
      cout<<" ----------  Summary phase fitting for cross-correlations ---------- " << endl;
      for(size_t i=0; i<v_acxd.size(); i++) {
	cout << "================================= TrackSet["<<i+1<<"] ================================="<<endl;
	v_acxd[i].PrintCxPhaseFitSummary(cout);
      }
      cout<<"=================================================================================="<<endl;
    }
    if (do_baselinefit || do_baselinesimplex)  {   // fit simultane des 6 cross-cor 
      string cxbofile = "cxb6_"+outfilename;
      string cxbckfile = "cxb6_"+checkfilename;
      CxBaselineFitter cxbfit(v_acxd, v_trk);
      //cxbfit.doSimplexMinimize();
      cxbfit.initFitParams();
      cxbfit.saveExpectedCx(cxbckfile);
      if (do_baselinesimplex) {
	cxbfit.doSimplexMinimize();
	cxbckfile="cxb6s_"+checkfilename;
	cxbfit.saveExpectedCx(cxbckfile);
	//      cxbfit.doCheck();
      }
      if (do_baselinefit) {
	cxbfit.dofit(cxbofile,fg_fixbaseline,fg_phi0only);
	cxbckfile = "cxb6f_"+checkfilename;
	cxbfit.saveExpectedCx(cxbckfile);
      }
    }
    {  // Resume resultat fits auto-correlations
      cout<<" ----------  Summary Dish diameter / pointing fit from Auto-Correlations ---------- " << endl;
      for(size_t i=0; i<v_acxd.size(); i++) {
	cout << "================================= TrackSet["<<i+1<<"] ================================="<<endl;
	v_acxd[i].PrintACFitSummary(cout);
      }
      cout<<"=================================================================================="<<endl;
    }
    if (do_cxfit)  {  //  Resume resultat fit cross-cor
      cout<<" ----------  Summary phase fitting for cross-correlations ---------- " << endl;
      for(size_t i=0; i<v_acxd.size(); i++) {
	cout << "================================= TrackSet["<<i+1<<"] ================================="<<endl;
	v_acxd[i].PrintCxPhaseFitSummary(cout);
      }
      cout<<"=================================================================================="<<endl;
    }

    Timer tm("trkacxfit");
  }  // End of try bloc 
  catch (PThrowable & exc) {
    cerr << " trkacxfit.cc: Catched Exception (PThrowable)" << (string)typeid(exc).name() 
         << " - Msg= " << exc.Msg() << endl;
    rc = 99;
  }
  catch (std::exception & e) {
    cerr << " trkacxfit.cc: Catched std::exception "  
         << " - what()= " << e.what() << endl;
    rc = 98;
  }
  catch (...) {
    ;    cerr << " trkacxfit.cc: some other exception (...) was caught ! " << endl;
    rc = 97;
  }
  cout << " ------------ END OF trkacxfit.cc  (Rc= " << rc << ") ------------ " << endl;
  return rc;
}


//-----------------------------------------------------------------------------------------
//-----------------------------------------------------------------------------------------


/* --Fonction-- */
int decode_args(int narg, char* arg[])
//! Decode program arguments 
{
  bool fglonghelp=false;
  if ((narg>1)&&(strcmp(arg[1],"-h")==0))  fglonghelp=true;
  if ((narg<2)||fglonghelp) {
    cout << "  trkacxfit : fit array geometry and   \n"
	 << "  Usage: trkacxfit [-options]  track_Set_1 [track_Set_2]  [track_Set_3 ...] \n"
	 << "  options: -inp def_input_path -out OutFileName -ckf CheckFileName \n" 
	 << "           -nocxf -cxfB -phifreq -docx6f -fixb -docx6s  \n"
	 <<"            -D dish_diameter -ngb  -prt PrintLevel  \n"<<endl;
    if (!fglonghelp) {
      cout << " trkacxfit -h  to get option description " << endl;
      return 2;
    }
    cout << "  -inp def_input_path : default directory path for input files (track_Set_J files) \n"
      	 << "    This path can be overwritten in each track_Set_J datacard file  \n"
      	 << "  -out OutFileName : Output file (text) to save fitted parameters vectors (def trkfit.txt) \n"
      	 << "  -ckf CheckFileName : Output file (PPF) to save expected signals (def chktrkfit.ppf)  \n"
	 << "    Output file names are constructed from OutFileName and CheckFileName prepending \n" 
	 << "  ac_JJ_ or cx_JJ_  where JJ=1...n is the track_set number \n"
      	 << "  -nocxf : Don't perform fit on cross-cor (AutoCor only fit) - default YES -> DO Cx fits) \n"
      	 << "  -cxfB : Fit the offset B when fitting cross-cor (default NO -> fix B0=complex(0,0)) \n"
	 << "  -phifreq : Use linearly varying phase with frequency model - default No fixed Phi=Phi0  \n"
	 <<"             Linear phi model : Phi(freq)=Phi_0 + a_Phi*(freq-1250)/250 \n"
	 << "  -docx6f : Perform simultaneous fit over 6 cross-corr (baseline and phases determination) - default NO  \n"
	 << "  -fixb : perform the previous fit (docx6f) with fixed baselines - default NO  \n"
	 << "  -docx6s : try simplex minimisation over 6 cross-corr (baseline and phases determination) - default NO  \n"
	 << "  -ngb : Use non gaussian beam profile (Bessel j1(angle) - default Gaussian beam)  \n"
	 << "  -D dish_diameter : define effective dish diameter (in meter D*eff , def=4.5) \n"
	 << "  -prt PrintLevel: specify print level \n" 
	 << "    track_Set_J (J=1..n) are the auto/cross correlation and track data sets  \n" 
	 << "    in datacard format [ @trk (multiple) @zenang cards @inpath] \n" 
	 << "    @trk AutoCrossVisDataTableFile Tstart,Tend Freq TrackFileName \n" 
	 << "    Filenames without .ppf extension, Tstart,Tend in minutes Freq in MHz for @trk cards \n" <<endl;
    return 1;
  }

  outfilename = "trkfit.txt";
  checkfilename = "chktrkfit.ppf";
  prtlevel=0;
  D_dish = 4.5;
  fggaussbeam=true;
  do_baselinefit=false;
  fg_fixbaseline=false;
  do_baselinesimplex=false;
  fg_cxB0=true;

  vector<string> lastargs;
  while (narg>1) {
    string fbo = arg[1];
    if (fbo=="-D") {  // dish size   (effective dish diameter) 
      if (narg<2) { cout << "trkacxfit/decode_args missing/bad argument, -h for help " << endl;  return -1; }
      D_dish=atof(arg[2]);   arg+=2; narg-=2;   lastargs.clear(); 
    }
    else if (fbo=="-ngb") {  // Use Non gaussian beam 
      fggaussbeam=false;   arg++;  narg--;  lastargs.clear(); 
    }
    else if (fbo=="-nocxf") {  // DON'T FIT Cross-Correlations 
      do_cxfit=false;   arg++;  narg--;  lastargs.clear(); 
    }
    else if (fbo=="-cxfB") {  // Determine B  when fitting Cross-cor
      fg_cxB0=false;   arg++;  narg--;  lastargs.clear(); 
    }
    else if (fbo=="-phifreq") {  // Use linear model for phase(frequency)  - If not Phi(freq) = Phi0  
      fg_phi0only=false;   arg++;  narg--;  lastargs.clear(); 
    }
    else if (fbo=="-docx6f") {  // Perform  simultaneous fit over 6 cross-corr 
      do_baselinefit=true;   arg++;  narg--;  lastargs.clear(); 
    }
    else if (fbo=="-fixb") {  // Perform  the previous docx6f fit with fixed baselines (phases only)
      fg_fixbaseline=true;   arg++;  narg--;  lastargs.clear(); 
    }
    else if (fbo=="-docx6s") {  // Perform  simultaneous fit over 6 cross-corr 
      do_baselinesimplex=true;   arg++;  narg--;  lastargs.clear(); 
    }
    else if (fbo=="-inp") {  // output file name 
      if (narg<2) { cout << "trkacxfit/decode_args missing/bad argument, -h for help " << endl;  return -1; }
      default_input_path=arg[2];   arg+=2; narg-=2;  lastargs.clear(); 
    }
    else if (fbo=="-out") {  // output file name 
      if (narg<2) { cout << "trkacxfit/decode_args missing/bad argument, -h for help " << endl;  return -1; }
      outfilename=arg[2];   arg+=2; narg-=2;  lastargs.clear(); 
    }
    else if (fbo=="-ckf") {  // output file name for checking (data & fitted 
      if (narg<2) { cout << "trkacxfit/decode_args missing/bad argument, -h for help " << endl;  return -1; }
      checkfilename=arg[2];   arg+=2; narg-=2;  lastargs.clear(); 
    }
    else if (fbo=="-prt") {  // print level 
      if (narg<2) { cout << "trkacxfit/decode_args missing/bad argument, -h for help " << endl;  return -1; }
      prtlevel=atoi(arg[2]);   arg+=2; narg-=2;  lastargs.clear(); 
    }
    else { arg++; narg--; lastargs.push_back(fbo); }
  }

  if (lastargs.size() < 1) {
    cout << "  trkacxfit/ERROR : No input track/data specified (trkacxfit -h for usage)  \n";
    return 9;
  }

  if (!do_cxfit) do_baselinefit=do_baselinesimplex=false;
  
  cout << " ------------------- trkacxfit/run parameters:"<<endl; 
  cout << " Beam: D_dish (initial guess)="<<D_dish<<" GaussianBeam ? "<<(fggaussbeam?"Yes":"No")<<endl;
  cout << " OutFileName= "<<outfilename<<" CheckFileName= "<<checkfilename<<endl;
  cout << " Perform fits on cross-correlations ? " << (do_cxfit?"Yes":"No")
       << "  Phase model : Phase(freq)= " << (fg_phi0only?"Phi0":"Phi0+aPhi*(freq-1250)/250") 
       << "  Fit B0 ? "<<(fg_cxB0?"No (B=(0,0)":"Yes")<<endl;
  cout << " Perform baseline fit ? " << (do_baselinefit?"Yes":"No")
       <<(fg_fixbaseline?"  With FIXED baselines":" (Phases & baselines)")
       <<" Simplex-Minimisation ? "<<(do_baselinesimplex?"Yes":"No")<<endl;
  cout << " --- TrackSetFiles: (NbFiles="<<lastargs.size()<<" default directory: "<<default_input_path<<")"<<endl;
  trksetfiles = lastargs;
  for (size_t i=0; i<trksetfiles.size(); i++) {
    cout << "["<<i<<"] : "<< trksetfiles[i]<<endl;
  }
  cout<<"---------------------------------------------------------------------------------"<<endl;
  return 0;
}