int IX(int i, int j) {
  return i+j*(N+2);
}

//////////////////////////////////////////////////////////////////////////////////////////

int[] XI(int k) {
  int i = k%(N+2);
  int j = int(k/float(N+2));
  int[] t = {
    i, j
  };
  return t;
}

//////////////////////////////////////////////////////////////////////////////////////////

void swap(float[] x0, float[] x) {
  float[] tmp = new float[x0.length];
  arrayCopy(x0, 0, tmp, 0, x0.length);
  arrayCopy(x, 0, x0, 0, x0.length);
  arrayCopy(tmp, 0, x, 0, x0.length);
}

//////////////////////////////////////////////////////////////////////////////////////////

void add_source(float[] x, float[] s, float dt) {
  for (int i=0; i<(N+2)*(N+2); i++) {
    x[i] += s[i]*dt;
  }
}


//////////////////////////////////////////////////////////////////////////////////////////

void add_force(float[] u, float[] fx, float dt) {
  for (int i=0; i<(N+2)*(N+2); i++) {
    u[i] += fx[i]*dt;
  }
}

//////////////////////////////////////////////////////////////////////////////////////////

void set_bnd(int b, float[] x)
{
  for (int i=1; i<N; i++) {
    for (int j=1; j<N; j++) {
      if (isFluid[IX(i,j)]==false) {
        vx[IX(i,j)]   = (vy[IX(i-1,j)]-vy[IX(i+1,j)])/h;
        vy[IX(i,j)]   = (vx[IX(i,j-1)]-vx[IX(i,j+1)])/h;       
        dens[IX(i,j)] = (0.5*(dens[IX(i,j-1)]-dens[IX(i,j+1)])+0.5*(dens[IX(i-1,j)]-dens[IX(i+1,j)]))/h;;
        s[IX(i,j)]    = 0;        
      }
    }
  }

  for (int i=1 ; i<=N ; i++) {
    if ( b==1 ) {
      x[IX(0, i)]   = -x[IX(1, i)]; 
      x[IX(N+1, i)] = -x[IX(N, i)];
    }    
    else {
      x[IX(N+1, i)] = x[IX(N, i)];
      x[IX(0, i)]   = x[IX(1, i)];
    }
    if ( b==2 ) {
      x[IX(i, 0)]   = -x[IX(i, 1)]; 
      x[IX(i, N+1)] = -x[IX(i, N)];
    }   
    else {
      x[IX(i, N+1)] = x[IX(i, N)];
      x[IX(i, 0)]   = x[IX(i, 1)];
    }
  }
  x[IX(0, 0 )]    = 0.5*(x[IX(1, 0 )]  + x[IX(0, 1)]);
  x[IX(0, N+1)]   = 0.5*(x[IX(1, N+1)] + x[IX(0, N )]);
  x[IX(N+1, 0 )]  = 0.5*(x[IX(N, 0 )]  + x[IX(N+1, 1)]);
  x[IX(N+1, N+1)] = 0.5*(x[IX(N, N+1)] + x[IX(N+1, N )]);
}

//////////////////////////////////////////////////////////////////////////////////////////

void diffuse(int b, float[] x, float[] x0, float diff, float dt)
{
  int i, j, k;
  float a=dt*diff*N*N;
  // (GAUSS-SEIDEL RELAXATION)
  for ( k=0 ; k<30 ; k++ ) {
    for ( i=1 ; i<=N ; i++ ) {
      for ( j=1 ; j<=N ; j++ ) {
        x[IX(i, j)] = (x0[IX(i, j)] + a*(x[IX(i-1, j)]+x[IX(i+1, j)]+
          x[IX(i, j-1)]+x[IX(i, j+1)]))/(1+4*a);
      }
    }
    set_bnd (b, x );
  }
}

//////////////////////////////////////////////////////////////////////////////////////////

void advect(int b, float[] d, float[] d0, float[] u, float[] v, float dt)
{
  int i0, j0, i1, j1;
  float x, y, s0, t0, s1, t1, dt0;
  dt0 = dt*N;
  for (int i=1 ; i<=N ; i++ ) {
    for (int j=1 ; j<=N ; j++ ) {
      x = i-dt0*u[IX(i, j)]; 
      y = j-dt0*v[IX(i, j)];
      if (x<0.5) x=0.5; 
      if (x>N+0.5) x=N+ 0.5; 
      i0=(int)x; 
      i1=i0+1;
      if (y<0.5) y=0.5; 
      if (y>N+0.5) y=N+ 0.5; 
      j0=(int)y; 
      j1=j0+1;
      s1 = x-i0; 
      s0 = 1-s1; 
      t1 = y-j0; 
      t0 = 1-t1;
      d[IX(i, j)] = s0*(t0*d0[IX(i0, j0)]+t1*d0[IX(i0, j1)])+
        s1*(t0*d0[IX(i1, j0)]+t1*d0[IX(i1, j1)]);
    }
  }
  set_bnd(b, d);
}


//////////////////////////////////////////////////////////////////////////////////////////

void project(float[] u, float[] v, float[] p, float[] div)
{
  int i, j, k;
  float h;
  h = w/N;
  for ( i=1 ; i<=N ; i++ ) {
    for ( j=1 ; j<=N ; j++ ) {
      div[IX(i, j)] = -0.5*h*(u[IX(i+1, j)]-u[IX(i-1, j)]+
        v[IX(i, j+1)]-v[IX(i, j-1)]);
      p[IX(i, j)] = 0;
    }
  }
  set_bnd(0, div); 
  set_bnd(0, p);
  for ( k=0 ; k<30 ; k++ ) {
    for ( i=1 ; i<=N ; i++ ) {
      for ( j=1 ; j<=N ; j++ ) {
        p[IX(i, j)] = (div[IX(i, j)]+p[IX(i-1, j)]+p[IX(i+1, j)]+
          p[IX(i, j-1)]+p[IX(i, j+1)])/4;
      }
    }
    set_bnd(0, p);
  }
  for ( i=1 ; i<=N ; i++ ) {
    for ( j=1 ; j<=N ; j++ ) {
      u[IX(i, j)] -= 0.5*(p[IX(i+1, j)]-p[IX(i-1, j)])/h;
      v[IX(i, j)] -= 0.5*(p[IX(i, j+1)]-p[IX(i, j-1)])/h;
    }
  }
  set_bnd(1, u); 
  set_bnd(2, v);
}

//////////////////////////////////////////////////////////////////////////////////////////

void calculateVelocity(float[] u, float[] v, float[] u0, float[] v0, float visc, float dt) {
  add_force(u, fx, dt);
  add_force(v, fy, dt); 
  add_source(u, u0, dt); 
  add_source(v, v0, dt);
  swap(u0, u); 
  swap(v0, v);
  diffuse(1, u, u0, visc, dt);   
  diffuse(2, v, v0, visc, dt);
  project(u, v, u0, v0);
  swap(u0, u); 
  swap(v0, v);
  advect(1, u, u0, u0, v0, dt); 
  advect(2, v, v0, u0, v0, dt);
  project(u, v, u0, v0);
}

//////////////////////////////////////////////////////////////////////////////////////////

void calculateDensity(float[] u, float[] v, float[] dens, float[] dens0, float[] s, float diff, float dt) {
  add_source( dens, s, dt ); 
  swap(dens0, dens);
  diffuse( 0, dens, dens0, diff, dt );  
  swap(dens0, dens);
  advect( 0, dens, dens0, u, v, dt );
}

//////////////////////////////////////////////////////////////////////////////////////////