const float epsilon = 0.001;

// Complex multiply
vec2 cmul(vec2 a, vec2 b) {
    return vec2(a.x * b.x - a.y * b.y, 2. * a.x * b.y);
}

// Advances j through one iteration of Julia fractal escape calculation.
vec2 julia_step(vec2 j) {
    vec2 ret = cmul(j, j);
    ret += vec2(sin(iTime*0.8) - 1.4, sin(iTime*0.73) * 1.8);
    return ret;
}

void mainImage(out vec4 fragColor, in vec2 fragCoord)
{
    vec2 uv = fragCoord / iResolution.xy;
    uv -= 0.5;
    uv /= vec2(iResolution.y/iResolution.x, 1.);
    uv *= 2.5;

    vec2 uv2 = uv;
    vec2 d = uv + vec2(epsilon);
    float traveled = 0.;
    float pTraveled = 0.;
    float traveled2 = 0.;
    float pTraveled2 = 0.;
    float uvLength = length(uv);

    int i = 0;
    while ( i < 100 && length(uv) < 4.) {
        i++;
        vec2 t = julia_step(uv);

        pTraveled2 = traveled2;
        pTraveled = traveled;
        vec2 s = t - uv2;
        float sLength = length(s);
        float inter = abs(sLength - uvLength);
        traveled += abs((length(t) - inter)/(uvLength - inter));
        traveled2 += abs((length(t) - inter)/(sLength - inter));
        uv = t;
    }

    // Reduces some of the terracing seen, especially in the blue color channel
    float overage = abs(length(uv) / 40.);
    if (overage < epsilon) {
        overage = epsilon;
    }

    for (int y = 0; y < 15; y++) {
        uv2 = julia_step(uv2);
        d = julia_step(d);
    }

    float dist = length(uv2 - d);
    float ldist = log(log(dist));

    float tia = traveled / float(i);
    float pTia = pTraveled / float(i - 1);
    float tia2 = traveled2 / float(i);
    float pTia2 = pTraveled2 / float(i - 1);

    fragColor = vec4(
        1. - 1./ldist,
        1. - (1. / tia + (pTia - tia) * fract(overage)),
        1. - (1. / tia2 + (pTia2 - tia2) * fract(overage)),
        1);
}