#version 120

//uniform sampler2D nv_t_diffuse;
uniform vec4 light_diffuse;
uniform vec4 light_specular;

//varying vec2 v_texcoord;
//varying vec3 v_normal;
//varying vec3 v_light_vector;
//varying vec3 v_view_vector;
varying vec3 v_world_pos;
varying vec3 v_camera_loc;
varying vec3 v_light_loc;
uniform vec3 center;
uniform float radius;
uniform mat4 nv_m_mvp;

vec3 mod289(vec3 x) {
  return x - floor(x * (1.0 / 289.0)) * 289.0;
}

vec4 mod289(vec4 x) {
  return x - floor(x * (1.0 / 289.0)) * 289.0;
}

vec4 permute(vec4 x) {
     return mod289(((x*34.0)+1.0)*x);
}

vec4 taylorInvSqrt(vec4 r)
{
  return 1.79284291400159 - 0.85373472095314 * r;
}

#define jitter 1.0 // smaller jitter gives more regular pattern

float snoise(vec3 v)
  { 
  const vec2  C = vec2(1.0/6.0, 1.0/3.0) ;
  const vec4  D = vec4(0.0, 0.5, 1.0, 2.0);

// First corner
  vec3 i  = floor(v + dot(v, C.yyy) );
  vec3 x0 =   v - i + dot(i, C.xxx) ;

// Other corners
  vec3 g = step(x0.yzx, x0.xyz);
  vec3 l = 1.0 - g;
  vec3 i1 = min( g.xyz, l.zxy );
  vec3 i2 = max( g.xyz, l.zxy );

  //   x0 = x0 - 0.0 + 0.0 * C.xxx;
  //   x1 = x0 - i1  + 1.0 * C.xxx;
  //   x2 = x0 - i2  + 2.0 * C.xxx;
  //   x3 = x0 - 1.0 + 3.0 * C.xxx;
  vec3 x1 = x0 - i1 + C.xxx;
  vec3 x2 = x0 - i2 + C.yyy; // 2.0*C.x = 1/3 = C.y
  vec3 x3 = x0 - D.yyy;      // -1.0+3.0*C.x = -0.5 = -D.y

// Permutations
  i = mod289(i); 
  vec4 p = permute( permute( permute( 
             i.z + vec4(0.0, i1.z, i2.z, 1.0 ))
           + i.y + vec4(0.0, i1.y, i2.y, 1.0 )) 
           + i.x + vec4(0.0, i1.x, i2.x, 1.0 ));

// Gradients: 7x7 points over a square, mapped onto an octahedron.
// The ring size 17*17 = 289 is close to a multiple of 49 (49*6 = 294)
  float n_ = 0.142857142857; // 1.0/7.0
  vec3  ns = n_ * D.wyz - D.xzx;

  vec4 j = p - 49.0 * floor(p * ns.z * ns.z);  //  mod(p,7*7)

  vec4 x_ = floor(j * ns.z);
  vec4 y_ = floor(j - 7.0 * x_ );    // mod(j,N)

  vec4 x = x_ *ns.x + ns.yyyy;
  vec4 y = y_ *ns.x + ns.yyyy;
  vec4 h = 1.0 - abs(x) - abs(y);

  vec4 b0 = vec4( x.xy, y.xy );
  vec4 b1 = vec4( x.zw, y.zw );

  //vec4 s0 = vec4(lessThan(b0,0.0))*2.0 - 1.0;
  //vec4 s1 = vec4(lessThan(b1,0.0))*2.0 - 1.0;
  vec4 s0 = floor(b0)*2.0 + 1.0;
  vec4 s1 = floor(b1)*2.0 + 1.0;
  vec4 sh = -step(h, vec4(0.0));

  vec4 a0 = b0.xzyw + s0.xzyw*sh.xxyy ;
  vec4 a1 = b1.xzyw + s1.xzyw*sh.zzww ;

  vec3 p0 = vec3(a0.xy,h.x);
  vec3 p1 = vec3(a0.zw,h.y);
  vec3 p2 = vec3(a1.xy,h.z);
  vec3 p3 = vec3(a1.zw,h.w);

//Normalise gradients
  vec4 norm = taylorInvSqrt(vec4(dot(p0,p0), dot(p1,p1), dot(p2, p2), dot(p3,p3)));
  p0 *= norm.x;
  p1 *= norm.y;
  p2 *= norm.z;
  p3 *= norm.w;

// Mix final noise value
  vec4 m = max(0.6 - vec4(dot(x0,x0), dot(x1,x1), dot(x2,x2), dot(x3,x3)), 0.0);
  m = m * m;
  return 42.0 * dot( m*m, vec4( dot(p0,x0), dot(p1,x1), 
                                dot(p2,x2), dot(p3,x3) ) );
}

// Permutation polynomial: (34x^2 + x) mod 289
vec3 permute(vec3 x) {
  return mod((34.0 * x + 1.0) * x, 289.0);
}

#define jitter2x2x2 0.8 // smaller jitter gives less errors in F2

vec2 cellular(vec3 P) {
#define K 0.142857142857 // 1/7
#define Ko 0.428571428571 // 1/2-K/2
#define K2 0.020408163265306 // 1/(7*7)
#define Kz 0.166666666667 // 1/6
#define Kzo 0.416666666667 // 1/2-1/6*2

	vec3 Pi = mod(floor(P), 289.0);
 	vec3 Pf = fract(P) - 0.5;

	vec3 Pfx = Pf.x + vec3(1.0, 0.0, -1.0);
	vec3 Pfy = Pf.y + vec3(1.0, 0.0, -1.0);
	vec3 Pfz = Pf.z + vec3(1.0, 0.0, -1.0);

	vec3 p = permute(Pi.x + vec3(-1.0, 0.0, 1.0));
	vec3 p1 = permute(p + Pi.y - 1.0);
	vec3 p2 = permute(p + Pi.y);
	vec3 p3 = permute(p + Pi.y + 1.0);

	vec3 p11 = permute(p1 + Pi.z - 1.0);
	vec3 p12 = permute(p1 + Pi.z);
	vec3 p13 = permute(p1 + Pi.z + 1.0);

	vec3 p21 = permute(p2 + Pi.z - 1.0);
	vec3 p22 = permute(p2 + Pi.z);
	vec3 p23 = permute(p2 + Pi.z + 1.0);

	vec3 p31 = permute(p3 + Pi.z - 1.0);
	vec3 p32 = permute(p3 + Pi.z);
	vec3 p33 = permute(p3 + Pi.z + 1.0);

	vec3 ox11 = fract(p11*K) - Ko;
	vec3 oy11 = mod(floor(p11*K), 7.0)*K - Ko;
	vec3 oz11 = floor(p11*K2)*Kz - Kzo; // p11 < 289 guaranteed

	vec3 ox12 = fract(p12*K) - Ko;
	vec3 oy12 = mod(floor(p12*K), 7.0)*K - Ko;
	vec3 oz12 = floor(p12*K2)*Kz - Kzo;

	vec3 ox13 = fract(p13*K) - Ko;
	vec3 oy13 = mod(floor(p13*K), 7.0)*K - Ko;
	vec3 oz13 = floor(p13*K2)*Kz - Kzo;

	vec3 ox21 = fract(p21*K) - Ko;
	vec3 oy21 = mod(floor(p21*K), 7.0)*K - Ko;
	vec3 oz21 = floor(p21*K2)*Kz - Kzo;

	vec3 ox22 = fract(p22*K) - Ko;
	vec3 oy22 = mod(floor(p22*K), 7.0)*K - Ko;
	vec3 oz22 = floor(p22*K2)*Kz - Kzo;

	vec3 ox23 = fract(p23*K) - Ko;
	vec3 oy23 = mod(floor(p23*K), 7.0)*K - Ko;
	vec3 oz23 = floor(p23*K2)*Kz - Kzo;

	vec3 ox31 = fract(p31*K) - Ko;
	vec3 oy31 = mod(floor(p31*K), 7.0)*K - Ko;
	vec3 oz31 = floor(p31*K2)*Kz - Kzo;

	vec3 ox32 = fract(p32*K) - Ko;
	vec3 oy32 = mod(floor(p32*K), 7.0)*K - Ko;
	vec3 oz32 = floor(p32*K2)*Kz - Kzo;

	vec3 ox33 = fract(p33*K) - Ko;
	vec3 oy33 = mod(floor(p33*K), 7.0)*K - Ko;
	vec3 oz33 = floor(p33*K2)*Kz - Kzo;

	vec3 dx11 = Pfx + jitter2x2x2*ox11;
	vec3 dy11 = Pfy.x + jitter2x2x2*oy11;
	vec3 dz11 = Pfz.x + jitter2x2x2*oz11;

	vec3 dx12 = Pfx + jitter2x2x2*ox12;
	vec3 dy12 = Pfy.x + jitter2x2x2*oy12;
	vec3 dz12 = Pfz.y + jitter2x2x2*oz12;

	vec3 dx13 = Pfx + jitter2x2x2*ox13;
	vec3 dy13 = Pfy.x + jitter2x2x2*oy13;
	vec3 dz13 = Pfz.z + jitter2x2x2*oz13;

	vec3 dx21 = Pfx + jitter2x2x2*ox21;
	vec3 dy21 = Pfy.y + jitter2x2x2*oy21;
	vec3 dz21 = Pfz.x + jitter2x2x2*oz21;

	vec3 dx22 = Pfx + jitter2x2x2*ox22;
	vec3 dy22 = Pfy.y + jitter2x2x2*oy22;
	vec3 dz22 = Pfz.y + jitter2x2x2*oz22;

	vec3 dx23 = Pfx + jitter2x2x2*ox23;
	vec3 dy23 = Pfy.y + jitter2x2x2*oy23;
	vec3 dz23 = Pfz.z + jitter2x2x2*oz23;

	vec3 dx31 = Pfx + jitter2x2x2*ox31;
	vec3 dy31 = Pfy.z + jitter2x2x2*oy31;
	vec3 dz31 = Pfz.x + jitter2x2x2*oz31;

	vec3 dx32 = Pfx + jitter2x2x2*ox32;
	vec3 dy32 = Pfy.z + jitter2x2x2*oy32;
	vec3 dz32 = Pfz.y + jitter2x2x2*oz32;

	vec3 dx33 = Pfx + jitter2x2x2*ox33;
	vec3 dy33 = Pfy.z + jitter2x2x2*oy33;
	vec3 dz33 = Pfz.z + jitter2x2x2*oz33;

	vec3 d11 = dx11 * dx11 + dy11 * dy11 + dz11 * dz11;
	vec3 d12 = dx12 * dx12 + dy12 * dy12 + dz12 * dz12;
	vec3 d13 = dx13 * dx13 + dy13 * dy13 + dz13 * dz13;
	vec3 d21 = dx21 * dx21 + dy21 * dy21 + dz21 * dz21;
	vec3 d22 = dx22 * dx22 + dy22 * dy22 + dz22 * dz22;
	vec3 d23 = dx23 * dx23 + dy23 * dy23 + dz23 * dz23;
	vec3 d31 = dx31 * dx31 + dy31 * dy31 + dz31 * dz31;
	vec3 d32 = dx32 * dx32 + dy32 * dy32 + dz32 * dz32;
	vec3 d33 = dx33 * dx33 + dy33 * dy33 + dz33 * dz33;

	// Sort out the two smallest distances (F1, F2)
#if 0
	// Cheat and sort out only F1
	vec3 d1 = min(min(d11,d12), d13);
	vec3 d2 = min(min(d21,d22), d23);
	vec3 d3 = min(min(d31,d32), d33);
	vec3 d = min(min(d1,d2), d3);
	d.x = min(min(d.x,d.y),d.z);
	return sqrt(d.xx); // F1 duplicated, no F2 computed
#else
	// Do it right and sort out both F1 and F2
	vec3 d1a = min(d11, d12);
	d12 = max(d11, d12);
	d11 = min(d1a, d13); // Smallest now not in d12 or d13
	d13 = max(d1a, d13);
	d12 = min(d12, d13); // 2nd smallest now not in d13
	vec3 d2a = min(d21, d22);
	d22 = max(d21, d22);
	d21 = min(d2a, d23); // Smallest now not in d22 or d23
	d23 = max(d2a, d23);
	d22 = min(d22, d23); // 2nd smallest now not in d23
	vec3 d3a = min(d31, d32);
	d32 = max(d31, d32);
	d31 = min(d3a, d33); // Smallest now not in d32 or d33
	d33 = max(d3a, d33);
	d32 = min(d32, d33); // 2nd smallest now not in d33
	vec3 da = min(d11, d21);
	d21 = max(d11, d21);
	d11 = min(da, d31); // Smallest now in d11
	d31 = max(da, d31); // 2nd smallest now not in d31
	d11.xy = (d11.x < d11.y) ? d11.xy : d11.yx;
	d11.xz = (d11.x < d11.z) ? d11.xz : d11.zx; // d11.x now smallest
	d12 = min(d12, d21); // 2nd smallest now not in d21
	d12 = min(d12, d22); // nor in d22
	d12 = min(d12, d31); // nor in d31
	d12 = min(d12, d32); // nor in d32
	d11.yz = min(d11.yz,d12.xy); // nor in d12.yz
	d11.y = min(d11.y,d12.z); // Only two more to go
	d11.y = min(d11.y,d11.z); // Done! (Phew!)
	return sqrt(d11.xy); // F1, F2
#endif
}

// Permutation polynomial: (34x^2 + x) mod 289
vec4 permute4(vec4 x) {
  return mod((34.0 * x + 1.0) * x, 289.0);
}

// Cellular noise, returning F1 and F2 in a vec2.
// Speeded up by using 2x2x2 search window instead of 3x3x3,
// at the expense of some pattern artifacts.
// F2 is often wrong and has sharp discontinuities.
// If you need a good F2, use the slower 3x3x3 version.
vec2 cellular2x2x2(vec3 P) {
#define K 0.142857142857 // 1/7
#define Ko 0.428571428571 // 1/2-K/2
#define K2 0.020408163265306 // 1/(7*7)
#define Kz 0.166666666667 // 1/6
#define Kzo 0.416666666667 // 1/2-1/6*2
	vec3 Pi = mod(floor(P), 289.0);
 	vec3 Pf = fract(P);
	vec4 Pfx = Pf.x + vec4(0.0, -1.0, 0.0, -1.0);
	vec4 Pfy = Pf.y + vec4(0.0, 0.0, -1.0, -1.0);
	vec4 p = permute4(Pi.x + vec4(0.0, 1.0, 0.0, 1.0));
	p = permute4(p + Pi.y + vec4(0.0, 0.0, 1.0, 1.0));
	vec4 p1 = permute4(p + Pi.z); // z+0
	vec4 p2 = permute4(p + Pi.z + vec4(1.0)); // z+1
	vec4 ox1 = fract(p1*K) - Ko;
	vec4 oy1 = mod(floor(p1*K), 7.0)*K - Ko;
	vec4 oz1 = floor(p1*K2)*Kz - Kzo; // p1 < 289 guaranteed
	vec4 ox2 = fract(p2*K) - Ko;
	vec4 oy2 = mod(floor(p2*K), 7.0)*K - Ko;
	vec4 oz2 = floor(p2*K2)*Kz - Kzo;
	vec4 dx1 = Pfx + jitter*ox1;
	vec4 dy1 = Pfy + jitter*oy1;
	vec4 dz1 = Pf.z + jitter*oz1;
	vec4 dx2 = Pfx + jitter*ox2;
	vec4 dy2 = Pfy + jitter*oy2;
	vec4 dz2 = Pf.z - 1.0 + jitter*oz2;
	vec4 d1 = dx1 * dx1 + dy1 * dy1 + dz1 * dz1; // z+0
	vec4 d2 = dx2 * dx2 + dy2 * dy2 + dz2 * dz2; // z+1

	// Sort out the two smallest distances (F1, F2)
#if 0
	// Cheat and sort out only F1
	d1 = min(d1, d2);
	d1.xy = min(d1.xy, d1.wz);
	d1.x = min(d1.x, d1.y);
	return sqrt(d1.xx);
#else
	// Do it right and sort out both F1 and F2
	vec4 d = min(d1,d2); // F1 is now in d
	d2 = max(d1,d2); // Make sure we keep all candidates for F2
	d.xy = (d.x < d.y) ? d.xy : d.yx; // Swap smallest to d.x
	d.xz = (d.x < d.z) ? d.xz : d.zx;
	d.xw = (d.x < d.w) ? d.xw : d.wx; // F1 is now in d.x
	d.yzw = min(d.yzw, d2.yzw); // F2 now not in d2.yzw
	d.y = min(d.y, d.z); // nor in d.z
	d.y = min(d.y, d.w); // nor in d.w
	d.y = min(d.y, d2.x); // F2 is now in d.y
	return sqrt(d.xy); // F1 and F2
#endif
}

float fbm_map( vec3 position )
{
//	float frequency  = 1.0;
//	float lacunarity = 2.0;
//	float gain       = 0.5;
	float fbm=snoise(1.0*position)
 		+ 0.5*snoise(2.0*position)
 		+ 0.25*snoise(4.0*position)
 		+ 0.125*snoise(8.0*position)
 		+ 0.0625*snoise(16.0*position)
 		;
	return fbm;
}

float crater_map( vec3 position, float cutoff, float inset )
{
	vec2 value = cellular2x2x2( position ); // 0..1
	float factor = smoothstep( 0.0, cutoff, value.x );
	float rmaxabs = 1.0 / max( 1.0 - inset, inset );
	return 1.0 - smoothstep( -0.05, 1.0, abs( factor - inset ) * rmaxabs ); 
}

float height_map( vec3 position )
{
	float crater1 = crater_map( 3*position, 0.5, 0.5 ); 
	float crater2 = crater_map( 1.7*position, 0.5, 0.8 ); 
	float crater  = 0.7 * crater2 + 0.3 * crater1;
	float fbm     = ( fbm_map( 6*position ) + 1.0 ) / 2.0;
	//fbm = 1.0;
	return 1.5*crater*(0.15*fbm+0.85) + fbm * 0.2;
}

vec3 color_map( vec3 position, float h )
{
	float fbm = ( fbm_map( 2*position ) + 1.0 ) / 2.0;
	//fbm = 1.0;
	return (fbm*0.5 + 0.5)*vec3(0.9,0.9,1.0);
}

vec3 self_ilum_map( vec3 position, float h )
{
//	float test     = fbm_map( 8*position );
//	float test = crater_map( 1.7*position, 0.5, 0.8 ); 
//	if ( test < -1.0 )
//		return vec3(1.0,0.0,0.0);
//	if ( test > 1.2 )
//		return vec3(0.0,1.0,0.0);
	return vec3(0.0,0.0,0.0);
}

float specular_map( vec3 position, float h )
{
	return 0.2*h;
}

vec3 normal_map( vec3 position, vec3 normal, vec3 tangent, vec3 bitangent )
{
	float step  = 0.001;
	float ystep = 0.1;
	vec3 n = ystep * normal;
	vec3 t = step * tangent;
	vec3 b = step * bitangent;
	vec3 prev_x = position - t;
	vec3 prev_z = position - b;
	vec3 next_x = position + t;
	vec3 next_z = position + b;
	float px = height_map( prev_x );
	float nx = height_map( next_x );
	float pz = height_map( prev_z );
	float nz = height_map( next_z );
	prev_x += px * n;
	prev_z += pz * n;
	next_x += nx * n;
	next_z += nz * n;
	vec3 nt = next_x - prev_x;
	vec3 nb = next_z - prev_z;
	return normalize( cross( nt, nb ) );
}


float world_depth( vec3 world_pos )
{
	vec4 v = nv_m_mvp * vec4( world_pos, 1.0 );
	v.z /= v.w;
	v.z = ( v.z + 1.0 ) * 0.5;
	return v.z;
}

vec4 sphere_isect( vec3 origin, vec3 ray, vec3 center, float r2 )
{
	vec3 sd = center - origin;
	float b = dot( ray, sd );
	float disc = b*b + r2 - dot(sd,sd);
	if ( disc > 0.0 )
	{
		float tnow = b - sqrt(disc);
		return vec4( origin + tnow * ray, tnow );
	}
	return vec4(0,0,0,0);
}

// Minnaert limb darkening diffuse term
float minnaert( vec3 L, vec3 Nf, float k) {
	float ndotl = max( 0.0, dot(L, Nf));
	return pow( ndotl, k);
}

// Ward isotropic specular term
float wardiso( vec3 Nf, vec3 Ln, vec3 Hn, float roughness, float ndotv ) {
	float ndoth = dot( Nf, Hn);
	float ndotl = dot( Nf, Ln);
	float tandelta = tan( acos(ndoth));
	return exp( -( pow( tandelta, 2.0) / pow( roughness, 2.0)))
		* (1.0 / sqrt( ndotl * ndotv ))
		* (1.0 / (4.0 * pow( roughness, 2.0)));
	}
	
float schlick( vec3 Nf, vec3 Vf, float ior, float ndotv ) {
	float kr = (ior-1.0)/(ior+1.0);
	kr *= kr;
	return kr + (1.0-kr)*pow( 1.0 - ndotv, 5.0);
}

void main(void) {
	vec3 view_vector  = normalize( v_world_pos - v_camera_loc );
	vec4 inter = sphere_isect( v_camera_loc, view_vector, center, radius*radius );
	//float height = height_map( inter.xyz );
	//vec3 position = inter.xyz + height*0.1*inter.xyz;
	vec3 position = inter.xyz;
	float height = height_map( position );
	vec3 light_vector = normalize( position - v_light_loc );
	view_vector       = normalize( position - v_camera_loc );
	vec3 normal       = normalize( position - center );
	vec3 tangent      = normalize( cross( normal, vec3(0.0,1.0,0.0) ) );
	vec3 bitangent    = cross( normal, tangent );

	normal       = normal_map( position, normal, tangent, bitangent );
	tangent      = normalize( cross( normal, vec3(0.0,1.0,0.0) ) );
	bitangent    = cross( normal, tangent );

	//vec3 diff_texel      = vec3( texture2D( nv_t_diffuse, v_texcoord ) );	
	vec3 diff_texel       = color_map( position, height );
	float specular_amount = specular_map( position, height );
	vec3 self_ilum_color  = self_ilum_map( position, height );
	vec3 ambient_color    = vec3 (0.1, 0.1, 0.1);
	float diffuse_amount  = 1.0 - specular_amount;

	vec3 reflect_vector     = reflect( light_vector, normal );
	float dot_prod_specular = dot( reflect_vector, -view_vector );
	float dot_prod_diffuse  = dot( -light_vector, normal );

	float diffuse_factor  = max( dot_prod_diffuse, 0.0 );
	float specular_factor = pow( max( dot_prod_specular, 0.0 ), 16.0 ); // 100.0

	float final_diffuse  = diffuse_amount * diffuse_factor;
	float final_specular = specular_amount * specular_factor;


	vec3 N = normalize(normal);
	vec3 V = normalize(view_vector);
	vec3 L = normalize(-light_vector);
	vec3 Vf = -V;
	float ndotv = dot(N, Vf);
	vec3 H = normalize(L+Vf);

	final_diffuse  = minnaert( L, N, 1.5) * diffuse_factor;
	//float fresnel  = schlick( N, V, 0.1, ndotv);
	//final_specular = wardiso( N, L, H, 0.1, ndotv) * fresnel;

	vec3 diffuse_color   = light_diffuse.xyz  * final_diffuse * diff_texel;
	vec3 specular_color  = light_specular.xyz * final_specular * diff_texel;

	if ( inter.w < 0.0 || inter.w > 0.0 )
		gl_FragColor = vec4( max( self_ilum_color, diffuse_color + specular_color + ambient_color ), 1.0 );
	else
		discard;
	gl_FragDepth = world_depth( inter.xyz );
	//gl_FragColor = vec4( 0.0, 1.0, 0.5, 0.8 );

}