
#include <vector>
#include <math.h>
#include <float.h>
#include <iostream>
#include "LatticeFunctions.hxx"
#include "CatomSim.hxx"


float firstSlice(string latticeType)
{
  if(latticeType == "cubic" || latticeType == "hexp" || latticeType == "fcc"||latticeType == "bcc"){
    return CatomSim::catom_radius;
  }

  return -1;
}


float nextSlice(string latticeType, float prevSlice)
{
  if(latticeType == "cubic" || latticeType == "hexp"){
    return prevSlice + CatomSim::catom_radius;
  }
  else if(latticeType =="fcc"){
    return prevSlice + (sqrt(2.0f)*CatomSim::catom_radius);
  }
  else if(latticeType =="bcc"){
    return prevSlice + (2*CatomSim::catom_radius/sqrt(3.0f));
  }
  return -1;
}


//returns the maximum possible number of neigbors
//a catom can have within a given lattice
int neighborCount(string latticeType){
  if(latticeType == "cubic") return 6;
  if(latticeType == "hexp") return 8;
  if(latticeType == "fcc") return 12;
  if(latticeType == "bcc") return 8;
  return -1;
}


//returns a vector of point3ds corresponding to the locations
//of origin's neigbors withing the given lattice
//NOTE: origin is assumed to be a point on the lattice
list<Point3D> latticeNeighbors(Point3D origin, string latticeType)
{
  float r = CatomSim::catom_radius;
  list<Point3D> result;
  origin = snapToLattice(origin, latticeType);

  if(latticeType == "cubic") {
    
    result.push_back(origin + Point3D(2*r, 0, 0));
    result.push_back(origin + Point3D(-2*r, 0, 0));
    result.push_back(origin + Point3D(0, 2*r, 0));
    result.push_back(origin + Point3D(0, -2*r, 0));
    result.push_back(origin + Point3D(0, 0, 2*r));
    result.push_back(origin + Point3D(0, 0, -2*r));
    
    
  }
  else if(latticeType == "hexp") {

    result.push_back(snapToLattice(origin + Point3D(0, 0, 2*r), "hexp"));
    result.push_back(snapToLattice(origin + Point3D(0, 0, -2*r), "hexp"));

    result.push_back(snapToLattice(origin + Point3D(2*r, 0, 0), "hexp"));
    result.push_back(snapToLattice(origin + Point3D(-2*r, 0, 0), "hexp"));

    float dx = cos(M_PI/3)*2*r;
    float dy = sin(M_PI/3)*2*r;

    result.push_back(snapToLattice(origin + Point3D(dx, dy, 0), "hexp"));
    result.push_back(snapToLattice(origin + Point3D(-dx, dy, 0), "hexp"));

    result.push_back(snapToLattice(origin + Point3D(dx, -dy, 0), "hexp"));
    result.push_back(snapToLattice(origin + Point3D(-dx, -dy, 0), "hexp"));

  }
  else if(latticeType == "fcc") {

    float d = (sqrt(2.0f)*r);
    
    result.push_back(snapToLattice(origin + Point3D(d, d, 0), "fcc"));
    result.push_back(snapToLattice(origin + Point3D(d, -d, 0), "fcc"));
    result.push_back(snapToLattice(origin + Point3D(-d, d, 0), "fcc"));
    result.push_back(snapToLattice(origin + Point3D(-d, -d, 0), "fcc"));

    result.push_back(snapToLattice(origin + Point3D(d, 0, d), "fcc"));
    result.push_back(snapToLattice(origin + Point3D(-d, 0, d), "fcc"));
    result.push_back(snapToLattice(origin + Point3D(0, d, d), "fcc"));
    result.push_back(snapToLattice(origin + Point3D(0, -d, d), "fcc"));

    result.push_back(snapToLattice(origin + Point3D(d, 0, -d), "fcc"));
    result.push_back(snapToLattice(origin + Point3D(-d, 0, -d), "fcc"));
    result.push_back(snapToLattice(origin + Point3D(0, d, -d), "fcc"));
    result.push_back(snapToLattice(origin + Point3D(0, -d, -d), "fcc"));

  }  
  else if(latticeType == "bcc") {

    float d = (2*r/sqrt(3.0f));
    
    
    result.push_back(snapToLattice(origin + Point3D(d, d, d), "bcc"));
    result.push_back(snapToLattice(origin + Point3D(d, -d, d), "bcc"));
    result.push_back(snapToLattice(origin + Point3D(-d, d, d), "bcc"));
    result.push_back(snapToLattice(origin + Point3D(-d, -d, d), "bcc"));
    
    result.push_back(snapToLattice(origin + Point3D(d, d, -d), "bcc"));
    result.push_back(snapToLattice(origin + Point3D(d, -d, -d), "bcc"));
    result.push_back(snapToLattice(origin + Point3D(-d, d, -d), "bcc"));
    result.push_back(snapToLattice(origin + Point3D(-d, -d, -d), "bcc"));

  }
  
  return result;
}

//This function takes an arbitrary point3d and returns
//a nearby point3D on the given lattice
Point3D snapToLattice(Point3D target, string latticeType)
{
  float r = CatomSim::catom_radius;
  float x = target.getX();
  float y = target.getY();
  float z = target.getZ();
  Point3D result(x,y,z);
  if(latticeType == "cubic") {
    result.setX(ceil(x/(2.0f*r))*(2.0f*r)-r);    
    result.setY(ceil(y/(2.0f*r))*(2.0f*r)-r);    
    result.setZ(ceil(z/(2.0f*r))*(2.0f*r)-r);  
  } 
  else if(latticeType == "hexp") {
    
    /*
      need to convert x(1,0) + y(0,1) (cartesian coordinates)
      into a(2*radius,0) + b(2*radius*cos(60), 2*radius*sin(60)) (hexagonal, catom-sized coordinates)
      from these equations it is apparent that
    */
    double b = target.getY()/(2*r*sin(M_PI/3));
    double a = (target.getX() - (b*2*r*cos(M_PI/3)))/(2*r);
    
    //these hexagonal coordinates must both be rounded so that they take up some integer number of catom-spaces:
    
    a = int(a+r/2*(a<0?-1:1));
    b = int(b+r/2*(b<0?-1:1));
    
    //now these even hexagonal coordinates can be converted back to cartesian:
    
    result.setX(a*2*r + b*2*r*cos(M_PI/3));
    result.setY(b*2*r*sin(M_PI/3));
    result.setZ(ceil(z/(2.0f*r))*(2.0f*r)-r);  
  }
  else if(latticeType == "fcc") {
    float d = (sqrt(2.0f)*r);
    
    result.setZ(ceil((z-r)/d-.0001)*d+r); 

    int zType = (int)ceil((z-r)/d-.0001) % 2;
    
    result.setX(ceil(x/(d)-.0001)*d);
    
    int xType = (int)ceil(x/(d)-.0001) % 2;
    
    result.setY(ceil(y/(2*d)-.0001)*2*d); 
    
    if((zType+ xType)%2){
      result.setY(result.getY()-d);
    }
  }
  else if(latticeType == "bcc") {
    float d = (2*r/sqrt(3.0f));
    
    result.setZ(ceil((z-r)/d-.0001)*d+r); 
    
    int zType = (int)ceil((z-r)/d-.0001) % 2;
    
    result.setX(ceil(x/(2*d)-.0001)*2*d);
    
    result.setY(ceil(y/(2*d)-.0001)*2*d); 
    
    if(zType){
      result.setY(result.getY()-d);
      result.setX(result.getX()-d);
    }
    
  }
  
  return result;
}
