/** applet No. 1055
 *
 * EM wave - finite difference time domain (FDTD) 3D
 * - multi thread  FDTD3D :: display - asynchronous
 *
 * Created by Ikeuchi Mitsuru on November 11 2006.
 * Copyright (c) 2006-2007 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2006.11.11  created
 * ver 0.0.2   2007.06.02  improved code
 *
 */

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

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

  // ------------------------------------ preset field -----------

  Thread th = null; // for run()-paint() thread
  FDTD3D fdtd = null;
  PlotFDTD3D pf3d = new PlotFDTD3D();

  // for event
  Choice ch_view;
  Button bt_reset, bt_startStop, bt_step;
  Scrollbar sc_lambda, sc_mag;

  // for off-paint buffer
  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int dgX, dgY, dgXb,dgYb; // mouse

  int sleepTime = 100;
  int dispMode = 0;
  int thCount = 0;


  // ------------------------------ applet main thread -----------

  public void init() {

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

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

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

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

    sc_lambda= new Scrollbar(Scrollbar.HORIZONTAL,8,10,8,42);
    sc_lambda.addAdjustmentListener(this);

    sc_mag= new Scrollbar(Scrollbar.HORIZONTAL,100,10,10,310);
    sc_mag.addAdjustmentListener(this);

    ch_view = new Choice();
    ch_view.add("Ez dens");
    ch_view.add("E vector");
    ch_view.add("H vector");
    ch_view.add("Ez wave");
    ch_view.addItemListener(this); 
    ch_view.select("Ez dens");

    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(bt_step);
    pnl.add(sc_lambda);
    pnl.add(sc_mag);
    pnl.add(ch_view);
    add(pnl,"North");

  }

  public void start() {
    if (fdtd == null) {
      fdtd = new FDTD3D();
      fdtd.start();
    }
    if (th == null) {
      th = new Thread(this);
      th.start();
    }
  }

  public void stop() {
    if (th != null) th = null;
    if (fdtd != null) fdtd = null;
  }


  // ---------------------------------- event listener -----------

  public void itemStateChanged(ItemEvent ev){
    if (ev.getSource() == ch_view){
      dispMode = ch_view.getSelectedIndex();
    }
  }

  public void actionPerformed(ActionEvent ev){
    if(ev.getSource() == bt_reset){ 
      fdtd.reset();
      thCount = 0;
    } else if (ev.getSource() == bt_startStop){
      if (fdtd.getStartSW()==0) {
        fdtd.setStartSW(1);
      } else {
        fdtd.setStartSW(0);
      }
    } else if (ev.getSource() == bt_step){
      if (fdtd.getStartSW()==0) {
        fdtd.setStepSW();
      }
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_lambda) { 
      fdtd.setLambda( (double)sc_lambda.getValue() );
    } else if(ev.getSource() == sc_mag) { 
      pf3d.setMag( 0.01*(double)(sc_mag.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;
    }
    pf3d.setViewDelta( (dgX-dgXb), (dgY-dgYb) );
  }



  // ========================= run() - paint() loop ==============

  public void run() {
    while (th != null) {
      try {
        Thread.sleep(sleepTime);
      } catch (InterruptedException e) { }
      thCount += 1;
      offPaint();
      repaint();
    }
  }

  public void update(Graphics g) {
    paint(g);
  }

  public void paint(Graphics g) {
    g.drawImage(imgOff,0,0,this);
  }

  void offPaint() {
    if (fdtd != null) {
      pf3d.plotGraph(gOff, fdtd, dispMode, thCount);
    }
  }
}



// =============================== draw graph 2D class ===========

class PlotFDTD3D {

  private Graphics gOff;
  private FDTD3D fdtd;

  private int NNx;
  private int NNy;
  private int NNz;
  private double dx = 1.0;
  private double dy = 1.0;
  private double dz = 1.0;
  private double xMax;
  private double yMax;
  private double zMax;
  private double cxx;
  private double cyy;
  private double czz;

  private double dispScale;
  private double viewScale = 1.0;
  private double scale;
  private double magn = 1.0;

  private int xImageLoc = 60;
  private int yImageLoc = 100;

  private double viewTheta = -15.0*3.14/180.0;
  private double viewFai = -72.0*3.14/180.0;
  private double cosTh = Math.cos(viewTheta); 
  private double sinTh = Math.sin(viewTheta);
  private double cosFi = Math.cos(viewFai);
  private double sinFi = Math.sin(viewFai);

  private double wkx[] = new double[8]; 
  private double wky[] = new double[8]; 
  private double wkz[] = new double[8]; 

  private double pwkx[] = new double[8]; 
  private double pwky[] = new double[8]; 
  private double pwkz[] = new double[8]; 

  private int boxp[][] = {
    {0,1},{0,2},{0,4},{1,3},{1,5},{2,3},{2,6},{4,5},{4,6},{3,7},{5,7},{6,7}
  };

  private int xpts[] = new int[200];
  private int ypts[] = new int[200];
  private double zval[] = new double[200];

  private int drawCount = 0;


  // ------------------------------- class constructor -----------

  PlotFDTD3D() {
  }


  // ----------------------------------- class methods -----------

  public void plotGraph(Graphics gg, FDTD3D fd, int dispMode, int thCount) {
    int fdtdCount,gb;

    gOff = gg;
    fdtd = fd;
    if (fd==null) return;
    fdtdCount = fdtd.getCount();
    if (fdtdCount==drawCount) return;
    drawCount = fdtdCount;
    setGlobal();

    gOff.setColor(Color.white);
    gOff.fillRect(0,0,630,420);

    if (dispMode==0) {
      densityPlot(0,2.0*magn,2); // Ez density

    } else if (dispMode==1) {
      densityPlot(1,20.0*magn,4);  // E vector

    } else if (dispMode==2) {
      densityPlot(2,20.0*magn,4); // H vector

    } else if (dispMode==3) {
      gridPlot(0, NNz/2, 20.0*magn); // Ez(@z=0) wave;

    } else if (dispMode==4) {
      ;

    }

    gOff.setColor(Color.black);
    gOff.drawString("t="+(int)(fdtd.getTime()+0.5)+" ",20,40);
    if (fdtd.getStartSW()==1) {
      gOff.drawString("start",630/6*1+20,40);
    } else {
      gOff.drawString("stop",630/6*1+20,40);
    }
    gOff.drawString("lambda="+(int)(fdtd.getLambda()+0.5)+" ",630/6*3+10,40);
    gOff.drawString("mag="+(int)(magn*100.0+0.5)+"%",630/6*4+10,40);

    gOff.drawString("grid=80x64x64",500,60);

    /*
    gb = 500;
    gOff.setColor(Color.green);
    gOff.drawString("vaccum",gb,325);
    gOff.setColor(Color.orange);
    gOff.drawString("metal",gb,340);
    gOff.setColor(Color.red);
    gOff.drawString("dielectrics",gb,355);
    gOff.setColor(Color.blue);
    gOff.drawString("ferromagnetic",gb,370);
    gOff.setColor(Color.magenta);
    gOff.drawString("diele+ferro",gb,385);
    gOff.setColor(Color.gray);
    gOff.drawString("resistor",gb,400);
    */

    gb = 500;
    gOff.setColor(Color.black);
    gOff.drawString("thread count",gb,340);
    gOff.drawString("FDTD = "+fdtd.getCount()+" ",gb,360);
    gOff.drawString("Plot = "+thCount+" ",gb,380);
  }

  private void setGlobal() {

    NNx = fdtd.getNNx();
    NNy = fdtd.getNNy();
    NNz = fdtd.getNNz();
    xMax = NNx*dx;
    yMax = NNy*dy;
    zMax = NNz*dz;
    cxx = xMax/2.0;
    cyy = yMax/2.0;
    czz = zMax/2.0;
    dispScale = (4.0*NNx)/xMax;
    scale = dispScale*viewScale;
  }


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

  private void densityPlot(int mode, double mag,int idiv) {
    int i,j,k,ix,iy,ir,fp;
    int iBegin, iEnd, iInc, jBegin, jEnd, jInc, kBegin, kEnd, kInc;
    double x,y,z,dens,col;

    fp = drawBox(0);

    iBegin = 0; iEnd = NNx; iInc = idiv;
    if (wkx[fp]>10.0) {
      iBegin = NNx-1; iEnd = 0; iInc = -idiv;
    }
    jBegin = 0; jEnd = NNy; jInc = idiv;
    if (wky[fp]>10.0) {
      jBegin = NNy-1; jEnd = 0; jInc = -idiv;
    }
    kBegin = 0; kEnd = NNz; kInc = idiv;
    if (wkz[fp]>10.0) {
      kBegin = NNz-1; kEnd = 0; kInc = -idiv;
    }

    for (i=iBegin; ((i<iEnd && iInc==idiv) || (i>=iEnd && iInc==-idiv)) ; i += iInc) {
      x = dx*i;
      for (j=jBegin; ((j<jEnd && jInc==idiv) || (j>=jEnd && jInc==-idiv)) ; j += jInc) {
        y = dy*j;
        for (k=kBegin; ((k<kEnd && kInc==idiv) || (k>=kEnd && kInc==-idiv)) ; k += kInc) {
          z = dz*k;

          if (mode==0) {
            dens = mag*fdtd.getEz(i,j,k);
            col = 0.33f;
            if (dens<0) {
              col = 0.01f;
              dens = -dens;
            }
            if (dens>0.999) dens = 0.999;
            ir = (int)(10.0*dens+0.9);
            if (ir>0) plotDot(x,y,z,ir, (float)col);
	
          } else if (mode==1) {
            plotEVector(x,y,z, 
              mag*fdtd.getEx(i,j,k),
              mag*fdtd.getEy(i,j,k),
              mag*fdtd.getEz(i,j,k) );

          } else if (mode==2) {
            plotHVector(x,y,z, 
              mag*fdtd.getHx(i,j,k),
              mag*fdtd.getHy(i,j,k),
              mag*fdtd.getHz(i,j,k) );
          }
        }
      }
    }
    drawBox(1);
  }

  private void plotDot(double x, double y, double z, int ir, float col) {
    int ix,iy;
    double ppx,ppy,ppz;

    ppy = cosFi*(y-cyy)+sinFi*(z-czz) + cyy; 
    ppz = -sinFi*(y-cyy)+cosFi*(z-czz) + czz;
    ppx = cosTh*(x-cxx)+sinTh*(ppz-czz) + cxx; 
    //ppz = -sinTh*(x-cxx)+cosTh*(ppz-czz) + czz;

    ix = (int)(scale*ppx)+xImageLoc;
    iy = (int)(scale*ppy)+yImageLoc;

    gOff.setColor(Color.getHSBColor(col,0.9f, 0.9f ));
    gOff.fillOval(ix-ir/2,iy-ir/2,ir,ir);
  }

  private void plotEVector(double x, double y, double z, double u, double v, double w) {
    double x2,y2,z2;
    float col;

    x2 = x + u;
    y2 = y + v;
    z2 = z + w;
    if (w>=0.0) {
      col = 0.33f;
    } else {
      col = 0.05f;
    }
    plotLine(x,y,z, x2,y2,z2, col);
  }

  private void plotHVector(double x, double y, double z, double u, double v, double w) {
    double x2,y2,z2;
    float col;

    x2 = x + u;
    y2 = y + v;
    z2 = z + w;
    if (v>=0.0) {
      col = 0.66f;
    } else {
      col = 0.85f;
    }
    plotLine(x,y,z, x2,y2,z2, col);
  }

  private void plotLine(double x, double y, double z, double x2, double y2, double z2, float col) {
    int ix,iy, ix2,iy2;
    double ppx,ppy,ppz, ppx2,ppy2,ppz2;

    ppy = cosFi*(y-cyy)+sinFi*(z-czz) + cyy; 
    ppz = -sinFi*(y-cyy)+cosFi*(z-czz) + czz;
    ppx = cosTh*(x-cxx)+sinTh*(ppz-czz) + cxx; 
    //ppz = -sinTh*(x-cxx)+cosTh*(ppz-czz) + czz;

    ix = (int)(scale*ppx)+xImageLoc;
    iy = (int)(scale*ppy)+yImageLoc;

    ppy2 = cosFi*(y2-cyy)+sinFi*(z2-czz) + cyy; 
    ppz2 = -sinFi*(y2-cyy)+cosFi*(z2-czz) + czz;
    ppx2 = cosTh*(x2-cxx)+sinTh*(ppz2-czz) + cxx; 
    //ppz2 = -sinTh*(x2-cxx)+cosTh*(ppz2-czz) + czz;

    ix2 = (int)(scale*ppx2)+xImageLoc;
    iy2 = (int)(scale*ppy2)+yImageLoc;

    gOff.setColor(Color.getHSBColor(col,0.9f, 0.9f ));
    gOff.drawLine(ix,iy,ix2,iy2);
  }

// draw box

  private int drawBox(int imode) {
    int i,gx,gy,g2x,g2y,farPoint;
    double sc;

    if (imode==0) {
      setBox();
      rotateBox();
    }
    sc = dispScale*viewScale;
    farPoint = findFarPoint();

    for (i=0; i<12; i++) {
      gx = (int)(sc*pwkx[boxp[i][0]])+xImageLoc;
      gy = (int)(sc*pwky[boxp[i][0]])+yImageLoc;
      g2x = (int)(sc*pwkx[boxp[i][1]])+xImageLoc;
      g2y = (int)(sc*pwky[boxp[i][1]])+yImageLoc;
      if (imode==0) {
        if (boxp[i][0]==farPoint || boxp[i][1]==farPoint) { 
          gOff.setColor(Color.lightGray);
          gOff.drawLine(gx,gy, g2x, g2y);
        }
      } else if (imode==1) {
        if (boxp[i][0]!=farPoint && boxp[i][1]!=farPoint) { 
          gOff.setColor(Color.darkGray);
          gOff.drawLine(gx,gy, g2x, g2y);
        }
      }
    }
    return farPoint;
  }

  private int findFarPoint() {
    int i,im;
    double m;

    im = 0; m = pwkz[im];
    for (i=0; i<8; i++) {
      if (pwkz[i]<m) {
        m = pwkz[i];
        im = i;
      }
    }
    return(im);
  }

  private void setBox() {
    int i,ix,iy,iz;
    double ballSize;

    i=0;
    ballSize = 0.0*xMax;

    for (ix=0; ix<2; ix++) {
      for (iy=0; iy<2; iy++) {
        for (iz=0; iz<2; iz++) {
          wkx[i] = -0.5*ballSize+(xMax+ballSize)*(double)(ix);
          wky[i] = -0.5*ballSize+(yMax+ballSize)*(double)(iy);
          wkz[i] = -0.5*ballSize+(zMax+ballSize)*(double)(iz);
          i += 1;
        }
      }
    }
  }

  private void rotateBox() {
    int i;
 
    for (i=0; i<8; i++) {
      pwky[i] = cosFi*(wky[i]-cyy)+sinFi*(wkz[i]-czz) + cyy; 
      pwkz[i] = -sinFi*(wky[i]-cyy)+cosFi*(wkz[i]-czz) + czz;
      pwkx[i] = cosTh*(wkx[i]-cxx)+sinTh*(pwkz[i]-czz) + cxx; 
      pwkz[i] = -sinTh*(wkx[i]-cxx)+cosTh*(pwkz[i]-czz) + czz;
    }  
  }

  //

  void gridPlot(int mode, int k, double mag) {
    int i,j,gbx,gby;
    int cx,cy,cz;
    double ff,z,px,py,pz,sc;

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

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

    for(j=0; j<NNy; j += 2) {
      for (i=0; i<NNx; i++) {
        ff = fdtd.getEz(i,j,k);
        z = mag*ff + 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;
        zval[i] = ff;
      }
      for (i=0; i<NNx-1; i++) {
        gOff.setColor(Color.getHSBColor(lineColor(i),0.6f,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++) {
        ff = fdtd.getEz(i,j,k);
        z = mag*ff + 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;
        zval[j] = ff;
      }
      for (j=0; j<NNy-1; j++) { 
        gOff.setColor(Color.getHSBColor(lineColor(j),0.6f,0.9f));
        gOff.drawLine(xpts[j],ypts[j],xpts[j+1],ypts[j+1]);
      }
    }
  }

  private float lineColor(int i) {
	
    if (zval[i]*zval[i+1]<0.0) return 0.66f;
    if (zval[i]+zval[i+1]>0.02) return 0.33f;
    if (zval[i]+zval[i+1]<-0.02) return 0.05f;
    return 0.66f;
  }


  // ------------------------------------- I/O methods -----------

  public void setViewDelta(double xDisplacement, double yDisplacement) {
    viewTheta += 0.5*3.14159/180.0*xDisplacement;
    viewFai += 0.5*3.14159/180.0*yDisplacement;
    cosTh = Math.cos(viewTheta); 
    sinTh = Math.sin(viewTheta);
    cosFi = Math.cos(viewFai);
    sinFi = Math.sin(viewFai);
  }
  public void setMag(double mg) { magn = mg; }
  public double getMag() { return magn; }
}



// ===================================== FDTD 3D class ===========

class FDTD3D extends Thread {

  private double omega = Math.PI/8.0;
  private double theta = 0.0;

  private double t = 0.0;
  private double dt = 0.5;
  private double dx = 1.0;
  private double dy = 1.0;
  private double dz = 1.0;

  private int NNx = 64+16;
  private int NNy = 48+16;
  private int NNz = 48+16;

  private double Ex[][][] = new double[NNx+1][NNy+1][NNz+1];
  private double Ey[][][] = new double[NNx+1][NNy+1][NNz+1];
  private double Ez[][][] = new double[NNx+1][NNy+1][NNz+1];
  private double Hx[][][] = new double[NNx+1][NNy+1][NNz+1];
  private double Hy[][][] = new double[NNx+1][NNy+1][NNz+1];
  private double Hz[][][] = new double[NNx+1][NNy+1][NNz+1];

  private double eps[][][] = new double[NNx+1][NNy+1][NNz+1];
  private double mue[][][] = new double[NNx+1][NNy+1][NNz+1];
  private double sgm[][][] = new double[NNx+1][NNy+1][NNz+1];

  private int resetSW = 0;
  private int startSW = 1;
  private int stepSW = 1;
  private int sleepTime = 1;
  private int stopSleepTime = 100;
  private int count = 0;


  // ------------------------------- class constructor -----------

  FDTD3D() {
    setInitialCondition();
  }


  // -------------------------------------- thread run -----------

  public void run() {
    while (true) {
      count += 1;
      timeEvolution();
      if (sleepTime>0) {
        try {
          Thread.sleep(sleepTime);
        } catch (InterruptedException e) { }
      }
    }
  }


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

  private void setInitialCondition() {

    count = 0;
    t = 0.0;
    theta = 0.0;

    clearField();
    setMaterials();
  }

  void clearField() {
    int i,j,k;

    for(i=0;i<NNx+1;i++) {
      for(j=0;j<NNy+1;j++) {
        for(k=0;k<NNz+1;k++) {
          Ex[i][j][k] = 0.0;
          Ey[i][j][k] = 0.0;
          Ez[i][j][k] = 0.0;
          Hx[i][j][k] = 0.0;
          Hy[i][j][k] = 0.0;
       	  Hz[i][j][k] = 0.0;
        }
      }
    }
  }


  void setMaterials() {
    int i,j,k;

    for(i=0;i<NNx+1;i++) {
      for(j=0;j<NNy+1;j++) {
        for(k=0;k<NNz+1;k++) {
          eps[i][j][k] = 1.0;
          mue[i][j][k] = 1.0;
       	  sgm[i][j][k] = 0.0;
        }
      }
    }

    for(i=1;i<8;i++) {
      for(j=1;j<NNy-1;j++) {
        for(k=1;k<NNz-1;k++) {
          eps[i][j][k] = 1.0;
          mue[i][j][k] = 1.0;
       	  sgm[i][j][k] = 0.3;
        }
      }
    }
    for(i=NNx-8;i<NNx-1;i++) {
      for(j=1;j<NNy-1;j++) {
        for(k=1;k<NNz-1;k++) {
          eps[i][j][k] = 1.0;
          mue[i][j][k] = 1.0;
       	  sgm[i][j][k] = 0.3;
        }
      }
    }
    for(i=1;i<NNx-1;i++) {
      for(j=1;j<8;j++) {
        for(k=1;k<NNz-1;k++) {
          eps[i][j][k] = 1.0;
          mue[i][j][k] = 1.0;
       	  sgm[i][j][k] = 0.3;
        }
      }
    }
    for(i=1;i<NNx-1;i++) {
      for(j=NNy-8;j<NNy-1;j++) {
        for(k=1;k<NNz-1;k++) {
          eps[i][j][k] = 1.0;
          mue[i][j][k] = 1.0;
       	  sgm[i][j][k] = 0.3;
        }
      }
    }
    for(i=1;i<NNx-1;i++) {
      for(j=1;j<NNy-1;j++) {
        for(k=1;k<8;k++) {
          eps[i][j][k] = 1.0;
          mue[i][j][k] = 1.0;
       	  sgm[i][j][k] = 0.3;
        }
      }
    }
    for(i=1;i<NNx-1;i++) {
      for(j=1;j<NNy-1;j++) {
        for(k=NNz-8;k<NNz-1;k++) {
          eps[i][j][k] = 1.0;
          mue[i][j][k] = 1.0;
       	  sgm[i][j][k] = 0.3;
        }
      }
    }

    for(i=0;i<NNx;i++) {
      for(j=0;j<NNy;j++) {
        eps[i][j][0] = 1.0e200;
        eps[i][j][NNz-1] = 1.0e200;
      }
    }
    for(i=0;i<NNx;i++) {
      for(k=0;k<NNz;k++) {
        eps[i][0][k] = 1.0e200;
        eps[i][NNy-1][k] = 1.0e200;
      }
    }
    for(j=0;j<NNy+1;j++) {
      for(k=0;k<NNz+1;k++) {
        eps[0][j][k] = 1.0e200;
        eps[NNx-1][j][k] = 1.0e200;
      }
    }
  }


  // ---------------------------------- time evolution -----------

  private void timeEvolution(){
    int i;

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

    if (startSW==1) {
      timeStep();

    } else {
      if (stepSW==1) {
        timeStep();
        stepSW = 0;
      }
      try {
        Thread.sleep(stopSleepTime);
      } catch (InterruptedException e) { }
    }
  }

  private void timeStep() {

    ExyzEvolution();
    generator();
    t += 0.5*dt;
    HxyzEvolution();
    t += 0.5*dt;
  }


  private void ExyzEvolution() {
    int i,j,k;

    for (i=0; i<NNx; i++) {
      for (j=1; j<NNy; j++) {
        for (k=1; k<NNz; k++) {
          Ex[i][j][k] = (
            (1.0-sgm[i][j][k]*dt/eps[i][j][k])*Ex[i][j][k] 
            +(dt/eps[i][j][k]/dy)*(Hz[i][j][k]-Hz[i][j-1][k])
            -(dt/eps[i][j][k]/dz)*(Hy[i][j][k]-Hy[i][j][k-1]) );
        }
      }
    }
    for (i=1; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        for (k=1; k<NNz; k++) {
          Ey[i][j][k] = (
            (1.0-sgm[i][j][k]*dt/eps[i][j][k])*Ey[i][j][k] 
            +(dt/eps[i][j][k]/dz)*(Hx[i][j][k]-Hx[i][j][k-1])
            -(dt/eps[i][j][k]/dx)*(Hz[i][j][k]-Hz[i-1][j][k]) );
        }
      }
    }
    for (i=1; i<NNx; i++) {
      for (j=1; j<NNy; j++) {
        for (k=0; k<NNz; k++) {
          Ez[i][j][k] = (
            (1.0-sgm[i][j][k]*dt/eps[i][j][k])*Ez[i][j][k] 
            +(dt/eps[i][j][k]/dx)*(Hy[i][j][k]-Hy[i-1][j][k])
            -(dt/eps[i][j][k]/dy)*(Hx[i][j][k]-Hx[i][j-1][k]) );
        }
      }
    }

    for (j=0; j<NNy; j++) {
      for (k=0; k<NNz; k++) {
        Ey[1][j][k] = 0.0;
        Ez[1][j][k] = 0.0;
      }
    }
    for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        Ex[i][j][1] = 0.0;
        Ey[i][j][1] = 0.0;
      }
    }
  }


  private void generator() {
    int k;
    double Ezt;

    theta += omega*dt;
    Ezt = Math.sin(theta);
    for (k=(NNz-16)/4+8; k<NNz/2; k++) {
      Ez[12][NNy/2][k] = +Ezt;
    }
  }


  private void HxyzEvolution() {
    int i,j,k;

    for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        for (k=1; k<NNz; k++) {
          Hx[i][j][k] = ( Hx[i][j][k] 
            -(dt/mue[i][j][k]/dy)*(Ez[i][j+1][k]-Ez[i][j][k])
            +(dt/mue[i][j][k]/dz)*(Ey[i][j][k+1]-Ey[i][j][k]) );
        }
      }
    }
    for (i=0; i<NNx; i++) {
      for (j=1; j<NNy; j++) {
        for (k=0; k<NNz; k++) {
          Hy[i][j][k] = ( Hy[i][j][k] 
            -(dt/mue[i][j][k]/dz)*(Ex[i][j][k+1]-Ex[i][j][k])
            +(dt/mue[i][j][k]/dx)*(Ez[i+1][j][k]-Ez[i][j][k]) );
        }
      }
    }
    for (i=1; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        for (k=0; k<NNz; k++) {
          Hz[i][j][k] = ( Hz[i][j][k] 
            -(dt/mue[i][j][k]/dx)*(Ey[i+1][j][k]-Ey[i][j][k])
            +(dt/mue[i][j][k]/dy)*(Ex[i][j+1][k]-Ex[i][j][k]) );
        }
      }
    }
  }


  // ------------------------------------- I/O methods -----------

  // control
  public void reset() { resetSW = 1; }
  public void setStartSW(int sw) { startSW = sw; }
  public int getStartSW() { return startSW; }
  public void setStepSW() { stepSW = 1; }
  public void setLambda(double lm) { omega = Math.PI/lm; }
  public double getLambda() { return Math.PI/omega; }

  public int getCount() { return count; }

  // get data
  public double getTime() { return t; }
  public int getNNx() { return NNx; }
  public int getNNy() { return NNy; }
  public int getNNz() { return NNz; }
  public double getEx(int i, int j, int k) { return Ex[i][j][k]; }
  public double getEy(int i, int j, int k) { return Ey[i][j][k]; }
  public double getEz(int i, int j, int k) { return Ez[i][j][k]; }
  public double getHx(int i, int j, int k) { return Hx[i][j][k]; }
  public double getHy(int i, int j, int k) { return Hy[i][j][k]; }
  public double getHz(int i, int j, int k) { return Hz[i][j][k]; }

  public double getSgm(int i, int j, int k) { return sgm[i][j][k]; }
  public double getEps(int i, int j, int k) { return eps[i][j][k]; }
  public double getMue(int i, int j, int k) { return mue[i][j][k]; }
}


/** finite differnce time domain ( FDTD ) 3D
 *
 * electro-magnetic field : Maxwell equation
 *
 * rot H = sgm E + eps dE/dt
 * rot E = -mue dH/dt
 *
 * in 3 dim system: Ex,Ey,Ez,Hx,Hy,Hz
 * dEx/dt = (sgm Ex/eps)+(1/eps)(dHz/dy-dHy/dz)
 * dEy/dt = (sgm Ey/eps)+(1/eps)(dHx/dz-dHz/dx)
 * dEz/dt = (sgm Ez/eps)+(1/eps)(dHy/dx-dHx/dy)
 * dHx/dt = -(1/mue)(dEz/dy-dEy/dz)
 * dHy/dt = -(1/mue)(dEx/dz-dEz/dx)
 * dHz/dt = -(1/mue)(dEy/dx-dEx/dy)
 *
 */