/* template QM2D                */
/* coded by Ikeuchi Mitsuru     */
/* ver 0.0.1   2004.05.08       */

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

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

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

  Choice cvw;
  Button bt_reset, bt_start, bt_stop;
  Scrollbar spx,spy;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 30;

  int dispMode = 0;
  int started = 1;

  double px0 = 24.0;
  double py0 = 8.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/256.0;
  double dy = 8.0/256.0;

  double dt = 0.5*dx*dx;

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


  int NNx = 256;
  int NNy = 256;
  int NN = NNx+1; /* max(NNx,NNy)+1 */
  int NNr = 200;

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

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

  double srh[] = new double[NNr+1];
  int xptr[] = new int[NNr+1];
  int yptr[] = new int[NNr+1];

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

  public void init() {

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

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

    cvw = new Choice();
    cvw.addItem("density"); cvw.addItem("grid ph"); cvw.addItem("grid V");
    cvw.addItem("grid ph+V"); cvw.addItem("phase"); cvw.addItem("cloud");
    cvw.addItemListener(this);

    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);
    
    spx= new Scrollbar(Scrollbar.HORIZONTAL,240,10,0,310);
    spx.addAdjustmentListener(this);

    spy= new Scrollbar(Scrollbar.HORIZONTAL,80,10,0,310);
    spy.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(spx);
    pnl.add(spy);
    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() == spx) { 
      px0 = 0.1*(double)(spx.getValue());
    } else if (ev.getSource() == spy) { 
      py0 = 0.1*(double)(spy.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() {
    double ev,ek;

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

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


    if (dispMode==0) {

      densityPlot(2.0);

    } else if (dispMode==1) {

      gridPlot(50.0, 0.0);

    } else if (dispMode==2) {

      gridPlot(0.0, 0.1);

    } else if (dispMode==3) {

      gridPlot(50.0, 0.1);

    } else if (dispMode==4) {

      phasePlot(4.0);

    } else if (dispMode==5) {

      cloudPlot();

    }

    gOff.setColor(Color.black);
    gOff.drawString("view",630/6*3+10,40);
    gOff.drawString("px0="+(int)(px0*10.0+0.5)/10.0+" au",630/6*4+10,40);
    gOff.drawString("py0="+(int)(py0*10.0+0.5)/10.0+" au",630/6*5+10,40);

    gOff.drawString("Lx="+(int)(NNx*dx*52.918+0.5)+" pm",630/6*4+10,60);
    gOff.drawString("Ly="+(int)(NNy*dy*52.918+0.5)+" pm",630/6*5+10,60);

    gOff.setColor(Color.blue);
    gOff.drawString("t="+(int)(t*2418.9+0.5)/100.0+" as",10,40);
    gOff.drawString("norm="+norm()+" ",630/6*1+10,40);

    gOff.drawString("V="+(int)(ev*100+0.5)/100.0+" ",630/6*0+10,55);
    gOff.drawString("K="+(int)(ek*100+0.5)/100.0+" ",630/6*1+10,55);
    gOff.drawString("K+V="+(int)((ev+ek)*100+0.5)/100.0+" ",630/6*2+10,55);

  }

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

  void densityPlot(double weight) {
    int i,j,gbx,gby,mag;
    double ph2;

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

    for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        ph2 = weight*(phRe[i][j]*phRe[i][j]+phIm[i][j]*phIm[i][j]);
        if (ph2>0.9) {
          ph2 = 0.9;
        }
        if (ph2>0.02) {
          gOff.setColor(Color.getHSBColor(0.15f, 0.95f, (float)(ph2+0.1)));
          gOff.fillRect(gbx+mag*i,gby+mag*j, mag, mag);
        } else if (vv[i][j]>0.0) {
          gOff.setColor(Color.getHSBColor(0.33f, 0.95f, 0.6f));
          gOff.fillRect(gbx+mag*i,gby+mag*j, mag, mag);
        }
      }
    }
  }

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

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

    gbx = 0+(int)(cx/sc);
    gby = -100+(int)(cy/sc);

    for(j=0; j<NNy; j += 8) {
      for (i=0; i<NNx; i++) { 
        ph2 = phRe[i][j]*phRe[i][j]+phIm[i][j]*phIm[i][j];
        v = vv[i][j];
        if (v>300.0) {
          v = 300;
        }
        z = phWeight*ph2 + vvWeight*v + 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.02) {
          icol[i] = 0.15;
        } else if (vv[i][j]>0.0) {
          icol[i] = 0.33;
        } else {
          icol[i] = 0.66;
        }
      }
      for (i=0; i<NNx-2; i++) {
        gOff.setColor(Color.getHSBColor((float)icol[i],0.9f,0.9f));
        gOff.drawLine(xpts[i],ypts[i],xpts[i+1],ypts[i+1]);
      }
    }

    for(i=0; i<NNx; i += 8) {
      for (j=0; j<NNy; j++) {
        ph2 = phRe[i][j]*phRe[i][j]+phIm[i][j]*phIm[i][j];
        v = vv[i][j];
        if (v>300.0) {
          v = 300;
        }
        z = phWeight*ph2 + vvWeight*v + 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.02) {
          icol[j] = 0.15;
        } else if (vv[i][j]>0.0) {
          icol[j] = 0.33;
        } else {
          icol[j] = 0.66;
        }
      }
      for (j=0; j<NNy-2; j++) { 
        gOff.setColor(Color.getHSBColor((float)icol[j],0.9f,0.9f));
        gOff.drawLine(xpts[j],ypts[j],xpts[j+1],ypts[j+1]);
      }
    }

  }

  void phasePlot(double weight) {
    int i,j,gbx,gby,mag;
    double ph2,r;

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

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

  void cloudPlot() {
    int i,j,m,gbx,gby,mag;
    double ph2,r;

    setsrh();

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

    for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        if (vv[i][j]>0.0) {
          gOff.setColor(Color.getHSBColor(0.33f, 0.95f, 0.6f));
          gOff.fillRect(gbx+mag*i,gby+mag*j, mag, mag);
        }
      }
    }

    for (m=0; m<NNr; m++) {
      i=xptr[m]; j=yptr[m];
      gOff.setColor(Color.getHSBColor(0.15f, 0.9f, 0.9f));
      gOff.fillRect(gbx+mag*i,gby+mag*j, mag, mag);
    }

  }

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

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

    p = 0.0;
    for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        p += (phRe[i][j]*phRe[i][j]+phIm[i][j]*phIm[i][j])*dx*dy;
      }
    }
    return ( p );
  }

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

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

  double meanKinetic() {
    int i,j;
    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++) {
        d2phRe = (phRe[i+1][j]+phRe[i-1][j]+phRe[i][j+1]+phRe[i][j-1]-4.0*phRe[i][j])/h2;
        d2phIm = (phIm[i+1][j]+phIm[i-1][j]+phIm[i][j+1]+phIm[i][j-1]-4.0*phIm[i][j])/h2;
        p += (phRe[i][j]*d2phRe+phIm[i][j]*d2phIm)*dx*dy;
      }
    }

    for (j=1; j<NNy-1; j++) {
      d2phRe = (phRe[1][j]+phRe[NNx-1][j]+phRe[0][j+1]+phRe[0][j-1]-4.0*phRe[0][j])/h2;
      d2phIm = (phIm[1][j]+phIm[NNx-1][j]+phIm[0][j+1]+phIm[0][j-1]-4.0*phIm[0][j])/h2;
        p += (phRe[0][j]*d2phRe+phIm[0][j]*d2phIm)*dx*dy;
    }

    for (j=1; j<NNy-1; j++) {
        d2phRe = (phRe[0][j]+phRe[NNx-2][j]+phRe[NNx-1][j+1]+phRe[NNx-1][j-1]-4.0*phRe[NNx-1][j])/h2;
        d2phIm = (phIm[0][j]+phIm[NNx-2][j]+phIm[NNx-1][j+1]+phIm[NNx-1][j-1]-4.0*phIm[i][j])/h2;
        p += (phRe[NNx-1][j]*d2phRe+phIm[NNx-1][j]*d2phIm)*dx*dy;
    }

    for (i=1; i<NNx-1; i++) {
      d2phRe = (phRe[i+1][0]+phRe[i-1][0]+phRe[i][1]+phRe[i][NNy-1]-4.0*phRe[i][0])/h2;
      d2phIm = (phIm[i+1][0]+phIm[i-1][0]+phIm[i][1]+phIm[i][NNy-1]-4.0*phIm[i][0])/h2;
      p += (phRe[i][j]*d2phRe+phIm[i][j]*d2phIm)*dx*dy;
    }

    for (i=1; i<NNx-1; i++) {
        d2phRe = (phRe[i+1][NNy-1]+phRe[i-1][NNy-1]+phRe[i][0]+phRe[i][NNy-2]-4.0*phRe[i][NNy-1])/h2;
        d2phIm = (phIm[i+1][NNy-1]+phIm[i-1][NNy-1]+phIm[i][0]+phIm[i][NNy-2]-4.0*phIm[i][NNy-1])/h2;
        p += (phRe[i][j]*d2phRe+phIm[i][j]*d2phIm)*dx*dy;
    }

    return ( -0.5*p );
  }

  void setsrh() {
    int i,j,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++) {
        s += (phRe[i][j]*phRe[i][j]+phIm[i][j]*phIm[i][j])*dx*dy;
        while (s>srh[m] && m<NNr-1) {
          xptr[m] = i; yptr[m] = j;
          m +=1;
        }
      }
    }

  }



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

  void setInitialCondition() {
    t = 0.0;
    setWave(2.0,2.0, 0.5, px0, py0);
    setPotential();

  }

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

  void setPotential() {
    int i,j;
    double wallPos, wallWidth, wallHight;
    double x,y;

    wallPos = 5.0;
    wallWidth = 0.2;
    wallHight = 0.5*24*24*1.1;

    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        x = i*dx;
        y = j*dy;
        if (x>wallPos && x<(wallPos+wallWidth)) {
          vv[i][j] = wallHight;
        } else {
          vv[i][j] = 0.0;
        }
      }
    }

    for (i=0;i<NNx;i++) {
      vv[i][0] = 10000.0;
      vv[i][NNy-1] = 10000.0;
    }
    for (j=0;j<NNy;j++) {
      vv[0][j] = 10000.0;
      vv[NNx-1][j] = 10000.0;
    }

  }

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

  void setWave(double xPos, double yPos, double waveWidth, double kx, double ky) {
    int i,j;
    double a,x,y,phAb,phPh;

    for (i=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        x = i*dx;
        y = j*dy;
        phAb = Math.exp(-((x-xPos)*(x-xPos)+(y-yPos)*(y-yPos))/(4.0*waveWidth*waveWidth) );
        phPh = kx*x+ky*y;
        phRe[i][j] = phAb*Math.cos(phPh);
        phIm[i][j] = phAb*Math.sin(phPh);
      }
    }
    for (i=1;i<NNx-1;i++) {
      phRe[i][0] = 0.0;
      phIm[i][0] = 0.0;
      phRe[i][NNy-1] = 0.0;
      phIm[i][NNy-1] = 0.0;
    }
    for (j=1;j<NNy-1;j++) {
      phRe[0][j] = 0.0;
      phIm[0][j] = 0.0;
      phRe[NNx-1][j] = 0.0;
      phIm[NNx-1][j] = 0.0;
    }

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

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

  void timeEvolution() {
    int i;

    if (started==1) {
      for (i=0; i<5; i++) {
        t = t + dt;
        phaseStep();
        kxStep();
        kyStep();
        phaseStep();
      }
    }
  }

/* ----------------------------- time step -------------- */
  
    void kxStep() {
    int i,j;

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

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

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

    xRe = phRe[i][j]; xIm = phIm[i][j];
    yRe = phRe[i2][j]; yIm = phIm[i2][j];
    wRe = alphaXRe*xRe-alphaXIm*xIm; wIm = alphaXRe*xIm+alphaXIm*xRe;
    zRe = betaXRe*yRe-betaXIm*yIm; zIm = betaXRe*yIm+betaXIm*yRe;
    phRe[i][j] = wRe + zRe; phIm[i][j] = wIm + zIm;
    wRe = alphaXRe*yRe-alphaXIm*yIm; wIm = alphaXRe*yIm+alphaXIm*yRe;
    zRe = betaXRe*xRe-betaXIm*xIm; zIm = betaXRe*xIm+betaXIm*xRe;
    phRe[i2][j] = wRe + zRe; phIm[i2][j] = wIm + zIm;

  }

    void kyStep() {
    int i,j;

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

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

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

    xRe = phRe[i][j]; xIm = phIm[i][j];
    yRe = phRe[i][j2]; yIm = phIm[i][j2];
    wRe = alphaYRe*xRe-alphaYIm*xIm; wIm = alphaYRe*xIm+alphaYIm*xRe;
    zRe = betaYRe*yRe-betaYIm*yIm; zIm = betaYRe*yIm+betaYIm*yRe;
    phRe[i][j] = wRe + zRe; phIm[i][j] = wIm + zIm;
    wRe = alphaYRe*yRe-alphaYIm*yIm; wIm = alphaYRe*yIm+alphaYIm*yRe;
    zRe = betaYRe*xRe-betaYIm*xIm; zIm = betaYRe*xIm+betaYIm*xRe;
    phRe[i][j2] = wRe + zRe; phIm[i][j2] = wIm + zIm;

  }

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

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

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

}