#ifdef GL_ES
precision mediump float;
#endif

uniform float time;
uniform vec2 mouse;
uniform vec2 resolution;

#define kPI05    1.5707963267949
#define kPI      3.1415926535898
#define kPI2     6.2831853071796

//
vec3 rotZ(vec3 p, float a) {
    vec2 sc = sin(vec2(a, a + kPI05));
    return vec3(
        p.x * sc.y + p.y * -sc.x,
        p.x * sc.x + p.y * sc.y,
        p.z
    );
}

vec3 lookAt(vec3 o, vec3 t, vec3 u, vec3 d) {
    vec3 f = normalize(o - t);
    vec3 s = normalize(cross(u, f));
    vec3 h = normalize(cross(f, s));
    return vec3(
        dot(vec3(s.x, h.x, f.x), d),
        dot(vec3(s.y, h.y, f.y), d),
        dot(vec3(s.z, h.z, f.z), d));
}

// Hash without Sine https://www.shadertoy.com/view/4djSRW
float hash11(float p) {
    p = fract(p * .1031);
    p *= p + 33.33;
    p *= p + p;
    return fract(p);
}
float hash13(vec3 p3) {
    p3  = fract(p3 * .1031);
    p3 += dot(p3, p3.yzx + 33.33);
    return fract((p3.x + p3.y) * p3.z);
}
vec2 hash22(vec2 p) {
    vec3 p3 = fract(vec3(p.xyx) * vec3(.1031, .1030, .0973));
    p3 += dot(p3, p3.yzx+33.33);
    return fract((p3.xx+p3.yz)*p3.zy);
}
vec3 hash33(vec3 p3) {
    p3 = fract(p3 * vec3(.1031, .1030, .0973));
    p3 += dot(p3, p3.yxz+33.33);
    return fract((p3.xxy + p3.yxx)*p3.zyx);
}

//
vec3 voro3(vec3 p, vec3 r) {
    vec3 ret;
    vec3 ip = floor(p);
    const int N = 1;
    float lmax = 10.0;
    vec3 minip;
    vec3 mintp;
    
    for(int iz = -N; iz <= N; iz++) {
        for(int iy = -N; iy <= N; iy++) {
            for(int ix = -N; ix <= N; ix++) {
                vec3 tp = ip + vec3(ivec3(ix, iy, iz));
                tp += (hash33(tp) - vec3(0.5)) * r;
                float l = length(tp - p);
                if(l < lmax) {
                    lmax = l;
                    minip = ip + vec3(ivec3(ix, iy, iz));
                    mintp = tp;
                }
            }
        }
    }
    ret.x = lmax;
    ret.z = hash13(minip);
    
    lmax = 10.0;
    vec3 minvp;
    for(int iz = -N; iz <= N; iz++) {
        for(int iy = -N; iy <= N; iy++) {
            for(int ix = -N; ix <= N; ix++) {
                vec3 tp = minip + vec3(ivec3(ix, iy, iz));
                tp += (hash33(tp) - vec3(0.5)) * r;
                vec3 vp = tp - mintp;
                float l = length(vp);
                if(l > 0.0) {
                    vec3 mp = (tp + mintp) * 0.5;
                    l = abs(dot(p - mp, vp / l));
                    if(l < lmax) {
                        lmax = l;
                        minvp = vp;
                    }
                }
            }
        }
    }
    ret.y = lmax;
    return ret;
}

float star(vec2 p) {
    float a = atan(p.y, p.x) / kPI2 + 0.5;
    a = abs(fract(a * 5.0) * 2.0 - 1.0) * kPI2 / 10.0;
    float r = length(p);
    vec2 q = cos(vec2(a, a - kPI05)) * r;
    float s = kPI / 180.0 * 55.0;
    vec2 n = cos(vec2(s, s - kPI05));
    return 0.5 - dot(q, n);
}

float flower(vec2 p, float w) {
    vec2 tt = vec2(1e2);
    for(float i = 0.0; i < 5.0; i+=1.0) {
        float a = kPI2 * i / 5.0;
        vec2 o = cos(vec2(a, a - kPI05));
        float to = length(p - o);
        tt = min(tt, vec2(to, abs(to - 1.1756)));
    }
    return w - (length(p) > 1.0 ? tt.x : tt.y);
}

//
float smin(float a, float b, float k) {
    k *= 4.0;
    float h = max(k - abs(a - b), 0.0) / k;
    return min(a, b) - h * h * k * 0.25;
}

float sphere(vec3 p, vec3 o, float r) {
    return length(p - o) - r;
}

vec2 tsphere(vec3 ro, vec3 rd, vec3 sc, float sr, out vec3 ret0, out vec3 ret1) {
    vec3 oc = ro - sc;
    float a = dot(rd, rd);
    float b = 2.0 * dot(rd, oc);
    float c = dot(oc, oc) - sr * sr;
    float d = b * b - 4.0 * a * c;
    if(d < 0.0) return vec2(-1.0);
    vec2 t = (vec2(-b) + vec2(-1.0, 1.0) * sqrt(d)) * 0.5 / a;
    ret0 = rd * t.x + ro;
    ret1 = rd * t.y + ro;
    return t;
}

float tplane(vec3 o, vec3 d, vec3 n, vec3 p, out vec3 ret) {
    vec3 po = p - o;
    float nd = dot(n, d);
    float t = dot(n, po) / nd;
    if(t < 0.0 || abs(nd) <= 1e-8) return -1.0;
    ret = d * t + o;
    return t;
}

// 
vec3 bg(vec3 p) {
    vec3 d = normalize(p);
    vec3 c = vec3(0.0, 0.0, 0.02);

    vec3 vr = voro3(d * 120.0, vec3(1.0));
    float f = 10.0 + abs(sin((vr.z * 24.0 + time) * 2.0)) * 20.0;
    vec3 s = mix(vec3(1.0, 0.6, 0.4), vec3(0.4, 0.9, 1.0), fract(vr.z * vr.z * 1e3));
    c += pow(vec3(1.0 - sqrt(vr.x)), vec3(10.0)) * f * s;

    vec4 hc = mix(vec4(0.1, 0.15, 0.05, 25.0),
                    vec4(0.0, 0.2, 0.8, 50.0),
                    smoothstep(-0.03, 0.03, d.y));
    c += pow(1.0-abs(d.y), hc.w) * hc.rgb;
    c += exp(d.y * -d.y * 3e4) * pow(-d.z, 12.0) * vec3(1.0);

    return c;
}

// scene
const float kDz = 4.0;
const float kIntro = 8.0;
const float kSeq = 16.0;

vec2 i2p(vec2 p, float i, inout float y) {
    if(i > kIntro) {
        float iseq = mod(i - kIntro, kSeq);
        if(iseq == 0.0) {
            p = vec2(0.0, -1.0);
            y *= 2.0;
        } else if(iseq > kSeq - 4.0) {
            p = vec2(p.x * 0.1, 2.0 - (kSeq - iseq));
        }
    }
    return p;
}

vec2 pos2(float t) {
    float i = floor(t);
    float f = fract(t);

    vec2 xx = vec2(-1.0, 1.0) * (fract(i * 0.5) * 4.0 - 1.0);
    float y = 4.0 * f * (1.0 - f);

    vec2 o0 = hash22(vec2(i) + vec2(0.0, 0.3456789));
    vec2 o1 = hash22(vec2(i + 1.0) + vec2(0.0, 0.3456789));

    o0 = o0 * 2.0 - vec2(1.0);
    o1 = o1 * 2.0 - vec2(1.0);

    vec2 p0 = vec2(xx.x, 0.0) + o0;
    vec2 p1 = vec2(xx.y, 0.0) + o1;

    p0 = i2p(p0, i, y);
    p1 = i2p(p1, i + 1.0, y);

    return mix(p0, p1, f) + vec2(0.0, y);
}

float map(vec3 p, vec3 tbif) {
    float tb = tbif.x;
    float ti = tbif.y;
    float tf = tbif.z;

    const float N = 16.0;
    const float r0 = 0.01;
    float d = sphere(p, vec3(pos2(tb), 0.0), r0);
    float dt = 0.0;
    for(float i = 1.0; i < N; i+=1.0) {
        float f = sqrt(i / (N - 1.0));
        f = 4.0 * f * (1.0 - f);
        dt += mix(0.02, 0.05, f);
        vec2 sp = pos2(tb - dt);
        float r = 0.1 * f + r0;
        float s = mix(0.03, 0.05, f);
        d = smin(d, sphere(p, vec3(sp, dt * kDz), r), s);
    }
    return d;
}

vec3 normal(vec3 p, vec3 t) {
    int id;
    vec2 e = vec2(1.0, -1.0) * 0.5773;
    const float eps = 0.00025;
    return normalize( e.xyy * map(p + e.xyy * eps, t) + 
					  e.yyx * map(p + e.yyx * eps, t) + 
					  e.yxy * map(p + e.yxy * eps, t) + 
					  e.xxx * map(p + e.xxx * eps, t) );
}

vec3 fxs(vec3 ro, vec3 rd, float maxt, vec3 tbif) {
    float tb = tbif.x;
    float ti = tbif.y;
    float tf = tbif.z;

    vec3 ret = vec3(0.0, 0.0, 0.0);

    for(float i = 0.0; i < 4.0; i+=1.0) {
        float pti = ti - i;
        if(pti < 0.0) break;

        float ft = (tb - pti);
        vec3 fp = vec3(pos2(pti), ft * kDz);
        if(pti < kIntro) maxt = 100.0;

        {
            vec3 hp;
            float ht = tplane(ro, rd, vec3(0.0, 1.0, 0.0), fp, hp);
            if(ht >= 0.0 && ht < maxt) {
                float r = ft * 1.2;
                float t = length(fp - hp);
                float f = smoothstep(0.95, 1.0, 1.0 - abs(t - r));
                f *= max(0.0, 1.0 - t * t * 1.0);
                ret += vec3(1.0, 1.0, 1.0) * f;
            }
        }

        if(pti < kIntro) continue;
        float scl = (mod(pti - kIntro, kSeq) == 0.0) ? 4.0 : 1.0;
        scl = (pti == kIntro) ? 1.0 : scl;

        for(float k = 0.0; k < 2.0; k+=1.0) {
            float fd = 1.0 - smoothstep(1.0, 1.5, ft);
            float sr = (0.5 * ft * ft + 1e-1) * (1.0 + k * 2.0) * scl;
            vec3 sp = fp + vec3(0.0, sr, 0.0);
            vec3 shp0, shp1;
            vec2 sht = tsphere(ro, rd, sp, sr, shp0, shp1);
            for(int j = 0; j < 2; j++) {
                float st = (j == 0) ? sht.x : sht.y;
                if(st < 0.0 || st > maxt) continue;
                vec3 lp = rd * st + ro - sp;

                float a = atan(lp.z, lp.x) + pti * 0.7 + k * 0.6;
                float r = acos(clamp(-lp.y / sr, -1.0, 1.0)) / kPI;

                r *= max(1.0, sr * 5.0 / scl);
                vec2 uv = cos(vec2(a, a - kPI05)) * r;
                // ret += vec3(fract(length(uv)));
                float tx = flower(uv, 0.04);
                ret += vec3(4.0, 3.0, 2.0) * smoothstep(-0.02, 0.04, tx) * fd;
            }
        }

        float strrng = fract(sin(pti + 1.7654321));
        for(float k = 0.0; k < 20.0; k+= 1.0) {
            strrng = fract(strrng * (1.23456 + strrng) * 100.0);
            float rnd = strrng;
            float a = (k + pti) * 2.4 + rnd;
            vec3 pp = fp;
            pp.xz += cos(vec2(a, a - kPI05)) * ft * (0.4 + rnd * 0.8) * scl;
            pp.y += ft * (0.9 + rnd * 0.5 - ft) * 2.0 * scl;

            vec3 hp;
            float ht = tplane(ro, rd, vec3(0.0, 0.0, -1.0), pp, hp);
            if(ht >= 0.0 && ht < maxt) {
                vec2 uv = (hp.xy - pp.xy) * 12.0 / scl;
                uv = rotZ(uv.xyy, k * 2.4).xy;
                float tx = star(uv);
                float fd = 1.0 - smoothstep(1.0, 1.5, ft);
                ret += vec3(1.5, 1.0, 0.6) * fd * smoothstep(-0.1, 0.1, tx);
            }
        }
    }

    return ret;
}

void main() {
    vec3 t;
    t.x = time * 120.0 / 60.0;
    t.y = floor(t.x);
    t.z = fract(t.x);

    vec3 color = vec3(0.0);
    vec2 sc = (gl_FragCoord.xy * 2.0 - resolution.xy) / min(resolution.x, resolution.y);

    vec3 ro = vec3(0.0, 0.0, 10.0);
    vec3 rd = normalize(vec3(sc.xy, -4.0));
    ro.xy += vec2(0.0, 0.4 + sin(time * 1.3) * 0.03) + sin(vec2(time) * vec2(0.7, 1.15)) * 0.04;
    vec3 lp = vec3(pos2(t.x - 0.6) * vec2(0.02, 0.002), sin(time * 0.2) * 0.05);
    rd = lookAt(ro, lp, vec3(0.0, 1.0, 0.0), rd);

    color = bg(rd);

    float td = 0.0;
    float blm = 0.0;
    for(int i = 0; i < 64; i++) {
        vec3 rp = rd * td + ro;
        float d = map(rp, t);
        if(d < 1e-2) {
            vec3 nv = normal(rp, t);
            color = mix(vec3(0.4, 0.4, 0.55), vec3(0.8), (nv.y * 0.5 + 0.5));
            color += vec3(0.2, 0.3, 0.0) * (1.0 - abs(nv.y));
            break;
        }
        blm += pow(max(0.0, 1.0 - d * 8.0), 10.0); 
        td += d;
    }
    color += vec3(0.1, 0.2, 0.2) * blm;

    if(t.y < kIntro) {
        color = vec3(0.0);
    } else {
        color += max(0.0, 1.0 - (t.x - kIntro) * 2.0);
    }

    color += fxs(ro, rd, td, t);
    gl_FragColor = vec4(pow(color, vec3(1.0/2.2)), 1.0);
}