//#extension GL_EXT_gpu_shader4 : require
#extension GL_EXT_gpu_shader4 : enable
//#extension GL_EXT_GPU_SHADER4 : require
//#version 140 // textureSize()&texelFetch() needs at least version 130, gl_InstanceID needs at least version 140

/****************************************************************************************************

	Code is based on

	.) Simple and Fast Fluid Flow Simulation on the GPU by Martin Guay, Fabrice Colin and Richard Egli
	available here:
	http://artis.imag.fr/Membres/Martin.Guay/

	.) http://http.developer.nvidia.com/GPUGems/gpugems_ch38.html

	.) http://http.developer.nvidia.com/GPUGems3/gpugems3_ch30.html	

****************************************************************************************************/




//#define OBSTACLES


varying vec2 v_tc0;

uniform sampler2D u_field_tex;
uniform sampler2D u_smokeDensity_tex;
uniform sampler2D u_obstacle_tex;

uniform sampler2D u_phi_n_hat_tex;
uniform sampler2D u_phi_n1_hat_tex;

uniform sampler2D u_vorticity_tex;

uniform sampler2D u_initialValue_tex;
uniform sampler2D u_base_tex;

uniform int u_callFunction;

uniform float u_velocityToQuantityTcScale;
uniform float u_vortConfinementScale;
uniform float u_dt;

uniform vec3 u_posAndRadius;
uniform float u_oldValScale;
uniform vec4 u_val;
uniform int u_forceMode;

uniform vec2 u_texSize;
uniform vec2 u_fieldTexStep; // 1.0 / textureSize( u_field_tex, 0 )
uniform vec2 u_smokeTexStep; // 1.0 / textureSize( u_smokeDensity_tex, 0 )

#define K  					0.01
#define viscosityFactor		0.1
#define gravityY			0.0

#define stabilityCondition	vec2( 3.0, 7.0 )


// helper functions to check for borders and obstacles
bool isBorder( vec2 tc, vec2 fieldStep )
{
	return ( tc.x <= fieldStep.x || tc.x > ( 1.0 - fieldStep.x ) || tc.y <= fieldStep.y || tc.y > ( 1.0 - fieldStep.y ) );
}

bool isSolidCell( vec2 cellTexCoords )
{
	return texture2D( u_obstacle_tex, cellTexCoords ).r < 0.8; // the smaller, the 'more  obstacle'
}

/*
bool isSolidCellTexOffset( vec2 cellTexCoords, ivec2 texOffset )
{
	// textureOffset requires #version 130 or later
	return texture2DOffset( u_obstacle_tex, cellTexCoords, texOffset ).r < 0.8; // the smaller, the 'more  obstacle'

	//http://developer.download.nvidia.com/opengl/specs/GL_EXT_gpu_shader4.txt
	//vec4 texture2DOffset(sampler2D sampler, vec2 coord, ivec2 offset [,float bias])
}
*/
bool isSolidCell( vec2 cellTexCoords, vec2 texOffset )
{
	return texture2D( u_obstacle_tex, cellTexCoords + texOffset ).r < 0.8; // the smaller, the 'more  obstacle'
}


void initDensityDomain()
{
	//gl_FragColor = vec4( texture2D( u_initialValue_tex, v_tc0 ).xxx, 1.0 ); // CORRECT
	gl_FragColor = texture2D( u_initialValue_tex, v_tc0 );
	//	gl_FragColor = vec4( 0.7, 0.2, 0.6, 1.0 );
}


float gaussian( vec2 pos, float radius )
{
  return exp( -dot(pos, pos) / radius );
}

void addForce()
{
	vec2 coords = v_tc0;

	float gauss = gaussian( ( coords - u_posAndRadius.xy ) * u_texSize, u_posAndRadius.z );

	if( u_forceMode == 0 ) {
		//gl_FragColor = u_oldValScale * texelFetch2D( u_base_tex, ivec2( coords * u_texSize ), 0 ) + u_val * gauss;
		gl_FragColor = u_oldValScale * texture2D( u_base_tex, floor( coords * u_texSize ) + 0.5 * ( 1.0 / u_texSize ), 0 ) + u_val * gauss;
	} else {
		// u_base_tex is now the field texture
		vec4 fieldC	= texture2D( u_base_tex, v_tc0 );
		/*
		vec3 fieldR	= texture2DOffset( u_base_tex, v_tc0, ivec2(+1, 0 ) ).xyz;
		vec3 fieldL	= texture2DOffset( u_base_tex, v_tc0, ivec2(-1, 0 ) ).xyz;
		vec3 fieldT	= texture2DOffset( u_base_tex, v_tc0, ivec2( 0,+1 ) ).xyz;
		vec3 fieldB	= texture2DOffset( u_base_tex, v_tc0, ivec2( 0,-1 ) ).xyz;
		*/
		vec2 stepSize = ( 1.0 / u_texSize );
		vec3 fieldR	= texture2D( u_base_tex, v_tc0 + vec2(+stepSize.x, 0 ) ).xyz;
		vec3 fieldL	= texture2D( u_base_tex, v_tc0 + vec2(-stepSize.x, 0 ) ).xyz;
		vec3 fieldT	= texture2D( u_base_tex, v_tc0 + vec2( 0,+stepSize.y ) ).xyz;
		vec3 fieldB	= texture2D( u_base_tex, v_tc0 + vec2( 0,-stepSize.y ) ).xyz;
		
		vec3 UdX = 0.5 * vec3( fieldR - fieldL );
		vec3 UdY = 0.5 * vec3( fieldT - fieldB );			
		vec2 DdX = vec2( UdX.z, UdY.z );

		float Udiv = UdX.x + UdY.y;
		fieldC.z -= u_dt * dot( vec3( DdX.xy, Udiv ), fieldC.xyz );

		// stability condition.
		//fieldC.z = clamp( fieldC.z, 0.1, 6.0 );
		fieldC.z = clamp( fieldC.z, stabilityCondition.x, stabilityCondition.y );

		gl_FragColor = u_oldValScale * fieldC + u_val * gauss;
	}
}


void singlePassSolver()
{	
	ivec2 directions[4];
	directions[0] = ivec2( 1, 0 );
	directions[1] = ivec2( 0, 1 );
	directions[2] = ivec2(-1, 0 );
	directions[3] = ivec2( 0,-1 );

	//vec2 texSize = vec2 ( textureSize( u_field_tex, 0 ) );
	//vec2 fieldStep = 1.0 / texSize;
	vec4 fieldC	= texture2D( u_field_tex, v_tc0 );

	bool isBoundary = texture2D( u_obstacle_tex, v_tc0 ).r < 0.8f;

	if ( isBoundary ) {
		fieldC.xy = vec2( 0.0 );
	}
		
	//if ( !isBoundary && !isBorder( v_tc0, fieldStep ) )	{
	if ( !isBoundary && !isBorder( v_tc0, u_fieldTexStep ) )	{
/*
		vec3 fieldR	= texture2DOffset( u_field_tex, v_tc0, ivec2( +1, 0 ) ).xyz;
		vec3 fieldL	= texture2DOffset( u_field_tex, v_tc0, ivec2( -1, 0 ) ).xyz;
		vec3 fieldT	= texture2DOffset( u_field_tex, v_tc0, ivec2(  0,+1 ) ).xyz;
		vec3 fieldB	= texture2DOffset( u_field_tex, v_tc0, ivec2(  0,-1 ) ).xyz;
*/
		vec3 fieldR	= texture2D( u_field_tex, v_tc0 + vec2( +u_fieldTexStep.x, 0 ) ).xyz;
		vec3 fieldL	= texture2D( u_field_tex, v_tc0 + vec2( -u_fieldTexStep.x, 0 ) ).xyz;
		vec3 fieldT	= texture2D( u_field_tex, v_tc0 + vec2(  0,+u_fieldTexStep.y ) ).xyz;
		vec3 fieldB	= texture2D( u_field_tex, v_tc0 + vec2(  0,-u_fieldTexStep.y ) ).xyz;

		// partial derivatives du/dx,du/dy
		vec3 UdX = 0.5 * ( fieldR - fieldL );
		vec3 UdY = 0.5 * ( fieldT - fieldB );
		
		vec2 DdX = vec2( UdX.z, UdY.z );
		
		// solve for velocity 	
		vec2 laplacian = fieldR.xy + fieldL.xy + fieldT.xy + fieldB.xy - 4.0 * fieldC.xy;
					
		vec2 PdX = ( K / u_dt ) * DdX;
		vec2 viscosityForce = viscosityFactor * laplacian;
		

		vec2 externalForces = vec2( 0.0 );
		if ( texture2D( u_smokeDensity_tex, v_tc0 ).x > 0.1f ) {
			externalForces += vec2( 0.0f, gravityY );			
		}
			
		// semi-lagrangian advection of velocity field itself
		//???? vec2 was = v_tc0 - u_dt * fieldC.xy * fieldStep;
		vec2 was = v_tc0 - u_dt * fieldC.xy * u_fieldTexStep;
		fieldC.xy = texture2D( u_field_tex, was ).xy;	
		fieldC.xy += u_dt * ( viscosityForce - PdX + externalForces );	


		// solve for density
		float Udiv = UdX.x + UdY.y;
		fieldC.z -= u_dt * dot( vec3( DdX.xy, Udiv ), fieldC.xyz );

		// stability condition.
		fieldC.z = clamp( fieldC.z, stabilityCondition.x, stabilityCondition.y );

		// boundary conditions.
		/*
		for ( int i = 0; i < 4; ++i )
		{
			if ( texture2DOffset( u_obstacle_tex, v_tc0, directions[i] ).x < 0.1 ) 
			{
				vec2 setToZero = ( 1.0 - abs( directions[i] ) );
				fieldC.xy *= setToZero;
			}
		}
		*/
	}	
		
	fieldC.w = 1.0;

	gl_FragColor = fieldC;
}


void updateDensity()
{	
	//vec2 densTexSize = vec2 ( textureSize( u_quantity_tex, 0 ) );
	//vec2 densFieldStep = 1.0 / densTexSize;

	float ink			= texture2D( u_smokeDensity_tex, v_tc0 ).x;
	float boundary		= texture2D( u_obstacle_tex, v_tc0 ).x;
	vec2 density = vec2( ink, boundary );

	if ( density.y > 0.8 )
	{					
		// semi-lagrangian advection of smoke particles
		//vec2 was 	= v_tc0 - u_dt * texture2D( u_field_tex, v_tc0 ).xy * densFieldStep;
		vec2 was 	= v_tc0 - u_dt * texture2D( u_field_tex, v_tc0 ).xy * u_smokeTexStep;
		density.x 	= texture2D( u_smokeDensity_tex, was ).x;	
	}

	gl_FragColor = vec4( density.xxx, 1.0 );
}


void updateDensityMacCormack(void)
{
	//vec2 densTexSize = vec2 ( textureSize( u_quantity_tex, 0 ) );
	//vec2 densFieldStep = 1.0 / densTexSize;


#ifdef OBSTACLES
	if ( isSolidCell( v_tc0 ) ) {
		out_color0 = vec2( 0.0 );
		return;
	}
#endif

	vec2 cellVelocity = texture2D( u_field_tex, v_tc0 ).xy;

	//vec2 backTrackedPos = v_tc0 - u_dt * cellVelocity * densFieldStep;
	vec2 backTrackedPos = v_tc0 - u_dt * cellVelocity * u_smokeTexStep;

	// limiter term
	float value[4];
	value[0] = texture2D( u_smokeDensity_tex, backTrackedPos + vec2( +0.5, +0.5 ) * u_smokeTexStep ).x;
	value[1] = texture2D( u_smokeDensity_tex, backTrackedPos + vec2( +0.5, -0.5 ) * u_smokeTexStep ).x;
	value[2] = texture2D( u_smokeDensity_tex, backTrackedPos + vec2( -0.5, +0.5 ) * u_smokeTexStep ).x;
	value[3] = texture2D( u_smokeDensity_tex, backTrackedPos + vec2( -0.5, -0.5 ) * u_smokeTexStep ).x;
	float phiMin = min( min( min( value[0], value[1] ), value[2] ), value[3] );
	float phiMax = max( max( max( value[0], value[1] ), value[2] ), value[3] );

	float phi_n1_hat = texture2D( u_phi_n1_hat_tex, backTrackedPos ).x;
	float phi_n = texture2D( u_smokeDensity_tex, v_tc0 ).x;
	float phi_n_hat = texture2D( u_phi_n_hat_tex, v_tc0 ).x;

	float result = phi_n1_hat + 0.5 * ( phi_n - phi_n_hat );

	gl_FragColor = vec4( vec3( clamp( result, phiMin, phiMax ) ), 1.0 );
}


void calculateVorticity()
{
	//vec2 densTexSize = vec2 ( textureSize( u_quantity_tex, 0 ) );
	//vec2 densFieldStep = 1.0 / densTexSize;
/*
	vec2 vR = texture2DOffset( u_field_tex, v_tc0, ivec2( +1,  0 ) ).xy;
	vec2 vL = texture2DOffset( u_field_tex, v_tc0, ivec2( -1,  0 ) ).xy;
	vec2 vT = texture2DOffset( u_field_tex, v_tc0, ivec2(  0, +1 ) ).xy;
	vec2 vB = texture2DOffset( u_field_tex, v_tc0, ivec2(  0, -1 ) ).xy;
*/
	vec2 vR = texture2D( u_field_tex, v_tc0 + vec2( +u_fieldTexStep.x,  0 ) ).xy;
	vec2 vL = texture2D( u_field_tex, v_tc0 + vec2( -u_fieldTexStep.x,  0 ) ).xy;
	vec2 vT = texture2D( u_field_tex, v_tc0 + vec2(  0, +u_fieldTexStep.y ) ).xy;
	vec2 vB = texture2D( u_field_tex, v_tc0 + vec2(  0, -u_fieldTexStep.y ) ).xy;

/*
	if ( isSolidCellTexOffset( v_tc0, ivec2( -1,  0 ) ) ) vL = vec2( 0.0 );
	if ( isSolidCellTexOffset( v_tc0, ivec2( +1,  0 ) ) ) vR = vec2( 0.0 );
	if ( isSolidCellTexOffset( v_tc0, ivec2(  0, -1 ) ) ) vB = vec2( 0.0 );
	if ( isSolidCellTexOffset( v_tc0, ivec2(  0, +1 ) ) ) vT = vec2( 0.0 );
*/

	if ( isSolidCell( v_tc0, vec2( -u_smokeTexStep.x,  0 ) ) ) vL = vec2( 0.0 );
	if ( isSolidCell( v_tc0, vec2( +u_smokeTexStep.x,  0 ) ) ) vR = vec2( 0.0 );
	if ( isSolidCell( v_tc0, vec2(  0, -u_smokeTexStep.y ) ) ) vB = vec2( 0.0 );
	if ( isSolidCell( v_tc0, vec2(  0, +u_smokeTexStep.y ) ) ) vT = vec2( 0.0 );

	float vorticity = 0.5 * ( ( vR.y - vL.y ) - ( vT.x - vB.x ) );

	//if ( isBorder( v_tc0, 1.0 / vec2( textureSize( u_field_tex, 0 ) ) ) )
	if ( isBorder( v_tc0, u_fieldTexStep ) )
		vorticity = 0.0;

	//if ( isSolidCell( v_tc0, ivec2(  0,  0 ) ) )
	if ( isSolidCell( v_tc0 ) )
		vorticity = 0.0;

	gl_FragColor = vec4( vec3( vorticity ), 1.0 );
	//gl_FragColor = vec4( vec3( vorticity*0.5 + 0.5 ), 1.0 ); // only for debugging, since vorticity may be negative => remap to positive range
}


void applyVorticity()
{
/*	
	float vortC = texture2D( u_vorticity_tex, v_tc0 ).r;
	float vortR = texture2DOffset( u_vorticity_tex, v_tc0, ivec2( +1,  0 ) ).r;
	float vortL = texture2DOffset( u_vorticity_tex, v_tc0, ivec2( -1,  0 ) ).r;
	float vortT = texture2DOffset( u_vorticity_tex, v_tc0, ivec2(  0, +1 ) ).r;
	float vortB = texture2DOffset( u_vorticity_tex, v_tc0, ivec2(  0, -1 ) ).r;
*/
	float vortC = texture2D( u_vorticity_tex, v_tc0 ).r;
	float vortR = texture2D( u_vorticity_tex, v_tc0 + vec2( +u_fieldTexStep.x,  0 ) ).r;
	float vortL = texture2D( u_vorticity_tex, v_tc0 + vec2( -u_fieldTexStep.x,  0 ) ).r;
	float vortT = texture2D( u_vorticity_tex, v_tc0 + vec2(  0, +u_fieldTexStep.y ) ).r;
	float vortB = texture2D( u_vorticity_tex, v_tc0 + vec2(  0, -u_fieldTexStep.y ) ).r;

/*
	if ( isSolidCellTexOffset( v_tc0, ivec2(  0,  0 ) ) ) vortC = 0.0;
	if ( isSolidCellTexOffset( v_tc0, ivec2( -1,  0 ) ) ) vortL = 0.0;
	if ( isSolidCellTexOffset( v_tc0, ivec2( +1,  0 ) ) ) vortR = 0.0;
	if ( isSolidCellTexOffset( v_tc0, ivec2(  0, -1 ) ) ) vortB = 0.0;
	if ( isSolidCellTexOffset( v_tc0, ivec2(  0, +1 ) ) ) vortT = 0.0;
*/


//!!!
	if ( isSolidCell( v_tc0 ) ) vortC = 0.0;
	if ( isSolidCell( v_tc0, vec2( -u_smokeTexStep.x,  0 ) ) ) vortL = 0.0;
	if ( isSolidCell( v_tc0, vec2( +u_smokeTexStep.x,  0 ) ) ) vortR = 0.0;
	if ( isSolidCell( v_tc0, vec2(  0, -u_smokeTexStep.y ) ) ) vortB = 0.0;
	if ( isSolidCell( v_tc0, vec2(  0, +u_smokeTexStep.y ) ) ) vortT = 0.0;

	vec2 force = 0.5 * vec2( abs( vortT ) - abs( vortB ), abs( vortR ) - abs( vortL ) );

	// safe normalize
	const float EPSILON = 2.4414e-4; // 2^-12
	float magSqr = max( EPSILON, dot( force, force ) ); 
	force = force * inversesqrt( magSqr ); 

	force *= u_vortConfinementScale * vortC * vec2( 1.0f, -1.0f );

	vec4 field = texture2D( u_field_tex, v_tc0 );
	vec2 velocity = field.xy + u_dt * force;

	//if( isBorder( v_tc0, 1.0 / vec2( textureSize( u_field_tex, 0 ) ) ) )
	//!! if( isBorder( v_tc0, u_fieldTexStep ) )
	//!! 	velocity = field.xy;

	gl_FragColor = vec4( velocity, field.zw );
}



//void main()
//{
//	gl_FragColor = vec4( 0.7, 0.2, 0.6, 1.0 );
//}

/*
void main()
{
	switch( u_callFunction ) {

		case 0:
			initDensityDomain();
			break;

		case 10:
			addForce();
			break;

		case 20:
			singlePassSolver();
			break;

		case 30:
			updateDensity();
			break;

		case 40:
			updateDensityMacCormack();
			break;

		case 50:
			calculateVorticity();
			break;

		case 60:
			applyVorticity();
			break;

		default:
			discard;
			break;
	}
}
*/

void main()
{
	if ( u_callFunction == 0 ) {
		initDensityDomain();
	} else if ( u_callFunction == 10 ) {
		addForce();
	} else if ( u_callFunction == 20 ) {
		singlePassSolver();
	} else if ( u_callFunction == 30 ) {
		updateDensity();
	} else if ( u_callFunction == 40 ) {
		updateDensityMacCormack();
	} else if ( u_callFunction == 50 ) {
		calculateVorticity();
	} else if ( u_callFunction == 60 ) {
		applyVorticity();
	} else {
		discard;
	}
}
