/* phase of wave function QED1D */
/* coded by Ikeuchi Mitsuru     */
/* ver 0.0.1   2004.04.25       */
/* ver 0.0.2   2005.11.23  bug in energy() fixed */

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

public class phaseQED1D extends Applet
  implements ActionListener, AdjustmentListener, Runnable {

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

  Button bt_reset, bt_start, bt_stop;
  Scrollbar sc_ww, sc_wh, spx;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 30;

  int started = 1;
  double px0 = 12.0;
  double wallPos = 7.0;
  double wallWidth = 0.3;
  double wallHight = 85.0;


  double t = 0.0;
  double dx = 16.0/512.0;

  double dt = 0.25*dx*dx;


  int NNx = 512+1;
  double phRe[] = new double[NNx];
  double phIm[] = new double[NNx];  

  double vv[] = new double[NNx];

  int NN = NNx+1;
  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];

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

  public void init() {

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

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

    bt_reset= new Button("reset");
    bt_reset.addActionListener(this);

    bt_start= new Button("start");
    bt_start.addActionListener(this); 

    bt_stop = new Button("stop");
    bt_stop.addActionListener(this);

    sc_ww= new Scrollbar(Scrollbar.HORIZONTAL,30,10,0,110);
    sc_ww.addAdjustmentListener(this);

    sc_wh= new Scrollbar(Scrollbar.HORIZONTAL,85,10,0,160);
    sc_wh.addAdjustmentListener(this);

    spx= new Scrollbar(Scrollbar.HORIZONTAL,120,10,0,210);
    spx.addAdjustmentListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,6,5,0));
    pnl.add(bt_reset);
    pnl.add(bt_start);
    pnl.add(bt_stop);
    pnl.add(sc_ww);
    pnl.add(sc_wh);
    pnl.add(spx);
    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 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());
      setInitialCondition();
    } else if (ev.getSource() == sc_ww) { 
      wallWidth = 0.01*(double)(sc_ww.getValue());
      setInitialCondition();
    } else if (ev.getSource() == sc_wh) { 
      wallHight = 1.0*(double)(sc_wh.getValue());
      setInitialCondition();
    }
  }

  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 i,gbx,gby;
    double v1,v2,ph2,y,z;

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

    gbx = 30;
    gby = 220;

    gOff.setColor(Color.green);
    for (i=1; i<NNx-2; i++) {
      v1 = vv[i]; v2 = vv[i+1];
	  gOff.drawLine(gbx+i,gby-(int)(v1), gbx+i, gby-(int)(v2));
    }

    gOff.setColor(Color.cyan);
	for (i=1; i<NNx-1; i++) {
      ph2 = phRe[i]*phRe[i]+phIm[i]*phIm[i];
      if (ph2*100.0>=1.0) {
	    gOff.drawLine(gbx+i,gby-(int)(100.0*ph2), gbx+i, gby);
      }
    }

    for (i=0; i<NNx; i++) { 
      y = 30.0*phRe[i];
      z = 30.0*phIm[i];
      xpts[i] = gbx+i-(int)(0.5*z);
      ypts[i] = gby+(int)(0.5*z)-(int)(y);
    }
    gOff.setColor(Color.getHSBColor(0.15f,0.9f,0.6f));
    for (i=0; i<NNx-2; i++) {
      gOff.drawLine(xpts[i],ypts[i],xpts[i+1],ypts[i+1]);
    }

    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.setColor(Color.black);
    gOff.drawString("Vwidth="+(int)(wallWidth*52.918+0.5)+" pm",630/6*3+10,40);
    gOff.drawString("Vhight="+(int)(wallHight+0.5)+" ",630/6*4+10,40);
    gOff.drawString("px0="+(int)(px0*10.0+0.5)/10.0+" au",630/6*5+10,40);
    gOff.drawString("vx0="+(int)(px0*2.79+0.5)/10.0+" %c",630/6*5+10,60);
    gOff.drawString("Lx="+(int)(NNx*dx*52.918+0.5)+" pm",630/6*5+10,80);

    gOff.setColor(Color.blue);
    gOff.drawString("E="+(int)(energy()*100.0+0.5)/100.0+" au",630/6*5+10,100);

    gOff.drawString("nLeft="+(int)(probabilityLeft()*10000.0+0.5)/100.0+" %",630/6*1+10,80);
    gOff.drawString("nRight="+(int)(probabilityRight()*10000.0+0.5)/100.0+" %",630/6*3+10,80);

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

  double norm() {
    int i;
    double p;

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

  double probabilityLeft() {
    int i;
    double x,p;

    p = 0.0;
    for (i=1;i<NNx-1;i++) {
      x = i*dx;
      if (x<wallPos) {
        p += (phRe[i]*phRe[i]+phIm[i]*phIm[i])*dx;
      }
    }
    return ( p );
  }

  double probabilityRight() {
    int i;
    double x,p;

    p = 0.0;
    for (i=1;i<NNx-1;i++) {
      x = i*dx;
      if (x>wallPos+wallWidth) {
        p += (phRe[i]*phRe[i]+phIm[i]*phIm[i])*dx;
      }
    }
    return ( p );
  }

  double energy() {
    int i;
    double s,hphRe,hphIm;

    s = 0.0;
    for (i=1; i<NNx-1; i++) {
      hphRe = (2.0*phRe[i]-phRe[i+1]-phRe[i-1])/(2.0*dx*dx)+vv[i]*phRe[i];
      hphIm = (2.0*phIm[i]-phIm[i+1]-phIm[i-1])/(2.0*dx*dx)+vv[i]*phIm[i];
      s += (phRe[i]*hphRe + phIm[i]*hphIm)*dx;
    }
    return ( s );
  }

  double kineticEnergy() {
    int i;
    double s,hphRe,hphIm;

    s = 0.0;
    for (i=1; i<NNx-1; i++) {
      hphRe = (2.0*phRe[i]-phRe[i+1]-phRe[i-1])/(2.0*dx*dx);
      hphIm = (2.0*phIm[i]-phIm[i+1]-phIm[i-1])/(2.0*dx*dx);
      s += (phRe[i]*hphRe + phIm[i]*hphIm)*dx;
    }
    return ( s );
  }

  double potentialEnergy() {
    int i;
    double s,hphRe,hphIm;

    s = 0.0;
    for (i=1; i<NNx-1; i++) {
      s += vv[i]*(phRe[i]*phRe[i] + phIm[i]*phIm[i])*dx;
    }
    return ( s );
  }

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

  void setInitialCondition() {
    t = 0.0;
    setWave(3.0, 0.5, px0);
    setPotential(wallPos, wallWidth, wallHight);
  }

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

  void setPotential(double vPos, double vWidth, double vHight) {
    int i;
    double x;

    for (i=1;i<NNx-1;i++) {
      x = i*dx;
      if (x>vPos && x<(vPos+vWidth)) {
        vv[i] = vHight;
      } else {
        vv[i] = 0.0;
      }
    }

  }

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

  void setWave(double wavePos, double waveWidth, double momentum) {
    int i;
    double a,x,phAb,phPh;

    a = Math.pow(2.0*Math.PI*waveWidth*waveWidth,-0.25);
    for (i=1;i<NNx-1;i++) {
      x = i*dx;
      phAb = a*Math.exp(-((x-wavePos)/(2.0*waveWidth))*((x-wavePos)/(2.0*waveWidth)));
      phPh = momentum*x;
      phRe[i] = phAb*Math.cos(phPh);
      phIm[i] = phAb*Math.sin(phPh);
    }
    phRe[0] = 0.0;
    phIm[0] = 0.0;
    phRe[NNx-1] = 0.0;
    phIm[NNx-1] = 0.0;

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

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

  void timeEvolution() {
    int i;

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

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

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

    for (i=1; i<NNx-1; i++) {
      bbRe[i] = 2.0;
      bbIm[i] = 4*dx*dx/dt;
      aaRe[i] = -2.0;
      aaIm[i] = 4*dx*dx/dt;
      bRe[i] = bbRe[i]*phRe[i]-bbIm[i]*phIm[i] - phRe[i+1] - phRe[i-1];
      bIm[i] = bbRe[i]*phIm[i]+bbIm[i]*phRe[i] - phIm[i+1] - phIm[i-1];
    }

    aaAb = aaRe[1]*aaRe[1]+aaIm[1]*aaIm[1];
	uRe[1] = aaRe[1]/aaAb;
	uIm[1] = -aaIm[1]/aaAb;
    phRe[1] = bRe[1]*uRe[1] - bIm[1]*uIm[1];
    phIm[1] = 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] = (bRe[i]-phRe[i-1])*uRe[i] - (bIm[i]-phIm[i-1])*uIm[i];
      phIm[i] = (bRe[i]-phRe[i-1])*uIm[i] + (bIm[i]-phIm[i-1])*uRe[i];
    }

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

  }

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

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

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

}