// CPP includes
#include <cstdlib>
#include <cstdio>
#include <iostream>
#include <cmath>
#include <vector>
#include <string>
#include <iostream>
#include <sstream>
#include <fstream>

// C includes
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

// ROOT includes
#include <TSystem.h>
#include <TROOT.h>
#include <TApplication.h>
#include <TRint.h>
#include <TDirectory.h>
#include <TFile.h>
#include <TChain.h>
#include <TTree.h>
#include <TNtuple.h>
#include <TBranch.h>
#include <TH1.h>
#include <TAxis.h>
#include <TCanvas.h>
#include <TColor.h>
#include <TPad.h>
#include <TStyle.h>
#include <TKey.h>
#include <TLegend.h>
#include <TText.h>
#include <TPaveText.h>
#include <TLine.h>
#include <TMarker.h>
#include <TLatex.h>
#include <TMath.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TF2.h>
#include <TString.h>
#include "TObjString.h"
#include "TRandom3.h"
#include "TGraph.h"
#include "TGraphSmooth.h"
#include "TSpline.h"

// BOOST
//#include <stdexcept>
#include <mpi.h>

// NANOFLANN & EIGEN
#include <nanoflann.hpp>
#include <Eigen/Dense>



#include </local/home/fbelloni/Software/Useful/PartialCSVParser/csv.h>
using namespace Eigen;
using namespace nanoflann;
using namespace std;




#define VERSION "v3.02-08-2016";

//#define IPM_Z_shift  -0.206//-0.372//-0.206 // -0.372//[m] //-0.206
void run(int, double , double , double , double ,double ); 

double IPM_Z_shift = -0.372;
// Other GLOBAL VARIABLES (these are used to read/store COMSOL file) 
  int nRow_file = 1030301;//6105475;//6105475;//12057692;//5256614;//1e+6+8;
  int nCol_file = 6;
  int nDim      = 3;
  Eigen::Matrix<float,Dynamic,Dynamic>  matCOMSOL(nRow_file, nCol_file); 
  Eigen::Matrix<float,Dynamic,Dynamic>  matREDUCED(nRow_file,nDim);
  typedef KDTreeEigenMatrixAdaptor< Eigen::Matrix<float,Dynamic,Dynamic> >  my_kd_tree_t;
  my_kd_tree_t   mat_index(nDim, matREDUCED);
  std::vector<float> query_pt(nDim);

void run(int tc, double ef, double ke, double sx, double sy,double sz){


  cout << " CIAO " << endl;  
   TH1::AddDirectory(kFALSE);
   TH2::AddDirectory(kFALSE);
   FILE* myfile;  
  
  
  //======================================================================================//
  //
  //                          LOAD THE COMSOL FILE 
  //
  //=======================================================================================//
   //nRow_file = counter_pt_inside;
   matCOMSOL.resize(nRow_file,nCol_file);
   double Com_X, Com_Y, Com_Z, Com_Ex, Com_Ey, Com_Ez;
   io::CSVReader<6> in2("/local/home/fbelloni/WORK/MATLAB2C++/COMSOL_EF/ESS_AUGUST_2018/IPM_30_372.csv");
   int cc=0;  
   int jj=0;
   while(in2.read_row(Com_X, Com_Y, Com_Z, Com_Ex, Com_Ey, Com_Ez)){
    matCOMSOL(cc,jj) =  roundf(Com_X * 10000) / 10000;//
    matCOMSOL(cc,jj+1) =roundf(Com_Y * 10000) / 10000;
    matCOMSOL(cc,jj+2) =roundf((Com_Z+IPM_Z_shift) * 10000) / 10000;
    matCOMSOL(cc,jj+3) = Com_Ex*20./30;
    matCOMSOL(cc,jj+4) = Com_Ey*20./30;
    matCOMSOL(cc,jj+5) = Com_Ez*20./30;
    cc++;
      //                                                                                              }
                                                                 }
}
 
 
 
 
 
 
 

//=====

int main(int argc, char **argv)
{
    
//   int ierr = MPI_Init(&argc, &argv);
//   int procid, numprocs;
//   ierr = MPI_Comm_rank(MPI_COMM_WORLD, &procid);
//   ierr = MPI_Comm_rank(MPI_COMM_WORLD, &numprocs);


  int    tc             = atoi(argv[1]);     // test charge
  double ef             = atof(argv[2]);   // electric field
  double ke             = atof(argv[3]);   // proton kinetic energy
  double sx             = atof(argv[4]);
  double sy             = atof(argv[5]);
  double sz             = atof(argv[6]);
  
      
  TRint *theApp = new TRint("App", &argc, argv, NULL, 0);
  run(tc, ef, ke, sx, sy,sz); 
  theApp->Run();
  
//   ierr = MPI_Finalize();

  return 0;
  

    }