#version 300 es
precision highp float;

uniform float iTime;
uniform vec2  iResolution;

uniform float ShadeMode; // value=1, min=0, max=1, step=1 (No, Yes)
uniform float ImageMode; // value=1, min=0, max=1, step=1 (No, Yes)
uniform float ImageXScale; // value=1, min=0, max=2, step=0.1
uniform float ImageYScale; // value=1, min=0, max=2, step=0.1
uniform float Front; // value=1, min=0, max=1, step=1 (No, Yes)
uniform float Animate; // value=1, min=0, max=1, step=1 (No, Yes)
uniform float Radius; // value=0.5, min=0, max=1, step=0.01
uniform float Normalised; // value=1, min=0, max=1, step=1 (No, Yes)
uniform float Julia; // value=1, min=0, max=1, step=1 (No, Yes)
uniform float Julia_c1; // value=0, min=-1.0, max=1.0, step=0.01
uniform float Julia_c2; // value=0, min=-1.0, max=1.0, step=0.01
uniform int MaxIter; // value=50, min=1, max=1000, step=1
uniform int AA; // value=4, min=1, max=9, step=1
uniform float Zoom; // value=1, min=0.01, max=9, step=0.01
uniform float Xoffset; // value=0, min=-1, max=1, step=0.01
uniform float Xpan; // value=0, min=-1, max=1, step=0.01
uniform float Yoffset; // value=0, min=-1, max=1, step=0.01
uniform float Ypan; // value=0, min=-1, max=1, step=0.01
uniform float light; // value=1.5, min=0, max=4, step=0.01
uniform float angle; // value=45, min=0, max=180, step=0.1
uniform float EscapeRadius; // value=4, min=0, max=100, step=0.01

float jc1;
float jc2;

out vec4 fragColor;

vec2 cMul(vec2 a, vec2 b) {
    return vec2( a.x*b.x -  a.y*b.y,a.x*b.y + a.y * b.x);
}

vec2 cInverse(vec2 a) {
    return  vec2(a.x,-a.y)/dot(a,a);
}
vec2 cDiv(vec2 a, vec2 b) {
    return cMul( a,cInverse(b));
}
vec2 cExp(in vec2 z){
    return vec2(exp(z.x)*cos(z.y),exp(z.x)*sin(z.y));
}

// https://thebookofshaders.com/10/
// Random function
float random (vec2 st) {
    return fract(sin(dot(st.xy,
                         vec2(12.9898,78.233)))*
        43758.5453123);
}

// Colour scheme
// Stolen from llemarie
// vec3 HSVtoRGB(float iter){
//     iter = mod(iter+iTime, float(MaxIter));
//     float value = 1. - float(iter)/float(MaxIter);
//     value *= value;
        
//     float h = value;
//     float s = 1.;
//     float v = (iter < float(MaxIter)) ? 1. - value*value : 0.;
    
//     float r, g, b, f, p, q, t;
//     int i;
//     i = int(h * 6.);
//     f = h * 6. - float(i);
//     p = v * (1. - s);
//     q = v * (1. - f * s);
//     t = v * (1. - (1. - f) * s);
//     switch (i % 6){
//         case 0: r = v, g = t, b = p; break;
//         case 1: r = q, g = v, b = p; break;
//         case 2: r = p, g = v, b = t; break;
//         case 3: r = p, g = q, b = v; break;
//         case 4: r = t, g = p, b = v; break;
//         case 5: r = v, g = p, b = q; break;
//     }
//     return vec3(r, g, b);
// }

vec3 palette(float t)
{
    //return vec3(0.5,0.5,0.5) + vec3(0.5,0.5,0.5)*cos(6.28318*(vec3(2.0,1.0,0.0)*t + vec3(0.5,0.2,0.25)));
    return vec3(0.5) + vec3(0.5)*cos(6.28318*(vec3(1.0)*t + vec3(0.,0.33,0.67) + iTime));
    //return vec3(0.5, 0.5, 0.5) + vec3(0.5, 0.5, 0.5)*cos(6.28318*(vec3(1.0, 1.0, 1.0)*t+2.*iTime*vec3(0.01, 0.10, 0.20)) );
}


// Put anything you like here. This is the image that is rendered in the image orbit traps
vec4 subrender(float xin, float yin){
    xin = xin/ImageXScale;
    yin = yin/ImageYScale;
    float sqr = sqrt((xin-0.5)*(xin-0.5)+(yin-0.5)*(yin-0.5));
    if(sqr<Radius){
        return vec4(sin(xin + iTime)*2. + 0.5,yin,sqr,1.0);
    }
    else{
        return vec4(0.0);
    }
}

vec4 render(float xin, float yin){
    vec4 final_col = vec4(0.0);
    for(int i = 0; i<AA;i++){
        float x = 0.0;
        float y = 0.0;
        float xx = 0.0;
        float yy = 0.0;
        float xy = 0.0;
        
        // sample points around the main pixel for anti-aliasing
        float r = random(vec2(xin + float(i),yin + float(i)))*0.0001;
        float c1 = (1./Zoom)*(xin + Xoffset + r + Xpan);
        float c2 = (1./Zoom)*(yin + Yoffset + r - Ypan);
        
        // Julia mode
        if(Julia>0.0){
            x = c1;
            y = c2;
            c1 = jc1;
            c2 = jc2;
        }
        
        vec2 dc = vec2(1, 0);
        vec2 der = dc;
        
        float h2 = light;
        float angle = angle;
        vec2 v = cExp(vec2(0.0, 45.*2.*3.1415/360.));
        v = vec2(sin(iTime),cos(iTime));
        
        int j = 0;
        // Flag keeps track of opqque pixels
        int flag = 0;
        
        vec4 col = vec4(palette(1.0), 1.0);
        
        while(j<MaxIter){
            xx = x*x;
            yy = y*y;
            if((xx + yy) > EscapeRadius){
                if(flag<1){
                    // Normalised/Smooth iteration count
                    if(Normalised>0.0){
                        float iter =float(j) - log2(log2(dot(vec2(x,y),vec2(x,y)))) + 4.0;
                        
                        // float logzn = log(xx+yy);
                        // float log4 = log(4.0);
                        // float log2 = 1.44269504089;
                        // float logzndiv4 = logzn/log4;
                        // float loglogzndiv4 = log(logzndiv4);
                        // float nu = loglogzndiv4 * 1.0/log2;
                        // nu = log(logzn/log(4.0))/log(2.0);
                        // float iter = float(j) + 1.0 - nu;
                        
                        col = vec4(palette(iter), 1.0);
                    }
                    else{
    
                        
                        col = vec4(palette(float(j)), 1.0);
                    }
                }
                
                break;
            }
            der = cMul(vec2(x,y), der) * 2.0 + dc;
            
            // TODO: Karatsuba style - not sure how much time it saves. https://mathr.co.uk/blog/2016-04-10_complex_squaring.html
            xy = x*y;
            x = xx - yy + c1;
            y = xy+xy + c2;
            
            // Image orbit trap
            if(ImageMode>0.0 && x>-1.0*ImageXScale && x<=1.0*ImageXScale && y > -1.0*ImageYScale && y<=1.0*ImageYScale){
                vec4 sr = subrender(x+0.5, y+0.5);
                if(sr.a>0.0){
                    col = sr;
                    if(Front>0.0){
                        break;
                    }
                    flag = 1;
                }
            }
            j++;
        }
        
        if(ShadeMode>0.0){
        
            vec2 u = cDiv(vec2(x,y), der);
            u = cDiv(u, abs(vec2(u.x, u.y)));
            float t = dot(u,v) + h2;
            t = t/(1. + h2);
            
            if(t<0.){
                t =0.;
            }
            
            col.xyz *= t;
        }
        
        final_col += col;
        
        

    }
    
    return final_col/float(AA);
}

void main() {
    // Screen coordinate (from [-aspect, -1] to [aspect, 1])
    vec2 q   = (2. * gl_FragCoord.xy - iResolution) / iResolution.y;
    
    
    jc1 = Julia_c1;
    jc2 = Julia_c2;
    if(Animate>0.0){
        jc1 = sin(Julia_c1 + iTime + 5.);
        jc2 = cos(Julia_c2 + iTime + 7.);
    }
    
    vec4 col = render(q.x, q.y);
    col.xyz = pow(col.xyz, vec3(0.9));
    
    // Output
    fragColor = col;
}