//created by Skye Adaire

#define tau32 6.2831853072
#define eps32 1e-10

float alpha(float x, float a, float b)
{
   return (x - a) / (b - a);
}

#define uclamp(x) clamp(x, 0.0, 1.0)
#define ualpha(x, a, b) uclamp(alpha(x, a, b))

#define Real float
#define Complex vec2

Complex H_conjugate(Complex h)
{
    return Complex(h[0], -h[1]);
}

Real H_sqnorm(Complex h)
{
    return dot(h, h);
}

Real H_norm(Complex h)
{
    return length(h);
}

Complex H_inverse(Complex h)
{
    return H_conjugate(h) / H_sqnorm(h);
}

Complex H_multiply(Complex lhs, Complex rhs)
{
	return Complex(
        lhs.x * rhs.x - lhs.y * rhs.y, 
        lhs.x * rhs.y + lhs.y * rhs.x);
}

Complex H_sin(Complex h)
{
	return Complex(sin(h.x) * cosh(h.y), cos(h.x) * sinh(h.y));
}

Real H_argument(Complex h)
{
   return atan(h.y, h.x);//[-pi, pi]
}

Real H_argument2(Complex h)
{
    Real angle = H_argument(h);
    return angle < Real(0) ? angle + tau32 : angle;//[0, tau]
}

Complex H_versor(Real angle)
{
    return Complex(cos(angle), sin(angle));
}

struct PolarComplex
{
    float norm;
    float argument;
};

PolarComplex H_toPolar(Complex h)
{
    return PolarComplex(H_norm(h), H_argument(h));
}

//the hyperplane centered at the origin with normal 0,0,1
bool intersectHyperplane(vec3 rayPosition, vec3 rayDirection, out float t)
{
    t = -rayPosition[2] / rayDirection[2];

    return t > 0.0 && !isinf(t);
}

//return the first positive solution along the ray
bool solveQuadraticIntersection(float a, float b, float c, out float t)
{
    if(abs(a) < eps32)
    {
        t = -c / b;
        return t > 0.0;
    }

	float discriminant = b * b - 4.0 * a * c;

    if(abs(discriminant) < eps32)
    {
        t = - b / (2.0 * a);
        return true;
    }
    else if(discriminant < 0.0)
    {
        return false;
    }
    else
	{
        float sqrtd = sqrt(discriminant);

        float t0 = (-b + sqrtd) / (2.0 * a);
        float t1 = (-b - sqrtd) / (2.0 * a);

        if(t1 < t0)
        {
            float tt = t0;
            t0 = t1;
            t1 = tt;
        }

        if(t0 > 0.0)
        {
            t = t0;
            return true;
        }

        if(t1 > 0.0)
        {
            t = t1;
            return true;
        }

        return false;
	}
}

//the hypersphere centered at the origin with radius 1
bool intersectHypersphere(vec3 rayPosition, vec3 rayDirection, out float t)
{
    float a = dot(rayDirection, rayDirection);
    float b = 2.0 * dot(rayDirection, rayPosition);
    float c = dot(rayPosition, rayPosition) - 1.0;

    return solveQuadraticIntersection(a, b, c, t);
}

//https://www.shadertoy.com/view/lsS3Wc
vec3 hsv2rgb( in vec3 c )
{
    vec3 rgb = clamp( abs(mod(c.x*6.0+vec3(0.0,4.0,2.0),6.0)-3.0)-1.0, 0.0, 1.0 );
    return c.z * mix( vec3(1.0), rgb, c.y);
}

//http://blog.hvidtfeldts.net/index.php/2012/03/lifted-domain-coloring/
vec3 liftedDomainColor(vec2 z)
{
    PolarComplex polar = H_toPolar(z);

    float magnitude = (1.0-1.0/pow(2.0,polar.norm)) * 0.9 + 0.1;
    float logradius = log(polar.norm);

    //black rings
    float fractlog = fract(logradius);
    float ringdist = min(abs(fractlog-0.5), fractlog > 0.5 ? 1.0-fractlog : fractlog);
    float ring = (1.0 - smoothstep(0.00, 0.02, ringdist)) * 0.8;

    //white rays
    float k = 12.0;
    float sectorsize = (tau32) / k;
    float anglemod = mod(polar.argument, sectorsize);
    float sectordist = anglemod > sectorsize/2.0 ? sectorsize-anglemod : anglemod;
    float raywidth = 0.02;
    float ray = (1.0 - smoothstep(0.0, raywidth, sectordist)) * 0.8;

    //infinity will be white
    float infinityFade = pow(magnitude,100000000.0);

    //growth ring shade
    float growth = (fractlog)*0.7 + 0.3;
    float darkening = uclamp(1.5*magnitude * (fractlog*0.5 + 0.5) + ray + infinityFade);

    float hue = polar.argument/tau32;
    float saturation = 1.0 - infinityFade;
    float value = darkening;
    
    vec3 color = hsv2rgb(vec3(hue, saturation, value));
    color = mix(color, vec3(1.0), darkening * ray);
    color = mix(color, vec3(0.0), darkening * (ring-infinityFade));

    return color;
}

mat3 rotationXY(float a)
{
    float c = cos(a);
    float s = sin(a);
    return mat3(c, s, 0, -s, c, 0, 0, 0, 1);
}

mat3 rotationXZ(float a)
{
    float c = cos(a);
    float s = sin(a);
    return mat3(c, 0, -s, 0, 1, 0, s, 0, c);
}

mat3 rotationYZ(float a)
{
    float c = cos(a);
    float s = sin(a);
    return mat3(1, 0, 0, 0, c, s, 0, -s, c);
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    float aspectRatio = iResolution.x / iResolution.y;
    vec2 uv = fragCoord/iResolution.xy;
    vec2 clip = uv * 2.0 - 1.0;
    vec2 unitSpacePosition = 0.5 * clip;
    vec2 ratioSpacePosition = vec2(aspectRatio, 1) * unitSpacePosition;

    mat3 viewTransform = rotationXY(-tau32 / 8.0) * rotationYZ(tau32/6.0);
    
    //view basis
    vec3 viewPosition = viewTransform * vec3(0, 0, 6);
    vec3 viewRight = viewTransform * vec3(1, 0, 0);
    vec3 viewUp = viewTransform * vec3(0, 1, 0);
    vec3 viewForward = viewTransform * vec3(0, 0, -1);
    
    //view ray
    vec3 frustumPoint = viewPosition - viewForward;
    vec3 screenRayPosition =
       viewPosition +
       viewRight * ratioSpacePosition.x +
       viewUp * ratioSpacePosition.y;
    vec3 screenRayDirection = normalize(screenRayPosition - frustumPoint);
       
    //transform sphere local to global
    float angle = iTime * 0.4;
    mat3 rotation = rotationXZ(angle);
    vec3 translation = vec3(0, 0, 1.0 + cos(iTime*0.4 + tau32/2.0));
    vec3 pole = translation + vec3(0, 0, 1);
    
    //transform sphere global to local
    mat3 rotationInv = rotationXZ(-angle);
    vec3 translationInv = -translation;
    vec3 poleInv = vec3(0, 0, 1);
    
    //transform the ray to local sphere space
    vec3 rpt = rotationInv * (translationInv + screenRayPosition);
    vec3 rdt = rotationInv * screenRayDirection;
    
    float t, tmin;
    int object = -1;
    vec3 position;
    
    //sphere in local space
    if(intersectHypersphere(rpt, rdt, t))
    {
        tmin = t;
        object = 1;
    }
    
    //plane in global space
    if(intersectHyperplane(screenRayPosition, screenRayDirection, t))
    {
        if(object == -1 || t < tmin)
        {
            tmin = t;
         	object = 2;
        }
    }
    
    vec3 color = vec3(0);
    
    if(object == 2)//we are on the global plane, project to the local sphere
    {
        vec3 rp = screenRayPosition + tmin * screenRayDirection;
        vec3 rd = normalize(rp - pole) * (length(rpt) < 1.0 ? 1.0 : -1.0);;
        
    	rpt = rotationInv * (translationInv + rp);
    	rdt = rotationInv * rd;
        
        intersectHypersphere(rpt, rdt, tmin);//must hit
        object = 1;
    }
    
    if(object == 1)//we are on the local sphere, project to the local plane
    {
        vec3 rp =  rpt + tmin * rdt;
        vec3 rd = normalize(rp - poleInv);

        float t;
        intersectHyperplane(rp, rd, t);//must hit
        
        vec2 localPlanePosition = (rp + t * rd).xy;
        
        //apply a complex transform to the plane
        vec2 z = localPlanePosition;
        vec2 zt = H_sin(z);

        color = liftedDomainColor(zt);
    }
    
    fragColor = vec4(color,1.0);
}