/* transmission coefficient QED2D */
/* coded by Ikeuchi Mitsuru   */
/* ver 0.0.1   2005.12.25     */

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

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

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

  Choice ch_view;
  Button bt_reset, bt_startStop;
  Scrollbar sc_efield, sc_hight, sc_px0, sc_py0;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 30;

  int resetSW = 0;
  int started = 1;
  int dispMode = 3;

  double vvXpos = 64.0;
  double vvWidth = 6.0;
  double vvhight = 0.2;
  double fieldStrength = 0.01;
  double px0 = 0.5;
  double py0 = 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/2.0;
  double dy = 1.0/2.0;

  double dt = 1.0*dx*dx;

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

  double phRe[][] = new double[NNx][NNy];
  double phIm[][] = new double[NNx][NNy];  

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

  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,400);
    setBackground(Color.white);

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

    ch_view = new Choice();
    ch_view.addItem("density");
    ch_view.addItem("phase");
    ch_view.addItem("momentum");
    ch_view.addItem("grid");
    ch_view.addItemListener(this);
    ch_view.select("grid");

    bt_reset= new Button("reset");
    bt_reset.addActionListener(this);

    bt_startStop = new Button("start/stop");
    bt_startStop.addActionListener(this); 

    sc_efield = new Scrollbar(Scrollbar.HORIZONTAL,100,10,10,1010);
    sc_efield.addAdjustmentListener(this);

    sc_hight = new Scrollbar(Scrollbar.HORIZONTAL,20,10,0,110);
    sc_hight.addAdjustmentListener(this);

    sc_px0 = new Scrollbar(Scrollbar.HORIZONTAL,50,10,0,110);
    sc_px0.addAdjustmentListener(this);

    sc_py0 = new Scrollbar(Scrollbar.HORIZONTAL,0,10,0,110);
    sc_py0.addAdjustmentListener(this);

    addMouseListener(this);  
    addMouseMotionListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,6,5,0));
	/* pnl.add(new Label("  ")); */
    pnl.add(bt_reset);
    pnl.add(bt_startStop);
    pnl.add(sc_px0);
    pnl.add(sc_hight);
    pnl.add(sc_efield);
    pnl.add(ch_view);
    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() == ch_view){
      dispMode = ch_view.getSelectedIndex();
    }
  }

  public void actionPerformed(ActionEvent ev){
    if(ev.getSource() == bt_reset){ 
      resetSW = 1;
    } else if (ev.getSource() == bt_startStop){
      if (started==0) {
        started = 1;
      } else {
        started = 0;
      }
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_efield) { 
      fieldStrength = 1.0e-4*(double)(sc_efield.getValue());
      setPotential(vvXpos,vvhight,fieldStrength);
    } else if (ev.getSource() == sc_hight) { 
      vvhight = 0.01*(double)(sc_hight.getValue());
      setPotential(vvXpos,vvhight,fieldStrength);
    } else if (ev.getSource() == sc_px0) { 
      px0 = 0.01*(double)(sc_px0.getValue());
    } else if (ev.getSource() == sc_py0) { 
      py0 = 0.01*(double)(sc_py0.getValue());
    }
  }

  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 gx;
    double kk,uu;

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


    if (dispMode==0) {
      densityPlot(5000.0, 0);

    } else if (dispMode==1) {
	  densityPlot(5000.0,1);

    } else if (dispMode==2) {
      momentumPlot(5000.0);

    } else if (dispMode==3) {
      gridPlot(20000.0);

    }

    gOff.setColor(Color.black);
    gOff.drawString("t="+(int)(t*10)/10.0+" au",10,40);
    gOff.drawString("t="+(int)(t*2.4189+0.5)/100.0+" fs",10,60);
    if (started==1) { 
      gOff.drawString("start",630/6*1+10,40);
	} else {
      gOff.drawString("stop",630/6*1+10,40);
	}
    gOff.drawString("px0="+(int)(px0*100.0+0.5)/100.0+" au",630/6*2+10,40);
    gOff.drawString("hight="+(int)(vvhight*100.0+0.5)/100.0+" au",630/6*3+10,40);
    gOff.drawString("(="+(int)(vvhight*27.2*100.0+0.5)/100.0+" eV)",630/6*3+10,60);
    gOff.drawString("field="+(int)(fieldStrength*1000.0+0.5)/1000.0+" au",630/6*4+10,40);
    gOff.drawString("(="+(int)(fieldStrength*51.0*1000+0.5)/1000.0+" GV/m)",630/6*4+10,60);
    gOff.drawString("view",630/6*5+10,40);

    gOff.drawString("Lx x Ly="+(int)(NNx*dx*52.918+0.5)+" x "+(int)(NNy*dy*52.918+0.5)+" pm",630/6*1+10,60);

    kk = kineticEnergy();
    uu = potentialEnergy();
    gOff.drawString("K="+(float)(kk)+" au",500,80);
    gOff.drawString("U="+(float)(uu)+" au",500,100);
    gOff.drawString("E="+(float)(kk+uu)+" au",500,120);

    gOff.drawString("nLeft="+(float)(nLeft())+" ",30,360);
    gOff.drawString("nRight="+(float)(nRight())+" ",180,360);
    gOff.drawString("norm="+norm()+" ",30,380);
  }

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

  void densityPlot(double weight, int imode) {
    int i,j;
    int gbx,gby,mag;
    double ph2,col;

    gbx = 30;
    gby = 70;
    mag = 1;
    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++) {
        if (vv[i][j]>0.0) {
          gOff.setColor(Color.getHSBColor(0.33f, 0.95f, 0.6f));
          gOff.fillRect(gbx+mag*i,gby+mag*j, mag, mag);
        }
        ph2 = weight*(phRe[i][j]*phRe[i][j]+phIm[i][j]*phIm[i][j]);
        if (ph2>0.01) {
          if (ph2>1.0) {
            ph2 = 1.0;
		  }
          if (imode==0) {
            gOff.setColor(Color.getHSBColor(0.15f, 0.95f, (float)(ph2)));
          } else if (imode==1) {
            gOff.setColor(Color.getHSBColor((float)(0.5+0.5*phRe[i][j]/Math.sqrt(ph2/weight)), 0.9f, (float)(ph2)));
          }
          gOff.fillRect(gbx+mag*i,gby+mag*j, mag, mag);
        }
      }
    }
  }

  void momentumPlot(double weight) {
    int i,j;
    int gbx,gby,mag;
    double amp,ph2,px,pxphRe,pxphIm,py,pyphRe,pyphIm;

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

    for (i=2; i<NNx-1; i+=4) {
      for (j=2; j<NNy-1; j+=4) {
        if (vv[i][j]>0.0) {
          gOff.setColor(Color.getHSBColor(0.33f, 0.95f, 0.6f));
          gOff.fillRect(gbx+i-2,gby+j-2, 4, 4);
        }
        ph2 = weight*(phRe[i][j]*phRe[i][j]+phIm[i][j]*phIm[i][j]);
        if (ph2>0.01) {
          pxphRe = (phIm[i+1][j]-phIm[i-1][j])/(2*dx);
          pxphIm = (-phRe[i+1][j]+phRe[i-1][j])/(2*dx);
          px = (phRe[i][j]*pxphRe + phIm[i][j]*pxphIm);
          pyphRe = (phIm[i][j+1]-phIm[i][j-1])/(2*dy);
          pyphIm = (-phRe[i][j+1]+phRe[i][j-1])/(2*dy);
          py = (phRe[i][j]*pyphRe + phIm[i][j]*pyphIm);

          if (px>0.0) {
            gOff.setColor(Color.getHSBColor(0.66f, 0.9f, 0.9f));
          } else {
            gOff.setColor(Color.getHSBColor(0.01f, 0.9f, 0.9f));
          }
		  gOff.drawLine(gbx+mag*i,gby+mag*j, gbx+mag*i+(int)(amp*px), gby+mag*j+(int)(amp*py));
        }
      }
    }
  }

  void gridPlot(double weight) {
    int i,j,gbx,gby;
    int cx,cy,cz;
    double z,px,py,pz,sc;
    double v1,v2,ph2;

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

    gbx = -100+(int)(cx/sc);
    gby = -40+(int)(cy/sc);
    for(j=0; j<NNy; j += 4) {
      for (i=0; i<NNx; i++) { 
        ph2 = phRe[i][j]*phRe[i][j]+phIm[i][j]*phIm[i][j];
        z = weight*ph2 + 0.005*weight*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 (weight*ph2>1.0) {
          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 = phRe[i][j]*phRe[i][j]+phIm[i][j]*phIm[i][j];
        z = weight*ph2 + 0.005*weight*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 (weight*ph2>1.0) {
          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]);
      }
    }
  }

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

  double norm() {
    int i,j;
    double p;

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

  double kineticEnergy() {
    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 = (phRe[i+1][j]+phRe[i-1][j]+phRe[i][j+1]+phRe[i][j-1]-4.0*phRe[i][j])/h2;
          d2phIm = (phIm[i+1][j]+phIm[i-1][j]+phIm[i][j+1]+phIm[i][j-1]-4.0*phIm[i][j])/h2;
          p += (phRe[i][j]*d2phRe+phIm[i][j]*d2phIm);
      }
    }
    return ( -0.5*p*dx*dy );
  }

  double potentialEnergy() {
    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]*(phRe[i][j]*phRe[i][j]+phIm[i][j]*phIm[i][j]);
      }
    }
    return ( p*dx*dy );
  }

  double nLeft() {
    int i,j;
    double p;

    p = 0.0;
    for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        if (i*dx<vvXpos) {
          p += phRe[i][j]*phRe[i][j]+phIm[i][j]*phIm[i][j];
        }
      }
    }
    return ( p*dx*dy );
  }

  double nRight() {
    int i,j;
    double p;

    p = 0.0;
    for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        if (i*dx>vvXpos+vvWidth) {
          p += phRe[i][j]*phRe[i][j]+phIm[i][j]*phIm[i][j];
        }
      }
    }
    return ( p*dx*dy );
  }

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

  void setInitialCondition() {
    t = 0.0;
    setWave(32.0,64.0, 7.0, px0, py0);
    setPotential(vvXpos,vvhight,fieldStrength);
  }

  void setPotential(double xpos, double hight, double eField) {
    int i,j;
    double x,y;

    for (i=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        x = i*dx;
        if (x>xpos) {
          vv[i][j] = hight-eField*(x-xpos);
          if (vv[i][j]<-0.1) vv[i][j] = -0.1;
        } else {
          vv[i][j] = 0.0;
        }
      }
    }
    vvWidth = setWidth(xpos);

  }

  double setWidth(double vPos) {
    int i;

    for (i=1;i<NNx-1;i++) {
      if (vv[i][NNy/2]<0.0) {
        return(i*dx-vPos);
      }
    }
    return ( NNx*dx-vPos );
  }


  void setWave(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;
        phRe[i][j] = phAb*Math.cos(phPh);
        phIm[i][j] = phAb*Math.sin(phPh);
      }
    }
    for (i=1;i<NNx-1;i++) {
      phRe[i][0] = 0.0;
      phIm[i][0] = 0.0;
      phRe[i][NNy-1] = 0.0;
      phIm[i][NNy-1] = 0.0;
    }
    for (j=1;j<NNy-1;j++) {
      phRe[0][j] = 0.0;
      phIm[0][j] = 0.0;
      phRe[NNx-1][j] = 0.0;
      phIm[NNx-1][j] = 0.0;
    }

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

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

  void timeEvolution() {
    int i;

    if (resetSW==1) {
      setInitialCondition();
      resetSW = 0;
      started = 1;
    }

    if (started==1) {
      for (i=0; i<2; i++) {
        t = t + dt;
        phaseStep();
        kxStep();
        kyStep();
        phaseStep();
      }
    }
  }

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

  void kxStep() {
    int i,j;
    double bbRe,bbIm,aaRe,aaIm,aaAb,auAb;

    bbRe = 2.0;
    bbIm = 4*dx*dx/dt;
    aaRe = -2.0;
    aaIm = 4*dx*dx/dt;
    aaAb = aaRe*aaRe+aaIm*aaIm;

    for (j=1; j<NNy-1; j++) {
      for (i=1; i<NNx-1; i++) {
        bRe[i] = bbRe*phRe[i][j]-bbIm*phIm[i][j]-phRe[i+1][j]-phRe[i-1][j];
        bIm[i] = bbRe*phIm[i][j]+bbIm*phRe[i][j]-phIm[i+1][j]-phIm[i-1][j];
      }

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

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

  void kyStep() {
    int i,j;
    double bbRe,bbIm,aaRe,aaIm,aaAb,auAb;

    bbRe = 2.0;
    bbIm = 4*dx*dx/dt;
    aaRe = -2.0;
    aaIm = 4*dx*dx/dt;
    aaAb = aaRe*aaRe+aaIm*aaIm;

    for (i=1; i<NNx-1; i++) {
      for (j=1; j<NNy-1; j++) {
        bRe[j] = bbRe*phRe[i][j]-bbIm*phIm[i][j]-phRe[i][j+1]-phRe[i][j-1];
        bIm[j] = bbRe*phIm[i][j]+bbIm*phRe[i][j]-phIm[i][j+1]-phIm[i][j-1];
      }

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

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

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

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

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

}