// Created by Patrik Colling - cyperus/2020 (https://www.youtube.com/user/cyperquantus)
// License Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.
//
// Parametric procedural 3d-object allowing smooth transitions between
// spheres, tori, torusknots.
//

#define AA

#define PI 3.14159265359
// Time modulation
#define PIiTime PI * iTime
#define sinTime02 sin(0.2*iTime)
#define sinTime03 sin(0.3*iTime)
#define sinTime05 sin(0.5*iTime)
#define sinTime07 sin(0.7*iTime)

// complex number transformations
// z_out = za*zb
vec2 cmul(vec2 za,vec2 zb){
    return za*mat2(zb.x,-zb.y,zb.yx);
    }
// z*e^ia
vec2 crot(vec2 z, float a){
    return cmul(z, vec2(cos(a),sin(a)));
    }
// z_out = z^q
vec2 cpowq(vec2 z, float q){
    float r = pow(length(z), q); float a=q* atan(z.y,z.x); return vec2(r*cos(a),r*sin(a));
    }
// z_out = (z*e^ia)^p-x0
vec2 ctfr(vec2 z,float x0, float p, float a){
	return cpowq( cmul(z, vec2(cos(-a),sin(-a))),p) - vec2(x0, 0.);
    }

// came
mat3 camerabase(vec3 co, vec3 ct, vec3 cup){
    // co	: camera origin point in worldspace
    // cup	: camera up direction vector in worldspace
    // ct	: camera target point in worldspace
	vec3 cw = normalize(ct - co);		// camera ponting direction
	vec3 cu = normalize(cross(cup, cw));// camera left right
	vec3 cv = normalize(cross(cw, cu));	// camera down up
    return mat3(cu,cv,cw); // return camera orhtogonal basis as matrix
}

vec3 cameraraydirection(vec2 uv, mat3 cam, float f){
    // uv : Viewport coordinates
    // cam : camera orhtogonal basis
    // f : focal length zoom-in: abs(f) ==> 0., zoom-out: abs(f) ==> +inf.
    return normalize(cam * vec3(f*uv,1.));
}

vec4 df(vec3 p)
{	///// MetaTorus-transformation

    // component swizzling: openGL => math coordinatesystem notation.
    p = p.zxy; //(z, x, y) => (x, y, z)

    //// (torus,sphere,2spheres)-transformation
    // cylindrical coordinate system
    float au = atan(p.y, p.x); // [-PI, +PI]
    float rxy = length(p.xy); // [0., +inf]
    // 3D => 2D-space: radial half plane in cylindrical coordinate system == half complex plane
    vec2 z = vec2(rxy, p.z);
    // julia, half complex plane => complex plane
	float shift0 = 1.4+(1.0+sinTime05); // [-inf, +inf]
    z = cmul(z,z)- vec2(shift0, 0.); // realaxis-translation, fraction == 2

    //// juliafractal-transformation
	/// julia fractal iteration 0
    const float fracu1 = 3.0; // int in [1,2,3,..]
    const float fracv1 = 2.0; // int in [1,2,3,..]
    const float twist1 = 0.0; // int in [...,-1,0,+1,...]
    // shiftX modulation with FourierSerie F(fracu1 * au)
    float shift1	= 2.0 * sinTime02 												// 0. order
        			+ 0.5 * cos(1. * (fracu1 * au) + (0.5 * PIiTime) );				// 1. order
            		//+ 0.25 * cos(2. * (fracu1 * au) + (0.5 +0.01 * PIiTime) ); 	// 2. order
    // torsionX must meet knot condition.
    float torsion1 = au * twist1 / fracv1 - 0.25 * PIiTime;
    z = ctfr(z, shift1, fracv1,  torsion1); // realaxis-translation,fraction,rotation
	/// julia fractal iteration 1
    const float fracu2 = 7.0;
    const float fracv2 = 3.0;
    const float twist2 = 5.0;
    float shift2 = 0.5 +1.4 * sinTime03 + 0.5 * cos(fracu2 * au +0.0- 0.02 * PIiTime);
    float torsion2 = au * twist2 / fracv2 + 0.05 * PIiTime;
    z = ctfr(z, shift2, fracv2, torsion2);

	/// estimated ray-step-length
	float d = log(length(z)); // :( but it works!
	return vec4(d , z, au * fracu1 * fracu2);
}

vec3 nor( vec3 p, float prec )
{
    vec2 e = vec2( prec, 0.0 );
    vec3 n = vec3(
		df(p+e.xyy).x - df(p-e.xyy).x,
		df(p+e.yxy).x - df(p-e.yxy).x,
		df(p+e.yyx).x - df(p-e.yyx).x );
    return normalize(n);
}

// from iq code
float softshadow( in vec3 ro, in vec3 rd, in float mint, in float tmax )
{
	float res = 1.0;
    float t = mint;
    for( int i=0; i<1; i++ )
    {
		float h = df( ro + rd*t ).x;
        res = min( res, 8.0*h/t );
        t += h*.25;
        if( h<0.001 || t>tmax ) break;
    }
    return clamp( res, 0., 1. );
}

// from iq code
float calcAO( in vec3 pos, in vec3 nor )
{
	float occ = 0.0;
    float sca = 1.0;
    // antialeasing
    for( int i=0; i<1; i++ )
    {
        float hr = 0.01 + 0.12*float(i)/4.0;
        vec3 aopos =  nor * hr + pos;
        float dd = df( aopos ).x;
        occ += -(dd-hr)*sca;
        sca *= 0.95;
    }
    return clamp( 1.0 - 3.0*occ, 0.0, 1.0 );
}

//--------------------------------------------------------------------------

// Grab all sky information for a given ray from camera
// from Dave Hoskins // https://www.shadertoy.com/view/Xsf3zX
vec3 GetSky(in vec3 rd, in vec3 sunDir, in vec3 sunCol)
{
	float sunAmount = max( dot( rd, sunDir), 0.0 );
	float v = pow(1.0-max(rd.y,0.0),6.);
	vec3  sky = mix(vec3(.1, .2, .3), vec3(.32, .32, .92), v);
	sky = sky + sunCol * sunAmount * sunAmount * .25;
	sky = sky + sunCol * min(pow(sunAmount, 800.0)*1.5, .3);
	return clamp(sky, 0.0, 1.0);
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec3 tot = vec3(0.0);
#ifdef AA
	vec2 rook[4];
    rook[0] = vec2( 1./8., 3./8.);
    rook[1] = vec2( 3./8.,-1./8.);
    rook[2] = vec2(-1./8.,-3./8.);
    rook[3] = vec2(-3./8., 1./8.);
    for( int n=0; n<4; ++n )
    {
    //// Pixel2Viewport-space
    vec2 o = rook[n];
    vec2 uv = (2.*(fragCoord+o)-iResolution.xy)/iResolution.y;
#else //AA
	vec2 uv = (2.*fragCoord-iResolution.xy)/iResolution.y;
#endif //AA

    // mouse
    vec2 m = (2.*iMouse.xy-iResolution.xy)/iResolution.y;

    //// camera, Viewport2World-space
    m *= 0.49*PI; float aa = m.x; float ab = -m.y;
    vec3 co = 4.*vec3( cos(ab)*sin(aa),sin(ab),cos(ab)*cos(aa)); // camera origin
	const vec3 cup = vec3(0,1,0);	    // camera up direction in Worldspace
	const vec3 ct = vec3(0,0,0);	    // camera target point in Worldspace
	mat3 cam = camerabase(co,ct,cup);   // cameradirection = cam[0];

    const float f = 0.6;
    vec3 rd =  cameraraydirection(uv,cam, f);

	// ray marching
    const float rayDistMax = 14.;	// maximal ray distance length
    const float rayiMax = 300.;		// maximal ray marching iterations
    const float rayF = 0.02;		// ray step multiplier
    const float rayEpsi = 0.0001;
	float rayStep = 1.; // ray step length
    float rayDist = 0.; // ray distance length
	for (float i=0.; i<rayiMax; i++){
		vec3 p = co + rd *rayDist; // point in 3D worldspace
		if (abs(rayStep)<rayEpsi || rayDist>rayDistMax) break;
		rayStep = df(p).x;
		rayDist += rayStep * rayF;
	}

	//// World2Color-space
    const vec3 ld = vec3(0.,1.,.5); // ligth direction
    const vec3 lc = vec3(0.4); // ligth color
    vec3 col;
    vec3 sky = GetSky(rd, ld, lc);

	if (rayDist<rayDistMax)  // intersection with object
	{
        vec3 p = co + rd * rayDist;
		vec3 n = nor(p, 0.001);	// compute normale
		vec4 mat = df(p);
        // uv-coords
        vec2 mat_uv;
        mat_uv.x = atan(mat[2],mat[1])/PI;
        mat_uv.y = mat[3]/PI;
        // texture color
        const float smoothness = 0.05; // TODO: (2) Calculation inside df()?
		col.rgb = vec3(smoothstep(-smoothness,smoothness,cos(PI*mat_uv.x)*cos(PI*mat_uv.y)));
		// lighting
		float occ = calcAO( p, n );
        float amb = clamp( 0.5+0.5*n.y, 0.0, 1.0 );
        float dif = clamp( dot( n, ld ), 0.0, 1.0 ) * (df(p+n*1.16).x);
        float spe = pow(clamp( dot( rd, ld ), 0.0, 1.0 ),16.0);
        float sss = df(p - n*0.001).x/0.01;
        // shading
        dif *= softshadow( p, ld, 0.1, 1. );
        vec3 brdf = vec3(0.0);
        brdf += 0.5*dif*vec3(1.00,0.90,0.60);
        brdf += 0.5*spe*vec3(0.8,0.60,0.20)*dif;
        brdf += 0.3*amb*vec3(0.40,0.60,0.40)*occ;
        brdf += 0.4;
        col.rgb *= brdf;
        // fog
        col.rgb = mix( col.rgb, sky, 1.0-exp( -0.02*rayDist*rayDist ) );
	}
	else
	{// no intersection => background
		col.rgb = sky;
	}

    tot += col;

#ifdef AA
    }
    tot /= 4.;
#endif
    // gamma correction
	fragColor = vec4( sqrt(0.4*tot), 1.0 );
}