#include<iostream>
#include<fstream>
#include<string.h>
#include<sstream>
#include<math.h>
#include<iomanip>
#include<stdlib.h>
#include<stdio.h>
#include<time.h>
#include<unistd.h>
#include "QC_MOM_stream_Tint.h" //Tvib not equal to Trot
//#include "QC_MOM.h" //Tvib neq Trot

using namespace std;

extern  void spline_init(double*, double*, int, double*);
extern  void spline_interpolation(double*, double*, double*, int, double, double*, double* );
void transform(double* pos_pol, double* pos_cart, double* mom_pol, double* mom_cart, double ma, double mb);
void transform_back(double* pos_pol, double* pos_cart, double* mom_pol, double* mom_cart, double ma, double mb);

/****************************************************************************/
  //READ DATA FROM INPUT FILE
  int rem_com(char* filename, char* streamstring, int string_length){ 
    const char com_B = '#';
    const char com_E = '\n';
    int pos = 0;
    char cc;
    ifstream inf(filename);
    while(inf.get(cc)&& pos < string_length-1){
      if(cc != com_B) 
	streamstring[pos++] = cc;
      else{
	while(cc != com_E && inf.get(cc));
	streamstring[pos++] = com_E;
      }
    }
    streamstring[pos] = 0;
    if(pos == string_length-1){
      cerr << "Buffer size exceeded !\n"; exit(0);
    }
    return(strlen(streamstring));
  }
  //END READING
  /****************************************************************************/

int main(int argc, char* argv[]){
  if(argc != 3){
    cerr << "need input_file(1), output_pref\n";
    exit(1);
  }
  double pos_cart[6]; double pos_pol[6]; //cartesian and relative positions
  double mom_cart[6]; double mom_pol[6]; //cartesian and relative momenta



  /**********************************************
   **************reading data from input********
   ***********************************************/
  int    n_trac     = 1;  // number of trajectories
  double zstart     = 13.228083; // z_start in a.u. 
  int    qvib       = 0; // vibrational state
  int    ql         = 0; // rotational quantum number l
  int    qm         = 0; // rotational quantum number l
  double EZ0        = 1.4; // velocity in Z (a.u.)
  char  potfile[1024];   //asymtotic potential file
//  sprintf(potfile,"gasphase.dat");
  double TSTREAM    = 0.; //Temperature of the beam 
  //FOR BEAM DISTRIBUTION
  int TON           = 0;
  double alpha      = 0.;
  double Estream    = 0.;
  int seed_nr       = 0;
  int seed_nr2      = 0;

  int buff_length  =  65536;
  char*  file_buff =  new char[buff_length];
  rem_com(argv[1], file_buff,  buff_length);

  istringstream ist(file_buff);

  ist >> n_trac  >> zstart >> qvib   >> ql  >> qm  >> EZ0 
      >> potfile >> TSTREAM >> TON >> alpha >> Estream >> seed_nr >> seed_nr2;
  delete file_buff; 
  ist.str(""); 

  srandom((int)seed_nr);

  double massCl  = 64626.865; //ATOMICMASS OF Cl^35
  double massH   = 1822.8885; //ATOMICMASS OF H^1
  double mu      = massH*massCl/(massH+massCl); //reduced mass
  EZ0 /= 27.211383;
  double p_Z0    = -sqrt(2.*EZ0*(massCl+massH));
  double vstream = sqrt(2.*Estream/(massCl+massH))*419382.9; 

  double ZSHIFT = 0.3132617835809899*5.97317529910482*3.9213310361935179 * 1.8897261; //bohrradius

  ofstream log_f; log_f.open("LOG_FILE.log");
  log_f << " INITIAL SLAB CONFIGURATIONS FOR HCl + Au(111) \n";
  log_f << "             NORMAL INCIDENCE                  \n";
  log_f << " Z_START " << zstart/1.8897261 << " Angstrom \n";
  log_f << " ZSHIFT  " << ZSHIFT/1.8897261 << " Angstrom \n";
  log_f << " SEED NUMBER   " << (int)seed_nr << "\n";
  log_f << " SEED NUMBER 2 " << (int)seed_nr2 << "\n";
  log_f << " NUMBER OF POSCARS CREATED " << n_trac << "\n";
  if(TON == 1){
    log_f << " BEAM IS SWITCHED ON; Trot .neq. (Tvib == Tnozzle) \n "; 
    log_f << " Tnozzle : " << TSTREAM << " K \n";
    log_f << " Tvib    : " << TSTREAM << " K \n";
    log_f << " Trot    : " << (-181.102 + 0.648034*TSTREAM) << " K \n";
    log_f << " Estream : " << Estream << " eV\n";
    log_f << " vstream : " << vstream << " m/s \n"; 
    log_f << " width a : " << alpha   << " m/s \n";
  }
  else{ log_f << " NO BEAM CHOSEN \n"; exit(1);}
  log_f.flush();


  //READ GASPHASE POTENTIAL AND INITIALIZE POSITIONS AND MOMENTA
  QC_MOM BLA;

  double* XSTART = new double[n_trac]; double* XMOM_START = new double[n_trac];
  double* YSTART = new double[n_trac]; double* YMOM_START = new double[n_trac];
  double* ZSTART = new double[n_trac]; double* ZMOM_START = new double[n_trac];
  double* RSTART = new double[n_trac]; double* RMOM_START = new double[n_trac];
  double* TSTART = new double[n_trac]; double* TMOM_START = new double[n_trac];
  double* PSTART = new double[n_trac]; double* PMOM_START = new double[n_trac];

  int* QVIB  = new int[n_trac]; int* QJ = new int[n_trac]; 
  int* QJ_MJ = new int[n_trac];

  ofstream out_init; out_init.open("initial_pos_momenta.info");
  double x, px;
  BLA.INITPOS(potfile, mu, TON, alpha, vstream, TSTREAM);
  ofstream init_read; init_read.open("init_vals.dat");
   char file_out[1024];
   ofstream poscar;

  for(int i = 0; i < n_trac; i++){
//    cout << "BLA " << i << "\n";
    BLA.init_pos_mom();
    BLA.getposr(&x, &px, 0); 
    XSTART[i] = x; XMOM_START[i] = 0.;
    BLA.getposr(&x, &px, 1);
    YSTART[i] = x; YMOM_START[i] = 0.;
    ZSTART[i] = zstart; 
    if(TON == 0) ZMOM_START[i] = p_Z0; 
    else if (TON == 1) ZMOM_START[i] = (massCl+massH)*BLA.get_vstream();
    BLA.getposr(&x, &px, 3);
    RSTART[i] = x; RMOM_START[i] = px;
    BLA.getposr(&x, &px, 4);
    TSTART[i] = x; TMOM_START[i] = px;
    BLA.getposr(&x, &px, 5);
    PSTART[i] = x; PMOM_START[i] = px;
    out_init << showpoint << setprecision(12) << fixed 
             << XSTART[i] << " " << YSTART[i] << " " << ZSTART[i] << " " 
             << RSTART[i] << " " << TSTART[i] << " " << PSTART[i] << " "
             << XMOM_START[i] << " " << YMOM_START[i] << " " << ZMOM_START[i] << " " 
             << RMOM_START[i] << " " << TMOM_START[i] << " " << PMOM_START[i] << "\n";
    QVIB[i]  = BLA.get_qvib();
    QJ[i]    = BLA.get_ql();
    QJ_MJ[i] = BLA.get_qm();
   }
//   srand((int)seed_nr2);

  ifstream inp_slab;   char slab_file[1024];  string slab_line;

  //generate an input to read
  for(int i = 0; i < n_trac; i++){
    pos_pol[0] = XSTART[i]; pos_pol[1] = YSTART[i]; pos_pol[2] = ZSTART[i];
    pos_pol[3] = RSTART[i]; pos_pol[4] = TSTART[i]; pos_pol[5] = PSTART[i];
    //LOAD START MOMENTA
    mom_pol[0] = XMOM_START[i]; mom_pol[1] = YMOM_START[i];
    mom_pol[2] = ZMOM_START[i]; mom_pol[3] = RMOM_START[i];
    mom_pol[4] = TMOM_START[i]; mom_pol[5] = PMOM_START[i];
    transform(pos_pol, pos_cart, mom_pol, mom_cart, massH, massCl);

    double Ekin = 0.;
    for(int ij = 0; ij < 3; ij++) Ekin += mom_cart[ij] * mom_cart[ij] * 0.5 / massH;
    for(int ij = 3; ij < 6; ij++) Ekin += mom_cart[ij] * mom_cart[ij] * 0.5 / massCl;
    for(int ij = 0; ij < 6; ij++)  pos_cart[ij] = pos_cart[ij]/1.8897261; //position in angstrom
    for(int ij = 0; ij < 3; ij++)  mom_cart[ij] = (mom_cart[ij]/massH)*21.876913; //velocity in Angstrom/fs
    for(int ij = 3; ij < 6; ij++)  mom_cart[ij] = (mom_cart[ij]/massCl)*21.876913; //velocity in Angstrom/fs 

    //MAKE SURE THAT THE Z-VALUES ARE SHIFTED ACCORDING  TO THE UPPERMOST LAYER POSITION OF THE SLAB 
    pos_cart[2] += ZSHIFT/1.8897261;
    pos_cart[5] += ZSHIFT/1.8897261; 

    init_read << i << " " << Ekin*27.211383 << " " << ZMOM_START[i]*ZMOM_START[i]*0.5/(massH+massCl) * 27.211383 << " " 
              <<  QVIB[i] << " " << QJ[i] << " " << QJ_MJ[i] << "\n";

    int dir   = random()%20 ;
    int slab  = random()%1000;
    dir = dir + 1;
    

    init_read << dir << " " << slab << "\n";
    for(int ij = 0; ij < 6; ij++) init_read << showpoint << setprecision(12) << fixed << setw(16) << pos_cart[ij] << " ";
    init_read << "\n";
    for(int ij = 0; ij < 6; ij++) init_read << showpoint << setprecision(12) << fixed << setw(16) << mom_cart[ij]  << " ";
    init_read << "\n";





    sprintf(file_out,"POSCAR%d.dat",i); 
    poscar.open(file_out);
    sprintf(slab_file,"/home/gernot/AuCl/SRP_PROJECT/SLAB_170K/VERLET_SLABS/P%d.dat",dir);
    inp_slab.open(slab_file);

    for(int kk = 0; kk < slab; kk++){
      for(int hh = 0; hh < 42; hh++){
        getline(inp_slab, slab_line);
      }
    }    
    
    for(int kk = 0; kk < 5; kk++){
      getline(inp_slab, slab_line);
      poscar << slab_line << "\n";
    }
    for(int kk = 0; kk < 4; kk++) getline(inp_slab, slab_line);
    poscar << "Au H Cl \n";
    poscar << "16 1 1 \n";
    poscar << "Selective dynamics\n";
    poscar << "Cartesian\n";
    for(int kk = 0; kk < 16; kk++){
      getline(inp_slab, slab_line);
      poscar << slab_line << "\n";  
    }
    poscar <<  showpoint << setprecision(12) << fixed << setw(16) << pos_cart[0]  << " ";
    poscar <<  showpoint << setprecision(12) << fixed << setw(16) << pos_cart[1]  << " "; 
    poscar <<  showpoint << setprecision(12) << fixed << setw(16) << pos_cart[2]  << "  T  T  T \n"; 
    poscar <<  showpoint << setprecision(12) << fixed << setw(16) << pos_cart[3]  << " "; 
    poscar <<  showpoint << setprecision(12) << fixed << setw(16) << pos_cart[4]  << " ";
    poscar <<  showpoint << setprecision(12) << fixed << setw(16) << pos_cart[5]  << "  T  T  T \n";
     for(int kk = 0; kk < 17; kk++){
      getline(inp_slab, slab_line);
      poscar << slab_line << "\n";
    }
    poscar <<  showpoint << setprecision(12) << fixed << setw(16) << mom_cart[0]  << " "; 
    poscar <<  showpoint << setprecision(12) << fixed << setw(16) << mom_cart[1]  << " ";
    poscar <<  showpoint << setprecision(12) << fixed << setw(16) << mom_cart[2]  << "\n";
    poscar <<  showpoint << setprecision(12) << fixed << setw(16) << mom_cart[3]  << " ";
    poscar <<  showpoint << setprecision(12) << fixed << setw(16) << mom_cart[4]  << " ";
    poscar <<  showpoint << setprecision(12) << fixed << setw(16) << mom_cart[5]  << "\n";
    poscar.close();
    inp_slab.close();

 
 


//    write_poscar(file_out, pos_cart, mom_cart);

  }
  out_init.close();

  double sum = 0.;
  for(int i = 0; i < n_trac; i++) sum += ZMOM_START[i]*ZMOM_START[i];
  log_f << "AVERAGE KINETIC ENERGY IN Z : " << 0.5*sum/(massH+massCl) * 27.211383/n_trac << " eV\n";
  log_f << "------------------------------------------\n";
  log_f << "EXPECTATION RO-VIBR. ENERGY VALUE \n";
  log_f << "ANALYTICAL VALUE: " << BLA.GET_EXPECT_INT_EXPECTED()*27211.383 << " meV \n";
  log_f << "NUMERICAL  VALUE: " << BLA.GET_EXPECT_INT()*27211.383/n_trac << " meV \n";



}//END PROGRAM