/* harmonic oscillator - diffusion Monte-Carlo method 2D */
/* coded by Ikeuchi Mitsuru   */
/* ver 0.0.1   2005.05.28     */

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

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

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

  Choice ch_view;
  Button bt_reset, bt_start, bt_stop;
  Scrollbar sc_val;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  Random rand = new Random();

  int sleepTime = 50;

  int dispMode = 0;
  int resetSW = 0;
  int started = 1;
  double kConstant = 0.5;


  double deltaTau = 0.1;
  double sqrtdt = Math.sqrt(deltaTau);
  double hbar = 1.0;

  double t = 0.0;
  double dt = 0.1;

  int NNr = 2000;
  int repStatus[] = new int[NNr];
  double replica[][] = new double[NNr][3];

  double eeRef0 = 0.0;
  double eeRef1 = eeRef0;
  double eeRef = eeRef0;

  int NN0 = 1000;
  int NN1 = NN0;

  int NN00 = NN0;

  int NNbx = 40;
  int NNby = 40;
  int box[][] = new int[NNbx][NNby];
  double meanBox[][] = new double[NNbx][NNby];
  double psi[][] = new double[NNbx][NNby];


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

  public void init() {

    resize(630,320);
    setBackground(Color.white);

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

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

    bt_start= new Button("start");
    bt_start.addActionListener(this); 

    bt_stop = new Button("stop");
    bt_stop.addActionListener(this);

    ch_view = new Choice();
    ch_view.addItem("replicas");
    ch_view.addItem("psi");
    ch_view.addItemListener(this);
    ch_view.select("replicas");

    sc_val= new Scrollbar(Scrollbar.HORIZONTAL,50,10,50,210);
    sc_val.addAdjustmentListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,6,5,0));
	pnl.add(bt_reset);
    pnl.add(bt_start);
    pnl.add(bt_stop);
    pnl.add(sc_val);
    pnl.add(new Label("  "));
    pnl.add(ch_view);
    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;
      started = 0;
    } else if(ev.getSource() == bt_start){ 
      started = 1;
    } else if(ev.getSource() == bt_stop){
      started = 0;
    }
		}

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if(ev.getSource() == sc_val) { 
      kConstant = 0.01*(double)(sc_val.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() {
    int gx,gy,gb;

    gx = 30; gy = 60;

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


    if (dispMode==0) {
      plotReplicas(gx,gy);
	} else if (dispMode==1) {
	  setPsi();
      plotMatrix(psi,gx, gy, 5.0, 0.15);
	}  else if (dispMode==2) {
	  ;
	} else if (dispMode==3) {
      ;
	}


    gOff.setColor(Color.black);
    gOff.drawString("t="+(int)(t*10+0.5)/10.0+" ",630/6*0+10,40);
    gOff.drawString("k="+kConstant+"",630/6*3+10,40);
    gOff.drawString("view",630/6*5+10,40);

    gb = 460;
    gOff.drawString("Nmax="+NNr+"",gb,100);
    gOff.drawString("d tau="+(float)(deltaTau)+"",gb,120);
    gOff.drawString("V(x)=k*(x^2+y^2)",gb,140);
	
    gOff.drawString("ER="+(float)(eeRef0)+"",gb,200);
    gOff.drawString("<ER>="+(float)(eeRef)+"",gb,220);
	gOff.drawString("N="+NN0+"",gb,240);

	gOff.drawString("E=1.0 at k=0.5",gb,300);
  }

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

  void plotReplicas(int gx, int gy) {
    int i,ix,iy;
    double x,sc,p1,p2;

    gOff.setColor(Color.black);
    gOff.fillRect(gx,gy,200,200);

    gOff.setColor(Color.green);
    for (i=0; i<NNr; i++) {
      if (repStatus[i]==1) {
        ix = 100+(int)(20.0*replica[i][0]+0.5);
        iy = 100+(int)(20.0*replica[i][1]+0.5);
		if (ix>=0 && ix<=200 && iy>=0 && iy<=200) {
          gOff.drawOval(gx+ix-2,gy+iy-2, 4, 4);
		}
	  }
    }
  }

  void plotMatrix(double mat[][], int gx, int gy, double mag, double col) {
    int i,j,ix,iy,ixs,iys;
    double f;

	ixs = 200/NNbx; iys = 200/NNby;
    for (i=0; i<NNbx; i++) {
      for (j=0; j<NNby; j++) {
        ix = gx+i*ixs;
        iy = gy+j*iys;
		f = mag*mat[i][j]; if (f>0.999) f=0.999;
		if (ix>=0 && ix<=300 && iy>=0 && iy<=300) {
		  gOff.setColor(Color.getHSBColor((float)col, 0.9f,(float)f));
          gOff.fillRect(ix,iy, ixs, iys);
		}
	  }
    }
  }

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

  void setInitialCondition() {
    int i,n0;

    t = 0.0;
    for (i=0; i<NNr; i++) {
	  repStatus[i] = 0;
	}

    n0 = NNr/2;
    for (i=0; i<n0; i++) {
	  repStatus[i] = 1;
	  replica[i][0] = 4.0*(Math.random()-0.5);
	  replica[i][1] = 4.0*(Math.random()-0.5);
	}
    setBox(0);

    NN0 = numberOfReplica();
	NN00 = NN0;
	eeRef0 = averagedPotential();
  }

/* ----------------------------- setBox/ psi -------------- */

  void setBox(int sw) {
    int i,ix,iy;

    for (ix=0; ix<NNbx; ix++) {
      for (iy=0; iy<NNby; iy++) {
	    box[ix][iy] = 0;
	  }
	}

    for (i=0; i<NNr; i++) {
	  if (repStatus[i]==1) {
	    ix = (int)(NNbx/2+replica[i][0]*NNbx/10.0+0.5);
		if (ix<0) ix=0;
		if (ix>NNbx-1) ix=NNbx-1;
		iy = (int)(NNby/2+replica[i][1]*NNby/10.0+0.5);
		if (iy<0) iy=0;
		if (iy>NNby-1) iy=NNby-1;
	    box[ix][iy] += 1;
	  }
	}

    if (sw==1) {
      for (ix=0; ix<NNbx; ix++) {
        for (iy=0; iy<NNby; iy++) {
	      meanBox[ix][iy] = 0.99*meanBox[ix][iy]+0.01*box[ix][iy];
	    }
	  }
	} else {
      for (ix=0; ix<NNbx; ix++) {
        for (iy=0; iy<NNby; iy++) {
	      meanBox[ix][iy] = box[ix][iy];
		}
	  }
	}
  }

  void setPsi() {
    int i,j;
    double s,s2;

    s2 = 0;
	for (i=0; i<NNbx; i++) {
	  for (j=0; j<NNby; j++) {
	    s2 += meanBox[i][j]*meanBox[i][j];
	  }
	}
    s = Math.sqrt(s2);

	for (i=0; i<NNbx; i++) {
	  for (j=0; j<NNby; j++) {
	    psi[i][j] = meanBox[i][j]/s;
	  }
	}
  }

  double averagedPotential() {
    int i,n;
    double s;

    s = 0.0; n = 0;
    for (i=0; i<NNr; i++) {
	  if (repStatus[i]==1) {
	    n += 1;
        s += potential(replica[i][0],replica[i][1]);
	  }
	}
	return(s/n);
  }

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

  void timeEvolution() {
    int i;

    if (resetSW==1) {
	  setInitialCondition();
	  resetSW = 0;
	}
    if (started==1) {
      t = t + dt;
      timeStep();
    }
  }

  void timeStep() {
    int i;

    t = t + dt;
	for (i=0; i<NNr; i++) {
	  if (repStatus[i]==1) {
        replica[i][0] += sqrtdt*gaussianRandom();
		replica[i][1] += sqrtdt*gaussianRandom();
	  }
	}

    for (i=0; i<NNr; i++) {
	  if (repStatus[i]==1) {
        birthDeathProcess(i);
	  }
	}
    for (i=0; i<NNr; i++) {
	  if (repStatus[i]==2) {
        repStatus[i] = 1;
	  }
	}

    NN1 = numberOfReplica();
    eeRef1 = eeRef0 + hbar/deltaTau*(1.0-(NN1+1.0)/(NN0+1.0)*(NN1+10.0*NNr)/(NN00+10.0*NNr));
	NN0 = NN1;
	eeRef0 = eeRef1;
	eeRef = 0.99*eeRef+0.01*eeRef0;

    setBox(1);
  }


  void birthDeathProcess(int i) {
    int mn;
    double xi,yi;

    xi = replica[i][0]; yi = replica[i][1];
    mn = (int)(weightFunction(xi,yi)+Math.random());
	if (mn>3) mn = 3;

    if (mn==0) {
      repStatus[i] = 0;
	} else if (mn==2) {
      birthReplica(i);
	} else if (mn==3) {
      birthReplica(i);
	  birthReplica(i);
	}
  }


  void birthReplica(int i) {
    int ip;

    ip = seekSpace(i);
	if (ip>=0) {
	  repStatus[ip] = 2;
	  replica[ip][0] = replica[i][0];
	  replica[ip][1] = replica[i][1];
    }
  }

  int seekSpace(int istart) {
    int i,ip;

    ip = -1;
	for (i=1; i<NNr; i++) {
	  if (repStatus[(i+istart)%NNr]==0) {
	    ip = (i+istart)%NNr;
		break;
	  }
	}
	return(ip);
  }

/* ----------------------------- utility methods -------------- */

  int numberOfReplica() {
    int i,n;
	
    n = 0;
	for (i=0; i<NNr; i++) {
	  if (repStatus[i]==1) n+=1;
	}
	return(n);
  }

  double weightFunction(double xn,double yn) {
    return(Math.exp(-(potential(xn,yn)-eeRef0)*deltaTau/hbar));
  }

  double gaussianRandom() {
    return(rand.nextGaussian());
  }

  double potential(double x, double y) {
    return(kConstant*(x*x+y*y));
  }


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

}