/** applet No. 1026
 *
 * template - real-coded lattice gas model 2D
 * - multi thread  rlg2d :: display - synchronous
 * ( get position and velocity --> cycle steal )
 *
 * Created by Ikeuchi Mitsuru on September 24 2006.
 * Copyright (c) 2006-2007 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2006.09.24  created
 * ver 0.0.2   2006.10.06  added DrawGraph2D class
 * ver 0.0.3   2006.12.07  improved code
 * ver 0.0.4   2007.05.23  improved code
 *
 */

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

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

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

  Thread th = null; // for run()-paint() thread
  RealLatticeGas2D rlg = null;
  Notice ntc = new Notice();
  DrawGraph2D dg2d = new DrawGraph2D();

  // for event
  Button bt_reset, bt_startStop, bt_step;
  Choice ch_view;
  Scrollbar sc_rlgSleep, sc_scale;

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

  int sleepTime = 100; // ms
  int dispMode = 2;

  int dgX, dgY, dgXb, dgYb; // mouse

  int thCount = 0;


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

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

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

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

    sc_scale= new Scrollbar(Scrollbar.HORIZONTAL,100,10,50,210);
    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(bt_step);
    pnl.add(new Label("  "));
    pnl.add(sc_scale);
    pnl.add(ch_view);
    add(pnl,"North");

  }

  public void start() {

    if (rlg == null) {
      rlg = new RealLatticeGas2D(ntc);
      rlg.start();
    }
    if (th == null) {
      th = new Thread(this);
      th.start();
    }
  }

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


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

  public void itemStateChanged(ItemEvent ev){
    if (ev.getSource() == ch_view) {
      dispMode = ch_view.getSelectedIndex();
    }
  }

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

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_scale) { 
      dg2d.setViewScale( 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.setShift( (dgX-dgXb), (dgY-dgYb) );
  }



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

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

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

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

  public void offPaint() {
    if (rlg != null) {
      dg2d.plotGraph(gOff, rlg, dispMode, thCount);
    }
  }

}



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

class DrawGraph2D {

  private Graphics gOff;
  private RealLatticeGas2D rlg;

  private double dispScale = (320.0/(320*1.0));
  private double viewScale = 1.0;
  private double scale = dispScale*viewScale;

  private int xImageLoc = 30;
  private int yImageLoc = 80;

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

  private int nsx = 0;
  private int nsy = 0;
  private double xMax = 0.0;
  private double yMax = 0.0;

  private int drawnCount = 0;


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

  DrawGraph2D() {
  }


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

  public void plotGraph(Graphics gg, RealLatticeGas2D rlg2d, int dispMode, int thCount) {
    int rlgCount;
    double temp;

    gOff = gg;
    rlg = rlg2d;
    if (rlg==null) return;
    rlgCount = rlg.getCount();
    if (rlgCount==drawnCount) return;
    drawnCount = rlgCount;

    nsx = rlg.getNsx();
    nsy = rlg.getNsy();
    xMax = rlg.getxMax();
    yMax = rlg.getyMax();

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

    xImageLoc = 30+(int)(scale*xShift);
    yImageLoc = 80+(int)(scale*yShift);

    if (dispMode==0) {
      ovalPlot(8000);

    } else if (dispMode==1) {
      velosityPlot(8000, 10.0);

    } else if (dispMode==2) {
      meanDensityPlot();

    } else if (dispMode==3) {
      meanFlowPlot(20.0);

    } else if (dispMode==4) {
      meanFieldPlot(1); // mode==1 temperature

    } else if (dispMode==5) {
      meanFieldPlot(2); // mode==2 pressure

    } else if (dispMode==6) {
      vSpacePlot(8000);

    } else if (dispMode==7) {
      ;
	
    }

    velocityDistributionPlot(410,200,0.015);

    gOff.setColor(Color.black);
    gOff.drawString("t="+(int)(rlg.getTime()*10.0+0.5)/10.0+" ",630/6*0+10,40);
    if (rlg.getStartSW()==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);
    gOff.drawString("Box="+rlg.getNsx()+"x"+rlg.getNsy()+"",630/6*0+10,60);
    gOff.drawString("N="+rlg.getNumberOfParticles()+"",630/6*1+10,60);
    temp = rlg.getSystemTemperature();
    gOff.drawString("T="+(int)(temp*1000.0+0.5)/1000.0+"",630/6*2+10,60);
    gOff.drawString("nue="+(int)(rlg.getViscosityMin(temp)*1000.0+0.5)/1000.0+" ",630/6*3+10,60);

    gOff.drawString("loop count",400,100);
    gOff.drawString("thread RLG ="+rlg.getCount()+"",400,120);
    gOff.drawString("thread disp ="+thCount+"",400,140);
  }


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

  private void ovalPlot(int Nsample) {
    int i,ir,ix,iy;

    rlg.setStatMode(1);

    drawWaku();
    drawBlock();

    ir = (int)(2.0*scale+0.5);
    for (i=0; i<Nsample; i++) {
      ix = xImageLoc + (int)(rlg.getxx(i)*scale+0.5);
      iy = yImageLoc + (int)(rlg.getyy(i)*scale+0.5);
      gOff.setColor(Color.getHSBColor(0.66f,0.9f,0.9f));
      gOff.drawOval(ix-ir/2,iy-ir/2, ir, ir);
    }
  }

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

    rlg.setStatMode(1);
    dt = rlg.getdt();

    drawWaku();
    drawBlock();

    for (i=0; i<Nsample; i++) {
      ix = xImageLoc + (int)(rlg.getxx(i)*scale+0.5);
      iy = yImageLoc + (int)(rlg.getyy(i)*scale+0.5);
	  vvx = rlg.getvx(i);
      ix2 = ix+(int)(vvx*dt*scale*mag+0.5);
      iy2 = iy+(int)(rlg.getvy(i)*dt*scale*mag+0.5);
      if (vvx>0.0) {
        gOff.setColor(Color.getHSBColor(0.6f,0.9f,0.7f));
      } else {
        gOff.setColor(Color.getHSBColor(0.01f,0.9f,0.7f));
      }
      gOff.drawLine(ix,iy,ix2,iy2);
    }
  }

  private void meanDensityPlot() {
    int ir,ic,jc,ix,iy;
    double den;

    rlg.setStatMode(1);

    ir = (int)(4*scale+0.999);

    for (ic=0; ic<nsx; ic+=4) {
      for (jc=0; jc<nsy; jc+=4) {
        den = 0.5*rlg.getrho(ic/4,jc/4); if (den>0.99) den = 0.99;
        ix = xImageLoc + (int)(ic*scale);
        iy = yImageLoc + (int)(jc*scale);
        gOff.setColor(Color.getHSBColor(0.3f,0.9f,(float)den));
        gOff.fillRect(ix,iy,ir,ir);
      }
    }

    drawWaku();
    drawBlock();
  }

  private void meanFlowPlot(double mag) {
    int ic,jc,ix,iy,ix2,iy2;
    double dt, mvx, mvy, den;

    rlg.setStatMode(2);

    dt = rlg.getdt();

    drawWaku();
    drawBlock();

    for (ic=0; ic<nsx; ic+=4) {
      for (jc=0; jc<nsy; jc+=4) {
        mvx = rlg.getxVelocity(ic/4,jc/4);
        mvy = rlg.getyVelocity(ic/4,jc/4); 
        den = 0.8*rlg.getrho(ic/4,jc/4); if (den>0.9) den = 0.9;

        ix = xImageLoc + (int)((ic+4/2)*scale);
        iy = yImageLoc + (int)((jc+4/2)*scale);
        ix2 = ix+(int)(mvx*dt*scale*mag+0.5);
        iy2 = iy+(int)(mvy*dt*scale*mag+0.5);
        if (mvx>=0.0) {
          gOff.setColor(Color.getHSBColor(0.6f,(float)den,0.7f));
        } else {
          gOff.setColor(Color.getHSBColor(0.01f,(float)den,0.7f));
        }
        gOff.drawLine(ix,iy,ix2,iy2);
      }
    }
  }

  private void meanFieldPlot(int mode) {
    int ic,jc,ix,iy,ir;
    double ke,d,den;

    rlg.setStatMode(3);

    ir = (int)(4*scale+0.999);

    for (ic=0; ic<nsx; ic+=4) {
      for (jc=0; jc<nsy; jc+=4) {
        den = rlg.getrho(ic/4,jc/4); if (den>0.9) den = 0.9;
        ke = rlg.getKEnergy(ic/4,jc/4);
        d = 0.0;
        if (mode==1) { // temperature
          d = 1.5*ke; if (d>0.99) d = 0.99;
        } else if (mode==2) { // pressure
          d = 1.5*(ke*den); if (d>0.99) d = 0.99;
        }
        ix = xImageLoc + (int)(ic*scale);
        iy = yImageLoc + (int)(jc*scale);
        gOff.setColor(Color.getHSBColor((float)(0.6-0.6*d),(float)den,0.9f));
        gOff.fillRect(ix,iy,ir,ir);
      }
    }

    drawWaku();
    drawBlock();
  }

  // utility

  private void drawBlock() {
    int ix,iy,iss,ca;
    double ss,rt;

    ss = scale*xMax/(double)nsx;
    iss = (int)(ss+0.999);
	
    for (ix=0; ix<nsx; ix++) {
      for (iy=0; iy<nsy; iy++) {
        ca = rlg.getCellAttribute(ix,iy);
        if (ca>0) {
          rt = rlg.getCellRefTemp(ca);
          if (rt>0.0) {
            gOff.setColor(Color.getHSBColor(0.85f,0.9f,(float)(0.2+0.8*rt)));
          } else {
            gOff.setColor(Color.getHSBColor(0.85f,0.3f,0.9f));
          }
          gOff.fillRect(xImageLoc+(int)(ss*ix+0.5),yImageLoc+(int)(ss*iy+0.5), iss, iss);
        }
      }
    }
  }

  private void drawWaku() {
    int xLen,yLen;

    xLen = (int)(scale*xMax+0.5);
    yLen = (int)(scale*yMax+0.5);
    gOff.setColor(Color.black);
    gOff.drawRect(xImageLoc,yImageLoc,xLen,yLen);
  }

  // velocity space plot

  private void vSpacePlot(int Nsample) {
    int i,ix,iy,bx,by;
    double mag;

    rlg.setStatMode(1);

    mag = 50.0;
    bx = 150+xImageLoc; by = 150+yImageLoc;
    for (i=0; i<Nsample; i++) {
      ix = bx+(int)(rlg.getvx(i)*mag);
      iy = by+(int)(rlg.getvy(i)*mag);
      gOff.setColor(Color.getHSBColor(0.66f,0.9f,0.9f));
      gOff.drawLine(ix,iy, ix,iy);
    }
    gOff.setColor(Color.black);
    gOff.drawLine(bx,by-150, bx,by+150);
    gOff.drawLine(bx-150,by, bx+150,by);
  }

  // velocity distribution plot

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

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

    // v distribution plot
    for (i=0; i<200; i++) {
      n1 = (int)(yMag*rlg.getVDistribution(i));
      if (n1>0) {
        gOff.setColor(Color.lightGray);
        gOff.drawLine(gx+i,gy+hoganHight,gx+i,gy+hoganHight-n1);
      }
    }

    // draw hogan
    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);

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


  // ------------------------------------- I/O methods -----------
  // control
  public void setShift(double xDisplacement, double yDisplacement) {
    xShift += (1.0/viewScale/dispScale)*xDisplacement;
    yShift += (1.0/viewScale/dispScale)*yDisplacement;
  }

  public void setViewScale(double sc) {
    viewScale = sc;
    scale = dispScale*viewScale;
  }

  public double getViewScale() { return viewScale; }
}



// ========== notice ( RLG free run and notice -> display ) ======

class Notice { 
  private int flag = 1;

  public synchronized void attention() { 
    notifyAll(); 
    flag = 0; 
  }

  public synchronized void getup() {
    while (flag==1) { 
      try { 
        wait(); 
      } catch (InterruptedException e) { } 
    } 
    flag = 1; 
  } 
}



// ============================== Real Lattice Gas 2Dl ===========

class RealLatticeGas2D extends Thread {

  private Notice note;
  private StealPositionThread sp;
  private StealVelocityThread sv;

  private double t = 0.0;
  private double dt = 1.0;
  private double dx = 1.0;
  private double dy = 1.0;

  int Nsx = 320;
  int Nsy = 320;
  private double xMax = Nsx*dx;
  private double yMax = Nsy*dy;
  private int NinCell = 4;

  private double gravity = 0.0;

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

  private int NNs = 20;
  int section[][][] = new int [Nsx][Nsy][NNs];
  private int cellAttribute[][] = new int [Nsx][Nsy];

  int density[][] = new int [Nsx/4][Nsy/4];
  double velocity[][][] = new double [Nsx/4][Nsy/4][2];
  double kEnergy[][] = new double [Nsx/4][Nsy/4];
  int vDistribution[][] = new int [200][4];

  private double refTemp[] = { 0.0, 0.0, 0.0, 0.0 };

  public int statMode = 2;

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


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

  RealLatticeGas2D(Notice ntc) {
    note = ntc;
    setInitialCondition();
  }


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

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


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

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

    t = 0.0;
    count = 0;

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

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


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

  private void timeEvolution() {
    int i;

    if (resetSW==1) {
      t = 0.0;
      setInitialCondition();
      startSW = 1;
      resetSW = 0;
    }
	
    if (startSW==1) {
      timeStep();

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

  void timeStep() {
    int i;

    t = t + dt;

    sv = new StealVelocityThread(this);
    sv.start();
    movement();
    try {
      sv.join();
    } catch(Exception ex) {}

    dividingSection();

    sp = new StealPositionThread(this);
    sp.start();
    collision();
    try {
      sp.join();
    } catch(Exception ex) {}

  }

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

      ca = cellAttributeAt(xx[i],yy[i]);
      if (ca>0) {
        xx[i] -= vx[i]*dt;
        yy[i] -= vy[i]*dt;
        r = 1.0;
        if (refTemp[ca]>0.0) {
          temp = 0.5*(vx[i]*vx[i]+vy[i]*vy[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];
      }
    }

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

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

    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;
    return( cellAttribute[ix][iy] );
  }

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

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

  void collisionInTheCell(int ic, int jc) {
    int i,k,n;
    double cTh,sTh,vxm,vym,vxs,vys;

    cTh = 0.0; sTh = 1.0;
    if (Math.random()<0.5) {
      cTh = 0.0; sTh = -1.0;
    }

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

    for(i=1;i<=n;i++) {
      k = section[ic][jc][i];
      vxs = vx[k]-vxm;
      vys = vy[k]-vym;
      vx[k] = vxm +cTh*vxs +sTh*vys;
      vy[k] = vym -sTh*vxs +cTh*vys;
    }
  }

  // dividing section
  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++) {
      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;
      if (iq<NNs) {
        section[i][j][0] = iq;
        section[i][j][iq] = ip;
      }
    }
  }


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

  private 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]);
    }
    return( tke );
  }

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

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

  private double viscosityMin(double temp) {
    int i;
    double rho;

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

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


  // ------------------------------------- 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 int getCount() { return count; }

  public void setStatMode(int sm) { statMode = sm; }

  // get data
  public int getNumberOfParticles() { return NNp; }
  public int getNsx() { return Nsx; }
  public int getNsy() { return Nsy; }
  public double getxMax() { return xMax; }
  public double getyMax() { return yMax; }

  public double getTime() { return t; }
  public double getdt() { return dt; }

  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 double getrho(int i, int j) { return density[i][j]/(16.0*NinCell); }
  public double getxVelocity(int i, int j) { return velocity[i][j][0]; }
  public double getyVelocity(int i, int j) { return velocity[i][j][1]; }
  public double getKEnergy(int i, int j) { return kEnergy[i][j]; }
  public int getCellAttribute(int i, int j) { return cellAttribute[i][j]; }

  public double getCellRefTemp(int ca) { return refTemp[ca]; }
  public int getVDistribution(int i) { return vDistribution[i][0]; }

  public double getSystemTemperature() { return systemTemperature(); }
  public double getViscosityMin(double temp) { return viscosityMin(temp); }
  public int getMaxSection() { return maxSection(); }

}



// ======================== position cycle steal class ===========

class StealPositionThread extends Thread {
  RealLatticeGas2D rlg;

  StealPositionThread(RealLatticeGas2D rg) {
    rlg = rg;
  }

  public void run() {
    setDensity();
  }

  void setDensity() {
    int i,j,ii,jj,n,nn;

    for (i=0; i<rlg.Nsx; i+=4) {
      for (j=0; j<rlg.Nsy; j+=4) {
        n = 0;
        for (ii=i; ii<i+4; ii++) {
          for (jj=j; jj<j+4; jj++) {
            nn = rlg.section[ii][jj][0];
            n += nn;
          }
        }
        rlg.density[i/4][j/4] = n;
      }
    }
  }

}



// ======================== velocity cycle steal class ===========

class StealVelocityThread extends Thread {
  RealLatticeGas2D rlg;

  StealVelocityThread(RealLatticeGas2D rg) {
    rlg = rg;
  }

  public void run() {
    setVelocityDistribution();
	if (rlg.statMode>1) meanDensity();
  }

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

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

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

  void meanDensity() {
    int i,j,ii,jj,n,nn,ip,k;
    double mvx,mvy,rv2;

    for (i=0; i<rlg.Nsx; i+=4) {
      for (j=0; j<rlg.Nsy; j+=4) {
        n = 0; mvx = 0.0; mvy = 0.0;
        for (ii=i; ii<i+4; ii++) {
          for (jj=j; jj<j+4; jj++) {
            nn = rlg.section[ii][jj][0];
            n += nn;
            if (rlg.statMode>1) {
              for (ip=1; ip<=nn; ip++) {
                k = rlg.section[ii][jj][ip];
                mvx += rlg.vx[k]; mvy += rlg.vy[k];
              }
            }
          }
        }
        if (rlg.statMode>1) {
          mvx = mvx/(n+1.0e-20);
          mvy = mvy/(n+1.0e-20);
          rlg.velocity[i/4][j/4][0] = mvx;
          rlg.velocity[i/4][j/4][1] = mvy;
          if (rlg.statMode==3) {
            rv2 = 0.0;
            for (ii=i; ii<i+4; ii++) {
              for (jj=j; jj<j+4; jj++) {
                nn = rlg.section[ii][jj][0];
                for (ip=1; ip<=nn; ip++) {
                  k = rlg.section[ii][jj][ip];
                  rv2 +=  (rlg.vx[k]-mvx)*(rlg.vx[k]-mvx)+(rlg.vy[k]-mvy)*(rlg.vy[k]-mvy);
                }
              }
            }
            rlg.kEnergy[i/4][j/4] = 0.5*rv2/(n+1.0e-20);
          }
        }
      }
    }
  }

}


/* real-coded lattice gas model 2D
 *
 * 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)(vi0-vm),  (theta = +-pai/2 : minmum viscosity)
 *
 */