// Author: Egon Pavlica <http://www.ses-ng.si/~pavlica/>

#include "Lattice.h"
#include "nrutil.h"
#include "nr.h"
#include <stdlib.h>
#include <time.h>//for random seed

//________________________________________
//  This class is for use to montecarlo
//  simulation of charge particle in a 
//  disordered molecular crystal
//  
//  it contains mechanism to get periodic
//  boundary conditions of lattice.
//
//  lattice contains integer coordinates
//  base cell size are floats
//  working space is by defaults 7x7x7 sites
//
//  Note:
//     - no bias here!!! see LSpace
//     - 1-based indexes of fspace(3D array)
//     - 0-based private index fdistances
//
//
// Logs:
//
//  $Log: Lattice.cc,v $
//  Revision 1.7  2003/06/26 14:46:16  cvs
//  found a bug in LSpace::GetEnergy() with fbias. Bug successfully resolved :)
//
//  Revision 1.6  2003/06/26 11:24:25  cvs
//  modified particle moving with LMiki and faster LDisorder operations
//
//  Revision 1.5  2003/06/26 09:27:11  cvs
//  added log fields...
//

ClassImp(Lattice)

 Lattice::Lattice(const int Na,const int Nb,const int Nc){
  //create float lattice [1..Na][1..Nb][1..Nc] in memory with nrc f3tensor
  //base cell size is 1x1x1
  SetBase(1.0,1.0,1.0);
  fNa=Na;fNb=Nb;fNc=Nc;
  fspace=f3tensor(1,fNa,1,fNb,1,fNc);
  if (fverbose) std::cout<<"Lattice created ["<<fNa*fa<<"A,"<<fNb*fb<<"A,"<<fNc*fc<<"A]n";
  fdistances=0;
}

 Lattice::~Lattice() {
  //destroy Lattice and particles and free memory space
  if (fdistances) {
    Tdelete(fdistances,0);
    fdistances=0;
  }
  if (fspace) {
    free_f3tensor(fspace,1,fNa,1,fNb,1,fNc);
    if (fverbose) std::cout<<"Lattice destroyed ["<<fNa<<","<<fNb<<","<<fNc<<"]n";
  }
}


 void Lattice::SetBase(const float a,const float b,const float c){
  //Set lattice base cell size
  fa=a;fb=b;fc=c;
};

 void Lattice::FillRandom(){
  //fill fspace with normal energy distribution of sigma=unity
  //use ScaleEnergy() to get diferent sigma
  static long seed=-time(NULL);
  if (!fspace) std::cerr<<"Lattice must be created first!!!n";
  else {
    if (fverbose) std::cout<<"Seed: "<<seed<<std::endl;
    for (int i=1;i<=fNa;i++) 
      for (int j=1;j<=fNb;j++) 
	for (int k=1;k<=fNc;k++) 
	  fspace[i][j][k]=gasdev(&seed);
    if (fverbose) std::cout<<"Lattice energies randomly setn";
  }
}

 void Lattice::ScaleEnergy(const float sigma){
  //scale energy values
  //sigma in eV!!!
  fsigma=sigma;
    for (int i=1;i<=fNa;i++) 
      for (int j=1;j<=fNb;j++) 
	for (int k=1;k<=fNc;k++) 
	  fspace[i][j][k]*=fsigma;
    
    if (fverbose) std::cout<<"Energies scaled with "<<fsigma<<" eVn";
}

 float Lattice::Distance(const int i,const int j,const  int k) const{
    //Calculates distance from the given lattice vector
    //this method has no control about limits since it
    // is intended to be used with small values.
  double sum;
  sum=i*i*fa*fa+j*j*fb*fb+k*k*fc*fc;
  return sqrt(sum);
}

 float Lattice::Distance(const int idx){
  //return distance from the given lattice index
  // is intended to be used with LDisorder::Step.
  if (!fdistances) PrecalculateDistances();
  return fdistances[idx];
}

 void Lattice::PrecalculateDistances(){
  //precalculates distances to speed up the calculation
  //index is calculated as in LDisorder::Step
  //WARNING!!!
  //set parameters here regarding LDisorder::SetWSpace
  //working space for example 7x7x7
  int a,b,c;
  a=b=c=7;//Working space
  int idx;
  int ac=a/2+1,bc=b/2+1,cc=c/2+1;
  
  if (fdistances) {
    Tdelete(fdistances,0);
    fdistances=0;
  }
  Tnew(fdistances,0,(a*b*c+1));
  
  for(int i=1;i<=a;i++)
    for(int j=1;j<=b;j++)
      for(int k=1;k<=c;k++){
	idx=b*c*(i-1)+c*(j-1)+k;
	fdistances[idx]=Distance(i-ac,j-bc,k-cc);
	//	if (fverbose) std::cout<<"i,j,k="<<i<<","<<j<<","<<k<<" idx="<<idx<<" dist:"<<fdistances[idx]<<"n";
      }
}

 float Lattice::GetEnergy(int i,int j,int k) const{
    //return the energy of the lattice at the current position.
    //periodic boundary condition is used
    //Details:
    //     This procedure returns the i,j,k is centered in the sublattice fWa,fWb,fWc 
    //     of latice fspace and periodic boudary condition is included 
    
    //periodic boundary condition
    if (i>fNa) i%=fNa;
    if (i<1) {i%=fNa;i+=fNa;} 
    if (j>fNb) j%=fNb;
    if (j<1) {j%=fNb;j+=fNb;} 
    if (k>fNc) k%=fNc;
    if (k<1) {k%=fNc;k+=fNc;}
    
    return fspace[i][j][k];
}

 float Lattice::GetEnergyDifference(int i,int j,int k,int di,int dj,int dk) const{
    //return the energy difference of the latice at the i,j,k position (e1) and
    //the lattice site in the i+di,j+dj,k+dk position (e2)
    // de=e2-e1
    //calls GetEnergy()
    //it is not faster than two GetEnergy()
  return GetEnergy(i+di,j+dj,k+dk)-GetEnergy(i,j,k);
}