#include<iostream>
#include<fstream>
#include<math.h>
#include<string>
#include<sstream>
#include<stdlib.h>
#include<iomanip>

using namespace std;

int main(void){

  ifstream xdat; xdat.open("XDATCAR");
  ofstream out; out.open("trajectory.dat", ios_base::app); 
  string line_inp;
  stringstream data_in;

  double L = 5.981689620000000;
  double LZ = 3.915749424;
  double ZMAX = 7.03;  //ANGSTROM
  double RMAX = 2.58; //ANGSTROM

  double ZSHIFT = 7.3374673808/(L*LZ); 

  double massH = 1.; double massCl = 35.453; 
  double massCOM = massH + massCl;

  double XH[3];  double XCl[3]; int count = 0;
  for(int i = 0; i < 7; i++) getline(xdat, line_inp);

  if(xdat.eof()) exit(1);


  while(true){
    getline(xdat, line_inp);
//    cout << line_inp <<"\n";
    for(int i = 0; i < 16; i++) getline(xdat, line_inp);
//    cout << line_inp <<"\n";
    xdat >> XH[0] >> XH[1] >> XH[2];
    xdat >> XCl[0] >> XCl[1] >> XCl[2]; 
//    cout << count << " " << XH[0] << " " << XH[1] << " " << XH[2] << " "
//                     << XCl[0] << " " << XCl[1] << " " << XCl[2] << "\n";
    count++;
    getline(xdat, line_inp);
    if(xdat.eof()) break;
 }
 count--;
 //cout << "START " << count << " " << XH[0] << " " << XH[1] << " " << XH[2] << " "
 //                    << XCl[0] << " " << XCl[1] << " " << XCl[2] << "\n";
 

 //CONVERT DIRECTS POSITIONS TO CARTESIANS
 XH[0] = (XH[0] + 0.5*XH[1])*L;
 XH[1] = sqrt(3)*0.5*XH[1]*L;
 XH[2] = (XH[2] - ZSHIFT)*L*LZ;

 XCl[0] = (XCl[0] + 0.5*XCl[1])*L;
 XCl[1] = sqrt(3)*0.5*XCl[1]*L;
 XCl[2] = (XCl[2] - ZSHIFT)*L*LZ;


 double Z = massH*XH[2]/massCOM + massCl*XCl[2]/massCOM;
 double RSAVE = 1000.;
 double XHSAVE[3]; 

 //cout << count << " " << XH[0] << " " << XH[1] << " " << XH[2] << " "
 //                     << XCl[0] << " " << XCl[1] << " " << XCl[2] << "\n";
   
 for(int i = -1; i <= 1; i++){
   double vshift = i*L;
   for(int j = -1; j <= 1; j++){
     double ushift = j*L;
     double x1 = ushift + 0.5*vshift;
     double y1 = sqrt(3)*0.5*vshift;

     double R = ((XH[0]+x1)-XCl[0])*((XH[0]+x1)-XCl[0])
              + ((XH[1]+y1)-XCl[1])*((XH[1]+y1)-XCl[1])
              + ((XH[2])-XCl[2])*((XH[2])-XCl[2]);
            R = sqrt(R);
     
     if(RSAVE > R){
       RSAVE = R; 
       XHSAVE[0] = XH[0] + x1;
       XHSAVE[1] = XH[1] + y1;
       XHSAVE[2] = XH[2];
     }           
   }

 }


 double rxy = (XHSAVE[0] - XCl[0])*(XHSAVE[0] - XCl[0]) 
           +  (XHSAVE[1] - XCl[1])*(XHSAVE[1] - XCl[1]); 
        rxy = sqrt(rxy);
 double tt  = rxy/RSAVE;
        tt = (XHSAVE[2] - XCl[2])/RSAVE;
        if(fabs(tt) > 1.0) tt = (int)tt;       
 tt = acos(tt)*180./M_PI;

 out << setw(8) << count << " " << setprecision(12) << fixed 
              << setw(15)  << RSAVE << " " 
              << setw(15)  << Z << " " 
              << setw(15)  << (massH*XHSAVE[0] + massCl*XCl[0])/massCOM << " " 
              << setw(15)  << (massH*XHSAVE[1] + massCl*XCl[1])/massCOM << " "  
              << setw(15)  << tt     << " "
              << setw(15)  << XH[0]  << " " 
              << setw(15)  << XH[1]  << " " 
              << setw(15)  << XH[2]  << " "
              << setw(15)  << XCl[0] << " " 
              << setw(15)  << XCl[1] << " " 
              << setw(15)  << XCl[2] << " ";
  

 if(RSAVE > RMAX){
  ofstream out2; out2.open("STOPCAR");
  out2 << "LSTOP = .TRUE.";
  out << " adsorbed ";
 }
 if(Z > ZMAX){
  ofstream out2; out2.open("STOPCAR");
  out2 << "LSTOP = .TRUE.";
  out << " scattered ";
 }



  
  


}