/* Class for treating splitted-cell (cell bank) with respect to given translation vector

   Written by k-fleak 

*/

#pragma warning( disable : 4786 )
#include <vector>
#include <map>
#include <math.h>

#include "../vector3D.hpp"
#include "../atom.cpp"

class CellBank{
private: 
  std::map<std::vector<int>, std::vector<int> > cellMap;
  std::vector<int> baseAtomIndex;
  std::vector<std::vector<int> > splitCellArray;
  double eps;
  std::vector<double> translationVec;
  void splitCell(std::vector<Atom>& atomArray, std::vector<int>& splitNum);
  void calcTranslationVec(std::vector<int>& splitNum);
  void calcEquivAtom(std::vector<Atom>& atomArray, std::vector<int>& splitNum);

public:
  CellBank(std::vector<Atom>& atomArray, std::vector<int>& splitNum);
  ~CellBank();
  std::vector<std::vector<int> > getCellBank();
};


CellBank::CellBank(std::vector<Atom>& atomArray, std::vector<int>& splitNum){
  eps = 0.001;
  splitCell(atomArray, splitNum);
  calcTranslationVec(splitNum);
  calcEquivAtom(atomArray, splitNum);
  
}


CellBank::~CellBank(){}

std::vector<std::vector<int> > CellBank::getCellBank(){
  return splitCellArray;
}

void CellBank::calcTranslationVec(std::vector<int>& splitNum){
  for (int i = 0; i < splitNum.size(); ++i){
    translationVec.push_back( double(1.0) / splitNum[i]);
  }
}

void CellBank::calcEquivAtom(std::vector<Atom>& atomArray, std::vector<int>& splitNum){
  std::map<std::vector<int>, std::vector<int> >::iterator itmap;
  std::vector<int>::const_iterator itAtom;
  int atom1, atom2;
  double diffx, diffy, diffz;
  double intDiffx, intDiffy, intDiffz;
  double eps2 = 0.00001;
  for (int i = 0; i < baseAtomIndex.size() ; ++i){
    int baseAtom = baseAtomIndex[i]; // Atom index in base-cell
    std::vector<int> vtmp;
    itmap = cellMap.begin();
    vtmp.push_back(baseAtom);
    while ( itmap != cellMap.end()){
      itAtom = (*itmap).second.begin();
      while (itAtom != (*itmap).second.end()){
	atom1 = baseAtom;
	atom2 = (*itAtom);
	diffx = (atomArray[atom1].getPosition().x - atomArray[atom2].getPosition().x) / translationVec[0] + eps2;
	diffy = (atomArray[atom1].getPosition().y - atomArray[atom2].getPosition().y) / translationVec[1] + eps2;
	diffz = (atomArray[atom1].getPosition().z - atomArray[atom2].getPosition().z) / translationVec[2] + eps2;
	if ( fabs(floor(diffx) - diffx) > fabs(ceil(diffx) - diffx) ){
	  intDiffx = ceil(diffx);
	}else{
	  intDiffx = floor(diffx);
	}
	if ( fabs(floor(diffy) - diffy) > fabs(ceil(diffy) - diffy) ){
	  intDiffy = ceil(diffy);
	}else{
	  intDiffy = floor(diffy);
	}
	if ( fabs(floor(diffz) - diffz) > fabs(ceil(diffz) - diffz) ){
	  intDiffz = ceil(diffz);
	}else{
	  intDiffz = floor(diffz);
	}
	
	if ( (fabs(diffx - intDiffx) < eps) and (fabs(diffy - intDiffy) < eps) and (fabs(diffz - intDiffz) < eps) )  {
	  // std::cout << "o  " << atom1 << " " << atom2 << "   ->   " << diffx << "  " << diffy << "  " << diffz << std::endl;
	  vtmp.push_back(atom2);
	}else{
	  // std::cout << "x  " << atom1 << " " << atom2 << "   ->   " << diffx << "  " << diffy << "  " << diffz << std::endl;
	}
	++itAtom;
      }
      ++itmap;
    }
    splitCellArray.push_back(vtmp);
  }

  // Debug
  /*
  std::vector<std::vector<int> >::iterator it0;
  std::vector<int>::iterator it1;
  it0 = splitCellArray.begin();
  while( it0 != splitCellArray.end()){
    it1 = (*it0).begin();
    while ( it1 != (*it0).end()){
      std::cout << (*it1) << "  ";
      ++it1;
    }
    std::cout << std::endl;
    ++it0;
  }
  */
  //

}

void CellBank::splitCell(std::vector<Atom>& atomArray, std::vector<int>& splitNum){
  std::vector<Atom>::iterator itAtom;
  itAtom = atomArray.begin();
  vector3D<double> pos;
  int atomIndex = 0;
  while ( itAtom != atomArray.end()){
    pos = (*itAtom).getPosition();
    std::vector<int> mtmp;
    mtmp.push_back(static_cast<int>(floor(splitNum[0] * pos.x + eps)));
    mtmp.push_back(static_cast<int>(floor(splitNum[1] * pos.y + eps)));
    mtmp.push_back(static_cast<int>(floor(splitNum[2] * pos.z + eps)));
    cellMap[mtmp].push_back(atomIndex);
    ++itAtom;
    ++atomIndex;
  }
  std::vector<int> baseCell(3, int(0));
  baseAtomIndex = cellMap[baseCell];
  cellMap.erase(baseCell);

  // Debug
  /*
 for (int i = 0; i < baseAtomIndex.size(); ++i){
    std::cout << baseAtomIndex[i] << "  " << std::endl;
  }
  */
}