shallow water equation c++ code

// minimal shallow water solver by Nils Thuerey
// see http://www.ntoken.com and http://www.matthiasmueller.info/realtimephysics/ for further info
// this code is released under the GPL see LICENSE.txt for details
#include 
#include 
#include 
#include 

// global parameters
const int SIZE = 100;

const int IMAGESTEPS = 2;

const int END_TIME = 25;

// single / double precision?
typedef float Real;

// use the turbulence model?
const bool SMAGORINSKY = false;

const int FLUID  = 0 NOSLIP = 1 VELOCITY  = 2;
const Real LDC_VELOCITY = 0.10;

// constants
const Real gravity = -10.;  // normal acceleration a_n
const Real dt = 10. / （Real）SIZE;       // time step size
const Real domainSize = 50.; // size of the domain
const Real dx = domainSize / （Real）SIZE;  
const Real dxInv = 1./dx;   // save some division later on

// impose global velocity
const Real GLOBAL_U = 0.;
const Real GLOBAL_V = 0.;

// main arrays velocity height and temporary storage
std::vector vel_x vel_y temp eta;

void writeImage（int s）;
void addRandomDrop（）;

Real interpolate（std::vector &array Real x Real y） {
	const int X = （int）x;
	const int Y = （int）y;
	const Real s1 = x - X;
	const Real s0 = 1.f - s1;
	const Real t1 = y - Y;
	const Real t0 = 1.f-t1;
	return  s0*（t0* array[X+SIZE*Y] + t1*array[X  +SIZE*（Y+1）] ）+
		s1*（t0*array[（X+1）+SIZE*Y]  + t1*array[（X+1）+SIZE*（Y+1）] ）;
}

void advectArray（std::vector &array int velocityType） { 
	for （int i=1;i		for （int j=1;j			const int index = i + j*SIZE;

			// distinguish different cases to interpolate the velocity
			Real u = GLOBAL_U v = GLOBAL_V; 
			switch（velocityType） {
				case 0: // heights
					u += （vel_x[index]+vel_x[index+1]） *0.5;
					v += （vel_y[index]+vel_y[index+SIZE]） *0.5;
					break;
				case 1: // x velocity
					u += vel_x[index];
					v += （vel_y[index]+vel_y[index+1]+vel_y[index+SIZE]+vel_y[index+SIZE+1]） *0.25;
					break;
				case 2: // y velocity
					u += （vel_x[index]+vel_x[index+1]+vel_x[index+SIZE]+vel_x[index+SIZE+1]） *0.25;
					v += vel_y[index];
					break;
				default: // invalid
					exit（1）;
			}

			// backtrace position
			Real srcpi = （Real）i - u * dt * dxInv;
			Real srcpj = （Real）j - v * dt * dxInv;

			// clamp range of accesses
			if（srcpi<0.） srcpi = 0.;
			if（srcpj<0.） srcpj = 0.;
			if（srcpi>SIZE-1.） srcpi = SIZE-1.;
			if（srcpj>SIZE-1.） srcpj = SIZE-1.;

			// interpolate source value
			temp[index] = interpolate（array srcpi srcpj）;
		}

	// copy back...
	for （int i=1;i		for （int j=1;j			const int index = i + j*SIZE;
			array[index] = temp[index];
		}
}


void updateHeight（） {
	// update heights
	for （int i=1;i		for （int j=1;j			const int index = i + j*SIZE;
			Real dh = -0.5 * eta[index] * dxInv * （
					 （vel_x[index+1]    - vel_x[index]） +
					 （vel_y[index+SIZE] - vel_y[index]） ）;

			eta[index] += dh * dt;
		}
}

void updateVelocities（