/** applet No. 1218
 *
 * sintering powder - kinetic Monte-Carlo Method 2D
 * - multi thread  KMC :: display - asynchronous
 *
 * Created by Ikeuchi Mitsuru on September 29 2007.
 * Copyright (c) 2007 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2007.09.29  created
 *
 */

import java.awt.*;
import java.awt.event.*;
import java.applet.*;

public class sinteringPowderMtKMC2D extends Applet
  implements ItemListener, ActionListener, AdjustmentListener, Runnable {

  // ------------------------------------ preset field -----------

  Thread th = null; // for run()-paint() thread
  KMC2D kmc = null;

 // for event
  Choice ch_view;
  Button bt_reset, bt_startStop;
  Scrollbar sc_add, sc_qqhop, sc_qqnb, sc_kT;

  // for off-paint buffer
  Dimension dim;
  Image imgOff;
  Graphics gOff;

  // display
  int dispWidth = 630;
  int dispHeight = 420;
  int sleepTime = 50;
  int dispMode = 0;
  int thCount = 0;

  // ------------------------------ applet main thread -----------

  public void init() {

    resize(dispWidth,dispHeight);
    setBackground(Color.white);

    dim = getSize();
    imgOff = createImage(dim.width,dim.height);
    gOff = imgOff.getGraphics();

    ch_view = new Choice();
    ch_view.add("0 yellow");
    ch_view.add("1 magenta");
    ch_view.add("2 green");
    ch_view.addItemListener(this);
    ch_view.select("1 magenta");

    bt_reset= new Button("reset");
    bt_reset.addActionListener(this);

    bt_startStop= new Button("start/stop");
    bt_startStop.addActionListener(this); 

    sc_add= new Scrollbar(Scrollbar.HORIZONTAL,20,10,5,110);
    sc_add.addAdjustmentListener(this);

    sc_qqhop= new Scrollbar(Scrollbar.HORIZONTAL,20,10,0,110);
    sc_qqhop.addAdjustmentListener(this);

    sc_qqnb= new Scrollbar(Scrollbar.HORIZONTAL,50,10,0,110);
    sc_qqnb.addAdjustmentListener(this);

    sc_kT= new Scrollbar(Scrollbar.HORIZONTAL,100,10,25,210);
    sc_kT.addAdjustmentListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,6,5,0));
    //pnl.add(new Label("  "));
    pnl.add(bt_reset);
    pnl.add(bt_startStop);
    pnl.add(new Label("  "));
    pnl.add(sc_qqhop);
    pnl.add(sc_qqnb);
    pnl.add(sc_kT);
    add(pnl,"North");
  }

  public void start() {
    if (kmc == null) {
      kmc = new KMC2D();
      kmc.start();
    }
    if (th == null) {
      th = new Thread(this);
      th.start();
    }
  }

  public void stop() {
    if (kmc != null) kmc = null;
    if (th != null) th = null;
  }


  // ---------------------------------- event listener -----------

  public void itemStateChanged(ItemEvent ev){
    if (ev.getSource() == ch_view) {
      dispMode = ch_view.getSelectedIndex();
    }
  }

  public void actionPerformed(ActionEvent ev){
    if(ev.getSource() == bt_reset){ 
      kmc.reset();
      thCount = 0;
    } else if (ev.getSource() == bt_startStop){
      if (kmc.getStartSW()==0) {
        kmc.setStartSW(1);
      } else {
        kmc.setStartSW(0);
      }
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_add) { 
      kmc.setAddingRate( 0.001*(double)(sc_add.getValue()) );
	  
    } else if (ev.getSource() == sc_qqhop) { 
      kmc.setQHopping( 0.01*(double)(sc_qqhop.getValue()) );

    } else if (ev.getSource() == sc_qqnb) { 
      kmc.setQNeighbor( 0.01*(double)(sc_qqnb.getValue()) );

    } else if (ev.getSource() == sc_kT) { 
      kmc.setKT( 0.001*(double)(sc_kT.getValue()) );

    }
  }



  // ========================= run() - paint() loop ==============

  public void run() {
    while (th != null) {
      thCount += 1;
      if (kmc != null) offPaint();
      repaint();
      try {
        Thread.sleep(sleepTime);
      } catch (InterruptedException e) { }
    }
  }

  public void update(Graphics g) {
    paint(g);
  }

  public void paint(Graphics g) {
    g.drawImage(imgOff,0,0,this);
  }


  // ---------------------------------------- offPaint -----------

  void offPaint() {
    int NNx, NNy;

    NNx = kmc.getNNx();
    NNy = kmc.getNNy();

    gOff.setColor(Color.white);
    gOff.fillRect(0,0,dim.width,dim.height);

    if (dispMode==0) {
      dispLattice(30, 80, 3);

    } else if (dispMode==1) {
      ;

    } else if (dispMode==2) {
      ;

    }

    gOff.setColor(Color.gray);
    gOff.drawRect(30,80,3*NNx,3*NNy);

    gOff.setColor(Color.black);
    gOff.drawString("t="+(int)(kmc.getTime()*1e12+0.5)/1000.0+" ns",10,40);
    if (kmc.getStartSW()==1) {
      gOff.drawString("start",dispWidth/6*1+20,40);
    } else {
      gOff.drawString("stop",dispWidth/6*1+20,40);
    }
    gOff.drawString("N="+kmc.getNumberOfAtom()+"",dispWidth/6*2+10,40);
    gOff.drawString("Qhop="+(int)(kmc.getQHopping()*100.0+0.5)/100.0+" eV",dispWidth/6*3+10,40);
    gOff.drawString("Qnbr="+(int)(kmc.getQNeighbor()*100.0+0.5)/100.0+" eV",dispWidth/6*4+10,40);
    gOff.drawString("kT="+(int)(kmc.getKT()*1000+0.5)/1000.0+" eV",dispWidth/6*5+10,40);
    gOff.drawString("  ="+(int)(kmc.getKT()*11600+0.5)+" K",dispWidth/6*5+10,60);

    gOff.setColor(Color.black);
    gOff.drawString("thread KMC = "+kmc.getCount()+" ",400,380);
    gOff.drawString("thread disp = "+thCount+" ",400,400);
  }


  // ------------------------------------ plot methods -----------

  void dispLattice(int gx, int gy, int span) {
    int i,j,ir,ia;
    double col;
    int NNx, NNy;

    NNx = kmc.getNNx();
    NNy = kmc.getNNy();

    ir = span;
    for(i=0;i<NNx;i++) {
      for(j=0;j<NNy;j++) {
        ia = kmc.getLattice(i,j);
        if (ia>=0) {
          if (kmc.getKind(ia)==0) {
            gOff.setColor(Color.getHSBColor(0.1f,0.6f,0.8f));
          } else {
            gOff.setColor(Color.getHSBColor(0.05f,0.6f,0.8f));
          }
          gOff.fillRect(gx+span*i,gy+span*j, ir, ir);
        }
      }
    }
  }

}



// ================== kinetic Monte-Carlo Method 2D ==============

class KMC2D extends Thread {

  private double t = 0.0;
  private double dt = 1.0;

  private double qqHopping = 0.2;
  private double qqNeighbor = 0.5;
  private double kT = 0.10;
  private double addingRate = 0.02;
  private double totalRate = 1.0e15;

  private int NNx = 100;
  private int NNy = 100;
  private int NNa = 0;
  private int maxNNa = 20000;
  private int NNN = 20000;

  private int lattice[][] = new int[NNx][NNy];
  private int atom[][] = new int[NNN][2];
  private int kind[] = new int[NNN];
  private double transition[] = new double[NNN*4];
  private double transitionAtom[] = new double[NNN];
  private double rate[] = new double[8];

  private int resetSW = 0;
  private int startSW = 1;
  private int stepSW = 1;
  private int changeFlag = 1;
  private int sleepTime = 1;
  private int stopSleepTime = 100;
  private int count = 0;


  // ------------------------------- class constructor -----------

  KMC2D() {
    setInitialCondition();
  }


  // -------------------------------------- thread run -----------

  public void run() {
    while (true) {
      count += 1;
      timeEvolution();
      if (sleepTime>0) {
        try {
          Thread.sleep(sleepTime);
        } catch (InterruptedException e) { }
      }
    }
  }


  // --------------------------- set initial condition -----------

  void setInitialCondition() {
    int ia,i,j;

    t = 0.0;
    count = 0;

    initLattice();

    ia = 0;
    for (i=0; i<3; i++) {
      for (j=0; j<3; j++) {
        ia = setPowder(ia, 20+29*i, 20+29*j, 14);
      }
	}
    NNa = ia;

    setAllTransition();
    changeFlag = 0;
  }

  void initLattice() {
    int i,j;

    for(i=0;i<NNx;i++) {
      for(j=0;j<NNy;j++) {
        lattice[i][j] = -1;
      }
    }
  }

  int setPowder(int ii, int ix, int iy, int nr) {
    int i,j,ia;

    ia = ii;
    for(i=ix-nr;i<=ix+nr;i++) {
      for(j=iy-nr;j<=iy+nr;j++) {
	    if ((i-ix)*(i-ix)+(j-iy)*(j-iy)<=nr*nr) {
          kind[ia] = 0;
          atom[ia][0] = i; atom[ia][1] = j;
          lattice[i][j] = ia;
          ia += 1;
		}
      }
    }
    return (ia);
  }


  void setAtom(int ia) {
    int i,j;

    do {
      i = (int)(0.99999999*Math.random()*NNx);
      j = (int)(0.99999999*Math.random()*NNy);
    } while (lattice[i][j]>= 0);
	
    atom[ia][0] = i; atom[ia][1] = j;
    lattice[i][j] = ia;
  }

  //

  void addAnAtom(int mxNNa) {
    int ia;

    if (NNa<mxNNa) {
      ia = NNa;
      setAtom(ia);
      NNa += 1;

      kind[ia] = 1;
      totalRate = setAllTransition();
      changeFlag = 0;
    }
  }


  // ---------------------------------- time evolution -----------

  private void timeEvolution(){
    int i;

    if (resetSW==1) {
      setInitialCondition();
      resetSW = 0;
    }

    if (startSW==1) {
      timeStep();

    } else {
      if (stepSW==1) {
        timeStep();
        stepSW = 0;
      }
      try {
        Thread.sleep(stopSleepTime);
      } catch (InterruptedException e) { }
    }
  }

  private void timeStep() {
    int i,nn;
    double deltat;

    //if (Math.random()<addingRate) addAnAtom(maxNNa);

    //nn = NNa/10+1; if (nn>1000) nn = 1000;
    nn=50;
    for(i=0;i<nn;i++) {
      deltat = kmcStep();
      t = t + deltat;
    }
    correctTotalRate();
  }

  private void correctTotalRate() {
    int i;
    double s;

    s = 0.0;
    for(i=0;i<NNa;i++) {
      s += transitionAtom[i];
    }
    totalRate = s;
  }

  private double kmcStep() {
    int i,ip;
    double r;

    if (changeFlag==1) {
      totalRate = setAllTransition();
      changeFlag = 0;
    }

    ip = selectTransion(totalRate);
    move(ip);

    r = (1.0-1.0e-12)*Math.random()+1.0e-12;
    return(-Math.log(r)/totalRate);
  }


  // set transition , and calc. sum of transition 

  private double setAllTransition() {
    int ia;
    double rsum;

    setRate();

    rsum = 0.0;
    for(ia=0;ia<NNa;ia++) {
      rsum += setTransition(ia);
    }
    return (rsum);
  }

  private double setTransition(int ia) {
    int i,j,nb;

    i = atom[ia][0]; j = atom[ia][1];
    nb = noOfNeighbor(i,j);

    transition[ia*4+0] = dirRate((i+1)%NNx, j, nb);
    transition[ia*4+1] = dirRate((i-1+NNx)%NNx, j, nb);
    transition[ia*4+2] = dirRate(i, (j+1)%NNy, nb);
    transition[ia*4+3] = dirRate(i, (j-1+NNy)%NNy, nb);
    transitionAtom[ia] = transition[ia*4+0]+transition[ia*4+1]+transition[ia*4+2]+transition[ia*4+3];
    return (transitionAtom[ia]);
  }

  private double dirRate(int ii, int jj, int nb) {
    int nb1;
    double rt;

    rt = 0.0;
    if (lattice[ii][jj]>=0) {
      rt = 0.0;
    } else {
      nb1 = noOfNeighbor(ii,jj)-1;
      if (nb<=nb1) {
        rt = rate[0];
      } else {
        rt = rate[nb-nb1];
      }
    }
    return(rt);
  }

  private int noOfNeighbor(int i,int j) {
    int nb;

    nb = 0;
    if (lattice[(i+1)%NNx][j]>=0) nb += 1;
    if (lattice[(i-1+NNx)%NNx][j]>=0) nb += 1;
    if (lattice[i][(j+1)%NNy]>=0) nb += 1;
    if (lattice[i][(j-1+NNy)%NNy]>=0) nb += 1;
    return( nb );
  }

  private void setRate() {
    rate[0] = transitionRate(qqHopping);
    rate[1] = transitionRate(qqHopping + qqNeighbor);
    rate[2] = transitionRate(qqHopping + 2.0*qqNeighbor);
    rate[3] = transitionRate(qqHopping + 3.0*qqNeighbor);
    rate[4] = transitionRate(qqHopping + 4.0*qqNeighbor);
  }

  private double transitionRate(double energy) {
    double h,nue0;

    h = 6.626e-34;
    nue0 = kT*1.602e-19/h;

    return( nue0*Math.exp(-energy/kT) );
  }

  /* ----- hermonic transition state theory -----

  double transitionRate(double energy) {
    double m,r,omega,nue,rate;

    m = 137.0*1.67e-27;
    r = 0.3*1.0e-9;
    omega = Math.sqrt(energy*1.6e-19/(m*r*r));
    nue = omega/(2.0*3.1416);//s

    rate = nue*Math.exp(-energy/kT);
    return( rate );
  }

  */


  // select transion

  private int selectTransion(double rsum) {
    int ia,iap;
    double rnd;

    rnd = rsum*Math.random();

    for (ia=0; ia<NNa; ia++) {
      rnd -= transitionAtom[ia];
      if (rnd<0.0) {
        rnd += transitionAtom[ia];
        for (iap=0; iap<4; iap++) {
          rnd -= transition[ia*4+iap];
          if (rnd<0.0) return (ia*4+iap);
        }
        return(ia*4+3);
      }
    }
    return((NNa-1)*4+3);
  }

  // move (update adatom[][],lattice[][])

  private void move(int ip) {
    int ia,id,i0,j0,i,j;
    double rs0,rs1;

    if (transition[ip]>0.0) {
      ia = ip/4;
      id = ip%4;
      i = atom[ia][0]; j = atom[ia][1];

      lattice[i][j] = -1;
      if (id==0) i = (i+1)%NNx;
      if (id==1) i = (i-1+NNx)%NNx;
      if (id==2) j = (j+1)%NNy;
      if (id==3) j = (j-1+NNy)%NNy;

      atom[ia][0] = i; atom[ia][1] = j;
      lattice[i][j] = ia;

      rs0 = rateAround(i,j);
      rs1 = setRateAround(i,j);
      totalRate = totalRate - rs0 + rs1;
    }
  }

  private double setRateAround(int i0, int j0) {
    int i,j,ii,jj;
    double srate;

    srate = 0.0;
    for (i=i0-2; i<=i0+2; i++) {
      ii = (i+NNx)%NNx;
      for (j=j0-2; j<=j0+2; j++) {
        jj =  (j+NNy)%NNy;
        if (lattice[ii][jj]>=0) {
          srate += setTransition(lattice[ii][jj]);
        }
      }
    }
    return ( srate );
  }

  private double rateAround(int i0, int j0) {
    int i,j,ii,jj;
    double srate;

    srate = 0.0;
    for (i=i0-2; i<=i0+2; i++) {
      ii = (i+NNx)%NNx;
      for (j=j0-2; j<=j0+2; j++) {
        jj =  (j+NNy)%NNy;
        if (lattice[ii][jj]>=0) {
          srate += transitionAtom[lattice[ii][jj]];
        }
      }
    }
    return ( srate );
  }


  // --------------------------------------- utilities -----------

  private int numberOfAtom() {
    int i,j,n;

    n = 0;
    for(i=0;i<NNx;i++) {
      for(j=0;j<NNy;j++) {
        if (lattice[i][j]>=0) n += 1;
      }
    }
    return(n);
  }


  // ------------------------------------- I/O methods -----------

  // control
  public void reset() { resetSW = 1; }
  public void setStartSW(int sw) { startSW = sw; }
  public int getStartSW() { return startSW; }
  public void setStepSW() { stepSW = 1; }
  public int getCount() { return count; }

  public void setAddingRate(double ar) { addingRate = ar; changeFlag = 1; }
  public double getAddingRate() { return addingRate; }
  public void setQHopping(double qh) { qqHopping = qh; changeFlag = 1; }
  public double getQHopping() { return qqHopping; }
  public void setQNeighbor(double qn) { qqNeighbor = qn; changeFlag = 1; }
  public double getQNeighbor() { return qqNeighbor; }
  public void setKT(double ev) { kT = ev; changeFlag = 1; }
  public double getKT() { return kT; }

  // get data
  public int getNNx() { return NNx; }
  public int getNNy() { return NNy; }
  public double getTime() { return t; }
  public int getNumberOfAtom() { return numberOfAtom(); }
  public int getLattice(int i, int j) { return lattice[i][j]; }
  public int getKind(int ia) { return kind[ia]; }

}


/* ----- kinetic Monte-Carlo algorithm -----
 * ( or the Bortz-Kalos-Liebowitz (BKL) algorithm )
 *
 * i-th transition rate in the system : ri
 * total transition rate : RN = sum(ri, {i=1,2,...,N})
 *
 * 1) t = 0
 *
 * 2) select i : randomly (proportional to ri/RN)
 *  calculate R(i) = sum(rj, {j=1,2,...,i})
 *  get a uniform random number u in [0,1]
 *  find i  R(i-1) < u RN <= R(i)
 *
 * 3) execute i-th transition
 *
 * 4) update the time
 *  get a uniform random numner u in (0,1]
 *  dt = -log(u)/R
 *  t = t + dt
 *
 * 5) goto 2)
 *
 * ---------- transition state -----
 *
 * transion rate k
 *  k = nue0*exp(-dE/kT)
 *  nue0 = kT/h , h: Plank constant
 *
 *
 * ----- hermonic transition state theory -----
 *
 * oscillation in hermonic potential 
 *  U(x) = k x^2, k = V/r^2 (<- U(r) = V)
 *
 * motion equation m a = - kx 
 *   omega = sqrt(k/m) = sqrt( V/(m*r^2) )
 *
 * transion rate = (omega/(2 pai))*exp(-V/kT)
 *
 * m = 137*1.67e-27; // Ba
 * V = 0.3eV = 0.3*1.6e-19;
 * r = 0.3e-9;
 *
 * omega/(2 pai) = 2.4e11 (1/s)
 *
 *
 */