#define AA 2  //set to 1 if too slow, increase if you got a phat GPU!

//the story is I have my own webgl browser thingy where I can interactively
//set all these parameters via javascript, but I am exporting them here
//so I can show you the interesting ones
float params[38] = float[38]
(0.,-0.131,0.,0.,1.1,0.,0.11,1.9,1.5, 3.,
 15.,0.035,-0.6,0.,-2.,0.,0.,4.,5., 4.,
 -0.20,0.6,0.,0.,0.,0.,0.,-2.,0.,0.,
 0.,0.,0.,2.,-0.01,0.01,0.,0.6);

const float tmax = 20.0;
vec3 color;
vec3 ray_origin;  
float inside_of_sphere( vec3 p, float r ) {
    return -( length(p) - r ); //negating the usual SDF gets us on the inside
}

float map(vec3 p) {

    // the guts of this came from Space Monolith shadertoy by zackpudil
    // but I parameterized it to death and changed the coloring and lighting
    // and found this sierpinski gasket hanging in the columns
    vec3 orig_p = p;

    p.xy = mod(p.xy + 1.0, 3.0) - 1.0;

    p.z = abs(p.z) - 0.75 - params[24];

    vec4 q = vec4(p, 1.0);
    
    float mscale = 2.*(1.+params[1]);
    float clamp_min = .5 - params[0];
    float clamp_max = 1. + params[0];

    color = vec3(0.);
    float color_radius = params[25];
    float num_iter = 0.;
    for(int i = 0; i < 100; i++) {

        if ( float(i) > ( params[10] ) ) break;

        float ilength = length(q.xyz);  //using the lagging q.xyz yield nice results

        q.xyz = abs(q.xyz) - vec3(0.3, 1.0, -0.0) + vec3(params[21], params[22], params[23]);

        if (params[29] == 0. ) 
            ilength = length(q.xyz - vec3( params[30], params[31], params[32]) );   


        q = mscale*q/clamp( pow(ilength,2.+params[28]), clamp_min, clamp_max) 
            - vec4(1.0+params[34], 0.0+params[35], 0.3+params[36], 0.0);

      
        // I have found that nice colors can be generated for iterated
        // systems that send points flying in 3d by keeping count of how 
        // many times an orbit passes through a big block of subspace
        // then using cosine to transform it - see down below
        if      ( q.x*q.y > color_radius ) { color.x ++ ; }
        else if ( q.y*q.z > color_radius ) { color.y ++ ; }
        else if ( q.z*q.x > color_radius ) { color.z ++ ; }
        

        num_iter ++;
    }

    color /= (1. + params[37] * num_iter);
    
    return max( abs(q.x + q.y + q.z)/q.w, inside_of_sphere(orig_p - ray_origin, params[13]) );

    //return ( length(q.xyz)/q.w );
}

float march(vec3 ro, vec3 rd) {

    float d = 1.0;
    float e, t = 0.0;

    ray_origin = ro;  //needed for inside of sphere

    float detail = .001 * (1.+params[20]); //of course detail increases as this variable decreases
    for(int i = 0; i < 200; i++) {
        if(abs(d) < e || t >= tmax) break;
        d = map(ro + rd*t);
        e = detail * (1.0 + t*4.0);  //tolerance for hitting an object increases with distance to reduce background noise
        t += d*(0.55 + 0.05*t);
    }
    
    return t;
}

vec3 normal(vec3 p) {
    vec2 h = vec2(0.01 * (1.+params[12]), 0.0);
    vec3 n = vec3(
        map(p + h.xyy) - map(p - h.xyy),
        map(p + h.yxy) - map(p - h.yxy),
        map(p + h.yyx) - map(p - h.yyx)
    );
    return normalize(n);
}


float hash2(vec2 n) {
    return fract(sin(dot(n, vec2(27.233, 71.989)))*43758.5453);
}

float calcAO( vec3 pos, vec3 nor ) {
    float detail = .001;
	float aodet=detail*8.;  //40.
	float totao = 0.0;
  	float sca = 25.0;
  	for( int aoi=0; aoi<40; aoi++ ) {
        if ( float(aoi) > params[33]) break;
		float hr = aodet*float(aoi*aoi);        
		vec3 aopos =  nor * hr + pos;
		float dd = map( aopos );
		totao += -(dd-hr)*sca;
		sca *= 0.7;
    	}
    	return clamp( 1.0 - 2.*totao, 0.2, 1.0 );
}



void mainImage( out vec4 fragColor, in vec2 fragCoord )
{


    vec3 tot = vec3(0.0);
    
    //anti aliasing taken from Soft Shadow Variation shadertoy by IQ
#if AA>1
    for( int m=0; m<AA; m++ )
    for( int n=0; n<AA; n++ )
    {
            // pixel coordinates
        vec2 o = vec2(float(m),float(n)) / float(AA) - 0.5;
        vec2 uv = (-iResolution.xy + 2.0*(fragCoord+o))/iResolution.y; 
    
#else    
        vec2 uv = (-iResolution.xy + 2.0*fragCoord)/iResolution.y;
#endif


    vec3 col = vec3(0.); //vec3(1)*step(0.997, hash2(uv));
    
    vec3 ro = vec3( 0., 0., max(min(2. + 3.*sin(iTime/3.), 4.),-.5) ); 
 
    float  focal_point = params[14];
    vec3 rd =  normalize(vec3(uv,focal_point));

    vec3 camera_direction = vec3( 0., 0., 1. );  
    
    vec3 ltDir = normalize ( camera_direction + vec3(0.,params[4]+2.*cos(iTime/2.),0.) );  

    float i = march(ro, rd);

    if(i < tmax) {

        vec3 pos = ro + rd*i;
        vec3 nor = normal(pos);

        col = vec3(0.0);

        float bounce = clamp( 1. + dot( rd, nor ), 0., 1. );

        col += (1.+params[8])*0.1*pow(bounce, 2.0);

        //specular
        col += (1.0+params[9]) * pow (max (0., dot (ltDir, reflect (rd, nor))), 32.*params[6]);

        float cx = cos(color.x*params[17]);
        float cy = cos(color.y*params[18]);
        float cz = cos(color.z*params[19]);

        if ( params[7] != 0. )
            col += .1 * params[7] * vec3( cx, cy, cz ); 

        //fade with distance
        col *= smoothstep(-8.0, -6.0*(1.+params[27]), -i);

        col *= calcAO( pos, nor );

        col = clamp( col, 0., 1.); 
    }

    col = pow(col, vec3(.4154545 * (1.+params[26])  ) );
    tot += col;
    
#if AA>1
    }
    tot /= float(AA*AA); //average color for anti aliasing
#endif
    
    fragColor = vec4(tot, 1.);
       
   
}