/** applet No. 1120
 *
 * plane wave - quantum electron dynamics 3D
 * - finite difference explicit method
 * - multi thread  PQM3D :: display - asynchronous
 *
 * Created by Ikeuchi Mitsuru on February 10 2007.
 * Copyright (c) 2007 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2007.02.10  created
 * ver 0.0.2   2007.06.16  improved code
 *
 */

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

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

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

  Thread th = null; // for run()-paint() thread
  PQM3D qed = null;
  DisplayQED3D dsp = new DisplayQED3D();

  // for event
  Choice ch_view;
  Button bt_reset, bt_startStop, bt_step;
  Scrollbar sc_px0, sc_py0, sc_pz0, sc_kPos;

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

  int dgX, dgY, dgXb,dgYb;

  int sleepTime = 50;
  int dispMode = 0;
  int thCount = 0;


  // ------------------------------ applet main thread -----------

  public void init() {

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

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

    ch_view = new Choice();
    ch_view.add("psi^2");
    ch_view.add("phase");
    ch_view.add("p flow");
    ch_view.add("cloud");
    ch_view.add("grid 2D");
    ch_view.add("psi^2 2D");
    ch_view.add("phase 2D");
    ch_view.addItemListener(this);
    ch_view.select("psi^2");

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

    sc_px0 = new Scrollbar(Scrollbar.HORIZONTAL,40,10,0,50);
    sc_px0.addAdjustmentListener(this);

    sc_py0 = new Scrollbar(Scrollbar.HORIZONTAL,20,10,0,50);
    sc_py0.addAdjustmentListener(this);

    sc_pz0 = new Scrollbar(Scrollbar.HORIZONTAL,20,10,0,50);
    sc_pz0.addAdjustmentListener(this);

    sc_kPos = new Scrollbar(Scrollbar.HORIZONTAL,32,10,0,73);
    sc_kPos.addAdjustmentListener(this);

    addMouseListener(this);  
    addMouseMotionListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(2,6,5,0));
    //pnl.add(new Label("  "));
    pnl.add(bt_reset);
    pnl.add(bt_startStop);
    pnl.add(bt_step);
    pnl.add(new Label("  "));
    pnl.add(new Label("  "));
    pnl.add(new Label("  "));

    pnl.add(sc_px0);
    pnl.add(new Label("  "));
    pnl.add(new Label("  "));
    pnl.add(new Label("  "));
    pnl.add(sc_kPos);
    pnl.add(ch_view);
    add(pnl,"North");
  }

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

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


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

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

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

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_px0) { 
      qed.setpx0( 0.1*(double)(sc_px0.getValue())-2.0 );
    } else if (ev.getSource() == sc_py0) { 
      qed.setpy0( 0.1*(double)(sc_py0.getValue())-2.0 );
    } else if (ev.getSource() == sc_pz0) { 
      qed.setpz0( 0.1*(double)(sc_pz0.getValue())-2.0 );
    } else if (ev.getSource() == sc_kPos) { 
      dsp.setkPos( sc_kPos.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;
    }
    dsp.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 (qed != null) {
      dsp.plotGraph(gOff, qed, dispMode, thCount);
    }
  }

}



// =============================== display QED3D class ===========

class DisplayQED3D {

  private Graphics gOff;
  private PQM3D qed;

  private double dt = PQM3D.dt;
  private int NNx = PQM3D.NNx;
  private int NNy = PQM3D.NNy;
  private int NNz = PQM3D.NNz;
  private double dx = PQM3D.dx;
  private double dy = PQM3D.dy;
  private double dz = PQM3D.dz;
  private double dv = dx*dy*dz;
  private double xMax = NNx*dx;
  private double yMax = NNy*dy;
  private double zMax = NNz*dz;
  private double cxx = xMax/2.0;
  private double cyy = yMax/2.0;
  private double czz = zMax/2.0;

  private double dispScale = 200.0/yMax;
  private double viewScale = 1.0;
  private double scale = dispScale*viewScale;

  private int xImageLoc = 50;
  private int yImageLoc = 150;
  
  private int kPos = 32;

  private double viewTheta = -22.5*3.1416/180.0;
  private double viewFai = -40.0*3.1416/180.0;
  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 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 int xpts[] = new int[500];
  private int ypts[] = new int[500];
  private float colz[] = new float[500];

  private int NNr = 1000;
  private double srh[] = new double[NNr+1];
  private short cloud[][][] = new short[NNx][NNy][NNz];

  private int drawnCount = 0;


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

  DisplayQED3D() {
  }


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

  public void plotGraph(Graphics gg, PQM3D qe, int dispMode, int thCount) {
    int qedCount,gb,i;
    double qqv,qqk;
    gOff = gg;
    qed = qe;
    if (qed==null) return;
    qedCount = qed.getCount();
    //if (qedCount==drawnCount) return;
    //drawnCount = qedCount;

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

    if (dispMode==0) { // density3D psi^2
      density3DPlot(0,2);

    } else if (dispMode==1) { // density3D phase
      density3DPlot(1,2);

    } else if (dispMode==2) { // density3D flow
      density3DPlot(2,4);

    } else if (dispMode==3) { // density3D cloud
      density3DPlot(3,1);

    } else if (dispMode==4) { // grid psi^2
      gridPlot(0, 1.0);

    } else if (dispMode==5) { // density2D psi^2
      densityPlot(0);

    } else if (dispMode==6) { // density2D phase
      densityPlot(1);

    } else if (dispMode==7) {
      ;

    }

    gOff.setColor(Color.black);
    gOff.drawString("px0="+(float)(qed.getpx0())+" ",630/6*0+10,70);
    //gOff.drawString("py0="+(float)(qed.getpy0())+" ",630/6*1+10,70);
    //gOff.drawString("pz0="+(float)(qed.getpz0())+" ",630/6*2+10,70);

    gOff.setColor(Color.black);
    gOff.drawString("2D z="+(kPos-32)*dz+"",630/6*4+10,70);
    gOff.drawString("view",630/6*5+20,70);

    gOff.setColor(Color.blue);
    gOff.drawString("t="+(int)(qed.getTime()*2418.9+0.5)/100.0+" as",10,90);
    gOff.drawString("norm="+qed.getNorm()+" ",120,90);
    gOff.setColor(Color.black);
    gOff.drawString("box = 32x16x16 au(1au=0.0529nm)",360,90);

    qqv = qed.getMeanPotential();
    qqk = qed.getMeanKinetic();
    gOff.setColor(Color.blue);
    gOff.drawString("<V>= "+(float)(qqv)+" ",500,120);
    gOff.drawString("<K>= "+(float)(qqk)+" ",500,140);
    gOff.drawString("<E>= "+(float)(qqv+qqk)+" ",500,160);

    gOff.setColor(Color.black);
    gOff.drawString("thread count",500,200);
    gOff.setColor(Color.blue);
    gOff.drawString("PQM3D= "+qed.getCount()+" ",500,220);
    gOff.drawString("Disp= "+thCount+" ",500,240);

    //gOff.drawString("theta= "+(float)(viewTheta/3.1416*180.0)+" ",500,380);
    //gOff.drawString("fai= "+(float)(viewFai/3.1416*180.0)+" ",500,400);
  }


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

  private void density3DPlot(int mode, int inc) {
  // mode: 0-psi^2 1-phase  2-flow
    int i,j,k,ir,fp;
    int iBegin,iEnd,iInc, jBegin,jEnd,jInc, kBegin,kEnd,kInc;
    double x,y,z,dens;
    double v[];

    if (mode==3) { // cloud
      setCloud();
    }

    fp = drawBox(0);
    iBegin = 0; iEnd = NNx-inc; iInc = inc;
    if (wkx[fp]>xMax/2.0) {
      iBegin = NNx-inc; iEnd = 0; iInc = -inc;
    }
    jBegin = 0; jEnd = NNy-inc; jInc = inc;
    if (wky[fp]>yMax/2.0) {
      jBegin = NNy-inc; jEnd = 0; jInc = -inc;
    }
    kBegin = 0; kEnd = NNz-inc; kInc = inc;
    if (wkz[fp]>zMax/2.0) {
      kBegin = NNz-inc; kEnd = 0; kInc = -inc;
    }

    for (i=iBegin; ((i<=iEnd && iInc>0) || (i>=iEnd && iInc<0)) ; i += iInc) {
      x = dx*i;
      for (j=jBegin; ((j<=jEnd && jInc>0) || (j>=jEnd && jInc<0)) ; j += jInc) {
        y = dy*j;
        for (k=kBegin; ((k<=kEnd && kInc>0) || (k>=kEnd && kInc<0)) ; k += kInc) {
          z = dz*k;

          if (qed.getvv(i,j,k)>0.0) {
            ir = (int)(3*inc*viewScale);
            plotOval(x,y,z, ir, 0.85f );
            continue;
          }

          if (mode==0) { // psi^2
            dens = 20000.0*dv*qed.getPh2(i,j,k); if (dens>0.999) dens = 0.999;
            ir = (int)(10.0*viewScale*dens+0.5);
            if (ir>0) plotOval(x,y,z, ir, (float)(0.66-0.66*dens) );

          } else if (mode==1) { // phase
            v = qed.getPhase(i,j,k);
            dens = 1000.0*dv*v[0]; if (dens>0.999) dens = 0.999;
            ir = (int)(3.3*viewScale*dens+0.5);
            if (ir>0) plotOval(x,y,z, ir, (float)v[1] );

          } else if (mode==2) { // flow
            plotFlow(i,j,k);

          } else if (mode==3) { // cloud
            if (cloud[i][j][k]>0) {
              ir = cloud[i][j][k]+2;
              plotOval(x,y,z, ir, 0.66f);
            }
          }
        }
      }
    }

    drawBox(1);
  }

  private void plotOval(double x, double y, double z, int ir, float col) {
    int ix, iy;
    double ppx,ppy,ppz;

    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)col, 0.9f,0.9f));
    gOff.fillOval(ix-ir/2,iy-ir/2,ir,ir);
  }

  private void plotFlow(int i, int j, int k) {
    int ix,iy, ix2,iy2, ir;
    double x,y,z, ppx,ppy,ppz, mag, dens;
    double v[];
 
    if ( i-3<0 || i+3>NNx || j-3<0 || j+3>NNy ||k-3<0 || k+3>NNz ) return;

    mag = 2000000.0;
    v = qed.getpp(i,j,k);

    x = i*dx; y = j*dy; z = k*dz;
    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;
    dens = 50000.0*dv*v[0]; if (dens>0.999) dens = 0.999;
    ir = (int)(2.0*viewScale*dens+0.9);
    if (ir==0) return;

    x += mag*dv*v[1]*dt; y += mag*dv*v[2]*dt; z += mag*dv*v[3]*dt;
    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;
    ix2 = (int)(scale*ppx)+xImageLoc;
    iy2 = (int)(scale*ppy)+yImageLoc;
 
    if (v[1]>0.0) {
      gOff.setColor(Color.getHSBColor(0.66f,(float)(0.3+0.7*dens),0.9f));
    } else {
      gOff.setColor(Color.getHSBColor(0.01f,(float)(0.3+0.7*dens),0.9f));
    }
    gOff.drawLine(ix,iy,ix2,iy2);

    gOff.fillOval(ix-ir/2,iy-ir/2,ir,ir);
  }

  //

  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.0*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;
    }  
  }

  //

  void gridPlot(int mode, double mag) {
    int i,j,k;
    double x,y,z,ppx,ppy,ppz;
    double zv;

    k = kPos;

    for(j=0; j<NNy; j += 2) {
      for (i=0; i<NNx; i++) {
        x = i*dx; y = j*dy;
        zv = getValue(i,j,k);
        z = mag*zv + (kPos-32)*dz + czz;
        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;
        xpts[i] = (int)(scale*ppx)+xImageLoc;
        ypts[i] = (int)(scale*ppy)+yImageLoc;
        if (qed.getvv(i,j,k)>0) {
          colz[i] = 0.85f;
        } else {
          colz[i] = colzv(zv);
        }
      }
      for (i=0; i<NNx-1; i++) {
        gOff.setColor(Color.getHSBColor(colz[i],0.6f,0.9f));
        gOff.drawLine(xpts[i],ypts[i],xpts[i+1],ypts[i+1]);
      }
    }

    for(i=0; i<NNx; i += 2) {
      for (j=0; j<NNy; j++) {
        x = i*dx; y = j*dy;
        zv = getValue(i,j,k);
        z = mag*zv + (kPos-32)*dz + czz;
        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;
        xpts[j] = (int)(scale*ppx)+xImageLoc;
        ypts[j] = (int)(scale*ppy)+yImageLoc;
        if (qed.getvv(i,j,k)>0) {
          colz[j] = 0.85f;
        } else {
          colz[j] = colzv(zv);
        }
      }
      for (j=0; j<NNy-1; j++) { 
        gOff.setColor(Color.getHSBColor(colz[j],0.6f,0.9f));
        gOff.drawLine(xpts[j],ypts[j],xpts[j+1],ypts[j+1]);
      }
    }
  }

  private double getValue(int i, int j, int k) {
    return 200000.0*dv*qed.getPh2(i,j,k)+1.0*qed.getvv(i,j,k);
  }

  private float colzv(double zv) {
    float col;

    col = 0.66f;
    if (zv>0.2) col = 0.15f;
    return col;
  }

  //

  private void densityPlot(int mode) {
    int i,j,k,ix,iy,ir;
    double x,y,z,zv,den;
    float col;
    double v[];

    k = kPos;
    ir = (int)(scale*dx+0.999);

    for(i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        x = i*dx; y = j*dy;
        ix = (int)(scale*x)+xImageLoc;
        iy = (int)(scale*y)+yImageLoc;
        if (qed.getvv(i,j,k)>0) {
          gOff.setColor(Color.getHSBColor(0.85f,0.6f,0.8f));
          gOff.fillRect(ix-ir/2,iy-ir/2,ir,ir);
        } else {
          if (mode==0) { // density
            zv = 30000.0*dv*qed.getPh2(i,j,k); if (zv>0.999) zv = 0.999;
            col = (float)(0.66-0.66*zv);
            den = 15.0*zv; if (den>0.999) den=0.999;
            gOff.setColor(Color.getHSBColor(col,0.6f,(float)den));
          } else if (mode==1) { // phase
            v = qed.getPhase(i,j,k);
            den = 1500.0*dv*v[0]; if (den>0.999) den = 0.999;
            gOff.setColor(Color.getHSBColor((float)v[1],0.9f,(float)den));
          }
          gOff.fillRect(ix-ir/2,iy-ir/2,ir,ir);
        }
      }
    }
  }

  //

  private void setCloud() {
    int i,j,k,m;
    double s;

    for (i=0;i<NNr;i++) {
      srh[i] = Math.random();
    }

    for (i=0;i<NNr;i++) {
      for (j=0;j<NNr-1;j++) {
        if (srh[j]>srh[j+1]) {
          s = srh[j+1]; srh[j+1] = srh[j]; srh[j] = s;
        }
      }
    }

    m = 0;
    s = 0.0;
    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        for (k=0;k<NNz;k++) {
          cloud[i][j][k] = 0;
          s += qed.getPh2(i,j,k)*dv;
          while (s>srh[m] && m<NNr) {
            m += 1;
            cloud[i][j][k] +=1;
          }
        }
      }
    }
  }


  // ------------------------------------- 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 setkPos(int kp) { kPos = kp; }
}



// ================ quantum electron dynamics 3D class ===========

class PQM3D extends Thread {

  private double t = 0.0;

  public static final double dx = 0.25;
  public static final double dy = 0.25;
  public static final double dz = 0.25;
  public static final double dt = 0.25*(dx*dx);

  public static final int NNx = 128;
  public static final int NNy = 64;
  public static final int NNz = 64;

  private double psi[][][][] = new double[NNx][NNy][NNz][2];
  private double vv[][][] = new double[NNx][NNy][NNz];
  private double wrk[][][][] = new double[NNx][NNy][NNz][2];

  private double epsilonX = dt/(dx*dx);
  private double alphaXRe = 0.5 * (1.0+Math.cos(epsilonX/2));
  private double alphaXIm = -0.5 * Math.sin(epsilonX/2);
  private double betaXRe = Math.sin(epsilonX/4)*Math.sin(epsilonX/4);
  private double betaXIm = 0.5 * Math.sin(epsilonX/2);
  
  private double epsilonY = dt/(dy*dy);
  private double alphaYRe = 0.5 * (1.0+Math.cos(epsilonY/2));
  private double alphaYIm = -0.5 * Math.sin(epsilonY/2);
  private double betaYRe = Math.sin(epsilonY/4)*Math.sin(epsilonY/4);
  private double betaYIm = 0.5 * Math.sin(epsilonY/2);

  private double epsilonZ = dt/(dz*dz);
  private double alphaZRe = 0.5 * (1.0+Math.cos(epsilonZ/2));
  private double alphaZIm = -0.5 * Math.sin(epsilonZ/2);
  private double betaZRe = Math.sin(epsilonZ/4)*Math.sin(epsilonZ/4);
  private double betaZIm = 0.5 * Math.sin(epsilonZ/2);

  private double px0 = 1.0;
  private double py0 = 0.0;
  private double pz0 = 0.0;
  private double slitWidth = 2.0;
  private double slitHight = 8.0;
  private double slitDistance = 5.0;

  private int resetSW = 0;
  private int startSW = 1;
  private int stepSW = 0;
  private int lossSW = 0;
  private int sleepTime = 1;
  private int stopSleepTime = 100;
  private int count = 0;


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

  PQM3D() {
    setInitialCondition();
  }


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

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


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

  private void setInitialCondition() {

    t = 0.0;
    count = 0;

    setPlaneWave(NNx*dx/4.0, 2.0, px0);
    setPotential(NNx*dx/2.0,NNy*dy/2.0,NNz*dz/2.0,2.0);
  }

  //

  void setPotential(double x0, double y0, double z0, double r0) {
    int i,j,k;
    double x,y,z;

    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        for (k=0;k<NNz;k++) {
          x = i*dx;
          y = j*dy;
          z = k*dz;
          vv[i][j][k] = 0.0;
          if ((x-x0)*(x-x0)+(y-y0)*(y-y0)+(z-z0)*(z-z0)<r0*r0) {
            vv[i][j][k] = 5.0;
          }
        }
      }
    }

  }

  //

  private void setWave(double xPos, double yPos, double zPos, double waveWidth, double kx, double ky, double kz) {
    int i,j,k;
    double a,x,y,z,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;
          phAb = Math.exp(-((x-xPos)*(x-xPos)+(y-yPos)*(y-yPos)+(z-zPos)*(z-zPos))/(4.0*waveWidth*waveWidth) );
          phPh = kx*x+ky*y+kz*z;
          psi[i][j][k][0] = phAb*Math.cos(phPh);
          psi[i][j][k][1] = phAb*Math.sin(phPh);
        }
      }
    }

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

    normalize(psi);
  }

  void setPlaneWave(double xPos, double waveWidth, double kx) {
    int i,j,k;
    double a,x,y,z,r,phAb,phPh;

    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        for (k=0;k<NNz;k++) {
          x = i*dx;
          y = j*dy;
          z = k*dz;
          phAb = Math.exp(-((x-xPos)*(x-xPos))/(4.0*waveWidth*waveWidth) );
          phPh = kx*x;
          psi[i][j][k][0] = phAb*Math.cos(phPh);
          psi[i][j][k][1] = phAb*Math.sin(phPh);
        }
      }
    }

    normalize(psi);
  }


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

  private void timeEvolution(){

    if (resetSW==1) {
      setInitialCondition();
      resetSW = 0;
    }

    if (startSW==1) {
      timeStep();

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

  private void timeStep() {
    int i;

    for (i=0; i<2; i++) {
      t = t + dt;
      phaseStep(psi,vv,0.5*dt);
      kxStep(psi,dt);
      kyStep(psi,dt);
      kzStep(psi,dt);
      phaseStep(psi,vv,0.5*dt);
    }
    if (lossSW==1) {
      steepestDescent(psi,vv,0.02);
    }
  }

  //

  void kxStep(double ph[][][][], double deltat) {
    int i,j,k;

    for(j=0;j<NNy;j++){
      for(k=0;k<NNz;k++){
        for(i=0;i<NNx-1;i+=2){
         substepX(ph,i,i+1,j,k);
        }
        for(i=1;i<NNx-1;i+=2){
          substepX(ph,i,i+1,j,k);
        }
        substepX(ph,NNx-1,0,j,k);

        substepX(ph,NNx-1,0,j,k);
        for(i=1;i<NNx-1;i+=2){
          substepX(ph,i,i+1,j,k);
        }
        for(i=0;i<NNx-1;i+=2){
          substepX(ph,i,i+1,j,k);
        }
      }
    }
  } 

  void substepX(double ph[][][][], int i,int i2, int j, int k) {
    double xRe,xIm,yRe,yIm,wRe,wIm,zRe,zIm;

    xRe = ph[i][j][k][0]; xIm = ph[i][j][k][1];
    yRe = ph[i2][j][k][0]; yIm = ph[i2][j][k][1];
    wRe = alphaXRe*xRe-alphaXIm*xIm; wIm = alphaXRe*xIm+alphaXIm*xRe;
    zRe = betaXRe*yRe-betaXIm*yIm; zIm = betaXRe*yIm+betaXIm*yRe;
    ph[i][j][k][0] = wRe + zRe; ph[i][j][k][1] = wIm + zIm;
    wRe = alphaXRe*yRe-alphaXIm*yIm; wIm = alphaXRe*yIm+alphaXIm*yRe;
    zRe = betaXRe*xRe-betaXIm*xIm; zIm = betaXRe*xIm+betaXIm*xRe;
    ph[i2][j][k][0] = wRe + zRe; ph[i2][j][k][1] = wIm + zIm;

  }

  void kyStep(double ph[][][][], double deltat) {
    int i,j,k;

    for(i=0;i<NNx;i++){
      for(k=0;k<NNz;k++){
        for(j=0;j<NNy-1;j+=2){
         substepY(ph,i,j,j+1,k);
        }
        for(j=1;j<NNy-1;j+=2){
          substepY(ph,i,j,j+1,k);
        }
        substepY(ph,i,NNy-1,0,k);

        substepY(ph,i,NNy-1,0,k);
        for(j=1;j<NNy-1;j+=2){
          substepY(ph,i,j,j+1,k);
        }
        for(j=0;j<NNy-1;j+=2){
          substepY(ph,i,j,j+1,k);
        }
      }
    }
  } 

  void substepY(double ph[][][][], int i,int j, int j2, int k) {
    double xRe,xIm,yRe,yIm,wRe,wIm,zRe,zIm;

    xRe = ph[i][j][k][0]; xIm = ph[i][j][k][1];
    yRe = ph[i][j2][k][0]; yIm = ph[i][j2][k][1];
    wRe = alphaYRe*xRe-alphaYIm*xIm; wIm = alphaYRe*xIm+alphaYIm*xRe;
    zRe = betaYRe*yRe-betaYIm*yIm; zIm = betaYRe*yIm+betaYIm*yRe;
    ph[i][j][k][0] = wRe + zRe; ph[i][j][k][1] = wIm + zIm;
    wRe = alphaYRe*yRe-alphaYIm*yIm; wIm = alphaYRe*yIm+alphaYIm*yRe;
    zRe = betaYRe*xRe-betaYIm*xIm; zIm = betaYRe*xIm+betaYIm*xRe;
    ph[i][j2][k][0] = wRe + zRe; ph[i][j2][k][1] = wIm + zIm;

  }

  void kzStep(double ph[][][][], double deltat) {
    int i,j,k;

    for(i=0;i<NNx;i++){
      for(j=0;j<NNy;j++){
        for(k=0;k<NNz-1;k+=2){
         substepZ(ph,i,j,k,k+1);
        }
        for(k=1;k<NNz-1;k+=2){
          substepZ(ph,i,j,k,k+1);
        }
        substepZ(ph,i,j,NNz-1,0);

        substepZ(ph,i,j,NNz-1,0);
        for(k=1;k<NNz-1;k+=2){
          substepZ(ph,i,j,k,k+1);
        }
        for(k=0;k<NNz-1;k+=2){
          substepZ(ph,i,j,k,k+1);
        }
      }
    }
  } 

  void substepZ(double ph[][][][], int i,int j, int k, int k2) {
    double xRe,xIm,yRe,yIm,wRe,wIm,zRe,zIm;

    xRe = ph[i][j][k][0]; xIm = ph[i][j][k][1];
    yRe = ph[i][j][k2][0]; yIm = ph[i][j][k2][1];
    wRe = alphaZRe*xRe-alphaZIm*xIm; wIm = alphaZRe*xIm+alphaZIm*xRe;
    zRe = betaZRe*yRe-betaZIm*yIm; zIm = betaZRe*yIm+betaZIm*yRe;
    ph[i][j][k][0] = wRe + zRe; ph[i][j][k][1] = wIm + zIm;
    wRe = alphaZRe*yRe-alphaZIm*yIm; wIm = alphaZRe*yIm+alphaZIm*yRe;
    zRe = betaZRe*xRe-betaZIm*xIm; zIm = betaZRe*xIm+betaZIm*xRe;
    ph[i][j][k2][0] = wRe + zRe; ph[i][j][k2][1] = wIm + zIm;

  }

  private void phaseStep(double ph[][][][], double vv[][][], double deltat) {
    int i,j,k;
    double th,cs,sn,phr,phi;

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

  // set loss

  void steepestDescent(double ph[][][][], double v[][][], double dump) {
    int i,j,k,ip,jp,kp,im,jm,km;
    double ee,d2phRe,d2phIm,h2;

    h2 = 2.0*dx*dx;
    ee = meanKinetic(ph)+meanPotential(ph,vv);

    for (i=0;i<NNx;i++) {
      ip = i+1; if (ip>=NNx) ip = 0;
      im = i-1; if (im<0) im = NNx-1;
      for (j=0;j<NNy;j++) {
        jp = j+1; if (jp>=NNy) jp = 0;
        jm = j-1; if (jm<0) jm = NNy-1;
        for (k=0;k<NNz;k++) {
          kp = k+1; if (kp>=NNz) kp = 0;
          km = k-1; if (km<0) km = NNz-1;
          wrk[i][j][k][0] = -(ph[ip][j][k][0]+ph[im][j][k][0]+ph[i][jp][k][0]+ph[i][jm][k][0]+ph[i][j][kp][0]+ph[i][j][km][0]-6.0*ph[i][j][k][0])/h2+(v[i][j][k]-ee)*ph[i][j][k][0];
          wrk[i][j][k][1] = -(ph[ip][j][k][1]+ph[im][j][k][1]+ph[i][jp][k][1]+ph[i][jm][k][1]+ph[i][j][kp][1]+ph[i][j][km][1]-6.0*ph[i][j][k][1])/h2+(v[i][j][k]-ee)*ph[i][j][k][1];
        }
      }
    }

    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        for (k=0;k<NNz;k++) {
          ph[i][j][k][0] -= dump*wrk[i][j][k][0];
          ph[i][j][k][1] -= dump*wrk[i][j][k][1];
        }
      }
    }

    normalize(ph);
  }


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

  private double norm(double ph[][][][]) {
    int i,j,k;
    double p;

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

  private void normalize(double ph[][][][]) {
    int i,j,k;
    double a;

    a = 1.0/Math.sqrt(norm(ph));
	
    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        for (k=0;k<NNz;k++) {
          ph[i][j][k][0] *= a;
          ph[i][j][k][1] *= a;
        }
      }
    }
  }

  private double meanPotential(double ph[][][][], double vv[][][]) {
    int i,j,k;
    double p;

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

  private double meanKinetic(double ph[][][][]) {
    int i,j,k,ip,jp,kp,im,jm,km;
    double p,d2phRe,d2phIm,h2;

    h2 = dx*dx;
    p = 0.0;
    for (i=0;i<NNx;i++) {
      ip = i+1; if (ip>=NNx) ip = 0;
      im = i-1; if (im<0) im = NNx-1;
      for (j=0;j<NNy;j++) {
        jp = j+1; if (jp>=NNy) jp = 0;
        jm = j-1; if (jm<0) jm = NNy-1;
        for (k=0;k<NNz;k++) {
          kp = k+1; if (kp>=NNz) kp = 0;
          km = k-1; if (km<0) km = NNz-1;
          d2phRe = (ph[ip][j][k][0]+ph[im][j][k][0]+ph[i][jp][k][0]+ph[i][jm][k][0]+ph[i][j][kp][0]+ph[i][j][km][0]-6.0*ph[i][j][k][0])/h2;
          d2phIm = (ph[ip][j][k][1]+ph[im][j][k][1]+ph[i][jp][k][1]+ph[i][jm][k][1]+ph[i][j][kp][1]+ph[i][j][km][1]-6.0*ph[i][j][k][1])/h2;
          p += (ph[i][j][k][0]*d2phRe+ph[i][j][k][1]*d2phIm);
        }
      }
    }
    return ( -0.5*p*dx*dy*dz );
  }


  // ------------------------------------- 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 int getCount() { return count; }

  public void setpx0(double p) { px0 = p; }
  public double getpx0() { return px0; }
  public void setpy0(double p) { py0 = p; }
  public double getpy0() { return py0; }
  public void setpz0(double p) { pz0 = p; }
  public double getpz0() { return pz0; }

  // get data
  public double getTime() { return t; }

  public double getPh2(int i, int j, int k) {
    return (psi[i][j][k][0]*psi[i][j][k][0]+psi[i][j][k][1]*psi[i][j][k][1]);
  }

  public double getvv(int i, int j, int k) { return vv[i][j][k]; }

  public double[] getPhase(int i, int j, int k) {
    double v[] = new double[2];

    v[0] = Math.sqrt(psi[i][j][k][0]*psi[i][j][k][0]+psi[i][j][k][1]*psi[i][j][k][1]);
    v[1] = 0.33+0.33*psi[i][j][k][0]/v[0];
    return v;
  }

  public double[] getpp(int i, int j, int k) {
    int ip,jp,kp,im,jm,km;
    double pRe,pIm;
    double v[] = new double[4];

    ip = i+1; if (ip>=NNx) ip = 0;
    im = i-1; if (im<0) im = NNx-1;
    jp = j+1; if (jp>=NNy) jp = 0;
    jm = j-1; if (jm<0) jm = NNy-1;
    kp = k+1; if (kp>=NNz) kp = 0;
    km = k-1; if (km<0) km = NNz-1;
	
    v[0] = psi[i][j][k][0]*psi[i][j][k][0]+psi[i][j][k][1]*psi[i][j][k][1];
    pRe = (psi[ip][j][k][1]-psi[im][j][k][1])/(2*dx);
    pIm = (-psi[ip][j][k][0]+psi[im][j][k][0])/(2*dx);
    v[1] = psi[i][j][k][0]*pRe + psi[i][j][k][1]*pIm;
    pRe = (psi[i][jp][k][1]-psi[i][jm][k][1])/(2*dy);
    pIm = (-psi[i][jp][k][0]+psi[i][jm][k][0])/(2*dy);
    v[2] = (psi[i][j][k][0]*pRe + psi[i][j][k][1]*pIm);
    pRe = (psi[i][j][kp][1]-psi[i][j][km][1])/(2*dz);
    pIm = (-psi[i][j][kp][0]+psi[i][j][km][0])/(2*dz);
    v[3] = (psi[i][j][k][0]*pRe + psi[i][j][k][1]*pIm);
    return v;
  }

  public double getNorm() { return norm(psi); }
  public double getMeanPotential() { return meanPotential(psi,vv); }
  public double getMeanKinetic() { return meanKinetic(psi); }
}