/* H2+ ion in The Box QMD3D     */
/* coded by Ikeuchi Mitsuru     */
/* ver 0.0.0   2004.05.05       */
/* ver 0.0.0.2   2004.05.23  improved reset code */

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

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

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

  Choice cvw;
  Button bt_reset, bt_start, bt_stop;
  Scrollbar sxx,svx;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 50;

  int dispMode = 5;
  int started = 1;
  double xx0Nuc = 3.0;
  double vx0Nuc = 0.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 dx = 8.0/32.0;
  double dy = 8.0/32.0;
  double dz = 8.0/32.0;
  double dt = 1.0*(dx*dx);

  double xField = 0.0;
  double yField = 0.0;
  double zField = 0.0;

  double xxNuc = 0.0;
  double yyNuc = 0.0;
  double zzNuc = 0.0;
  double vxNuc = 0.0;
  double vyNuc = 0.0;
  double vzNuc = 0.0;
  double ffxNuc = 0.0;
  double ffyNuc = 0.0;
  double ffzNuc = 0.0;
  double fbxNuc = 0.0;
  double fbyNuc = 0.0;
  double fbzNuc = 0.0;

  double x1Field = 0.0;
  double y1Field = 0.0;
  double z1Field = 0.0;

  double xx2Nuc = 0.0;
  double yy2Nuc = 0.0;
  double zz2Nuc = 0.0;
  double vx2Nuc = 0.0;
  double vy2Nuc = 0.0;
  double vz2Nuc = 0.0;
  double ffx2Nuc = 0.0;
  double ffy2Nuc = 0.0;
  double ffz2Nuc = 0.0;
  double fbx2Nuc = 0.0;
  double fby2Nuc = 0.0;
  double fbz2Nuc = 0.0;

  double x2Field = 0.0;
  double y2Field = 0.0;
  double z2Field = 0.0;

  double xProtonField = 0.0;
  double yProtonField = 0.0;
  double zProtonField = 0.0;

  int NNx = 32+1;
  int NNy = 32+1;
  int NNz = 32+1;
  int NN = NNx+1; /* max(NNx,NNy,NNz)+1 */

  double phRe[][][] = new double[NNx][NNy][NNz];
  double phIm[][][] = new double[NNx][NNy][NNz];  

  double vv[][][] = new double[NNx][NNy][NNz];

  double aaRe[] = new double[NN];
  double aaIm[] = new double[NN];
  double bbRe[] = new double[NN];
  double bbIm[] = new double[NN];
  double bRe[] = new double[NN];
  double bIm[] = new double[NN];
  double uRe[] = new double[NN];
  double uIm[] = new double[NN];

  int xpts[] = new int[NN];
  int ypts[] = new int[NN];
  double icol[] = new double[NN];

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

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

  public void init() {

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

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

    cvw = new Choice();
    cvw.addItem("dens 3D"); cvw.addItem("phase 3D");
    cvw.addItem("grid ph2"); cvw.addItem("grid V");cvw.addItem("grid ph2+V");
    cvw.addItem("dens 2D"); cvw.addItem("phase 2D");
    cvw.addItemListener(this);
    cvw.select("dens 2D");

    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);
    
    sxx= new Scrollbar(Scrollbar.HORIZONTAL,30,10,20,55);
    sxx.addAdjustmentListener(this);

    svx= new Scrollbar(Scrollbar.HORIZONTAL,0,10,0,30);
    svx.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(cvw);
    pnl.add(sxx);
    pnl.add(svx);
    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() == cvw){
      dispMode = cvw.getSelectedIndex(); 
    }
  }

  public void actionPerformed(ActionEvent ev){
    if(ev.getSource() == bt_reset){ 
      t = 0.0;
      setInitialCondition();
      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() == sxx) { 
      xx0Nuc = 0.1*(double)(sxx.getValue());
    } else if (ev.getSource() == svx) { 
      vx0Nuc = 1000.0*(double)(svx.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 gb;
    double ev,ek;

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

	   ev = meanPotential();
    ek = meanKinetic();

    if (dispMode==0) {

      boxPlot();
      dens3DPlot(160.0,0);
      boxPlot2();

    } else if (dispMode==1) {

      boxPlot();
      dens3DPlot(160.0,1);
      boxPlot2();

    } else if (dispMode==2) {

      gridPlot(NNz/2,50.0,0.0);

    } else if (dispMode==3) {

      gridPlot(NNz/2,0.0,2.0);

    } else if (dispMode==4) {

      gridPlot(NNz/2,50.0,2.0);

    } else if (dispMode==5) {

      densPlot(20.0,0);

    } else if (dispMode==6) {

      densPlot(20.0,1);

    }

    gb = 400;
    gOff.setColor(Color.black);
    gOff.drawString("view",630/6*3+10,40);
    gOff.drawString("xx0="+(int)(xx0Nuc*10.0)/10.0+" au",630/6*4+10,40);
    gOff.drawString("vx0="+(int)(vx0Nuc/100.0)/10.0+" km/s",630/6*5+10,40);

    gOff.setColor(Color.blue);
    gOff.drawString("t="+(int)(t*2418.9+0.5)/100.0+" as",630/6*0+10,40);
    gOff.drawString("norm="+norm()+" ",630/6*1+10,40);
    gOff.drawString("xx ="+(int)(xxNuc*100.0)/100.0+" ",630/6*4+10, 60);
    gOff.drawString("vx ="+(int)(vxNuc/10)/100.0+" ",630/6*5+10,60);
    gOff.drawString("xx2 ="+(int)(xx2Nuc*100.0)/100.0+" ",630/6*4+10, 80);
    gOff.drawString("vx2 ="+(int)(vx2Nuc/10)/100.0+" ",630/6*5+10,80);

    gOff.setColor(Color.blue);
    gOff.drawString("<V>="+ev+" au",gb,100);
    gOff.drawString("<K>="+ek+" au",gb,120);
    gOff.drawString("<K>+<V>="+(ek+ev)+" au",gb,140);
    gOff.drawString("x1Field ="+x1Field+" au",gb,160);
    gOff.drawString("x2Field ="+x2Field+" au",gb,180);

    gOff.setColor(Color.black);
    gOff.drawString("box="+(int)(NNx*dx)+"x"+ +(int)(NNy*dy)+"x"+(int)(NNz*dz)+" au(1 au = 52.918 pm)",gb,200);
    gOff.drawString("V(r)= -1.0/r (r>=0.1)",gb,220);
    gOff.drawString("    = -10.0 (r<0.1)",gb,240);
    gOff.drawString("phi(r,t=0) = A exp(-r)  (1s state)",gb,260);

  }

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

  void densPlot(double weight, int mode) {
    int i,j,k;
    int gbx,gby,mag,ix,iy,ir;
    double ph2,col;

    gbx = 30;
    gby = 70;
    mag = 6;
    gOff.setColor(Color.black);
    gOff.fillRect(gbx-1,gby-1,mag*NNx+1,mag*NNy+1);

    k = NNz/2;
    for (i=1; i<NNx-1; i++) {
      for (j=1; j<NNy-1; j++) {
        if (vv[i][j][k]>0.0) {
          gOff.setColor(Color.getHSBColor(0.33f, 0.95f, 0.6f));
          gOff.fillRect(gbx+mag*i,gby+mag*j, mag, mag);
        }
        ph2 = weight*(phRe[i][j][k]*phRe[i][j][k]+phIm[i][j][k]*phIm[i][j][k]);
        if (ph2>0.1) {
          if (ph2>1.0) {
            ph2 = 1.0;
          }
          if (mode==0) {
            gOff.setColor(Color.getHSBColor(0.15f, 0.95f, (float)(ph2)));
          } else if (mode==1) {
            gOff.setColor(Color.getHSBColor((float)(0.5+0.5*phRe[i][j][k]/Math.sqrt(ph2/weight)), 0.9f, (float)(ph2)));
          }
          gOff.fillRect(gbx+mag*i,gby+mag*j, mag, mag);
        }
      }
    }

    ix = gbx+(int)(mag*xxNuc/dx)+mag/2;
    iy = gby+(int)(mag*yyNuc/dy)+mag/2;
    ir = mag*2;
    gOff.setColor(Color.getHSBColor(0.01f,0.9f,0.8f));
    gOff.fillOval(ix-ir/2,iy-ir/2,ir,ir);

    ix = gbx+(int)(mag*xx2Nuc/dx)+mag/2;
    iy = gby+(int)(mag*yy2Nuc/dy)+mag/2;
    ir = mag*2;
    gOff.setColor(Color.getHSBColor(0.85f,0.9f,0.8f));
    gOff.fillOval(ix-ir/2,iy-ir/2,ir,ir);

  }

  void gridPlot(int zPos,double phWeight, double vvWeight) {
    int i,j,k,gbx,gby;
    int cx,cy,cz;
    double z,px,py,pz,sc;
    double ph2;

    cx = NNx/2;
    cy = NNy/2;
    cz = 0;
    sc = 8.0;

    gbx = 30+(int)(cx/sc);
    gby = 50+(int)(cy/sc);
    k = zPos;

    for(j=0; j<NNy; j += 1) {
      for (i=0; i<NNx; i++) { 
        ph2 = phRe[i][j][k]*phRe[i][j][k]+phIm[i][j][k]*phIm[i][j][k];
        z = phWeight*ph2 + vvWeight*vv[i][j][k] + cz;
        py = cosFi*(j-cy)+sinFi*(z-cz) + cy; 
        pz = -sinFi*(j-cy)+cosFi*(z-cz) + cz;
        px = cosTh*(i-cx)+sinTh*(pz-cz) + cx; 
        pz = -sinTh*(i-cx)+cosTh*(pz-cz) + cz;
        xpts[i] = (int)(sc*px)+gbx;
        ypts[i] = (int)(sc*py)+gby;
        if (ph2>0.01) {
          icol[i] = 0.15;
        } else if (vv[i][j][k]>0.0) {
          icol[i] = 0.33;
        } else {
          icol[i] = 0.66;
        }
      }
      for (i=0; i<NNx-1; i++) {
        gOff.setColor(Color.getHSBColor((float)icol[i],0.9f,0.8f));
        gOff.drawLine(xpts[i],ypts[i],xpts[i+1],ypts[i+1]);
      }
    }

    for(i=0; i<NNx; i += 1) {
      for (j=0; j<NNy; j++) {
        ph2 = phRe[i][j][k]*phRe[i][j][k]+phIm[i][j][k]*phIm[i][j][k];
        z = phWeight*ph2 + vvWeight*vv[i][j][k] + cz;
        py = cosFi*(j-cy)+sinFi*(z-cz) + cy; 
        pz = -sinFi*(j-cy)+cosFi*(z-cz) + cz;
        px = cosTh*(i-cx)+sinTh*(pz-cz) + cx; 
        pz = -sinTh*(i-cx)+cosTh*(pz-cz) + cz;
        xpts[j] = (int)(sc*px)+gbx;
        ypts[j] = (int)(sc*py)+gby;
        if (ph2>0.01) {
          icol[j] = 0.15;
        } else if (vv[i][j][k]>0.0) {
          icol[j] = 0.33;
        } else {
          icol[j] = 0.66;
        }
      }
      for (j=0; j<NNy-1; j++) { 
        gOff.setColor(Color.getHSBColor((float)icol[j],0.9f,0.8f));
        gOff.drawLine(xpts[j],ypts[j],xpts[j+1],ypts[j+1]);
      }
    }

  }

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

    gbx = 60;
    gby = 80;
    mag = 6.0;
    cxx = mag*(NNx/2.0); 
    cyy = mag*(NNy/2.0);
    czz = mag*(NNz/2.0);

    pointPlot(xxNuc,yyNuc,zzNuc,0.01);
    pointPlot(xx2Nuc,yy2Nuc,zz2Nuc,0.85);

    for (i=0; i<NNx; i+=1) {
      for(j=0; j<NNy; j +=1) {
        for (k=0; k<NNz; k+=1) {
          if (vv[i][j][k]>0.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;
            ir = 2;
            ix = (int)(ppx)+gbx;
            iy = (int)(ppy)+gby;
            gOff.setColor(Color.getHSBColor(0.33f, 0.9f,0.9f));
            gOff.fillOval(ix-ir/2,iy-ir/2,ir,ir);
          }
          dens =  phRe[i][j][k]*phRe[i][j][k]+phIm[i][j][k]*phIm[i][j][k];
          ir = (int)(weight*dens);
          if (ir>0) {
            if (ir>10) {
              ir = 10;
            }
            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;
            col = 0.0;
            if (mode==0) {
              col = 0.66 -0.05*dens;
              if (col<0.0) {
                col = 0.0;
              }
            } else if (mode==1) {
              col = 0.5+0.5*phRe[i][j][k]/Math.sqrt(dens);
            }
            gOff.setColor(Color.getHSBColor((float)col, 0.9f,0.9f));
            gOff.fillOval(ix-ir/2,iy-ir/2,ir,ir);
          }
        }
      }
    }
  }

  void pointPlot(double px,double py,double pz,double col) {
    int i,j,k,ix,iy,ir;
    int gbx,gby;
    double cxx,cyy,czz;
    double mag,z,ppx,ppy,ppz;

    gbx = 60;
    gby = 80;
    mag = 6.0;
    cxx = mag*(NNx/2.0); 
    cyy = mag*(NNy/2.0);
    czz = mag*(NNz/2.0);

    i = (int)(px/dx); j = (int)(py/dy); k = (int)(pz/dz);
    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;
    ir = 6;
    ix = (int)(ppx)+gbx;
    iy = (int)(ppy)+gby;
    gOff.setColor(Color.getHSBColor((float)col, 0.9f,0.9f));
    gOff.fillOval(ix-ir/2,iy-ir/2,ir,ir);

  }


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

    gbx = 60;
    gby = 70;
    mag = 6.0;
    cxx = mag*(NNx/2.0); 
    cyy = mag*(NNy/2.0);
    czz = mag*(NNz/2.0);

    ip = 0;
    for (i=0; i<NNx; i+=NNx-1) {
      for(j=0; j<NNy; j+=NNy-1) {
        for (k=0; k<NNz; k+=NNz-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]);
      }
    }
  }

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

  double norm() {
    int i,j,k;
    double p;

    p = 0.0;
    for (i=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        for (k=1;k<NNz-1;k++) {
          p += (phRe[i][j][k]*phRe[i][j][k]+phIm[i][j][k]*phIm[i][j][k])*dx*dy*dz;
        }
      }
    }
    return ( p );
  }

  double meanPotential() {
    int i,j,k;
    double p;

    p = 0.0;
    for (i=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        for (k=1;k<NNz-1;k++) {
          p += vv[i][j][k]*(phRe[i][j][k]*phRe[i][j][k]+phIm[i][j][k]*phIm[i][j][k])*dx*dy*dz;
        }
      }
    }
    return ( p );
  }

  double meanKinetic() {
    int i,j,k;
    double p,d2phRe,d2phIm,h2;

    h2 = dx*dx;
    p = 0.0;
    for (i=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        for (k=1;k<NNz-1;k++) {
          d2phRe = (phRe[i+1][j][k]+phRe[i-1][j][k]+phRe[i][j+1][k]+phRe[i][j-1][k]+phRe[i][j][k-1]+phRe[i][j][k+1]-6.0*phRe[i][j][k])/h2;
          d2phIm = (phIm[i+1][j][k]+phIm[i-1][j][k]+phIm[i][j+1][k]+phIm[i][j-1][k]+phIm[i][j][k-1]+phIm[i][j][k+1]-6.0*phIm[i][j][k])/h2;
          p += (phRe[i][j][k]*d2phRe+phIm[i][j][k]*d2phIm)*dx*dy*dz;
        }
      }
    }
    return ( -0.5*p );
  }


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

  void setInitialCondition() {
    t = 0.0;
    xxNuc = xx0Nuc;
    yyNuc = 4.0;
    zzNuc = 4.0;
    vxNuc = vx0Nuc;
    vyNuc = 0.0;
    vzNuc = 0.0;
    ffxNuc = 0.0;
    ffyNuc = 0.0;
    ffzNuc = 0.0;
    fbxNuc = 0.0;
    fbyNuc = 0.0;
    fbzNuc = 0.0;

    xx2Nuc = 5.0;
    yy2Nuc = 4.0;
    zz2Nuc = 4.0;
    vx2Nuc = 0.0;
    vy2Nuc = 0.0;
    vz2Nuc = 0.0;
    ffx2Nuc = 0.0;
    ffy2Nuc = 0.0;
    ffz2Nuc = 0.0;
    fbx2Nuc = 0.0;
    fby2Nuc = 0.0;
    fbz2Nuc = 0.0;

    setWave(xxNuc,yyNuc,zzNuc, 0.0, 0.0, 0.0);
    setPotential(xxNuc,yyNuc,zzNuc,xx2Nuc,yy2Nuc,zz2Nuc);

  }

/* ----------------------------- set potential -------------- */

  void setPotential(double x0,double y0,double z0,double x20,double y20,double z20) {
    int i,j,k;
    double x,y,z,r,r2;

    for (i=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        for (k=1;k<NNz-1;k++) {
          x = i*dx;
          y = j*dy;
          z = k*dz;
          r = Math.sqrt((x-x0)*(x-x0)+(y-y0)*(y-y0)+(z-z0)*(z-z0));
          if (r<0.1) {
            r = 0.1;
          }
          r2 = Math.sqrt((x-x20)*(x-x20)+(y-y20)*(y-y20)+(z-z20)*(z-z20));
          if (r2<0.1) {
            r2 = 0.1;
          }
          vv[i][j][k] = -1.0/r -1.0/r2;

        }
      }
    }
  }

/* ----------------------------- set wave packet -------------- */

  void setWave(double xPos, double yPos, double zPos, double kx, double ky, double kz) {
    int i,j,k;
    double a,x,y,z,r,phAb,phPh;

    for (i=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        for (k=1;k<NNz-1;k++) {
          x = i*dx;
          y = j*dy;
          z = k*dz;
          r = Math.sqrt((x-xPos)*(x-xPos)+(y-yPos)*(y-yPos)+(z-zPos)*(z-zPos));
          phAb = Math.exp( -r );
          phPh = kx*x+ky*y+kz*z;
          phRe[i][j][k] = phAb*Math.cos(phPh);
          phIm[i][j][k] = phAb*Math.sin(phPh);
        }
      }
    }

    for (i=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        phRe[i][j][0] = 0.0;
        phIm[i][j][0] = 0.0;
        phRe[i][j][NNz-1] = 0.0;
        phIm[i][j][NNz-1] = 0.0;
      }
    }
    for (i=1;i<NNx-1;i++) {
      for (k=1;k<NNz-1;k++) {
        phRe[i][0][k] = 0.0;
        phIm[i][0][k] = 0.0;
        phRe[i][NNy-1][k] = 0.0;
        phIm[i][NNy-1][k] = 0.0;
      }
    }
    for (j=1;j<NNy-1;j++) {
      for (k=1;k<NNz-1;k++) {
        phRe[0][j][k] = 0.0;
        phIm[0][j][k] = 0.0;
        phRe[NNx-1][j][k] = 0.0;
        phIm[NNx-1][j][k] = 0.0;
      }
    }

    a = Math.sqrt(norm());
	
    for (i=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        for (k=1;k<NNz-1;k++) {
          phRe[i][j][k] = phRe[i][j][k]/a;
          phIm[i][j][k] = phIm[i][j][k]/a;
        }
      }
    }

  }

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

  void timeEvolution() {
    int i;

    if (started==-1) {
      t = 0.0;
      setInitialCondition();
      started = 0;
    } else if (started==1) {
      for (i=0; i<10; i++) {
        t = t + dt;
        phaseStep();
        kxStep();
        kyStep();
        kzStep();
        phaseStep();
      }
      moveAtom(10.0*dt);
    }

  }

  void moveAtom(double dtau) {
    double mNuc = 1.67e-27;
    double dtSI;

    dtSI = dtau*2.4189e-17;
    xxNuc += (vxNuc*dtSI+0.5*dtSI*dtSI*fbxNuc/mNuc)/5.2918e-11;
    yyNuc += (vyNuc*dtSI+0.5*dtSI*dtSI*fbyNuc/mNuc)/5.2918e-11;
    zzNuc += (vzNuc*dtSI+0.5*dtSI*dtSI*fbzNuc/mNuc)/5.2918e-11;
    xx2Nuc += (vx2Nuc*dtSI+0.5*dtSI*dtSI*fbx2Nuc/mNuc)/5.2918e-11;
    yy2Nuc += (vy2Nuc*dtSI+0.5*dtSI*dtSI*fby2Nuc/mNuc)/5.2918e-11;
    zz2Nuc += (vz2Nuc*dtSI+0.5*dtSI*dtSI*fbz2Nuc/mNuc)/5.2918e-11;
    fbxNuc = ffxNuc;
    fbyNuc = ffyNuc;
    fbzNuc = ffzNuc;
    fbx2Nuc = ffx2Nuc;
    fby2Nuc = ffy2Nuc;
    fbz2Nuc = ffz2Nuc;

    setFieldAt(xxNuc,yyNuc,zzNuc);
    x1Field=xField; y1Field=yField; z1Field=zField;
    setProtonField();
    ffxNuc = (4.3597e-18/5.2918e-11)*(x1Field+0.6*xProtonField);
    ffyNuc = (4.3597e-18/5.2918e-11)*(y1Field+0.6*yProtonField);
    ffzNuc = (4.3597e-18/5.2918e-11)*(z1Field+0.6*zProtonField);
    setFieldAt(xx2Nuc,yy2Nuc,zz2Nuc);
    x2Field=xField; y2Field=yField; z2Field=zField;
    ffx2Nuc = (4.3597e-18/5.2918e-11)*(x2Field-0.6*xProtonField);
    ffy2Nuc = (4.3597e-18/5.2918e-11)*(y2Field-0.6*yProtonField);
    ffz2Nuc = (4.3597e-18/5.2918e-11)*(z2Field-0.6*zProtonField);

    vxNuc += 0.5*dtSI*(fbxNuc+ffxNuc)/mNuc;
    vyNuc += 0.5*dtSI*(fbyNuc+ffyNuc)/mNuc;
    vzNuc += 0.5*dtSI*(fbzNuc+ffzNuc)/mNuc;
    vx2Nuc += 0.5*dtSI*(fbx2Nuc+ffx2Nuc)/mNuc;
    vy2Nuc += 0.5*dtSI*(fby2Nuc+ffy2Nuc)/mNuc;
    vz2Nuc += 0.5*dtSI*(fbz2Nuc+ffz2Nuc)/mNuc;
    setPotential(xxNuc,yyNuc,zzNuc,xx2Nuc,yy2Nuc,zz2Nuc);

  }

  void setFieldAt(double xPos,double yPos,double zPos) {
    int i,j,k;
    double sx,sy,sz,x,y,z,r,r2,p;

    sx = 0.0; sy = 0.0;sz = 0.0;
    for (i=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        for (k=1;k<NNz-1;k++) {
          x = i*dx;
          y = j*dy;
          z = k*dz;
          r2 = (x-xPos)*(x-xPos)+(y-yPos)*(y-yPos)+(z-zPos)*(z-zPos);
          r = Math.sqrt(r2);
          if (r>0.1) {
            p = 1.0/(r*r*r)*(phRe[i][j][k]*phRe[i][j][k] + phIm[i][j][k]*phIm[i][j][k])*dx*dy*dz;
            sx += (x-xPos)*p;
            sy += (y-yPos)*p;
            sz += (z-zPos)*p;
          }
        }
      }
    }
    xField = sx;
    yField = sy;
    zField = sz;

  }

  void setProtonField() {
    double x,y,z,r;

    x = xxNuc-xx2Nuc;
    y = yyNuc-yy2Nuc;
    z = zzNuc-zz2Nuc;
    r = Math.sqrt(x*x+y*y+z*z);

    xProtonField = x/(r*r*r);
    yProtonField = y/(r*r*r);
    zProtonField = z/(r*r*r);
  }

/* ----------------------------- time step -------------- */

  void kxStep() {
    int i,j,k;
    double a,aaAb,auAb;

    a = 4.0*dy*dy/dt;
    for (k=1; k<NNz-1; k++) {
      for (j=1; j<NNy-1; j++) {
        for (i=1; i<NNx-1; i++) {
          bbRe[i] = 2.0;
          bbIm[i] = a;
          aaRe[i] = -2.0;
          aaIm[i] = a;
          bRe[i] = bbRe[i]*phRe[i][j][k]-bbIm[i]*phIm[i][j][k] - phRe[i+1][j][k] - phRe[i-1][j][k];
          bIm[i] = bbRe[i]*phIm[i][j][k]+bbIm[i]*phRe[i][j][k] - phIm[i+1][j][k] - phIm[i-1][j][k];
        }

        aaAb = aaRe[1]*aaRe[1]+aaIm[1]*aaIm[1];
        uRe[1] = aaRe[1]/aaAb;
        uIm[1] = -aaIm[1]/aaAb;
        phRe[1][j][k] = bRe[1]*uRe[1] - bIm[1]*uIm[1];
        phIm[1][j][k] = bIm[1]*uRe[1] + bRe[1]*uIm[1];

        for (i=2; i<NNx-1; i++) {
          auAb = (aaRe[i]-uRe[i-1])*(aaRe[i]-uRe[i-1])+(aaIm[i]-uIm[i-1])*(aaIm[i]-uIm[i-1]);
          uRe[i] = (aaRe[i]-uRe[i-1])/auAb;
          uIm[i] = -(aaIm[i]-uIm[i-1])/auAb;
          phRe[i][j][k] = (bRe[i]-phRe[i-1][j][k])*uRe[i] - (bIm[i]-phIm[i-1][j][k])*uIm[i];
          phIm[i][j][k] = (bRe[i]-phRe[i-1][j][k])*uIm[i] + (bIm[i]-phIm[i-1][j][k])*uRe[i];
        }

        for (i=NNx-3; i>=1; i--) {
          phRe[i][j][k] -= phRe[i+1][j][k]*uRe[i] - phIm[i+1][j][k]*uIm[i];
          phIm[i][j][k] -= phRe[i+1][j][k]*uIm[i] + phIm[i+1][j][k]*uRe[i];
        }
      }
    }
  }

  void kyStep() {
    int i,j,k;
    double a,aaAb,auAb;

    a = 4.0*dy*dy/dt;
    for (i=1; i<NNx-1; i++) {
      for (k=1; k<NNz-1; k++) {
        for (j=1; j<NNy-1; j++) {
          bbRe[j] = 2.0;
          bbIm[j] = a;
          aaRe[j] = -2.0;
          aaIm[j] = a;
          bRe[j] = bbRe[j]*phRe[i][j][k]-bbIm[j]*phIm[i][j][k] - phRe[i][j+1][k] - phRe[i][j-1][k];
          bIm[j] = bbRe[j]*phIm[i][j][k]+bbIm[j]*phRe[i][j][k] - phIm[i][j+1][k] - phIm[i][j-1][k];
        }

        aaAb = aaRe[1]*aaRe[1]+aaIm[1]*aaIm[1];
        uRe[1] = aaRe[1]/aaAb;
        uIm[1] = -aaIm[1]/aaAb;
        phRe[i][1][k] = bRe[1]*uRe[1] - bIm[1]*uIm[1];
        phIm[i][1][k] = bIm[1]*uRe[1] + bRe[1]*uIm[1];

        for (j=2; j<NNy-1; j++) {
          auAb = (aaRe[j]-uRe[j-1])*(aaRe[j]-uRe[j-1])+(aaIm[j]-uIm[j-1])*(aaIm[j]-uIm[j-1]);
          uRe[j] = (aaRe[j]-uRe[j-1])/auAb;
          uIm[j] = -(aaIm[j]-uIm[j-1])/auAb;
          phRe[i][j][k] = (bRe[j]-phRe[i][j-1][k])*uRe[j] - (bIm[j]-phIm[i][j-1][k])*uIm[j];
          phIm[i][j][k] = (bRe[j]-phRe[i][j-1][k])*uIm[j] + (bIm[j]-phIm[i][j-1][k])*uRe[j];
        }

        for (j=NNy-3; j>=1; j--) {
          phRe[i][j][k] -= phRe[i][j+1][k]*uRe[j] - phIm[i][j+1][k]*uIm[j];
          phIm[i][j][k] -= phRe[i][j+1][k]*uIm[j] + phIm[i][j+1][k]*uRe[j];
        }
      }
    }
  }

  void kzStep() {
    int i,j,k;
    double a,aaAb,auAb;

    a = 4.0*dy*dy/dt;
    for (i=1; i<NNx-1; i++) {
      for (j=1; j<NNy-1; j++) {
        for (k=1; k<NNz-1; k++) {
          bbRe[k] = 2.0;
          bbIm[k] = a;
          aaRe[k] = -2.0;
          aaIm[k] = a;
          bRe[k] = bbRe[k]*phRe[i][j][k]-bbIm[k]*phIm[i][j][k] - phRe[i][j][k+1] - phRe[i][j][k-1];
          bIm[k] = bbRe[k]*phIm[i][j][k]+bbIm[k]*phRe[i][j][k] - phIm[i][j][k+1] - phIm[i][j][k-1];
        }

        aaAb = aaRe[1]*aaRe[1]+aaIm[1]*aaIm[1];
        uRe[1] = aaRe[1]/aaAb;
        uIm[1] = -aaIm[1]/aaAb;
        phRe[i][j][1] = bRe[1]*uRe[1] - bIm[1]*uIm[1];
        phIm[i][j][1] = bIm[1]*uRe[1] + bRe[1]*uIm[1];

        for (k=2; k<NNz-1; k++) {
          auAb = (aaRe[k]-uRe[k-1])*(aaRe[k]-uRe[k-1])+(aaIm[k]-uIm[k-1])*(aaIm[k]-uIm[k-1]);
          uRe[k] = (aaRe[k]-uRe[k-1])/auAb;
          uIm[k] = -(aaIm[k]-uIm[k-1])/auAb;
          phRe[i][j][k] = (bRe[k]-phRe[i][j][k-1])*uRe[k] - (bIm[k]-phIm[i][j][k-1])*uIm[k];
          phIm[i][j][k] = (bRe[k]-phRe[i][j][k-1])*uIm[k] + (bIm[k]-phIm[i][j][k-1])*uRe[k];
        }

        for (k=NNz-3; k>=1; k--) {
          phRe[i][j][k] -= phRe[i][j][k+1]*uRe[k] - phIm[i][j][k+1]*uIm[k];
          phIm[i][j][k] -= phRe[i][j][k+1]*uIm[k] + phIm[i][j][k+1]*uRe[k];
        }
      }
    }

  }

  void phaseStep() {
    int i,j,k;
    double th,cs,sn,phr,phi;

    for (i=1; i<NNx-1; i++) {
      for (j=1; j<NNy-1; j++) {
        for (k=1; k<NNz-1; k++) {
          th = 0.5*dt*vv[i][j][k];
          cs = Math.cos(th); sn = Math.sin(th);
          phr = phRe[i][j][k];
          phi = phIm[i][j][k];
          phRe[i][j][k] = cs*phr+sn*phi;
          phIm[i][j][k] = cs*phi-sn*phr;
        }
      }
    }
  }

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

}