/** applet No. 1280
 *
 * wall - wind tunnel
 * - MAC(Marker and Cell) method 2D
 * - multi thread  MAC2D :: display - asynchronous
 *
 * Created by Ikeuchi Mitsuru on April 06 2008.
 * Copyright (c) 2008 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2008.04.06  created
 *
 */

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

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

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

	Thread th = null; // for run()-paint() thread
	Thread th_mac = null; // for MAC thread
	MarkerAndCell2D mac = new MarkerAndCell2D();
	DrawGraph2D dg2d = new DrawGraph2D(mac);

	// for event
	private Choice ch_view;
	private Button bt_reset, bt_startStop;
	private Scrollbar sc_Reynolds,sc_vx0;
	private int dgX, dgY, dgXb, dgYb; // mouse

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

	private int sleepTime = 50;
	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_mac == null) {
			th_mac = new Thread(mac);
			th_mac.start();
		}
	}

	public void stop() {
		if (th != null) th = null;
		if (th_mac != null) th_mac = 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_Reynolds);
		pnl.add(sc_vx0);
		pnl.add(new Label("  "));
		pnl.add(ch_view);
		add(pnl,"North");
	}

	// choice
	private void setChoice() {
		ch_view = new Choice();
		ch_view.add("flow");
		ch_view.add("pressure");
		ch_view.add("flow+P");
		ch_view.add("grid P");
		ch_view.add("grid v,P");
		ch_view.add("flow line");
		ch_view.addItemListener(this);
		ch_view.select("flow line");
	}

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

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


	// scrollbar
	private void setScrollbar() {
		sc_Reynolds = new Scrollbar(Scrollbar.HORIZONTAL,400,10,10,2010);
		sc_Reynolds.addAdjustmentListener(this);

		sc_vx0 = new Scrollbar(Scrollbar.HORIZONTAL,100,10,0,110);
		sc_vx0.addAdjustmentListener(this);
	}

	public void adjustmentValueChanged(AdjustmentEvent ev){
		if (ev.getSource() == sc_Reynolds) { 
			mac.setReynoldsNumber( (double)(sc_Reynolds.getValue()) );
		} else if (ev.getSource() == sc_vx0) { 
			mac.setFanFlow( 0.01*(double)(sc_vx0.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.setViewDelta( (dgX-dgXb), (dgY-dgYb) );
	}




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

	public void run() {
		while (th != null) {
			dg2d.plotGraph2d(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 MarkerAndCell2D mac;
	private Graphics gOff;
	private int dispWidth = 630;
	private int dispHeight = 460;

	private int dispMode = 5;
	private double dispScale = 4.0;
	private double viewScale = 1.0;
	private double scale = dispScale*viewScale;

	private int nnx;
	private int nny;
	private int xImageLoc = 20;
	private int yImageLoc = 80;

	private double viewTheta = -15.0*3.14/180.0;
	private double viewFai = -72.0*3.14/180.0;
	private double cosTh = Math.cos(viewTheta); 
	private double sinTh = Math.sin(viewTheta);
	private double cosFi = Math.cos(viewFai);
	private double sinFi = Math.sin(viewFai);

	private int xpts[] = new int[200];
	private int ypts[] = new int[200];
	private double fpts[] = new double[200];

	private int drawnCount = 0;


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

	DrawGraph2D(MarkerAndCell2D mac2d) {
		mac = mac2d;
	}


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

	public void plotGraph2d(int thCount) {
		int macCount,gb;
		double xBox, yBox;

		macCount = mac.getCount();
		if (macCount==drawnCount) return;
		drawnCount = macCount;

		nnx = mac.getNNx();
		nny = mac.getNNy();

		gOff.setColor(Color.white);
		gOff.fillRect(0,0,dispWidth,dispHeight);

		if (dispMode==0) {
			boundaryPlot();
			flowPlot(20.0);

		} else if (dispMode==1) {
			pressurePlot(5.0);
			boundaryPlot();

		} else if (dispMode==2) {
			pressurePlot(5.0);
			boundaryPlot();
			flowPlot(20.0);

		} else if (dispMode==3) {
			gridPressurePlot(10.0);

		} else if (dispMode==4) {
			gridFlowPlot(10.0, 3.0);

		} else if (dispMode==5) {
			boundaryPlot();
			flowPlot(20.0);
			flowLinePlot();

		} else if (dispMode==6) {
			;

		}

		gOff.setColor(Color.blue);
		gOff.drawString("t="+(int)(mac.getTime()*100.0+0.5)/100.0+" ",10,40);
		if (mac.getStartSW()==0) {
			gOff.drawString("stopped",630/6*1+10,40);
		} else {
			gOff.drawString("started",630/6*1+10,40);
		}
		gOff.drawString("Re="+(int)(mac.getReynoldsNumber()+0.5)+" ",630/6*2+10,40);
		gOff.drawString("fan="+(int)(mac.getFanFlow()*100.0+0.5)/100.0+" ",630/6*3+10,40);
		gOff.drawString("view",630/6*5+10,40);

		gOff.setColor(Color.black);
		xBox = mac.getNNx()*mac.getdx(); yBox = mac.getNNy()*mac.getdy();
		gOff.drawString("wind tunnel="+(int)(xBox*10.0+0.5)/10.0+"x"+(int)(yBox*10.0+0.5)/10.0+" ",10,60);

		//gOff.setColor(Color.blue);
		//gOff.drawString("pError="+(float)(mac.getPressureError())+" ",200,60);
		//gOff.setColor(Color.black);
		//gOff.drawString("thread MAC2D = "+mac.getCount()+" ",350,60);
		//gOff.drawString("disp = "+thCount+" ",500,60);
	}


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

	private void flowPlot(double vmag) {
		int i,j,ix,iy,ix2,iy2;
		double vmx,vmy;

		for (i=0; i<nnx; i++) {
			for (j=0; j<nny; j++) {
				ix = xImageLoc + (int)(scale*(i+0.5));
				iy = yImageLoc + (int)(scale*(j+0.5));

				vmx = mac.getvx(i,j); vmy = mac.getvy(i,j); 
				ix2 = ix+(int)(viewScale*vmag*vmx);
				iy2 = iy+(int)(viewScale*vmag*vmy);
				if (vmx>0) {
					gOff.setColor(Color.getHSBColor(0.70f,0.9f,0.8f));
				} else {
					gOff.setColor(Color.getHSBColor(0.85f,0.9f,0.8f));
				}
				gOff.drawLine(ix,iy, ix2,iy2);
			}
		}
	}

	private void pressurePlot(double mag) {
		int i,j,ix,iy;
		double c,absc,p0;

		p0 = mac.getMeanPressure();
		for (i=0; i<nnx; i++) {
			for (j=0; j<nny; j++) {
				ix = xImageLoc + (int)(scale*i);
				iy = yImageLoc + (int)(scale*j);
				c = mag*(mac.getPressure(i,j) - p0);
				if (c>0.999) c = 0.999;
				if (c<-0.999) c = -0.999;
				absc = c; if (absc<0) absc = -absc;
				gOff.setColor(Color.getHSBColor((float)(0.33-0.33*c), (float)(0.1+0.3*absc), 0.8f));
				gOff.fillRect(ix,iy, 4,4);
			}
		}
	}

	//

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

		for (i=0; i<nnx; i++) {
			for (j=0; j<nny; j++) {
				if (mac.getCellKind(i,j)==1) {
					ix = xImageLoc + (int)(scale*i);
					iy = yImageLoc + (int)(scale*j);
					gOff.setColor(Color.gray);
					gOff.fillRect(ix,iy, 4,4);
				}
			}
		}
	}

	//

	private void flowLinePlot() {
		int i;
		double sx,sy;

		gOff.setColor(Color.red);
		for (i=0; i<9; i++) {
			flowLine(1.0, 0.5+0.4*i, 500);
		}
	}

	private void flowLine(double x0, double y0, int tstep) {
		int i,j,ix,iy,ix2,iy2,px,py;
		double dt,dx,dy, sx,sy;

		dt = 2.0*mac.getdt();
		dx = mac.getdx(); dy = mac.getdy();

		sx = x0; sy = y0;
		for (i=0; i<tstep; i++) {
			ix = xImageLoc + (int)(scale*sx/dx+0.5);
			iy = yImageLoc + (int)(scale*sy/dy+0.5);
			px = (int)(sx/dx+0.5);if (px>=nnx) break;
			py = (int)(sy/dy+0.5);if (py>=nny) break;
			sx += mac.getvx(px,py)*dt;
			sy += mac.getvy(px,py)*dt;
			ix2 = xImageLoc + (int)(scale*sx/dx+0.5);
			iy2 = yImageLoc + (int)(scale*sy/dy+0.5);
			gOff.setColor(Color.red);
			gOff.drawLine(ix,iy, ix2,iy2);
		}
	}

	// grid pressure plot

	private void gridPressurePlot(double weight) {
		int i,j,k,gbx,gby;
		int cx,cy,cz;
		double z,f,px,py,pz,col,p0;

		cx = nnx/2;
		cy = nny/2;
		cz = 0;

		gbx = 20+(int)(cx/scale);
		gby = 60+(int)(cy/scale);

		for(j=0; j<nny; j += 1) {
			for (i=0; i<nnx; i++) {
				if (mac.getCellKind(i,j)==0) {
					f = weight*mac.getPressure(i,j);
				} else {
					f = 0.0;
				}
				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)(scale*px)+gbx;
				ypts[i] = (int)(scale*py)+gby;
				fpts[i] = f;
			}
			for (i=0; i<nnx-1; i++) {
				if (fpts[i]+fpts[i+1]>0.1) {
					gOff.setColor(Color.getHSBColor(0.10f,0.7f,0.9f));
				} else if (fpts[i]+fpts[i+1]<-0.1) {
					gOff.setColor(Color.getHSBColor(0.66f,0.7f,0.9f));
				} else {
					gOff.setColor(Color.getHSBColor(0.33f,0.7f,0.9f));
				}
				gOff.drawLine(xpts[i],ypts[i],xpts[i+1],ypts[i+1]);
			}
		}

		for(i=0; i<nnx; i += 1) {
			for (j=0; j<nny; j++) {
				if (mac.getCellKind(i,j)==0) {
					f = weight*mac.getPressure(i,j);
				} else {
					f = 0.0;
				}
				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)(scale*px)+gbx;
				ypts[j] = (int)(scale*py)+gby;
				fpts[j] = f; 
			}
			for (j=0; j<nny-1; j++) { 
				if (fpts[j]+fpts[j+1]>0.1) {
					gOff.setColor(Color.getHSBColor(0.10f,0.7f,0.9f));
				} else if (fpts[j]+fpts[j+1]<-0.1) {
					gOff.setColor(Color.getHSBColor(0.66f,0.7f,0.9f));
				} else {
					gOff.setColor(Color.getHSBColor(0.33f,0.7f,0.9f));
				}
				gOff.drawLine(xpts[j],ypts[j],xpts[j+1],ypts[j+1]);
			}
		}
	}

	// grid flow plot

	private void gridFlowPlot(double pweight, double fweight) {
		int i,j,k,gbx,gby;
		int cx,cy,cz;
		double z,f,px,py,pz,px2,py2,pz2,ivx,ivy,col,p0;
		double vmx,vmy;

		cx = nnx/2;
		cy = nny/2;
		cz = 0;

		gbx = 20+(int)(cx/scale);
		gby = 60+(int)(cy/scale);

		for(j=0; j<nny; j += 1) {
			for (i=0; i<nnx; i++) {
				if (mac.getCellKind(i,j)==0) {
					f = pweight*mac.getPressure(i,j);
				} else {
					f = 0.0;
				}
				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;

				vmx = mac.getvx(i,j); vmy = mac.getvy(i,j);
				ivx = i-cx+fweight*vmx; ivy = j-cy+fweight*vmy;
				py2 = cosFi*(ivy)+sinFi*(z-cz) + cy; 
				pz2 = -sinFi*(ivy)+cosFi*(z-cz) + cz;
				px2 = cosTh*(ivx)+sinTh*(pz2-cz) + cx; 
				pz2 = -sinTh*(ivx)+cosTh*(pz2-cz) + cz;

				if (f>=0.0) {
					if (vmx>0.0) {
						gOff.setColor(Color.getHSBColor(0.45f,0.7f,0.9f));
					} else {
						gOff.setColor(Color.getHSBColor(0.01f,0.7f,0.9f));
					}
				} else {
					if (vmx>0.0) {
						gOff.setColor(Color.getHSBColor(0.66f,0.7f,0.9f));
					} else {
						gOff.setColor(Color.getHSBColor(0.85f,0.7f,0.9f));
					}
				}
				gOff.drawLine( (int)(scale*px)+gbx, (int)(scale*py)+gby, (int)(scale*px2)+gbx, (int)(scale*py2)+gby );
			}
		}
	}


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

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

	public void setViewDelta(double xDisplacement, double yDisplacement) {
		viewTheta += 0.5*3.14159/180.0*xDisplacement;
		viewFai += 0.5*3.14159/180.0*yDisplacement;
		cosTh = Math.cos(viewTheta); 
		sinTh = Math.sin(viewTheta);
		cosFi = Math.cos(viewFai);
		sinFi = Math.sin(viewFai);
	}

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

	public void setDispMode(int m) { dispMode = m; }

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




// ======================================== MAC 2D class ===========

class MarkerAndCell2D implements Runnable {

	private double t = 0.0;
	private double dt = 0.0025;
	private double dx = 1.0/20.0;
	private double dy = 1.0/20.0;

	private double ReynoldsNumber = 400.0; 
	private double alpha = 3.0;
	private double fanFlow = 1.0;
	private double PressureReference = 1.0;
	private double pressureMean = 1.0;
	private double pressureError = 0.0;
	private double pressureErrorLimit = 0.001;  
	private double SORomega = 1.5;
	private int SORmaxLoop = 200;

	private int NNx = 140;
	private int NNy = 84;
	private byte cellKind[][] = new byte[NNx][NNy];
	private double vx[][] = new double[NNx][NNy];
	private double vy[][] = new double[NNx][NNy];
	private double nvx[][] = new double[NNx][NNy];
	private double nvy[][] = new double[NNx][NNy];
	private double source[][] = new double[NNx][NNy];
	private double pressure[][] = new double[NNx][NNy];
	private double readPressure[][] = new double[NNx][NNy];

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

	MarkerAndCell2D() {
		setInitialCondition();
	}


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

	public void run() {
		while (true) {
			timeEvolution();
			count += 1;

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


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

	private void setInitialCondition() {
		int i,j;
	
		t = 0.0;
		count = 0;

		setBoundary();

		for (i=0; i<NNx; i++) {  
			for (j=0; j<NNy; j++) {
				if (cellKind[i][j]==1) continue;
				pressure[i][j] = PressureReference;
				vx[i][j] = fanFlow;
				vy[i][j] = 0.0;
			}
		}
	}

	private void setBoundary() {
		int i,j;

		for (i=0; i<NNx; i++) {  
			for (j=0; j<NNy; j++) {
				cellKind[i][j] = 0;
				if (j<2 || j>=NNy-2) { // wall of wind tunnel
					cellKind[i][j] = 1;
					pressure[i][j] = PressureReference;
					vx[i][j] = 0.0;
					vy[i][j] = 0.0;
				} else if ( i>=40 && i<=43 && j>=2*NNy/3 ) { // wall
					cellKind[i][j] = 1;
					pressure[i][j] = PressureReference;
					vx[i][j] = 0.0;
					vy[i][j] = 0.0;
				}
			}
		}
	}


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

	private void timeEvolution() {

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

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

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

	private void timeStep() {
		t = t + dt;
		vEvolution();
		setFan(fanFlow);
		setNextPressure();
	}

	private void vEvolution() {
		int i,ip,ipp,im,imm,j,jp,jpp,jm,jmm;
		double absui,absvi, dpx,d2ux,d2uy,udux,vduy, dpy,d2vx,d2vy,udvx,vdvy;

		for (i=0; i<NNx; i++) {
			ip = i+1; if (ip>=NNx) ip = 0;
			ipp = ip+1; if (ipp>=NNx) ipp = 0;
			im = i-1; if (im<0) im = NNx-1;
			imm = im-1; if (imm<0) imm = NNx-1;	  
			for (j=0; j<NNy; j++) {
				if (cellKind[i][j]==1) continue;
				jp = j+1; if (jp>=NNy) jp = 0;
				jpp = jp+1; if (jpp>=NNy) jpp = 0;
				jm = j-1; if (jm<0) jm = NNy-1;
				jmm = jm-1; if (jmm<0) jmm = NNy-1;

				absui = vx[i][j]; if (absui<0.0) absui = -absui;
				absvi = vy[i][j]; if (absvi<0.0) absvi = -absvi;

				dpx = (pressure[ip][j]-pressure[im][j])/(2.0*dx);
				d2ux = (vx[ip][j]+vx[im][j]-2.0*vx[i][j])/(dx*dx);
				d2uy = (vx[i][jp]+vx[i][jm]-2.0*vx[i][j])/(dy*dy);

				udux = vx[i][j]*(-vx[ipp][j]+8.0*vx[ip][j]-8.0*vx[im][j]+vx[imm][j])/(12.0*dx) +
					alpha*absui*(vx[ipp][j]-4.0*vx[ip][j]+6*vx[i][j]-4.0*vx[im][j]+vx[imm][j])/(12.0*dx);
				vduy = vy[i][j]*(-vx[i][jpp]+8.0*vx[i][jp]-8.0*vx[i][jm]+vx[i][jmm])/(12.0*dy) +
					alpha*absvi*(vx[i][jpp]-4.0*vx[i][jp]+6*vx[i][j]-4.0*vx[i][jm]+vx[i][jmm])/(12.0*dy);
				nvx[i][j] = vx[i][j]+dt*(-udux-vduy-dpx+(d2ux+d2uy)/ReynoldsNumber);

				dpy = (pressure[i][jp]-pressure[i][jm])/(2.0*dy);
				d2vx = (vy[ip][j]+vy[im][j]-2.0*vy[i][j])/(dx*dx);
				d2vy = (vy[i][jp]+vy[i][jm]-2.0*vy[i][j])/(dy*dy);

				udvx = vx[i][j]*(-vy[ipp][j]+8.0*vy[ip][j]-8.0*vy[im][j]+vy[imm][j])/(12.0*dx) +
					alpha*absui*(vy[ipp][j]-4.0*vy[ip][j]+6*vy[i][j]-4.0*vy[im][j]+vy[imm][j])/(12.0*dx);
				vdvy = vy[i][j]*(-vy[i][jpp]+8.0*vy[i][jp]-8.0*vy[i][jm]+vy[i][jmm])/(12.0*dy) +
					alpha*absvi*(vy[i][jpp]-4.0*vy[i][jp]+6*vy[i][j]-4.0*vy[i][jm]+vy[i][jmm])/(12.0*dy);
				nvy[i][j] = vy[i][j]+dt*(-udvx-vdvy-dpy+(d2vx+d2vy)/ReynoldsNumber);
			}
		}

		for (i=0; i<NNx; i++) {  
			for (j=0; j<NNy; j++) {
				if (cellKind[i][j]==1) continue;
				vx[i][j] = nvx[i][j];
				vy[i][j] = nvy[i][j];
			}
		}
	}

	private void setFan(double v0) {
		int j;
	
		for (j=2; j<NNy-2; j++) {
			vx[0][j] = v0;
			vy[0][j] = 0.0+0.01*(Math.random()-0.5);
			vx[1][j] = v0;
			vy[1][j] = 0.0+0.01*(Math.random()-0.5);
		}
	}

	private void setNextPressure() {
		int i,j,ir;

		setSource();
	
		for (ir=0; ir<SORmaxLoop; ir++) {
			pressureError = poisson(SORomega);
			if (pressureError<pressureErrorLimit) break;
		}

		pressureMean = meanPressure();
		for (i=0; i<NNx; i++) {  
			for (j=0; j<NNy; j++) {
				pressure[i][j] -= pressureMean;
				readPressure[i][j] = pressure[i][j];
			}
		}
	}

	private double meanPressure() {
		int i,j,n;
		double s;

		s = 0.0; n = 0;
		for (i=0; i<NNx; i++) {
			for (j=0; j<NNy; j++) {
				if (cellKind[i][j]==0) {
					s += pressure[i][j];
					n += 1;
				}
			}
		}
		return(s/n);
	}

	private void setSource() {
		int i,ip,im,j,jp,jm;
		double dux,dvx,duy,dvy;

		for (i=0; i<NNx; i++) {
			ip = i+1; if (ip>=NNx) ip = 0;
			im = i-1; if (im<0) im = NNx-1;	  
			for (j=0; j<NNy; j++) {
				if (cellKind[i][j]==1) continue;
				jp = j+1; if (jp>=NNy) jp = 0;
				jm = j-1; if (jm<0) jm = NNy-1;

				dux = (vx[ip][j]-vx[im][j])/(2.0*dx);
				dvx = (vy[ip][j]-vy[im][j])/(2.0*dx);
				duy = (vx[i][jp]-vx[i][jm])/(2.0*dy);
				dvy = (vy[i][jp]-vy[i][jm])/(2.0*dy);
				source[i][j] = -(dux+dvy)/dt +(dux*dux+2.0*dvx*duy+dvy*dvy);
			}
		}
	}

	private double poisson(double omega) {
		int i,ip,im,j,jp,jm,ir,np;
		double dp,absdp,dpmax,p;

		dpmax = 0.0;
		for (i=0; i<NNx; i++) {
			ip = i+1; if (ip>=NNx) ip = 0;
			im = i-1; if (im<0) im = NNx-1;	  
			for (j=0; j<NNy; j++) {
				jp = j+1; if (jp>=NNy) jp = 0;
				jm = j-1; if (jm<0) jm = NNy-1;
				if (cellKind[i][j]==0) {
					dp = 0.25*(pressure[ip][j]+pressure[im][j]+pressure[i][jp]+pressure[i][jm]-4.0*pressure[i][j]+(dx*dx)*source[i][j]);
					pressure[i][j] += omega*dp;
					absdp = dp; if (absdp<0) absdp = -absdp;
					if (absdp>dpmax) dpmax = absdp;
				} else {
					p = 0.0; np = 0;
					if (cellKind[ip][j]==0) { p += pressure[ip][j]; np++; };
					if (cellKind[im][j]==0) { p += pressure[im][j]; np++; };
					if (cellKind[i][jp]==0) { p += pressure[i][jp]; np++; };
					if (cellKind[i][jm]==0) { p += pressure[i][jm]; np++; };
					if (np>0) {
						pressure[i][j] = p/np;
					} else {
						pressure[i][j] = PressureReference;
					}
				}
			}
		}
		return ( dpmax );
	}


	// ------------------------------------- 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 setReynoldsNumber(double re) { ReynoldsNumber = re; }
	public double getReynoldsNumber() { return ReynoldsNumber; }

	public void setFanFlow(double vx0) { fanFlow = vx0; }
	public double getFanFlow() { return fanFlow; }

	// get data
	public double getTime() { return t; }
	public double getdt() { return 0.01; }
	public double getdx() { return dx; }
	public double getdy() { return dy; }
	public int getNNx() { return NNx; }
	public int getNNy() { return NNy; }

	public double getvx(int i, int j) { return vx[i][j]; }
	public double getvy(int i, int j) { return vy[i][j]; }
	public double getPressure(int i, int j) { return readPressure[i][j]; }
	public int getCellKind(int i, int j) { return ((int)cellKind[i][j]); }
	public double getMeanPressure() { return pressureMean; }
	public double getPressureError() { return pressureError; } 
}


/*
 * MAC(marker and cell) method 2D (regular grid)
 * 
 * imcompressible fluid
 * equation of continuity
 *   div(v)=0   (I)
 * conservation of momentum (Navier-Stokes)
 *   dv/dt + (v.nabra)v = -grad(p)+(1/Re)div(grad(v))   (II)
 *
 * finite difference
 * dv/dt --> {v(t+dt) - v(t)}/dt
 * dvx/dx --> {vx[i+1][j] - vx[i-1][j]}/(2dx),
 * dvx/dy --> {vx[i][j+1] - vx[i][j-1]}/(2dy),
 * ...
 * dvx2/dx2 --> {vx[i+1][j] + vx[i-1][j] - 2vx[i][j]}/(dx^2),
 *
 * Kawamura-Kuwahara scheme alpha = 3,
 *   u df/dx --> u[i][j]*dr + alpha*|u[i][j]|*ds,
 *     dr = {-f[i+2][j]+8f[i+1][j]-8f[i-1][j]+f[i-2][j]}/(12dx),
 *     ds = {f[i+2][j]-4f[i+1][j]+6f[i][j]-4f[i-1][j]+f[i-2][j]}/(12dx),
 *
 * pressure <-- (I),(II)
 * div(grad(p)) = dD/dt -{(dvx/dx)^2+2(dvx/dy)(dvy/dx)+(dvy/dy)^2}
 *   D = div(v) = dvx/dx+dvy/dy
 *
 *
 * time evolution (vx,vy,p)
 *
 * (1) vx,vy,p --> next vx,vy
 *    v(t+dt) = v + dt(-v(v.nabra)v -grad(p)+(1/Re)div(grad(v)))
 *
 * (2) next vx,vy --> p(next)
 *    div(grad(p)) = -s,
 *      s = -dD/dt +{(dvx/dx)^2+2(dvx/dy)(dvy/dx)+(dvy/dy)^2},
 *      D = div(v).
 *    solve Poissons equation 
 *
 * goto (1)
 *
 */