#version 330 core

out vec4 FragColor;
in vec2 vUV;

uniform vec3 iResolution;
uniform float iTime;
uniform float uAlpha;
uniform sampler2D iChannel0;

#define ITR 120
#define FAR 400.
#define PRECISION 0.00012
#define BASE_STRIDE .28
#define FAR_STRIDE 0.3
#define MIN_STEP .001

const float logvl = 1. + MIN_STEP;
#define time iTime
#define pi 3.1415926535

mat2 mm2(in float a)
{
    float c = cos(a), s = sin(a);
    return mat2(c, -s, s, c);
}

float sdSphere(vec3 p, float s)
{
    return length(p) - s;
}

float sdBox(vec3 p, vec3 b)
{
    vec3 d = abs(p) - b;
    return min(max(d.x, max(d.y, d.z)), 0.0) + length(max(d, 0.0));
}

float sdCylinder(vec3 p, vec2 h)
{
    vec2 d = abs(vec2(length(p.xz), p.y)) - h;
    return min(max(d.x, d.y), 0.0) + length(max(d, 0.0));
}

float opS(float d1, float d2)
{
    return max(-d2, d1);
}

vec2 opU(vec2 d1, vec2 d2)
{
    return (d1.x < d2.x) ? d1 : d2;
}

vec3 opRep(vec3 p, vec3 c)
{
    vec3 q = p;
    if (abs(c.x) > 1e-6) q.x = mod(p.x, c.x) - 0.5 * c.x;
    if (abs(c.y) > 1e-6) q.y = mod(p.y, c.y) - 0.5 * c.y;
    if (abs(c.z) > 1e-6) q.z = mod(p.z, c.z) - 0.5 * c.z;
    return q;
}

float hash31(vec3 p)
{
    return fract(sin(dot(p, vec3(127.1, 311.7, 74.7))) * 43758.5453);
}

vec3 hash33(vec3 p)
{
    float n = sin(dot(p, vec3(7, 157, 113)));
    return fract(vec3(2097152, 262144, 32768) * n) * 2. - 1.;
}

#define ANG_TO_WORLD 0.1
#define DNA_RADIUS 10.72 * ANG_TO_WORLD
#define RISE 3.32 * ANG_TO_WORLD
#define OFFSET 21.44 * ANG_TO_WORLD
#define ROTATION_PER_BP -0.180315633 / ANG_TO_WORLD
#define BP_APOTHEM 3.84170 * ANG_TO_WORLD
#define BP_WIDTH 20.01596 * ANG_TO_WORLD / 2.0

vec3 opDNATwist(vec3 p, float twistFactor)
{
    float c = cos(ROTATION_PER_BP * p.y * twistFactor);
    float s = sin(ROTATION_PER_BP * p.y * twistFactor);
    mat2 m = mat2(c, -s, s, c);
    return vec3(m * p.xz, p.y);
}

float mapDNA(vec3 pos, float ySplit, float twistFactor, float basePairY)
{
    float splitWidth = 0.2;
    float splitOffset = exp(-pow(splitWidth * (pos.y - ySplit), 2.0)) * .9;
    float direction = atan(pos.z, pos.x);

    float dnaBasePairs1 = sdBox(
        opRep(
            opDNATwist(pos - vec3(0, -RISE / 2.0, 0), twistFactor)
            + vec3(BP_APOTHEM, BP_WIDTH / 2., basePairY)
            + vec3(0, splitOffset, 0),
            vec3(0, 0, RISE)
        ),
        vec3(.07, BP_WIDTH / 2., .1)
    );

    float dnaBasePairs2 = sdBox(
        opRep(
            opDNATwist(pos - vec3(0, -RISE / 2.0, 0), twistFactor)
            + vec3(BP_APOTHEM, -BP_WIDTH / 2., basePairY)
            - vec3(0, splitOffset, 0),
            vec3(0, 0, RISE)
        ),
        vec3(.07, BP_WIDTH / 2., .1)
    );

    vec2 strands = opU(
        vec2(
            sdBox(
                opDNATwist(pos - vec3(cos(direction) * splitOffset, 0.0, sin(direction) * splitOffset), twistFactor) - vec3(0, DNA_RADIUS, 0),
                vec3(.225, .125, 1e6)
            ),
            1e6
        ),
        vec2(
            sdBox(
                opDNATwist(pos + vec3(0, OFFSET, 0) - vec3(cos(direction) * splitOffset, 0.0, sin(direction) * splitOffset), twistFactor) - vec3(0, DNA_RADIUS, 0.0),
                vec3(.225, .125, 1e6)
            ),
            1e6
        )
    );

    vec2 completeBP = opU(
        vec2(dnaBasePairs1, 40.0),
        vec2(dnaBasePairs2, 40.0)
    );

    return opU(strands, completeBP).x;
}

#define REPEAT vec3(6, 0, 6)

float mapProtein(vec3 pos, float yPos)
{
    vec3 add = vec3(0, yPos, 0);
    float proteins = opS(
        sdSphere(opRep(pos, REPEAT) + add, DNA_RADIUS * 2.2),
        sdCylinder(opRep(pos, REPEAT) + add, vec2(DNA_RADIUS * 2., DNA_RADIUS * 2.4))
    );
    return proteins;
}

float hash(vec2 p) { return fract(sin(1.0 + dot(p, vec2(127.1, 311.7))) * 43758.545); }
vec2 sincos(float x) { return vec2(sin(x), cos(x)); }

float map(in vec3 pos)
{
    vec2 id = floor((pos.xz) / REPEAT.x);
    float ph = hash(id + 113.1);
    float ve = hash(id);

    vec3 p1 = pos;
    p1.xz = mod(pos.xz, REPEAT.x) - REPEAT.x / 2.;
    p1.xz += 0.5 * cos(2.0 * ve * iTime + (pos.y + ph) * vec2(0.53, 0.32) - vec2(1.57, 0.0));

    float direction = -(mod(floor((p1.x - REPEAT.x / 2.) / REPEAT.x), 2.) * 2.0 - 1.0);
    float proteinPos = cos((id.x / 5. + id.y / 10.) * time) * 30. + time * 2.0;

    float res = mapDNA(
        id.x == 0. && id.y == 0. ? pos - REPEAT.x / 2. : p1,
        id.x == 0. && id.y == 0. ? time * 2.0 : proteinPos,
        1.,
        direction * time * 0.5
    );

    float proteins = mapProtein(
        p1 + REPEAT / 2.,
        id.x == 0. && id.y == 0. ? time * 2.0 : -proteinPos
    );

    return opU(vec2(res, 1.0), vec2(proteins, 1.0)).x;
}

float shadow(in vec3 ro, in vec3 rd, float mint, float maxt)
{
    for (float t = mint; t < maxt;)
    {
        float h = map(ro + rd * t);
        if (h < 0.001) return 0.0;
        t += h;
    }
    return 1.0;
}

float bisect(in vec3 ro, in vec3 rd, in float near, in float far)
{
    float mid = 0.;
    float sgn = sign(map(rd * near + ro));
    for (int i = 0; i < 6; i++)
    {
        mid = (near + far) * .5;
        float d = map(rd * mid + ro);
        if (abs(d) < PRECISION) break;
        d * sgn < 0. ? far = mid : near = mid;
    }
    return (near + far) * .5;
}

float intersect(in vec3 ro, in vec3 rd)
{
    float t = 0.;
    float d = map(rd * t + ro);
    float sgn = sign(d);
    float told = t;
    bool doBisect = false;

    for (int i = 0; i <= ITR; i++)
    {
        if (abs(d) < PRECISION || t > FAR) break;
        else if (i == ITR) t = 1000.;

        if (sign(d) != sgn)
        {
            doBisect = true;
            break;
        }

        told = t;
        if (d > 1.) t += d * FAR_STRIDE;
        else t += log(abs(d) + logvl) * BASE_STRIDE;
        d = map(rd * t + ro);
    }

    if (doBisect) t = bisect(ro, rd, told, t);
    return t;
}

float getGray(vec3 p)
{
    return p.x * 0.299 + p.y * 0.587 + p.z * 0.114;
}

mat2 rot(float th)
{
    float cs = cos(th), sn = sin(th);
    return mat2(cs, -sn, sn, cs);
}

vec3 firePalette(float i)
{
    float T = 1600. + 1300. * i;
    vec3 L = vec3(
        5.0,
        1.0,
        (texture(iChannel0, vec2(.1, .2)).x + 0.01) * 10.
    );
    L = pow(L, vec3(5.0)) * (exp(1.43876719683e5 / (T * L)) - 1.0);
    return 1.0 - exp(-5e8 / L);
}

vec3 normal(const in vec3 p)
{
    vec2 e = vec2(-1., 1.) * 0.01;
    return normalize(
        e.yxx * map(p + e.yxx) +
        e.xxy * map(p + e.xxy) +
        e.xyx * map(p + e.xyx) +
        e.yyy * map(p + e.yyy)
    );
}

float trig3(in vec3 p)
{
    p = cos(p * 2. + (cos(p.yzx) + 1. + time * 4.) * 1.57);
    return dot(p, vec3(0.1666)) + 0.5;
}

float trigNoise3D(in vec3 p)
{
    const mat3 m3RotTheta = mat3(0.25, -0.866, 0.433, 0.9665, 0.25, -0.2455127, -0.058, 0.433, 0.899519) * 1.5;
    float res = 0.;

    float t = trig3(p * pi);
    p += (t - time * 0.25);
    p = m3RotTheta * p;
    res += t;

    t = trig3(p * pi);
    p += (t - time * 0.25) * 0.7071;
    p = m3RotTheta * p;
    res += t * 0.7071;

    t = trig3(p * pi);
    res += t * 0.5;

    return res / 2.2071;
}

float getMist(in vec3 ro, in vec3 rd, in vec3 lp, in float t)
{
    float mist = 0.;
    ro += rd * t / 8.;

    for (int i = 0; i < 4; i++)
    {
        float sDi = length(lp - ro) / FAR;
        float sAtt = min(1. / (1. + sDi * 0.25 + sDi * sDi * 0.05), 1.);
        mist += trigNoise3D(ro * 2.) * sAtt;
        ro += rd * t / 4.;
    }

    return clamp(mist / 2. + hash31(ro) * 0.1 - 0.05, 0., 1.);
}

float softShadow(in vec3 ro, in vec3 rd, in float start, in float end, in float k)
{
    float shade = 1.0;
    const int maxIterationsShad = 12;
    float dist = start;

    for (int i = 0; i < maxIterationsShad; i++)
    {
        float h = map(ro + rd * dist);
        shade = min(shade, k * h / dist);
        dist += clamp(h, 0.0001, .2);
        if (abs(h) < 0.001 || dist > end) break;
    }

    return min(max(shade, 0.) + 0.4, 1.0);
}

float calculateAO(vec3 p, vec3 n)
{
    const float AO_SAMPLES = 5.0;
    float r = 0.0, w = 1.0, d;

    for (float i = 1.0; i < AO_SAMPLES + 1.1; i++)
    {
        d = i / AO_SAMPLES;
        r += w * (d - map(p + n * d));
        w *= 0.5;
    }

    return 1.0 - clamp(r, 0.0, 1.0);
}

void mainImage(out vec4 fragColor, in vec2 fragCoord)
{
    vec2 uv = (fragCoord.xy - iResolution.xy * 0.5) / iResolution.y;
    vec3 ro = vec3(REPEAT.x / 2. + 0.00001, time * 2., REPEAT.z / 2.);
    vec3 rd = normalize(vec3(uv, 0.5));

    mat2 m2 = rot(time * 0.15);
    mat2 m2_2 = rot(time * 0.055);
    mat2 m2_3 = rot(time * 0.075);
    rd.xz *= m2;
    rd.xy *= m2_2;
    rd.yz *= m2_3;

    vec3 lp = vec3(0., 0., FAR) + ro + rd * 10.;
    float bgShade = pow(max(dot(rd, normalize(lp - ro)), 0.) * 0.5 + 0.5, 4.);
    vec3 bc = mix(vec3(.0, .0, .05), vec3(1.), bgShade);
    vec3 sc = bc;

    float t = intersect(ro, rd);

    if (t < FAR)
    {
        vec3 sp = ro + rd * t;
        vec3 sn = normal(sp);
        vec3 objCol = vec3(1.0, 0.5, 0.7);
        vec3 ld = lp - sp;
        float lDist = max(length(ld), 0.001);
        ld /= lDist;
        lDist /= FAR;
        float sAtten = min(1. / (1. + lDist * 0.125 + lDist * lDist * 0.05), 1.);
        float shad = softShadow(sp, ld, 0.05, FAR, 8.);
        float ao = calculateAO(sp, sn);
        float diff = max(dot(sn, ld), 0.);
        float spec = pow(max(dot(reflect(-ld, sn), -rd), 0.0), 8.);
        sc = (objCol * (diff + 0.5) + spec) * sAtten;
        sc = min(sc, 1.) * shad * ao;
    }

    float fog = min(1.0 / (1. + (t * 0.25) + (t * t * 0.025)), 1.);
    sc = mix(bc, sc, fog);

    vec3 sc2 = firePalette(getGray(sc));
    float fadeFactor = min(1.0 / (1. + t), 1.);
    sc = mix(sc, sc2, fadeFactor * 0.34 + 0.66);

    float mist = getMist(ro, rd, lp, t);
    vec3 fogCol = mix(vec3(bgShade * 0.8 + 0.2) * vec3(1., 0.85, 0.6), sc, mist * fog);
    sc = sc * 0.65 + fogCol * 0.35;

    fragColor = vec4(clamp(sc, 0., 1.), 1.0);
}

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