#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_stat.h" //Tvib not equal to 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);

void write_poscar(char* file, double* pos_cart, double* speed);

/****************************************************************************/
  //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

  time_t seed_time;
  time(&seed_time);
  srandom((int)seed_time);

  /**********************************************
   **************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 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;
  delete file_buff; 
  ist.str(""); 

  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; 



  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 << " 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);}


  //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];
  double x, px;
  BLA.INITPOS(potfile, mu, TON, alpha, vstream, TSTREAM);

  int stat_th[91];
  for(int i = 0; i < 91; i++) stat_th[i] = 0;

  for(int i = 0; i < n_trac; i++){
    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;

    int count = TSTART[i] *90./M_PI;
    stat_th[count]++;

  }

  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";


  ofstream out_t; out_t.open("theta_stat.dat");
  for(int i = 0; i < 91; i++) out_t << i << " " << (double)stat_th[i]/n_trac << "\n";




}//END PROGRAM