/* template - direct simulation Monte Carlo method duasi-2D (3D) */
/* coded by Ikeuchi Mitsuru   */
/* ver 0.0.1   2005.11.19     */

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

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

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

  Button bt_reset, bt_startStop;
  Choice ch_kind, ch_kind2, ch_col, ch_view;
  Scrollbar sc_ratio, sc_temp, sc_scale;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 50;

  double t = 0.0;
  double dt = 1.0e-7;/* s */

  double dx = 1.0e-3;
  double dy = 1.0e-3;
  double dz = 1.0e-3;
  double vvc = dx*dy*dz;
  int Nsx = 40;
  int Nsy = 40;
  int Nsz = 3;
  double xMax = Nsx*dx;
  double yMax = Nsy*dy;
  double zMax = Nsz*dz;
  double zDepth = yMax;

  double ffn = 3.0e10;
  double sigmaMax = 4.0*3.1416*(3.0e-10)*(3.0e-10);
  double vvmax2 = 2000.0;


  int dispMode = 0;
  int resetSW = 0;
  int started = 1;
  int incRate = 1;
  double dispScale = (300.0/yMax);
  double scale = 1.0;
  double scrVar = 50.0;
  int collisionMode = 1;
  int kindOfGas = 18; /* Xe */
  int kindOfGas2 = 2; /* He */
  int kind1Ratio = 50;
  double adsorptionEnergy = 0.1; /* eV */
  double setTemp = 300.0;

  int dispWidth = 630;
  int dispHeight = 400;

  int dgX, dgY, dgXb,dgYb;
  double cx = 0.5*xMax;
  double cy = 0.5*yMax;
  double cz = 0.5*zMax;
  double theta = -20.0*3.14/180.0;
  double fai = -30.0*3.14/180.0;
  double dtheta = 0.0*3.14/180.0;
  double pai = 3.1415926536;
  double cosTh = Math.cos(theta); 
  double sinTh = Math.sin(theta);
  double cosFi = Math.cos(fai);
  double sinFi = Math.sin(fai);

  int NNp = 40*40*3*30;
  int kind[] = new int [NNp];

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

  double gasData[][] = {
    /* mass(x1.66e-27 kg),  d(x1.0e-10 m), HSBcolor */
    { 29.0 , 4.19, 0.00 }, /*  0 dummy */
    {  2.02, 2.29, 0.01 }, /*  1 H2  */
    {  4.0 , 2.33, 0.05 }, /*  2 He  */
    { 16.05, 4.83, 0.95 }, /*  3 CH4 */
    { 17.04, 5.94, 0.90 }, /*  4 NH3 */
    { 20.18, 2.77, 0.10 }, /*  5 Ne  */
    { 28.01, 4.19, 0.85 }, /*  6 CO  */

    { 28.02, 4.17, 0.30 }, /*  7 N2  */
    { 30.01, 4.20, 0.80 }, /*  8 NO  */
    { 29.0 , 4.19, 0.33 }, /*  9 air */
    { 32.0 , 4.07, 0.40 }, /* 10 O2  */
    { 36.46, 5.76, 0.75 }, /* 11 HCl */
    { 39.95, 4.17, 0.15 }, /* 12 Ar  */

    { 44.01, 5.62, 0.70 }, /* 13 CO2 */
    { 44.02, 5.71, 0.65 }, /* 14 N2O */
    { 64.06, 7.16, 0.60 }, /* 15 SO2 */
    { 70.9 , 6.98, 0.55 }, /* 16 Cl2 */
    { 83.8 , 4.76, 0.20 }, /* 17 Kr  */
    {131.3 , 5.74, 0.25 }  /* 18 Xe  */
  };

  String gasStr[] = { 
    " ","H2","He","CH4","NH3","Ne","CO",
    "N2","NO","air","O2","HCl","Ar",
	"CO2","N2O","SO2","Cl2","Kr","Xe"
  };

  int section[][][][] = new int [Nsx][Nsy][Nsz][200];
  int cellAttribute[][][] = new int [Nsx][Nsy][Nsz];

  double density[][][] = new double [Nsx][Nsy][Nsz];
  double temperature[][][] = new double [Nsx][Nsy][Nsz];
  double pressure[][][] = new double [Nsx][Nsy][Nsz];

  double vSpace[][][] = new double [20][20][20];

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

  int NplotMax = 12000;
  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 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}
  };

  int xpts[] = new int[100];
  int ypts[] = new int[100];

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

  public void init() {

    resize(dispWidth,dispHeight);
    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_kind = new Choice();
    ch_kind.addItem("H2"); ch_kind.addItem("He"); ch_kind.addItem("CH4");
    ch_kind.addItem("NH3"); ch_kind.addItem("Ne"); ch_kind.addItem("CO");
    ch_kind.addItem("N2"); ch_kind.addItem("NO"); ch_kind.addItem("air");
    ch_kind.addItem("O2"); ch_kind.addItem("HCl"); ch_kind.addItem("Ar");
    ch_kind.addItem("CO2"); ch_kind.addItem("N2O"); ch_kind.addItem("SO2");
    ch_kind.addItem("Cl2"); ch_kind.addItem("Kr"); ch_kind.addItem("Xe");
    ch_kind.addItemListener(this);
    ch_kind.select("Xe");

    ch_kind2 = new Choice();
    ch_kind2.addItem("H2"); ch_kind2.addItem("He"); ch_kind2.addItem("CH4");
    ch_kind2.addItem("NH3"); ch_kind2.addItem("Ne"); ch_kind2.addItem("CO");
    ch_kind2.addItem("N2"); ch_kind2.addItem("NO"); ch_kind2.addItem("air");
    ch_kind2.addItem("O2"); ch_kind2.addItem("HCl"); ch_kind2.addItem("Ar");
    ch_kind2.addItem("CO2"); ch_kind2.addItem("N2O"); ch_kind2.addItem("SO2");
    ch_kind2.addItem("Cl2"); ch_kind2.addItem("Kr"); ch_kind2.addItem("Xe");
    ch_kind2.addItemListener(this);
    ch_kind2.select("He");

    ch_col = new Choice();
    ch_col.addItem("null");
    ch_col.addItem("rigid sphere");
    ch_col.addItemListener(this);
    ch_col.select("rigid sphere");

    ch_view = new Choice();
    ch_view.addItem("sample 1000");
    ch_view.addItem("sample 8000");
    ch_view.addItem("veloc. 1000");
    ch_view.addItem("veloc. 8000");
	ch_view.addItem("mass flow");
    ch_view.addItem("v space");
	ch_view.addItem("density");
    ch_view.addItem("temperature");
    ch_view.addItem("pressure");
    ch_view.addItem("n mean 2D");
    ch_view.addItem("T mean 2D");
    ch_view.addItem("P mean 2D");
    ch_view.addItem("n grid 2D");
    ch_view.addItem("T grid 2D");
    ch_view.addItem("P grid 2D");

    ch_view.addItemListener(this);
    ch_view.select("sample 1000");

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

    sc_temp= new Scrollbar(Scrollbar.HORIZONTAL,300,10,100,610);
    sc_temp.addAdjustmentListener(this);

    sc_scale= new Scrollbar(Scrollbar.HORIZONTAL,100,10,50,1010);
    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(ch_kind);
    pnl.add(ch_kind2);
    pnl.add(sc_ratio);
    pnl.add(sc_temp);
    pnl.add(sc_scale);
    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 itemStateChanged(ItemEvent ev){
    if (ev.getSource() == ch_kind) {
      kindOfGas = 1 + ch_kind.getSelectedIndex();
	  resetSW = 1;
    } else if (ev.getSource() == ch_kind2) {
      kindOfGas2 = 1 + ch_kind2.getSelectedIndex();
	  resetSW = 1;
    } else if (ev.getSource() == ch_view) {
      dispMode = ch_view.getSelectedIndex();
    } else if (ev.getSource() == ch_col) {
      collisionMode = ch_col.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_ratio) { 
      kind1Ratio = sc_ratio.getValue();
      resetSW = 1;
    } else if (ev.getSource() == sc_temp) { 
      setTemp = 1.0*sc_temp.getValue();
      setGasTemperature(setTemp,kindOfGas);
	  setGasTemperature(setTemp,kindOfGas2);
    } else if (ev.getSource() == sc_scale) { 
      scale = 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;
    }
    theta += 0.5*3.14/180.0*(dgX-dgXb);
    fai += 0.5*3.14/180.0*(dgY-dgYb);
    cosTh = Math.cos(theta); 
    sinTh = Math.sin(theta);
    cosFi = Math.cos(fai);
    sinFi = Math.sin(fai);
  }

  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 gx, gbx, gby;
    double sc;

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

    sc = scale*dispScale;
    gbx = 50;
    gby = 60;


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

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

	} else if (dispMode==2) {
	  velocity3DPlot(gbx, gby, sc, 1000, 10.0);

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

	} else if (dispMode==4) {
	  flow3DPlot(gbx, gby, sc, 0.1);

	} else if (dispMode==5) {
	  vSpace3DPlot(gbx+40, gby+60, sc, 0.02);

	} else if (dispMode==6) {
	  setDensity();
	  field3DPlot(density, gbx, gby, sc, 5.0e-22);

	} else if (dispMode==7) {
	  setDensity();
	  setTemperature();
	  field3DPlot(temperature, gbx, gby, sc, 0.001);

	} else if (dispMode==8) {
	  setDensity();
	  setTemperature();
	  setPressure();
	  field3DPlot(pressure, gbx, gby, sc, 0.1);

	} else if (dispMode==9) {
      setDensity();
	  field2DPlot(density, gbx, gby, sc, 5.0e-22);

	} else if (dispMode==10) {
	  setDensity();
	  setTemperature();
	  field2DPlot(temperature, gbx, gby, sc, 0.001);

	} else if (dispMode==11) {
	  setDensity();
	  setTemperature();
	  setPressure();
	  field2DPlot(pressure, gbx, gby, sc, 0.1);

	} else if (dispMode==12) {
      setDensity();
	  grid2DPlot(density, 10.0*5.0e-22);

	} else if (dispMode==13) {
	  setDensity();
	  setTemperature();
	  grid2DPlot(temperature, 20.0*0.001);

	} else if (dispMode==14) {
	  setDensity();
	  setTemperature();
	  setPressure();
	  grid2DPlot(pressure, 20.0*0.1);

	} else if (dispMode==15) {
	  ;

	}

    velocityDistributionPlot(400, 200, 0.02, kindOfGas, kindOfGas2);


    gOff.setColor(Color.black);
    gOff.drawString("t="+(int)(t*1.0e6*1000.0+0.5)/1000.0+" us",630/6*0+10,40);
    gOff.drawString(gasStr[kindOfGas]+"("+kind1Ratio+") : "+gasStr[kindOfGas2]+"("+(100-kind1Ratio)+")",630/6*2+5,40);

    gOff.drawString("setT="+(int)(setTemp*10.0+0.5)/10.0+" K",630/6*3+10,40);
    gOff.drawString("scale="+(int)(scale*100.0)/100.0+" ",630/6*4+10,40);
    gOff.drawString("view",630/6*5+10,40);

    gx = 430;
    gOff.drawString("Box="+(float)(xMax*1000.0)+"x"+(float)(yMax*1000.0)+"x"+(float)(zMax*1000.0)+"(mm)",gx,60);
    gOff.drawString("N="+NNp+" particles",gx,80);
    gOff.drawString("FN = "+(float)ffn+" /particless",gx,100);
    gOff.drawString("max Nc="+maxSection()+"",gx,120);

    gOff.drawString("sgmMax="+(float)sigmaMax+"(m^2)",gx,140);
    gOff.drawString("vmax="+(float)maxSpeed()+"(m/s)",gx,160);

    gOff.drawString("alive=",gx,180);
    gOff.setColor(Color.getHSBColor((float)gasData[kindOfGas][2],0.9f,0.7f));
    gOff.drawString(" "+totalNumberOf(kindOfGas)+", ",gx+40,180);
    gOff.setColor(Color.getHSBColor((float)gasData[kindOfGas2][2],0.9f,0.7f));
    gOff.drawString(" "+totalNumberOf(kindOfGas2)+" ",gx+110,180);
  }

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

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

    drawBox(gx, gy, sc, 0);

    nn = setPlotData(Nsample);

    for (i=0; i<nn; i++) {
	  j = srtz[0][i];
      if (j<Nsample && kind[j]!=0 || j>=Nsample) {
        ir = (int)(ballSize*plotSize[j]*scale*(0.8+0.2*rotPlotz[j]/zDepth));
        gOff.setColor(Color.getHSBColor((float)(plotColor[j]),0.9f,(float)(plotShade[j]+0.15*rotPlotz[j]/zDepth)));
        gOff.fillOval(gx+(int)(rotPlotx[j]*sc+0.5)-ir/2,gy+(int)(rotPloty[j]*sc+0.5)-ir/2, ir, ir);
      }
    }
    drawBox(gx, gy, sc, 1);
  }

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

    for (i=0; i<Nsample; i++) {
	  aki = kind[i]; if (aki<0) aki = -aki;
        plotx[i] = xx[i]; ploty[i] = yy[i]; plotz[i] = zz[i];
        plotColor[i] = gasData[aki][2];
        plotShade[i] = 0.8; if (kind[i]<0) plotShade[i] = 0.3;
        plotSize[i] = gasData[aki][1];
    }

    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] = 8.0;
            plotColor[i] = 0.85;
			if (cellAttribute[ix][iy][iz]==1) plotShade[i] = 0.7;
			if (cellAttribute[ix][iy][iz]==2) plotShade[i] = 0.1;
            if (cellAttribute[ix][iy][iz]==3) plotShade[i] = 0.3;
		    i += 1;
		  }
        }
      }
    }
    nn = i;

    rotate(nn);
    zSort(nn,nn);

    return (nn);
  }


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

    drawBox(gx, gy, sc, 0);

    nn = setPlotData(Nsample);

    for (i=0; i<nn; i++) {
	  j = srtz[0][i];
	  if (j<Nsample) {
	    if (kind[j]!=0) {
	      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.getHSBColor((float)(gasData[kind[j]][2]),0.9f,0.7f));
          gOff.drawLine(gx+(int)(rotPlotx[j]*sc+0.5),gy+(int)(rotPloty[j]*sc+0.5), gx+(int)(rotp2x*sc+0.5),gy+(int)(rotp2y*sc+0.5));
        }
      } else {
        ir = (int)(1.0*plotSize[j]*scale*(0.8+0.2*rotPlotz[j]/zDepth));
        gOff.setColor(Color.getHSBColor((float)(plotColor[j]),0.9f,(float)(plotShade[j]+0.15*rotPlotz[j]/zDepth)));
        gOff.fillOval(gx+(int)(rotPlotx[j]*sc+0.5)-ir/2,gy+(int)(rotPloty[j]*sc+0.5)-ir/2, ir, ir);
	  }
    }
    drawBox(gx, gy, sc, 1);
  }


  void flow3DPlot(int gx, int gy, double sc,double mag) {
    int i,j,ix,iy,iz,jx,jy,jz,ip,jp,nn,nnj,ir;
	double p2x,p2y,p2z,rotp2x,rotp2y,rotp2z,vmx,vmy,vmz,m;

    drawBox(gx, gy, sc, 0);

    i = 0;
    for (ix=0; ix<Nsx; ix++) {
      for (iy=0; iy<Nsy; iy++) {
        for (iz=0; iz<Nsz; iz++) {
          plotx[i] = (ix+0.5)*dx;
          ploty[i] = (iy+0.5)*dy;
          plotz[i] = (iz+0.5)*dz;
		  if (cellAttribute[ix][iy][iz]>0) {
			plotSize[i] = 8.0;
            plotColor[i] = 0.85;
			if (cellAttribute[ix][iy][iz]==1) plotShade[i] = 0.7;
			if (cellAttribute[ix][iy][iz]==2) plotShade[i] = 0.1;
            if (cellAttribute[ix][iy][iz]==3) plotShade[i] = 0.3;
		  }
		  i += 1;
        }
      }
    }
    nn = i;

    rotate(nn);
    zSort(nn,nn);

    for (i=0; i<nn; i++) {
	  j = srtz[0][i];

	  jx = j/(Nsy*Nsz); jy = (j-Nsy*Nsz*jx)/Nsz; jz = (j-Nsy*Nsz*jx)%Nsz;
	  if (cellAttribute[jx][jy][jz]==0) {
	    vmx = 0.0; vmy = 0.0; vmz = 0.0;
	    nnj = section[jx][jy][jz][0];
	    if (nnj>0) {
          for (ip=1; ip<=nnj; ip++) {
            jp = section[jx][jy][jz][ip];
			m = gasData[kind[jp]][0];
            vmx += m*vx[jp]; vmy += m*vy[jp]; vmz += m*vz[jp];
          }

	      p2x = plotx[j]+mag*vmx*dt;
	      p2y = ploty[j]+mag*vmy*dt;
	      p2z = plotz[j]+mag*vmz*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 (vmx>0) {
            gOff.setColor(Color.blue);
          }
          gOff.drawLine(gx+(int)(rotPlotx[j]*sc+0.5),gy+(int)(rotPloty[j]*sc+0.5), gx+(int)(rotp2x*sc+0.5),gy+(int)(rotp2y*sc+0.5));
        }
      } else if (cellAttribute[jx][jy][jz]>0) {
        ir = (int)(1.0*plotSize[j]*scale*(0.8+0.2*rotPlotz[j]/zDepth));
        gOff.setColor(Color.getHSBColor((float)(plotColor[j]),0.9f,(float)(plotShade[j]+0.15*rotPlotz[j]/zDepth)));
        gOff.fillOval(gx+(int)(rotPlotx[j]*sc+0.5)-ir/2,gy+(int)(rotPloty[j]*sc+0.5)-ir/2, ir, ir);
	  }
    }

    drawBox(gx, gy, sc, 1);
  }


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

    sc = 1.2*scl;
    drawBox(gx, gy, sc, 2);

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

    for (i=0; i<NNp; i++) {
      ix = (int)((vx[i]+1000.0)/100.0);
      iy = (int)((vy[i]+1000.0)/100.0);
      iz = (int)((vz[i]+1000.0)/100.0);
	  if (ix>=0 && ix<20 && iy>=0 &&iy<20 && iz>=0 && iz<20) {
	    vSpace[ix][iy][iz] += 1.0;
	  }
    }

    i = 0;
    for (ix=0; ix<30; ix++) {
      for (iy=0; iy<20; iy++) {
        for (iz=0; iz<20; iz++) {
          plotx[i] = (ix+0.5)*dx;
          ploty[i] = (iy+0.5)*dy;
          plotz[i] = (iz+0.5)*dz;
		  i += 1;
        }
      }
    }
    nn = 8000;

    rotate(nn);
    zSort(nn,nn);

    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)(6.0*ff);
	  if (ir>0) {
	    gOff.setColor(Color.getHSBColor((float)(0.66-0.66*ff),0.8f,(float)(0.3+0.4*rotPlotz[j]/zDepth)));
        gOff.fillOval(gx+(int)(rotPlotx[j]*sc+0.5)-ir/2,gy+(int)(rotPloty[j]*sc+0.5)-ir/2, ir, ir);
      }
	}

    drawBox(gx, gy, sc, 3);
  }

/*---- field 3D plot ---- */

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

    drawBox(gx, gy, sc, 0);

    i = 0;
    for (ix=0; ix<Nsx; ix++) {
      for (iy=0; iy<Nsy; iy++) {
        for (iz=0; iz<Nsz; iz++) {
          plotx[i] = (ix+0.5)*dx;
          ploty[i] = (iy+0.5)*dy;
          plotz[i] = (iz+0.5)*dz;
		  if (cellAttribute[ix][iy][iz]==0) {
		    plotSize[i] = 0.1*section[ix][iy][iz][0]; if (plotSize[i]>10.0) plotSize[i]=10.0;
		    ff = mag*fld[ix][iy][iz]; if (ff>0.999) ff = 0.999;
			plotColor[i] = 0.66-0.66*ff;
			plotShade[i] = 0.8;
		  } else if (cellAttribute[ix][iy][iz]>0) {
			plotSize[i] = 8.0;
            plotColor[i] = 0.85;
			if (cellAttribute[ix][iy][iz]==1) plotShade[i] = 0.7;
			if (cellAttribute[ix][iy][iz]==2) plotShade[i] = 0.1;
            if (cellAttribute[ix][iy][iz]==3) plotShade[i] = 0.3;
		  }
		  i += 1;
        }
      }
    }
    nn = i;

    rotate(nn);
    zSort(nn,nn);

    for (i=0; i<nn; i++) {
	  j = srtz[0][i];
      ir = (int)(1.0*plotSize[j]*scale*(0.8+0.2*rotPlotz[j]/zDepth));
	  if (ir>0) {
        gOff.setColor(Color.getHSBColor((float)(plotColor[j]),0.9f,(float)(plotShade[j]+0.15*rotPlotz[j]/zDepth)));
        gOff.fillOval(gx+(int)(rotPlotx[j]*sc+0.5)-ir/2,gy+(int)(rotPloty[j]*sc+0.5)-ir/2, ir, ir);
      }
    }

    drawBox(gx, gy, sc, 1);
  }

/* set field[][] */

  void setDensity() {
    int ix,iy,iz;

    for (ix=0; ix<Nsx; ix++) {
      for (iy=0; iy<Nsy; iy++) {
        for (iz=0; iz<Nsz; iz++) {
          density[ix][iy][iz] = section[ix][iy][iz][0]*ffn/vvc;
        }
      }
    }
  }

  void setTemperature() {
    int ix,iy,iz,ip,j,n;
    double vmx,vmy,vmz,tke;

    for (ix=0; ix<Nsx; ix++) {
      for (iy=0; iy<Nsy; iy++) {
        for (iz=0; iz<Nsz; iz++) {
		  n = section[ix][iy][iz][0];
		  temperature[ix][iy][iz] = 0.0;
		  if (n>10) {
		    vmx = 0.0; vmy = 0.0; vmz = 0.0;
		    for (ip=1; ip<=section[ix][iy][iz][0]; ip++) {
			  j = section[ix][iy][iz][ip];
			  vmx += vx[j]; vmy += vy[j]; vmz += vz[j];
		    }
            vmx = vmx/n; vmy = vmy/n; vmz = vmz/n;

		    tke = 0.0;
		    for (ip=1; ip<=section[ix][iy][iz][0]; ip++) {
			  j = section[ix][iy][iz][ip];
              tke += 0.5*(gasData[kind[j]][0]*1.66e-27)*((vx[j]-vmx)*(vx[j]-vmx)+(vy[j]-vmy)*(vy[j]-vmy)+(vz[j]-vmz)*(vz[j]-vmz));
		    }
			tke = tke/n;
			temperature[ix][iy][iz] = tke/(1.5*1.38e-23);
		  }
        }
      }
    }
  }

  void setPressure() {
    int ix,iy,iz,ip,n;
    double kB;

    kB = 1.38e-23;
    for (ix=0; ix<Nsx; ix++) {
      for (iy=0; iy<Nsy; iy++) {
        for (iz=0; iz<Nsz; iz++) {
		  pressure[ix][iy][iz] = density[ix][iy][iz]*kB*temperature[ix][iy][iz];
        }
      }
    }
  }


/* -------------------------------  box Box  -------------- */

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

    setWaku(mode);

    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 (mode%2==0) {
        if (boxp[i][0]==farPoint || boxp[i][1]==farPoint) { 
          gOff.setColor(Color.lightGray);
          gOff.drawLine(gx,gy, g2x, g2y);
        }
	  } else if (mode%2==1) {
        if (boxp[i][0]!=farPoint && boxp[i][1]!=farPoint) { 
          gOff.setColor(Color.black);
          gOff.drawLine(gx,gy, g2x, g2y);
        }
	  }
    }
  }

  void setWaku(int mode) {
    int i,ix,iy,iz;
    double a,llx,lly,llz;

    a = 0.5*dx;
	if (mode<2) {
      llx = xMax; lly = yMax; llz = zMax;
	} else {
      llx = 20.0*dx; lly = 20.0*dx; llz = 20.0*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+(llx+a)*(double)(ix);
          wky[i] = -0.5*a+(lly+a)*(double)(iy);
          wkz[i] = -0.5*a+(llz+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);
  }

/* ------------------------------- 3D display utilities -------------- */

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

    for (i=0; i<Nmax; i++) {
      srtz[0][i] = i;
      srtz[1][i] = (int)(32767.0*rotPlotz[i]/zMax);
	}
    for (i=0; i<Nmax; i++) {
      for (j=0; j<Ntimes-1; j++) {
        if (srtz[1][j]>srtz[1][j+1]) {
          w = srtz[0][j]; srtz[0][j]=srtz[0][j+1]; srtz[0][j+1]=w;
          w = srtz[1][j]; srtz[1][j]=srtz[1][j+1]; srtz[1][j+1]=w;
        }
      }
    }    
  }

  void rotate(int Nmax) {
    int i;

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

/*------------------- 2D plot --------------- */

  void field2DPlot(double fld[][][], int gx, int gy, double sc, double mag) {
    int ix,iy,iz,iss,ca;
    double ss,f;

    ss = sc*xMax/(double)Nsx;
    iss = (int)(ss+0.999);
    for (ix=0; ix<Nsx; ix++) {
      for (iy=0; iy<Nsy; iy++) {
		f = meanField(fld,ix,iy)*mag;
		if (f>0.999) f = 0.999;
        gOff.setColor(Color.getHSBColor((float)(0.66-0.66*f),0.9f,0.9f));
		ca = cellAttribute[ix][iy][Nsz/2];
		if (ca>0) {
          if (ca==1) {
            gOff.setColor(Color.getHSBColor(0.85f,0.9f,0.9f));
          } else if (ca==2) {
            gOff.setColor(Color.getHSBColor(0.85f,0.9f,0.2f));
		  } else if (ca==3) {
            gOff.setColor(Color.getHSBColor(0.85f,0.9f,0.5f));
		  }
		}
        gOff.fillRect(gx+(int)(ss*ix+0.5),gy+(int)(ss*iy+0.5), iss, iss);
      }
    }
  }

  void grid2DPlot(double fld[][][], double weight) {
    grid2DPlot11(fld,0.0);
    grid2DPlot11(fld,weight);
  }


  void grid2DPlot11(double fld[][][], double weight) {
    int i,j,k,gbx,gby;
    int cx,cy,cz;
    double z,f,px,py,pz,sc,col;

    col = 0.85; if (weight>0.0) col = 0.66;

    cx = Nsx/2;
    cy = Nsy/2;
    cz = 0;
    sc = 8.0;

    gbx = 50+(int)(cx/sc);
    gby = 60+(int)(cy/sc);

    for(j=0; j<Nsy; j += 1) {
      for (i=0; i<Nsx; i++) {
	    f = 0.0;
        if (weight>0.0) f = weight*meanField(fld,i,j);
        z = f + 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<Nsx-1; i++) {
        gOff.setColor(Color.getHSBColor((float)col,0.9f,0.9f));
        gOff.drawLine(xpts[i],ypts[i],xpts[i+1],ypts[i+1]);
      }
    }

    for(i=0; i<Nsx; i += 1) {
      for (j=0; j<Nsy; j++) {
	    f = 0.0;
        if (weight>0.0) f = weight*meanField(fld,i,j);
        z = f + 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<Nsy-1; j++) { 
        gOff.setColor(Color.getHSBColor((float)col,0.9f,0.9f));
        gOff.drawLine(xpts[j],ypts[j],xpts[j+1],ypts[j+1]);
      }
    }
  }

  double meanField(double fld[][][], int ix, int iy) {
    int iz;
    double f;

    f = 0.0;
    for (iz=0; iz<Nsz; iz++) {
      f += fld[ix][iy][iz];
    }
    return( (f/Nsz) );
  }

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

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

    setVelocityDistribution(kind1, kind2);

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

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

    for (i=0; i<200; i++) {
      n1 = (int)(yMag*vDistribution[1][i]);
      n2 = (int)(yMag*vDistribution[2][i]);
      c1 = gasData[kindOfGas][2];
      c2 = gasData[kindOfGas2][2];
      cm = 0.5*(c1+c2);
      if (n1>=n2 && n1>0) {
        gOff.setColor(Color.getHSBColor((float)cm,0.9f,0.6f));
        gOff.drawLine(gx+i,gy+hoganHight,gx+i,gy+hoganHight-n2);
        gOff.setColor(Color.getHSBColor((float)c1,0.9f,0.9f));
        gOff.drawLine(gx+i,gy+hoganHight-n2-1,gx+i,gy+hoganHight-n1);
      } else if (n1<n2 && n2>0)  {
        gOff.setColor(Color.getHSBColor((float)cm,0.9f,0.6f));
        gOff.drawLine(gx+i,gy+hoganHight,gx+i,gy+hoganHight-n1);
        gOff.setColor(Color.getHSBColor((float)c2,0.9f,0.9f));
        gOff.drawLine(gx+i,gy+hoganHight-n1-1,gx+i,gy+hoganHight-n2);
      }
    }
    drawHogan(gx,gy,hoganWidth,hoganHight,xDiv,yDiv);
    gOff.setColor(Color.black);
    gOff.drawString("0",gx-5,gy+hoganHight+15);
    gOff.drawString("2000",gx+hoganWidth-20,gy+hoganHight+15);
    gOff.drawString("velocity distribution (m/s)",gx+20,gy+hoganHight+30);
  }

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

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

    for(i=0;i<NNp;i++) {
      if (kind[i]>0) {
        v = Math.sqrt(vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i]);
        iv = (int)(v/10.0);
        if (iv<200) {
          vDistribution[0][iv] += 1;
          if (kind[i]==kind1) {
            vDistribution[1][iv] += 1;
          } else if (kind[i]==kind2) {
            vDistribution[2][iv] += 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 -------- */

  int totalNumberOf(int gasKind) {
    int i,n;

    n = 0;
    for(i=0;i<NNp;i++) {
      if (kind[i]==gasKind) n += 1;
    }
    return(n);
  }

  double totalKineticEnergy() {
    int i;
    double tke;

    tke = 0.0;
    for(i=0;i<NNp;i++) {
	  if (kind[i]>0) {
        tke += 0.5*(gasData[kind[i]][0]*1.66e-27)*(vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i]);
      }
    }
    return( ffn*tke );
  }


/*-------------------------- set initial condition ------------ */
/*  cellAttribute -  0:space, 1:wall, 2:absorb, 3:adsorb        */
/*  kind>0 (free molec.), kind<0 (adsorbed), kind==0 (absorbed) */

  void setInitialCondition() {
    int i,ix,iy,iz,ir2;

    t = 0.0;

    for(ix=0;ix<Nsx;ix++) {
      for(iy=0;iy<Nsy;iy++) {
        for(iz=0;iz<Nsz;iz++) {
          cellAttribute[ix][iy][iz] = 0;

          ir2 = (iy-Nsy/2)*(iy-Nsy/2);
          if (ix==30 && ir2>10) {
            cellAttribute[ix][iy][iz] = 1;
          }
          if (ix==Nsx-1) {
		    cellAttribute[ix][iy][iz] = 2;
          }
        }
      }
    }

    for (i=0; i<NNp; i++) {
	  if (i%100<kind1Ratio) {
        kind[i] = kindOfGas;
	  } else {
	    kind[i] = kindOfGas2;
	  }
      xx[i] = xMax*Math.random()*0.74;
      yy[i] = yMax*Math.random();
      zz[i] = zMax*Math.random();
      vx[i] = 1000.0*(Math.random()-0.5);
      vy[i] = 1000.0*(Math.random()-0.5);
      vz[i] = 1000.0*(Math.random()-0.5);
    }

    setSigmaMax();
    setGasTemperature(setTemp,kindOfGas);
    setGasTemperature(setTemp,kindOfGas2);
  }

  void setSigmaMax() {
    int i;
    double d,dMax;

    dMax = 0.0;
    for(i=0;i<NNp;i++) {
      if (kind[i]>0) {
        d = gasData[kind[i]][1];
        if (d>dMax) dMax = d;
      }
    }
    sigmaMax = 1.0e-20*pai*dMax*dMax;
  }

  void setGasTemperature(double temp0, int knd) {
    int i,n;
    double m,tke,tmp,r;

    m = gasData[knd][0]*1.66e-27;
    n = 0; tke = 0.0;
    for(i=0;i<NNp;i++) {
      if (kind[i]==knd) {
	    n += 1;
        tke += 0.5*m*(vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i]);
	  }
    }
	tmp = tke/n/1.38e-23;
	r = Math.sqrt(temp0/(tmp+1.0e-20));
    for(i=0;i<NNp;i++) {
      if (kind[i]==knd) {
	    vx[i] = r*vx[i];
	    vy[i] = r*vy[i];
		vz[i] = r*vz[i];
	  }
    }
  }

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

  void timeEvolution() {
    int i;

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


  void timeStep() {

    t = t + dt;

    movement();
    dividingSection();
    if (collisionMode>0) {
      interaction();
    }
  }

/*-------------------------------  movement   -------------*/

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

    rr = 1.0;

    for (i=0; i<NNp; i++) {
      if (kind[i]>0) {
        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) {
          if (ca==1) {
            wallCell(i);
          } else if (ca==2) {
            absorbCell(i);
          } else if (ca==3) {
            adsorbCell(i);
          }
        }
      } else if (kind[i]<0) {
        if (Math.random()<1.0*Math.exp(-adsorptionEnergy/0.025)) {
		  kind[i] = -kind[i];
		  vx[i] = -vx[i]; vy[i] = -vy[i]; vz[i] = -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] = zz[i]+zMax;
      }
      if (zz[i] > zMax) {
        zz[i] = zz[i]-zMax;
      }
/*
      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)(Nsx*x/xMax); if (ix>=Nsx) ix = Nsx-1;
    iy = (int)(Nsy*y/yMax); if (iy>=Nsy) iy = Nsy-1;
    iz = (int)(Nsz*z/zMax); if (iz>=Nsz) iz = Nsz-1;
    return( cellAttribute[ix][iy][iz] );
  }

  void wallCell(int i) {
    int at100, at010, at001;

    xx[i] -= vx[i]*dt;
    yy[i] -= vy[i]*dt;
    zz[i] -= vz[i]*dt;

    at100 = cellAttributeAt(xx[i]+vx[i]*dt, yy[i], zz[i]);
    at010 = cellAttributeAt(xx[i], yy[i]+vy[i]*dt, zz[i]);
    at001 = cellAttributeAt(xx[i], yy[i], zz[i]+vz[i]*dt);
    if (at100==1) vx[i] = -vx[i];
    if (at010==1) vy[i] = -vy[i];
    if (at001==1) vz[i] = -vz[i];
	if (at100==0 && at010==0 && at001==0) {
      vx[i] = -vx[i]; vy[i] = -vy[i]; vz[i] = -vz[i];
	}
  }

  void absorbCell(int i) {
    kind[i] = 0;
  }

  void adsorbCell(int i) {
    kind[i] = -kind[i];
    xx[i] -= vx[i]*dt;
    yy[i] -= vy[i]*dt;
    zz[i] -= vz[i]*dt;
  }

/*-------------------------------  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++) {
      if (kind[ip]>0) {
        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;
		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 );
  }

/*-------------------------------  interaction   -------------*/

  void interaction() {
    int ic,jc,kc,ii,ir,i,j,nn,immc,pi,pj;
    double vmax,mmc;

    vvmax2 = 2.0*maxSpeed();

    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) {
            mmc = 0.5*nn*nn*ffn*sigmaMax*vvmax2*dt/vvc;
            immc = (int)mmc;
            for(ii=0;ii<=mmc;ii++) {
              collisionInTheCell(ic, jc, kc);
            }
            if ( Math.random()<(mmc-immc) ) {
              collisionInTheCell(ic, jc, kc);
		    }
          }
        }
      }
    }
  }

  double maxSpeed() {
    int i;
    double v2,vm2;

    vm2 = 0.0;
    for(i=0;i<NNp;i++) {
      v2 = vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i];
      if (v2>vm2) vm2 = v2;
    }
    return( Math.sqrt(vm2) );
  }

  void collisionInTheCell(int ic, int jc, int kc) {
    int i,j,pi,pj,nnc;

    nnc = section[ic][jc][kc][0];
    do {
      i = (int)(nnc*Math.random()); if (i>=nnc) i = nnc-1;
      j = (int)(nnc*Math.random()); if (j>=nnc) j = nnc-1;
    } while (i==j);
    pi = section[ic][jc][kc][i+1]; pj = section[ic][jc][kc][j+1];
    collision(pi,pj);
  }

  void collision(int i, int j) {
    double vrel,sgm,a,b,vcmx,vcmy,vcmz,th,costh,sinth,phi,cosphi,sinphi;

    vrel = Math.sqrt((vx[i]-vx[j])*(vx[i]-vx[j])+(vy[i]-vy[j])*(vy[i]-vy[j])+(vz[i]-vz[j])*(vz[i]-vz[j]));
	sgm = 1.0e-20*pai*0.25*(gasData[kind[i]][1]+gasData[kind[j]][1])*(gasData[kind[i]][1]+gasData[kind[j]][1]);
    if ( vrel/vvmax2 > Math.random() && (sgm/sigmaMax) > Math.random() ) {
      a = gasData[kind[i]][0]/(gasData[kind[i]][0]+gasData[kind[j]][0]);
      b = 1.0 - a;
      vcmx = a*vx[i]+b*vx[j]; vcmy = a*vy[i]+b*vy[j]; vcmz = a*vz[i]+b*vz[j];
      th = pai*(2.0*Math.random()-1);
      costh = Math.cos(th); sinth = Math.sin(th);
      cosphi = 2.0*Math.random()-1.0; 
      sinphi = Math.sqrt(1.0-cosphi*cosphi);
      vx[i] = vcmx + vrel*costh*sinphi*b;
      vx[j] = vcmx - vrel*costh*sinphi*a;
      vy[i] = vcmy + vrel*sinth*sinphi*b;
      vy[j] = vcmy - vrel*sinth*sinphi*a;
      vz[i] = vcmz + vrel*cosphi*b;
      vz[j] = vcmz - vrel*cosphi*a;
    }
  }


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

}

/* - memo - direct simulation Monte Carlo 3D
 *
 *  cell size; Lc = 0.001 (m)
 *  mean free path; lambda ~ 0.002 (m)
 *  No. of representative particle in the cell; Nc ~ 30 (particles)
 *  maximum relative speed; Vmax ~ 1000*2 (m/s)
 *  time interval; dt = 0.2*Lc/Vmax = 1.0e-7 (s)
 *  cross section of the collision; sgmMax = pai*(4.19e-10)^2 = 5.5e-19 (m^2)
 *  number density; n = Nc*FN/Lc^3, lambda = 1/(n*sgm)
 *    n = 1/(lambda*sgmMax) = 1/(0.002*5.5e-19) = 9.1e20 (particles/m^3)
 *  represent numbers of the real molecule; FN
 *    FN = n*Lc^3/Nc = 9.1e20*1e-9/30 = 3.0e10 (particles)
 *  the number of collisions in the cell; Mc
 *    Mc = 0.5*Nc*Nc*FN*sgmMax*Vmax*dt/Lc^3 ~ 1.5(times)
 *    collision occer only if (V/Vmax>random[0-1] and sgm/sgmMax>random[0-1])
 *
 */