/** applet No. 1034
 *
 * H2 molecule - diffusion Monte-Carlo method
 * - multi thread  DMC :: display - asynchronous
 *
 * Created by Ikeuchi Mitsuru on October 07 2006.
 * Copyright (c) 2006-2007 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2006.10.07  created
 * ver 0.0.2   2006.11.30  improved code
 * ver 0.0.3   2006.12.19  bug in densityPlot() fixed
 * ver 0.0.4   2007.05.30  improved code
 *
 */

import java.awt.*;
import java.awt.event.*;
import java.applet.*;
import java.util.*; // Random class

public class H2MtDMC extends Applet
  implements MouseListener, MouseMotionListener,
    ItemListener, ActionListener, AdjustmentListener, Runnable {

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

  Thread th = null; // for run()-paint() thread
  DiffusionMC dmc = null;
  DrawGraph3d dp3d = new DrawGraph3d( 20.0, 20.0, 20.0 );

  // for event
  Choice ch_average, ch_view;
  Button bt_reset, bt_startStop, bt_step;
  Scrollbar sc_val;

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

  int sleepTime = 50;

  int dispMode = 0;

  int dgX, dgY, dgXb, dgYb; // mouse 

  int thCount = 0;


  // -------------------------------- applet main loop -----------

  public void init() {

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

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

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

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

    bt_step= new Button("step");
    bt_step.addActionListener(this); 

    ch_average = new Choice();
    ch_average.add("100");
    ch_average.add("1000");
    ch_average.add("10000");
    ch_average.addItemListener(this);
    ch_average.select("100");

    ch_view = new Choice();
    ch_view.add("replicas");
    ch_view.add("density");
    ch_view.add("density/2");
    ch_view.addItemListener(this);
    ch_view.select("replicas");

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

    addMouseListener(this);  
    addMouseMotionListener(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(bt_step);
    pnl.add(sc_val);
    pnl.add(ch_average);
    pnl.add(ch_view);
    add(pnl,"North");

  }

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

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


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

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

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

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if(ev.getSource() == sc_val) { 
      dmc.setppDistance( 0.01*(double)(sc_val.getValue()) );
    }
  }

  public void mousePressed(MouseEvent ev){ }
  public void mouseReleased(MouseEvent ev){
    dgXb = 0; dgYb = 0;
    dgX = 0; dgY = 0;
  }
  public void mouseClicked(MouseEvent ev){ }
  public void mouseEntered(MouseEvent ev){ } 
  public void mouseExited (MouseEvent ev){ }

  public void mouseMoved  (MouseEvent ev){ }

  public void mouseDragged(MouseEvent ev){
    dgXb = dgX; dgYb = dgY;
    dgX=ev.getX();
    dgY=ev.getY();
    if (dgXb==0 && dgYb==0) {
      dgXb = dgX; dgYb = dgY;
    }
    dp3d.setViewDelta( (dgX-dgXb), (dgY-dgYb) );
  }



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

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

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

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

  public void offPaint() {
    if (dmc != null) {
      dp3d.plotGraph(gOff, dmc, dispMode, thCount);
    }
  }

}



// =============================== draw graph 3D class ===========

class DrawGraph3d {

  private Graphics gOff;
  private DiffusionMC dmc;

  private int xImageLoc = 60;
  private int yImageLoc = 100;

  private double viewTheta = -15.0*3.14/180.0;
  private double viewFai = -72.0*3.14/180.0;
  private double dtheta = 0.0*3.14/180.0;
  private double pai = 3.1415926536;

  private double cosTh = Math.cos(viewTheta); 
  private double sinTh = Math.sin(viewTheta);
  private double cosFi = Math.cos(viewFai);
  private double sinFi = Math.sin(viewFai);

  private int nMax = 0;
  private double xMax = 20.0;
  private double yMax = 20.0;
  private double zMax = 20.0;
  private double cxx = xMax/2.0;
  private double cyy = yMax/2.0;
  private double czz = zMax/2.0;

  private double dispScale = 250.0/xMax;
  private double viewScale = 1.0;
  private double scale = dispScale*viewScale;

  private int NN = 10000;

  private short pID[] = new short[NN];

  private double px[] = new double[NN]; 
  private double py[] = new double[NN]; 
  private double pz[] = new double[NN]; 

  private int srtz[][] = new int[2][NN];

  private double wkx[] = new double[8]; 
  private double wky[] = new double[8]; 
  private double wkz[] = new double[8]; 

  private double pwkx[] = new double[8]; 
  private double pwky[] = new double[8]; 
  private double pwkz[] = new double[8]; 

  private int boxp[][] = {
   {0,1},{0,2},{0,4},{1,3},{1,5},{2,3},{2,6},{4,5},{4,6},{3,7},{5,7},{6,7}
  };

  private double particleColor[] = { 0.66, 0.33, 0.01, 0.05 };

  private int drawnCount = 0;


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

  DrawGraph3d(double xMaxSize, double yMaxSize, double zMaxSize) {
    xMax = xMaxSize; yMax = yMaxSize; zMax = zMaxSize;
    cxx = xMax/2.0; cyy = yMax/2.0; czz = zMax/2.0;
  }


  // ----------------------------------- class methods -----------

  public void plotGraph(Graphics gg, DiffusionMC dmc3d, int dispMode, int thCount) {
    int dmcCount,gb;

    gOff = gg;
    dmc = dmc3d;
    if (dmc==null) return;
    dmcCount = dmc.getCount();
    if (dmcCount==drawnCount) return;
    drawnCount = dmcCount;

    gOff.setColor(Color.white);
    gOff.fillRect(0,0,630,420);

    if (dispMode==0) {
      replicaPlot(10.0);

    } else if (dispMode==1) {
      densityPlot(0);

    } else if (dispMode==2) {
      densityPlot(1);

    } else if (dispMode==3) {
      ;
    }

    gOff.setColor(Color.black);
    gOff.drawString("tau="+(int)(dmc.getTime()*10+0.5)/10.0+" ",630/6*0+10,40);
    if (dmc.getStartSW()==0) {
      gOff.drawString("stop",630/6*1+10,40);
    } else {
      gOff.drawString("start",630/6*1+10,40);
    }
	gOff.drawString("R="+(int)(dmc.getppDistance()*100.0+0.5)/100.0+" au",630/6*3+10,40);
    gOff.drawString("average N",630/6*4+10,40);
    //gOff.drawString("view",630/6*5+10,40);

    gb = 400;
    gOff.drawString("Box = 20 x 20 x 20 au",gb,80);
    gOff.drawString("Nmax="+dmc.getMaxOfReplica()+"",gb,100);
    gOff.drawString("V(r1,r2,R1,R2) =",gb,140);
    gOff.drawString("     sum(-1/|ri-Rj| ,{i,j=1,2})",gb,155);
    gOff.drawString("     +1/|r1-r2| ",gb,170);

    gOff.drawString("N="+dmc.getNumberOfReplica()+"",gb,240);
    gOff.drawString("ER="+(float)(dmc.getReferenceEnergy())+"",gb,260);
    gOff.drawString("<ER>="+(float)(dmc.getAveragedReferenceEnergy())+"",gb,280); 
    gOff.drawString("Etotal="+(float)(dmc.getTotalEnergy())+"",gb,300); 	
    gOff.drawString("(Etotal=-1.165au at dPP=1.39au)",gb,320);

    gOff.drawString("thread count",gb,360);
    gOff.drawString("DMC ="+dmc.getCount()+"",gb,380);
    gOff.drawString("disp ="+thCount+"",gb,400);
  }


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

  void densityPlot(int mode) {
    int i,j,k,ix,iy,ir,fp;
    int iBegin, iEnd, iInc, jBegin, jEnd, jInc, kBegin, kEnd, kInc;
    double x,y,z,ppx,ppy,ppz,dens,psi;

    fp = drawBox(0);

    iBegin = 0; iEnd = 40; iInc = 1;
    if (wkx[fp]>10.0) {
      iBegin = 39; iEnd = 0; iInc = -1;
    }
    jBegin = 0; jEnd = 40; jInc = 1; if (mode==1) jEnd = 20;
    if (wky[fp]>10.0) {
      jBegin = 39; jEnd = 0; jInc = -1; if (mode==1) jBegin = 19;
    }
    kBegin = 0; kEnd = 40; kInc = 1;
    if (wkz[fp]>10.0) {
      kBegin = 39; kEnd = 0; kInc = -1;
    }

    for (i=iBegin; ((i<iEnd && iInc>0) || (i>=iEnd && iInc<0)) ; i += iInc) {
      x = 0.5*i;
      for (j=jBegin; ((j<jEnd && jInc>0) || (j>=jEnd && jInc<0)) ; j += jInc) {
        y = 0.5*j;
        for (k=kBegin; ((k<kEnd && kInc>0) || (k>=kEnd && kInc<0)) ; k += kInc) {
          z = 0.5*k;
          psi = dmc.getPsiDensity(i,j,k);
          dens = 0.000001*psi*psi;
          if (dens>0.999) dens = 0.999;
          ir = (int)(10.0*dens+0.5);
          if (ir>0) {
            ppy = cosFi*(y-cyy)+sinFi*(z-czz) + cyy; 
            ppz = -sinFi*(y-cyy)+cosFi*(z-czz) + czz;
            ppx = cosTh*(x-cxx)+sinTh*(ppz-czz) + cxx; 
            ppz = -sinTh*(x-cxx)+cosTh*(ppz-czz) + czz;

            ix = (int)(scale*ppx)+xImageLoc;
            iy = (int)(scale*ppy)+yImageLoc;
			
            gOff.setColor(Color.getHSBColor((float)(0.66-0.66*dens), 0.9f,0.9f));
            gOff.fillOval(ix-ir/2,iy-ir/2,ir,ir);
          }
        }
      }
    }

    drawBox(1);
  }

  //

  public void replicaPlot(double size) {
    int nn,i,j,k,jj,gx,gy,gz,g2x,g2y;
    double sz,rt,rg;

    setData();
    rotate();
    zSort();

    drawBox(0);

    for (i=0; i<nMax; i++) {
      j = srtz[0][i];
      gx = (int)(scale*px[j])+xImageLoc;
      gy = (int)(scale*py[j])+yImageLoc;
      gz = (int)(size*(0.4*pz[j]/zMax+0.6));

      gOff.setColor(Color.getHSBColor((float)(particleColor[pID[j]]), 0.95f,(float)(0.4*pz[j]/zMax+0.3)));
      gOff.drawOval(gx-gz/2,gy-gz/2, gz, gz);
    }

    drawBox(1);
  }

  // 

  private void setData() {
    int i,nn,np,idisp,ip;
    double cx,cy,cz;

    cx = xMax/2.0; cy = yMax/2.0; cz = zMax/2.0;

    nn = dmc.getMaxOfReplica(); np = dmc.getNumberOfParticle();
    idisp = 0;
    for (i=0; i<nn; i++) {
      if (dmc.getRepStatus(i)>0) { // if replica is alive
        for (ip=0; ip<np; ip++) {
          px[idisp] = cx + dmc.getReplica(i,ip,0);
          py[idisp] = cy + dmc.getReplica(i,ip,1);
          pz[idisp] = cz + dmc.getReplica(i,ip,2);
          pID[idisp] = (short)ip;
          idisp += 1;
        }
      }
    }
    nMax = idisp;
  }

  //

  private void rotate() {
    int i;
    double x,y,z;

    for (i=0; i<nMax; i++) {
	  x = px[i]; y = py[i]; z = pz[i];
      py[i] = cosFi*(y-cyy)+sinFi*(z-czz) + cyy; 
      pz[i] = -sinFi*(y-cyy)+cosFi*(z-czz) + czz;
      px[i] = cosTh*(x-cxx)+sinTh*(pz[i]-czz) + cxx; 
      pz[i] = -sinTh*(x-cxx)+cosTh*(pz[i]-czz) + czz;
    } 
  }

  //

  private void zSort() {
    int nn,i,j,w;

    nn = nMax;
    for (i=0; i<nn; i++) {
      srtz[0][i] = i;
      srtz[1][i] = (int)(32767.0*pz[i]/zMax);
    }
    qSort(0,nn-1);    
  }

  private void qSort(int i,int j) {
    int k,p;

    if(i==j) return;
    p = i+1;
    while (p<=j && srtz[1][i]==srtz[1][p]) p++;
    if (p<=j) {
      if (srtz[1][i]>=srtz[1][p]) p = i;
      k = qSortPartition(i,j,srtz[1][p]);
      qSort(i,k-1);
      qSort(k,j);
    }
  }

  private int qSortPartition(int i,int j, int x) {
    int l,r,w;

    l = i; r = j;
    while (l<=r) {
      while (l<=j && srtz[1][l]<x) l++;
      while (r>=i && srtz[1][r]>=x) r--;
      if (l>r) break;
      w = srtz[0][l]; srtz[0][l] = srtz[0][r]; srtz[0][r] = w;
      w = srtz[1][l]; srtz[1][l] = srtz[1][r]; srtz[1][r] = w;
      l++; r--;
    }
    return ( l );
  }

  //

  private int drawBox(int imode) {
    int i,gx,gy,g2x,g2y,farPoint;
    double sc;

    if (imode==0) {
      setBox();
      rotateBox();
    }
    sc = dispScale*viewScale;
    farPoint = findFarPoint();

    for (i=0; i<12; i++) {
      gx = (int)(sc*pwkx[boxp[i][0]])+xImageLoc;
      gy = (int)(sc*pwky[boxp[i][0]])+yImageLoc;
      g2x = (int)(sc*pwkx[boxp[i][1]])+xImageLoc;
      g2y = (int)(sc*pwky[boxp[i][1]])+yImageLoc;
      if (imode==0) {
        if (boxp[i][0]==farPoint || boxp[i][1]==farPoint) { 
          gOff.setColor(Color.lightGray);
          gOff.drawLine(gx,gy, g2x, g2y);
        }
      } else if (imode==1) {
        if (boxp[i][0]!=farPoint && boxp[i][1]!=farPoint) { 
          gOff.setColor(Color.darkGray);
          gOff.drawLine(gx,gy, g2x, g2y);
        }
      }
    }
    return farPoint;
  }

  private int findFarPoint() {
    int i,im;
    double m;

    im = 0; m = pwkz[im];
    for (i=0; i<8; i++) {
      if (pwkz[i]<m) {
        m = pwkz[i];
        im = i;
      }
    }
    return(im);
  }

  private void setBox() {
    int i,ix,iy,iz;
    double ballSize;

    i=0;
    ballSize = 0.03*xMax;

    for (ix=0; ix<2; ix++) {
      for (iy=0; iy<2; iy++) {
        for (iz=0; iz<2; iz++) {
          wkx[i] = -0.5*ballSize+(xMax+ballSize)*(double)(ix);
          wky[i] = -0.5*ballSize+(yMax+ballSize)*(double)(iy);
          wkz[i] = -0.5*ballSize+(zMax+ballSize)*(double)(iz);
          i += 1;
        }
      }
    }
  }

  private void rotateBox() {
    int i;
 
    for (i=0; i<8; i++) {
      pwky[i] = cosFi*(wky[i]-cyy)+sinFi*(wkz[i]-czz) + cyy; 
      pwkz[i] = -sinFi*(wky[i]-cyy)+cosFi*(wkz[i]-czz) + czz;
      pwkx[i] = cosTh*(wkx[i]-cxx)+sinTh*(pwkz[i]-czz) + cxx; 
      pwkz[i] = -sinTh*(wkx[i]-cxx)+cosTh*(pwkz[i]-czz) + czz;
    }  
  }


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

  public void setViewDelta(double xDisplacement, double yDisplacement) {
    viewTheta += 0.5*3.14159/180.0*xDisplacement;
    viewFai += 0.5*3.14159/180.0*yDisplacement;
    cosTh = Math.cos(viewTheta); 
    sinTh = Math.sin(viewTheta);
    cosFi = Math.cos(viewFai);
    sinFi = Math.sin(viewFai);
  }

  public void setScale(double sc) {
    viewScale = sc;
	scale = dispScale*viewScale;
  }
}



// ====================== diffusion Monte-Carlo method ===========

class DiffusionMC extends Thread {

  private Random rand = new Random();
  private double t = 0.0;
  private double dt = 0.1;
  private double deltaTau = 0.1;
  private double sqrtdt = Math.sqrt(deltaTau);
  private double hbar = 1.0;

  private int NNrep = 2000; // max number of replicant
  private int NNp = 2; // number of particle (in 3D)
  private int NNdim = 3; // 3D
  private int repStatus[] = new int[NNrep];
  private double replica[][][] = new double[NNrep][NNp][NNdim];

  private int Ndiv = 40;
  private float psiDensity[][][] = new float[Ndiv][Ndiv][Ndiv];

  private double referenceEnergy = 0.0;
  private double averagedReferenceEnergy = referenceEnergy;

  private int NN0 = 1000;
  private int NN1 = NN0;
  private int NN00 = NN0;

  private int sleepTime = 5;
  private int stopSleepTime = 100;

  private int resetSW = 0;
  private int startSW = 1;
  private int stepSW = 0;
  private int densSW = 1;

  private double dPP = 1.39;

  private int averageCount = 100;
  private double averageFactor = 1.0/((double)averageCount);


  private int dmcCount = 0;


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

  DiffusionMC() {
    setInitialCondition();
  }


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

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


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

  private void setInitialCondition() {
    int i,jp,jd,n0;

    t = 0.0;
    dmcCount = 0;

    for (i=0; i<NNrep; i++) {
      repStatus[i] = 0;
    }

    n0 = NNrep/2;
    for (i=0; i<n0; i++) {
      repStatus[i] = 1;
      for (jp=0; jp<NNp; jp++) {
        for (jd=0; jd<NNdim; jd++) {
		  replica[i][jp][jd] = 4.0*(Math.random()-0.5);
          //replica[i][jp][jd] = particle[jp][3+jd]+1.0/particle[jp][2]*(Math.random()-0.5);
        }
      }
    }

    NN0 = numberOfReplica();
    NN00 = NN0;
    referenceEnergy = averagedPotential();

    clearPsiDensity();
  }

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

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


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

  private void timeEvolution() {
    int i;

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

  private void timeStep() {
    int i,jp,jd;
    double q;

    t = t + dt;

    // (1) replica walk
    for (i=0; i<NNrep; i++) {
      if (repStatus[i]==1) {
        for (jp=0; jp<NNp; jp++) {
          q = 1.0; //q = 1.0/particle[jp][2];
          for (jd=0; jd<NNdim; jd++) {
            replica[i][jp][jd] += q*sqrtdt*gaussianRandom();
          }
        }
      }
    }

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

    // (3) count
    NN1 = numberOfReplica();

    // (4) energy
    referenceEnergy += hbar/deltaTau*(1.0-(NN1+1.0)/(NN0+1.0)*(NN0+10.0*NNrep)/(NN00+10.0*NNrep));
    NN0 = NN1;

    averagedReferenceEnergy = (1.0-averageFactor)*averagedReferenceEnergy+averageFactor*referenceEnergy;

    if (densSW==1) {
      setPsiDensity();
    }
  }

  private void birthDeathProcess(int i) {
    int mn;

    mn = (int)(weightFunction(i)+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);
    }
  }

  private void birthReplica(int i) {
    int ip,jp,jd;

    ip = seekSpace(i);
    if (ip>=0) {
      repStatus[ip] = 2;
      for (jp=0; jp<NNp; jp++) {
        for (jd=0; jd<NNdim; jd++) {
          replica[ip][jp][jd] = replica[i][jp][jd];
        }
      }
    }
  }

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

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


  private void setPsiDensity() {
    int i,ip,ix,iy,iz;
	
    for (i=0; i<NNrep; i++) {
      if (repStatus[i]==1) {
        for (ip=0; ip<NNp; ip++) {
          ix = (int)(Ndiv*(replica[i][ip][0]+10.0)/20.0);
          iy = (int)(Ndiv*(replica[i][ip][1]+10.0)/20.0);
          iz = (int)(Ndiv*(replica[i][ip][2]+10.0)/20.0);
          if (ix>=0 && ix<Ndiv && iy>=0 && iy<Ndiv && iz>=0 && iz<Ndiv ) {
            psiDensity[ix][iy][iz] += 1.0f;
          }
        }
      }
    }
    decayPsiDensity();
  }

  private void decayPsiDensity() {
    int ix,iy,iz;
    float dc;

    dc = 0.999f;
    for (ix=0; ix<Ndiv; ix++) {
      for (iy=0; iy<Ndiv; iy++) {
        for (iz=0; iz<Ndiv; iz++) {
          psiDensity[ix][iy][iz] *= dc;
        }
      }
    }
  }

  private void clearPsiDensity() {
    int ix,iy,iz;
	
    for (ix=0; ix<Ndiv; ix++) {
      for (iy=0; iy<Ndiv; iy++) {
        for (iz=0; iz<Ndiv; iz++) {
          psiDensity[ix][iy][iz] = 0.0f;
        }
      }
    }
  }


  // ----------------------------------------- utility -----------

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

  private double weightFunction(int i) {
    return(Math.exp(-(potential(i)-referenceEnergy)*deltaTau/hbar));
  }

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

  /*
  private double gaussianRandom() { // Box-Muller method
    double rnd;
    rnd = Math.sqrt( -2.0*Math.log(Math.random()) );
    return rnd*Math.sin(2.0*Math.PI*Math.random());
  }
  */

  private double potential(int i) {
    double x1,y1,z1,x2,y2,z2;
    double rr,r12,r22,r32,r42,r52;

    x1 = replica[i][0][0]; y1 = replica[i][0][1]; z1 = replica[i][0][2];
    x2 = replica[i][1][0]; y2 = replica[i][1][1]; z2 = replica[i][1][2];
    rr = dPP;
    r12 = x1*x1+y1*y1+(z1-0.5*rr)*(z1-0.5*rr);
    if (r12<1.0e-10) r12=1.0e-10;
    r22 = x1*x1+y1*y1+(z1+0.5*rr)*(z1+0.5*rr);
    if (r22<1.0e-10) r22=1.0e-10;
    r32 = x2*x2+y2*y2+(z2-0.5*rr)*(z2-0.5*rr);
    if (r32<1.0e-10) r32=1.0e-10;
    r42 = x2*x2+y2*y2+(z2+0.5*rr)*(z2+0.5*rr);
    if (r42<1.0e-10) r42=1.0e-10;
    r52 = (x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2);
    if (r52<1.0e-10) r52=1.0e-10;
    return(-1.0/Math.sqrt(r12)-1.0/Math.sqrt(r22)-1.0/Math.sqrt(r32)-1.0/Math.sqrt(r42)+1.0/Math.sqrt(r52));
  }


  // ------------------------------------- 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 void setDensSW(int sw) { densSW = sw; }
  public void setppDistance(double d) { dPP = d; }
  public double getppDistance() { return dPP; }
  public void setAverageCount(int ac) {
    averageCount = 100;
    if (ac==1) averageCount = 1000;
    if (ac==2) averageCount = 10000;
    averageFactor = 1.0/((double)averageCount);
  }

  // get data
  public int getMaxOfReplica() { return NNrep; }
  public int getNumberOfReplica() { return NN0; }
  public int getNumberOfParticle() { return NNp; } 

  public double getTime() { return t; }
  public double getReplica(int i,int ip,int id) { return replica[i][ip][id]; }
  public int getRepStatus(int i) { return repStatus[i]; }
  public double getPsiDensity(int i, int j, int k) { return psiDensity[i][j][k]; }

  public double getReferenceEnergy() { return referenceEnergy; }
  public double getAveragedReferenceEnergy() { return averagedReferenceEnergy; }
  public double getTotalEnergy() { return (averagedReferenceEnergy+1.0/dPP); }

  public int getCount() { return dmcCount; }

}

/* diffusion Monte Carlo method
 *
 * I. Kosztin et. al.; e-print arXiv:physics/9702023v1 (1997)
 * "Introduction to the Diffusion Monte Carlo Method"
 *
 * 
 * DMC algolithm
 *
 * (1) replica walk
 *     x -> x + sqrt(dtau)*GaussianRandom
 *
 * (2) branch: replica birth-death (branching) processes
 *     weight function W
 *       W = exp( -(V-ER)*dtau/h ) --> Wi ~ 1-(Vi-Ei)*dtau/h
 *
 *     m = int( wi + random[0-1] )
 *     m == 0 death
 *     m == 1 alive
 *     m == 2 alive + birth one replica
 *     m >3   alive + birth two replica
 *
 * (3) count: number of replicas N1
 *
 * (4) energy calc
 *     ER(k+1) = ER(k) + (h/dtau)(1-N1/N0)
 *        (N0 reference number)
 *
 * goto (1)
 *
 */