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


using namespace std;

void transform(double* pos_pol, double* pos_cart);

/****************************************************************************/
  //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 != 2){
    cerr << "need input_file(1)\n";
    exit(1);
  }

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

  istringstream ist(file_buff);

  int nx = 4;
  int ny = 0;
  double theta = 0.;
  double phi = 0.;

  ist >> nx >> ny >> theta >> phi;



  ifstream inp; inp.open("POSCAR_INP");

  double Lx, Ly, Lz;
  double val;  

  int N = 24;

  char lines[1024];

  double X[N]; double Y[N]; double Z[N];

  inp >> lines >> lines >> lines >> lines >> lines >> lines >> lines;
  inp >> val;
  inp >> Lx >> val >> val;
  inp >> val >> Ly >> val;
  inp >> val >> val >> Lz;
  inp >> lines >> lines;
  inp >> val >> val;
  inp >> lines >> lines;
  inp >> lines;
  for(int i = 0; i < N; i++){
    inp >> X[i] >> Y[i] >> Z[i] >> lines >> lines >> lines;
  }




  double pos_pol[6]; double pos_cart[6];
  ofstream out; 
  char file[1024];


  pos_pol[0] = (Lx* nx)/9.;
  pos_pol[1] = (Ly* ny)/8.;
  pos_pol[4] = theta * M_PI/180.;
  pos_pol[5] = phi   * M_PI/180.;


  int count = 0;  





  for(int i1 = 0; i1 < 22; i1++){
    pos_pol[3] = 0.4 + 0.1 * i1; //H-H distance 
    for(int i2 = 0; i2 < 45; i2++){
      if(i2 <= 34){ 
        pos_pol[2] = -1.0 + i2*0.15;  // molecule-surface distance
      }
      else if(i2 > 34){
        pos_pol[2] = 4.1 + (i2-34)*0.25;
      }
    

    transform(pos_pol, pos_cart);

    sprintf(file,"POSCAR_%d.dat", count);
    out.open(file);


    out << "H2 + Pt(211) 4 layer slab \n";
    out << " 1.0\n";
    out << setprecision(15) << showpoint << Lx << " " << 0.0 << " " << 0.0 << "\n";
    out << setprecision(15) << showpoint << 0.0 << " " << Ly << " " << 0.0 << "\n";
    out << setprecision(15) << showpoint << 0.0 << " " << 0.0 << " " << Lz << "\n";
    out << "Pt H \n";
    out << "24 2 \n";
    out << "Selective dynamics\n";
    out << "Direct\n";
    for(int i = 0; i < N; i++){
      out << showpoint << setprecision(15) << setw(18) << X[i] << " " << setw(18) << Y[i]  << " " << setw(18) << Z[i] << " F F F \n";
    }
    out << showpoint << setprecision(15) << setw(18) << pos_cart[0]/Lx << " " << setw(18) << pos_cart[1]/Ly << " " << setw(18) << pos_cart[2]/Lz + Z[0] << " F F F\n";
    out << showpoint << setprecision(15) << setw(18) << pos_cart[3]/Lx << " " << setw(18) << pos_cart[4]/Ly << " " << setw(18) << pos_cart[5]/Lz + Z[0] << " F F F\n";
    count++;
    out.close();
   }
 }

}

void transform(double* pos_pol, double* pos_cart){

  double chi = 20.0*M_PI/180.;
  


  double tt = pos_pol[4];
  double pp = pos_pol[5];
  double cos_theta = cos(tt)*cos(chi) + cos(pp)*sin(tt) *sin(chi);
  if(fabs(cos_theta) > 1.0) cos_theta = (int)cos_theta;
  double theta = acos(cos_theta);  
  double sin_theta = sin(theta);
  double cos_phi = sin(tt)*cos(pp)*cos(chi) - cos(tt)*sin(chi);
  if(fabs(sin_theta) > 1e-10) cos_phi /= sin_theta;
  if(fabs(cos_phi) > 1.0) cos_phi = (int)cos_phi;
  double phi = acos(cos_phi);
 
  if(sin(pp) < 0.) phi = 2.*M_PI - phi;


  cout << "PHI " << phi*180./M_PI << "   THETA   " << theta*180./M_PI  << " R " << pos_pol[3] << " Z " <<  pos_pol[2] << "\n";


  pos_cart[0] = pos_pol[0] + 0.5*pos_pol[3]*sin(theta)*cos(phi);
  pos_cart[3] = pos_pol[0] - 0.5*pos_pol[3]*sin(theta)*cos(phi);
 
  pos_cart[1] = pos_pol[1] + 0.5*pos_pol[3]*sin(theta)*sin(phi);
  pos_cart[4] = pos_pol[1] - 0.5*pos_pol[3]*sin(theta)*sin(phi);

  pos_cart[2] = pos_pol[2] + 0.5*pos_pol[3]*cos(theta);
  pos_cart[5] = pos_pol[2] - 0.5*pos_pol[3]*cos(theta);

}