/* Kronig-Penney Model                 */
/* periodic steepest descent method 1D */
/* coded by Ikeuchi Mitsuru            */
/* ver 0.0.1   2004.10.25              */

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

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

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

  Choice cpt,cvw;
  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 vWidth = 20.0;
  double vHight = 1.0;
  double energy = 0.5;
  int potentialMode = 1;
  int dispMode = 0;

  double t = 0.0;
  double dt = 1.0/256.0;
  double dx = 1.0/8.0;

  int NNx = 256;
  double ph[] = new double[NNx];
  double vv[] = new double[NNx];

  double sdEnergy[] = new double[20];
  double sdState[][] = new double[20][NNx];
  double wrk[] = new double[NNx];

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

  public void init() {

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

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


    cpt = new Choice();
    cpt.addItem("well"); cpt.addItem("parabolic"); cpt.addItem("triangle");
    cpt.addItemListener(this);

    cvw = new Choice();
    cvw.addItem("wave func"); cvw.addItem("density"); cvw.addItem("both");
    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);

    sc_ww= new Scrollbar(Scrollbar.HORIZONTAL,20,10,0,42);
    sc_ww.addAdjustmentListener(this);

    sc_wh= new Scrollbar(Scrollbar.HORIZONTAL,100,10,0,210);
    sc_wh.addAdjustmentListener(this);

    spx= new Scrollbar(Scrollbar.HORIZONTAL,500,10,0,2010);
    spx.addAdjustmentListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,6,5,0));
    pnl.add(new Label("  "));
    pnl.add(sc_ww);
    pnl.add(sc_wh);
    pnl.add(new Label("  "));
    pnl.add(new Label("  "));
    pnl.add(new Label("  "));
    add(pnl,"North");
	
    initSD(20);
    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() == cpt){
      potentialMode = 1+cpt.getSelectedIndex();
      setInitialCondition();
    } else if (ev.getSource() == cvw){
      dispMode = cvw.getSelectedIndex();
    }
  }

  public void actionPerformed(ActionEvent ev){
    if(ev.getSource() == bt_reset){ 
      started = -1;
    } 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) { 
      energy = 0.001*(double)(spx.getValue())-1.0;
    } else if (ev.getSource() == sc_ww) { 
      vWidth = (double)(sc_ww.getValue());
      setInitialCondition();
    } else if (ev.getSource() == sc_wh) { 
      vHight = 0.01*(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,istate,gbx,gby,gby2;
    double v1,v2,ph2,y,z,p1,p2;

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

    gbx = 30;
    gby = 280;

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

    for (istate=0;istate<20; istate++) {
      gOff.setColor(Color.getHSBColor( 0.66f, 0.6f, 0.9f ));
      gby2 = gby-(int)(100.0*sdEnergy[istate])-1;
      for (i=1; i<NNx-2; i++) {
        v1 = 40.0*sdState[istate][i]; v2 = 40.0*sdState[istate][i+1];
        gOff.drawLine(gbx+i,gby2-(int)(v1), gbx+i, gby2-(int)(v2));
      }
      gOff.setColor(Color.blue);
      gOff.drawString("E["+istate+"]="+(float)sdEnergy[istate]+"",400,80+12*istate);
    }


    gOff.setColor(Color.black);
    gOff.drawString("a/(a+b)="+(int)(vWidth*10.0)/10.0+"/32",630/6*1,40);
    gOff.drawString("hight="+(int)(vHight*1000.0)/1000.0+" au",630/6*2+20,40);
    gOff.drawString("V(x)=0(x<=a) else V(x)=Vh ",400,40);
	gOff.drawString("a+b=4au, L=32au ",400,60);
	
  }

/*------------------------- <i|j> integral ----------*/

  double intersection(double f[], double g[]) {
    int i;
    double s;

    s = 0.0;
    for (i=0; i<NNx; i++) {
      s += f[i]*g[i]*dx;
    }
    return( s );
  }


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

  void setInitialCondition() {
    setPotential();
  }

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

  void setPotential() {
    int i,j;
    int vw;

    vw = (int)vWidth;
    for (i=0;i<8;i++) {
      for (j=0;j<32;j++) {
	    if (j>=vw) {
          vv[i*32+j] = vHight;
        } else {
          vv[i*32+j] = 0.0;
        }
      }
    }
  }

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

  void timeEvolution() {
    int i;

    if (started==1) {

      SD(20, vv, 100, 0.01);
      sortState(20);

    } else if (started==-1) {
      t = 0.0;
      setInitialCondition();
      started = 0;
    }
  }

  void sortState(int maxState) {
    int ist,i;
    double w;

    for (ist=0; ist<maxState-1; ist++) {
      if (sdEnergy[ist]>sdEnergy[ist+1]+0.001) {
        for (i=0;i<NNx;i++) {
          w = sdState[ist][i];
          sdState[ist][i] = sdState[ist+1][i];
		  sdState[ist+1][i] = w;
        }
        w = sdEnergy[ist]; sdEnergy[ist]= sdEnergy[ist+1]; sdEnergy[ist+1]= w;
      }
    }
  }

/* ----------------- Gram-Schmidt --------- */

  void GramSchmidt(int stateMax) {
    int i,ist,istate;
    double s;

    normalize(sdState[0]);

    for (istate=1; istate<stateMax; istate++) {
      for (ist=0; ist<istate; ist++) {
        s = innerProduct(sdState[ist],sdState[istate]);
        for (i=0;i<NNx;i++) {
          sdState[istate][i] -= s*sdState[ist][i];
        }
      }
      normalize(sdState[istate]);
    }
  }

  double innerProduct(double f[],double g[]) {
    int i;
    double s;

    s = 0.0;
    for (i=0;i<NNx;i++) {
      s += f[i]*g[i]*dx;
    }
    return ( s );
  }

/* --------------------- steepest descent method -------------- */

  void SD(int stateMax, double v[], int iterMax, double dump) {
    int i,ist;

    for (i=0; i<iterMax; i++) {
	  for (ist=0; ist<stateMax; ist++) {
        sdEnergy[ist] = steepestDescent(sdState[ist], v, dump);
      }
	  GramSchmidt(stateMax);
    }
  }

  void initSD(int stateMax) {
    int i,ist;
    double pai = 3.14159265;

    for (ist=0; ist<stateMax; ist++) {
      for (i=0; i<NNx; i++) {
        sdState[ist][i] = Math.sin((ist+1.0)*i*dx*pai/16.0);
      }
	  normalize(sdState[ist]);
    }
  }

  double steepestDescent(double psi[], double v[], double dump) {
    int i;
    double h2,ei;

    h2 = 2.0*dx*dx;
    ei = energyOf(psi);
    for (i=1; i<NNx-1; i++) {
      wrk[i] = (2.0*psi[i]-psi[i+1]-psi[i-1])/h2 + (v[i]-ei)*psi[i];
    }
    wrk[0] = (2.0*psi[0]-psi[1]-psi[NNx-1])/h2 + (v[0]-ei)*psi[0];
    wrk[NNx-1] = (2.0*psi[NNx-1]-psi[0]-psi[NNx-2])/h2 + (v[NNx-1]-ei)*psi[NNx-1];
    for (i=0; i<NNx; i++) {
      psi[i] = psi[i]-dump*wrk[i];
    }
    normalize(psi);
    return (ei);
  }

  double energyOf(double psi[]) {
    int i;
    double h2,s;

    h2 = 2.0*dx*dx;
    s = 0.0;
    for (i=1; i<NNx-1; i++) {
      s += psi[i]*((2.0*psi[i]-psi[i+1]-psi[i-1])/h2 + vv[i]*psi[i])*dx;
    }
    s += psi[0]*((2.0*psi[0]-psi[1]-psi[NNx-1])/h2 + vv[0]*psi[0])*dx;
    s += psi[NNx-1]*((2.0*psi[NNx-1]-psi[0]-psi[NNx-2])/h2 + vv[NNx-1]*psi[NNx-1])*dx;
    return ( s );
  }

  void normalize(double ph[]) {
    int i;
    double s,a;

    s = 0.0;
    for (i=0; i<NNx; i++) {
      s += ph[i]*ph[i]*dx;
    }
    a = Math.sqrt(1.0/s);
    for (i=0; i<NNx; i++) {
      ph[i] = a*ph[i];
    }
  }

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

}