/* lattice Boltzmann method     */
/* coded by Ikeuchi Mitsuru     */
/* ver 0.0.1   2003.12.20       */
/* ver 0.0.2   2003.12.23  bug(setTemp) fixed */
/* ver 0.0.3   2004.01.23  bugs in scattering() fixed */
/* ver 0.0.4   2004.02.24  improved pannel layout */
/* ver 0.0.5   2004.03.15  bug(offPaint) fixed */

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

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

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

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

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 50;


  double omega = 0.50;
  double scale = 1.0;
  int started = 1;
  int dispMode = 0;

  int dgX, dgY, dgXb,dgYb;

  int cx = (int)(50.0/scale);
  int cy = (int)(50.0/scale);
  int cz = (int)(50.0/scale);

  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 dt = 1.0;
  double dx = 1.0;
  double dy = 1.0;

  int tpp = 0;
  int tnx = 1;

  int NNx = 100;
  int NNy = 100;
  double ccx[] = { 0.0, 1.0, 1.0, 0.0,-1.0,-1.0,-1.0, 0.0, 1.0 };
  double ccy[] = { 0.0, 0.0,-1.0,-1.0,-1.0, 0.0, 1.0, 1.0, 1.0 };
  double cc2[] = { 0.0, 1.0, 2.0, 1.0, 2.0, 1.0, 2.0, 1.0, 2.0 };

  int kind[][] = new int[NNx][NNy];
  double v9[][][][] = new double[2][NNx][NNy][9];
  double feq[] = new double[9];
  double dfn[] = new double[9];
  double rho[][] = new double[NNx][NNy];
  double ux[][] = new double[NNx][NNy];
  double uy[][] = new double[NNx][NNy];
  double flow[][] = new double[NNx][NNy];
  double temp[][] = new double[NNx][NNy];

  int NN = 100+2;
  int xpts[] = new int[NN];
  int ypts[] = new int[NN];


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

  public void init() {

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

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

    cvw = new Choice();
    cvw.addItem("flow"); cvw.addItem("dens");
    cvw.addItem("vx-vy"); cvw.addItem("flowE");
    cvw.addItemListener(this);

    scr= new Scrollbar(Scrollbar.HORIZONTAL,50,10,10,110);
    scr.addAdjustmentListener(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);

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

    addMouseListener(this);  
    addMouseMotionListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,6,5,0));
    pnl.add(cvw);
    pnl.add(bt_reset);
    pnl.add(bt_start);
    pnl.add(bt_stop);
    pnl.add(new Label("  "));
    pnl.add(scr);
    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 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 itemStateChanged(ItemEvent ev){
    if (ev.getSource() == cvw) {
      dispMode = cvw.getSelectedIndex();
      repaint(); 
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if(ev.getSource() == scr) { 
      omega = 0.01*(double)(scr.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 i,gb;

    vDistribution();

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

    if (dispMode==0) {

      densPlot(rho, 0.66, 0.25);
      flowPlot(40.0);
      gOff.setColor(Color.black);
      gOff.drawString("density",100,265);

    } else if (dispMode==1) {

      gOff.setColor(Color.green);
      linePlot(rho,20.0);
      gOff.drawString("density",30,70);

    } else if (dispMode==2) {

      gOff.setColor(Color.blue);
      xLinePlot(ux,50.0);
      gOff.drawString("ux",30,70);
      gOff.setColor(Color.red);
      yLinePlot(uy,50.0);
      gOff.drawString("uy",60,70);

    } else if (dispMode==3) {

      setFlowEnergy();
      densPlot(flow, 0.33, 6.0);
      flowPlot(40.0);
      gOff.setColor(Color.black);
      gOff.drawString("flow energy",100,265);

    }


    gb = 400;
    gOff.setColor(Color.black);
    gOff.drawRect(gb,50,90,200);
    gOff.setColor(Color.pink);
    for (i=0; i<9; i++) {
      gOff.fillRect(gb+10*i+1,250-(int)dfn[i]/24, 8, (int)dfn[i]/24);
      gOff.drawString(""+i,gb+10*i,265);
    }

    gOff.setColor(Color.black);
    gOff.drawString("view",630/6*0+10,40);
    gOff.drawString("omega="+(int)(omega*100+0.5)/100.0+"",630/6*5+10,40);
    gb = 500;
    gOff.setColor(Color.black);
    gOff.drawString("t="+(int)(t+0.5)+"(times)",500,70);
    for (i=0; i<9; i++) {
      gOff.drawString("f"+i+"="+(int)(dfn[i]*10000/(NNx*NNy))+"/10000",500,110+i*15);
    }

  }

/*-----------------------------------  plot utility  ----------------*/

  void densPlot(double fn[][], double col, double sc) {
    int i,j,ix,iy;
    int gbx,gby;
    double s;

    gbx = 30;
    gby = 50;
    for(i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        ix = i*2+gbx;
        iy = j*2+gby;
        s = sc*fn[i][j];
        gOff.setColor(Color.getHSBColor((float)(col-s),0.95f,0.95f));
        gOff.fillRect(ix,iy,2,2);
      }
    }
  }

  void flowPlot(double mag) {
    int i,j,ix,iy,ix2,iy2;
    int gbx,gby;

    gbx = 30;
    gby = 50;
    for(i=5/2; i<NNx; i+=5) {
      for (j=5/2; j<NNy; j+=5) {
        ix = i*2+gbx;
        iy = j*2+gby;
        ix2 = ix + (int)(mag*ux[i][j]);
        iy2 = iy + (int)(mag*uy[i][j]);

        gOff.setColor(Color.red);
        gOff.drawLine(ix,iy, ix2,iy2);
      }
    }
  }


  void linePlot(double fn[][], double mag) {
    xLinePlot(fn,mag);
    yLinePlot(fn,mag);
  }

  void xLinePlot(double fn[][], double mag) {
    int i,j;
    int gbx,gby;
    double z,ppx,ppy,ppz;

    gbx = 50;
    gby = 70;

    for(j=0; j<NNy; j += 5) {
      for (i=0; i<NNx; i++) {
        z = mag*fn[i][j]+cz;
        ppy = cosFi*(j-cy)+sinFi*(z-cz) + cy; 
        ppz = -sinFi*(j-cy)+cosFi*(z-cz) + cz;
        ppx = cosTh*(i-cx)+sinTh*(ppz-cz) + cx; 
        ppz = -sinTh*(i-cx)+cosTh*(ppz-cz) + cz;
        xpts[i] = (int)(scale*2*ppx)+gbx;
        ypts[i] = (int)(scale*2*ppy)+gby;
      }
      for (i=0; i<NNx-2; i++) {
        gOff.drawLine(xpts[i],ypts[i],xpts[i+1],ypts[i+1]);
      }
    }
  }

 void yLinePlot(double fn[][], double mag) {
    int i,j;
    int gbx,gby;
    double z,ppx,ppy,ppz;

    gbx = 50;
    gby = 70;

    for(i=0; i<NNx; i += 5) {
      for (j=0; j<NNy; j++) { 
        z = mag*fn[i][j]+cz;
        ppy = cosFi*(j-cy)+sinFi*(z-cz) + cy; 
        ppz = -sinFi*(j-cy)+cosFi*(z-cz) + cz;
        ppx = cosTh*(i-cx)+sinTh*(ppz-cz) + cx; 
        ppz = -sinTh*(i-cx)+cosTh*(ppz-cz) + cz;
        xpts[j] = (int)(scale*2*ppx)+gbx;
        ypts[j] = (int)(scale*2*ppy)+gby;
      }
      for (j=0; j<NNy-2; j++) { 
        gOff.drawLine(xpts[j],ypts[j],xpts[j+1],ypts[j+1]);
      }
    }
  }

/*-------------------------------------  statictics  ----------*/


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

    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        for (k=0;k<9;k++) {
          dfn[k] = 0.0;
        }
      }
    }

    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        for (k=0;k<9;k++) {
          dfn[k] += v9[tpp][i][j][k];
        }
      }
    }
  }


  void setFlowEnergy() {
    int i,j;

    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        flow[i][j] = 0.5*rho[i][j]*(ux[i][j]*ux[i][j]+uy[i][j]*uy[i][j]);
      }
    }
  }



  void setTemp() {
    int i,j,k;
    double s,vx,vy;

    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        s = 0.0;
        for (k=0;k<9;k++) {
          vx = ccx[k]-ux[i][j];
          vy = ccy[k]-uy[i][j];
          s += 0.5*v9[tpp][i][j][k]*(vx*vx+vy*vy);
        }
        if (rho[i][j]>0.0) {
          temp[i][j] = s/rho[i][j];
        } else {
          temp[i][j] = 0.0;
        }
      }
    }
  }



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

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

    t = 0.0;

    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        kind[i][j]=0;
        rho[i][j] = 0.0;
        ux[i][j] = 0.0;
        uy[i][j] = 0.0;

        if ((i-30)*(i-30)+(j-50)*(j-50)<20*20) {
          v9[tpp][i][j][0] = 0.14*3.0;
          v9[tpp][i][j][1] = 0.06*3.0;
          v9[tpp][i][j][2] = 0.02*3.0;
          v9[tpp][i][j][3] = 0.06*3.0;
          v9[tpp][i][j][4] = 0.02*3.0;
          v9[tpp][i][j][5] = 0.06*3.0;
          v9[tpp][i][j][6] = 0.02*3.0;
          v9[tpp][i][j][7] = 0.06*3.0;
          v9[tpp][i][j][8] = 0.02*3.0;
        } else {
          v9[tpp][i][j][0] = 0.14;
          v9[tpp][i][j][1] = 0.06;
          v9[tpp][i][j][2] = 0.02;
          v9[tpp][i][j][3] = 0.06;
          v9[tpp][i][j][4] = 0.02;
          v9[tpp][i][j][5] = 0.06;
          v9[tpp][i][j][6] = 0.02;
          v9[tpp][i][j][7] = 0.06;
          v9[tpp][i][j][8] = 0.02;
        }

      }
    }

  }



/* ----------------------  lattice Boltzmann time evolution  ----- */

  void timeEvolution() {
    if (started==1) {

      freeFlight();
      scattering(omega);

      t = t + dt;
      tpp = tnx;
      tnx = (tnx+1)%2;

    }
  }

  void freeFlight() {
    int i,j,ip,im,jp,jm;

    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        ip = (i+1)%NNx; im = (i-1+NNx)%NNx;
        jp = (j+1)%NNy; jm = (j-1+NNy)%NNy;
        v9[tnx][i][j][0] = v9[tpp][i][j][0];
        v9[tnx][i][j][1] = v9[tpp][im][j][1];
        v9[tnx][i][j][2] = v9[tpp][im][jp][2];
        v9[tnx][i][j][3] = v9[tpp][i][jp][3];
        v9[tnx][i][j][4] = v9[tpp][ip][jp][4];
        v9[tnx][i][j][5] = v9[tpp][ip][j][5];
        v9[tnx][i][j][6] = v9[tpp][ip][jm][6];
        v9[tnx][i][j][7] = v9[tpp][i][jm][7];
        v9[tnx][i][j][8] = v9[tpp][im][jm][8];
      }
    }

  }

  void scattering(double omg) {
    int i,j,k;
    double m;

    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        if (kind[i][j]==0) {
          setfeq(i,j);
          for (k=0;k<9;k++) {
            v9[tnx][i][j][k] -= omg*(v9[tnx][i][j][k]-feq[k]);
          }
        } else if (kind[i][j]==1) {
          v9[tnx][i][j][0] = 0.0;
          m = v9[tnx][i][j][1]; 
          v9[tnx][i][j][1]=v9[tnx][i][j][5];
          v9[tnx][i][j][5]=m;
          m = v9[tnx][i][j][2]; 
          v9[tnx][i][j][2]=v9[tnx][i][j][6];
          v9[tnx][i][j][6]=m;
          m = v9[tnx][i][j][3]; 
          v9[tnx][i][j][3]=v9[tnx][i][j][7];
          v9[tnx][i][j][7]=m;
          m = v9[tnx][i][j][4]; 
          v9[tnx][i][j][4]=v9[tnx][i][j][8];
          v9[tnx][i][j][8]=m;
        }
      }
    }

  }

  void setfeq(int i, int j) {
    int k;
    double rh,vx,vy,v2,c1v,c2v,c3v,c4v,c5v,c6v,c7v,c8v;
    double w0,w1,w2;
/*
    double ccx[] = { 0.0, 1.0, 1.0, 0.0,-1.0,-1.0,-1.0, 0.0, 1.0 };
    double ccy[] = { 0.0, 0.0,-1.0,-1.0,-1.0, 0.0, 1.0, 1.0, 1.0 };
*/

    w0 = 4.0/9.0;
    w1 = 1.0/9.0;
    w2 = 1.0/36.0;

    rh = density(i,j);
    rho[i][j] = rh;
    if (rh>0) {
      vx = rhvx(i,j)/rh;
      vy = rhvy(i,j)/rh;
      ux[i][j] = vx;
      uy[i][j] = vy;
      v2 = vx*vx+vy*vy;

      c1v = vx;
      c2v = vx-vy;
      c3v = -vy;
      c4v = -vx-vy;
      c5v = -vx;
      c6v = -vx+vy;
      c7v = vy;
      c8v = vx+vy;

      feq[0] = w0*rh*(1.0-1.5*v2);
      feq[1] = w1*rh*(1.0+3.0*c1v+4.5*c1v*c1v-1.5*v2);
      feq[2] = w2*rh*(1.0+3.0*c2v+4.5*c2v*c2v-1.5*v2);
      feq[3] = w1*rh*(1.0+3.0*c3v+4.5*c3v*c3v-1.5*v2);
      feq[4] = w2*rh*(1.0+3.0*c4v+4.5*c4v*c4v-1.5*v2);
      feq[5] = w1*rh*(1.0+3.0*c5v+4.5*c5v*c5v-1.5*v2);
      feq[6] = w2*rh*(1.0+3.0*c6v+4.5*c6v*c6v-1.5*v2);
      feq[7] = w1*rh*(1.0+3.0*c7v+4.5*c7v*c7v-1.5*v2);
      feq[8] = w2*rh*(1.0+3.0*c8v+4.5*c8v*c8v-1.5*v2);

    } else {
      feq[0] = 0.0;
      feq[1] = 0.0;
      feq[2] = 0.0;
      feq[3] = 0.0;
      feq[4] = 0.0;
      feq[5] = 0.0;
      feq[6] = 0.0;
      feq[7] = 0.0;
      feq[8] = 0.0;
    }

  }

  double rhvx(int i, int j) {
    int k;
    double s;

    s = 0.0;
    for (k=0;k<9;k++) {
      s += ccx[k]*v9[tnx][i][j][k];
    }
    return( s );
  }

  double rhvy(int i, int j) {
    int k;
    double s;

    s = 0.0;
    for (k=0;k<9;k++) {
      s += ccy[k]*v9[tnx][i][j][k];
    }
    return( s );
  }

  double density(int i, int j) {
    int k;
    double s;

    s = 0.0;
    for (k=0;k<9;k++) {
      s += v9[tnx][i][j][k];
    }
    return( s );
  }


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

}