/*
 * template - real-coded lattice gas model 3D
 *
 * Created by Ikeuchi Mitsuru on Apr 04 2006.
 * Copyright (c) 2006 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2006.04.09  created
 *
 */

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

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

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

  Button bt_reset, bt_startStop;
  Choice ch_view;
  Scrollbar sc_scale;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 50;

  double t = 0.0;
  double dt = 1.0;

  double dx = 1.0;
  double dy = 1.0;
  double dz = 1.0;

  int Nsx = 80;
  int Nsy = 48;
  int Nsz = 48;
  double xMax = Nsx*dx;
  double yMax = Nsy*dy;
  double zMax = Nsz*dz;
  int NinCell = 4;

  int dispMode = 3;
  int resetSW = 0;
  int started = 1;
  int incRate = 1;
  double gravity = 0.0;

  double dispScale = (240.0/yMax);
  double viewScale = 1.0;

  int dragSW = 2;
  int dgX, dgY, dgXb, dgYb;
  double xShift = 0.0;
  double yShift = 0.0;

  double cx = 0.5*xMax;
  double cy = 0.5*yMax;
  double cz = 0.5*zMax;
  double viewTheta = -15.0*3.14/180.0;
  double viewFai = 18.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);

  int NNp = Nsx*Nsy*Nsz*NinCell*2/3;
  double xx[] = new double [NNp];
  double yy[] = new double [NNp];
  double zz[] = new double [NNp];
  double vx[] = new double [NNp];
  double vy[] = new double [NNp];
  double vz[] = new double [NNp];

  int NNs = 20;
  int section[][][][] = new int [Nsx][Nsy][Nsz][NNs];
  int cellAttribute[][][] = new int [Nsx][Nsy][Nsz];
  double refTemp[] = { 0.0, 0.0, 0.0, 0.0 };

  int vDistribution[][] = new int [200][3];

  int NplotMax = 32000;
  double plotx[] = new double [NplotMax];
  double ploty[] = new double [NplotMax];
  double plotz[] = new double [NplotMax];
  double plotColor[] = new double [NplotMax];
  double plotShade[] = new double [NplotMax];
  double plotSize[] = new double [NplotMax];
  double rotPlotx[] = new double [NplotMax];
  double rotPloty[] = new double [NplotMax];
  double rotPlotz[] = new double [NplotMax];

  int vSpace[][][] = new int [20][20][20];

  int srtz[][] = new int[2][NplotMax];

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

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

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

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

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

    ch_view = new Choice();
    ch_view.addItem("sample 8k");
    ch_view.addItem("velocity 8k");
    ch_view.addItem("density");
    ch_view.addItem("flow");
    ch_view.addItem("temperature");
    ch_view.addItem("pressure");
    ch_view.addItem("v space");
    ch_view.addItem("v distr.");
    ch_view.addItemListener(this);
    ch_view.select("flow");

    sc_scale= new Scrollbar(Scrollbar.HORIZONTAL,100,10,50,510);
    sc_scale.addAdjustmentListener(this);

    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(new Label("  "));
    pnl.add(new Label("  "));
    pnl.add(sc_scale);
    pnl.add(ch_view);
    add(pnl,"North");

    setWaku();
    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 itemStateChanged(ItemEvent ev){
    if (ev.getSource() == ch_view) {
      dispMode = ch_view.getSelectedIndex();
    }
  }

  public void actionPerformed(ActionEvent ev){
    if(ev.getSource() == bt_reset){ 
      resetSW = 1;
    } else if (ev.getSource() == bt_startStop){
      if (started==0) {
        started = 1;
      } else {
        started = 0;
      }
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_scale) { 
      viewScale = 0.01*sc_scale.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;
    }
	
    if (dragSW==1) {
      xShift += (1.0/viewScale/dispScale)*(dgX-dgXb);
      yShift += (1.0/viewScale/dispScale)*(dgY-dgYb);
    } else if (dragSW==2) {
      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 gbx, gby;
    double sc,temp;

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

    sc = viewScale*dispScale;
    gbx = 60+(int)(sc*xShift);
    gby = 120+(int)(sc*yShift);

    if (dispMode==0) {
      oval3DPlot(gbx,gby,sc,8000);

    } else if (dispMode==1) {
      velocity3DPlot(gbx,gby,sc, 8000, 2.0);

    } else if (dispMode==2) {
      field3DPlot(0,gbx,gby,sc,1.0);

    } else if (dispMode==3) {
      flow3DPlot(gbx,gby,sc,10.0);

    } else if (dispMode==4) {
      field3DPlot(1,gbx,gby,sc,1.5);

    } else if (dispMode==5) {
      field3DPlot(2,gbx,gby,sc,0.3);

    } else if (dispMode==6) {
      vSpace3DPlot(gbx,gby,sc,0.01);

    } else if (dispMode==7) {
      velocityDistributionPlot(80,100,0.01);

    } else if (dispMode==8) {
      ;

    }

    //velocityDistributionPlot(410,200,0.01);


    gOff.setColor(Color.black);
    gOff.drawString("t="+(int)(t*10.0+0.5)/10.0+" ",630/6*0+10,40);
    if (started==1) {
      gOff.drawString("started",630/6*1+10,40);
    } else {
      gOff.drawString("stopped",630/6*1+10,40);
    }
    gOff.drawString("scale="+(int)(viewScale*100.0+0.5)/100.0+" ",630/6*4+10,40);
    gOff.drawString("view",630/6*5+10,40);

    gOff.setColor(Color.blue);
    temp = systemTemperature();
    gOff.drawString("Box="+Nsx+"x"+Nsy+"",630/6*0+10,60);
    gOff.drawString("N="+NNp+"",630/6*1+10,60);
    gOff.drawString("T="+(int)(temp*1000.0+0.5)/1000.0+"",630/6*2+10,60);
    gOff.drawString("nue="+(int)(viscosity(temp)*1000.0+0.5)/1000.0+" ",630/6*3+10,60);
    gOff.drawString("maxNc="+maxSection()+" ",630/6*5+10,60);

//    gOff.drawString("theta="+(float)(viewTheta*180/3.1416)+" ",500,100);
//    gOff.drawString("fai="+(float)(viewFai*180/3.1416)+" ",500,120);
	
  }


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

  void oval3DPlot(int gx, int gy, double sc,int Nsample) {
    int i,j,ix,iy,ir,nn;

    nn = setPlotData(Nsample);
    drawBackBox(gx,gy,sc);

    for (i=0; i<nn; i++) {
      j = srtz[0][i];
      ix = gx+(int)(rotPlotx[j]*sc+0.5);
      iy = gy+(int)(rotPloty[j]*sc+0.5);
      ir = (int)(2.0*plotSize[j]*viewScale*(0.8+0.2*rotPlotz[j]/zMax));
      gOff.setColor(Color.getHSBColor((float)(plotColor[j]),0.9f, (float)(plotShade[j]+0.5*rotPlotz[j]/zMax)));
      gOff.drawOval(ix-ir/2,iy-ir/2, ir, ir);
    }
    drawFrontBox(gx,gy,sc);
  }

  int setPlotData(int Nsample) {
    int i,ix,iy,iz,nn,aki;

    for (i=0; i<Nsample; i++) {
      plotx[i] = xx[i]; ploty[i] = yy[i]; plotz[i] = zz[i];
      plotColor[i] = 0.33;
      plotShade[i] = 0.3;
      plotSize[i] = 4.0;
    }

    for (ix=0; ix<Nsx; ix++) {
      for (iy=0; iy<Nsy; iy++) {
        for (iz=0; iz<Nsz; iz++) {
          if (cellAttribute[ix][iy][iz]>0) {
            plotx[i] = (ix+0.5)*dx;
            ploty[i] = (iy+0.5)*dy;
            plotz[i] = (iz+0.5)*dz;

            plotSize[i] = 2.0;
            plotColor[i] = 0.85;
            plotShade[i] = 0.3;
            i += 1;
          }
        }
      }
    }
    nn = i;

    rotate(viewTheta, viewFai, nn);
    zSort(nn);

    return (nn);
  }


  void velocity3DPlot(int gx, int gy, double sc, int Nsample, double mag) {
    int i,j,ix,iy,ir,nn;
    double p2x,p2y,p2z,rotp2x,rotp2y,rotp2z;

    nn = setPlotData(Nsample);
    drawBackBox(gx,gy,sc);

    for (i=0; i<nn; i++) {
      j = srtz[0][i];
      ix = gx+(int)(rotPlotx[j]*sc+0.5);
      iy = gy+(int)(rotPloty[j]*sc+0.5);
      if (j<Nsample) {
        p2x = xx[j]+mag*vx[j]*dt;
        p2y = yy[j]+mag*vy[j]*dt;
        p2z = zz[j]+mag*vz[j]*dt;
        rotp2y = cosFi*(p2y-cy)+sinFi*(p2z-cz) + cy; 
        rotp2z = -sinFi*(p2y-cy)+cosFi*(p2z-cz) + cz;
        rotp2x = cosTh*(p2x-cx)+sinTh*(rotp2z-cz) + cx; 
        rotp2z = -sinTh*(p2x-cx)+cosTh*(rotp2z-cz) + cz;
        gOff.setColor(Color.red); if (vx[j]>0) gOff.setColor(Color.blue);
        gOff.drawLine(ix,iy, gx+(int)(rotp2x*sc+0.5),gy+(int)(rotp2y*sc+0.5));
      } else {
        ir = (int)(2.0*plotSize[j]*viewScale*(0.8+0.2*rotPlotz[j]/zMax));
        gOff.setColor(Color.getHSBColor((float)(plotColor[j]),0.9f, (float)(plotShade[j]+0.5*rotPlotz[j]/zMax)));
        gOff.drawOval(ix-ir/2,iy-ir/2, ir, ir);
      }
    }
    drawFrontBox(gx,gy,sc);
  }


  void field3DPlot(int imode, int gx, int gy, double sc, double mag) {
    int i,j,ix,iy,iz,ir,nn;
    double ff,gg;

    i = 0;
    for (ix=0; ix<Nsx; ix+=2) {
      for (iy=0; iy<Nsy; iy+=2) {
        for (iz=0; iz<Nsz; iz+=2) {
          plotx[i] = (ix+1)*dx;
          ploty[i] = (iy+1)*dy;
          plotz[i] = (iz+1)*dz;
          if (cellAttribute[ix][iy][iz]==0) {
            ff = (0.6/NinCell)*meanDensity(ix,iy,iz,2); if (ff>0.999) ff = 0.999;
            plotSize[i] = 8.0*ff;
            if (imode==0) {
              plotColor[i] = 0.33;
              plotShade[i] = 0.3;
            } else if (imode==1) {
              gg = mag*meanField(1,ix,iy,iz,2); if (gg>0.999) gg = 0.999;
              plotColor[i] = 0.66-0.66*gg;
              plotShade[i] = 0.3;
            } else if (imode==2) {
              gg = mag*meanField(2,ix,iy,iz,2); if (gg>0.999) gg = 0.999;
              plotColor[i] = 0.66-0.66*gg;
              plotShade[i] = 0.3;
            }
          } else if (cellAttribute[ix][iy][iz]>0) {
            plotSize[i] = 4.0;
            plotColor[i] = 0.85;
            plotShade[i] = 0.3;
          }
          i += 1;
        }
      }
    }
    nn = i;
    rotate(viewTheta, viewFai, nn);
    zSort(nn);

    drawBackBox(gx,gy,sc);

    for (i=0; i<nn; i++) {
      j = srtz[0][i];
      ix = gx+(int)(rotPlotx[j]*sc+0.5);
      iy = gy+(int)(rotPloty[j]*sc+0.5);
      ir = (int)(2.0*plotSize[j]*viewScale*(0.8+0.2*rotPlotz[j]/zMax));
	  if (ir>0) {
        gOff.setColor(Color.getHSBColor((float)(plotColor[j]), 0.9f,(float)(plotShade[j]+0.5*rotPlotz[j]/zMax)));
        gOff.drawOval(ix-ir/2,iy-ir/2, ir, ir);
      }
    }
    drawFrontBox(gx,gy,sc);
  }

  double meanDensity(int ix, int iy, int iz,int len) {
    int i,j,k,n,nn;

    n = 0;
    for (i=ix; i<ix+len; i++) {
      for (j=iy; j<iy+len; j++) {
        for (k=iz; k<iz+len; k++) {
          nn = section[i][j][k][0];
          n += nn;
        }
      }
    }
    return( n/(1.0*len*len*len) );
  }

  double meanField(int imode, int ix, int iy, int iz,int len) {
    int i,j,k,n,nn,ip,iq;
    double mvx,mvy,mvz,ke,rho;

    n = 0;
    mvx = 0.0; mvy = 0.0; mvz = 0.0;
    for (i=ix; i<ix+len; i++) {
      for (j=iy; j<iy+len; j++) {
        for (k=iz; k<iz+len; k++) {
          nn = section[i][j][k][0];
          n += nn;
          for (ip=1; ip<nn; ip++) {
            iq = section[i][j][k][ip];
            mvx += vx[iq]; mvy += vy[iq]; mvz += vz[iq];
          }
        }
      }
    }
	mvx = mvx/(n+1.0e-20); mvy = mvy/(n+1.0e-20); mvz = mvz/(n+1.0e-20);
    rho = n/(1.0*len*len*len);

	ke = 0.0;
    for (i=ix; i<ix+len; i++) {
      for (j=iy; j<iy+len; j++) {
        for (k=iz; k<iz+len; k++) {
        nn = section[i][j][k][0];
          for (ip=1; ip<nn; ip++) {
            iq = section[i][j][k][ip];
            ke += (vx[iq]-mvx)*(vx[iq]-mvx)+(vy[iq]-mvy)*(vy[iq]-mvy)+(vz[iq]-mvz)*(vz[iq]-mvz);
          }
        }
      }
    }
    ke = 0.5*ke/(n+1.0e-20);
    if (imode==1) return( ke );
    if (imode==2) return( ke*rho );
    return(-1.0);
  }

/*------------------------------ flow 3D plot ----------*/

  void flow3DPlot(int gx, int gy, double sc, double mag) {
    int i,j,ix,iy,iz,ir,n,nn,nnn,ii,jj,kk,ip,iq,igx,igy;
    double rho,p2x,p2y,p2z,rotp2x,rotp2y,rotp2z;
    double mvx,mvy,mvz;

    i = 0;
    for (ix=0; ix<Nsx; ix+=2) {
      for (iy=0; iy<Nsy; iy+=2) {
        for (iz=0; iz<Nsz; iz+=2) {
          plotx[i] = (ix+1)*dx;
          ploty[i] = (iy+1)*dy;
          plotz[i] = (iz+1)*dz;
          if (cellAttribute[ix][iy][iz]==0) {
            plotSize[i] = 0.0;
            plotColor[i] = 0.33;
            plotShade[i] = 0.3;
          } else if (cellAttribute[ix][iy][iz]>0) {
            plotSize[i] = 4.0;
            plotColor[i] = 0.85;
            plotShade[i] = 0.3;
          }
        i += 1;
        }
      }
    }
    nn = i;
    rotate(viewTheta, viewFai, nn);
    zSort(nn);

    drawBackBox(gx,gy,sc);

    for (i=0; i<nn; i++) {
      j = srtz[0][i];
      igx = gx+(int)(rotPlotx[j]*sc+0.5);
      igy = gy+(int)(rotPloty[j]*sc+0.5);
      if (plotSize[j]<1.0) {
        ix = (int)(Nsx*plotx[j]/xMax); if (ix>=Nsx) ix = Nsx-1;
        iy = (int)(Nsy*ploty[j]/yMax); if (iy>=Nsy) iy = Nsy-1;
        iz = (int)(Nsz*plotz[j]/zMax); if (iz>=Nsz) iz = Nsz-1;
        ix = (ix/2)*2; iy = (iy/2)*2; iz = (iz/2)*2;
		
        n = 0;
        mvx = 0.0; mvy = 0.0; mvz = 0.0;
        for (ii=ix; ii<ix+2; ii++) {
          for (jj=iy; jj<iy+2; jj++) {
            for (kk=iz; kk<iz+2; kk++) {
              nnn = section[ii][jj][kk][0];
              n += nnn;
              for (ip=1; ip<nnn; ip++) {
                iq = section[ii][jj][kk][ip];
                mvx += vx[iq]; mvy += vy[iq]; mvz += vz[iq];
              }
            }
          }
        }
        mvx = mvx/(n+1.0e-20); mvy = mvy/(n+1.0e-20); mvz = mvz/(n+1.0e-20);
        rho = 2.0*n/(NinCell*2*2*2); if (rho>=1.0) rho = 1.0;

        p2x = plotx[j]+mag*mvx*dt;
        p2y = ploty[j]+mag*mvy*dt;
        p2z = plotz[j]+mag*mvz*dt;
        rotp2y = cosFi*(p2y-cy)+sinFi*(p2z-cz) + cy; 
        rotp2z = -sinFi*(p2y-cy)+cosFi*(p2z-cz) + cz;
        rotp2x = cosTh*(p2x-cx)+sinTh*(rotp2z-cz) + cx; 
        rotp2z = -sinTh*(p2x-cx)+cosTh*(rotp2z-cz) + cz;
        gOff.setColor(Color.red);
        gOff.setColor(Color.getHSBColor(0.01f,(float)(0.9*rho), (float)(plotShade[j]+0.5*rotPlotz[j]/zMax)));
        if (mvx>0) {
          gOff.setColor(Color.getHSBColor(0.66f,(float)(0.9*rho), (float)(plotShade[j]+0.5*rotPlotz[j]/zMax)));
        }
		if (n>0) {
         gOff.drawLine(igx,igy, gx+(int)(rotp2x*sc+0.5),gy+(int)(rotp2y*sc+0.5));
        }

      } else {
        ir = (int)(2.0*plotSize[j]*viewScale*(0.8+0.2*rotPlotz[j]/zMax));
        gOff.setColor(Color.getHSBColor((float)(plotColor[j]),0.9f, (float)(plotShade[j]+0.5*rotPlotz[j]/zMax)));
        gOff.drawOval(igx-ir/2,igy-ir/2, ir, ir);
      }
    }

    drawFrontBox(gx,gy,sc);
  }

/*--------------------------- velocity space plot ---------*/

  void vSpace3DPlot(int gx, int gy, double sc, double mag) {
    int i,ix,iy,iz,j,jx,jy,jz,ir,nn;
    double ff;

    for (ix=0; ix<20; ix++) {
      for (iy=0; iy<20; iy++) {
        for (iz=0; iz<20; iz++) {
          vSpace[ix][iy][iz] = 0;
        }
      }
    }

    for (i=0; i<NNp; i++) {
      ix = (int)((vx[i]+2.0)*5.0);
      if (ix<0) ix = 0;
      if (ix>19) ix = 19;
      iy = (int)((vy[i]+2.0)*5.0);
      if (iy<0) iy = 0;
      if (iy>19) iy = 19;
      iz = (int)((vz[i]+2.0)*5.0);
      if (iz<0) iz = 0;
      if (iz>19) iz = 19;
	  vSpace[ix][iy][iz] += 1;
    }

    i = 0;
    for (ix=0; ix<30; ix++) {
      for (iy=0; iy<20; iy++) {
        for (iz=0; iz<20; iz++) {
          plotx[i] = (2*(ix-10)+Nsx/2)*dx;
          ploty[i] = (2*(iy-10)+Nsy/2)*dy;
          plotz[i] = (2*(iz-10)+Nsz/2)*dz;
          i += 1;
        }
      }
    }
    nn = 20*20*20;

    rotate(viewTheta, viewFai, nn);
    zSort(nn);

    drawAxis(gx,gy,sc, 0);
    drawAxis(gx,gy,sc, 1);
    drawAxis(gx,gy,sc, 2);

    for (i=0; i<nn; i++) {
      j = srtz[0][i];
      jx = j/400; jy = (j-400*jx)/20; jz = (j-400*jx)%20;
      ff = mag*vSpace[jx][jy][jz]; if (ff>0.999) ff = 0.999;
      ir = (int)(3.0*ff+0.5);
      if (ir>0) {
        gOff.setColor(Color.getHSBColor((float)(0.66-0.63*ff),0.8f,(float)(0.3+0.4*rotPlotz[j]/zMax)));
        gOff.fillOval(gx+(int)(rotPlotx[j]*sc+0.5)-ir/2,gy+(int)(rotPloty[j]*sc+0.5)-ir/2, ir, ir);
      }
    }
  }

  void drawAxis(int gx, int gy, double sc, int imode) {
    double px,py,pz,px2,py2,pz2;
    double rotpx,rotpy,rotpz,rotpx2,rotpy2,rotpz2;

    px = Nsx/2*dx; py = Nsy/2*dy; pz = Nsz/2*dz;
    if (imode==0) {
      px2 = Nsx/2*dx+40; py2 = Nsy/2*dy; pz2 = Nsz/2*dz;
    } else if (imode==1) {
      px2 = Nsx/2*dx; py2 = Nsy/2*dy-40; pz2 = Nsz/2*dz;
    } else {
      px2 = Nsx/2*dx; py2 = Nsy/2*dy; pz2 = Nsz/2*dz+40;
    }

    rotpy = cosFi*(py-cy)+sinFi*(pz-cz) + cy; 
    rotpz = -sinFi*(py-cy)+cosFi*(pz-cz) + cz;
    rotpx = cosTh*(px-cx)+sinTh*(rotpz-cz) + cx; 
    rotpz = -sinTh*(px-cx)+cosTh*(rotpz-cz) + cz;

    rotpy2 = cosFi*(py2-cy)+sinFi*(pz2-cz) + cy; 
    rotpz2 = -sinFi*(py2-cy)+cosFi*(pz2-cz) + cz;
    rotpx2 = cosTh*(px2-cx)+sinTh*(rotpz2-cz) + cx; 
    rotpz2 = -sinTh*(px2-cx)+cosTh*(rotpz2-cz) + cz;

    gOff.setColor(Color.red);
	gOff.drawLine(gx+(int)(rotpx*sc+0.5),gy+(int)(rotpy*sc+0.5),gx+(int)(rotpx2*sc+0.5),gy+(int)(rotpy2*sc+0.5));
  }

/* -------------------------------  box waku  -------------- */

  void drawBackBox(int gbx, int gby, double sc) {
    int i,gx,gy,g2x,g2y,farPoint;

    for (i=0; i<8; i++) {
      pwky[i] = cosFi*(wky[i]-cy)+sinFi*(wkz[i]-cz) + cy; 
      pwkz[i] = -sinFi*(wky[i]-cy)+cosFi*(wkz[i]-cz) + cz;
    } 
    for (i=0; i<8; i++) {
      pwkx[i] = cosTh*(wkx[i]-cx)+sinTh*(pwkz[i]-cz) + cx; 
      pwkz[i] = -sinTh*(wkx[i]-cx)+cosTh*(pwkz[i]-cz) + cz;
    }  

    farPoint = findMin();
    for (i=0; i<12; i++) {
      gx = (int)(sc*pwkx[boxp[i][0]])+gbx;
      gy = (int)(sc*pwky[boxp[i][0]])+gby;
      g2x = (int)(sc*pwkx[boxp[i][1]])+gbx;
      g2y = (int)(sc*pwky[boxp[i][1]])+gby;
      if (boxp[i][0]==farPoint || boxp[i][1]==farPoint) { 
        gOff.setColor(Color.lightGray);
        gOff.drawLine(gx,gy, g2x, g2y);
      } 
    }
  }

  void drawFrontBox(int gbx, int gby, double sc) {
    int i,gx,gy,g2x,g2y,farPoint;

    farPoint = findMin();
    for (i=0; i<12; i++) {
      gx = (int)(sc*pwkx[boxp[i][0]])+gbx;
      gy = (int)(sc*pwky[boxp[i][0]])+gby;
      g2x = (int)(sc*pwkx[boxp[i][1]])+gbx;
      g2y = (int)(sc*pwky[boxp[i][1]])+gby;
      if (boxp[i][0]!=farPoint && boxp[i][1]!=farPoint) { 
        gOff.setColor(Color.black);
        gOff.drawLine(gx,gy, g2x, g2y);
      }
    }
  }

  void setWaku() {
    int i,ix,iy,iz;
    double a;

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

  int findMin() {
    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);
  }

/* ------------------------------- display utilities -------------- */

  void zSort(int Nmax) {
    int i,j,w;

    for (i=0; i<Nmax; i++) {
      srtz[0][i] = i;
      srtz[1][i] = (int)(32767.0*rotPlotz[i]/zMax);
    }
    qSort(0,Nmax-1);    
  }

  void qSort(int i,int j) {
    int k,p;

    if(i==j) return;
	p = i+1;
    while (p<=j && srtz[1][i]==srtz[1][p]) p++;
    if (p<=j) {
      if (srtz[1][i]>=srtz[1][p]) p = i;
      k = qSortPartition(i,j,srtz[1][p]);
      qSort(i,k-1);
      qSort(k,j);
    }
  }

  int qSortPartition(int i,int j, int x) {
    int l,r,w;

    l = i; r = j;
    while (l<=r) {
      while (l<=j && srtz[1][l]<x) l++;
      while (r>=i && srtz[1][r]>=x) r--;
      if (l>r) break;
      w = srtz[0][l]; srtz[0][l] = srtz[0][r]; srtz[0][r] = w;
      w = srtz[1][l]; srtz[1][l] = srtz[1][r]; srtz[1][r] = w;
      l++; r--;
    }
    return ( l );
  }

  void rotate(double th,double fi, int Nmax) {
    int i;
    double cosTh,sinTh,cosFi,sinFi;

    cosTh = Math.cos(th); sinTh = Math.sin(th);
    cosFi = Math.cos(fi); sinFi = Math.sin(fi);
    for (i=0; i<Nmax; i++) {
      rotPloty[i] = cosFi*(ploty[i]-cy)+sinFi*(plotz[i]-cz) + cy; 
      rotPlotz[i] = -sinFi*(ploty[i]-cy)+cosFi*(plotz[i]-cz) + cz;

      rotPlotx[i] = cosTh*(plotx[i]-cx)+sinTh*(rotPlotz[i]-cz) + cx; 
      rotPlotz[i] = -sinTh*(plotx[i]-cx)+cosTh*(rotPlotz[i]-cz) + cz;
    }  
  }

/* ------------------- velocity distribution plot --------- */

  void velocityDistributionPlot(int gx, int gy, double yMag) {
    int i,n1,n2;
    int hoganWidth,hoganHight,xDiv,yDiv;
    double c1,c2,cm;

    setVelocityDistribution();

    hoganWidth = 200; hoganHight = 160; xDiv = 50; yDiv = 80;

    for (i=0; i<200; i++) {
      n1 = (int)(yMag*vDistribution[i][0]);
      if (n1>0) {
        gOff.setColor(Color.lightGray);
        gOff.drawLine(gx+i,gy+hoganHight,gx+i,gy+hoganHight-n1);
      }
    }

    drawHogan(gx,gy,hoganWidth,hoganHight,xDiv,yDiv);
    gOff.setColor(Color.black);
    gOff.drawString("0",gx-5,gy+hoganHight+15);
    gOff.drawString("2",gx+hoganWidth/2-5,gy+hoganHight+15);
    gOff.drawString("4",gx+hoganWidth-5,gy+hoganHight+15);
    gOff.drawString("velocity distribution",gx+20,gy+hoganHight+35);
  }

  void setVelocityDistribution() {
    int i,iv;
    double v;

    for(i=0;i<200;i++) {
      vDistribution[i][0] = 0;
    }

    for(i=0;i<NNp;i++) {
      v = Math.sqrt(vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i]);
      iv = (int)(50.0*v);
      if (iv<200) {
        vDistribution[iv][0] += 1;
      }
    }
  }

  void drawHogan(int gx, int gy, int hoganWidth, int hoganHight, int xDiv, int yDiv) {
    int i;

    gOff.setColor(Color.gray);
    for (i=0; i<hoganWidth; i+=xDiv) {
      gOff.drawLine(gx+i,gy, gx+i, gy+hoganHight);
    }
    for (i=0; i<hoganHight; i+=yDiv) {
      gOff.drawLine(gx,gy+hoganHight-i, gx+hoganWidth, gy+hoganHight-i);
    }
    gOff.setColor(Color.black);
    gOff.drawRect(gx,gy,hoganWidth,hoganHight);
  }

/*-------------------------- statistics -------- */

  double totalKineticEnergy() {
    int i;
    double tke;

    tke = 0.0;
    for(i=0;i<NNp;i++) {
      tke += 0.5*(vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i]);
    }
    return( tke );
  }

  double systemTemperature() {
    int i;
    double tke;

    tke = totalKineticEnergy();
    return( tke/NNp );
  }

  double viscosity(double temp) {
    int i;
    double rho;

    rho = NNp/(1.0*Nsx*Nsy*Nsz);
    return( 1.0/12.0+0.5*temp*(rho+1-Math.exp(-rho)/(rho-1+Math.exp(-rho))) );
  }

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

  void setInitialCondition() {
    int i,ic,jc,kc;
    double x,y,z;

    t = 0.0;

    for (ic=0; ic<Nsx; ic++) {
      for (jc=0; jc<Nsy; jc++) {
        for (kc=0; kc<Nsz; kc++) {
          cellAttribute[ic][jc][kc] = 0;
          if (ic>=Nsx-Nsx/3 && ic<=Nsx-Nsx/3+3 && (jc-Nsy/2)*(jc-Nsy/2)+(kc-Nsz/2)*(kc-Nsz/2)>30) {
            cellAttribute[ic][jc][kc] = 1;
          }
        }
      }
    }

    i = 0;
    while ( i<NNp ) {
      x = xMax*Math.random()*0.67;
      y = yMax*Math.random();
	  z = zMax*Math.random();
      if (cellAttributeAt(x,y,z)==0) {
        xx[i] = x;
        yy[i] = y;
        zz[i] = z;
        vx[i] = 1.6*(Math.random()-0.5);
        vy[i] = 1.6*(Math.random()-0.5);
        vz[i] = 1.6*(Math.random()-0.5);
        i += 1;
      }
    }
  }

/* ----------------------------- timeEvolution -------------- */

  void timeEvolution() {
    int i;

    if (resetSW==1) {
      t = 0.0;
      setInitialCondition();
      started = 0;
      resetSW = 0;
    }
	
    if (started==1) {
      for (i=0; i<incRate; i++) {
        timeStep();
      }
    }
  }

  void timeStep() {
    int i;

    t = t + dt;

    movement();
    collision();
  }

  void movement() {
    int i,ca;
    double rr,r,temp;

    for (i=0; i<NNp; i++) {
      vy[i] += gravity*dt;
      xx[i] += vx[i]*dt;
      yy[i] += vy[i]*dt;
      zz[i] += vz[i]*dt;

      ca = cellAttributeAt(xx[i],yy[i],zz[i]);
      if (ca>0) {
        xx[i] -= vx[i]*dt;
        yy[i] -= vy[i]*dt;
		zz[i] -= vz[i]*dt;
        r = 1.0;
        if (refTemp[ca]>0.0) {
          temp = 0.5*(vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i]);
          r = Math.sqrt(refTemp[ca]/temp)*(1.0-1.0*(Math.random()-0.5));
        } else {
          r = 1.0;
        }
        vx[i] = -r*vx[i];
        vy[i] = -r*vy[i];
        vz[i] = -r*vz[i];
      }
    }

    rr = 1.0;
    for (i=0; i<NNp; i++) {
      if (xx[i] < 0.0) { 
        xx[i] = 0.0; vx[i] = -rr*vx[i]; vy[i] = rr*vy[i]; vz[i] = rr*vz[i];
      }
      if (xx[i] > xMax) {
        xx[i] = xMax; vx[i] = -rr*vx[i]; vy[i] = rr*vy[i]; vz[i] = rr*vz[i];
      }
      if (yy[i] < 0.0) { 
        yy[i] = 0.0; vx[i] = rr*vx[i]; vy[i] = -rr*vy[i]; vz[i] = rr*vz[i];
      }
      if (yy[i] > yMax) {
        yy[i] = yMax; vx[i] = rr*vx[i]; vy[i] = -rr*vy[i]; vz[i] = rr*vz[i];
      }
      if (zz[i] < 0.0) { 
        zz[i] = 0.0; vx[i] = rr*vx[i]; vy[i] = rr*vy[i]; vz[i] = -rr*vz[i];
      }
      if (zz[i] > zMax) {
        zz[i] = zMax; vx[i] = rr*vx[i]; vy[i] = rr*vy[i]; vz[i] = -rr*vz[i];
      }
    }
  }

  int cellAttributeAt(double x, double y, double z) {
    int ix,iy,iz;

    ix = (int)(x/dx); if (ix>=Nsx) ix = Nsx-1;
    if (ix<0) ix = 0;
    iy = (int)(y/dy); if (iy>=Nsy) iy = Nsy-1;
    if (iy<0) iy = 0;
    iz = (int)(z/dz); if (iz>=Nsz) iz = Nsz-1;
    if (iz<0) iz = 0;
    return( cellAttribute[ix][iy][iz] );
  }

  void collision() {
    int ic,jc,kc,nn,ip,k;

    dividingSection();

    for(ic=0;ic<Nsx;ic++) {
      for(jc=0;jc<Nsy;jc++) {
        for(kc=0;kc<Nsz;kc++) {
          nn = section[ic][jc][kc][0];
          if (nn>1) {
            if (cellAttribute[ic][jc][kc]==0) {
              collisionInTheCell(ic,jc,kc);
            }
          }
        }
      }
    }
  }

  void collisionInTheCell(int ic, int jc, int kc) {
    int i,ip,n;
    double thx,thy,thz,cx,sx,cy,sy,cz,sz;
    double vxm,vym,vzm,vxs,vys,vzs;

    thx = 3.14159*(2.0*Math.random()-1.0);
    thy = 3.14159*(2.0*Math.random()-1.0);
    thz = 3.14159*(2.0*Math.random()-1.0);
    cx = Math.cos(thx); sx = Math.sin(thx);
    cy = Math.cos(thy); sy = Math.sin(thy);
    cz = Math.cos(thz); sz = Math.sin(thz);

    n = section[ic][jc][kc][0];
    vxm = 0.0; vym = 0.0; vzm = 0.0;
    for(i=1;i<=n;i++) {
      ip = section[ic][jc][kc][i];
      vxm += vx[ip]; vym += vy[ip]; vzm += vz[ip];
    }
    vxm = vxm/n; vym = vym/n; vzm = vzm/n;

    for(i=1;i<=n;i++) {
      ip = section[ic][jc][kc][i];
      vxs = vx[ip]-vxm;
      vys = vy[ip]-vym;
      vzs = vz[ip]-vzm;
      vx[ip] = vxm +cy*cz*vxs +(sx*sy*cz-cx*sz)*vys +(cx*sy*cz+sx*sz)*vzs;
      vy[ip] = vym +cy*sz*vxs +(sx*sy*sz+cx*cz)*vys +(cx*sy*sz-sx*cz)*vzs;
      vz[ip] = vzm -sy*vxs +sx*cy*vys +cx*cy*vzs;
    }
  }


/*-----  dividing section  -----*/

  void dividingSection() {
    int i,j,k,ip,iq;

    for(i=0;i<Nsx;i++) {
      for(j=0;j<Nsy;j++) {
        for(k=0;k<Nsz;k++) {
          section[i][j][k][0] = 0;
        }
      }
    }

    for(ip=0;ip<NNp;ip++) {
      i = (int)(Nsx*xx[ip]/xMax); if (i>=Nsx) i = Nsx-1;
      j = (int)(Nsy*yy[ip]/yMax); if (j>=Nsy) j = Nsy-1;
      k = (int)(Nsz*zz[ip]/zMax); if (k>=Nsz) k = Nsz-1;
      iq = section[i][j][k][0]+1;
      if (iq<NNs) {
        section[i][j][k][0] = iq;
        section[i][j][k][iq] = ip;
      }
    }
  }

  int maxSection() {
    int i,j,k,m;

    m = 0;
    for(i=0;i<Nsx;i++) {
      for(j=0;j<Nsy;j++) {
        for(k=0;k<Nsz;k++) {
          if (section[i][j][k][0]>m) m = section[i][j][k][0];
        }
      }
    }
    return ( m );
  }

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

}

/* real-coded lattice gas model 3D
 *
 * method
 *
 * (1) movement
 *  r = r + v*dt
 *
 * (2) collision - all particles collide each other in the same cell
 *  vm = mean v in the cell
 *  vi = vm + rotate(theta,fai)(vi0-vm)
 *
 */