/* applet No. 996
 *
 * adatom on metal 
 * - kinetic Monte-Carlo simulation
 *
 * Created by Ikeuchi Mitsuru on August 13 2006.
 * Copyright (c) 2006 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2006.08.13  created
 *
 */

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

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

/* -------------------------------------- set global ------ */

  Choice ch_view;
  Button bt_reset, bt_startStop;
  Scrollbar sc_nn, sc_qqhop, sc_qqnb, sc_kT;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int dispWidth = 630;
  int dispHeight = 420;

  int sleepTime = 50;

  double t = 0.0;
  double dt = 1.0;

  int resetSW = 0;
  int started = 1;
  int dispMode = 0;

  //double qqAdsorption = 0.0;
  double qqHopping = 0.2;
  double qqNeighbor = 0.3;
  double kT = 0.1;


  int NNx = 32;
  int NNy = 32;
  int NNa = 300;

  int lattice[][] = new int[NNx][NNy];
  int adatom[][] = new int[600][2];
  double transition[] = new double[600*4];
  double rate[] = new double[8];

/* ---------------------------------- applet control -------- */

  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.addItem("adatom");
    ch_view.addItem("potential");
    ch_view.addItemListener(this);
    ch_view.select("adatom");

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

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

    sc_nn= new Scrollbar(Scrollbar.HORIZONTAL,300,10,1,610);
    sc_nn.addAdjustmentListener(this);

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

    sc_qqnb= new Scrollbar(Scrollbar.HORIZONTAL,30,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(sc_nn);
    pnl.add(sc_qqhop);
    pnl.add(sc_qqnb);
    pnl.add(sc_kT);
    add(pnl,"North");

    setInitialCondition();
  }

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

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

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

  public void actionPerformed(ActionEvent ev){
    if(ev.getSource() == bt_reset){ 
      resetSW = 1;
    } else if (ev.getSource() == bt_startStop){
      if (started==0) {
        started = 1;
      } else {
        started = 0;
      }
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_nn) { 
      NNa = sc_nn.getValue();
	  resetSW = 1;
    } else if (ev.getSource() == sc_qqhop) { 
      qqHopping = 0.01*(double)(sc_qqhop.getValue());
    } else if (ev.getSource() == sc_qqnb) { 
      qqNeighbor = 0.01*(double)(sc_qqnb.getValue());
    } else if (ev.getSource() == sc_kT) { 
      kT = 0.001*(double)(sc_kT.getValue());
    }
  }

  public void run() {
    while (th != null) {
      try {
        timeEvolution();
        offPaint();
        repaint();
        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() {

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

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

    } else if (dispMode==1) {
      ;

    } else if (dispMode==2) {
      ;

    }

    gOff.setColor(Color.black);
    gOff.drawString("t="+(int)(t*1e14+0.5)/100.0+" ps",10,40);
    if (started==0) {
      gOff.drawString("stopped",dispWidth/6*1+10,40);
    } else {
      gOff.drawString("started",dispWidth/6*1+10,40);
	}
    gOff.drawString("N="+NNa+"",dispWidth/6*2+10,40);
    gOff.drawString("Qhop="+(int)(qqHopping*100.0+0.5)/100.0+" eV",dispWidth/6*3+10,40);
    gOff.drawString("Qnbr="+(int)(qqNeighbor*100.0+0.5)/100.0+" eV",dispWidth/6*4+10,40);
    gOff.drawString("kT="+(int)(kT*1000+0.5)/1000.0+" eV",dispWidth/6*5+10,40);
    gOff.drawString("  ="+(int)(kT*11600+0.5)+" K",dispWidth/6*5+10,60);

    //gOff.setColor(Color.blue);
    //gOff.drawString("nn="+numberOfAdatom()+"",dispWidth/6*2+10,60);
  }

/*-------------------------- dispLattice --------- */

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

    ir = span-2;
    for(i=0;i<NNx;i++) {
      for(j=0;j<NNy;j++) {
	    if (lattice[i][j]>=0) {
          gOff.setColor(Color.getHSBColor(0.1f,0.4f,0.8f));
		} else {
		  gOff.setColor(Color.getHSBColor(0.66f,0.4f,0.8f));
		}
        gOff.fillOval(gx+span*i-ir/2,gy+span*j-ir/2, ir, ir);
      }
	}
  }

/*-------------------------- statistics ------------- */

  int numberOfAdatom() {
    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);
  }

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

  void setInitialCondition() {
    int ia;

    t = 0.0;

    initLattice();
	for(ia=0;ia<NNa;ia++) {
      setAdatom(ia);
	}
  }

  void initLattice() {
    int i,j;

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

  void setAdatom(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);
	
    adatom[ia][0] = i; adatom[ia][1] = j;
    lattice[i][j] = ia;
  }

/* ----------------------------- timeEvolution -------------- */

  void timeEvolution() {
    if (resetSW==1) {
      setInitialCondition();
      resetSW = 0;
      started = 1;
    }

    if (started==1) {
      timeStep();
    }
  }

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

    nn = NNa/10 + 1;

    for(i=0;i<nn;i++) {
      deltat = kmcStep();
      t = t + deltat;
    }
  }

  double kmcStep() {
    int ip;
    double rsum,r;

    rsum = setAllTransition();
    ip = selectTransion(rsum);
    move(ip);

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

/* ---------- set transition , and calc. sum of transition ---------- */

  double setAllTransition() {
    int ia;
    double rsum;

    setRate();

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

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

    i = adatom[ia][0]; j = adatom[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);

    return ( transition[ia*4+0] + transition[ia*4+1] + transition[ia*4+2] + transition[ia*4+3] );
  }

  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);
  }

  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 );
  }


  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);
  }

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

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

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

/* ---------- select transion ---------- */

  int selectTransion(double rsum) {
    int ip,ipmax;
    double rnd;

    ipmax = 4*NNa;
    rnd = rsum*Math.random();

    for (ip=0; ip<ipmax; ip++) {
	  rnd -= transition[ip];
	  if (rnd<0.0) return(ip);
	}
    return(ipmax-1);
  }

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

  void move(int ip) {
    int ia,id,i,j;

    ia = ip/4;
    id = ip%4;
    i = adatom[ia][0]; j = adatom[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;
    adatom[ia][0] = i; adatom[ia][1] = j;
    lattice[i][j] = ia;
  }

/* ----------------------------- end of applet -------------- */

}

/* ----- 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)
 *
 * ----------
 *
 * transion rate k
 *  k = nue0*exp(-dE/kT)
 *  nue0 = kT/h , h: Plank constant
 *
 */