/* tunnel effect - QED1D       */
/* coded by Ikeuchi Mitsuru    */
/* ver 0.0.1   2006.01.09      */

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

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

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

  Button bt_reset, bt_startStop;
  Scrollbar sc_wwidth, sc_whight, sc_spx;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 30;

  int resetSW = 0;
  int started = 1;

  double px0 = 0.55;
  double wallPos = 120.0;
  double wallWidth = 10.0;
  double wallHight = 0.165;

  double t = 0.0;
  double dx = 1.0/2.0;

  double dt = 0.5*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 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,400);
    setBackground(Color.white);

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

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

    bt_startStop= new Button("start/stop");
    bt_startStop.addActionListener(this); 

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

    sc_whight= new Scrollbar(Scrollbar.HORIZONTAL,165,10,0,510);
    sc_whight.addAdjustmentListener(this);

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

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,6,5,0));
/*
    pnl.add(new Label("  "));
*/
    pnl.add(bt_reset);
    pnl.add(bt_startStop);
    pnl.add(new Label("  "));
    pnl.add(sc_whight);
    pnl.add(sc_wwidth);
    pnl.add(sc_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;
	  resetSW = 1;
      setInitialCondition();
      started = 0;
    } else if (ev.getSource() == bt_startStop){
      if (started==0) {
        started = 1;
      } else {
        started = 0;
      }
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_spx) { 
      px0 = 0.01*(double)(sc_spx.getValue());
      setInitialCondition();
    } else if (ev.getSource() == sc_wwidth) { 
      wallWidth = 1.0*(double)(sc_wwidth.getValue());
      setInitialCondition();
    } else if (ev.getSource() == sc_whight) { 
      wallHight = 0.001*(double)(sc_whight.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;

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

    gbx = 30;
    gby = 240;

    plotPotential(gbx, gby, 400.0);
    plotProbability(gbx, gby, 2000.0);
    plotPhaseLine(gbx, gby, 200.0);
	//plotPhase(gbx, gby, 200.0);

    gOff.setColor(Color.black);
    gOff.drawString("0 ",gbx-5,gby+20);
    gOff.drawString(" "+NNx*dx+" au",gbx+NNx-5,gby+20);
    gOff.drawString("(="+(int)(NNx*dx*52.918/10.0+0.5)/100.0+" nm)",gbx+NNx-5,gby+40);

    gOff.setColor(Color.blue);
    gOff.drawString("t="+(int)(t*24.189/1000.0*100+0.5)/100.0+" fs",10,40);

    if (started==1) {
      gOff.drawString("started",630/6*1+10,40);
    } else {
      gOff.drawString("stopped",630/6*1+10,40);
    }

    gOff.drawString("hight="+(int)(wallHight*1000.0+0.5)/1000.0+" ",630/6*3+10,40);
    gOff.drawString("(="+(int)(wallHight*27.2*10.0+0.5)/10.0+" eV)",630/6*3+10,60);

    gOff.drawString("width="+(int)(wallWidth+0.5)+" ",630/6*4+10,40);
    gOff.drawString("(="+(int)(wallWidth*52.918/10.0+0.5)/100.0+" nm)",630/6*4+10,60);
	
    gOff.drawString("px0="+(int)(px0*1000.0+0.5)/1000.0+" ",630/6*5+10,40);
    gOff.drawString("K="+(int)(kineticEnergy()*1000.0+0.5)/1000.0+" au",630/6*5+10,60);
    gOff.drawString("(="+(int)(kineticEnergy()*27.2*100.0+0.5)/100.0+" eV)",630/6*5+10,80);
    gOff.drawString("V="+(int)(potentialEnergy()*1000.0+0.5)/1000.0+" au",630/6*5+10,100);
    gOff.drawString("E="+(int)(energy()*1000.0+0.5)/1000.0+" au",630/6*5+10,120);
    gOff.drawString("(="+(int)(energy()*27.2*100.0+0.5)/100.0+" eV)",630/6*5+10,140);

    gOff.setColor(Color.blue);
    gOff.drawString("nLeft="+(int)(probabilityLeft()*10000.0+0.5)/100.0+" %",60,360);
    gOff.drawString("nRight="+(int)(probabilityRight()*10000.0+0.5)/100.0+" %",260,360);
    gOff.drawString("norm="+norm()+" ",60,380);

  }


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

  void plotProbability(int gbx, int gby, double mag) {
    int i;
    double ph2;

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

  void plotPotential(int gbx, int gby, double mag) {
    int i;
    double v1,v2;

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

  void plotPhaseLine(int gbx, int gby, double mag) {
    int i;
    double y,z;

    for (i=0; i<NNx; i++) { 
      y = mag*phRe[i];
      z = mag*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]);
    }
  }

  void plotPhase(int gbx, int gby, double mag) {
    int i;
    double y,z;

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

/*-------------------------- 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(50.0, 10.0, px0);
    setPotential(wallPos, wallHight, wallWidth);
  }

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

  void setPotential(double vPos, double vHight, double vWidth) {
    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 (resetSW==1) {
      setInitialCondition();
      resetSW = 0;
    }

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

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

  void kxStep() {
    int i;
    double bbRe,bbIm,aaRe,aaIm,aaAb,auAb;

    bbRe = 2.0;
    bbIm = 4*dx*dx/dt;
    aaRe = -2.0;
    aaIm = 4*dx*dx/dt;
    aaAb = aaRe*aaRe+aaIm*aaIm;

    for (i=1; i<NNx-1; i++) {
      bRe[i] = bbRe*phRe[i]-bbIm*phIm[i] - phRe[i+1] - phRe[i-1];
      bIm[i] = bbRe*phIm[i]+bbIm*phRe[i] - phIm[i+1] - phIm[i-1];
    }

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

}