/* quantum well QED1D         */
/* coded by Ikeuchi Mitsuru   */
/* ver 0.0.1   2004.04.29     */

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

public class quantumWellQED1D extends Applet
  implements ActionListener, ItemListener, AdjustmentListener, Runnable {

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

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

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 30;

  int started = 1;
  int quantumNumber = 8;
  double wellPos = 4.0;
  double wellWidth = 8.0;
  double wellHight = 100.0;

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

  double dt = 0.25*dx*dx;


  int NNx = 32*16+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);

    cqn = new Choice();
    cqn.addItem("1"); cqn.addItem("2"); cqn.addItem("3");
    cqn.addItem("4"); cqn.addItem("5"); cqn.addItem("6");
    cqn.addItem("7"); cqn.addItem("8"); cqn.addItem("9"); 
    cqn.addItem("10"); cqn.addItem("11"); cqn.addItem("12");
    cqn.addItem("13"); cqn.addItem("14"); cqn.addItem("15");
    cqn.addItem("16"); 
    cqn.addItemListener(this);
    cqn.select("8");

    sc_ww= new Scrollbar(Scrollbar.HORIZONTAL,80,10,20,110);
    sc_ww.addAdjustmentListener(this);

    sc_wh= new Scrollbar(Scrollbar.HORIZONTAL,100,10,0,310);
    sc_wh.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(cqn);
    pnl.add(sc_ww);
    pnl.add(sc_wh);
    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 itemStateChanged(ItemEvent ev){
    if (ev.getSource() == cqn) {
      quantumNumber = 1+cqn.getSelectedIndex();
      t = 0.0;
      setInitialCondition();
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_wh) { 
      wellHight = 1.0*(double)(sc_wh.getValue());
      setInitialCondition();
    } else if (ev.getSource() == sc_ww) { 
      wellWidth = 0.1*(double)(sc_ww.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 = 50;
    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*200.0>=1.0) {
	    gOff.drawLine(gbx+i,gby-(int)(200.0*ph2), gbx+i, gby);
      }
    }

    for (i=0; i<NNx; i++) { 
      y = 100.0*phRe[i];
      z = 100.0*phIm[i];
      xpts[i] = gbx+i-(int)(0.5*z);
      ypts[i] = gby+(int)(0.5*z)-(int)(y);
    }
    for (i=0; i<NNx-2; i++) {
      gOff.setColor(Color.getHSBColor(0.15f,0.9f,(float)(0.5+0.3*phIm[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.drawString("K="+(int)(kineticEnergy()*100.0+0.5)/100.0+" au",630/6*0+10,60);
    gOff.drawString("V="+(int)(potentialEnergy()*100.0+0.5)/100.0+" au",630/6*1+10,60);
    gOff.drawString("E="+(int)(energy()*100.0+0.5)/100.0+" au",630/6*2+10,60);
    gOff.drawString("E(theory)="+theoreticalEnergy(wellWidth, quantumNumber)+" au",630/6*1+10,80);
	
    gOff.setColor(Color.black);
    gOff.drawString("quantum No.",630/6*3+10,40);
    gOff.drawString("width="+(int)(wellWidth*52.918+0.5)+" pm",630/6*4,40);
    gOff.drawString("hight="+(int)(wellHight*10+0.5)/10.0+" au",630/6*5,40);
    gOff.drawString("Lx="+(int)(NNx*dx*52.918+0.5)+" pm",630/6*5,60);

  }


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

  double theoreticalEnergy(double len, double qNumber) {

    return ( (2.0*Math.PI)*(2.0*Math.PI)/(8.0*len*len)*qNumber*qNumber );
  }

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

  void setInitialCondition() {
    t = 0.0;
    setWave(quantumNumber);
    setPotential();
  }

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

  void setPotential() {
    int i;
    double x;

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

  }

/* --------------------------- set wave function -------------- */

  void setWave(int n) {
    int i;
    double a,x,phAb,phPh;

    for (i=1;i<NNx-1;i++) {
      x = i*dx;
      if (x<wellPos) {
        phRe[i] = 0.0;
        phIm[i] = 0.0;
      } else if (x>=wellPos && x<=wellPos+wellWidth) {
        phAb = Math.sin(n*Math.PI*(x-wellPos)/wellWidth);
        phPh = 0.0;
		phRe[i] = phAb*Math.cos(phPh);
        phIm[i] = phAb*Math.sin(phPh);
      } else {
        phRe[i] = 0.0;
        phIm[i] = 0.0;
      }
    }
    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<50; i++) {
        t = t + dt;
        phaseStep();
        kxStep();
        phaseStep();
      }
    }
  }

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

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

    r = 4.0*dx*dx/dt;
    for (i=1; i<NNx-1; i++) {
      bbRe[i] = 2.0;
      bbIm[i] = r;
      aaRe[i] = -2.0;
      aaIm[i] = r;
      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 -------------- */

}