/** applet No. 1106
 *
 * adsorption - direct simulation Monte Carlo method 2D
 * - multi thread  MorseMD :: display - asynchronous
 *
 * Created by Ikeuchi Mitsuru on January 14 2007.
 * Copyright (c) 2007 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2007.01.14  created
 * ver 0.0.2   2007.06.16  improved code
 *
 */

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

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

  // ------------------------------------ preset field -----------

  Thread th = null; // for run()-paint() thread
  DSMC2D dsmc = null;
  DrawGraph2D dg2d = new DrawGraph2D();

  // for event
  Button bt_reset, bt_startStop, bt_step;
  Choice ch_kind, ch_kind2, ch_col, ch_view;
  Scrollbar sc_temp, sc_adsorption, sc_adsorpBT, sc_scale;

  // for off-paint buffer
  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int dgX, dgY, dgXb, dgYb;

  int sleepTime = 50;
  int dispMode = 3;
  int thCount = 0;


  // ------------------------------ applet main thread -----------

  public void init() {

    resize(630,460);
    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); 

    bt_step = new Button("step");
    bt_step.addActionListener(this);

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

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

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

    ch_view = new Choice();
    ch_view.add("particle");
    ch_view.add("collision");
    ch_view.add("velocity");
    ch_view.add("sample 1000");
    ch_view.add("veloc. 1000");
    ch_view.add("density");
    ch_view.add("mass flow");
    ch_view.add("v space");
    ch_view.add("temperature");
    ch_view.add("pressure");
    ch_view.addItemListener(this);
    ch_view.select("sample 1000");

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

    sc_adsorption = new Scrollbar(Scrollbar.HORIZONTAL,100,10,0,160);
    sc_adsorption.addAdjustmentListener(this);

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

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

    addMouseListener(this);  
    addMouseMotionListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(2,6,5,0));
    //pnl.add(new Label("  "));
    pnl.add(bt_reset);
    pnl.add(bt_startStop);
    pnl.add(bt_step);
    pnl.add(new Label("  "));
    pnl.add(new Label("  "));
    pnl.add(new Label("  "));

    pnl.add(sc_adsorption);
    pnl.add(sc_adsorpBT);
    pnl.add(new Label("  "));
    pnl.add(sc_temp);
    pnl.add(sc_scale);
    pnl.add(ch_view);
    add(pnl,"North");
  }


  public void start() {
    if (dsmc == null) {
      dsmc = new DSMC2D();
      dsmc.start();
    }
    if (th == null) {
      th = new Thread(this);
      th.start();
    }
  }

  public void stop() {
    if (th != null) th = null;
    if (dsmc != null) dsmc = null;
  }


  // ---------------------------------- event listener -----------

  public void itemStateChanged(ItemEvent ev){
    if (ev.getSource() == ch_kind) {
      dsmc.setKind1( 1 + ch_kind.getSelectedIndex() );
    } else if (ev.getSource() == ch_kind2) {
      dsmc.setKind2( 1 + ch_kind2.getSelectedIndex() );
    } else if (ev.getSource() == ch_view) {
      dispMode = ch_view.getSelectedIndex();
    } else if (ev.getSource() == ch_col) {
      dsmc.setCollisionMode( ch_col.getSelectedIndex() );
    }
  }

  public void actionPerformed(ActionEvent ev){
    if(ev.getSource() == bt_reset){ 
      dsmc.reset();
      thCount = 0;
    } else if (ev.getSource() == bt_startStop){
      if (dsmc.getStartSW()==0) {
        dsmc.setStartSW(1);
      } else {
        dsmc.setStartSW(0);
      }
    } else if (ev.getSource() == bt_step){
      if (dsmc.getStartSW()==0) {
        dsmc.setStepSW();
      }
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_temp) { 
      dsmc.setContTemp( 1.0*sc_temp.getValue() );

    } else if (ev.getSource() == sc_adsorption) { 
      dsmc.setAdsorptionEnergy( 0.001*sc_adsorption.getValue() );

    } else if (ev.getSource() == sc_adsorpBT) { 
      dsmc.setAdsorpBaseTemp( 1.0*sc_adsorpBT.getValue() );

    } else if (ev.getSource() == sc_scale) { 
      dg2d.setScale( 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;
    }
    dg2d.setViewDelta( (dgX-dgXb), (dgY-dgYb) );
  }



  // ========================= run() - paint() loop ==============

  public void run() {
    while (th != null) {
      try {
        Thread.sleep(sleepTime);
      } catch (InterruptedException e) { }
      thCount += 1;
      offPaint();
      repaint();
    }
  }

  public void update(Graphics g) {
    paint(g);
  }

  public void paint(Graphics g) {
    g.drawImage(imgOff,0,0,this);
  }

  private void offPaint() {
    if (dsmc != null) {
      dg2d.plotGraph2D(gOff, dsmc, dispMode, thCount);
    }
  }

}



// =============================== draw graph 2D class ===========

class DrawGraph2D {

  private Graphics gOff;
  private DSMC2D dsmc;

  private int NNp;
  private int Nsx;
  private int Nsy;
  private double xMax;
  private double yMax;

  private double dispScale;
  private double scale = 1.0;

  private double xShift = 0.0;
  private double yShift = 0.0;

  private int vDistribution[][] = new int [200][3];
  private int vSpace[][] = new int [250][250];

  private int drawnCount = 0;


  // ------------------------------- class constructor -----------


  DrawGraph2D() {
  }

  // ----------------------------------- class methods -----------

  public void plotGraph2D(Graphics gg, DSMC2D adsmc, int dispMode, int thCount) {
    int gx, gbx, gby;
    double sc;
    int dsmcCount;

    gOff = gg;
    dsmc = adsmc;

    if (dsmc==null) return;
    dsmcCount = dsmc.getCount();
    if (dsmcCount==drawnCount) return;
    drawnCount = dsmcCount;

    NNp = dsmc.getNNp();
    Nsx = dsmc.getNsx();
    Nsy = dsmc.getNsy();
    xMax = dsmc.getxMax();
    yMax = dsmc.getyMax();
    dispScale = 300.0/yMax;

    gOff.setColor(Color.white);
    gOff.fillRect(0,0,630,460);

    sc = scale*dispScale;
    gbx = 20+(int)(sc*xShift);
    gby = 100+(int)(sc*yShift);

    if (dispMode==0) {
      dsmc.setStatMode(0);
      drawBlock(gbx,gby,sc);
      ovalPlot(gbx,gby,sc,NNp);

    } else if (dispMode==1) {
      dsmc.setStatMode(0);
      drawBlock(gbx,gby,sc);
      markPlot(gbx,gby,sc);

    } else if (dispMode==2) {
      dsmc.setStatMode(0);
      drawBlock(gbx,gby,sc);
      velosityPlot(gbx,gby,sc,NNp,20.0);

    } else if (dispMode==3) {
      dsmc.setStatMode(0);
      drawBlock(gbx,gby,sc);
      ovalPlot(gbx,gby,sc,1000);

    } else if (dispMode==4) {
      dsmc.setStatMode(0);
      drawBlock(gbx,gby,sc);
      velosityPlot(gbx,gby,sc,1000,20.0);

    } else if (dispMode==5) {
      dsmc.setStatMode(1);
      densityPlot(gbx,gby,sc,1.0);
      drawBlock(gbx,gby,sc);

    } else if (dispMode==6) {
      dsmc.setStatMode(2);
      drawBlock(gbx,gby,sc);
      massflowPlot(gbx,gby,sc,5.0);

    } else if (dispMode==7) {
      dsmc.setStatMode(0);
      vSpacePlot(gbx,gby,0.1);

    } else if (dispMode==8) {
      dsmc.setStatMode(1);
      temperaturePlot(gbx,gby,sc,0.001);
      drawBlock(gbx,gby,sc);

    } else if (dispMode==9) {
      dsmc.setStatMode(1);
      pressurePlot(gbx,gby,sc,1.0e8);
      drawBlock(gbx,gby,sc);

    } else if (dispMode==10) {
      ;

    }

    vDistributionPlot(420,220,0.1);

    gOff.setColor(Color.black);
    gOff.drawString("adsQ="+(int)(dsmc.getAdsorptionEnergy()*1000.0+0.5)/1000.0+"eV",630/6*0+10,70);
    gOff.drawString("baseT="+(int)(dsmc.getAdsorpBaseTemp()*10.0+0.5)/10.0+" K",630/6*1+10,70);
    gOff.drawString("contT="+(int)(dsmc.getContTemp()*10.0+0.5)/10.0+" K",630/6*3+10,70);
    gOff.drawString("scale="+(int)(scale*100.0+0.5)+"%",630/6*4+10,70);
    gOff.drawString("view",630/6*5+10,70);

    gOff.setColor(Color.black);
    gx = 430;
    gOff.drawString("t="+(int)(dsmc.getTime()*1.0e6*1000.0+0.5)/1000.0+" us",gx,100);
	gOff.drawString("Lx="+(float)(xMax*1000.0)+"(mm), Ly="+(float)(yMax*1000.0)+"(mm)",gx,120);
    gOff.drawString("N="+NNp+", FN="+(float)dsmc.getffn()+"",gx,140);
    gOff.drawString("Etot="+(float)totalKineticEnergy()+"(J)",gx,160);

    gOff.setColor(Color.black);
    gOff.drawString("adsorped atom ="+ (dsmc.getNNp()-numberOfPartickes(dsmc.getKind1()))+" ",gx,180);
/*
    gOff.setColor(Color.black);
    gOff.drawString("alive =",gx,180);
    gOff.setColor(Color.getHSBColor(dsmc.getColorOfKind(dsmc.getKind1()),0.9f,0.6f));
    gOff.drawString(" "+numberOfPartickes(dsmc.getKind1())+",",gx+60,180);
    gOff.setColor(Color.getHSBColor(dsmc.getColorOfKind(dsmc.getKind2()),0.9f,0.6f));
    gOff.drawString(" "+numberOfPartickes(dsmc.getKind2())+" ",gx+120,180);
*/
    gOff.setColor(Color.black);
    gOff.drawString("thread DSMC = "+dsmc.getCount()+" ",20,440);
    gOff.drawString("thread disp = "+thCount+" ",220,440);
  }


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

  private void ovalPlot(int gx, int gy, double sc,int Nsample) {
    int i,ir,knd;

    ir = 2;
    for (i=0; i<Nsample; i++) {
      knd = dsmc.getKind(i);
      if (knd>0) {
        gOff.setColor(Color.getHSBColor(dsmc.getColor(i),0.9f,0.9f));
        gOff.drawOval(gx+(int)(dsmc.getxx(i)*sc+0.5)-ir/2,gy+(int)(dsmc.getyy(i)*sc+0.5)-ir/2, ir, ir);
      } else if (knd<0) {
        gOff.setColor(Color.getHSBColor(dsmc.getColorOfKind(-knd),0.9f,0.25f));
        gOff.drawOval(gx+(int)(dsmc.getxx(i)*sc+0.5)-ir/2,gy+(int)(dsmc.getyy(i)*sc+0.5)-ir/2, ir, ir);
      }
    }
  }

  private void markPlot(int gx, int gy, double sc) {
    int i,ir;

    ir = 2;
    for (i=0; i<NNp; i++) {
      if (dsmc.getKind(i)!=0) {
        if (dsmc.getMark(i)==0) {
          gOff.setColor(Color.green);
        } else {
          gOff.setColor(Color.red);
        }
        gOff.drawOval(gx+(int)(dsmc.getxx(i)*sc+0.5)-ir/2,gy+(int)(dsmc.getyy(i)*sc+0.5)-ir/2, ir, ir);
      }
    }
  }

  private void velosityPlot(int gx, int gy, double sc, int Nsample, double mag) {
    int i,ix,iy,ix2,iy2;
    double dt;

    dt = dsmc.getdt();
    for (i=0; i<Nsample; i++) {
      if (dsmc.getKind(i)>0) {
        ix = gx+(int)(dsmc.getxx(i)*sc+0.5);
        iy = gy+(int)(dsmc.getyy(i)*sc+0.5);
        if (dsmc.getvx(i)>0.0) {
          gOff.setColor(Color.getHSBColor(0.66f,0.8f,0.8f));
        } else {
          gOff.setColor(Color.getHSBColor(0.01f,0.8f,0.8f));
        }

        ix2 = ix+(int)(dsmc.getvx(i)*dt*sc*mag+0.5);
        iy2 = iy+(int)(dsmc.getvy(i)*dt*sc*mag+0.5);
        gOff.drawLine(ix,iy,ix2,iy2);
      }
    }
  }

  //

  private void drawBlock(int gx, int gy, double sc) {
    int ix,iy,iss,ca;
    double ss;

    gOff.setColor(Color.black);
    gOff.drawRect(gx,gy,(int)(sc*xMax+0.5),(int)(sc*yMax+0.5));

    ss = sc*xMax/(double)Nsx;
    iss = (int)(ss+0.999);
	
    for (ix=0; ix<Nsx; ix++) {
      for (iy=0; iy<Nsy; iy++) {
        ca = dsmc.getCellAttribute(ix,iy);
        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);
        }
      }
    }
  }

  // field plot

  private void densityPlot(int gx, int gy, double sc, double mag) {
    int ix,iy,iss;
    double ss,d,a;

	a = dsmc.getffn()/((xMax/Nsx)*(yMax/Nsy));
    ss = sc*xMax/(double)Nsx;
    iss = (int)(ss+0.999);
    for (ix=0; ix<Nsx; ix++) {
      for (iy=0; iy<Nsy; iy++) {
        d = 0.025*mag*dsmc.getDensity(ix,iy)/a; if (d>0.999) d = 0.999;
        gOff.setColor(Color.getHSBColor(0.33f,0.8f,(float)(d)));
        gOff.fillRect(gx+(int)(ss*ix+0.5),gy+(int)(ss*iy+0.5), iss, iss);
      }
    }
  }

  private void temperaturePlot(int gx, int gy, double sc, double mag) {
    int ix,iy,iss;
    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 = dsmc.getTemperature(ix,iy)*mag; if (f>0.999) f = 0.999;
        gOff.setColor(Color.getHSBColor((float)(0.66-0.66*f),0.9f,0.9f));
        gOff.fillRect(gx+(int)(ss*ix+0.5),gy+(int)(ss*iy+0.5), iss, iss);
      }
    }
  }

  private void pressurePlot(int gx, int gy, double sc, double mag) {
    int ix,iy,iss;
    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 = dsmc.getPressure(ix,iy)*mag; if (f>0.999) f = 0.999;
        gOff.setColor(Color.getHSBColor((float)(0.66-0.66*f),0.9f,0.9f));
        gOff.fillRect(gx+(int)(ss*ix+0.5),gy+(int)(ss*iy+0.5), iss, iss);
      }
    }
  }


  private void massflowPlot(int gx, int gy, double sc,double mag) {
    int ix,iy,iss,igx,igy,igx2,igy2;
    double ss,mvx,mvy,dt;
	
    ss = sc*xMax/(double)Nsx;
    iss = (int)(ss+0.999);
    dt = dsmc.getdt();

    for (ix=0; ix<Nsx; ix++) {
      for (iy=0; iy<Nsy; iy++) {
        mvx = dsmc.getmfx(ix,iy); mvy = dsmc.getmfy(ix,iy);

        igx = gx+(int)(ss*(ix+0.5)+0.5);
        igy = gy+(int)(ss*(iy+0.5)+0.5);

        if (mvx>0.0) {
          gOff.setColor(Color.getHSBColor(0.66f,0.8f,0.8f));
        } else {
          gOff.setColor(Color.getHSBColor(0.01f,0.8f,0.8f));
        }
        igx2 = igx+(int)(mvx*dt*sc*mag+0.5);
        igy2 = igy+(int)(mvy*dt*sc*mag+0.5);
        gOff.drawLine(igx,igy,igx2,igy2);
      }
    }
  }


  private void fieldPlot(double fld[][], int gx, int gy, double sc, double mag) {
    int ix,iy,iss;
    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 = 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));
        gOff.fillRect(gx+(int)(ss*ix+0.5),gy+(int)(ss*iy+0.5), iss, iss);
      }
    }
  }

  // velocity space plot

  private void vSpacePlot(int gx, int gy, double mag) {
    int i,j;
    double f;

    setVSpace();

    gOff.setColor(Color.getHSBColor(0.66f,0.9f,0.9f));
    gOff.fillRect(gx,gy,250,250);

    for (i=0; i<250; i++) {
      for (j=0; j<250; j++) {
        f = vSpace[i][j]*mag; if (f>0.999) f = 0.999;
        if (f>0.0) {
          gOff.setColor(Color.getHSBColor((float)(0.66-0.66*f),0.9f,0.9f));
          gOff.fillRect(gx+i,gy+j,1,1);
        }
      }
    }
    gOff.setColor(Color.yellow);
    gOff.drawLine(gx,gy+125,gx+250,gy+125);
    gOff.drawLine(gx+125,gy,gx+125,gy+250);
  }

  private void setVSpace() {
    int i,j,ip;
    double v;

    for(i=0;i<250;i++) {
      for(j=0;j<250;j++) {
        vSpace[i][j] = 0;
      }
    }

    for(ip=0;ip<NNp;ip++) {
      if (dsmc.getKind(ip)>0) {
        i = (int)((dsmc.getvx(ip)+1250.0)/10.0+0.5);
        if (i<0) { 
          i = 0;
        } else if (i>249) {
          i = 249;
        }
        j = (int)((dsmc.getvy(ip)+1250.0)/10.0+0.5);
        if (j<0) { 
          j = 0;
        } else if (j>249) {
          j = 249;
        }

        vSpace[i][j] += 1;
      }
    }
  }

  // velocity distribution 

  private void vDistributionPlot(int gx, int gy, double yMag) {
    int i,n1,n2;
    double c1,c2,cm;

    setVelocityDistribution();

    for (i=0; i<200-1; i++) {
      n1 = (int)(yMag*vDistribution[i][1]);
      n2 = (int)(yMag*vDistribution[i][2]);
      c1 = dsmc.getColorOfKind(dsmc.getKind1());
      c2 = dsmc.getColorOfKind(dsmc.getKind2());
      cm = 0.5*(c1+c2);
      if (n1>=n2) {
        gOff.setColor(Color.getHSBColor((float)cm,0.9f,0.6f));
        gOff.drawLine(gx+i,gy+200,gx+i,gy+200-n2);
        gOff.setColor(Color.getHSBColor((float)c1,0.9f,0.9f));
        gOff.drawLine(gx+i,gy+200-n2-1,gx+i,gy+200-n1);
      } else {
        gOff.setColor(Color.getHSBColor((float)cm,0.9f,0.6f));
        gOff.drawLine(gx+i,gy+200,gx+i,gy+200-n1);
        gOff.setColor(Color.getHSBColor((float)c2,0.9f,0.9f));
        gOff.drawLine(gx+i,gy+200-n1-1,gx+i,gy+200-n2);
      }
    }
    drawHogan(gx, gy, 200, 200, 50, 100);
    gOff.setColor(Color.black);
    gOff.drawString("0",gx-5,gy+220);
    gOff.drawString("1000",gx+100-15,gy+220);
    gOff.drawString("(m/s)",gx+200-30,gy+220);
    gOff.drawString("velocity distribution",gx+20,gy+235);
  }

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

  private void setVelocityDistribution() {
    int i,iv,knd;
    double vx,vy,v;

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

    for(i=0;i<NNp;i++) {
      knd = dsmc.getKind(i);
      if (knd>0) {
        vx = dsmc.getvx(i); vy = dsmc.getvy(i);
        v = Math.sqrt(vx*vx+vy*vy);
        if (v>1999.0) v = 1999.0;
        iv = (int)(v/10.0);
        vDistribution[iv][0] += 1;
        if (knd==dsmc.getKind1()) {
          vDistribution[iv][1] += 1;
        } else if (knd==dsmc.getKind2()) {
          vDistribution[iv][2] += 1;
        }
      }
    }
  }


  // ----------------------------------------- utility -----------

  private double totalKineticEnergy() {
    int i;
    double tke;

    tke = 0.0;
    for(i=0;i<NNp;i++) {
      if (dsmc.getKind(i)>0) {
        tke += dsmc.getKineticEnergy(i);
      }
    }
    return( dsmc.getffn()*tke );
  }

  private int numberOfPartickes(int knd) {
    int i,n;

    n = 0;
    for(i=0;i<NNp;i++) {
      if (dsmc.getKind(i)==knd) {
        n += 1;
      }
    }
    return( n );
  }


  // ------------------------------------- I/O methods -----------

  public void setViewDelta(double xs, double ys) {
    xShift += (1.0/scale/dispScale)*xs;
    yShift += (1.0/scale/dispScale)*ys;
  }
  public void setScale(double scl) { scale = scl; }
}



// =========== direct simulation Monte Carlo method 2D ===========

class DSMC2D extends Thread {

  private double t = 0.0;
  private static final double dt = 1.0e-7; // s

  private static final double dx = 1.0e-3; // m 
  private static final double dy = 1.0e-3; // m 

  private static final int Nsx = 40;
  private static final int Nsy = 40;
  private double xMax = Nsx*dx;
  private double yMax = Nsy*dy;

  private double ffn = 2.0e4;
  private double sigmaMax = 4.0*3.0e-10;
  private double vvmax2 = 2000.0;

  private double adsorptionEnergy = 0.1; // eV
  private double adsorpBaseTemp = 300.0; // K
  private int collisionMode = 1;
  private int statMode = 1;
  private int kindOfGas1 = 12; // Ar
  private int kindOfGas2 = 12; // dummy
  private double contTemp = 300.0;

  private int NNp = Nsx*Nsy*30;
  private int kind[] = new int [NNp];
  private double xx[] = new double [NNp];
  private double yy[] = new double [NNp];
  private double vx[] = new double [NNp];
  private double vy[] = new double [NNp];

  private int mark[] = new int [NNp];

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

    { 28.02, 4.17, 0.15 }, /*  7 N2  */
    { 30.01, 4.20, 0.25 }, /*  8 NO  */
    { 29.0 , 4.19, 0.35 }, /*  9 air */
    { 32.0 , 4.07, 0.45 }, /* 10 O2  */
    { 36.46, 5.76, 0.55 }, /* 11 HCl */
    { 39.95, 4.17, 0.30 }, /* 12 Ar  */

    { 44.01, 5.62, 0.70 }, /* 13 CO2 */
    { 44.02, 5.71, 0.75 }, /* 14 N2O */
    { 64.06, 7.16, 0.80 }, /* 15 SO2 */
    { 70.9 , 6.98, 0.85 }, /* 16 Cl2 */
    { 83.8 , 4.76, 0.40 }, /* 17 Kr  */
    {131.3 , 5.74, 0.50 }  /* 18 Xe  */
  };

  private static final String molecStr[] = {
    "H2","He","CH4","NH3","Ne","CO","N2","NO","air",
    "O2","HCl","Ar","CO2","N2O","SO2","Cl2","Kr","Xe"
  };
  private int maxKind = 18;

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

  double density[][] = new double [Nsx][Nsy];
  double temperature[][] = new double [Nsx][Nsy];
  double pressure[][] = new double [Nsx][Nsy];
  double mfx[][] = new double [Nsx][Nsy];
  double mfy[][] = new double [Nsx][Nsy];

  private int resetSW = 0;
  private int startSW = 1;
  private int stepSW = 1;
  private int sleepTime = 1;
  private int stopSleepTime = 100;
  private int count = 0;


  // ------------------------------- class constructor -----------

  DSMC2D() {
    setInitialCondition();
  }


  // -------------------------------------- thread run -----------

  public void run() {
    while (true) {
      count += 1;
      timeEvolution();
      if (sleepTime>0) {
        try {
          Thread.sleep(sleepTime);
        } catch (InterruptedException e) { }
      }
    }
  }


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

  // cellAttribute -  0:space, 1:wall, 2:absorb, 3:adsorb

  private void setInitialCondition() {
    int i,j,ix,iy;

    t = 0.0;

    for(i=0;i<Nsx;i++) {
      for(j=0;j<Nsy;j++) {
        cellAttribute[i][j] = 0; // space

        if (i==0 && (j>Nsy/2-10 && j<Nsy/2+10)) {
          cellAttribute[i][j] = 3; // adsorb
        }
      }
    }

    for (i=0; i<NNp; i++) {
	  kind[i] = kindOfGas1;
      xx[i] = xMax*(Math.random()*0.95+0.045);
      yy[i] = yMax*Math.random();

      vx[i] = 1000.0*(Math.random()-0.5);
      vy[i] = 1000.0*(Math.random()-0.5);
    }

    setSigmaMax();
    setTemperature(contTemp,kindOfGas1);
    setTemperature(contTemp,kindOfGas2);
  }

  private void setSigmaMax() {
    int i;
    double s,sMax;

    sMax = 0.0;
    for(i=0;i<NNp;i++) {
      if (kind[i]>0) {
        s = gasData[kind[i]][1]*1.0e-10;
        if (s>sMax) sMax = s;
      }
    }
    sigmaMax = 2.0*sMax;
  }

  private void setTemperature(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]);
      }
    }
	if (n==0) return;
    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];
      }
    }
  }


  // ---------------------------------- time evolution -----------

  private void timeEvolution(){
    int i;

    if (resetSW==1) {
      setInitialCondition();
      resetSW = 0;
    }

    if (startSW==1) {
      timeStep();

    } else {
      if (stepSW==1) {
        timeStep();
        stepSW = 0;
      }
      try {
        Thread.sleep(stopSleepTime);
      } catch (InterruptedException e) { }
    }
  }

  private void timeStep() {
    int i;

    t = t + dt;

    movement();
    dividingSection();

    for(i=0;i<NNp;i++) {
      mark[i] = 0;
    }
    if (collisionMode>0) {
      interaction();
    }

    if (statMode==1) {
      setField();
    } else if (statMode==2) {
      setMeanFlow();
    }
  }

  // movement

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

        ca = cellAttributeAt(xx[i],yy[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/(adsorpBaseTemp/11600.0))) {
          kind[i] = -kind[i];
          vx[i] = -vx[i]; vy[i] = -vy[i];
        }
      }
    }

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

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

    ix = (int)(Nsx*x/xMax); if (ix>=Nsx) ix = Nsx-1;
    iy = (int)(Nsy*y/yMax); if (iy>=Nsy) iy = Nsy-1;
    return( cellAttribute[ix][iy] );
  }

  private void wallCell(int i) {
    int at10, at01;

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

    at10 = cellAttributeAt(xx[i]+vx[i]*dt, yy[i]);
    at01 = cellAttributeAt(xx[i], yy[i]+vy[i]*dt);
    if (at10==1) vx[i] = -vx[i];
    if (at01==1) vy[i] = -vy[i];
    if (at10==0 && at01==0) {
      vx[i] = -vx[i]; vy[i] = -vy[i];
    }
  }

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

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

  // dividing section

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

    for(i=0;i<Nsx;i++) {
      for(j=0;j<Nsy;j++) {
        section[i][j][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;
        iq = section[i][j][0]+1;
        section[i][j][0] = iq;
        section[i][j][iq] = ip;
      }
    }
  }

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

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

  // interaction

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

    vvmax2 = 2.0*maxSpeed();

    vvc = dx*dy;
    for(ic=0;ic<Nsx;ic++) {
      for(jc=0;jc<Nsy;jc++) {
        nn = section[ic][jc][0];
        if (nn>1) {
          mmc = 0.5*nn*nn*ffn*sigmaMax*vvmax2*dt/vvc;
          immc = (int)mmc;
          for(ii=0;ii<immc;ii++) {
            collisionInTheCell(ic, jc);
          }
          if ( Math.random()<(mmc-immc) ) {
            collisionInTheCell(ic, jc);
          }
        }
      }
    }
  }

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

    nnc = section[ic][jc][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][i+1]; pj = section[ic][jc][j+1];
    collision(pi,pj);
  }


  private void collision(int i, int j) {
    double vrel,a,b,vcmx,vcmy,th,costh,sinth;

    vrel = Math.sqrt((vx[i]-vx[j])*(vx[i]-vx[j])+(vy[i]-vy[j])*(vy[i]-vy[j]));
    if ( (vrel/vvmax2)> Math.random() && ((gasData[kind[i]][1]+gasData[kind[j]][1])/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];
      th = 2.0*3.14159265*Math.random();
      costh = Math.cos(th); sinth = Math.sin(th);
      vx[i] = vcmx + vrel*costh*b;
      vx[j] = vcmx - vrel*costh*a;
      vy[i] = vcmy + vrel*sinth*b;
      vy[j] = vcmy - vrel*sinth*a;
      mark[i] = 1; mark[j] = 1;
    }
  }

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

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


  // ----------------------------------------- utility -----------

  private void setField() {
    int i,j,k,ip,n;
    double ke,vmx,vmy,vvc,temp;

    vvc = dx*dy;
    for(i=0;i<Nsx;i++) {
      for(j=0;j<Nsy;j++) {
        n = section[i][j][0];
        density[i][j] = 0.9*density[i][j] + 0.1*n*ffn/vvc;

        temp = 0.0;
        if (n>10) {
          vmx = 0.0; vmy = 0.0;
          for(ip=1;ip<=section[i][j][0];ip++) {
            k = section[i][j][ip];
            vmx += vx[k]; vmy += vy[k];
          }
          vmx = vmx/n;
          vmy = vmy/n;

          ke = 0.0; 
          for(ip=1;ip<=section[i][j][0];ip++) {
            k = section[i][j][ip];
            ke += 0.5*(gasData[kind[k]][0]*1.66e-27)*((vx[k]-vmx)*(vx[k]-vmx)+(vy[k]-vmy)*(vy[k]-vmy));
          }
          temp = ke/n/(1.0*1.38e-23);
        }
        temperature[i][j] = 0.9*temperature[i][j] + 0.1*temp;

        pressure[i][j] =  0.9*pressure[i][j] + 0.1*density[i][j]*1.38e-23*temperature[i][j];
      }
    }
  }

  private void setMeanFlow() {
    int ix,iy,ip,ii;
    double m,mvx,mvy;

    for (ix=0; ix<Nsx; ix++) {
      for (iy=0; iy<Nsy; iy++) {
        mvx = 0.0; mvy = 0.0;
        for (ip=1; ip<=section[ix][iy][0]; ip++) {
          ii = section[ix][iy][ip];
          m = gasData[kind[ii]][0]/29.0;
          mvx += m*vx[ii]; mvy += m*vy[ii];
        }
        mfx[ix][iy] = 0.9*mfx[ix][iy] + 0.1*mvx;
        mfy[ix][iy] = 0.9*mfy[ix][iy] + 0.1*mvy;
      }
    }
  }


  // ------------------------------------- I/O methods -----------

  // control
  public void reset() { resetSW = 1; }
  public void setStartSW(int sw) { startSW = sw; }
  public int getStartSW() { return startSW; }
  public void setStepSW() { stepSW = 1; }
  public void setStatMode(int m) { statMode = m; }
  public int getCount() { return count; }

  public String getMolecStr(int knd) { return molecStr[knd]; }
  public void setKind1(int knd) { kindOfGas1 = knd; resetSW = 1; }
  public int getKind1() { return kindOfGas1; }
  public void setKind2(int knd) { kindOfGas2 = knd; resetSW = 1; }
  public int getKind2() { return kindOfGas2; }

  public void setAdsorptionEnergy(double qq) { adsorptionEnergy = qq; }
  public double getAdsorptionEnergy() { return adsorptionEnergy; }
  public void setAdsorpBaseTemp(double tt) { adsorpBaseTemp = tt; }
  public double getAdsorpBaseTemp() { return adsorpBaseTemp; }

  public void setCollisionMode(int cm) { collisionMode = cm; }
  public int getCollisionMode() { return collisionMode; }
  public void setContTemp(double tt) {
    contTemp = tt;
    setTemperature(contTemp,kindOfGas1);
    setTemperature(contTemp,kindOfGas2);
  }
  public double getContTemp() { return contTemp; }

  // get data
  public double getTime() { return t; }
  public double getdt() { return dt; }
  public int getNNp() { return NNp; }
  public int getNsx() { return Nsx; }
  public int getNsy() { return Nsy; }
  public double getxMax() { return xMax; }
  public double getyMax() { return yMax; }
  public double getffn() { return ffn; }

  public int getMaxKind() { return maxKind; }
  public double getMassOfKind(int knd) { return gasData[knd][0]*1.66e-27; }
  public float getColorOfKind(int knd) { return (float)gasData[knd][2]; }

  public int getKind(int i) { return kind[i]; }
  public double getxx(int i) { return xx[i]; }
  public double getyy(int i) { return yy[i]; }
  public double getvx(int i) { return vx[i]; }
  public double getvy(int i) { return vy[i]; }
  public float getColor(int i) { return (float)gasData[kind[i]][2]; }
  public int getMark(int i) { return mark[i]; }
  public double getKineticEnergy(int i) {
    return 0.5*gasData[kind[i]][0]*1.66e-27*(vx[i]*vx[i]+vy[i]*vy[i]);
  }

  public int getCellAttribute(int ix, int iy) { return cellAttribute[ix][iy]; }
  public double getDensity(int ix, int iy) { return density[ix][iy]; }
  public double getTemperature(int ix, int iy) { return temperature[ix][iy]; }
  public double getPressure(int ix, int iy) { return pressure[ix][iy]; }
  public double getmfx(int ix, int iy) { return mfx[ix][iy]; }
  public double getmfy(int ix, int iy) { return mfy[ix][iy]; }
}


/** direct simulation Monte Carlo method 2D
 *
 * dt (s) ; time division
 * N (1/m^3) ; number density
 * sgm (m^2) ; cross section
 * lambda (m) ; mean free path
 * Lc (m) ; cell length
 * Vc (m^3) ; cell volume
 * Nc ; number of particles in the cell ( 20 - 40 )
 * FN = real particles / dsmc particle
 * Vmax (m/s) = max(|Vi-Vj| = maximum particle velosity x 2
 * Mc ; number of collisions in dt (s) period
 * Mc=0.5*Nc*Nc*FN*sgm*Vmax*dt/Vc
 *
 * -- numerical example --
 * 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; sgm = 4R = 4*3.0e-10 (m)
 * number density; n = Nc*FN/Lc^2, lambda = 1/(n*sgm)
 *  n = 1/(lambda*sgm) = 1/(0.002*1.2e-9) = 4.0e11 (particles/m^2)
 * represent numbers of the real molecule; FN
 *  FN = n*Lc^2/Nc = 4e11*1e-6/30 = 1.3e4 --> 2.0e4 (particles)
 * the number of collisions in the cell; Mc
 *  Mc = 0.5*Nc*Nc*FN*sgm*Vmax*dt/Lc^2 ~ 4(times)
 *
 * -- method --
 *
 * 1. dividing region ( Lc < lambda )
 *
 * 2. put dsmc particles ~Nc in the cell
 *
 * 3. time evolution
 *
 * (1) movement r = r + v*dt
 *
 * (2) dividing section 
 *
 * (3) collision
 * 3.1 try collision Mc=0.5*Nc*Nc*FN*sgm*Vmax*dt/Vc 
 * 3.2 select i,j particle in the cell randomly
 *   only if |vi-Vj|/Vmax > random[0-1]
 * 3.3 determine particle velocity vi , vj
 *  - energy, momentum conservation
 *  - angle - random
 * 
 */