#version 330 core

uniform vec3 iResolution;
uniform float iTime;
uniform float uAlpha;
uniform sampler2D iChannel0;
uniform float uFireBaseRadius;
uniform float uFireAmplitude;
uniform float uFireColorRed;
uniform float uFireColorGreen;
uniform float uFireColorBlue;
uniform float uMainIterations;
uniform float uFrequencyIterations;
uniform float uMandelbrotIterations;
uniform float uFrequencyInitial;
uniform float uRotationScale;
uniform float uTimeScale;
uniform float uFrequencyScale;
uniform float uDepthScale;

out vec4 FragColor;

// stable tanh for vec4
vec4 tanh4(vec4 x) {
    vec4 e = exp(clamp(2.0 * x, vec4(-20.0), vec4(20.0)));
    return (e - vec4(1.0)) / (e + vec4(1.0));
}

vec4 fire(vec3 p) {
    float s, i, d, n, T = iTime;
    vec4 o = vec4(0.0);
    p += cos(p.z + T + p.yzx * 0.5) * uFireAmplitude;
    s = uFireBaseRadius - length(p.xy);
    // rotate XY by time-based matrix using cos; original used mat2(cos(T+vec4(0,33,11,0))) which is odd; approximate
    p.xy = mat2(cos(T), -sin(T), sin(T), cos(T)) * p.xy;
    for (n = uFrequencyInitial; n < uFrequencyIterations; n += n)
        s -= abs(dot(sin(p.z + uFrequencyScale * T + p * n), vec3(0.5))) / n;
    s = (0.01 + abs(s) * 0.15);
    o += 2.0 / s;
    return (vec4(uFireColorRed, uFireColorGreen, uFireColorBlue, 1.0) * o * o);
}

void mainImage(out vec4 o, vec2 u) {
    float d = 0.5 * texelFetch(iChannel0, ivec2(u) % 1024, 0).a;
    float i = 0.0;
    float s = 0.0;
    float w = 0.0;
    float l = 0.0;
    float T = iTime;
    vec3 q, p = iResolution;

    u = (u - p.xy / 2.0) / p.y;
    for(o = vec4(0.0); i++ < uMainIterations; o += d / s + uTimeScale * fire(q)) {
        q = p = vec3(u * d, d + T * uTimeScale);
        // alternate rotation
        float cosR = cos(uRotationScale * T + p.z * uDepthScale);
        float sinR = sin(uRotationScale * T + p.z * uDepthScale);
        p.xy = mat2(cosR, sinR, -sinR, cosR) * p.xy;
        p *= 0.5;
        p.xy += 1.5;
        w = 0.75;
        for (int ii = 0; ii++ < int(uMandelbrotIterations); ) {
            p *= l = 3.0 / dot(p = sin(p), p);
            w *= l;
        }
        d += s = 0.002 + abs(length(p) / w);
    }
    o = mix(o, o.zyxw, smoothstep(0.2, 1.0, dot(u, u)));
    // avoid division by zero
    float lenu = max(length(u), 1e-6);
    float denom = max(d * lenu, 1e-6);
    o = tanh4(max(o / 200000.0 / denom, vec4(0.0)));
}

void main() {
    mainImage(FragColor, gl_FragCoord.xy);
    FragColor.a *= uAlpha;
}
