/** applet No. 1053
 *
 * powder scattering - finite difference time domain (FDTD) 2D
 * - multi thread  FDTD2D :: display - asynchronous
 * - float precision - grid 800 x 800
 * - material table 
 *
 * Created by Ikeuchi Mitsuru on November 05 2006.
 * Copyright (c) 2006-2007 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2006.11.05  created
 * ver 0.0.2   2007.06.02  improved code
 *
 */

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

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

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

  Thread th = null; // for run()-paint() thread
  FDTD2D fdtd = null;
  DrawGraph2D dg2d = new DrawGraph2D();

  // 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 = 50;
  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,32,10,8,110);
    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 wave");
    ch_view.add("Ez density");
    ch_view.add("Hxy field");
    ch_view.add("Ez + Hxy");
    ch_view.addItemListener(this); 
    ch_view.select("Ez wave");

    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 FDTD2D();
      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) { 
      dg2d.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;
    }
    dg2d.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) {
      dg2d.plotGraph2D(gOff, fdtd, dispMode, thCount);
    }
  }
}



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

class DrawGraph2D {

  private Graphics gOff;
  private FDTD2D fdtd;

  private int NNx;
  private int NNy;
  private double magn = 1.0;

  double viewTheta = -15.0*3.14/180.0;
  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 int xpts[] = new int[2200];
  private int ypts[] = new int[2200];

  private int drawnCount = 0;


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

  DrawGraph2D() {
  }


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

  public void plotGraph2D(Graphics gg, FDTD2D fd, int dispMode, int thCount) {
    int fdtdCount,gb;

    gOff = gg;
    fdtd = fd;
    if (fd==null) return;
    NNx = fdtd.getNNx(); NNy = fdtd.getNNy();

    fdtdCount = fdtd.getCount();
    if (fdtdCount==drawnCount) return;
    drawnCount = fdtdCount;

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

    if (dispMode==0) { // Ez grid plot
      gridPlot(15.0*magn);

    } else if (dispMode==1) { // Ez density
      densityPlot(15.0*magn, 1);

    } else if (dispMode==2) { // Hxy field
      densityPlot(15.0*magn, 2);

    } else if (dispMode==3) { // Ez + Hxy
      densityPlot(15.0*magn, 3);

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

    gb = 500;
    gOff.setColor(Color.black);
    gOff.drawString("grid=800x800",gb,60);

    gOff.drawString("thread count",gb,80);
    gOff.drawString("FDTD = "+fdtd.getCount()+" ",gb,100);
    gOff.drawString("disp = "+thCount+" ",gb,120);

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

  }


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

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

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

    gbx = 60;
    gby =-10;

    for(j=0; j<NNy; j += 16) {
      for (i=0; i<NNx; i++) { 
        z = mag*fdtd.getEz(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;
      }
      for (i=0; i<NNx-2; i++) {
        if (fdtd.getSgm(i,j)>0.0) {
          gOff.setColor(Color.getHSBColor(colorOf(i,j),0.9f,0.5f));
        } else {
          gOff.setColor(Color.getHSBColor(colorOf(i,j),0.9f,0.9f));
        }
        gOff.drawLine(xpts[i],ypts[i],xpts[i+1],ypts[i+1]);
      }
    }

    for(i=0; i<NNx; i += 16) {
      for (j=0; j<NNy; j++) { 
        z = mag*fdtd.getEz(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;
      }
      for (j=0; j<NNy-2; j++) { 
        if (fdtd.getSgm(i,j)>0.0) {
          gOff.setColor(Color.getHSBColor(colorOf(i,j),0.9f,0.5f));
        } else {
          gOff.setColor(Color.getHSBColor(colorOf(i,j),0.9f,0.9f));
        }
        gOff.drawLine(xpts[j],ypts[j],xpts[j+1],ypts[j+1]);
      }
    }
  }

  private void densityPlot(double mag, int mode) {
    int i,j,gx,gy,ix,iy,ir,ic;
    double cs,cb,sc;

    gx = 30; gy = 50;
    ic = 4;
    sc = 0.4;
    ir = (int)(sc*ic+0.99999);

    for(i=0+ic/2; i<NNx-ic; i += ic) {
      for (j=0+ic/2; j<NNy-ic; j += ic) { 
        cs = 0.3+0.6/(30.0*fdtd.getSgm(i,j)+1.0); if (cs>0.999) cs = 0.999;
        if (cs<0.0) cs = 0.0;
        cb = 0.5+0.25*fdtd.getEz(i,j); if (cb>0.999) cb = 0.999;
        if (cb<0.0) cb = 0.0;
        if (mode==1 || mode==3) {
          gOff.setColor(Color.getHSBColor( colorOf(i,j), (float)cs, (float)cb ));
        } else if (mode==2) {
          gOff.setColor(Color.getHSBColor( colorOf(i,j), 0.3f, 0.9f ));
        }
        ix = (int)(gx+sc*i);
        iy = (int)(gy+sc*j);
        gOff.fillRect(ix,iy,ir,ir);
      }
    }

    if (mode==2 || mode==3) {
      for(i=0+ic; i<NNx-ic*2; i += ic*2) {
        for (j=0+ic; j<NNy-ic*2; j += ic*2) {
          if (fdtd.getHy(i,j)>0) {
            gOff.setColor(Color.red);
          } else {
            gOff.setColor(Color.blue);
          }
          ix = (int)(gx+sc*(i+ic));
          iy = (int)(gy+sc*(j+ic));
          gOff.drawLine(ix,iy,
            ix+(int)(mag*fdtd.getHx(i,j)),
            iy+(int)(mag*fdtd.getHy(i,j)) );
        }
      }
    }
  }

  private float colorOf(int i, int j) {
    float c;

    if (fdtd.getMetalQ(i,j)) return 0.12f;
    c = 0.33f;
    if (fdtd.getEps(i,j)>1.0) {
      if (fdtd.getMue(i,j)>1.0) {
        c = 0.85f;
      } else {
        c = 0.01f;
      }
    } else {
      if (fdtd.getMue(i,j)>1.0) c = 0.66f;
    }
    if (fdtd.getSgm(i,j)>0.0) c = c-0.05f;
    return c;
  }


  // ------------------------------------- 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 2D class ===========

class FDTD2D extends Thread {

  private float t = 0.0f;
  private float dt = 1.0f;

  private float omega = (float)(Math.PI/32.0);
  private float theta = 0.0f;

  private int NNx = 801;
  private int NNy = 801;

  private float Ez[][] = new float[NNx+1][NNy+1];
  private float Hx[][] = new float[NNx+1][NNy+1];
  private float Hy[][] = new float[NNx+1][NNy+1];

  private byte mat[][] = new byte[NNx+1][NNy+1];

  private float matTable[][] = {
    // eps, mue, sgm
    { 1.0f, 1.0f, 0.0f }, // 0-vacuum or air
    { 2.0f, 1.0f, 0.0f }, // 1-dielectric
    { 1.0f, 4.0f, 0.0f }, // 2-magnetic
    { 2.0f, 4.0f, 0.0f }, // 3-dielectric + magnetic
    { 1.0f, 1.0f, 0.05f }, // 4-registor
    { 1.0e30f, 1.0f, 0.0f }, // 5-metal
    { 1.0f, 1.0f, 0.3f } // 6-wall
  };

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

  FDTD2D() {
    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.0f;
    theta = 0.0f;

    clearField();
    setMaterial();
  }

  private void clearField() {
    int i,j;

    for(i=0;i<NNx+1;i++) {
      for(j=0;j<NNy+1;j++) {
       	Ez[i][j] = 0.0f;
       	Hx[i][j] = 0.0f;
       	Hy[i][j] = 0.0f;
      }
    }
  }

  private void setMaterial() {
    int i,j;

    for(i=0;i<NNx+1;i++) {
      for(j=0;j<NNy+1;j++) {
       	mat[i][j] = 0; // vacuum
      }
    }

    for(i=0;i<5;i++) {
      for(j=0;j<6;j++) {
        setPowder( i*80+300, j*100+50*(i%2)+100, 30);
      }
    }

  }

  private void setPowder(int i0, int j0, int ir) {
    int i,j,iBeg,iEnd,jBeg,jEnd;

    iBeg = i0-ir; if (iBeg<0) iBeg=0;
    iEnd = i0+ir; if (iEnd>=NNx) iEnd=NNx-1;
    jBeg = j0-ir; if (jBeg<0) jBeg=0;
    jEnd = j0+ir; if (jEnd>=NNy) jEnd=NNy-1;

    for(i=iBeg;i<=iEnd;i++) {
      for(j=jBeg;j<=jEnd;j++) {
        if ((i-i0)*(i-i0)+(j-j0)*(j-j0)<ir*ir) {
          mat[i][j] = 1; // dielectric ( eps = 2.0 )
        }
      }
    }
  }

  // ---------------------------------- 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(){

    t = t + 0.5f*dt;
    EzEvolution();
    EzGenerator();
    HxyEvolution();

    t = t + 0.5f*dt;
    EzEvolution();
    EzGenerator();
    HxyEvolution();
  }

  private void EzGenerator() {
    int j;
    float Ezt;

    theta += omega*0.5f*dt;
    Ezt = (float)Math.sin(theta);
    if (theta<3.1416f*2.0f) {
      for(j=0;j<NNy;j++) {
        Ez[0][j]=Ezt;
      }
    }
  }

  private void EzEvolution() {
    int i,j,mij;
    float sqrt2v2,d2,sgm,eps;

    sqrt2v2 = 0.707f*0.5f;
    d2 = dt*0.5f;

    for(i=0;i<NNx;i++) {
      for(j=0;j<NNy;j++) {
        mij = mat[i][j];
        if (matTable[mij][0]<1.0e20f) {
          sgm = matTable[mij][2]; eps = matTable[mij][0];
          Ez[i][j] += (-sgm*d2*Ez[i][j] + (sqrt2v2/eps)*((Hy[i+1][j]-Hy[i][j])-(Hx[i][j+1]-Hx[i][j])));
        }
      }
    }
  }

  private void HxyEvolution() {
    int i,j;
    float sqrt2v2,c;

    sqrt2v2 = 0.707f*0.5f;
    c = (float)((1.0-sqrt2v2)/(1.0+sqrt2v2));

    for(i=0; i<NNx; i++) {
      Hx[i][0] = c*Hx[i][0] + Hx[i][1];
    }
    for(i=0; i<NNx; i++) {
      Hx[i][NNy] = c*Hx[i][NNy] + Hx[i][NNy-1];
    }
    for(j=0; j<NNy; j++) {
      Hy[0][j] = c*Hy[0][j] + Hy[1][j];
    }
    for(j=0; j<NNy; j++) {
      Hy[NNx][j] = c*Hy[NNx][j] + Hy[NNx-1][j];
    }

    for(i=0;i<NNx;i++) {
      for(j=1;j<NNy;j++) {
       	Hx[i][j] -= sqrt2v2/matTable[mat[i][j]][1]*(Ez[i][j]-Ez[i][j-1]);
      }
    }

    for(i=1;i<NNx;i++) {
      for(j=0;j<NNy;j++) {
       	Hy[i][j] += sqrt2v2/matTable[mat[i][j]][1]*(Ez[i][j]-Ez[i-1][j]);
      }
    }

    for(i=0;i<NNx;i++) {
      Hx[i][0] -= c*Hx[i][1];
    }
    for(i=0;i<NNx;i++) {
      Hx[i][NNy] -= c*Hx[i][NNy-1];
    }
    for(j=0;j<NNy;j++) {
      Hy[0][j] -= c*Hy[1][j];
    }
    for(j=0;j<NNy;j++) {
      Hy[NNx][j] -= c*Hy[NNx-1][j];
    }
  }


  // ------------------------------------- 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 = (float)(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 double getEz(int i, int j) { return Ez[i][j]; }
  public double getHx(int i, int j) { return Hx[i][j]; }
  public double getHy(int i, int j) { return Hy[i][j]; }

  public double getSgm(int i, int j) { return matTable[mat[i][j]][2]; }
  public double getEps(int i, int j) { return matTable[mat[i][j]][0]; }
  public double getMue(int i, int j) { return matTable[mat[i][j]][1]; }
  public boolean getMetalQ(int i, int j) { return (matTable[mat[i][j]][0]>1.0e20f); }
}


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