/* porous media flow - lattice Boltzmann method 2D9V */
/* coded by Ikeuchi Mitsuru   */
/* ver 0.0.1   2005.02.26     */

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

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

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

  Button bt_reset, bt_start, bt_stop;
  Choice ch_view;
  Scrollbar sc_omega, sc_fan;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 50;


  double omega = 1.7;
  double scale = 1.0;
  double fan = 0.6;
  int started = 1;
  int dispMode = 0;

  int dgX, dgY, dgXb,dgYb;

  int cx = (int)(300.0/2.0/scale);
  int cy = (int)(300.0/2.0/scale);
  int cz = (int)(300.0/2.0/scale);

  double pai = 3.1415926536;
  double viewTheta = -15.0*pai/180.0;
  double viewFai = -72.0*pai/180.0;
  double dtheta = 0.0*pai/180.0;

  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 = 600;
  int NNy = 300;
  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 = NNx+2;
  int xpts[] = new int[NN];
  int ypts[] = new int[NN];


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

  public void init() {

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

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

    sc_omega= new Scrollbar(Scrollbar.HORIZONTAL,170,10,1,209);
    sc_omega.addAdjustmentListener(this);

    sc_fan= new Scrollbar(Scrollbar.HORIZONTAL,60,10,0,110);
    sc_fan.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);

    ch_view = new Choice();
    ch_view.addItem("flow"); ch_view.addItem("dens");
    ch_view.addItem("vx-vy");
    ch_view.addItemListener(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(sc_omega);
    pnl.add(sc_fan);
    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 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() == ch_view) {
      dispMode = ch_view.getSelectedIndex();
      repaint(); 
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if(ev.getSource() == sc_omega) { 
      omega = 0.01*(double)(sc_omega.getValue());
    } else if(ev.getSource() == sc_fan) { 
      fan = 0.01*(double)(sc_fan.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, gbx, gby;
    double nue,vmean;

    vDistribution();

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

    if (dispMode==0) {

      flowPlot(200.0);

    } else if (dispMode==1) {

      densPlot(rho, 0.66, 0.25);
      flowPlot(200.0);
	  
    } else if (dispMode==2) {

      gOff.setColor(Color.blue);
      xLinePlot(ux,200.0);
      gOff.drawString("ux",100,360);
      gOff.setColor(Color.red);
      yLinePlot(uy,200.0);
      gOff.drawString("uy",150,360);

    } else if (dispMode==3) {

      ;

    } else if (dispMode==4) {

      ;

    }


    gOff.setColor(Color.black);
    gOff.drawString("t="+(int)(t+0.5)+" (times)",630/6*0+10,40);
    gOff.drawString("omega="+(int)(omega*100+0.5)/100.0+" ",630/6*3+10,40);
    gOff.drawString("fan="+(int)(fan*100+0.5)+"%",630/6*4+10,40);
    gOff.drawString("view",630/6*5+10,40);
    nue = (1.0/omega-0.5)/3.0;
	vmean = vFan();
	gOff.drawString("nue=(1/omega-1/2)/3="+(int)(nue*1000+0.5)/1000.0+" ",630/6*0+10,60);
    gOff.drawString("Re=v*D/nue="+(int)(vmean*1000+0.5)/1000.0+"*40/"+(int)(nue*1000+0.5)/1000.0+" = "+(int)(vmean*40.0/nue+0.5)+" ",630/6*3+10,60);

/*
    gbx = 460; gby = 60;
    gOff.setColor(Color.black);
    gOff.drawRect(gbx,gby,90,200);
    gOff.setColor(Color.pink);
    for (i=0; i<9; i++) {
      gOff.fillRect(gbx+10*i+1,gby+200-(int)dfn[i]/24, 8, (int)dfn[i]/24);
    }
    gOff.setColor(Color.black);
    for (i=0; i<9; i++) {
	  gOff.drawString(""+i,gbx+10*i,gby+215);
      gOff.drawString("f"+i+"="+(int)(dfn[i]*10000/(NNx*NNy))+"",gbx+100,gby+40+i*15);
    }
    gOff.drawString("   /10000",gbx+10*i,gby+40+9*15);
*/
  }

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

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

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

  }

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

    gbx = 10;
    gby = 80;

	gOff.setColor(Color.blue);
    for(i=100; i<450; i++) {
      for (j=0; j<NNy; j++) {
	    if (kind[i][j]==1) {
          ix = i*1+gbx;
          iy = j*1+gby;
          gOff.fillRect(ix,iy,1,1);
		}
      }
    }

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

        if (ix>ix2) {
          gOff.setColor(Color.blue);
		} else {
          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 = 30;
    gby = 100;

    for(j=0; j<NNy; j += 5) {
      for (i=0; i<NNx-20; 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*1*ppx)+gbx;
        ypts[i] = (int)(scale*1*ppy)+gby;
      }
      for (i=0; i<NNx-2-20; 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 = 30;
    gby = 100;

    for(i=0; i<NNx-20; 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*1*ppx)+gbx;
        ypts[j] = (int)(scale*1*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;
        }
      }
    }
  }

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

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

/*-------------------------- 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;

        v9[tpp][i][j][0] = 0.24;
        v9[tpp][i][j][1] = 0.06+0.03*fan;
        v9[tpp][i][j][2] = 0.02+0.01*fan;
        v9[tpp][i][j][3] = 0.06;
        v9[tpp][i][j][4] = 0.02-0.01*fan;
        v9[tpp][i][j][5] = 0.06-0.03*fan;
        v9[tpp][i][j][6] = 0.02-0.01*fan;
        v9[tpp][i][j][7] = 0.06;
        v9[tpp][i][j][8] = 0.02+0.01*fan;
      }
    }

    for (i=150;i<400;i+=51) {
      for (j=0+25*(i%2);j<NNy+25;j+=60) {
        putBall(i, j, 20);
      }
    }
  }

  void putBall(int ix, int iy, int ir) {
    int i,j,k,i0,i1,j0,j1;

    i0 = ix-ir; if (i0<0) i0 = 0;
	i1 = ix+ir; if (i1>=NNx) i1 = NNx-1;
    j0 = iy-ir; if (j0<0) j0 = 0;
	j1 = iy+ir; if (j1>=NNy) j1 = NNy-1;
    for (i=i0;i<=i1;i++) {
      for (j=j0;j<=j1;j++) {
		if ((i-ix)*(i-ix)+(j-iy)*(j-iy)<=ir*ir  ) {
          kind[i][j] = 1;
          for (k=0;k<9;k++) {
            v9[tpp][i][j][k] = 0.0;
          }
        }
	  }
	}
  }

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

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

      freeFlight();
      scattering(omega);
      setFan();

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

  void setFan() {
    int i,j;

    i = 0;
    for (j=0;j<NNy;j++) {
        v9[tnx][i][j][0] = 0.24;
        v9[tnx][i][j][1] = 0.06+0.03*fan;
        v9[tnx][i][j][2] = 0.02+0.01*fan;
        v9[tnx][i][j][3] = 0.06;
        v9[tnx][i][j][4] = 0.02-0.01*fan;
        v9[tnx][i][j][5] = 0.06-0.03*fan;
        v9[tnx][i][j][6] = 0.02-0.01*fan;
        v9[tnx][i][j][7] = 0.06;
        v9[tnx][i][j][8] = 0.02+0.01*fan;
    }
  }

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

}