/* Creative Commons Licence Attribution-NonCommercial-ShareAlike
   phreax 2022
*/

#define PI 3.141592
#define TAU (2.*PI)
#define SIN(x) (sin(x)*.5+.5)
#define PHI 1.618033988749895
#define hue(v) ( .6 + .6 * cos( 6.3*(v) + vec3(0,23,21) ) )
//#define KALEIDOSCOPIC
//#define POLARKALEI
#define THICKNESS .01
#define AA 0
#define BUMP_EPS 0.004


float tt, g_mat;

mat2 rot(float a) { return mat2(cos(a), -sin(a), sin(a), cos(a)); }

float fadeInOut(float t1, float t2, float fadeTime) {
    return smoothstep(t1, t1+fadeTime, iTime)-smoothstep(t2-fadeTime, t2, iTime);
}

// zucconis spectral palette https://www.alanzucconi.com/2017/07/15/improving-the-rainbow-2/
vec3 bump3y (vec3 x, vec3 yoffset)
{
    vec3 y = 1. - x * x;
    y = clamp((y-yoffset), vec3(0), vec3(1));
    return y;
}


vec3 spectral_zucconi6(float x) {
    x = fract(x);
    const vec3 c1 = vec3(3.54585104, 2.93225262, 2.41593945);
    const vec3 x1 = vec3(0.69549072, 0.49228336, 0.27699880);
    const vec3 y1 = vec3(0.02312639, 0.15225084, 0.52607955);
    const vec3 c2 = vec3(3.90307140, 3.21182957, 3.96587128);
    const vec3 x2 = vec3(0.11748627, 0.86755042, 0.66077860);
    const vec3 y2 = vec3(0.84897130, 0.88445281, 0.73949448);
    return bump3y(c1 * (x - x1), y1) + bump3y(c2 * (x - x2), y2) ;
}

// Amazing fractal texture from jarble https://www.shadertoy.com/view/csl3zl
vec2 triangle_wave(vec2 a){

    vec2 a2 = vec2(1.,0.5),

    a1 = a-a2;

    return abs(fract((a1)*(a2.x+a2.y))-.5);
}


float rect( vec2 p, vec2 b, float r ) {
    vec2 d = abs(p) - (b - r);
    return length(max(d, 0.)) + min(max(d.x, d.y), 0.) - r;
}


vec3 kalei(vec3 p) {
    float iter = 3.;
    p.yz *= rot(tt);
    p.xy = abs(p.xy)+.9;

    #ifdef POLARKALEI
    p.xy = vec2(length(p.xy)*.8, atan(p.y, p.x));
    #endif

    float s= 1.;
    for(float i=0.; i< iter; i++) {
        p = abs(p);
        p.xz *= rot(i/iter*PI+.2*tt);
        p.xz += .2 - .1*i + sin(p.y*.2) + .2*exp(.01+p.z*.2);
        p.z -= .2*sin(p.z*.5*i/iter*2.*PI);
        p -= .2;

    }

    return clamp(p, -1e5, 1e5);

}

// fold space for toroid
vec3 transform(vec3 p) {
    p.yz *= rot(PI*.25 + .1*tt);
    p.xz *= rot(.2*tt);

    #ifdef KALEIDOSCOPIC
    p = kalei(p);
    #endif

    float r = 1.9;
    vec2 cp = vec2(length(p.xz)-r, p.y);


    float rev = 2.5;
    float a = atan(p.z, p.x);


    cp *= rot(rev*a+.3*tt);
    cp= abs(cp) - mix(.1, 1., SIN(tt));
    cp *= rot(.5*tt);


    return vec3(cp, p.z);
}


float map(vec3 p) {

    vec3 bp = p;
    float edge = 0.01;

    vec2 cp = transform(p).xy;
    float dr = rect(cp.xy, vec2(.8, THICKNESS), edge);

    // g_mat = db < dr ? 1. : 0.;
    float d = dr;
    return .2*d;
}


vec3 getNormal(vec3 p) {

    vec2 eps = vec2(0.001, 0.0);
    return normalize(vec3(map(p + eps.xyy) - map(p - eps.xyy),
                          map(p + eps.yxy) - map(p - eps.yxy),
                          map(p + eps.yyx) - map(p - eps.yyx)
                         )
                     );
}

float bumpSurf3D( in vec3 p){

    p.z += .3*tt;
    p = abs(mod(p*2., 2.*0.125)-0.0125);

    float x = min(p.x,min(p.z, p.y))/0.03125;

    return clamp(x, 0., 1.);


}

// Standard function-based bump mapping function (from Shane)
vec3 doBumpMap(in vec3 p, in vec3 nor, float bumpfactor){

    const float eps = BUMP_EPS;
    float ref = bumpSurf3D(p);
    vec3 grad = vec3( bumpSurf3D(vec3(p.x-eps, p.y, p.z))-ref,
                      bumpSurf3D(vec3(p.x, p.y-eps, p.z))-ref,
                      bumpSurf3D(vec3(p.x, p.y, p.z-eps))-ref )/eps;

    grad -= nor*dot(nor, grad);

    return normalize( nor + bumpfactor*grad );

}

vec2 raymarch(vec3 ro, vec3 rd, float steps) {

    float mat = 0.,
          t   = 0.,
          d   = 0.;
    vec3 p = ro;
    for(float i=.0; i<steps; i++) {

        d = map(p);
        mat = g_mat;  // save global material

        if(abs(d) < 0.0001 || t > 100.) break;

        t += d;
        p += rd*d;
    }

    return vec2(t, mat);
}

float n21(vec2 p) {
      return fract(sin(dot(p, vec2(524.423,123.34)))*3228324.345);
}

float noise(vec2 n) {
    const vec2 d = vec2(0., 1.0);
    vec2 b = floor(n);
    vec2 f = mix(vec2(0.0), vec2(1.0), fract(n));
    return mix(mix(n21(b), n21(b + d.yx), f.x), mix(n21(b + d.xy), n21(b + d.yy), f.x), f.y);
}

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

    vec3 tot = vec3(0);
    #if AA>0

    for(int m=0; m<AA; m++)
    for(int n=0; n<AA; n++) {
        vec2 off = vec2(float(m), float(n));
  	    vec2 uv = (fragCoord+off - .5*iResolution.xy)/iResolution.y;
    #else
   	vec2 uv = (fragCoord - .5*iResolution.xy)/iResolution.y;
    #endif

    tt = iTime + 17.;
    vec3 ro = vec3(uv*6.,-4.),
          rd = vec3(0,0,.6),
          lp = vec3(2., 2., -10),
          lp2 = vec3(-2., -3., -5);

    vec3 col;

    float mat = 0.,
          t   = 0.,
          d   = 0.;


    vec2 e = vec2(0.0035, -0.0035);

    // background color
    vec3 c1 = vec3(0.000,1.000,0.835);
    vec3 c2 = vec3(1.000,0.494,0.078);

    // light color
    vec3 lc1 = vec3(0.961,0.816,0.737);
    vec3 lc2 = vec3(0.588,0.992,0.945);


    // currently only one pass
    for(float i = 0.; i < 1.; i++) {
        float steps = i > 0. ? 50. : 250.;
        vec2 rm = raymarch(ro, rd, steps);
        mat = rm.y;

        vec3 p = ro + rm.x*rd;

        vec3 n = normalize( e.xyy*map(p+e.xyy) + e.yyx*map(p+e.yyx) +
                                e.yxy*map(p+e.yxy) + e.xxx*map(p+e.xxx));

        vec3 pt = transform(p);
        vec3 np = transform(n);
        n = doBumpMap(pt, n, .004);


        if(rm.x < 50.) {

            vec3 l = normalize(lp-p);
            vec3 l2 = normalize(lp2-p);
            float dif = max(dot(n, l), .0);
            float dif2 = max(dot(n, l), .0);
            float spe = pow(max(dot(reflect(-rd, n), -l), .0),40.);

            float sss = smoothstep(0., 1., map(p + l * .4)) / .4;
            vec3 n2 = n;
            n2.xy += noise(p.xy) * .5 - .025;
            n2 = normalize(n2);
            float height = atan(n2.y, n2.x);

            vec3 iri = spectral_zucconi6(height*1.11)*smoothstep(.8, .2, abs(n2.z));

            col += dif*lc1*.6 + .8*dif2*lc2 + .8*iri + .1;

            if(mat == 0.) {
                n += .5*texture(iChannel1, n.xy*100.).rgb;
                rd = reflect(rd, n);

                vec3 refl = texture(iChannel0, rd).rgb;

                refl *= mix(vec3(1), spectral_zucconi6(n.x*n.y*3.), .6); // reflect rainbows too
                col = mix(col, refl, .7);

                float height = atan(n.y, n.x);

            }

        } else {
            col =  mix(c1-.5, c2, (.9-pow(dot(uv, uv), .3)))*.2+.7; // background

        }

    }


    col *= mix(.2, 1., (1.5-pow(dot(uv, uv), .5))); // vignette
    col = pow(col, vec3(.7)); // gamma

    tot += col;
    #if AA>0
    }
    tot /= float(AA*AA);
    #endif

    fragColor = vec4(tot, 1.0);
}