/* variational method to seek eigen value and eigen function 1D */
/* coded by Ikeuchi Mitsuru     */
/* ver 0.0.1   2004.07.24       */

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

public class variationalMethod 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, sc1,sc2,sc3,sc4,sc5,sc6;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 30;

  int started = 1;
  double px0 = 12.0;
  double vWidth = 8.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/16.0;

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

  double cc[] = new double[100];
  double eigenEnergy[] = new double[20];

/* ----------------------------- 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,800,10,50,1210);
    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);

    sc1= new Scrollbar(Scrollbar.HORIZONTAL,200,10,0,210);
    sc1.addAdjustmentListener(this);

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

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

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

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

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

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,6,5,0));
    pnl.add(cpt);
    pnl.add(sc_ww);
    pnl.add(sc_wh);
    pnl.add(new Label("  "));
    pnl.add(new Label("  "));
    pnl.add(new Label("  "));
    add(pnl,"North");

    Panel pnl2 = new Panel();
    pnl2.setLayout(new GridLayout(1,6,5,0));
    pnl2.add(sc1);
    pnl2.add(sc2);
    pnl2.add(sc3);
    pnl2.add(sc4);
    pnl2.add(sc5);
    pnl2.add(sc6);
    add(pnl2,"South");

    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 = 0.01*(double)(sc_ww.getValue());
      setInitialCondition();
    } else if (ev.getSource() == sc_wh) { 
      vHight = -0.01*(double)(sc_wh.getValue());
      setInitialCondition();
    } else if (ev.getSource() == sc1) { 
      cc[1] = 0.01*(double)(sc1.getValue())-1.0;
    } else if (ev.getSource() == sc2) { 
      cc[2] = 0.01*(double)(sc2.getValue())-1.0;
    } else if (ev.getSource() == sc3) { 
      cc[3] = 0.01*(double)(sc3.getValue())-1.0;
    } else if (ev.getSource() == sc4) { 
      cc[4] = 0.01*(double)(sc4.getValue())-1.0;
    } else if (ev.getSource() == sc5) { 
      cc[5] = 0.01*(double)(sc5.getValue())-1.0;
    } else if (ev.getSource() == sc6) { 
      cc[6] = 0.01*(double)(sc6.getValue())-1.0;
    }
  }

  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,gby2,gb;
    double v1,v2,p1,p2,ve;

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

    gbx = 30;
    gby = 100;

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

    if (dispMode==0) {
      trialPh(eigenEnergy[1]);
      normalize();
      gOff.setColor(Color.getHSBColor( 0.01f, 0.9f, 0.9f ));
      gby2 = gby-(int)(100.0*eigenEnergy[1]);
      for (i=1; i<NNx-2; i++) {
        p1 = 800.0*ph[i]; p2 = 800.0*ph[i+1];
        gOff.drawLine(gbx+i,gby2-(int)(p1), gbx+i, gby2-(int)(p2));
      }
    }

    setFourier(ph);
    ve = kineticEnergy(ph);
    gOff.setColor(Color.getHSBColor( 0.66f, 0.5f, 0.9f ));
    gby2 = gby-(int)(100.0*ve);
    for (i=1; i<NNx-2; i++) {
      p1 = 800.0*ph[i]; p2 = 800.0*ph[i+1];
      gOff.drawLine(gbx+i,gby2-(int)(p1), gbx+i, gby2-(int)(p2));
    }

    gOff.setColor(Color.black);
    gOff.drawString("potential",10,40);
    gOff.drawString("width="+(int)(vWidth*1000.0)/1000.0+" au",630/6*1+10,40);
    gOff.drawString("depth="+(int)(vHight*1000.0)/1000.0+" au",630/6*2+10,40);

    gb = 630/6*3+10;
    gOff.setColor(Color.red);
    gOff.drawString("E1="+(float)eigenEnergy[1]+" ",gb,40);
    gOff.setColor(Color.blue);
    gOff.drawString("E=<u|H|u>/<u|u>="+(float)ve+" ",gb,60);

    gOff.drawString("u(x) = c1 sin(1 k x)  , k=pai/L ",gb,100);
    gOff.drawString("    +c2 sin(2 k x) ",gb,115);
    gOff.drawString("    +c3 sin(3 k x) ",gb,130);
    gOff.drawString("    +c4 sin(4 k x) ",gb,145);
    gOff.drawString("    +c5 sin(5 k x) ",gb,160);
    gOff.drawString("    +c6 sin(6 k x) ",gb,175);
	
    for (i=1; i<=6; i++) {
      gOff.drawString("c"+i+"="+(int)(cc[i]*100.0)/100.0+" ",630/6*(i-1)+10,300);
    }
  }

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

  double kineticEnergy(double ph[]) {
    int i;
    double s,hph,sn;

    s = 0.0;
    sn = 0.0;
    for (i=1; i<NNx-1; i++) {
      hph = (2.0*ph[i]-ph[i+1]-ph[i-1])/(2.0*dx*dx)+vv[i]*ph[i];
      s += ph[i]*hph*dx;
      sn += ph[i]*ph[i]*dx;
    }
    return ( s/(sn+1.0e-20) );
  }

  void setFourier(double ph[]) {
    int i,j;
    double x,pi;

    pi = 3.14159265;
    for (i=1; i<NNx-1; i++) {
      ph[i] = 0.0;
      x = i*dx;
      for (j=1; j<=6; j++) {
        ph[i] += cc[j]*Math.sin(j*pi*x/16.0);
      }
    }
    normalize();
  }

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

  void setInitialCondition() {
    cc[1] = 1.0;
    setPotential(potentialMode, vWidth, vHight);
  }

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

  void setPotential(int mode, double vw, double vh) {
    int i;
    double x,x0,a;

    x0 = 8.0;
    if (mode==1) {

      for (i=1;i<NNx-1;i++) {
        x = i*dx;
        if (x>x0-0.5*vw && x<(x0+0.5*vw)) {
          vv[i] = vh;
        } else {
          vv[i] = 0.0;
        }
      }

    } else if (mode==2) {

	  a = -vh*4.0/(vw*vw);
      for (i=1;i<NNx-1;i++) {
        x = i*dx;
        if (x>x0-0.5*vw && x<(x0+0.5*vw)) {
          vv[i] = a*(x-x0)*(x-x0)+vh;
        } else {
          vv[i] = 0.0;
        }
      }

    } else if (mode==3) {

      for (i=1;i<NNx-1;i++) {
        x = i*dx;
        if (x>x0-0.5*vw && x<=x0) {
          vv[i] = 2.0*((0.5*vw-(x0-x))/vw)*vh;
        } else if (x>x0 && x<=x0+0.5*vw) {
          vv[i] = 2.0*((0.5*vw-(x-x0))/vw)*vh;
        } else {
          vv[i] = 0.0;
        }
      }

    } else if (mode==4) {

      ;

    }

  }

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

  void timeEvolution() {
    int i;

    if (started==1) {

      getEigenValue(vHight);
      trialPh(eigenEnergy[1]);
      normalize();

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

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

  int getEigenValue(double energyLow) {
    int i;
    double er = 1.0e-14;

    for (i=0; i<10; i++) {
      eigenEnergy[i] = 0.0;
    }

    eigenEnergy[1] = seekEigenValue(energyLow,er);
    for (i=2; i<10; i++) {
      eigenEnergy[i] = seekEigenValue(eigenEnergy[i-1]+0.001,er);
      if (eigenEnergy[i]>=0) break;
    }
    return(i-1);
  }

  double seekEigenValue(double energyLow,double error) {
    double eh,el,em, de;
    double phLow,phHigh,phMid;

    de = 0.001;
    el = energyLow;
    phLow = trialPh(el);
    phHigh = phLow;
    while (el<0.0 && phLow*phHigh>0.0) {
      el = el + de;
      phHigh = trialPh(el);
    }
    if (el>=0.0){
      return(0.0);
    }

    eh = el;
    el = el - de;
    em = el;
    while (eh-el>error) {
	  em = (el+eh)/2.0;
      phMid = trialPh(em);
      if (phMid*phLow>0.0) {
        el = em;
      } else {
        eh = em;
      }
    }
    return( em );
  }

  double trialPh(double energy) {
    int i;
    double h2;

    ph[0] = 0.0;
    ph[1] = dx;
    h2 = dx*dx;
    for (i=1; i<NNx-1; i++) {
      ph[i+1] = 2.0*ph[i]-ph[i-1]+2.0*(vv[i]-energy)*ph[i]*h2;
    }
    return( ph[NNx-1] );
  }

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

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

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

}