/** applet No. 1118
 *
 * H ion scatter - quantum electron dynamics 3D
 * - multi thread  QED3D :: display - asynchronous
 *
 * Created by Ikeuchi Mitsuru on February 04 2007.
 * Copyright (c) 2007 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2007.02.04  created
 * ver 0.0.2   2007.06.16  improved code
 *
 */

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

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

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

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

  // for event
  Choice ch_view;
  Button bt_reset, bt_startStop, bt_step, bt_loss;
  Scrollbar sc_px0, sc_sWidth, sc_sDistance, 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 psi");
    ch_view.add("psi^2 2D");
    ch_view.add("phase 2D");
    ch_view.add("grid V");
    ch_view.add("grid psi+V");
    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);

    bt_loss = new Button("loss on/off");
    bt_loss.addActionListener(this);

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

    sc_sWidth = new Scrollbar(Scrollbar.HORIZONTAL,20,10,0,60);
    sc_sWidth.addAdjustmentListener(this);

    sc_sDistance = new Scrollbar(Scrollbar.HORIZONTAL,50,10,0,110);
    sc_sDistance.addAdjustmentListener(this);

    sc_kPos = new Scrollbar(Scrollbar.HORIZONTAL,36,10,0,81);
    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(bt_loss);
    pnl.add(new Label("  "));
    pnl.add(sc_kPos);
    pnl.add(ch_view);
    add(pnl,"North");
  }

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


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

  public void stop() {
    if (th != null) th = null;
    if (qed != null) qed = 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){ 
      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();
      }
    } else if (ev.getSource() == bt_loss){
      if (qed.getLossSW()==0) {
        qed.setLossSW(1);
      } else {
        qed.setLossSW(0);
      }
    }
  }

  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_sWidth) { 
      qed.setSlitWidth( 0.1*(double)(sc_sWidth.getValue()) );
    } else if (ev.getSource() == sc_sDistance) { 
      qed.setSlitDistance( 0.1*(double)(sc_sDistance.getValue()) );
    } 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 QED3D qed;

  private double dt = QED3D.dt;
  private int NNx = QED3D.NNx;
  private int NNy = QED3D.NNy;
  private int NNz = QED3D.NNz;
  private double dx = QED3D.dx;
  private double dy = QED3D.dy;
  private double dz = QED3D.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 = 36;

  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, QED3D 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(1.0, 0.0);

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

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

    } else if (dispMode==7) { // grid potential
      gridPlot(0.0, 1.0);

    } else if (dispMode==8) { // grid psi^2 + potential
      gridPlot(1.0, 1.0);

    } else if (dispMode==9) {
      ;

    }

    gOff.setColor(Color.black);
    gOff.drawString("px0="+(float)(qed.getpx0())+" ",630/6*0+10,70);
    if (qed.getLossSW()==1) {
      gOff.drawString("loss = ON",630/6*2+10,70);
    } else {
      gOff.drawString("loss = OFF",630/6*2+10,70);
    }
    gOff.setColor(Color.black);
    gOff.drawString("2D z="+(kPos-36)*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 = "+(int)(NNx*dx*10.0)/10.0+" x "+(int)(NNy*dy*10.0)/10.0+" x "+(int)(NNz*dz*10.0)/10.0+" au",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("QED3D= "+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;
          }
*/
          dens =  qed.getvv(i,j,k);
		  if (i%4==0 && j%4==0 && k%4==0 && dens<-0.2) {
            ir = (int)(-dens*5.0*viewScale);
            plotOval(x,y,z, ir, (float)(0.75-0.01*dens) );
            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 = 500000.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 = 20000.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(double ph2Mag, double vvMag) {
    int i,j,k;
    double x,y,z,ppx,ppy,ppz;
    double ph2,zv;

    k = kPos;

    for(j=0; j<NNy; j += 2) {
      for (i=0; i<NNx; i++) {
        x = i*dx; y = j*dy;
		ph2 = 100000.0*dv*qed.getPh2(i,j,k);
		zv = 3.0*qed.getvv(i,j,k);
        z = ph2Mag*ph2+vvMag*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(ph2);
        }
      }
      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;
		ph2 = 100000.0*dv*qed.getPh2(i,j,k);
		zv = 3.0*qed.getvv(i,j,k);
        z = ph2Mag*ph2+vvMag*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(ph2);
        }
      }
      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 float colzv(double ph2) {
    float col;

    col = 0.66f;
    if (ph2>0.1) 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 = 10000.0*dv*qed.getPh2(i,j,k); if (zv>0.999) zv = 0.999;
            col = (float)(0.66-0.66*zv);
            den = 5.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 = 500.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 QED3D 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 = 1.0*(dx*dx);

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

  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 int NN = 200;
  private double bRe[] = new double[NN];
  private double bIm[] = new double[NN];
  private double uRe[] = new double[NN];
  private double uIm[] = new double[NN];

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

  QED3D() {
    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;

    setWave(8.0,9.0,9.0, 2.0, px0,py0,pz0);
    setPotential(20.0,9.0,9.0);
  }

  //

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

    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.2) {
            r = 0.2;
          }
          vv[i][j][k] = -1.0/r;

        }
      }
    }
  }

  //

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


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

  //

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

    a = 4.0*dy*dy/dt;
    aaAb = 4.0+a*a;
    for (k=1; k<NNz-1; k++) {
      for (j=1; j<NNy-1; j++) {
        for (i=1; i<NNx-1; i++) {
          bRe[i] = 2.0*ph[i][j][k][0]-a*ph[i][j][k][1] - ph[i+1][j][k][0] - ph[i-1][j][k][0];
          bIm[i] = 2.0*ph[i][j][k][1]+a*ph[i][j][k][0] - ph[i+1][j][k][1] - ph[i-1][j][k][1];
        }

        uRe[1] = -2.0/aaAb;
        uIm[1] = -a/aaAb;
        ph[1][j][k][0] = bRe[1]*uRe[1] - bIm[1]*uIm[1];
        ph[1][j][k][1] = bIm[1]*uRe[1] + bRe[1]*uIm[1];

        for (i=2; i<NNx-1; i++) {
          auAb = (-2.0-uRe[i-1])*(-2.0-uRe[i-1])+(a-uIm[i-1])*(a-uIm[i-1]);
          uRe[i] = (-2.0-uRe[i-1])/auAb;
          uIm[i] = -(a-uIm[i-1])/auAb;
          ph[i][j][k][0] = (bRe[i]-ph[i-1][j][k][0])*uRe[i] - (bIm[i]-ph[i-1][j][k][1])*uIm[i];
          ph[i][j][k][1] = (bRe[i]-ph[i-1][j][k][0])*uIm[i] + (bIm[i]-ph[i-1][j][k][1])*uRe[i];
        }

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

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

    a = 4.0*dy*dy/dt;
    aaAb = 4.0+a*a;
    for (i=1; i<NNx-1; i++) {
      for (k=1; k<NNz-1; k++) {
        for (j=1; j<NNy-1; j++) {
          bRe[j] = 2.0*ph[i][j][k][0]-a*ph[i][j][k][1] - ph[i][j+1][k][0] - ph[i][j-1][k][0];
          bIm[j] = 2.0*ph[i][j][k][1]+a*ph[i][j][k][0] - ph[i][j+1][k][1] - ph[i][j-1][k][1];
        }

        uRe[1] = -2.0/aaAb;
        uIm[1] = -a/aaAb;
        ph[i][1][k][0] = bRe[1]*uRe[1] - bIm[1]*uIm[1];
        ph[i][1][k][1] = bIm[1]*uRe[1] + bRe[1]*uIm[1];

        for (j=2; j<NNy-1; j++) {
          auAb = (-2.0-uRe[j-1])*(-2.0-uRe[j-1])+(a-uIm[j-1])*(a-uIm[j-1]);
          uRe[j] = (-2.0-uRe[j-1])/auAb;
          uIm[j] = -(a-uIm[j-1])/auAb;
          ph[i][j][k][0] = (bRe[j]-ph[i][j-1][k][0])*uRe[j] - (bIm[j]-ph[i][j-1][k][1])*uIm[j];
          ph[i][j][k][1] = (bRe[j]-ph[i][j-1][k][0])*uIm[j] + (bIm[j]-ph[i][j-1][k][1])*uRe[j];
        }

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

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

    a = 4.0*dy*dy/dt;
    aaAb = 4.0+a*a;
    for (i=1; i<NNx-1; i++) {
      for (j=1; j<NNy-1; j++) {
        for (k=1; k<NNz-1; k++) {
          bRe[k] = 2.0*ph[i][j][k][0]-a*ph[i][j][k][1] - ph[i][j][k+1][0] - ph[i][j][k-1][0];
          bIm[k] = 2.0*ph[i][j][k][1]+a*ph[i][j][k][0] - ph[i][j][k+1][1] - ph[i][j][k-1][1];
        }

        uRe[1] = -2.0/aaAb;
        uIm[1] = -a/aaAb;
        ph[i][j][1][0] = bRe[1]*uRe[1] - bIm[1]*uIm[1];
        ph[i][j][1][1] = bIm[1]*uRe[1] + bRe[1]*uIm[1];

        for (k=2; k<NNz-1; k++) {
          auAb = (-2.0-uRe[k-1])*(-2.0-uRe[k-1])+(a-uIm[k-1])*(a-uIm[k-1]);
          uRe[k] = (-2.0-uRe[k-1])/auAb;
          uIm[k] = -(a-uIm[k-1])/auAb;
          ph[i][j][k][0] = (bRe[k]-ph[i][j][k-1][0])*uRe[k] - (bIm[k]-ph[i][j][k-1][1])*uIm[k];
          ph[i][j][k][1] = (bRe[k]-ph[i][j][k-1][0])*uIm[k] + (bIm[k]-ph[i][j][k-1][1])*uRe[k];
        }

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

  }

  private void phaseStep(double ph[][][][], double vv[][][], double deltat) {
    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 = 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;
    double ee,d2phRe,d2phIm,h2;

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

    for (i=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        for (k=1;k<NNz-1;k++) {
          wrk[i][j][k][0] = -(ph[i+1][j][k][0]+ph[i-1][j][k][0]+ph[i][j+1][k][0]+ph[i][j-1][k][0]+ph[i][j][k+1][0]+ph[i][j][k-1][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[i+1][j][k][1]+ph[i-1][j][k][1]+ph[i][j+1][k][1]+ph[i][j-1][k][1]+ph[i][j][k+1][1]+ph[i][j][k-1][1]-6.0*ph[i][j][k][1])/h2+(v[i][j][k]-ee)*ph[i][j][k][1];
        }
      }
    }

    for (i=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        for (k=1;k<NNz-1;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=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        for (k=1;k<NNz-1;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=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        for (k=1;k<NNz-1;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=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        for (k=1;k<NNz-1;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;
    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 = (ph[i+1][j][k][0]+ph[i-1][j][k][0]+ph[i][j+1][k][0]+ph[i][j-1][k][0]+ph[i][j][k+1][0]+ph[i][j][k-1][0]-6.0*ph[i][j][k][0])/h2;
          d2phIm = (ph[i+1][j][k][1]+ph[i-1][j][k][1]+ph[i][j+1][k][1]+ph[i][j-1][k][1]+ph[i][j][k+1][1]+ph[i][j][k-1][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 void setLossSW(int sw) { lossSW = sw; }
  public int getLossSW() { return lossSW; }
  public int getCount() { return count; }

  public void setpx0(double p) { px0 = p; }
  public double getpx0() { return px0; }
  public void setSlitWidth(double w) { slitWidth = w; }
  public double getSlitWidth() { return slitWidth; }
  public void setSlitDistance(double d) { slitDistance = d; }
  public double getSlitDistance() { return slitDistance; }

  // 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) {
    double pRe,pIm;
    double v[] = new double[4];

    v[0] = psi[i][j][k][0]*psi[i][j][k][0]+psi[i][j][k][1]*psi[i][j][k][1];
    pRe = (psi[i+1][j][k][1]-psi[i-1][j][k][1])/(2*dx);
    pIm = (-psi[i+1][j][k][0]+psi[i-1][j][k][0])/(2*dx);
    v[1] = psi[i][j][k][0]*pRe + psi[i][j][k][1]*pIm;
    pRe = (psi[i][j+1][k][1]-psi[i][j-1][k][1])/(2*dy);
    pIm = (-psi[i][j+1][k][0]+psi[i][j-1][k][0])/(2*dy);
    v[2] = (psi[i][j][k][0]*pRe + psi[i][j][k][1]*pIm);
    pRe = (psi[i][j][k+1][1]-psi[i][j][k-1][1])/(2*dz);
    pIm = (-psi[i][j][k+1][0]+psi[i][j][k-1][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); }
}