//?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（