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

#include "LDisorder.h"

#include "nrutil.h"
#include "nr.h"
#include <time.h>
#include <stdlib.h>
#include <iostream>

//________________________________________
//  Simulation of hopping in dissordered 
// sistem
//
// an exaple of simulation is in file 
// final.cc. Full functionallity of this
// class is implemented in file hoppclient.cc
//
// with GetN() we get the time index (0..fmax-1)
// if the value is 0, the simulation didn't
// even started, if its value is equal to value
// returned by GetMax(), than the transient time 
// is longer than recording time, since the time 
// index could be used also to calculcate transient time
//
// The real time could be obtained also with 
// GetTTime(). If value is -1 than particle is trapped
// and has not traveled to oposite electrode.
// if value is 0, than particle exited space in first jump
//
// with GetKTrack() and with GetETrack() we
// get recorded tracks of electron...
// first track is particle position index (k axis base cell index)
// (if someone needs z coordinate the index should
// be multiplied with LDataHeader::fbasez) vs time
// second track is particle energy (eV) vs time
// every array field correspond to time index. Realtime
// is calculated using LDataHeader::fresolution
// array size is defined with LDataHeader::fmax 
//
// main function is Walk() and its clone WalkFast() 
// with almost no output bla bla.
// for 2D simulation use e.g. SetWSpace(7,1,7) or SetWSpace(1,7,7)
// for 1D simulation use e.g. SetWSpace(1,1,7)
// ...before starting simulation
//
//
// Test(mode) is for testing purposes
//
// Note:
//   - fprob zero-indexed
//   - fktrack zero-indexed
//   - fetrack zero-indexed
//   - dont forget to SetTimeResolution()
//
// Logs:
//
//  $Log: LDisorder.cc,v $
//  Revision 1.12  2004/07/13 09:57:57  cvs
//  found bug in tof normalization. solution is to multiply histogram with  tof->Scale(data->fvt->GetEntries()/(Double_t)data->fN/(header->rpmax-header->rpmin)/header->fresolution); and then the tof is normalized! This bug was corrected and commited to cvs.
//
//  Revision 1.11  2003/12/15 15:54:17  cvs
//  i think i found a bug in SaveProperties(int), when saving old value.
//
//  Revision 1.10  2003/11/11 17:33:03  cvs
//  LServer upgraded
//
//  Revision 1.9  2003/08/20 10:26:50  cvs
//  Added scripts directory with two utilities now:
//  modify_header which modify existing header and send it database
//  deliver_header program and script which deliver header to clients
//
//  libLComm was made ROOT compatible and GetBias(z) was added to LDataHeader
//
//  Revision 1.8  2003/07/25 11:36:26  cvs
//  simulation of IV measurements is done with iv program
//  some minority fixes to LDisorder - corrected ft to match the recorded position
//
//  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 13:26:05  cvs
//  some bugs removed + some test methods implemented in LDisorder
//
//  Revision 1.5  2003/06/26 11:24:25  cvs
//  modified particle moving with LMiki and faster LDisorder operations
//
//  Revision 1.4  2003/06/26 09:27:11  cvs
//  added log fields...
//

ClassImp(LDisorder)

 LDisorder::LDisorder():LMiki(){
  //empty constructor
  fe=0;
  fspace=0;
  fprob=0;
  fetrack=0;
  fktrack=0;
  ft=0.0;
  fN=0;
  fprob=0;
  SetTimeResolution(1.0);
  SetMax(10000);
};

 LDisorder::LDisorder(LSpace* space):LMiki(space){
  //default constructor. Default working space set to 7x7x7
  fe=0;
  fspace=0;
  fprob=0;
  fetrack=0;
  ftime=0.0;
  fktrack=0;
  ft=0.0;
  fN=0;
  fprob=0;
  SetTimeResolution(1.0);
  SetMax(10000);
};

 LDisorder::~LDisorder(){
  //default destructor
  if (fprob) Tdelete(fprob,0);
  if (fktrack) Tdelete(fktrack,0);
  if (fetrack) Tdelete(fetrack,0);
}

 void LDisorder::SetWSpace(const int a,const int b,const int c){
  //set working space size-lattice sites that are included in probablility calculation
  //prefered odd numbers!!
  fWa=a,fWb=b;fWc=c;
  fNprob=fWa*fWb*fWc;
  if (fprob) Tdelete(fprob,0);
  Tnew(fprob,0,fNprob);//for storing probability
  fprob[0]=0.0; //this must be set, used in Walk time calculation
  fWa2=fWa/2,fWb2=fWb/2,fWc2=fWc/2;
  fWa2c=fWa/2+1,fWb2c=fWb/2+1,fWc2c=fWc/2+1;
  idxexclude=fWb*fWc*fWa2+fWc*fWb2+fWc2+1;
  if (fWa2==fWa2c) std::cerr<<"Error: SetWSpace(a) - odd number needed!n"; 
  if (fWb2==fWb2c) std::cerr<<"Error: SetWSpace(b) - odd number needed!n"; 
  if (fWc2==fWc2c) std::cerr<<"Error: SetWSpace(c) - odd number needed!n"; 
}

 void LDisorder::SetMax(const long max){
  //Set maximum number of records
  //must be called before Walk() and similar methods
  //which need fktrack and fetrack
  if (!fktrack) fmax=max;
  else std::cerr<<"Error: SetMax() called too late, ignored!n"; 
}

 void LDisorder::SetTimeResolution(const float delta){
  //Set time resolution of records
  //must be called before Walk() and similar methods,
  //where the particle track is recorded
  if (!fN) ftime=delta;
  else std::cerr<<"Error: SetTimeResolution() called too late, ignored!n"; 
}

 void LDisorder::Init(){
  //get fe,fspace,fwspace and fmonitor from LMiki class
  
  if (!fe) fe=GetParticle();
  if (!fspace) fspace=GetSpace();
  if (!fktrack) Tnew(fktrack,0,fmax-1);
  Tclear(fktrack,0,fmax-1);
  if (!fetrack) Tnew(fetrack,0,fmax-1);
  Tclear(fetrack,0,fmax-1);
  ft=0; //time is measured from zero
  fN=0;
  if (!fprob) SetWSpace(7,7,7);
  //this was in Probability
  // boltzkT=ftemperature*1.3806503e-23;// J 
  boltzkTeV=8.617342e-5*ftemperature;// eV 
}

 void LDisorder::Walk(){
  //do a disorder hopping simulation of the particle
  //for more information get the report of this project.
  //for further references get the article:
  //     H.Baessler:phys.stat.sol.B 175, 15, 1993.
  //
  //Note:
  //   If monitor is set, then it outputs to monitor. 
  //   this is nice but VERY slow
  
  Init();//linking to other classes
  //physical properties here
  //I set them elsewhere
  if (fverbose) std::cout<<"Disorder Hopping simulationn";
  //initial savings
  SaveProperties(1);//first position

  if (!fspace->IsInside(fe)) return;

  while (1){
    if (fN>=fmax) break;//out of time

    if (ft<=fe->time) {
      SaveProperties(1);
      continue;
    }

    if (!Step()) break;//out of space
  }
  ft-=ftime;//since the ft points to the time where should store next
    
  //some verbosity
  if (fN<fmax) 
    if (fverbose) 
      std::cout<<"Particle exited the space at "<<fe->i<<","<<fe->j<<","<<fe->k<<std::endl; 
    else {
      if (fverbose) 
	std::cout<<"Particle timeout"<<std::endl;
    }
}


 int LDisorder::Where(const float a) const {
  //return position k in array, which satisfies arr[k-1]<a<=arr[k]
  //used with fprob
  unsigned long ju=fNprob+1,jm,jl=0;
  
  //if(a>fprob[fNprob-1]) return fNprob;//obsolete
  if(a<=fprob[1]) return 1;

  while(ju-jl>1){
    jm=(ju+jl)>>1;//compute midpoint 
    if(a>fprob[jm]) jl=jm; else ju=jm;
  }
  return ju;
}

 long LDisorder::initseed() const{
  //inits seed
  long seed=-time(NULL);
  //if (fverbose) std::cout<<"Seed: "<<seed<<std::endl;
  return seed;
}

inline
 int LDisorder::Step(){
  //A single step in dissorder hopping simulation
  //time is set as exponentially distributed deviate with average value of inverse of probability of jump
  static long seed=initseed();//for random number generator

  int idx;//unity index of lattice positions

  int i0=fe->i-fWa2c;
  int j0=fe->j-fWb2c;
  int k0=fe->k-fWc2c;

  float e1=fe->GetEnergy();//current energy 
  float e2,r,sum=0.0;

  //calculate probabilities
  for(int i=1;i<=fWa;i++)
    for(int j=1;j<=fWb;j++)
      for(int k=1;k<=fWc;k++){
	idx=fWb*fWc*(i-1)+fWc*(j-1)+k;
	if (idx!=idxexclude){//exclude sumation over the site of origin

	  e2=fspace->GetEnergy(i0+i,j0+j,k0+k);
	  r=fspace->Distance(idx);//to speed up
	  //r=fspace->Distance(i-fWa2c,j-fWb2c,k-fWc2c);

	  sum+=Probability(r,e1,e2);
	  fprob[idx]=sum;
	} else fprob[idx]=sum;  //so will never be selected in Where
      }
  idx=Where(sum*ran1(&seed));//where to jump?
  //inverse tranformation
  int x,z,y;
  idx--;
  x=idx/(fWb*fWc);
  y=(idx-x*fWb*fWc)/fWc;
  z=idx-x*fWb*fWc-y*fWc+1;
  idx++;x++;y++;

  //  fe->time+=(expdev(&seed)*0.166667*sum*exp(2.0*fgamma));
  //fe->time+=(expdev(&seed2)*0.166667*exp(2.0*fgamma));
  fe->time+=(expdev(&seed)*sum/(fprob[idx]-fprob[idx-1]));//no scale

  int answ=MoveTo(x-fWa2c,y-fWb2c,z-fWc2c);//move particle

  if (fverbose) if (e1>fe->GetEnergy()) std::cout<<"."; else std::cout<<"*";//down the energy and up the energy
  if (fverbose) std::cout<<"Time: "<<fe->time<<" Energy: "<<fe->GetEnergy()<<" eV Position: "<<fe->i<<","<<fe->j<<","<<fe->k<<std::endl;

  return answ;
}

inline
 float LDisorder::Probability(const float r,const float e1,const float e2) const{
  //Miller-Abrahams probabiliry for jump of length r 
  //from energy e1 to energy e2 (in eV). fgamma, and ftemperature 
  //are used too. Here it is assumed, that electron 
  //orbitals are spherical. If not, the r must be vector!
  //this is for HOLES!!!!!
  float bla1=exp(-2.0*fgamma*r);
  //this is in Init now
  // boltzkT=ftemperature*1.3806503e-23;// J 
  // boltzkTeV=8.617342e-5*ftemperature;// eV 
  //FOR HOLES
  if (e1>e2) return bla1*exp((e2-e1)/boltzkTeV);
  return bla1;
  //FOR ELECTRONS
  //  if (e2>e1) return bla1*exp((e1-e2)/boltzkTeV);
  //  return bla1;
}

inline
 void LDisorder::SaveProperties(const int save){
  //save particle properties (position,energy) and increment fN and ft
  //if save=1 it save position else it copy previous
  if (save) 
    SaveProperties(fe->k,fe->GetEnergy());
  else
    SaveProperties(fktrack[fN-1],fetrack[fN-1]);
}

inline
 void LDisorder::SaveProperties(const float position,const float energy){
  //save particle properties (position,energy) and increment fN and ft
  fktrack[fN]=position;
  fetrack[fN]=energy;
  fN++;
  ft+=ftime;

}


 void LDisorder::WalkFast(){
  //like LDisorder::Walk but optimized
  //no LMonitor connectivity and no verbose
  Init();
  SaveProperties(1);
  if (!fspace->IsInside(fe)) return;
  while (1){
    if (fN>=fmax) break;//out of time

    if (ft<=fe->time) {
      SaveProperties(1);
      continue;
    }

    if (!StepFast()) break;//out of space
  }
  ft-=ftime;//since the ft points to the time where should store next

}

inline
 int LDisorder::StepFast(){
  //like LDisorder:Step but optimized with no LMonitor connectivity
  //to use with WalkFast
  static long seed=initseed();
  int idx;

  int i0=fe->i-fWa2c;
  int j0=fe->j-fWb2c;
  int k0=fe->k-fWc2c;

  float e1=fe->GetEnergy();//current energy but no gradient
  float e2,r,sum=0.0;

  //calculate probabilities
  for(int i=1;i<=fWa;i++)
    for(int j=1;j<=fWb;j++)
      for(int k=1;k<=fWc;k++){
	idx=fWb*fWc*(i-1)+fWc*(j-1)+k;
	if (idx!=idxexclude){//exclude sumation over the site of origin

	  e2=fspace->GetEnergy(i0+i,j0+j,k0+k);
	  r=fspace->Distance(idx);//to speed up
	  //r=fspace->Distance(i-fWa2c,j-fWb2c,k-fWc2c);

	  sum+=Probability(r,e1,e2);
	  fprob[idx]=sum;
	} else fprob[idx]=sum;  //so will never be selected in Where
      }

  float a;a=sum*ran1(&seed);
  unsigned long jm,jl=0;
  idx=fNprob+1;
  if(a>fprob[fNprob-1]) idx=fNprob;
  else if(a<=fprob[1]) idx=1;
  else 
    while(idx-jl>1){
      jm=(idx+jl)>>1;//compute midpoint 
      if(a>fprob[jm]) jl=jm; else idx=jm;
    }
  int x,z,y;
  idx--;

  x=idx/(fWb*fWc);
  y=(idx-x*fWb*fWc)/fWc;
  z=idx-x*fWb*fWc-y*fWc+1;
  idx++;x++;y++;
  fe->time+=(expdev(&seed)*sum/(fprob[idx]-fprob[idx-1]));//no scale
  return MoveTo(x-fWa2c,y-fWb2c,z-fWc2c);//move particle
}

void
 LDisorder::Test(int testmode){
  //this is for testing only
  //tests are:
  //1 - particle is moving in k direction until it exits space
  //2 - particle is moving in -k direction until it exits space

  Init();
  std::cout<<"ft:"<<ft<<" t:"<<fe->time<<" k:"<<fe->k<<" e:"<<fe->GetEnergy()<<" fN"<<fN<<std::endl;
  SaveProperties(1);
  while(fspace->IsInside(fe)){
    if (fN>=fmax) break;
    switch (testmode){
    case 1: if (!MoveTo(0,0,1)) return;
      break;
    case 2: if (!MoveTo(0,0,-1)) return;
      break;
    }
    fe->time+=ftime;
    std::cout<<"ft:"<<ft<<" t:"<<fe->time<<" k:"<<fe->k<<" e:"<<fe->GetEnergy()<<" fN"<<fN<<std::endl;
    SaveProperties(1);
  }
}