/* harmonic oscillator - steepest descent 2D */
/* coded by Ikeuchi Mitsuru   */
/* ver 0.0.1   2004.08.06     */

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

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

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

  Choice cvw;
  Button bt_reset, bt_start, bt_stop;
  Scrollbar spx,spy;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 30;

  int dispMode = 3;
  int started = 1;


  double kx0 = 2.0;
  double ky0 = 2.0;


  double t = 0.0;

  double dx = 8.0/64.0;
  double dy = 8.0/64.0;

  double dt = 2.0*dx*dx;

  int NNx = 64+1;
  int NNy = 64+1;
  int NN = NNx+1; /* max(NNx,NNy)+1 */

  double vv[][] = new double[NNx][NNy];

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

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

  public void init() {

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

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

    cvw = new Choice();
    cvw.addItem("density"); cvw.addItem("psi");
    cvw.addItemListener(this);
    cvw.select("density");

    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);
    
    spx= new Scrollbar(Scrollbar.HORIZONTAL,200,10,50,210);
    spx.addAdjustmentListener(this);

    spy= new Scrollbar(Scrollbar.HORIZONTAL,200,10,50,210);
    spy.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(new Label("  "));
    pnl.add(spx);
    pnl.add(spy);
    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 itemStateChanged(ItemEvent ev){
    if (ev.getSource() == cvw){
      dispMode = cvw.getSelectedIndex();
    }
  }

  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) { 
      kx0 = 0.01*(double)(spx.getValue());
      setPotential();
    } else if (ev.getSource() == spy) { 
      ky0 = 0.01*(double)(spy.getValue());
      setPotential();
    }
  }


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


    for (i=0; i<10; i++) {
      gbx = 10+120*(i%5);
      gby = 50+130*(i/5);
      densPlot2D50(sdState[i], gbx, gby);
	  gOff.drawString("E["+i+"]="+(float)(sdEnergy[i])+" ",gbx+10,gby+130);
    }


    gOff.setColor(Color.black);
    gOff.drawString("steepest descent method",30,45);

    gOff.setColor(Color.blue);
    gOff.drawString("kx="+(int)(kx0*100.0+0.5)/100.0+" ",630/6*4+10,40);
    gOff.drawString("ky="+(int)(ky0*100.0+0.5)/100.0+" ",630/6*5+10,40);
  }

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


  void densPlot2D50(double psi[][], int gbx, int gby) {
    int i,j;
    double ph2;

    for (i=7; i<57; i++) {
      for (j=7; j<57; j++) {
        ph2 = 2.0*(psi[i][j]*psi[i][j])+0.0004*vv[i][j];
        if (ph2>0.9) {
          ph2 = 0.9;
        }
        if (ph2>=0.01) {
          gOff.setColor(Color.getHSBColor(
            (float)(0.15+0.0004*vv[i][j]), 0.95f, (float)(ph2+0.05)));
          gOff.fillRect(gbx+i*2,gby+j*2, 2, 2);
        }
      }
    }

  }


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

  void setInitialCondition() {

    setPotential();
    initSD(10);
  }

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

  void setPotential() {
    int i,j;
    double x,y,x0,y0;

    x0 = 4.0; y0 = 4.0;
    for (i=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        x = i*dx;
        y = j*dy;
        vv[i][j] = kx0*(x-x0)*(x-x0)+ky0*(y-y0)*(y-y0);
      }
    }
  }


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

  void timeEvolution() {
    int i;

    if (started==1) {

      SD(10, vv, 10, 0.005);
    }
  }

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

  void GramSchmidt(int stateMax) {
    int i,j,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=1; i<NNx-1; i++) {
          for (j=1; j<NNy-1; j++) {
            sdState[istate][i][j] -= s*sdState[ist][i][j];
          }
        }
      }
      normalize(sdState[istate]);
    }
  }

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

    s = 0.0;
    for (i=1; i<NNx-1; i++) {
      for (j=1; j<NNy-1; j++) {
        s += f[i][j]*g[i][j]*dx*dy;
      }
    }
    return ( s );
  }

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

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

/* --------------------- 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,j,ist;
    double pai = 3.14159265;

    for (ist=0; ist<stateMax; ist++) {
      for (i=1; i<NNx-1; i++) {
        for (j=1; j<NNx-1; j++) {
          sdState[ist][i][j] = Math.random()-0.5;
        }
      }
	  normalize(sdState[ist]);
    }
  }

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

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

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

    h2 = 2.0*dx*dx;
    s = 0.0; sn=0.0;
    for (i=1; i<NNx-1; i++) {
      for (j=1; j<NNx-1; j++) {
        s += psi[i][j]*((4.0*psi[i][j]-psi[i+1][j]-psi[i-1][j]-psi[i][j-1]-psi[i][j+1])/h2 + v[i][j]*psi[i][j])*dx*dy;
        sn += psi[i][j]*psi[i][j]*dx*dy;
      }
    }
    return ( s/sn );
  }

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

}