/* 4th-order laplacian - periodic steepest descent 2D */
/* coded by Ikeuchi Mitsuru     */
/* ver 0.0.1   2005.04.29       */

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

public class periodicO4SD2D extends Applet
  implements MouseListener, MouseMotionListener,
    ItemListener, ActionListener, AdjustmentListener, Runnable {

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

  Choice cvw;
  Button bt_reset, bt_start, bt_stop;
  Scrollbar sc_wr,sc_wd;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 30;

  int dispMode = 1;
  int objState = 1;
  int started = 1;


  double wellRadious = 2.0;
  double wellDepth = 1.0;


  int dgX, dgY, dgXb,dgYb;
  double viewTheta = -15.0*3.14/180.0;
  double viewFai = -72.0*3.14/180.0;
  double dtheta = 0.0*3.14/180.0;
  double pai = 3.1415926536;

  double cosTh = Math.cos(viewTheta); 
  double sinTh = Math.sin(viewTheta);
  double cosFi = Math.cos(viewFai);
  double sinFi = Math.sin(viewFai);

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

  int xpts[] = new int[NN];
  int ypts[] = new int[NN];
  double icol[] = new double[NN];

/* ----------------------------- 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("multi"); cvw.addItem("ph-0"); cvw.addItem("ph-1");
    cvw.addItem("ph-2"); cvw.addItem("ph-3"); cvw.addItem("ph-4");
    cvw.addItem("ph-5"); cvw.addItem("ph-6"); cvw.addItem("ph-7"); 
    cvw.addItem("ph-8"); cvw.addItem("ph-9"); cvw.addItem("V");
	cvw.addItemListener(this);
    cvw.select("ph-0");

    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_wr= new Scrollbar(Scrollbar.HORIZONTAL,200,10,0,410);
    sc_wr.addAdjustmentListener(this);

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

    addMouseListener(this);  
    addMouseMotionListener(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(cvw);
    pnl.add(sc_wr);
    pnl.add(sc_wd);
    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() == sc_wr) { 
      wellRadious = 0.01*(double)(sc_wr.getValue());
      setPotential();
    } else if (ev.getSource() == sc_wd) { 
      wellDepth = 0.01*(double)(sc_wd.getValue());
      setPotential();
    }
  }

  public void mousePressed(MouseEvent ev){ }
  public void mouseReleased(MouseEvent ev){
    dgXb = 0; dgYb = 0;
    dgX = 0; dgY = 0;
  }
  public void mouseClicked(MouseEvent ev){ }
  public void mouseEntered(MouseEvent ev){ } 
  public void mouseExited (MouseEvent ev){ }

  public void mouseMoved  (MouseEvent ev){ }

  public void mouseDragged(MouseEvent ev){
    dgXb = dgX; dgYb = dgY;
    dgX=ev.getX();
    dgY=ev.getY();
    if (dgXb==0 && dgYb==0) {
      dgXb = dgX; dgYb = dgY;
    }
    viewTheta += 0.5*3.14/180.0*(dgX-dgXb);
    viewFai += 0.5*3.14/180.0*(dgY-dgYb);
    cosTh = Math.cos(viewTheta); 
    sinTh = Math.sin(viewTheta);
    cosFi = Math.cos(viewFai);
    sinFi = Math.sin(viewFai);
  }


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

    if (dispMode==0) {
	
      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);
      }
	  
    } else if (dispMode>=1 && dispMode<11) {
	
	  gridPlot(sdState[dispMode-1],100.0,10.0);
	  for (i=0; i<10; i++) {
	    gOff.setColor(Color.blue);
		if (i==dispMode-1) {
		  gOff.setColor(Color.red);
		}
        gOff.drawString("E["+i+"]="+(float)(sdEnergy[i])+" ",400,i*20+80);
      }
	
    } else if (dispMode==11) {
	
      gridPlot(sdState[dispMode-1],0.0,10.0);
	  gOff.setColor(Color.blue);
	  for (i=0; i<10; i++) {
        gOff.drawString("E["+i+"]="+(float)(sdEnergy[i])+" ",400,i*20+80);
      }
	  
    } else if (dispMode==12) {
	  ;
    }

    gOff.setColor(Color.black);
    gOff.drawString("steepest descent method (periodic)",30,45);
    gOff.drawString("laplacian --> 4th-order approx.",30,65);

    gOff.setColor(Color.blue);
	gOff.drawString("view",630/6*3+10,40);
    gOff.drawString("radious="+(int)(wellRadious*100.0+0.5)/100.0+" ",630/6*4+10,40);
    gOff.drawString("depth="+(int)(wellDepth*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);
        }
      }
    }

  }

  void gridPlot(double ph[][],double phWeight,double vvWeight) {
    int i,j,gbx,gby;
    int cx,cy,cz;
    double z,px,py,pz,sc;
    double ph2;

    cx = NNx/2;
    cy = NNy/2;
    cz = 0;
    sc = 4.0;

    gbx = 60+(int)(cx/sc);
    gby = 60+(int)(cy/sc);

    for(j=0; j<NNy; j += 4) {
      for (i=0; i<NNx; i++) { 
        ph2 = phWeight*ph[i][j]*ph[i][j];
        z = ph2 + vvWeight*vv[i][j] + cz;
        py = cosFi*(j-cy)+sinFi*(z-cz) + cy; 
        pz = -sinFi*(j-cy)+cosFi*(z-cz) + cz;
		px = cosTh*(i-cx)+sinTh*(pz-cz) + cx; 
        pz = -sinTh*(i-cx)+cosTh*(pz-cz) + cz;
        xpts[i] = (int)(sc*px)+gbx;
        ypts[i] = (int)(sc*py)+gby;
		if (ph2>0.5) {
          icol[i] = 0.15;
        } else if (vv[i][j]>0.0) {
		  icol[i] = 0.33;
        } else {
          icol[i] = 0.66;
		}
	  }
      for (i=0; i<NNx-2; i++) {
		gOff.setColor(Color.getHSBColor((float)icol[i],0.9f,0.9f));
		gOff.drawLine(xpts[i],ypts[i],xpts[i+1],ypts[i+1]);
      }
	}

	for(i=0; i<NNx; i += 4) {
	  for (j=0; j<NNy; j++) {
		ph2 = phWeight*ph[i][j]*ph[i][j];
        z = ph2 + vvWeight*vv[i][j] + cz;
		py = cosFi*(j-cy)+sinFi*(z-cz) + cy; 
		pz = -sinFi*(j-cy)+cosFi*(z-cz) + cz;
		px = cosTh*(i-cx)+sinTh*(pz-cz) + cx; 
         pz = -sinTh*(i-cx)+cosTh*(pz-cz) + cz;
		xpts[j] = (int)(sc*px)+gbx;
		ypts[j] = (int)(sc*py)+gby;
		if (ph2>0.5) {
		icol[j] = 0.15;
		} else if (vv[i][j]>0.0) {
		icol[j] = 0.33;
		} else {
		icol[j] = 0.66;
		}
	  }
	  for (j=0; j<NNy-2; j++) { 
		gOff.setColor(Color.getHSBColor((float)icol[j],0.9f,0.9f));
		gOff.drawLine(xpts[j],ypts[j],xpts[j+1],ypts[j+1]);
	  }
	}

  }

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

  void setInitialCondition() {

    setPotential();
    initSD(10);
  }

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

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

    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;
	    r = Math.sqrt((x-x0)*(x-x0)+(y-y0)*(y-y0));
		if (r<=wellRadious) {
          vv[i][j] = -wellDepth;
		} else {
          vv[i][j] = 0.0;
		}
      }
    }
  }


/* ----------------------------- 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=0; i<NNx; i++) {
          for (j=0; j<NNy; 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=0; i<NNx; i++) {
      for (j=0; j<NNy; 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=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        s += ph[i][j]*ph[i][j]*dx*dy;
      }
    }
    a = Math.sqrt(1.0/s);
    for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; 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] = steepestDescentO4(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=0; i<NNx; i++) {
        for (j=0; j<NNx; j++) {
          sdState[ist][i][j] = Math.random()-0.5;
        }
      }
	  normalize(sdState[ist]);
    }
  }

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

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

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

    h2 = 2.0*dx*dx;
    s = 0.0; sn=0.0;
    for (i=0; i<NNx; i++) {
      ip = i+1; if (ip>=NNx) ip = 0;
      im = i-1; if (im<0) im = NNx-1;
      for (j=0; j<NNy; j++) {
        jp = j+1; if (jp>=NNy) jp = 0;
        jm = j-1; if (jm<0) jm = NNy-1;
        s += psi[i][j]*((4.0*psi[i][j]-psi[ip][j]-psi[im][j]-psi[i][jp]-psi[i][jm])/h2 + v[i][j]*psi[i][j])*dx*dy;
        sn += psi[i][j]*psi[i][j]*dx*dy;
      }
    }
    return ( s/sn );
  }

/* ----------------------------- 4th order laplacian -------------- */

  double steepestDescentO4(double psi[][], double v[][], double dump) {
    int i,ip,ip2,im,im2,j,jp,jp2,jm,jm2;
    double h2,ei;

    h2 = 12.0*2.0*dx*dx;
    ei = energyO4Of(psi);
    for (i=0; i<NNx; i++) {
	  ip = (i+1)%NNx; ip2=(i+2)%NNx; im = (i+NNx-1)%NNx; im2 = (i+NNx-2)%NNx;
      for (j=0; j<NNy; j++) {
        jp = (j+1)%NNy; jp2=(j+2)%NNy; jm = (j+NNy-1)%NNy; jm2 = (j+NNy-2)%NNy;
        wrk[i][j] = (60.0*psi[i][j]-16.0*psi[ip][j]-16.0*psi[im][j]-16.0*psi[i][jp]-16.0*psi[i][jm]+psi[ip2][j]+psi[im2][j]+psi[i][jp2]+psi[i][jm2])/h2 + (vv[i][j]-ei)*psi[i][j];
      }
    }
    for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        psi[i][j] = psi[i][j]-dump*wrk[i][j];
      }
    }
    normalize(psi);
    return (ei);
  }

  double energyO4Of(double psi[][]) {
    int i,ip,ip2,im,im2,j,jp,jp2,jm,jm2;
    double h2,dxy,s;

    h2 = 12.0*2.0*dx*dx;
	dxy = dx*dy;
    s = 0.0;
    for (i=0; i<NNx; i++) {
	  ip = (i+1)%NNx; ip2=(i+2)%NNx; im = (i+NNx-1)%NNx; im2 = (i+NNx-2)%NNx;
      for (j=0; j<NNy; j++) {
        jp = (j+1)%NNy; jp2=(j+2)%NNy; jm = (j+NNy-1)%NNy; jm2 = (j+NNy-2)%NNy;
        s += psi[i][j]*((60.0*psi[i][j]-16.0*psi[ip][j]-16.0*psi[im][j]-16.0*psi[i][jp]-16.0*psi[i][jm]+psi[ip2][j]+psi[im2][j]+psi[i][jp2]+psi[i][jm2])/h2 + vv[i][j]*psi[i][j])*dxy;
      }
    }
    return ( s );
  }

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

}