uniform vec3 camPos;
uniform vec3 lightDir;
uniform int iterations;
uniform float power;

const float BOUNDING_RADIUS_SQ = 100.0;

bool intersectBoundingSphere(vec3 origin, vec3 direction, out vec2 t) {
	float b = dot(origin, direction);
	float c = dot(origin, origin) - 1.25*1.25;
	float d = b*b - c;

	float sqrtD = sqrt( d );
	t.x = -b - d;
	t.y = -b + d;

	return d > 0.0;
}

bool iterateMandelbulb(vec3 position, out vec4 orbitTrap, out float potential) {
	vec3 cz = position;
	vec3 c = position;

	orbitTrap = vec4(100.0);
	for(int i = 1; i <= iterations; i++) {
		//Faster version of regular Mandelbulb by Iñiqo Quilez
		float x = cz.x;
		float y = cz.y;
		float z = cz.z;

		float x2 = x*x;
		float y2 = y*y;
		float z2 = z*z;

		float x4 = x2*x2;
		float y4 = y2*y2;
		float z4 = z2*z2;

		float k3 = x2 + z2;
		float k2 = inversesqrt(k3*k3*k3*k3*k3*k3*k3*k3);
		float k1 = x4 + y4 + z4 - 6.0*y2*z2 - 6.0*x2*y2 + 2.0*z2*x2;
		float k4 = x2 - y2 + z2;

		cz.x = + 64.0*x*y*z*(x2 - z2)*k4*(x4 - 6.0*x2*z2 + z4)*k1*k2;
		cz.y = - 16.0*y2*k3*k4*k4 + k1*k1;
		cz.z = - 8.0*y*k4*(x4*x4 - 28.0*x4*x2*z2 + 70.0*x4*z4 - 28.0*x2*z2*z4 + z4*z4)*k1*k2;

		/*float wr = sqrt(dot(cz,cz));
		float wo = acos(cz.y/wr);
		float wi = atan(cz.x,cz.z);
		wr = pow( wr, power );
		wo = wo * power;
		wi = wi * power;
		cz.x = wr*sin(wo)*sin(wi);
		cz.y = wr*cos(wo);
		cz.z = wr*sin(wo)*cos(wi);*/

		float distance = dot(cz, cz);
		orbitTrap = min(orbitTrap, vec4(cz*cz, distance));

		if(distance > BOUNDING_RADIUS_SQ) {
			potential = 0.5*log(distance)/pow(power, float(i));

			return false;
		}

		cz += c;
	}

	potential = 0.0;

	return true;
}

bool raymarch(vec3 origin, vec3 direction, out vec3 normal, out vec4 color) {
	vec2 intersectionT;
	if( !intersectBoundingSphere( origin, direction, intersectionT ) ) {
		return false;
	}

	if(intersectionT.x < 0.001) {
		intersectionT.x = 0.001;
	}

	vec3 gradient;
	for( float t = intersectionT.x; t < intersectionT.y; ) {
		float potential;

		vec3 position = origin + direction*t;

		if(iterateMandelbulb(position, color, potential)) {
			normal = normalize( gradient );

			return true;
		}

		float minimumStep = clamp(0.00001*t*0.7, 0.00001, 0.005);
		float epsilon = minimumStep*0.1;

		float potential1, potential2, potential3;

		vec4 dummyColor;

		iterateMandelbulb(position + vec3(epsilon, 0.0, 0.0), dummyColor, potential1);
		iterateMandelbulb(position + vec3(0.0, epsilon, 0.0), dummyColor, potential2);
		iterateMandelbulb(position + vec3(0.0, 0.0, epsilon), dummyColor, potential3);

		gradient = vec3(potential1 - potential, potential2 - potential, potential3 - potential);

		float stepSize = 0.5*potential*epsilon/length(gradient);

		if(stepSize < minimumStep) {
			normal = normalize( gradient );

			return true;
		}
		t += stepSize;
	}

	return false;
}

void main(void) {
	vec3 direction = normalize( vec3( gl_TexCoord[ 0 ] ) );
	vec3 position  = camPos;

	vec4 color;
	vec3 normal;

	if(raymarch(position, direction, normal, color)) {
		vec3 lighting = vec3(color.w);

        lighting = 1.5*(lighting*0.15 + 0.85*sqrt(lighting));

		gl_FragColor = vec4(lighting, 1.0);
	} else {
		gl_FragColor = vec4(0.2, 0.2, 0.2, 1.0);
	}
}
