//
// STL_particles.cpp
//    input STL, output particles, FCC lattice
//    (c) MIT CBA Neil Gershenfeld 6/4/20
// todo
//    test need for perturbation
//
// includes
//
#include <iostream>
#include <cmath>
#include <vector>
#include <algorithm>
#include <cstring>

using namespace std;
//
// variables
//
float xmin,xmax,ymin,ymax,zmin,zmax,d;
float *ax[3],*ay[3],*az[3],*sx[2],*sy[2];
//
// coordinate functions
//
float xpos(int x,int y,int z) {
   float dx;
   if (y%2 == 0) {
      if (z%3 == 0) dx = 0;
      else if (z%3 == 1) dx = 0.5;
      else dx = 0;
      }
   else {
      if (z%3 == 0) dx = 0;
      else if (z%3 == 1) dx = -0.5;
      else dx = 0;
      }
   return xmin+(x+(y%2)/2.0+dx)*d;
   }
float ypos(int y,int z) {
   float dy;
   if (z%3 == 0) dy = 0;
   else if (z%3 == 1) dy = 1/(2*sqrt(3));
   else dy = 1/sqrt(3);
   return ymin+(sqrt(3.0)*y/2.0+dy)*d;
   }
float zpos(int z) {
   return zmin+(sqrt(2.0/3.0)*z)*d;
   }
//
// main
//
int main(int argc, char** argv) {
   if (argc == 1) {
      cerr << "command line: STL_particles [arguments] < input.stl | particle sink " << endl;
      cerr << "   nlattice=(lattice size between x limits)" << endl;
      return 1;
      }
   int nlattice;
   if (nullptr == strstr(argv[1],"nlattice=")) {
      cerr << "STL_particles argument not recognized" << endl;
      return 1;
      }
   else
      nlattice = stoi(argv[1]+9);
   //
   // read number of triangles
   //
   int ntriangles;
   cin.seekg(80);
   cin.read((char*)&ntriangles,4);
   cerr << "STL_particles:" << endl;
   cerr << "   input " << ntriangles << " triangles" << endl;
   //
   // read triangles
   //
   ax[0] = new float[ntriangles];
   ax[1] = new float[ntriangles];
   ax[2] = new float[ntriangles];
   ay[0] = new float[ntriangles];
   ay[1] = new float[ntriangles];
   ay[2] = new float[ntriangles];
   az[0] = new float[ntriangles];
   az[1] = new float[ntriangles];
   az[2] = new float[ntriangles];
   sx[0] = new float[ntriangles];
   sx[1] = new float[ntriangles];
   sy[0] = new float[ntriangles];
   sy[1] = new float[ntriangles];
   xmin = 1e10; xmax = -1e10;
   ymin = 1e10; ymax = -1e10;
   zmin = 1e10; zmax = -1e10;
   float x0,y0,z0,x1,y1,z1,x2,y2,z2;
   float p = 1e-6; // perturbation to remove deneracies
   for (int t = 0; t < ntriangles; ++t) {
      //
      // skip normal
      //
      cin.seekg(3*4,ios::cur);
      //
      // read vertices
      //
      cin.read((char*)&x0,4);
      x0 += p;
      if (x0 > xmax) xmax = x0;
      if (x0 < xmin) xmin = x0;
      cin.read((char*)&y0,4);
      y0 += p;
      if (y0 > ymax) ymax = y0;
      if (y0 < ymin) ymin = y0;
      cin.read((char*)&z0,4);
      z0 += p;
      if (z0 > zmax) zmax = z0;
      if (z0 < zmin) zmin = z0;
      cin.read((char*)&x1,4);
      x1 += p;
      if (x1 > xmax) xmax = x1;
      if (x1 < xmin) xmin = x1;
      cin.read((char*)&y1,4);
      y1 += p;
      if (y1 > ymax) ymax = y1;
      if (y1 < ymin) ymin = y1;
      cin.read((char*)&z1,4);
      z1 += p;
      if (z1 > zmax) zmax = z1;
      if (z1 < zmin) zmin = z1;
      cin.read((char*)&x2,4);
      x2 += p;
      if (x2 > xmax) xmax = x2;
      if (x2 < xmin) xmin = x2;
      cin.read((char*)&y2,4);
      y2 += p;
      if (y2 > ymax) ymax = y2;
      if (y2 < ymin) ymin = y2;
      cin.read((char*)&z2,4);
      z2 += p;
      if (z2 > zmax) zmax = z2;
      if (z2 < zmin) zmin = z2;
      //
      // skip attribute
      //
      cin.seekg(2,ios::cur);
      //
      // sort by z and save
      //
      int i0,i1,i2;
      if ((z2 >= z1) && (z2 >= z0))
         if (z1 >= z0) {i0 = 0; i1 = 1; i2 = 2;}
         else {i0 = 1; i1 = 0; i2 = 2;}
      else if ((z1 >= z2) && (z1 >= z0))
         if (z2 >= z0) {i0 = 0; i1 = 2; i2 = 1;}
         else {i0 = 1; i1 = 2; i2 = 0;}
      else
         if (z2 >= z1) {i0 = 2; i1 = 0; i2 = 1;}
         else {i0 = 2; i1 = 1; i2 = 0;}
      ax[i0][t] = x0; ay[i0][t] = y0; az[i0][t] = z0;
      ax[i1][t] = x1; ay[i1][t] = y1; az[i1][t] = z1;
      ax[i2][t] = x2; ay[i2][t] = y2; az[i2][t] = z2;
      }
   cerr << "   xmin/max: " << xmin << ' ' << xmax << endl;
   cerr << "   ymin/max: " << ymin << ' ' << ymax << endl;
   cerr << "   zmin/max: " << zmin << ' ' << zmax << endl;
   d = (xmax-xmin)/(nlattice-1.0);
   int nx = nlattice;
   int ny = (ymax-ymin)/(d*sqrt(3.0)/2.0);
   int nz = (zmax-zmin)/(d*sqrt(2.0/3.0));
   //
   // loop over z slices
   //
   vector<float> xcross,px,py,pz;
   int nsegments;
   for (int iz = 0; iz < nz; ++iz) {
      float z = zpos(iz);
      //
      // loop over triangles to find slice crossings
      //
      nsegments = 0;
      for (int t = 0; t < ntriangles; ++t) {
         if ((az[2][t] > z) && (az[0][t] < z)) {
            //
            // crossing found, check side
            //
            if (az[1][t] > z) {
               float f10 = (z-az[0][t])/(az[1][t]-az[0][t]);
               x0 = ax[0][t]+f10*(ax[1][t]-ax[0][t]);
               y0 = ay[0][t]+f10*(ay[1][t]-ay[0][t]);
               }
            else {
               float f21 = (z-az[1][t])/(az[2][t]-az[1][t]);
               x0 = ax[1][t]+f21*(ax[2][t]-ax[1][t]);
               y0 = ay[1][t]+f21*(ay[2][t]-ay[1][t]);
               }
            float f20 = (z-az[0][t])/(az[2][t]-az[0][t]);
            x1 = ax[0][t]+f20*(ax[2][t]-ax[0][t]);
            y1 = ay[0][t]+f20*(ay[2][t]-ay[0][t]);
            //
            // sort y and save segment
            //
            if (y1 >= y0) {
               sx[0][nsegments] = x0; sy[0][nsegments] = y0;
               sx[1][nsegments] = x1; sy[1][nsegments] = y1;
               }
            else {
               sx[0][nsegments] = x1; sy[0][nsegments] = y1;
               sx[1][nsegments] = x0; sy[1][nsegments] = y0;
               }
            ++nsegments;
            }
         }
      //
      // loop over slice rows
      //
      for (int iy = 0; iy < ny; ++iy) {
         float y = ypos(iy,iz);
         //
         // loop over slice segments to find crossings
         //
         for (int s = 0; s < nsegments; ++s)
            if ((sy[1][s] >= y) && (sy[0][s] < y)) {
               float f = (y-sy[0][s])/(sy[1][s]-sy[0][s]);
               x0 = sx[0][s]+f*(sx[1][s]-sx[0][s]);
               xcross.push_back(x0);
               }
         //
         // sort crossings
         //
         sort(xcross.begin(),xcross.end());
         //
         // loop over slice columns
         //
         int pointer = 0;
         bool inside = false;
         for (int ix = 0; ix < nx; ++ix) {
            if (xcross.size() > 0) {
               float x = xpos(ix,iy,iz);
               if (x > xcross[pointer]) {
                  //
                  // passed pointer, toggle side
                  //
                  inside = !inside;
                  ++pointer;
                  if (pointer == xcross.size())
                     break;
                  }
               if (inside) {
                  //
                  // inside, save point
                  //
                  px.push_back(x);
                  py.push_back(y);
                  pz.push_back(z);
                  }
               }
            }
         xcross.clear();
         }
      }
   int number = px.size();
   cerr << "   output " << number << " particles" << endl;
   //
   // output particles
   //
   cout << "type:" << endl;
   cout << "FCC" << endl;
   cout << "format:" << endl;
   cout << "xyz float32" << endl;
   cout << "xmin:" << endl;
   cout << xmin << endl;
   cout << "xmax:" << endl;
   cout << xmax << endl;
   cout << "ymin:" << endl;
   cout << ymin << endl;
   cout << "ymax:" << endl;
   cout << ymax << endl;
   cout << "zmin:" << endl;
   cout << zmin << endl;
   cout << "zmax:" << endl;
   cout << zmax << endl;
   cout << "links per particle:" << endl;
   cout << 12 << endl;
   cout << "number:" << endl;
   cout << number << endl;
   cout << "pitch:" << endl;
   cout << d << endl;
   cout << "particles:" << endl;
   for (int i = 0; i < number; ++i) {
      cout.write((char*)&px[i],4);
      cout.write((char*)&py[i],4);
      cout.write((char*)&pz[i],4);
      }
   cout << "end:" << endl;
   }