/* steepest descent - time dependent Kohn-Sham 2D  */
/* coded by Ikeuchi Mitsuru     */
/* ver 0.0.1   2004.11.20       */

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

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

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

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

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 30;

  int dispMode = 0;
  int dispObject = 0;
  int started = 1;
  int lossFlag = 0;

  double px1 = 2.0;
  double py1 = 0.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 = 1.0/4.0;
  double dy = 1.0/4.0;

  double dt = 1.0*dx*dx;
  double lz = 8.0;

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

  int NNe = 2;

  int spin[] = { 1,1,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0 };
  double energy[] = new double[NNe];
  double psi[][][][] = new double[NNe][NNx][NNy][2];
  double wrk[][][] = new double[NNx][NNy][2];

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

  double vvext[][] = new double[NNx][NNy];
  double vvh[][] = new double[NNx][NNy];
  double vvx[][] = new double[NNx][NNy];
  double vvc[][] = new double[NNx][NNy];
  
  double rho[][] = new double[NNx][NNy];

  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];
  double icol[] = new double[NN];

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

  public void init() {

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

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

    che = new Choice();
    che.addItem("0"); che.addItem("1");
    che.addItemListener(this);

    cvw = new Choice();
    cvw.addItem("rho"); cvw.addItem("grid-i");
    cvw.addItem("dens-i"); cvw.addItem("phase-i");
    cvw.addItem("Vext"); cvw.addItem("VH x100");
    cvw.addItem("Vx x100"); cvw.addItem("Vc x100");
    cvw.addItem("Veff");
    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);

    bt_loss = new Button("loss");
    bt_loss.addActionListener(this);
    
    spx= new Scrollbar(Scrollbar.HORIZONTAL,70,10,0,110);
    spx.addAdjustmentListener(this);

    spy= new Scrollbar(Scrollbar.HORIZONTAL,50,10,0,110);
    spy.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(che);
    pnl.add(cvw);
    pnl.add(bt_loss);
    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() == che){
      dispObject = che.getSelectedIndex();
    } 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;
    } else if(ev.getSource() == bt_loss){
	  if (lossFlag==0) {
	    lossFlag = 1;
	  } else {
	    lossFlag = 0;
	  }
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == spx) { 
      px1 = 0.1*(double)(spx.getValue())-5.0;
    } else if (ev.getSource() == spy) { 
      py1 = 0.1*(double)(spy.getValue())-5.0;
    }
  }

  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 ie, gb, gbx,gby;
    double ev0,ek0,ev1,ek1;

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

    ev0 = meanPotential(psi[0]);
    ek0 = meanKinetic(psi[0]);

    ev1 = meanPotential(psi[1]);
    ek1 = meanKinetic(psi[1]);

    if (dispMode==0) {

      gridPlot(rho,800.0, 0.2,0);

    } else if (dispMode==1) {

      gridPsiPlot(dispObject,100.0, 0.2);

    } else if (dispMode==2) {

      densPlot(psi[dispObject],20.0, 0);

    } else if (dispMode==3) {
	
      densPlot(psi[dispObject],20.0, 1);

    } else if (dispMode==4) {
	
      gridPlot(vvext,0.2, 0.0,1);

    } else if (dispMode==5) {
	
      gridPlot(vvh,20.0, 0.0,1);

    } else if (dispMode==6) {
	
      gridPlot(vvx,20.0, 0.0,1);

    } else if (dispMode==7) {
	
      gridPlot(vvc,20.0, 0.0,1);

    } else if (dispMode==8) {
	
      gridPlot(vv,0.2, 0.0,1);

    }



    gOff.setColor(Color.black);
    gOff.drawString("object",630/6*3+10,40);
    gOff.drawString("view",630/6*4+10,40);
	if (lossFlag==1) {
      gOff.drawString("loss=ON",630/6*5+10,40);
    } else {
      gOff.drawString("loss=OFF",630/6*5+10,40);
    }
    gOff.setColor(Color.blue);
    gOff.drawString("t="+(int)(t*2418.9+0.5)/100.0+" as",10,40);

    gb = 330;
	gOff.setColor(Color.blue);
    gOff.drawString("Lx="+(int)(NNx*dx+0.5)+", Ly="+(int)(NNy*dy+0.5)+" au",gb,60);

    for (ie=0;ie<NNe;ie++) {
	  gbx = 330+150*(ie/2); gby = 100+80*(ie%2);
	  gOff.drawString("E"+ie+"="+(float)(energy[ie])+" ",gb,gby);
      gOff.drawString("<x"+ie+">="+(float)(meanPosX(psi[ie]))+" ",gb,gby+20);
      gOff.drawString("<y"+ie+">="+(float)(meanPosY(psi[ie]))+" ",gb,gby+40);
	}

	for (ie=0; ie<NNe; ie++) {
	  gOff.drawString("|<"+ie+"|"+dispObject+">|="+(float)(innerNorm(psi[ie],psi[dispObject]))+" ",gb,260+15*ie);
	}
  }

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

  void gridPsiPlot(int ie, 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 = 30+(int)(cx/sc);
    gby = 50+(int)(cy/sc);

    for(j=0; j<NNy; j += 2) {
      for (i=0; i<NNx; i++) { 
        ph2 = psi[ie][i][j][0]*psi[ie][i][j][0]+psi[ie][i][j][1]*psi[ie][i][j][1];
        z = phWeight*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.002) {
          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-1; 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 += 2) {
      for (j=0; j<NNy; j++) {
        ph2 = psi[ie][i][j][0]*psi[ie][i][j][0]+psi[ie][i][j][1]*psi[ie][i][j][1];
        z = phWeight*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.002) {
          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-1; j++) { 
        gOff.setColor(Color.getHSBColor((float)icol[j],0.9f,0.9f));
        gOff.drawLine(xpts[j],ypts[j],xpts[j+1],ypts[j+1]);
      }
    }
  }

  void densPlot(double ph[][][], double weight, int mode) {
    int i,j;
    int gbx,gby,mag;
    double ph2,col;

    gbx = 30;
    gby = 70;
    mag = 3;
    gOff.setColor(Color.black);
    gOff.fillRect(gbx-1,gby-1,mag*NNx+1,mag*NNy+1);

    for (i=1; i<NNx-1; i++) {
      for (j=1; j<NNy-1; j++) {

        ph2 = weight*(ph[i][j][0]*ph[i][j][0]+ph[i][j][1]*ph[i][j][1]);
        if (ph2>0.01) {
          if (ph2>1.0) {
            ph2 = 1.0;
          }
          if (mode==0) {
            gOff.setColor(Color.getHSBColor(0.15f, 0.95f, (float)(ph2)));
          } else if (mode==1) {
            gOff.setColor(Color.getHSBColor((float)(0.5+0.5*ph[i][j][0]/Math.sqrt(ph2/weight)), 0.9f, (float)(ph2)));
          }
          gOff.fillRect(gbx+mag*i,gby+mag*j, mag, mag);
        }
      }
    }
  }

  void gridPlot(double fn[][],double phWeight, double vvWeight,int mode) {
    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 = 30+(int)(cx/sc);
    gby = 50+(int)(cy/sc);

    for(j=0; j<NNy; j += 2) {
      for (i=0; i<NNx; i++) { 
        ph2 = fn[i][j];
        z = phWeight*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.001) {
          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-1; i++) {
        if (mode==0) {
          gOff.setColor(Color.getHSBColor((float)icol[i],0.9f,0.9f));
        } else {
          gOff.setColor(Color.getHSBColor(0.66f,0.9f,0.9f));
        }
        gOff.drawLine(xpts[i],ypts[i],xpts[i+1],ypts[i+1]);
      }
    }

    for(i=0; i<NNx; i += 2) {
      for (j=0; j<NNy; j++) {
        ph2 = fn[i][j];
        z = phWeight*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.001) {
          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-1; j++) { 
        if (mode==0) {
          gOff.setColor(Color.getHSBColor((float)icol[j],0.9f,0.9f));
        } else {
          gOff.setColor(Color.getHSBColor(0.66f,0.9f,0.9f));
        }
        gOff.drawLine(xpts[j],ypts[j],xpts[j+1],ypts[j+1]);
      }
    }
  }

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

  double norm(double ph[][][]) {
    int i,j;
    double p;

    p = 0.0;
    for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        p += (ph[i][j][0]*ph[i][j][0]+ph[i][j][1]*ph[i][j][1])*dx*dy;
      }
    }
    return ( p );
  }

  double meanPotential(double ph[][][]) {
    int i,j;
    double p;

    p = 0.0;
    for (i=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        p += vv[i][j]*(ph[i][j][0]*ph[i][j][0]+ph[i][j][1]*ph[i][j][1]);
      }
    }
    return ( p*dx*dy );
  }

  double meanKinetic(double ph[][][]) {
    int i,j;
    double p,d2phRe,d2phIm,h2;

    h2 = dx*dx;
    p = 0.0;
    for (i=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        d2phRe = (ph[i+1][j][0]+ph[i-1][j][0]+ph[i][j+1][0]+ph[i][j-1][0]-4.0*ph[i][j][0])/h2;
        d2phIm = (ph[i+1][j][1]+ph[i-1][j][1]+ph[i][j+1][1]+ph[i][j-1][1]-4.0*ph[i][j][1])/h2;
        p += (ph[i][j][0]*d2phRe+ph[i][j][1]*d2phIm);
      }
    }
    return ( -0.5*p*dx*dy );
  }

  double meanPosX(double ph[][][]) {
    int i,j;
    double p;

    p = 0.0;
    for (i=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        p += i*dx*(ph[i][j][0]*ph[i][j][0]+ph[i][j][1]*ph[i][j][1]);
      }
    }
    return ( p*dx*dy-8.0 );
  }

  double meanPosY(double ph[][][]) {
    int i,j;
    double p;

    p = 0.0;
    for (i=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        p += j*dy*(ph[i][j][0]*ph[i][j][0]+ph[i][j][1]*ph[i][j][1]);
      }
    }
    return ( p*dx*dy-8.0 );
  }


/*-------------------------- set loss energy ----------*/

  void setLoss() {
    int ie;

    for (ie=0;ie<NNe;ie++) {
      steepestDescent(psi[ie], vv, 0.02);
	}
  }

  void steepestDescent(double ph[][][], double v[][], double dump) {
    int i,j,k;
    double ee,d2phRe,d2phIm,h2;

    h2 = 2.0*dx*dx;
    ee = meanKinetic(ph)+meanPotential(ph);

    for (i=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
	    wrk[i][j][0] = -(ph[i+1][j][0]+ph[i-1][j][0]+ph[i][j+1][0]+ph[i][j-1][0]-4.0*ph[i][j][0])/h2+(v[i][j]-ee)*ph[i][j][0];
        wrk[i][j][1] = -(ph[i+1][j][1]+ph[i-1][j][1]+ph[i][j+1][1]+ph[i][j-1][1]-4.0*ph[i][j][1])/h2+(v[i][j]-ee)*ph[i][j][0];
      }
    }

    for (i=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
		ph[i][j][0] -= dump*wrk[i][j][0];
		ph[i][j][1] -= dump*wrk[i][j][1];
      }
    }

    normalize(ph);
  }

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

  void GramSchmidt(int stateMax) {
    int i,j,istate,ist;
    double pRe,pIm;

    normalize(psi[0]);

    for (istate=1; istate<stateMax; istate++) {
      for (ist=0; ist<istate; ist++) {
        pRe = 0.0; pIm = 0.0;
        for (i=1;i<NNx-1;i++) {
          for (j=1;j<NNy-1;j++) {
			pRe += (psi[ist][i][j][0]*psi[istate][i][j][0] + psi[ist][i][j][1]*psi[istate][i][j][1])*dx*dy;
			pIm += (psi[ist][i][j][0]*psi[istate][i][j][1] - psi[ist][i][j][1]*psi[istate][i][j][0])*dx*dy;
		  }
        }
        for (i=1;i<NNx-1;i++) {
          for (j=1;j<NNy-1;j++) {
			psi[istate][i][j][0] -= (pRe*psi[ist][i][j][0]-pIm*psi[ist][i][j][1]);
			psi[istate][i][j][1] -= (pRe*psi[ist][i][j][1]+pIm*psi[ist][i][j][0]);
		  }
        }
	  }
      normalize(psi[istate]);
    }
  }

  double innerNorm(double f[][][], double g[][][]) {
    int i,j;
    double pRe,pIm;

    pRe = 0.0;
    pIm = 0.0;
    for (i=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        pRe += (f[i][j][0]*g[i][j][0] + f[i][j][1]*g[i][j][1])*dx*dy;
        pIm += (f[i][j][0]*g[i][j][1] - f[i][j][1]*g[i][j][0])*dx*dy;
      }
    }
  
    return ( Math.sqrt(pRe*pRe+pIm*pIm) );
  }

  void normalize(double f[][][]) {
    int i,j,k;
	double a;
  
    a = Math.sqrt(norm(f));
    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
		f[i][j][0] = f[i][j][0]/a;
		f[i][j][1] = f[i][j][1]/a;
      }
    }
  }

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

  void setInitialCondition() {
    t = 0.0;
    setWave(psi[0], 8.0,8.0, 0.85, 0.0, 0.0);
    setWave(psi[1],12.0,4.0, 0.85, 0.0, 0.0);
    GramSchmidt(NNe);
    setPotential();
    poisson(5000);
  }

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

  void setPotential() {
    int i,j;

    for (i=0; i<NNx; i++) {
      for (j=1;j<NNy-1;j++) {
        vv[i][j] = 0.0;
        vvext[i][j] = 0.0;
        vvh[i][j] = 0.0;
        vvx[i][j] = 0.0;
        vvc[i][j] = 0.0;
      }
    }
    setVext();
  }

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

    x0 = 8.0;
    y0 = 8.0;
    k = 0.25;

    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        x = i*dx;
        y = j*dy;
        vvext[i][j] = k*((x-x0)*(x-x0)+(y-y0)*(y-y0));
      }
    }
  }

/* ----------------------------- set wave packet -------------- */

  void setWave(double ph[][][], double xPos, double yPos, double waveWidth, double kx, double ky) {
    int i,j;
    double a,x,y,phAb,phPh;

    for (i=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        x = i*dx;
        y = j*dy;
        phAb = Math.exp(-((x-xPos)*(x-xPos)+(y-yPos)*(y-yPos))/(4.0*waveWidth*waveWidth) );
        phPh = kx*x+ky*y;
        ph[i][j][0] = phAb*Math.cos(phPh);
        ph[i][j][1] = phAb*Math.sin(phPh);
      }
    }
    for (i=1;i<NNx-1;i++) {
      ph[i][0][0] = 0.0;
      ph[i][0][1] = 0.0;
      ph[i][NNy-1][0] = 0.0;
      ph[i][NNy-1][1] = 0.0;
    }
    for (j=1;j<NNy-1;j++) {
      ph[0][j][0] = 0.0;
      ph[0][j][1] = 0.0;
      ph[NNx-1][j][0] = 0.0;
      ph[NNx-1][j][1] = 0.0;
    }

    a = Math.sqrt(norm(ph));
	
    for (i=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        ph[i][j][0] = ph[i][j][0]/a;
        ph[i][j][1] = ph[i][j][1]/a;
      }
    }
  }

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

  void timeEvolution() {
    int i,ie;

    if (started==1) {

      for (i=0; i<1; i++) {
	  
        t = t + 0.5*dt;
        for (ie=0; ie<NNe; ie++) {
          kxStep(psi[ie],0.5*dt);
          kyStep(psi[ie],0.5*dt);
        }

        setVeff();

        for (ie=0; ie<NNe; ie++) {
          phaseStep(psi[ie],dt);
        }

        t = t + 0.5*dt;
        for (ie=0; ie<NNe; ie++) {
          kxStep(psi[ie],0.5*dt);
          kyStep(psi[ie],0.5*dt);
        }
      }

      for (ie=0; ie<NNe; ie++) {
		energy[ie]=meanKinetic(psi[ie])+meanPotential(psi[ie]);
	  }

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

    if (lossFlag==1) {
      setLoss();
      GramSchmidt(NNe);
    }

  }


/* ----------------- effective potential <-- electron density ------- */

  void setVeff() {
    int i,j;

    setRho(NNe);
    poisson(50);
    setVcx();
    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        vv[i][j] = vvext[i][j]+vvh[i][j]+vvx[i][j]+vvc[i][j];
      }
    }
  }

  void setRho(int istateMax) {
    int i,j,ie;
    double g;

    g = 1.0;

    for (i=1; i<NNx-1; i++) {
      for (j=1; j<NNy-1; j++) {
        rho[i][j] = 0.0;
        for (ie=0; ie<istateMax; ie++) {
          rho[i][j] += g*(psi[ie][i][j][0]*psi[ie][i][j][0]+psi[ie][i][j][1]*psi[ie][i][j][1])/lz;
        }
      }
    }
  }

  void poisson(int iterMax) {
    int i,j,iter;
    double h2;

    h2 = 2.0*3.1416*dx*dx;
    for (iter=1; iter<iterMax; iter++) {
      for (i=1; i<NNx-1; i++) {
        for (j=1; j<NNy-1; j++) {
          vvh[i][j] = vvh[i][j]+0.33*(vvh[i+1][j]+vvh[i-1][j]+vvh[i][j+1]+vvh[i][j-1]-4.0*vvh[i][j]+h2*rho[i][j]);
        }
      }
    }
  }

  /* LDA :  J. P. Perdew and A. Zunger; Phys. Rev., B23, 5048 (1981) */

  void setVcx() {
    int i,j;
    double rh,c1,drho;

    c1 = -0.984745022;
    drho = 0.01;
    for (i=1; i<NNx-1; i++) {
      for (j=1; j<NNy-1; j++) {
        rh = rho[i][j];
        vvx[i][j] = c1*Math.pow(rh,0.33333333);
        vvc[i][j] =  eeCorrelation(rh)+rh*(eeCorrelation(rh+drho)-eeCorrelation(rh))/drho;
      }
    }
  }

  double eeCorrelation(double rh) {
    double r,ec;

    r = 0.6204/(Math.pow(rh,0.33333333)+1.0e-20);
    if (r>=1.0) {
      ec = -0.1423/(1.0+1.0529*Math.sqrt(r)+0.3334*r);
    } else {
      ec = -0.0480-0.0116*r+(0.0311+0.0020*r)*Math.log(r);
    }
    return ( ec );
  }


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

  void kxStep(double ph[][][],double deltat) {
    int i,j;
    double a,aaAb,auAb;

    a = 4.0*dx*dx/deltat;
    for (j=1; j<NNy-1; j++) {
      for (i=1; i<NNx-1; i++) {
        bbRe[i] = 2.0;
        bbIm[i] = a;
        aaRe[i] = -2.0;
        aaIm[i] = a;
        bRe[i] = bbRe[i]*ph[i][j][0]-bbIm[i]*ph[i][j][1] - ph[i+1][j][0] - ph[i-1][j][0];
        bIm[i] = bbRe[i]*ph[i][j][1]+bbIm[i]*ph[i][j][0] - ph[i+1][j][1] - ph[i-1][j][1];
      }

      aaAb = aaRe[1]*aaRe[1]+aaIm[1]*aaIm[1];
      uRe[1] = aaRe[1]/aaAb;
      uIm[1] = -aaIm[1]/aaAb;
      ph[1][j][0] = bRe[1]*uRe[1] - bIm[1]*uIm[1];
      ph[1][j][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;
        ph[i][j][0] = (bRe[i]-ph[i-1][j][0])*uRe[i] - (bIm[i]-ph[i-1][j][1])*uIm[i];
        ph[i][j][1] = (bRe[i]-ph[i-1][j][0])*uIm[i] + (bIm[i]-ph[i-1][j][1])*uRe[i];
      }

      for (i=NNx-3; i>=1; i--) {
        ph[i][j][0] -= ph[i+1][j][0]*uRe[i] - ph[i+1][j][1]*uIm[i];
        ph[i][j][1] -= ph[i+1][j][0]*uIm[i] + ph[i+1][j][1]*uRe[i];
      }
    }
  }

  void kyStep(double ph[][][],double deltat) {
    int i,j;
    double a,aaAb,auAb;

    a = 4.0*dx*dx/deltat;
    for (i=1; i<NNx-1; i++) {
      for (j=1; j<NNy-1; j++) {
        bbRe[j] = 2.0;
        bbIm[j] = a;
        aaRe[j] = -2.0;
        aaIm[j] = a;
        bRe[j] = bbRe[j]*ph[i][j][0]-bbIm[j]*ph[i][j][1]-ph[i][j+1][0]-ph[i][j-1][0];
        bIm[j] = bbRe[j]*ph[i][j][1]+bbIm[j]*ph[i][j][0]-ph[i][j+1][1]-ph[i][j-1][1];
      }

      aaAb = aaRe[1]*aaRe[1]+aaIm[1]*aaIm[1];
      uRe[1] = aaRe[1]/aaAb;
      uIm[1] = -aaIm[1]/aaAb;
      ph[i][1][0] = bRe[1]*uRe[1] - bIm[1]*uIm[1];
      ph[i][1][1] = bIm[1]*uRe[1] + bRe[1]*uIm[1];
      for (j=2; j<NNy-1; j++) {
        auAb = (aaRe[j]-uRe[j-1])*(aaRe[j]-uRe[j-1])+(aaIm[j]-uIm[j-1])*(aaIm[j]-uIm[j-1]);
        uRe[j] = (aaRe[j]-uRe[j-1])/auAb;
        uIm[j] = -(aaIm[j]-uIm[j-1])/auAb;
        ph[i][j][0] = (bRe[j]-ph[i][j-1][0])*uRe[j] - (bIm[j]-ph[i][j-1][1])*uIm[j];
        ph[i][j][1] = (bRe[j]-ph[i][j-1][0])*uIm[j] + (bIm[j]-ph[i][j-1][1])*uRe[j];
      }

      for (j=NNy-3; j>=1; j--) {
        ph[i][j][0] -= ph[i][j+1][0]*uRe[j] - ph[i][j+1][1]*uIm[j];
        ph[i][j][1] -= ph[i][j+1][0]*uIm[j] + ph[i][j+1][1]*uRe[j];
      }
    }
  }


  void phaseStep(double ph[][][],double deltat) {
    int i,j;
    double th,cs,sn,phr,phi;

    for (i=1; i<NNx-1; i++) {
      for (j=1; j<NNy-1; j++) {
        th = deltat*vv[i][j];
        cs = Math.cos(th); sn = Math.sin(th);
        phr = ph[i][j][0];
        phi = ph[i][j][1];
        ph[i][j][0] = cs*phr+sn*phi;
        ph[i][j][1] = cs*phi-sn*phr;
      }
    }
  }

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

}