#define MAXSTEP 100
#define MAXDIST 100.
#define MINSURF 0.0001
#define pi 3.14159
#define ITER 21.

mat2 Rot(float b)
{
    float a = b*pi/180.;
    float s = sin(a);
    float c = cos(a);
    return mat2(c,-s,s,c);
}
vec2 DE(vec3 pos) {

    pos.xz*= Rot(-25.);
    pos.yz*=Rot(-40.);


    float Iterations = ITER;
    float Bailout = 40.;
    float Power = 8.;

    vec3 trap = vec3(0,0,0);
    float minTrap = 1e10;

	vec3 z = pos;
	float dr = 1.;
	float r = 0.0;

    float orb = 0.;
	for (float i = 0.; i < Iterations ; i++) {
        orb = i;
		r = length(z);
		if (r > Bailout) break;

        minTrap = min(minTrap, z.z);

		float theta = acos(z.z/r);
		float phi = atan(z.y,z.x);
		dr =  pow(r, Power-1.0)*Power*dr + 1.0;

		float zr = pow(r,Power);
		theta = theta*Power;
		phi = phi*Power;

        z = zr*vec3(cos(theta)*cos(phi), cos(theta)*sin(phi), sin(theta));
		z+=pos;
	}
	return vec2(0.5*log(r)*r/dr,orb);
}
vec2 rayMarch(vec3 ro, vec3 rd) {
    float surfDist = DE(ro).x/800.; //better surface
    float dO = 0.;
    int i = 0;
    float iit =float(MAXSTEP-1);
    vec2 t;
    do
    {
        vec3 p = ro+rd*dO;

        t = DE(p);
        dO+=t.x;
        if(dO>MAXDIST || t.x < surfDist) { iit = float(i);break; }
        i++;
    }while(i<MAXSTEP);

    return vec2(iit,t.y);
}

vec3 GetNormal(vec3 p) {
	float d = DE(p).x;
    vec2 e = vec2(.001, 0);

    vec3 n = d - vec3(DE(p-e.xyy).x,DE(p-e.yxy).x,DE(p-e.yyx).x);

    return normalize(n);
}
vec3 colors(float orbit) {
    float irit = ITER;
    vec3 color = vec3(1.,1.,1.);
    if(orbit < irit)
        color = vec3(1.,0.,0.);
    if(orbit < irit*0.875)
        color = vec3(1.,1.,0.);
    if(orbit < irit*0.75)
        color = vec3(0.,1.,0.);
    if(orbit < irit*0.625)
        color = vec3(0.,1.,1.);

    if(orbit < irit*0.5)
        color = vec3(0.,0.,1.);
    if(orbit < irit*0.375)
        color = vec3(1.,0.,1.);

    if(orbit < irit*0.25)
        color = vec3(0.5,1.,0.);
    if(orbit < irit*0.125)
        color = vec3(1.,0.5,0.);
    return color;
}
float zoomEX(float c, float start, float end,float speed) { //c is curve between 0 and 1, start is starting position, end is ending position
    float b = end;
    float a = start-end;
    return a*pow(c,iTime*speed)+b; //y=a*c^x+b

}//exponential zoom
void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    //Normalized coordinates (0 - 1)
    vec2 uv = fragCoord/iResolution.xy*2.-1.;

    float ttime = 1.0035394;
    vec3 ro = vec3(0,0,zoomEX(0.4,-2.,-2.+ttime,0.5));
    //vec3 ro = vec3(0,0,-2.+ttime);
    vec3 rd = normalize(vec3(uv.x,uv.y,1.));
    vec2 d = rayMarch(ro,rd);
    vec3 norm = GetNormal(ro+rd*d.x);
    // Output to screen
    vec3 p = ro+rd*d.x;
    vec3 bColor = colors(d.y);
    fragColor = vec4(bColor*(1.-pow(d.x/float(MAXSTEP),0.9)),1.0);
}