// PROBLEM APPEARS TO BE THAT THE SCALED COORDS DON'T RUN FROM -.5 to +.5
// Something is slightly wrong with the transformations.
#include<iostream>
#include<fstream>
#include<vector>
#include<ctime>
#include<cmath>
#include<cstdlib>
#include<algorithm>
#include"global.h"
#define PI 3.1415926536
using namespace std;

void boxptrsetup(int i, int j, int k, int l, int m, int n);

int main (int argc, const char *argv[]) {
   N = atoi(argv[1]);
   Nframes = atoi(argv[2]);
   R = atoi(argv[3]);   S = atoi(argv[4]);   T = atoi(argv[5]);
   const char *infname = argv[6];

   if (R < 3 || S < 3 || T < 3) {
           cout << "Need to use at least three boxes in each direction for code to work.  No speedup unless use at least 4 boxes in some direction." << endl;
           exit(EXIT_FAILURE);
           }

   x.resize(N); y.resize(N);     z.resize(N);             
   ix.resize(N);    iy.resize(N);    iz.resize(N);
   type.resize(N);

// matrices for scaling and unscaling coords -- needed to properly apply triclinic PBCs / min-img convention
   double h[3][3], hinv[3][3];
   double det; // determinant of h
   double scaledxq, scaledyq, scaledzq, scaledxr, scaledyr, scaledzr;
   double scaleddelx, scaleddely, scaleddelz;
   double scaledxlo, scaledylo, scaledzlo, scaledxhi, scaledyhi, scaledzhi;
   double xlo_bound, xhi_bound, ylo_bound, yhi_bound, zlo_bound, zhi_bound;
    
// various helper variables
   double delx, dely, delz, delsq, del;
   int f, i, j, k, l, m, n, o, p, q, r, t, u, v;
   int index, dum;
   int tox,toy,toz;    // used to place atoms in bins
   double Rd,Sd,Td;    // double versions of R,S,T
   Rd = R;      Sd = S;         Td = T;
// LAMMPS tilt factors: see https://lammps.sandia.gov/doc/dump.html and https://docs.lammps.org/Howto_triclinic.html
   double xy, xz, yz;

    
// for measuring density fluctuations
   vector<double> rho;      rho.resize(Nframes);    // overall number density
   vector<double> rhosq;    rhosq.resize(Nframes);
   double cellvol, rholocal, phi;
    
// for voronoi analyses -- use only if voronoi volumes are written to the config.dump file
   double vorov;  int nface;
//   vector<double> avgvorov;     avgvorov.resize(Nframes);
//   vector<double> avgvorovsq;   avgvorovsq.resize(Nframes);

// for contact analyses
   vector<double> Nc;  Nc.resize(Nframes);

// NEXT UP: add scaled avg. contact distance?
   vector<double> scavgcontdist;       scavgcontdist.resize(Nframes);

// input, output files
   char sponge[150];
   ifstream fin;   	fin.open(infname);
   if (!fin.is_open()) { cout << "Error: couldn't open configs file" << endl; exit(1); }
   ofstream fout;  

// timers
   time_t tstart, tend, tmid;
   tstart = time(0);

// Here do initial box setup
   vector <int> boxx, boxy, boxz;
   boxx.resize(N);		boxy.resize(N);		boxz.resize(N);
// set up bins
   vector< vector< vector<box> > > supercell(R,vector< vector<box> >(S,vector<box>(T)));
   vector< vector< vector<int> > > boxpops(R,vector< vector<int> >(S,vector<int>(T)));

// link bins; each box has a pointer to each of its neighbors, 3D-pbcs assumed
    for (i = 0; i < R; i++)  for (j = 0; j < S; j++) for (k = 0; k < T; k++) {
// each box has a pointer to each of its neighbors in 3D
       for (l = -1; l < 2; l++) for (m = -1; m < 2; m++) for (n = -1; n < 2; n++) {
          boxptrsetup(i,j,k,l,m,n);
          supercell[i][j][k].neighbor[l+1][m+1][n+1] = &supercell[pointvec[0]][pointvec[1]][pointvec[2]];
       }
    }

// main loop over frames
      for (f = 0; f < Nframes; f++) {
         tmid = time(0);
 //        cout << "Frame = " << f << ", time = " << difftime(tmid, tstart) << endl;
// read in frame header
         for (i = 0; i < 5; i++) fin.getline(sponge,150);
// note edits made for sheared-triclinic-cell format
// DO THESE NEED TO BE CHANGED TO XLO_BOUND, ETC?  YES!  See https://docs.lammps.org/Howto_triclinic.html
//         fin >> xlo >> xhi >> xy >> ylo >> yhi >> xz >> zlo >> zhi >> yz;
         fin >> xlo_bound >> xhi_bound >> xy >> ylo_bound >> yhi_bound >> xz >> zlo_bound >> zhi_bound >> yz;
         for (i = 0; i < 2; i++) fin.getline(sponge,150);
 
// Set cell edge lengths, cell volume, avg. number density
// First define corrected box bounds; see https://docs.lammps.org/Howto_triclinic.html
        vector<double> shiftxbounds{0,xy,xz,xy+xz};
        sort(shiftxbounds.begin(),shiftxbounds.end());
        xlo = xlo_bound - shiftxbounds[0];
        xhi = xhi_bound - shiftxbounds[3];
        ylo = ylo_bound - min(0.0,yz);
        yhi = yhi_bound - max(0.0,yz);
        zlo = zlo_bound;
        zhi = zhi_bound;
          
        DX = xhi - xlo;         DY = yhi - ylo;		DZ = zhi - zlo;
 //       cout << DX << " " << xhi_bound - xlo_bound << " " << DY << " " << yhi_bound - ylo_bound << " " << DZ << " " << zhi_bound - zlo_bound << endl;
         cellvol = (DX/R)*(DY/S)*(DZ/T);
         rho[f] = N/(DX*DY*DZ);
// estimated phi assumes 50:50 binary 1:1.4 mixture
          phi = rho[f]*.5*(PI/6)*(1 + 1.4*1.4*1.4);
 

// Next set cell side vectors -- used for PBCs: see https://docs.lammps.org/Howto_triclinic.html
// Are these the correct vectors?
// Vertical orientation of the vectors a, b, c, i.e. making them the COLUMN vectors of the shape matrix
         h[0][0] = DX;      h[0][1] = xy;       h[0][2] = xz;
         h[1][0] = 0;       h[1][1] = DY;       h[1][2] = yz;
         h[2][0] = 0;       h[2][1] = 0;        h[2][2] = DZ;
         cout << "Double-precision h matrix: " << endl;
         cout.precision(15);
         cout << DX << " " << xy << " " << xz << endl << 0 << " " << DY << " " << yz << endl << 0 << " " << 0 << " " << DZ << endl;
          
         det = h[0][0]*(h[1][1]*h[2][2] - h[2][1]*h[1][2]) -  h[0][1]*(h[1][0]*h[2][2] - h[2][0]*h[1][2]) +  h[0][2]*(h[1][0]*h[2][1] - h[2][0]*h[1][1]);
          phi = N*(PI/6)*(1 + 1.4*1.4*1.4)/det;
//         cout << "Frame = " << f << " phi = " << phi << " time = " << difftime(tmid, tstart) << endl;
         hinv[0][0] = h[1][1]*h[2][2] - h[2][1]*h[1][2];
         hinv[0][1] = h[0][2]*h[2][1] - h[0][1]*h[2][2];
         hinv[0][2] = h[0][1]*h[1][2] - h[0][2]*h[1][1];
         hinv[1][0] = h[1][2]*h[2][0] - h[1][0]*h[2][2];
         hinv[1][1] = h[0][0]*h[2][2] - h[0][2]*h[2][0];
         hinv[1][2] = h[1][0]*h[0][2] - h[0][0]*h[1][2];
         hinv[2][0] = h[1][0]*h[2][1] - h[2][0]*h[1][1];
         hinv[2][1] = h[2][0]*h[0][1] - h[0][0]*h[2][1];
         hinv[2][2] = h[0][0]*h[1][1] - h[1][0]*h[0][1];
         for (i = 0; i < 3; i++) for (j = 0; j < 3; j++) hinv[i][j] /= det;
// reads in positions for this frame -- tailored to match Kai's sheared-cohesive-grain dump files
         for (i = 0; i < N; i++) {
            fin >> index;   // >> dum;
            index -= 1;
            fin >> type[index] >> x[index] >> y[index] >> z[index] >> ix[index] >> iy[index] >> iz[index] >> vorov >> nface;
//            avgvorov[f] += vorov/N;
//            avgvorovsq[f] += vorov*vorov/N;
         }
         fin.getline(sponge,150);

// first step of putting atoms in boxes: set atoms' boxx,boxy,boxz, box populations
// code needed to be modified to handle triclinic boxes, use scaled coordinates
// Q: DO THESE NEED TO BE XLO_BOUND/XHI_BOUND, YLO_BOUND/YHI_BOUND, ZLO_BOUND/ZHI_BOUND instead?
// ARE THE SCALED COORDINATES ALL IN THE RANGE (-.5, .5)?
// NO!!! For example, one config. gives "Scaled box bounds: -0.500497 0.500772 -0.500138 0.500138 -0.5 0.5"
// Fixed one bug, now gives "Scaled box bounds: -0.500634 0.500634 -0.500138 0.500138 -0.5 0.5"
// Scaled box is symmetric about the origin but slightly too large along the x, y directions
         scaledxlo = hinv[0][0]*xlo + hinv[0][1]*ylo + hinv[0][2]*zlo;
         scaledylo = hinv[1][0]*xlo + hinv[1][1]*ylo + hinv[1][2]*zlo;
         scaledzlo = hinv[2][0]*xlo + hinv[2][1]*ylo + hinv[2][2]*zlo;
         scaledxhi = hinv[0][0]*xhi + hinv[0][1]*yhi + hinv[0][2]*zhi;
         scaledyhi = hinv[1][0]*xhi + hinv[1][1]*yhi + hinv[1][2]*zhi;
         scaledzhi = hinv[2][0]*xhi + hinv[2][1]*yhi + hinv[2][2]*zhi;
         cout << "Scaled box bounds: xlo " << scaledxlo << " xhi " << scaledxhi << " ylo " << scaledylo << " yhi " << scaledyhi << " zlo " << scaledzlo << " zhi " << scaledzhi << endl;
         for (i = 0; i < R; i++)  for (j = 0; j < S; j++) for (k = 0; k < T; k++) boxpops[i][j][k] = 0;
         for (q = 0; q < N; q++) {
            scaledxq = hinv[0][0]*x[q] + hinv[0][1]*y[q] + hinv[0][2]*z[q];
            scaledyq = hinv[1][0]*x[q] + hinv[1][1]*y[q] + hinv[1][2]*z[q];
            scaledzq = hinv[2][0]*x[q] + hinv[2][1]*y[q] + hinv[2][2]*z[q];
             if (scaledxq < -.5 || scaledxq > .5 || scaledyq < -.5 || scaledyq > .5 || scaledzq < -.5 || scaledzq > .5) {
//                cout << "q = " << q << " scaled coords " << scaledxq << " " << scaledyq << " " << scaledzq << endl;
//                exit(1);
             }
            tox = static_cast <int> ((scaledxq-scaledxlo)/(1.0/Rd));
            if (tox < 0) { tox = R-1; x[u] += DX; }
            else if (tox >= R) { tox = 0; x[u] -= DX; }
            boxx[q] = tox;
            toy = static_cast <int> ((scaledyq-scaledylo)/(1.0/Sd));
            if (toy < 0) { toy = S-1; y[u] += DY; }
            else if (toy >= S) { toy = 0;  y[u] -= DY; }
            boxy[q] = toy;
            toz = static_cast <int> ((scaledzq-scaledzlo)/(1.0/Td));
            if (toz < 0) { toz = T-1;  z[u] += DZ; }
            else if (toz >= T) { toz = 0;  z[u] -= DZ; }
            boxz[q] = toz;
            boxpops[tox][toy][toz] += 1;
         }

// set up atom lists of length boxpops
         for (i = 0; i < R; i++) for (j = 0; j < S; j++) for (k = 0; k < T; k++) {
            supercell[i][j][k].atomlist = new int[boxpops[i][j][k]];
            supercell[i][j][k].pos = 0;
         }
// second step of putting atoms in boxes
         for (u = 0; u < N; u++) {
            supercell[boxx[u]][boxy[u]][boxz[u]].atomlist[supercell[boxx[u]][boxy[u]][boxz[u]].pos] = u;
            supercell[boxx[u]][boxy[u]][boxz[u]].pos += 1;
         }

          
// Calculate density fluctuations
        for (i = 0; i < R; i++)  for (j = 0; j < S; j++) for (k = 0; k < T; k++) {
            rholocal = boxpops[i][j][k]/cellvol;
            rhosq[f] += rholocal*rholocal/(R*S*T);
        }
                    
// Calculate Z.  Assumes type=1 atoms have radius .5 and type 2 have radius .7
// PBC routine follows Tuckerman textbook
    for (i = 0; i < R; i++) for (j = 0; j < S; j++) for (k = 0; k < T; k++) for (o = 0; o < boxpops[i][j][k]; o++) {
        q = supercell[i][j][k].atomlist[o];   // index of the atom I'm currently mapping contacts for
// calculate scaled coordinates
        scaledxq = hinv[0][0]*x[q] + hinv[0][1]*y[q] + hinv[0][2]*z[q];
        scaledyq = hinv[1][0]*x[q] + hinv[1][1]*y[q] + hinv[1][2]*z[q];
        scaledzq = hinv[2][0]*x[q] + hinv[2][1]*y[q] + hinv[2][2]*z[q];
// loop over neighbors
        for (l = 0; l < 3; l++) for (m = 0; m < 3; m++) for (n = 0; n < 3; n++) for (p = 0; p < supercell[i][j][k].neighbor[l][m][n]->pos; p++) {
            r = supercell[i][j][k].neighbor[l][m][n]->atomlist[p];  // index of the atom testing whether in range of atom q
            if (r > q) {
// calculate scaled coordinates
                scaledxr = hinv[0][0]*x[r] + hinv[0][1]*y[r] + hinv[0][2]*z[r];
                scaledyr = hinv[1][0]*x[r] + hinv[1][1]*y[r] + hinv[1][2]*z[r];
                scaledzr = hinv[2][0]*x[r] + hinv[2][1]*y[r] + hinv[2][2]*z[r];
                scaleddelx = scaledxr - scaledxq;
                scaleddely = scaledyr - scaledyq;
                scaleddelz = scaledzr - scaledzq;
// min. image. conv in scaled coords.
                if (scaleddelx > 1.0/2) scaleddelx -= 1;
                else if (scaleddelx < -1.0/2) scaleddelx += 1;
                if (scaleddely > 1.0/2) scaleddely -= 1;
                else if (scaleddely < -1.0/2) scaleddely += 1;
                if (scaleddelz > 1.0/2) scaleddelz -= 1;
                else if (scaleddelz < -1.0/2) scaleddelz += 1;
// now convert back to unscaled coords to obtain the correct "del"
                delx = h[0][0]*scaleddelx + h[0][1]*scaleddely + h[0][2]*scaleddelz;
                dely = h[1][0]*scaleddelx + h[1][1]*scaleddely + h[1][2]*scaleddelz;
                delz = h[2][0]*scaleddelx + h[2][1]*scaleddely + h[2][2]*scaleddelz;
                delsq = delx*delx + dely*dely + delz*delz;
                del = sqrt(delsq);
// contact-counting
               if (del < 1.4) {
                  if (type[q] == 2 && type[r] == 2) Nc[f] += 2;
                  else if ((type[q] == 1 && type[r] == 2) || (type[q] == 2 && type[r] == 1)) if (del < 1.2) Nc[f] += 2;
                  else if (type[q] == 1 && type[r] == 1) if (del < 1) Nc[f] += 2;
               }
            }  // end of if (r > q) conditional
        } // end of the l-m-n-p loop
    } // end of the i-j-k-o loop


// delete atoms in linked cells
       for (i = 0; i < R; i++) for (j = 0; j < S; j++) for (k = 0; k < T; k++) delete [] supercell[i][j][k].atomlist;

   }  // end of f-loop
  fin.close();


// output Z
    char ofname[50];
    sprintf(ofname,"Zvsstrain");
    fout.open(ofname);
    for (f = 0; f < Nframes; f++) fout << f << " " << (1.0*Nc[f])/N << endl;
    fout.close();

    
// timer
   tend = time(0);
   cout << "Execution took " << difftime(tend, tstart) << " seconds." <<  endl;
   
   return 0;
}


// assignment of pointers to neighboring boxes (respecting PBCs)  
void boxptrsetup(int i, int j, int k, int l, int m, int n) {
   int o, p, q;
// PBCs expressed here: (l, m, n) is the vector from the pointer
    if ((i + l) == -1) o = R-1;
    else if ((i + l) == R) o = 0;				// handles x rollover
    else o = i + l;							
    pointvec[0] = o;                  
                   
    if ((j + m) == -1) p = S-1;
    else if ((j + m) == S) p = 0;				// handles y rollover
    else p = j + m;
    pointvec[1] = p;                                  											
    if ((k + n) == -1) q = T-1;					// handles z rollover
    else if ((k + n) == T) q = 0;                             								
    else q = k + n;                              									
    pointvec[2] = q;  

    return;                        
}