/* harmonic oscillator - diffusion Monte-Carlo method 1D */
/* 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 harmonicOscillatorDMC1D extends Applet
  implements ActionListener, AdjustmentListener, Runnable {

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

  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 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 = 500;
  int NN1 = NN0;

  int NN00 = NN0;

  int NNb = 100;
  int box[] = new int[NNb];
  double meanBox[] = new double[NNb];
  double psi[] = new double[NNb];


/* ----------------------------- 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);
    
    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(new Label("  "));
    pnl.add(new Label("  "));
    pnl.add(sc_val);
    add(pnl,"North");

    setInitialCondition();
    repaint();
  }

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

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

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

    plotReplica(gx,gy,10);
    plotMeanBox(gx,gy,2.0);
    drawPotential(gx,gy,10.0);
/*
    plotPsi(gx,gy, 500.0);
*/

    gOff.setColor(Color.black);
    gOff.drawString("k="+kConstant+"",630/6*5+10,40);
    gOff.drawString("t="+(int)(t*10+0.5)/10.0+" ",630/6*0+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",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 = 0.5au at k=0.5",gb,280);

  }

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

  void drawPotential(int gx, int gy, double mag) {
    int ix;
    double x,sc,p1,p2;

    gOff.setColor(Color.red);
    sc = 5.0/200.0;
    for (ix=0; ix<400-1; ix++) {
      p1 = mag*potential(sc*(ix-200));
      p2 = mag*potential(sc*(ix+1-200));
      gOff.drawLine(gx+ix,gy+200-(int)(p1), gx+ix+1, gy+200-(int)(p2));
    }
  }

  void plotReplica(int gx, int gy, int iy) {
    int i,ix;

    for (i=0; i<NNr; i++) {
      if (repStatus[i]==1) {
        ix = gx+(int)(200.0+replica[i][0]*(200.0/5.0)+0.5);
        gOff.setColor(Color.blue);
        gOff.drawLine(ix,gy+200+2,ix,gy+200+2+iy);
	  }
    }
  }

  void plotBox(int gx, int gy, double mag) {
    int i,iy,a;

    a = 400/NNb;
    gOff.setColor(Color.blue);
    for (i=0; i<NNb; i++) {
      iy = (int)(mag*box[i]+0.5);
      if (box[i]>0) {
        gOff.fillRect(gx+i*a-a/2,gy+200-iy,a,iy);
      }
    }
  }

  void plotMeanBox(int gx, int gy, double mag) {
    int i,iy,a;

    a = 400/NNb;
    gOff.setColor(Color.green);
    for (i=0; i<NNb; i++) {
      iy = (int)(mag*meanBox[i]+0.5);
      if (meanBox[i]>0.5) {
        gOff.fillRect(gx+i*a-a/2,gy+200-iy,a,iy);
      }
    }
  }

  void plotPsi(int gx, int gy, double mag) {
    int ix,a;
    double x,sc,p1,p2;

    setPsi();

    a = 400/NNb;
    gOff.setColor(Color.orange);
    for (ix=0; ix<NNb-1; ix++) {
      p1 = mag*psi[ix];
      p2 = mag*psi[ix+1];
      gOff.drawLine(gx+a*ix,gy+200-(int)(p1), gx+a*(ix+1), gy+200-(int)(p2));
    }
  }

/*-------------------------- 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);
    }
    setBox(0);

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

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

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

    for (ix=0; ix<NNb; ix++) {
	  box[ix] = 0;
	}

    for (i=0; i<NNr; i++) {
	  if (repStatus[i]==1) {
	    ix = (int)(50.0+replica[i][0]*10+0.5);
		if (ix<0) ix=0;
		if (ix>NNb-1) ix=NNb-1;
	    box[ix] += 1;
	  }
	}

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

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

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

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

/* ----------------------------- 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();
      }
    }

    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;

    xi = replica[i][0];
    mn = (int)(weightFunction(xi)+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];
    }
  }

  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) {
    return(Math.exp(-(potential(xn)-eeRef0)*deltaTau/hbar));
  }

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

/*
  double potential(double x) {
    double v;

    v = 0.0;
    if (x<-1.5 || x>1.5) v = 15.0;

    return(v);
  }
*/

  double potential(double x) {
    return(kConstant*x*x);
  }


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

}