/* expansion LBGK3D27V - D3Q27 model */
/* coded by Ikeuchi Mitsuru           */
/* ver 0.0.1   2004.01.24             */

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

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

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

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

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 50;

  double omega = 0.5;
  double fan = 0.5;
  int started = 1;
  int dispMode = 0;
  int viewk = 20;

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

  int tpp = 0;
  int tnx = 1;

  int NNx = 40;
  int NNy = 40;
  int NNz = 40;
  double ccx[] = { 0.0, 1.0,-1.0, 0.0, 0.0, 0.0, 0.0, 1.0,-1.0, 1.0,-1.0, 1.0,-1.0,  1.0,-1.0, 0.0, 0.0, 0.0, 0.0, 1.0,-1.0, 1.0,-1.0, 1.0,-1.0,  1.0,-1.0 };
  double ccy[] = { 0.0, 0.0, 0.0, 1.0,-1.0, 0.0, 0.0, 1.0,-1.0,-1.0, 1.0, 0.0, 0.0,  0.0, 0.0, 1.0,-1.0,  1.0,-1.0, 1.0, 1.0,-1.0,-1.0, 1.0, 1.0, -1.0,-1.0 };
  double ccz[] = { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0,-1.0, 0.0, 0.0, 0.0, 0.0, 1.0,-1.0, -1.0, 1.0, 1.0,-1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0,-1.0,-1.0, -1.0,-1.0 };
  double cc2[] = { 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0 };

  int kind[][][] = new int[NNx][NNy][NNz];
  double fv[][][][][] = new double[2][NNx][NNy][NNz][27];
  double feq[] = new double[27];
  double dfn[] = new double[27];
  double rho[][][] = new double[NNx][NNy][NNz];
  double ux[][][] = new double[NNx][NNy][NNz];
  double uy[][][] = new double[NNx][NNy][NNz];
  double uz[][][] = new double[NNx][NNy][NNz];
  double flow[][][] = new double[NNx][NNy][NNz];
  double temp[][][] = new double[NNx][NNy][NNz];

  int xpts[] = new int[NNx+NNy+2];
  int ypts[] = new int[NNx+NNy+2];


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

  public void init() {

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

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

    scr= new Scrollbar(Scrollbar.HORIZONTAL,50,10,10,200);
    add(new Label("omega=")); add(scr);
    scr.addAdjustmentListener(this);

    scf= new Scrollbar(Scrollbar.HORIZONTAL,50,10,0,110);
    add(new Label("fan=")); add(scf);
    scf.addAdjustmentListener(this);

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

    bt_start= new Button("start"); add(bt_start);
    bt_start.addActionListener(this); 

    bt_stop = new Button("stop");  add(bt_stop);
    bt_stop.addActionListener(this);

    addMouseListener(this);  
    addMouseMotionListener(this);

    cvw = new Choice();
    cvw.addItem("flow"); cvw.addItem("depth");
    cvw.addItem("mesh");
    add(new Label("view=")); add(cvw);
    cvw.addItemListener(this);

    sck= new Scrollbar(Scrollbar.HORIZONTAL,20,10,0,49);
    add(new Label("depth=")); add(sck);
    sck.addAdjustmentListener(this);

    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){
    int m;

    if (ev.getSource() == cvw) {
      dispMode = cvw.getSelectedIndex();
      repaint(); 
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == scr) { 
      omega = 0.01*(double)(scr.getValue());
    } else if (ev.getSource() == scf) { 
      fan = 0.01*(double)(scf.getValue());
    } else if (ev.getSource() == sck) { 
      viewk = sck.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) {

      uPlot();

    } else if (dispMode==1) {

      rhouPlot(0.66, viewk);
      gOff.drawString("rho,ux,uy, at z="+viewk+"/"+NNz+" ",30,280);
      gb = 430;
      gOff.setColor(Color.black);
      gOff.drawRect(gb,50,80,200);
      gOff.setColor(Color.pink);
      for (i=0; i<27; i++) {
        gOff.fillRect(gb+3*i+1,250-(int)dfn[i]/200, 2, (int)dfn[i]/200);
      }

    } else if (dispMode==2) {

      if (viewk != 0) {
        gOff.setColor(Color.gray);
        linePlot(rho, 5.0, 0);
      }
      gOff.setColor(Color.green);
      linePlot(rho, 5.0, viewk);
      if (viewk != NNz-1) {
        gOff.setColor(Color.gray);
        linePlot(rho, 5.0, NNz-1);
      }
      gOff.drawString("rho, at z="+viewk+"/"+NNz+" ",30,280);

    }


    gb = 520;
    gOff.setColor(Color.black);
    gOff.drawString("omega="+(int)(omega*100+0.5)/100.0+"(1/t)",gb,50);
    gOff.drawString("fan="+(int)(fan*100+0.5)+"(%)",gb,65);
    gOff.drawString("t="+(int)(t+0.5)+"(times)",gb,80);
/*
    gOff.setColor(Color.pink);
    for (i=0; i<27; i++) {
      gOff.drawString("f"+i+"="+(int)(dfn[i]*10000/(NNx*NNy*NNz))+"/10000",gb,95+i*11);
    }
*/
  }

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

  void uPlot() {
    int i,j,k,ix,iy,ir,ix2,iy2;
    int gbx,gby;
    double cxx,cyy,czz;
    double z,ppx,ppy,ppz;
    double mag,mm;

    gbx = 40;
    gby = 60;
    mag = 5.0;
    mm = 40.0;
    cxx = mag*(NNx/2.0); 
    cyy = mag*(NNy/2.0);
    czz = mag*(NNz/2.0);

    for (i=0; i<NNx; i+=2) {
      for(j=0; j<NNy; j +=2) {
        for (k=0; k<NNz; k+=2) {

          z = k*mag;
          ppy = cosFi*(j*mag-cyy)+sinFi*(z-czz) + cyy; 
          ppz = -sinFi*(j*mag-cyy)+cosFi*(z-czz) + czz;
          ppx = cosTh*(i*mag-cxx)+sinTh*(ppz-czz) + cxx; 
          ppz = -sinTh*(i*mag-cxx)+cosTh*(ppz-czz) + czz;
          ix = (int)(ppx)+gbx;
          iy = (int)(ppy)+gby;
          ir = (int)(2.0*rho[i][j][k]-0.01);
          gOff.setColor(Color.getHSBColor(0.33f,0.95f,0.95f));

          if (kind[i][j][k]==1) {
            ir = 10;
            gOff.setColor(Color.blue);
          }
          gOff.fillOval(ix-ir/2,iy-ir/2,ir,ir);

          if (i%4==0 && j%4==0 && k%4==0) {
            z = (k+mm*uz[i][j][k])*mag;
            ppy = cosFi*((j+mm*uy[i][j][k])*mag-cyy)+sinFi*(z-czz) + cyy; 
            ppz = -sinFi*((j+mm*uy[i][j][k])*mag-cyy)+cosFi*(z-czz) + czz;
            ppx = cosTh*((i+mm*ux[i][j][k])*mag-cxx)+sinTh*(ppz-czz) + cxx; 
            ppz = -sinTh*((i+mm*ux[i][j][k])*mag-cxx)+cosTh*(ppz-czz) + czz;

            ix2 = (int)(ppx)+gbx;
            iy2 = (int)(ppy)+gby;

            gOff.setColor(Color.getHSBColor((float)(0.01),0.95f,(float)(0.8*j/NNy+0.18)));
            gOff.drawLine(ix,iy, ix2,iy2);
          }

        }
      }
    }

  }

  void rhouPlot(double col, int kk) {
    int i,j,ix,iy,ix2,iy2;
    int gbx,gby;
    double s;

    gbx = 20;
    gby = 50;
    for(i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        s = rho[i][j][kk]-1.0;
        s = 0.33 - 1.0*s;
        while (s<0.0) {
          s += 1.0;
        }
        ix = i*5+gbx;
        iy = j*5+gby;
        gOff.setColor(Color.getHSBColor((float)(s),0.95f,0.95f));
        if (kind[i][j][kk]==1) {
            gOff.setColor(Color.blue);
        }
        gOff.fillRect(ix-2,iy-2,5,5);
      }
    }

    for(i=0; i<NNx; i+=2) {
      for (j=0; j<NNy; j+=2) {
        ix = i*5+gbx;
        iy = j*5+gby;
        ix2 = ix + (int)(100.0*ux[i][j][kk]);
        iy2 = iy + (int)(100.0*uy[i][j][kk]);
        gOff.setColor(Color.red);
        gOff.drawLine(ix,iy, ix2,iy2);
      }
    }

  }


  void linePlot(double fn[][][], double hight, int kk) {
    xLinePlot(fn,hight,kk);
    yLinePlot(fn,hight,kk);
  }

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

    mag = 5.0;
    cx = mag*(NNx/2.0); 
    cy = mag*(NNy/2.0);
    cz = mag*(NNz/2.0);

    gbx = 50;
    gby = 30;

    for(j=0; j<NNy; j += 2) {
      for (i=0; i<NNx; i++) {
        z = hight*20.0*(fn[i][j][kk]-1.0)+4*kk;
        if (kind[i][j][kk]==1) {
          z = 4*kk;
        }
        ppy = cosFi*(mag*j-cy)+sinFi*(z-cz) + cy; 
        ppz = -sinFi*(mag*j-cy)+cosFi*(z-cz) + cz;
        ppx = cosTh*(mag*i-cx)+sinTh*(ppz-cz) + cx; 
        ppz = -sinTh*(mag*i-cx)+cosTh*(ppz-cz) + cz;
        xpts[i] = (int)(ppx)+gbx;
        ypts[i] = (int)(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 hight, int kk) {
    int i,j;
    int gbx,gby;
    double mag,cx,cy,cz,z,ppx,ppy,ppz;

    mag = 5.0;
    cx = mag*(NNx/2.0); 
    cy = mag*(NNy/2.0);
    cz = mag*(NNz/2.0);

    gbx = 50;
    gby = 30;

    for(i=0; i<NNx; i += 2) {
      for (j=0; j<NNy; j++) { 
        z = hight*20.0*(fn[i][j][kk]-1.0)+4*kk;
        if (kind[i][j][kk]==1) {
          z = 4*kk;
        }
        ppy = cosFi*(mag*j-cy)+sinFi*(z-cz) + cy; 
        ppz = -sinFi*(mag*j-cy)+cosFi*(z-cz) + cz;
        ppx = cosTh*(mag*i-cx)+sinTh*(ppz-cz) + cx; 
        ppz = -sinTh*(mag*i-cx)+cosTh*(ppz-cz) + cz;
        xpts[j] = (int)(ppx)+gbx;
        ypts[j] = (int)(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,iv;

    for (iv=0;iv<27;iv++) {
      dfn[iv] = 0.0;
    }
    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        for (k=0;k<NNz;k++) {
          for (iv=0;iv<27;iv++) {
            dfn[iv] += fv[tpp][i][j][k][iv];
          }
        }
      }
    }
  }


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

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


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

    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        for (k=0;k<NNz;k++) {
          s = 0.0;
          for (iv=0;iv<27;iv++) {
            vx = ccx[iv]-ux[i][j][k];
            vy = ccy[iv]-uy[i][j][k];
            vz = ccz[iv]-uz[i][j][k];
            s += 0.5*fv[tpp][i][j][k][iv]*(vx*vx+vy*vy+vz*vz);
          }
          if (rho[i][j][k]>0.0) {
            temp[i][j][k] = s/rho[i][j][k];
          } else {
            temp[i][j][k] = 0.0;
          }
        }
      }
    }
  }


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

  void setInitialCondition() {
    int i,j,k,iv;
    double w0,w1,w2,w3;

    w0 = 8.0/27.0;
    w1 = 2.0/27.0;
    w2 = 1.0/54.0;
    w3 = 1.0/216.0;

    t = 0.0;

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

          fv[tpp][i][j][k][0] = w0;

          fv[tpp][i][j][k][1] = w1;
          fv[tpp][i][j][k][2] = w1;
          fv[tpp][i][j][k][3] = w1;
          fv[tpp][i][j][k][4] = w1;
          fv[tpp][i][j][k][5] = w1;
          fv[tpp][i][j][k][6] = w1;

          fv[tpp][i][j][k][7] = w2;
          fv[tpp][i][j][k][8] = w2;
          fv[tpp][i][j][k][9] = w2;
          fv[tpp][i][j][k][10] = w2;
          fv[tpp][i][j][k][11] = w2;
          fv[tpp][i][j][k][12] = w2;
          fv[tpp][i][j][k][13] = w2;
          fv[tpp][i][j][k][14] = w2;
          fv[tpp][i][j][k][15] = w2;
          fv[tpp][i][j][k][16] = w2;
          fv[tpp][i][j][k][17] = w2;
          fv[tpp][i][j][k][18] = w2;

          fv[tpp][i][j][k][19] = w3;
          fv[tpp][i][j][k][20] = w3;
          fv[tpp][i][j][k][21] = w3;
          fv[tpp][i][j][k][22] = w3;
          fv[tpp][i][j][k][23] = w3;
          fv[tpp][i][j][k][24] = w3;
          fv[tpp][i][j][k][25] = w3;
          fv[tpp][i][j][k][26] = w3;

          if ((i-15)*(i-15)+(j-15)*(j-15)+(k-20)*(k-20)<10*10) {
            for (iv=0;iv<27;iv++) {
              fv[tpp][i][j][k][iv] = fv[tpp][i][j][k][iv]*1.5;
            }
          }
        }
      }
    }

  }



/* ----------------------  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,k,ip,im,jp,jm,kp,km;

    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        for (k=0;k<NNz;k++) {
          ip = (i+1)%NNx; im = (i-1+NNx)%NNx;
          jp = (j+1)%NNy; jm = (j-1+NNy)%NNy;
          kp = (k+1)%NNz; km = (k-1+NNz)%NNz;
          fv[tnx][i][j][k][0] = fv[tpp][i][j][k][0];
          fv[tnx][i][j][k][1] = fv[tpp][im][j][k][1];
          fv[tnx][i][j][k][2] = fv[tpp][ip][j][k][2];
          fv[tnx][i][j][k][3] = fv[tpp][i][jm][k][3];
          fv[tnx][i][j][k][4] = fv[tpp][i][jp][k][4];
          fv[tnx][i][j][k][5] = fv[tpp][i][j][km][5];
          fv[tnx][i][j][k][6] = fv[tpp][i][j][kp][6];
          fv[tnx][i][j][k][7] = fv[tpp][im][jm][k][7];
          fv[tnx][i][j][k][8] = fv[tpp][ip][jp][k][8];
          fv[tnx][i][j][k][9] = fv[tpp][im][jp][k][9];
          fv[tnx][i][j][k][10] = fv[tpp][ip][jm][k][10];
          fv[tnx][i][j][k][11] = fv[tpp][im][j][km][11];
          fv[tnx][i][j][k][12] = fv[tpp][ip][j][kp][12];
          fv[tnx][i][j][k][13] = fv[tpp][im][j][kp][13];
          fv[tnx][i][j][k][14] = fv[tpp][ip][j][km][14];
          fv[tnx][i][j][k][15] = fv[tpp][i][jm][km][15];
          fv[tnx][i][j][k][16] = fv[tpp][i][jp][kp][16];
          fv[tnx][i][j][k][17] = fv[tpp][i][jm][kp][17];
          fv[tnx][i][j][k][18] = fv[tpp][i][jp][km][18];
          fv[tnx][i][j][k][19] = fv[tpp][im][jm][km][19];
          fv[tnx][i][j][k][20] = fv[tpp][ip][jm][km][20];
          fv[tnx][i][j][k][21] = fv[tpp][im][jp][km][21];
          fv[tnx][i][j][k][22] = fv[tpp][ip][jp][km][22];
          fv[tnx][i][j][k][23] = fv[tpp][im][jm][kp][23];
          fv[tnx][i][j][k][24] = fv[tpp][ip][jm][kp][24];
          fv[tnx][i][j][k][25] = fv[tpp][im][jp][kp][25];
          fv[tnx][i][j][k][26] = fv[tpp][ip][jp][kp][26];
        }
      }
    }

  }

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

    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        for (k=0;k<NNz;k++) {
          if (kind[i][j][k]==0) {

            setfeq(i,j,k);
            for (iv=0;iv<27;iv++) {
              fv[tnx][i][j][k][iv] -= omg*(fv[tnx][i][j][k][iv]-feq[iv]);
            }

          } else if (kind[i][j][k]==1) {

            fv[tnx][i][j][k][0] = 0.0;

            m = fv[tnx][i][j][k][1]; 
            fv[tnx][i][j][k][1]=fv[tnx][i][j][k][2];
            fv[tnx][i][j][k][2]=m;

            m = fv[tnx][i][j][k][3]; 
            fv[tnx][i][j][k][3]=fv[tnx][i][j][k][4];
            fv[tnx][i][j][k][4]=m;

            m = fv[tnx][i][j][k][5]; 
            fv[tnx][i][j][k][5]=fv[tnx][i][j][k][6];
            fv[tnx][i][j][k][6]=m;

            m = fv[tnx][i][j][k][7]; 
            fv[tnx][i][j][k][7]=fv[tnx][i][j][k][8];
            fv[tnx][i][j][k][8]=m;

            m = fv[tnx][i][j][k][9]; 
            fv[tnx][i][j][k][9]=fv[tnx][i][j][k][10];
            fv[tnx][i][j][k][10]=m;

            m = fv[tnx][i][j][k][11]; 
            fv[tnx][i][j][k][11]=fv[tnx][i][j][k][12];
            fv[tnx][i][j][k][12]=m;

            m = fv[tnx][i][j][k][13]; 
            fv[tnx][i][j][k][13]=fv[tnx][i][j][k][14];
            fv[tnx][i][j][k][14]=m;

            m = fv[tnx][i][j][k][15]; 
            fv[tnx][i][j][k][15]=fv[tnx][i][j][k][16];
            fv[tnx][i][j][k][16]=m;

            m = fv[tnx][i][j][k][17]; 
            fv[tnx][i][j][k][17]=fv[tnx][i][j][k][18];
            fv[tnx][i][j][k][18]=m;

            m = fv[tnx][i][j][k][19]; 
            fv[tnx][i][j][k][19]=fv[tnx][i][j][k][26];
            fv[tnx][i][j][k][26]=m;

            m = fv[tnx][i][j][k][20]; 
            fv[tnx][i][j][k][20]=fv[tnx][i][j][k][25];
            fv[tnx][i][j][k][25]=m;

            m = fv[tnx][i][j][k][21]; 
            fv[tnx][i][j][k][21]=fv[tnx][i][j][k][24];
            fv[tnx][i][j][k][24]=m;

            m = fv[tnx][i][j][k][22]; 
            fv[tnx][i][j][k][22]=fv[tnx][i][j][k][23];
            fv[tnx][i][j][k][23]=m;

          }
        }
      }
    }

  }


  void setfeq(int i, int j, int k) {
    double rh,vx,vy,vz,v2;
    double c1v,c2v,c3v,c4v,c5v,c6v,c7v,c8v,c9v,c10v,c11v,c12v,c13v,c14v,c15v,c16v,c17v,c18v;
    double c19v,c20v,c21v,c22v,c23v,c24v,c25v,c26v;
    double w0,w1,w2,w3;

    w0 = 8.0/27.0;
    w1 = 2.0/27.0;
    w2 = 1.0/54.0;
    w3 = 1.0/216.0;

    setRho(i,j,k);
    rh = rho[i][j][k];
    if (rh>0) {

      vx = ux[i][j][k];
      vy = uy[i][j][k];
      vz = uz[i][j][k];
      v2 = vx*vx+vy*vy+vz*vz;

      c1v = vx;
      c2v = -vx;
      c3v = vy;
      c4v = -vy;
      c5v = vz;
      c6v = -vz;

      c7v =  vx+vy;
      c8v = -vx-vy;
      c9v =  vx-vy;
      c10v = -vx+vy;
      c11v =  vx+vz;
      c12v = -vx-vz;
      c13v =  vx-vz;
      c14v = -vx+vz;
      c15v =  vy+vz;
      c16v = -vy-vz;
      c17v =  vy-vz;
      c18v = -vy+vz;

      c19v =  vx+vy+vz;
      c20v = -vx+vy+vz;
      c21v =  vx-vy+vz;
      c22v = -vx-vy+vz;
      c23v =  vx+vy-vz;
      c24v = -vx+vy-vz;
      c25v =  vx-vy-vz;
      c26v = -vx-vy-vz;

      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] = w1*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] = w1*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] = w1*rh*(1.0+3.0*c6v+4.5*c6v*c6v-1.5*v2);

      feq[7] = w2*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);
      feq[9] = w2*rh*(1.0+3.0*c9v+4.5*c9v*c9v-1.5*v2);
      feq[10] = w2*rh*(1.0+3.0*c10v+4.5*c10v*c10v-1.5*v2);
      feq[11] = w2*rh*(1.0+3.0*c11v+4.5*c11v*c11v-1.5*v2);
      feq[12] = w2*rh*(1.0+3.0*c12v+4.5*c12v*c12v-1.5*v2);
      feq[13] = w2*rh*(1.0+3.0*c13v+4.5*c13v*c13v-1.5*v2);
      feq[14] = w2*rh*(1.0+3.0*c14v+4.5*c14v*c14v-1.5*v2);
      feq[15] = w2*rh*(1.0+3.0*c15v+4.5*c15v*c15v-1.5*v2);
      feq[16] = w2*rh*(1.0+3.0*c16v+4.5*c16v*c16v-1.5*v2);
      feq[17] = w2*rh*(1.0+3.0*c17v+4.5*c17v*c17v-1.5*v2);
      feq[18] = w2*rh*(1.0+3.0*c18v+4.5*c18v*c18v-1.5*v2);

      feq[19] = w3*rh*(1.0+3.0*c7v+4.5*c7v*c7v-1.5*v2);
      feq[20] = w3*rh*(1.0+3.0*c8v+4.5*c8v*c8v-1.5*v2);
      feq[21] = w3*rh*(1.0+3.0*c9v+4.5*c9v*c9v-1.5*v2);
      feq[22] = w3*rh*(1.0+3.0*c10v+4.5*c10v*c10v-1.5*v2);
      feq[23] = w3*rh*(1.0+3.0*c11v+4.5*c11v*c11v-1.5*v2);
      feq[24] = w3*rh*(1.0+3.0*c12v+4.5*c12v*c12v-1.5*v2);
      feq[25] = w3*rh*(1.0+3.0*c13v+4.5*c13v*c13v-1.5*v2);
      feq[26] = w3*rh*(1.0+3.0*c14v+4.5*c14v*c14v-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;
      feq[9] = 0.0;
      feq[10] = 0.0;
      feq[11] = 0.0;
      feq[12] = 0.0;
      feq[13] = 0.0;
      feq[14] = 0.0;
      feq[15] = 0.0;
      feq[16] = 0.0;
      feq[17] = 0.0;
      feq[18] = 0.0;
      feq[19] = 0.0;
      feq[20] = 0.0;
      feq[21] = 0.0;
      feq[22] = 0.0;
      feq[23] = 0.0;
      feq[24] = 0.0;
      feq[25] = 0.0;
      feq[26] = 0.0;
    }

  }


  void setRho(int i, int j, int k) {
    int iv;
    double rh,rux,ruy,ruz;

    rh = 0.0; rux = 0.0; ruy = 0.0; ruz = 0.0;
    for (iv=0;iv<27;iv++) {
      rh += fv[tnx][i][j][k][iv];
      rux += ccx[iv]*fv[tnx][i][j][k][iv];
      ruy += ccy[iv]*fv[tnx][i][j][k][iv];
      ruz += ccz[iv]*fv[tnx][i][j][k][iv];
        
    }
    rho[i][j][k] = rh;
    if (rh>0.0) {
      ux[i][j][k] = rux/rh;
      uy[i][j][k] = ruy/rh;
      uz[i][j][k] = ruz/rh;
    } else {
      ux[i][j][k] = 0.0;
      uy[i][j][k] = 0.0;
      uz[i][j][k] = 0.0;
    }

  }

/*
  void setFan() {
    int i,j,k;
    double w0,w1,w2,w3;

    w0 = 8.0/27.0;
    w1 = 2.0/27.0;
    w2 = 1.0/54.0;
    w3 = 1.0/216.0;

    i = 0;
    for (j=0;j<NNy;j++) {
      for (k=0;k<NNz;k++) {
        fv[tnx][i][j][k][0] = w0;

        fv[tnx][i][j][k][1] = w1*(1.0+fan);
        fv[tnx][i][j][k][2] = w1*(1.0-fan);
        fv[tnx][i][j][k][3] = w1;
        fv[tnx][i][j][k][4] = w1;
        fv[tnx][i][j][k][5] = w1;
        fv[tnx][i][j][k][6] = w1;

        fv[tnx][i][j][k][7] = w2;
        fv[tnx][i][j][k][8] = w2;
        fv[tnx][i][j][k][9] = w2;
        fv[tnx][i][j][k][10] = w2;
        fv[tnx][i][j][k][11] = w2;
        fv[tnx][i][j][k][12] = w2;
        fv[tnx][i][j][k][13] = w2;
        fv[tnx][i][j][k][14] = w2;
        fv[tnx][i][j][k][15] = w2;
        fv[tnx][i][j][k][16] = w2;
        fv[tnx][i][j][k][17] = w2;
        fv[tnx][i][j][k][18] = w2;

        fv[tnx][i][j][k][19] = w3;
        fv[tnx][i][j][k][20] = w3;
        fv[tnx][i][j][k][21] = w3;
        fv[tnx][i][j][k][22] = w3;
        fv[tnx][i][j][k][23] = w3;
        fv[tnx][i][j][k][24] = w3;
        fv[tnx][i][j][k][25] = w3;
        fv[tnx][i][j][k][26] = w3;
      }
    }
  }
*/

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

}