/** applet No. 1262
 *
 * ball convection
 * - another RLG algorithm
 * - real-coded lattice gas model 2D
 * - multi thread  KRLG2D :: display - asynchronous
 * - float precision
 *
 * Created by Ikeuchi Mitsuru on February 24 2008.
 * Copyright (c) 2008 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2008.02.24  created
 * ver 0.0.2   2008.04.29  DrawGraph2D.setDisplay() changed
 *
 */

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

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

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

	Thread th = null; // for run()-paint() thread
	Thread th_rlg = null; // for RLG2D() thread
	ARLG2D rlg = new ARLG2D();
	DrawGraph2D dg2d = new DrawGraph2D(rlg);

	// for event
	private Button bt_reset, bt_startStop, bt_step;
	private Choice ch_view;
	private Scrollbar sc_temp1, sc_temp2, sc_scale;
	private int dgX, dgY, dgXb, dgYb; // mouse

	// for off-paint buffer
	private Image imgOff;
	private Graphics gOff;

	private int sleepTime = 100; // ms
	private int thCount = 0;


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

	public void init() {
		Dimension dim;

		resize(dg2d.getDispWidth(),dg2d.getDispHeight());
		setBackground(Color.white);

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

		setGUIPannel();
	}


	public void start() {

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


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




	// ===================================== GUI control ===========

	private void setGUIPannel() {

		setChoice();
		setButton();
		setScrollbar();
		setMouse();

		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(sc_temp1);
		pnl.add(sc_temp2);
		pnl.add(sc_scale);
		pnl.add(ch_view);
		add(pnl,"North");
	}


	// choice
	private void setChoice() {
		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("temperature");
		dg2d.setDispMode( 4 );
	}

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


	// button
	private void setButton() {
		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); 
	}

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


	// scrollbar
	private void setScrollbar() {
		sc_temp1= new Scrollbar(Scrollbar.HORIZONTAL,100,10,20,110);
		sc_temp1.addAdjustmentListener(this);

		sc_temp2= new Scrollbar(Scrollbar.HORIZONTAL,20,10,20,110);
		sc_temp2.addAdjustmentListener(this);

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

	public void adjustmentValueChanged(AdjustmentEvent ev){
		if (ev.getSource() == sc_temp1) { 
			rlg.setRefTemp1( 0.01*sc_temp1.getValue() );
		} else if (ev.getSource() == sc_temp2) { 
			rlg.setRefTemp2( 0.01*sc_temp2.getValue() );
		} else if (ev.getSource() == sc_scale) { 
			dg2d.setViewScale( 0.01*sc_scale.getValue() );
		}
	}


	// mouse
	private void setMouse() {
		addMouseListener(this);  
		addMouseMotionListener(this);
	}

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

			dg2d.plotGraph(thCount);
			repaint();
			thCount += 1;

			try {
				Thread.sleep(sleepTime);
			} catch (InterruptedException e) { }
		}
	}


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


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

}




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

class DrawGraph2D {

	private ARLG2D rlg;
	private Graphics gOff;
	private int dispWidth = 630;
	private int dispHeight = 460;

	private int dispMode = 4; //temperature
	private double dispScale = (320.0/240.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 double vDistribution[] = new double [200];

	private int drawnCount = 0;


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

	DrawGraph2D(ARLG2D rlg2d) {
		rlg = rlg2d;
	}


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

	public void plotGraph(int thCount) {
		int rlgCount;
		double temp;

		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,dispWidth,dispHeight);

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

		if (dispMode==0) {
			rlg.setTempSW(0);
			ovalPlot(8000);

		} else if (dispMode==1) {
			rlg.setTempSW(0);
			velosityPlot(8000, 10.0);

		} else if (dispMode==2) {
			rlg.setTempSW(0);
			meanDensityPlot();

		} else if (dispMode==3) {
			rlg.setTempSW(0);
			meanFlowPlot(50.0);

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

		} else if (dispMode==5) {
			rlg.setTempSW(1);
			meanFieldPlot(2); // mode==2 pressure

		} else if (dispMode==6) {
			rlg.setTempSW(0);
			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("temp1="+(int)(rlg.getRefTemp1()*100.0+0.5)/100.0+" ",630/6*2+10,40);
		gOff.drawString("temp2="+(int)(rlg.getRefTemp2()*100.0+0.5)/100.0+" ",630/6*3+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.getNNp()+"",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("max sec="+rlg.getMaxSection()+" ",630/6*5+10,60);

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


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

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

		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;

		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;

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

		for (ic=0; ic<nsx; ic+=4) {
			for (jc=0; jc<nsy; jc+=4) {
				den = 0.5*rlg.meanDensity(ic,jc); 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;
		double v[];

		dt = rlg.getdt();

		drawWaku();
		drawBlock();

		for (ic=0; ic<nsx; ic+=4) {
			for (jc=0; jc<nsy; jc+=4) {
				v = rlg.meanVelocity(ic,jc);
				den = 0.8*v[0]; if (den>0.9) den = 0.9;
				mvx = v[1];
				mvy = v[2]; 

				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 (mvy>=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 den,ke,d;

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

		for (ic=0; ic<nsx; ic+=4) {
			for (jc=0; jc<nsy; jc+=4) {
				den = 0.8*rlg.meanDensity(ic,jc); if (den>0.99) den = 0.99;
				ke = 1.0*rlg.meanTemperature(ic,jc);

				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;

		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;

		setVelocityDistribution();

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

		// v distribution plot
		for (i=0; i<200; i++) {
			n1 = (int)(yMag*vDistribution[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);
	}


	private void setVelocityDistribution() {
		int i,iv,NNp;

		NNp = rlg.getNNp();

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

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


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

	public void setDisplay(Graphics gg) {
		gOff = gg;
	}

	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 void setDispMode(int m) { dispMode = m; }

	public int getDispWidth() { return dispWidth; }
	public int getDispHeight() { return dispHeight; }
}




// ================================= Real Lattice Gas 2D ===========

class ARLG2D implements Runnable {

	private float t = 0.0f;
	private float dt = 1.0f;
	private float dx = 1.0f;
	private float dy = 1.0f;

	private int Nsx = 240;
	private int Nsy = 240;
	private float xMax = Nsx*dx;
	private float yMax = Nsy*dy;
	private int NinCell = 8;

	private float gravity = 0.001f;

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

	private byte cellAttribute[][] = new byte [Nsx][Nsy];
	private byte cellSigma[][] = new byte [Nsx][Nsy];
	private short cellNumber[][] = new short [Nsx][Nsy];
	private float cellmvx[][] = new float [Nsx][Nsy];
	private float cellmvy[][] = new float [Nsx][Nsy];
	private float cellTemp[][] = new float [Nsx][Nsy];

	private float refTemp[] = { 0.0f, 1.0f, 0.2f, 0.0f };

	private int sectionReadyFlag = 0;

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


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

	ARLG2D() {
		setInitialCondition();
	}


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

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


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

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

		t = 0.0f;
		count = 0;

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

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


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

	private void timeEvolution() {
		int i;

		if (resetSW==1) {
			t = 0.0f;
			setInitialCondition();
			startSW = 1;
			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();
		collision();
	}


	// (1) move
	private void movement() {
		int i,ca;
		float 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.0f;
				if (refTemp[ca]>0.0f) {
					temp = 0.5f*(vx[i]*vx[i]+vy[i]*vy[i]);
					r = (float)Math.sqrt(refTemp[ca]/temp)*(1.0f-1.0f*(floatRandom()-0.5f));
				} else {
					r = 1.0f;
				}
				vx[i] = -r*vx[i];
				vy[i] = -r*vy[i];
			}
		}

		rr = 1.0f;
		for (i=0; i<NNp; i++) {
			if (xx[i] < 0.0f) { 
				xx[i] = 0.0f; 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.0f) { 
				yy[i] = 0.0f; 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(float x, float 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] );
	}


	// (2) collision
	private void collision() {
		int i,ix,iy,n;
		float cTh,sTh,vmx,vmy,vxs,vys;

		setCell();

		for(i=0;i<NNp;i++) {
			ix = (int)(xx[i]/dx); if (ix>=Nsx) ix = Nsx-1;
			iy = (int)(yy[i]/dy); if (iy>=Nsy) iy = Nsy-1;
			n = cellNumber[ix][iy];
			if (n<2) continue;

			cTh = 0.0f; sTh = 1.0f;
			if (cellSigma[ix][iy]==1) {
				cTh = 0.0f; sTh = -1.0f;
			}
			vmx = cellmvx[ix][iy]/n;
			vmy = cellmvy[ix][iy]/n;
	  
			vxs = vx[i]-vmx;
			vys = vy[i]-vmy;
			vx[i] = vmx +cTh*vxs +sTh*vys;
			vy[i] = vmy -sTh*vxs +cTh*vys;
		}
	}


	// 

	private void setCell() {
		int i,ix,iy,n;
		float m,vmx,vmy,vxs,vys;

		sectionReadyFlag = 0;

		for(ix=0;ix<Nsx;ix++) {
			for(iy=0;iy<Nsy;iy++) {
				cellSigma[ix][iy] = 0; if (Math.random()>0.5) cellSigma[ix][iy] = 1;
				cellNumber[ix][iy] = 0;
				cellmvx[ix][iy] = 0.0f;
				cellmvy[ix][iy] = 0.0f;
				cellTemp[ix][iy] = 0.0f;
			}
		}

		for(i=0;i<NNp;i++) {
			ix = (int)(xx[i]/dx); if (ix>=Nsx) ix = Nsx-1;
			iy = (int)(yy[i]/dy); if (iy>=Nsy) iy = Nsy-1;
			cellNumber[ix][iy] += 1;
			cellmvx[ix][iy] += vx[i];
			cellmvy[ix][iy] += vy[i];
		}

		if (tempSW==1) {
			for(i=0;i<NNp;i++) {
				ix = (int)(xx[i]/dx); if (ix>=Nsx) ix = Nsx-1;
				if (ix<0) ix = 0;
				iy = (int)(yy[i]/dy); if (iy>=Nsy) iy = Nsy-1;
				if (iy<0) iy = 0;

				n = cellNumber[ix][iy];
				if (n>0) {
					vmx = cellmvx[ix][iy]/cellNumber[ix][iy];
					vmy = cellmvy[ix][iy]/cellNumber[ix][iy];
					vxs = vx[i]-vmx; vys = vy[i]-vmy;
					cellTemp[ix][iy] += 0.5f*(vxs*vxs+vys*vys);
				}
			}
		}

		sectionReadyFlag = 1;
	}


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

	private float floatRandom() { return (float)Math.random(); }


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

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


	private float systemTemperature() {
		return( totalKineticEnergy()/NNp );
	}


	private float viscosityMin(double temp) {
		float rho;

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


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

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


	// mean field

	public double meanDensity(int ix, int iy) {
		int i,j,n;

		cellDataReady();
		n = 0;
		for (i=ix; i<ix+4; i++) {
			for (j=iy; j<iy+4; j++) {
				n += cellNumber[i][j];
			}
		}
		return n/(16.0*NinCell);
	}


	public double[] meanVelocity(int ix, int iy) {
		int i,j,n;
		float m,mvx,mvy; 
		double v[] = new double [4];

		cellDataReady();
		n = 0; m = 0.0f; mvx = 0.0f; mvy = 0.0f;
		for (i=ix; i<ix+4; i++) {
			for (j=iy; j<iy+4; j++) {
				n += cellNumber[i][j];
				mvx += cellmvx[i][j];
				mvy += cellmvy[i][j];
			}
		}
		if (n>0) {
			v[0] = n/(16.0*NinCell); v[1] = mvx/n; v[2] = mvy/n;
		} else {
			v[0] = 0.0; v[1] = 0.0; v[2] = 0.0;
		}
		return v;
	}


	public double meanTemperature(int ix, int iy) {
		int i,j,n;
		float s; 

		cellDataReady();
		n = 0; s = 0.0f;
		for (i=ix; i<ix+4; i++) {
			for (j=iy; j<iy+4; j++) {
				n += cellNumber[i][j];
				s += cellTemp[i][j];
			}
		}

		if (n==0) return 0.0;
		return (double)(s/n);
	}


	//

	private void cellDataReady() {
		while(sectionReadyFlag==0) {
			try {
				Thread.sleep(1);
			} catch (InterruptedException e) { }
		}
	}


	// ------------------------------------- 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 setTempSW(int sw) { tempSW = sw; }
	public int getTempSW() { return tempSW; }
	public int getCount() { return count; }

	public void setRefTemp1(double t1) { refTemp[1] = (float)t1; }
	public double getRefTemp1() { return refTemp[1]; }
	public void setRefTemp2(double t2) { refTemp[2] = (float)t2; }
	public double getRefTemp2() { return refTemp[2]; }

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

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

	public double getxx(int i) { return (double)xx[i]; }
	public double getyy(int i) { return (double)yy[i]; }
	public double getvx(int i) { return (double)vx[i]; }
	public double getvy(int i) { return (double)vy[i]; }
	public double getv2(int i) { return (double)(vx[i]*vx[i]+vy[i]*vy[i]); }

	public int getCellAttribute(int i, int j) { return cellAttribute[i][j]; }
	public double getCellRefTemp(int ca) { return (double)refTemp[ca]; }
	public double getSystemTemperature() { return (double)systemTemperature(); }
	public double getViscosityMin(double temp) { return (double)viscosityMin(temp); }
	public int getMaxSection() { return maxSection(); }

}

/* 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)
 *
 */