/*
 * rotor in the wind tunnel - MAC(Marker and Cell) method 2D
 *
 * coded by Ikeuchi Mitsuru
 * ver 0.0.1   2006.03.05  created
 *
 */

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

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

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

  Choice ch_view;
  Button bt_reset, bt_startStop;
  Scrollbar sc_Reynolds,sc_vx0,sc_vr0;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 50;

  int resetSW = 0;
  int started = 1;
  int dispMode = 5;

  double t = 0.0;
  double dt = 0.01;
  double dx = 1.0/20.0;
  double dy = 1.0/20.0;
  int NNx = 140;
  int NNy = 84;

  double ReynoldsNumber = 400.0; 
  double alpha = 3.0;
  double vx0 = 1.0;
  double vr0 = 1.0;
  double pressure0 = 1.0;
  double pError = 0.0;
  double dispScale = (400.0/(NNx*dx));
  double viewScale = 1.0;

  int dragSW = 2;
  int dgX, dgY, dgXb, dgYb;
  double xShift = 0.0;
  double yShift = 0.0;

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


  int cellKind[][] = new int[NNx][NNy];
  double vx[][] = new double[NNx][NNy];
  double vy[][] = new double[NNx][NNy];
  double nvx[][] = new double[NNx][NNy];
  double nvy[][] = new double[NNx][NNy];
  double source[][] = new double[NNx][NNy];
  double pressure[][] = new double[NNx][NNy];

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

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

  public void init() {
    resize(630,420);
    setBackground(Color.white);

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

    ch_view = new Choice();
    ch_view.addItem("flow");
    ch_view.addItem("pressure");
    ch_view.addItem("flow+P");
    ch_view.addItem("grid P");
    ch_view.addItem("grid v,P");
    ch_view.addItem("flow line");
    ch_view.addItemListener(this);
    ch_view.select("flow line");

    bt_reset= new Button("reset");
    bt_reset.addActionListener(this);

    bt_startStop= new Button("start/stop");
    bt_startStop.addActionListener(this); 
    
    sc_Reynolds = new Scrollbar(Scrollbar.HORIZONTAL,400,10,50,410);
    sc_Reynolds.addAdjustmentListener(this);

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

    sc_vr0 = new Scrollbar(Scrollbar.HORIZONTAL,100,10,0,110);
    sc_vr0.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(sc_Reynolds);
    pnl.add(sc_vx0);
    pnl.add(sc_vr0);
    pnl.add(ch_view);
    add(pnl,"North");

    setInitialCondition();
  }

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

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

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

  public void actionPerformed(ActionEvent ev){
    if(ev.getSource() == bt_reset){ 
      resetSW = 1;
    } else if (ev.getSource() == bt_startStop){
      if (started==0) {
        started = 1;
      } else {
        started = 0;
      }
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_Reynolds) { 
      ReynoldsNumber = 1.0*(double)(sc_Reynolds.getValue());
    } else if (ev.getSource() == sc_vx0) { 
      vx0 = 0.01*(double)(sc_vx0.getValue());
    } else if (ev.getSource() == sc_vr0) { 
      vr0 = 0.01*(double)(sc_vr0.getValue());
	  setRotor();
    }
  }

  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;
    }
	
    if (dragSW==1) {
      xShift += (1.0/viewScale/dispScale)*(dgX-dgXb);
      yShift += (1.0/viewScale/dispScale)*(dgY-dgYb);
    } else if (dragSW==2) {
      viewTheta += 0.5*3.14/180.0*(dgX-dgXb);
      viewFai += 0.5*3.14/180.0*(dgY-dgYb);
      cosTh = Math.cos(viewTheta); 
      sinTh = Math.sin(viewTheta);
      cosFi = Math.cos(viewFai);
      sinFi = Math.sin(viewFai);
    }
  }

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

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

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

/* ----------------------------- offPaint -------------------- */

  void offPaint() {
    int gx,gy;
	
	gx = 20;
	gy = 80;

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

    if (dispMode==0) {
      boundaryPlot(gx,gy);
      flowPlot(gx,gy, 20.0);

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

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

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

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

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

    } else if (dispMode==6) {
      ;

    }

    gOff.setColor(Color.blue);
    gOff.drawString("t="+(int)(t*100.0+0.5)/100.0+" ",10,40);
    if (started==0) {
      gOff.drawString("stopped",630/6*1+10,40);
	} else {
      gOff.drawString("started",630/6*1+10,40);
	}
    gOff.drawString("Re="+(int)(ReynoldsNumber+0.5)+" ",630/6*2+10,40);
    gOff.drawString("fan="+(int)(vx0*100.0+0.5)/100.0+" ",630/6*3+10,40);
    gOff.drawString("rotate="+(int)(vr0*100.0+0.5)/100.0+" ",630/6*4+10,40);
    gOff.drawString("view",630/6*5+10,40);

    gOff.setColor(Color.black);
    gOff.drawString("Box="+(int)(NNx*dx*10.0+0.5)/10.0+"x"+(int)(NNy*dy*10.0+0.5)/10.0+" ",10,60);
/*
    gOff.setColor(Color.blue);
    gOff.drawString("dpMax="+(float)(pError)+" ",150,60);
*/
  }

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

  void flowPlot(int gx, int gy, double vmag) {
    int i,j,ix,iy,ix2,iy2;
    double sc,viewScale;

    sc = 4.0; viewScale=1.0;
    for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        ix = (int)(sc*(i+0.5))+gx;
        iy = (int)(sc*(j+0.5))+gy;

        ix2 = ix+(int)(viewScale*vmag*vx[i][j]);
        iy2 = iy+(int)(viewScale*vmag*vy[i][j]);
	    if (vx[i][j]>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);
      }
    }
  }

  void pressurePlot(int gx, int gy, double mag) {
    int i,j,ix,iy;
    double sc,viewScale,c,absc,p0;

    p0 = meanPressure();
    sc = 4.0; viewScale=1.0;
    for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        ix = (int)(sc*i)+gx;
        iy = (int)(sc*j)+gy;
		c = mag*(pressure[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);
      }
    }
  }

  void boundaryPlot(int gx, int gy) {
    int i,j,ix,iy;
    double sc;

    sc = 4.0;
    for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
	    if (cellKind[i][j]==1) {
          ix = (int)(sc*i)+gx;
          iy = (int)(sc*j)+gy;
	      gOff.setColor(Color.gray);
          gOff.fillRect(ix,iy, 4,4);
		}
      }
    }
  }

  void flowLinePlot(int gx, int gy) {
    int i,j,ix,iy,ix2,iy2,px,py;
    double sc,viewScale,sx,sy;

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

  void flowLine(int gx, int gy, double x0, double y0, int tstep) {
    int i,j,ix,iy,ix2,iy2,px,py;
    double sc,viewScale,sx,sy;

    sc = 4.0; viewScale=1.0;
	sx = x0; sy = y0;
    for (i=0; i<tstep; i++) {
      ix = (int)(sc*sx/dx+0.5)+gx;
      iy = (int)(sc*sy/dy+0.5)+gy;
	  px = (int)(sx/dx+0.5);if (px>=NNx) break;
	  py = (int)(sy/dy+0.5);if (py>=NNy) break;
      sx += vx[px][py]*dt;
	  sy += vy[px][py]*dt;
      ix2 = (int)(sc*sx/dx+0.5)+gx;
      iy2 = (int)(sc*sy/dy+0.5)+gy;
	  gOff.setColor(Color.red);
	  gOff.drawLine(ix,iy, ix2,iy2);
    }
  }

/* -------------------------- grid pressure plot -------------- */

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

    cx = NNx/2;
    cy = NNy/2;
    cz = 0;
    sc = 4.0;

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

    for(j=0; j<NNy; j += 1) {
      for (i=0; i<NNx; i++) {
	    if (cellKind[i][j]==0) {
          f = weight*pressure[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)(sc*px)+gbx;
        ypts[i] = (int)(sc*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 (cellKind[i][j]==0) {
          f = weight*pressure[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)(sc*px)+gbx;
        ypts[j] = (int)(sc*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 -------------- */

  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,sc,col,p0;

    cx = NNx/2;
    cy = NNy/2;
    cz = 0;
    sc = 4.0;

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

    for(j=0; j<NNy; j += 1) {
      for (i=0; i<NNx; i++) {
	    if (cellKind[i][j]==0) {
          f = pweight*pressure[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;

        ivx = i-cx+fweight*vx[i][j]; ivy = j-cy+fweight*vy[i][j];
        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 (vx[i][j]>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 (vx[i][j]>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)(sc*px)+gbx,(int)(sc*py)+gby,(int)(sc*px2)+gbx,(int)(sc*py2)+gby);
      }
    }
  }

/*------------------------------------- statistics --------- */

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

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

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

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

  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) {
	      cellKind[i][j] = 1;
          pressure[i][j] = pressure0;
          vx[i][j] = 0.0;
          vy[i][j] = 0.0;
        } else if ( (i-50)*(i-50)+(j-NNy/2)*(j-NNy/2)<90 ) {
	      cellKind[i][j] = 1;
          pressure[i][j] = pressure0;
          vx[i][j] = -1.0*vr0*(j-NNy/2)/10.0;
          vy[i][j] = 1.0*vr0*(i-50)/10.0;
		}
      }
    }
  }

  void setRotor() {
    int i,j;

    for (i=0; i<NNx; i++) {  
      for (j=0; j<NNy; j++) {
	    if ( (i-50)*(i-50)+(j-NNy/2)*(j-NNy/2)<90 ) {
          vx[i][j] = -1.0*vr0*(j-NNy/2)/10.0;
          vy[i][j] = 1.0*vr0*(i-50)/10.0;
		}
      }
    }
  }

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

  void timeEvolution() {

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

    if (started==1) {
      t = t + dt;
      vEvolution();
	  setFan(vx0);
	  setNextPressure(100);
    }
  }

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

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

  void setNextPressure(int iter) {
    int i,j,ir;
    double p0;

    setSource();
	
    for (ir=0; ir<iter; ir++) {
      pError = poisson(1.5);
	}

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

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

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

    setSource();

    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] = pressure0;
          }
        }
      }
    }
    return ( dpmax );
  }


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

}

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