/* hydrogen H2+ ion - diffusion Monte-Carlo method 3D */
/* coded by Ikeuchi Mitsuru   */
/* ver 0.0.1   2005.05.29     */

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

public class H2IonDMC3D extends Applet
  implements MouseListener, MouseMotionListener,
    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 dPP = 2.0;
  double boxSize = 20.0;


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

  int dgX, dgY, dgXb,dgYb;
  double viewTheta = -15.0*3.14/180.0;
  double viewFai = -72.0*3.14/180.0;
  double dtheta = 0.0*3.14/180.0;
  double pai = 3.1415926536;

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

  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 = 20;
  int NNby = 20;
  int NNbz = 20;
  int box[][][] = new int[NNbx][NNby][NNbz];
  double meanBox[][][] = new double[NNbx][NNby][NNbz];
  double psi[][][] = new double[NNbx][NNby][NNbz];

  int pxBox[] = new int[8];
  int pyBox[] = new int[8];
  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}
  };

  int xpts[] = new int[NNbx+1];
  int ypts[] = new int[NNbx+1];

/* ----------------------------- 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,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(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) { 
      dPP = 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;
    }
    viewTheta += 0.5*3.14/180.0*(dgX-dgXb);
    viewFai += 0.5*3.14/180.0*(dgY-dgYb);
    cosTh = Math.cos(viewTheta); 
    sinTh = Math.sin(viewTheta);
    cosFi = Math.cos(viewFai);
    sinFi = Math.sin(viewFai);
  }


  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;
    double mag;

    gx = 60; gy = 80;
    mag = 10.0;

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


    if (dispMode==0) {

      plot3DReplicas(gx,gy,mag);

	} else if (dispMode==1) {

	  setPsi();
	  plot3DMatrix(psi, 20.0,gx,gy,mag);

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


    gOff.setColor(Color.black);
    gOff.drawString("tau="+(int)(t*10+0.5)/10.0+" ",630/6*0+10,40);
    gOff.drawString("R="+(int)(dPP*100.0+0.5)/100.0+" au",630/6*3+10,40);
    gOff.drawString("view",630/6*5+10,40);

    gb = 400;
    gOff.drawString("L="+(int)(boxSize*10.0+0.5)/10.0+"",gb,80);
    gOff.drawString("Nmax="+NNr+"",gb,100);
    gOff.drawString("d tau="+(float)(deltaTau)+"",gb,120);
    gOff.drawString("V(r)=-1/|r-R/2|-1/|r+R/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("Etotal=ER+1/R="+(float)(eeRef+1.0/dPP)+"",gb,280);    
	gOff.drawString("(Etotal=-0.597au at R=2.0au)",gb,300);
  }

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

  void plot3DReplicas(int gx, int gy,double mag) {
    int i,ix,iy,ir;
    double cxx,cyy,czz,x,y,z,ppx,ppy,ppz;

    boxPlot(gx,gy,mag);

    cxx = mag*(NNbx/2.0); 
    cyy = mag*(NNby/2.0);
    czz = mag*(NNbz/2.0);

    gOff.setColor(Color.green);
    for (i=0; i<NNr; i++) {
      if (repStatus[i]==1) {
        x = cxx+(NNbx/boxSize)*mag*replica[i][0]; 
		y = cyy+(NNby/boxSize)*mag*replica[i][1];
		z = czz+(NNbz/boxSize)*mag*replica[i][2];
		
        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)(ppx)+gx;
		iy = (int)(ppy)+gy;
        ir = (int)(0.5+0.6*mag);
        gOff.setColor(Color.getHSBColor((float)(0.33), 0.9f,0.9f));
        gOff.drawOval(ix-ir/2,iy-ir/2,ir,ir);
	  }
    }
	
    boxPlot2();
  }


  void plot3DMatrix(double mat[][][], double weight, int gbx, int gby, double mag) {
    int i,j,k,ix,iy,ir;
    double cxx,cyy,czz;
    double z,ppx,ppy,ppz;
    double dens,col;

    boxPlot(gbx,gby,mag);

    cxx = mag*(NNbx/2.0); 
    cyy = mag*(NNby/2.0);
    czz = mag*(NNbz/2.0);

    for (i=1; i<NNbx-1; i+=2) {
      for(j=1; j<NNby-1; j +=2) {
        for (k=1; k<NNbz-1; k+=2) {

          dens = weight*mat[i][j][k];
		  if (dens>0.999) dens = 0.999;
          ir = (int)(1.2*mag*dens+0.5);
          if (ir>0) {
            z = k*mag;
            ppy = cosFi*(j*mag-cyy)+sinFi*(z-czz) + cyy; 
            ppz = -sinFi*(j*mag-cyy)+cosFi*(z-czz) + czz;
            ppx = cosTh*(i*mag-cxx)+sinTh*(ppz-czz) + cxx; 
            ppz = -sinTh*(i*mag-cxx)+cosTh*(ppz-czz) + czz;

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

    boxPlot2();
  }


  void boxPlot(int gbx, int gby, double mag) {
    int i,j,k,ip,i1,i2;
    double cxx,cyy,czz;
    double z,ppx,ppy,ppz;
    double dens;

    cxx = mag*(NNbx/2.0); 
    cyy = mag*(NNby/2.0);
    czz = mag*(NNbz/2.0);

    ip = 0;
    for (i=0; i<NNbx; i+=NNbx-1) {
      for(j=0; j<NNby; j+=NNby-1) {
        for (k=0; k<NNbz; k+=NNbz-1) {
          z = k*mag;
          ppy = cosFi*(j*mag-cyy)+sinFi*(z-czz) + cyy; 
          ppz = -sinFi*(j*mag-cyy)+cosFi*(z-czz) + czz;
          ppx = cosTh*(i*mag-cxx)+sinTh*(ppz-czz) + cxx; 
          ppz = -sinTh*(i*mag-cxx)+cosTh*(ppz-czz) + czz;

          pxBox[ip] = (int)(ppx)+gbx;
          pyBox[ip] = (int)(ppy)+gby;
          ip += 1;
        }
      }
    }
    gOff.setColor(Color.gray);
    for (i=0; i<12; i++) {
      i1= boxp[i][0]; i2= boxp[i][1];
      if (i1==0) {
        gOff.setColor(Color.gray);
        gOff.drawLine(pxBox[i1],pyBox[i1],pxBox[i2],pyBox[i2]);
      }
    }
  }

  void boxPlot2() {
    int i,i1,i2;

    for (i=0; i<12; i++) {
      i1= boxp[i][0]; i2= boxp[i][1];
      if (i1 != 0) {
        gOff.setColor(Color.black);
        gOff.drawLine(pxBox[i1],pyBox[i1],pxBox[i2],pyBox[i2]);
      }
    }
  }

/*-------------------------- 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);
      replica[i][2] = 4.0*(Math.random()-0.5);
	}
    setBox(0);

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

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

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

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

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

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

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

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

	for (i=0; i<NNbx; i++) {
	  for (j=0; j<NNby; j++) {
        for (k=0; k<NNbz; k++) {
	      psi[i][j][k] = meanBox[i][j][k]/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],replica[i][2]);
	  }
	}
	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();
        replica[i][2] += 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,zi;

    xi = replica[i][0]; yi = replica[i][1]; zi = replica[i][2];
    mn = (int)(weightFunction(xi,yi,zi)+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];
	  replica[ip][2] = replica[i][2];
    }
  }

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

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

  double potential(double x, double y, double z) {
    double rr,r12,r22;

    rr = dPP;
    r12 = x*x+y*y+(z-0.5*rr)*(z-0.5*rr);
	if (r12<1.0e-10) r12=1.0e-10;
    r22 = x*x+y*y+(z+0.5*rr)*(z+0.5*rr);
	if (r22<1.0e-10) r22=1.0e-10;
    return(-1.0/Math.sqrt(r12)-1.0/Math.sqrt(r22));
  }


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

}