precision highp float;
uniform vec2 resolution;
uniform vec2 mouse;
uniform float time;
uniform sampler2D backbuffer;
// Created by evilryu
// License Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.

#define FAR 100.

vec3 sun_dir = normalize(vec3(0., 1., 100.));
//vec3 blackHolePos = vec3(30, 40, 180.);
vec3 blackHolePos = vec3(-10, 40, 160.);
float blackHoleRadius = 22.;
float blackHoleDiscRadius = 120.;
//float blackHoleDiscRadius = 220.;
float blackHoleTilt = .2;
float myTime,fade = 0.;
float sandSpeed = 0.;
float fogDensity = 0.;
float renderBH = 1.;
float hackFlicker = 0.;


#define MAT_BUILDING 0.
#define MAT_WINDOW 1.
#define MAT_SAND 2.
#define MAT_DISC 3.
#define MAT_BACKGROUND -1.


#define PI 3.1415926535

float hash12(vec2 p)
{
    p = fract(p * vec2(5.3983, 5.4427));
    p += dot(p.yx, p.xy + vec2(21.5351, 14.3137));
    return fract(p.x * p.y * 95.4337);
}

vec2 hash22(vec2 p)
{
    p = fract(p * vec2(5.3983, 5.4427));
    p += dot(p.yx, p.xy +  vec2(21.5351, 14.3137));
    return fract(vec2(p.x * p.y * 95.4337, p.x * p.y * 97.597));
}

float hash13(vec3 p)
{
    p = fract(p * vec3(5.3983, 5.4472, 6.9371));
    p += dot(p.yzx, p.xyz + vec3(21.5351, 14.3137, 15.3219));
    return fract(p.x * p.y * p.z * 95.4337);
}


// Robust hash without sine
vec4 hash4(vec4 p)
{
    p = fract(p * .1031);
    p *= p + 33.33;
    p *= p + p;
    return fract(p);
}

float noise(vec3 p)
{
    vec3 s = vec3(7, 157, 113);
    vec3 ip = floor(p);
    vec4 h = vec4(0, s.yz, s.y + s.z) + dot(ip, s);
    p -= ip;
    p = p * p * (3. - 2. * p);
    
    vec4 h1 = hash4(h);
    vec4 h2 = hash4(h + s.x);
    vec4 res = mix(h1, h2, p.x);
    res.xy = mix(res.xz, res.yw, p.y);
    return mix(res.x, res.y, p.z);
}


float noise(vec2 x)
{
    vec2 p = floor(x);
    vec2 w = fract(x);
    vec2 u = w * w * w * (w * (w * 6.0 - 15.0) + 10.0);

    float a = hash12(p + vec2(0, 0));
    float b = hash12(p + vec2(1, 0));
    float c = hash12(p + vec2(0, 1));
    float d = hash12(p + vec2(1, 1));

    return (a + (b - a) * u.x + (c - a) * u.y + (a - b - c + d) * u.x * u.y);
}

float fbm3d(vec3 p)
{
    mat3 m = mat3(.0, .8, .6,
        -.8, .36, -.48,
        -.6, -.48, .64);
    float f = 0., s = .5;
    for (int i = 0; i < 4; ++i)
    {
        f += s * noise(p);
        p = p * m * 2.01;
        s *= .5;
    }

    return f;
}

void rot(inout vec2 p, float a) {
    p = cos(a) * p + sin(a) * vec2(p.y, -p.x);
}

void opU(inout vec2 obj, float b, float matB)
{
    obj = obj.x < b ? obj : vec2(b, matB);
}

//https://mercury.sexy/hg_sdf/
void opStairU(inout vec2 obj, float b, float r, float n, float matB)
{
    float s = r / n;
    float u = b - r;
    obj = vec2(min(min(obj.x, b), .5 * (u + obj.x + abs(mod(u - obj.x + s, 2. * s) - s))),
        obj.x < b ? obj.y : matB);
}

void opCut(inout vec2 obj, float b, float matB)
{
    obj = obj.x > -b ? obj : vec2(-b, matB);
}

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

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

vec2 polarRep(vec2 p, float n)
{
    if(n == 0.) return p;
    n = PI * .5 / n;
    float a = mod(atan(p.y, p.x) + n / 2., n) - n / 2.;
    p = length(p) * vec2(cos(a), sin(a));
    return .5 * (p + p - vec2(1, 0));
}

float ellipsoid(vec3 p, vec3 r)
{
    return (length(p / r) - 1.) * min(min(r.x, r.y), r.z);
}

float cone( vec3 p, vec2 c, float h )
{
  float q = length(p.xz);
  return max(dot(c.xy,vec2(q,p.y)),-h-p.y);
}

float disc(vec3 p, float r) 
{
	float l = length(p.xz) - r;
	return l < .0 ? abs(p.y) : length(vec2(p.y, l));
}

vec3 palette(float t, vec3 a, vec3 b, vec3 c, vec3 d)
{
    return a+b*cos(2.*PI*(c*t+d));
}

vec3 firecolor(float f)
{
    f=f*f*(3.0-2.0*f);
    return  min(vec3(f+.8, f*f*1.4+.1, f*f*f*.7)*f, 1.0);
}

vec2 ufo(vec3 p)
{
    vec2 res = vec2(0., 0.);
    vec3 q = p;
    q.y = -q.y;
    float d1 = cone(q, vec2(.5, .6), .75);
    res.x = d1;

    p.y -= 1.;
    float d2 = vcylinder(p, .3, 1.);

    opStairU(res, d2, 0.2, 3., 0.);


    q = p;
    q.y -= 1.;
    d1 = cone(q, vec2(.5,.6), .7);
    opStairU(res, d1, 0.1, 3., 0.);

    q.y -= 2.8;
    d1 = cone(q, vec2(1.,.05), 3.);
    opStairU(res, d1, .3, 3., 0.);
    return res;
}

vec2 perl(vec3 p, vec4 rnd, vec2 fl)
{
    //p *=5.;
    float width = 1.3;//rnd.x * 0.4 + 0.5;

    vec3 q = p;
    p.y -= width;
    vec2 obj = vec2(ellipsoid(p, vec3(width*.4, width*.35, width * .4)), 0.);

    //rot(q.xz, .6);
    q.xz=polarRep(q.xz, 1.);
    q.x+=.3 - width * 0.1;
    rot(q.xy, .2);

    float d1 = vcylinder(q, width, width * 0.08);
    opU(obj, d1, 0.);

    q = p;
    for(int i = 0; i < 4; ++i)
    {
        q.y -= width * .8;
        vec3 qq= q;
    
        qq.xz=polarRep(qq.xz, 1.);
        //qq.x += log(width*1.8+1.) * .35;
        qq.x += .4-width*.05;
        float d2 = vcylinder(qq, width, width*0.08);
        opU(obj, d2, 0.);

        //if(rnd.w > 0.5)
        {
            q.y -= width * .8+.1;

            d1 = ellipsoid(q, vec3(width*.4,width*.3,width*.4));
        }
        #if 0
        else
        {
            q.y -= width - float(i) * 0.1;
            
            d1 = ufo(q*(2./width)).x/(2./width);//
        }
        #endif
        opU(obj, d1, 0.);
    }
    

    q.y -= width * .5;
    d1 = vcylinder(q, width, width * 0.08);
    opU(obj, d1, 0.);
    d1 = vcylinder(q, width * 0.15, width * 0.1);
    opStairU(obj, d1, width * 0.1, 2., 0.);
    

    q.y -= width;
    //if(rnd.w > 0.5)
    {
        d1 = length(q) - width * .15;
    }
    #if 0
    else
    {
        d1 = ufo(q*(3./width)).x/(3./width);
    }
    #endif

    opU(obj, d1, 0.);
    
    q.y-=width * 1.5;
    d1 = cone(q, vec2(.9, .05), width * 1.5);
    opU(obj, d1, 0.);
    //obj.x/=5.;
    return obj;
}


vec2 building(vec3 p, vec4 rnds, vec2 fl)
{
    float rep = floor(pow(rnds.x,3.) * 6.);

    vec2 rnds2 = hash22(rnds.zw);
    float layer = floor(rnds2.x * 6.);
    
    float width = clamp(rnds.z*(1.+rep), 0.3, 1.);// + mix(2., 0., rnds.y);
    float height = max(0.3, 1. * rnds.w) + (1.-width) * .2;

    if(rep > 0.)
    {
        width *= .5;    
    }

    width = max(0.3, width);

    vec3 q = p;

    //rot(q.xz, params0.x);

    q.xz = polarRep(q.xz, rep);
    if(rep > 0.) q.x += width * .6;

    float bs = width * 0.4;

    vec3 size0 = vec3(bs, height * 1.2, bs);

    q.y -= size0.y;
    float d0 = box(q, size0);
    vec2 res = vec2(d0, 0.);

    vec3 size1 = vec3(bs * .9, height * .48, bs * .9);
    float lasty = size0.y, cury = size1.y;
    for(int i = 0; i < 6; ++i)
    {
        if(i >= int(layer)) break;
        q.y -= lasty + cury;
        float d1 = box(q, size1);
        opStairU(res, d1, bs * .1, 2., 0.);
        lasty = cury;

        size1.xz *= 0.9;
        size1.y *= 0.5;
        //size1 = vec3(size1.x-.1, size1.y*.5, size1.z - .1);
        size1 = max(vec3(bs*.2, height*.2, bs * .2), size1);
        cury = size1.y;
    }

    if(layer > 2. && rep == 0. && rnds.w > .95)
    {
        q.y -= lasty;
        vec3 qq = q;
        //rot(qq.xz, PI/4.);
        //float d1 = octahedron(qq, size1.x);
        vec2 res2 = ufo(qq / (size1.x*1.5));
        opStairU(res, res2.x * (size1.x*1.5), size1.x*0.5, 2., res2.y);
    }
    else if(layer >= 1. && rep < 2.)
    {
        q.y -= size1.y;
        float d1 = ellipsoid(q, vec3(size1.x*.9));
        opStairU(res, d1, size1.x * .2, 2., 0.);
    }
    else
    {
        q.y -= lasty; 
        opCut(res, box(q, vec3(size1.x, size1.x * .5, size1.z)), 0.);     
    }
    return  res;
}

vec2 map(vec3 p)
{
    vec3 q = p;
    vec2 fl = floor(q.xz);
    vec2 fr = fract(q.xz);

    vec4 rnds;
    rnds.xy = hash22(fl);
    rnds.zw = hash22(rnds.xy);
    q.xz = fract(q.xz) - .5;
    //q.xz -= 0.5;
    float dist = length(fl);
    float heightScale = max(.5, (2. - dist/40.));
    vec2 res=vec2(1e9, MAT_BACKGROUND);
    if(fl == vec2(0,0))
    {
        res = perl(q, vec4(0), fl);
    }
    else
    {
        if(dist < 80.)
        {
            res = building(vec3(q.x, q.y + sin(myTime+rnds.x * 2.)*0., q.z), rnds * vec4(1,1,1,heightScale), fl);
        }
       
    }
    
    q = p;
    q.x = mod(q.x-20., 50.) - 25.;
    opCut(res, box(q, vec3(1,100,100)), 0.);

    // q = p;
    // q.z = mod(q.z, 20. ) - 10.;
    // opCut(res, box(q, vec3(100,100,1)), 0.);
    
    opU(res, p.y, 1.);
    return res;
}

void blackholeRot(inout vec3 p, float flip)
{
    float t = -blackHoleTilt * flip;
    rot(p.yz, t*1.);
    rot(p.xy, -t);  
}

void blackholeTransform(inout vec3 p, float flip)
{
    p-=blackHolePos;
    blackholeRot(p, flip);
}

vec2 mapDisc(vec3 p)
{
    blackholeTransform(p, 1.);
    vec2 res = vec2(disc(p, blackHoleDiscRadius), MAT_DISC);
    return res;
}


vec3 getNormal(vec3 p, float t)
{
    vec3 eps = vec3(0.001 * t, 0, 0);

    return normalize(vec3(map(p + eps.xzy).x - map(p - eps.xyz).x,
        map(p + eps.yxz).x - map(p - eps.yxz).x,
        map(p + eps.zyx).x - map(p - eps.zyx).x));
}

float sabs(float x, float k)
{
    return sqrt(x * x + k) - 0.1;
}

vec2 vsabs(vec2 p, float k)
{
    return vec2(sabs(p.x, k), sabs(p.y, k));
}

float gerstner(vec2 p, float freq, vec2 phase, float sharpness)
{
    p += phase;
    p *= freq;
    p += noise(p) * 2. - 1.;
    vec2 wave = (1. - vsabs(sin(p), .005 - .005 * sharpness));
    return pow(wave.x + wave.y, 1.2);
}

mat2 m2 = mat2(0.80, 0.60,
    -0.60, 0.80);
    
float sand(vec3 pos)
{
    vec2 p = pos.xz;

    p *= 0.2;

    float amp = 1.;
    float freq = .8;
    float sharpness = 0.;

    vec2 offset = vec2(.5 + myTime * sandSpeed);

    amp *= .5;

    float h = sin(p.x*.3 + 2. * cos(p.y * .2)) * .2;

    for (int i = 0; i < 2; ++i)
    {
        h += amp * gerstner(p, freq, offset, sharpness);
        h += amp * gerstner(p, freq * 1.2, -offset * .5 + vec2(7.), sharpness);
        amp *= .3;

        freq *= 2.7;
        p = m2 * p;
        offset = m2 * .6 * offset;
        sharpness = .5;
    }
    return h;
}


float eps = 0.001;

vec3 sandNormal(vec3 pos, float t)
{
    vec2 e = vec2(eps * t, 0.0);
    return normalize(vec3(sand(pos - e.xyy) - sand(pos + e.xyy),
        2.0 * e.x,
        sand(pos - e.yyx) - sand(pos + e.yyx)));
}

// fast one for non rough sand
vec2 intersectsand(vec3 ro, vec3 rd)
{
    float nearT = 0.01;
    float farT = FAR;

    // if the end of the ray is above sand, return
    vec3 p = ro + farT * rd;
    float farHeight = p.y - sand(p);
    if (farHeight > 0.0) return vec2(farT, MAT_BACKGROUND);

    p = ro + nearT * rd;
    float nearHeight = p.y - sand(p);

    for (int i = 0; i < 256; i++)
    {
        // bisect
        float midT = mix(nearT, farT, nearHeight / (nearHeight - farHeight));

        p = ro + midT * rd;
        float midHeight = p.y - sand(p);

        // height of the midpoint is close to the sand enough 
        if (abs(midHeight) < 0.001) return vec2(midT, MAT_SAND);

        // update the bisection range based on current height
        if (midHeight < 0.0)
        {
            farT = midT;
            farHeight = midHeight;
        }
        else
        {
            nearT = midT;
            nearHeight = midHeight;
        }

        if (abs(nearT - farT) < 0.001) break;
    }


    return vec2(nearT, MAT_BACKGROUND);
}


float schlick(float f0, float NoV)
{
    return f0 + (1. - f0) * pow(1.0 - NoV, 5.);
}

// r is roughness
float directSpecular(vec3 V, vec3 L, vec3 N, float r)
{
    vec3 H = normalize(V + L);
    float NoH = max(0., dot(N, H));
    float NoV = max(0., dot(N, V));
    float NoL = max(0., dot(N, L));
    float r2 = r * r;
    float a = NoH * r2;
    float k = r2 / ((1.0 - NoH * NoH) + a * a);
    float d = k * k * (1.0 / PI);

    float v = 0.5 / (NoL * (NoV * (1. - r2) + r2) + NoV * (NoL * (1. - r2) + r2));

    float f = schlick(.02, NoV);

    return d * f;
}


float shadow(vec3 sro, vec3 srd)
{
    float s = 1., t = .001;
    for (int i = 0; i < 10; i+=1)
    {
        vec2 h = map(sro + srd * t);
        s = min(s, 16. * h.x / t);
        if (h.x < 1e-4) break;
        t += clamp(h.x*.5 + .5*hash13(sro+t*srd), .1, 2.);
    }

    return clamp(s, 0., 1.);
}


float mediumDensity(vec3 p, float h, float t)
{
    float fade = 10.;
    p.y += (noise(p*.3))*5.;
	return exp(-0.015 * p.y * p.y) * (1.0 - smoothstep(h-fade, h, p.y));
}

vec4 renderVolShadow(vec3 ro, vec3 rd, float t)
{
    float t0 = 0.,  st = min(t, 100.) / 24.;
    vec4 col = vec4(0);

    float h = 15.;
    if(ro.y > h)
    {
        t0 = abs(rd.y) < 0.0001 ? 0. : (ro.y - h) / (-rd.y);
    }
 
    t0 += hash12(gl_FragCoord.xy);

    for(int i = 0; i < 24; ++i)
    {
        vec3 p = ro + t0 * rd;
        float d = fbm3d(p*1.+ vec3(myTime * 1., myTime, -myTime*10.)*.5) * mediumDensity(p, h, t0) * fogDensity;

        float vis = shadow(p, sun_dir);
        col.rgb += d * st * vis* vec3(204,228,235)/255.* vec3(118,75,22)/255. * exp(-col.a);

        col.a = d * st;

        t0 +=st;
    }

    return col;
}

vec3 renderMist(vec3 ro, vec3 rd, float t)
{
    float t0 = 1000., st = 1., k = 0.;
    vec3 col = vec3(0);
    for(int i = 0; i < 32; ++i)
    {
        vec3 p = ro + t0 * rd;
        float density = fbm3d(p * .006) * .4;
        density *= density;
        col+=density * st * exp(-k) * palette(smoothstep(0.,.35, density), vec3(15.,11.,11.)/255.,
                                                vec3(218.,181.,181.)/255.,
                                                vec3(254.,254.,254.)/255.,
                                                vec3(194.,120.,0.)/255. ) * 2. ;
        k+=density * st;
        t0 += st;
    }
    col /= 16.;
    return col;
}


vec3 starLayer(vec2 uv, float r)
{
    vec3 col = vec3(0);
    vec2 fl = floor(uv);
    vec2 fr = fract(uv) -.5;

    float rnd = noise(fl);

    vec2 rndPos = (hash22(fl)-.5) * .3;
    r = (.5+.5*hash12(fl)) * r;

    if(rnd > 0.9)
    {
        float d = length(fr-rndPos) - r;
        
        d = exp(-d * 50.);
        
        float h = hash12(fl + 13.0);
        vec3 starCol = mix(vec3(1.0, 0.8, 0.5), vec3(0.6, 0.8, 1.0), h);
        col = starCol * d * 10. * (0.8 + 0.4 * hash12(fl + 7.0));

        float twinkleFreq = 2.0 + 8.0 * hash12(fl + 11.0);
        float twinklePhase = hash12(fl + 22.0) * 6.28;
        float twinkle = 0.7 + 0.6 * sin(myTime * .3 * twinkleFreq + twinklePhase);
        col *= twinkle;
    }

    return col;
}

vec2 toSperical(vec3 v)
{
    float theta = acos(clamp(v.y / length(v), -1.0, 1.0)) / PI;
    float phi   = atan(v.z, v.x) / PI * .5 + .5; 
    return vec2(theta, phi);
}

vec3 renderSky(vec3 ro, vec3 bendRd, vec3 rd, float t)
{
    vec3 dir = normalize(blackHolePos - ro);
    
    float distToBH = length(blackHolePos - ro);
    float shadowRadius = blackHoleRadius * 3.35; 
    float sinTheta = min(1.0, shadowRadius / max(0.1, distToBH));
    float cosTheta = sqrt(1.0 - sinTheta * sinTheta);

    float starFade= 1.;

    vec2 uv = toSperical(bendRd);
    uv.y*=1.5;

    // turn off some stars around the black hole for reducing flickering
    if(dot(rd, dir) > cosTheta && hackFlicker > 0.)
    {
        uv*=.7;
        starFade = 0.;
    }

    vec3 col = vec3(0);
    
    float baseFreq = 220.0; 

    // Layer 1
    col += starLayer(uv * baseFreq + vec2(0., time*.1), 0.001) * starFade * 0.5;
    
    // Layer 2
    vec2 uv2 = uv; rot(uv2, 1.1);
    col += starLayer((uv2 + vec2(0.13, 0.42)) * (baseFreq * 1.1), 0.0005) * 0.8 * starFade; 
    
    // Layer 3
    vec2 uv3 = uv; rot(uv3, 2.4);
    col += starLayer((uv3 + vec2(0.55, 0.91)) * (baseFreq * 0.95) - vec2(time * .1, 0.), 0.0005) * 0.6 * starFade;

    // Layer 4
    vec2 uv4 = uv; rot(uv4, -0.7);
    col += starLayer((uv4 + vec2(0.82, 0.26)) * (baseFreq * 1.05), 0.001) * 0.5 * starFade;

    // Layer 5
    vec2 uv5 = uv; rot(uv5, 3.8);
    col += starLayer((uv5 + vec2(0.37, 0.74)) * (baseFreq * 1.15), 0.001) * 0.4 * starFade;
    col += starLayer(uv5 * 150.+vec2(time*.2,0.), 0.005) * .8;

    col += renderMist(ro, bendRd, t);


    // trick for proper blending with the blackhole
    shadowRadius = blackHoleRadius * 2.0; 
    sinTheta = min(1.0, shadowRadius / max(0.1, distToBH));
    cosTheta = sqrt(1.0 - sinTheta * sinTheta);
    
    if(dot(rd, dir) > cosTheta)
        col = vec3(0);


    //  a little glow
    col += firecolor(.12).zyx * pow(max(0.,dot(bendRd, dir)),12.);

    return col;
}

vec2 intersect(vec3 ro, vec3 rd)
{
    float t = 0.;
    vec2 res = vec2(1e9, 0);

    int loopCount = 0;
    for (int i = 0; i < 256; i++)
    {
        vec3 p = ro + t * rd;
        res = map(p);

        if (res.x < 0.002 || t > FAR)
        {
            break;
        }

        // voxel
        float dx = -fract(p.x);
        if (rd.x > 0.)
        {
            dx = 1. - fract(p.x);
        }

        float dz = -fract(p.z);
        if (rd.z > 0.)
        {
            dz = 1. - fract(p.z);
        }

        float distToVoxelEdge = min(fract(dx / rd.x), fract(dz / rd.z)) + 0.1;

        distToVoxelEdge = max(0.1, distToVoxelEdge);

        t += min(res.x, distToVoxelEdge);
        loopCount++;
    }

    if (t > FAR || loopCount == 256){ res.y = MAT_BACKGROUND; t = FAR;}

    return vec2(t, res.y);
}

vec2 intersectBlackHole(inout vec3 ro, inout vec3 rd)
{
    float t =.1;
    vec2 res = vec2(1e9, 0);
    int loopCount = 0;
    ro += t*rd;

    bool hit = false;

    for (int i = 0; i < 256; i++)
    {
        res = mapDisc(ro);

        vec3 q = ro;
        blackholeTransform(q, 1.); 

        bool yeah = length(q) > blackHoleRadius *0.95;
        if (yeah&& (res.x < .0002 * t || t > FAR+1500.))
        {
            break;
        }

        if(!yeah)
        {
            res.x = max(res.x, 2.5 * hash13(ro));
        }
        
        t += res.x;

        
        //vec3 d = q - blackHolePos;
        float r = length(q);
        float bend = 3.  / (r*r);
        vec3 n = q / r;
        blackholeRot(n,  -1.);
        rd = normalize(rd - n * (bend * res.x*4.5));
        ro += rd * res.x;

        loopCount++;
    }

    if (t > FAR+1500. || loopCount == 256) 
    {
        res.y = MAT_BACKGROUND; 
    }
    return vec2(t, res.y);
}


vec3 window(vec2 p, inout vec3 n, float t)
{
    p*=20.;
    n = vec3(0);
    vec2 fl = floor(p);
    vec2 fr = fract(p);
    
    fr -= .5;
    
    float d = length(fr) - .1;
    
    vec2 rnds = hash22(fl);
    vec2 rnds1 = hash22(fl + vec2(1.,0.));
    vec2 rnds2 = hash22(fl + vec2(0.,1.));
    vec2 rnds3 = hash22(fl + vec2(1.,1.));

    float oddRow = 2.* fract(.5 * fl.x);

    vec2 rndsl = hash22(floor(p*.5));
    
       
    d = max(fr.x, fr.y);
    d = smoothstep(.2, .1, d); 

    float c = abs(sin(myTime) * .1);

    if(rnds.x > .5 * c && rnds.y < .5*c)
    {
        //emission = 20.* vec3(rnds.xxx * .2) * (1.-step(.35, max(fr.x,fr.y))) * oddRow;
    }
	
    // Analytic gradient normal offset
    float edge   = max(abs(fr.x), abs(fr.y));
    float tEdge  = clamp((edge - 0.35) / (0.50 - 0.35), 0.0, 1.0);
    // float dSmoothStep = (edge > 0.35 && edge < 0.50) ? (6.0 * tEdge * (1.0 - tEdge)) / (0.50 - 0.35) : 0.0;

    // vec2 gradDir = (abs(fr.x) > abs(fr.y)) ? 
    //                     vec2(sign(fr.x), 0.0) : (abs(fr.y) > abs(fr.x)) ?     
    //                             vec2(0.0, sign(fr.y)) : normalize(vec2(sign(fr.x), sign(fr.y))); // tie case

    // //float heightField = 0.08 * smoothstep(0.35, 0.50, edge);
    // vec2 gradHeightField = 0.08 * dSmoothStep * gradDir;

    // n += 1.* vec3(-gradHeightField.x, -gradHeightField.y, 0.0) * 1. * oddRow;

    return mix(vec3(.3,.3,.4), vec3(.1,.1,.2), (1.-tEdge) * oddRow);
}



vec3 getMat(vec3 p, inout vec3 n, float id, float t)
{
    vec3 col = vec3(.5, .5, .55)*.2;

    if(id == MAT_BUILDING)
    {
        col = vec3(.3,.3,.4);
        n = getNormal(p, t);
        if(n.y < 0.01)
        {
            vec2 uv = p.yz;
            if(abs(n.x) < abs(n.z))
            {
                uv=p.yx;
            }

            vec3 noff;
            col = window(uv, noff, t);

            if(abs(n.x) > abs(n.z))
            {
                noff.xyz = noff.zyx;
            }
            n -= .5*noff;
            n=normalize(n);
        }
    }
    else if (id == MAT_WINDOW)
    {
        n = getNormal(p, t);

        col = map(p).x * vec3(1000.);
        p.xz = fract(p.xz) - .5;
        col *= smoothstep(-.5, -.4, p.x);
        col *= smoothstep(-.5, -.4, p.z);
    }
    else if(id == MAT_SAND)
    {
        n = sandNormal(p, t);
        col = vec3(1., .9, .5) * (.8 + .2*hash12(p.xz * 50.));
    }
    else if(id == MAT_DISC)
    {
        vec3 q = p;

        blackholeTransform(q, 1.);
        
        vec2 uv = vec2(atan(q.z, q.x)/ PI * .5 + .5 , 
                    (length(q) - blackHoleRadius)/blackHoleDiscRadius);
        float f = fbm3d(vec3(uv.x * 50.+myTime, uv.y*50., q.y));
       // f *= smoothstep(-blackHoleDiscRadius/4., blackHoleDiscRadius/8.,(blackHoleDiscRadius - length(q.xz)));
        float intensity = smoothstep(0., .5, .7 - uv.y);

        col = vec3(f, intensity, 0.);
    }

    return col;
}


void getCamera(out vec3 ro, out vec3 ta, vec2 uv)
{
    sandSpeed = .1;

    if(myTime < 10.)
    {
        vec3 ro0 = vec3(46.87, 12.17, -99.36);
        vec3 ta0 = vec3(46.97, 12.23, -100.35);

        vec3 ro1 = vec3(72.43, 12.17, -96.76);
        vec3 ta1 = vec3(72.53, 12.23, -97.75);

        float t = myTime/10.;
        fogDensity = mix(0., .1, t);
        sandSpeed = smoothstep(.5,1.,t) * .05;
        ro = mix(ro0, ro1, smoothstep(0., 1., t));
        ta = mix(ta0, ta1, smoothstep(0., 1., t));

        fade = smoothstep(0.,2., myTime)-smoothstep(8.,10.,myTime);
    }
    else if(myTime < 20.)
    {
        vec3 ro0 = vec3(-27.41, 3.92, -89.24);
        vec3 ta0 = vec3(-26.41, 3.9, -89.28);

        vec3 ro1 = vec3(-25.84, 3.44,-33.98);
        vec3 ta1 = vec3(-24.85,3.51, -33.83);

        float t = myTime/20.;
        ro = mix(ro0, ro1, smoothstep(0., 1., t));
        ta = mix(ta0, ta1, smoothstep(0., 1., t));
        fogDensity=.4;
        renderBH = 0.;
        fade = smoothstep(10.,12., myTime)-smoothstep(18.,20.,myTime);
    }
     else if(myTime < 40.)
    {
        vec3 ro0 = vec3(16.38, 9.21, -10.13);
        vec3 ta0 = vec3(15.54, 9.21, -9.59);;

        vec3 ro1 = vec3(-13.14, 9.21, -15.71);
        vec3 ta1 = vec3(-12.54,9.20, -14.91);

        float t = max(0.,(myTime-20.)/20.);
        ro = mix(ro0, ro1, smoothstep(0., 1., t));
        ta = mix(ta0, ta1, smoothstep(0., 1., t));
        fogDensity=.4;
        blackHoleTilt = .3;
        hackFlicker = 1.;

        fade = smoothstep(20.,20., myTime)-smoothstep(38.,40.,myTime);
    }
    else if(myTime < 60.)
    {
        vec3 ro0 = vec3(-6.78, 12.96, -12.83);
        vec3 ta0 = vec3(-6.78, 13.09, -11.83);

        vec3 ro1 = vec3(-15.39, 12.24, -39.41);
        vec3 ta1 = vec3(-15.31, 12.37, -38.42);
        
        float t = clamp((myTime - 40.) / 20., 0., 1.);
        //float t_mapped = mix(0.8, 1.0, smoothstep(0., 1., t));
        
        ro = mix(ro0, ro1, t);
        ta = mix(ta0, ta1, t);
        fogDensity=.4;
        hackFlicker = 1.;
        fade = smoothstep(40.,40., myTime)-smoothstep(59.,60.,myTime);
    }
    else if(myTime < 80.)
    {
        blackHolePos = vec3(30, 40, 180.);
        blackHoleDiscRadius = 220.;
        
        vec3 ro0 = vec3(-26.79, 9.10, -12.66);
        vec3 ta0 = vec3(-26., 9.2, -12.05);

        vec3 ro1 = vec3(-24.7, 9.74, -22.62);
        vec3 ta1 = vec3(-23.9, 9.84, -21.93);
        
        float t = clamp((myTime - 60.) / 20., 0., 1.);
        
        ro = mix(ro0, ro1, smoothstep(0., 1., t));
        ta = mix(ta0, ta1, smoothstep(0., 1., t));
        fogDensity=.4;
        fade = smoothstep(60.,60., myTime)-smoothstep(79.,80.,myTime);

    }
    else if(myTime < 90.)
    {
        vec3 ro0 = vec3(-5.05, 4.39, -24.4);
        vec3 ta0 = vec3(-5.1, 4.41, -23.4);

        vec3 ro1 = vec3(-5.12, 3.83, 6.11);
        vec3 ta1 = vec3(-5.11, 3.9, 7.1);
        
        float t = clamp((myTime - 80.) / 10., 0., 1.);
        
        ro = mix(ro0, ro1, smoothstep(0., 1., t));
        ta = mix(ta0, ta1, smoothstep(0., 1., t));
        fogDensity=.4;
        blackHoleTilt = .3;
        fade = smoothstep(80.,82., myTime)-smoothstep(89.,90.,myTime);
    }
    else if(myTime < 110.)
    {
        vec3 ro0 = vec3(-11.25, 11.34, 40.25);
        vec3 ta0 = vec3(-11.6, 11.07, 41.15);

        vec3 ro1 = vec3(-0.49, 21.96,67.23);
        vec3 ta1 = vec3(0.04, 22.31, 68.);
        
        float t = clamp((myTime - 90.) / 20., 0., 1.);
        
        ro = mix(ro0, ro1, smoothstep(0., 1., t));
        ta = mix(ta0, ta1, smoothstep(0., 1., t));
        fogDensity=.4;
        blackHoleTilt = .3;
        fade = 1.;
    }
    else 
    {
        fogDensity = .4;
        blackHoleTilt = .2;
        // ro = rayOrg;
        // ta = rayTarget;

        ro = vec3(-0.49, 21.96,67.23);
        ta = vec3(0.04, 22.31, 68.);
        fade = 1.;
    }
}


vec3 scene(vec3 ro, vec3 rd)
{

    vec3 bg = vec3(.5);

    vec2 res = intersect(ro, rd);

    vec2 res2 = intersectsand(ro, rd);
    if(res2.x < res.x) res = res2;

    if(res.y == MAT_BACKGROUND) res.x = 1000.;

    vec3 bhro = ro, bhrd = rd;
    
    if(renderBH > 0.)
    {
        res2 = intersectBlackHole(bhro, bhrd);
        if(res2.x < res.x) res = res2;
    }

    vec3 col = renderSky(ro, bhrd, rd, res.x).rgb;
    vec3 pos = ro + res.x * rd;
    vec3 n;


    if (res.y > MAT_BACKGROUND)
    {
       if(res.y == MAT_DISC)
        {
            float Tr = 1.0;                  
            vec3 accumDisc = vec3(0.0); 

            // for strengthen the disc intensity close to the center
            vec3 acc2 = vec3(0);

            vec3 m = getMat(bhro, n, res.y, res.x);

            vec3 c = vec3(firecolor(m.x).zyx);
            
            accumDisc += Tr * m.x * c;
            Tr *= (1.0 - m.x);

            acc2 += firecolor(m.y*.8).zyx;
            

            // disc is transparent
            for (int k = 0; k < 1; ++k)
            {
                if(Tr <= 1e-3) 
                    break;
                bhro += 1.*bhrd;
                res2 = intersectBlackHole(bhro, bhrd);
                if (res2.y != MAT_DISC) break;

                m = getMat(bhro, n, MAT_DISC, res2.x);
                c = vec3(firecolor(m.x).zyx);
                accumDisc += Tr * m.x * c;
                Tr *= (1.0 - m.x);
                acc2 += firecolor(m.y * 0.8).zyx;
            }

            // Composite over background seen through remaining transmittance
            col = col * Tr + accumDisc + acc2;
            col.z *=.45;
        }
       
        else
        {
            vec3 m = getMat(pos, n, res.y, res.x);

            float dif = max(0., dot(n, sun_dir));
            float bac = max(0., dot(n, vec3(-sun_dir.x, -sun_dir.y, 0.)));
            float amb = max(0., n.y);
            float spe = directSpecular(-rd, sun_dir, n, .5);

            // specular cased flickering for most of the myTime
            if(!(myTime >=80.&&myTime<90.)&&res.y == MAT_BUILDING)
            {
                spe=0.;
            }

            col = ((dif+spe) + bac * .5) * m  + amb * vec3(0., .5, 1.) * .06 + vec3(1)*bac;
            col *= .5;
        }
    }

    vec4 volShadow = renderVolShadow(ro, rd, res.x);
    col = col * exp(-volShadow.a) + volShadow.rgb;

    return col;
}

vec3 tonemap(vec3 x) 
{
    const float a=2.51, b=0.03, c=2.43, d=0.59, e=0.14;
    return (x*(a*x+b))/(x*(c*x+d)+e);
}

void main()
{
    vec2 uv = gl_FragCoord.xy / resolution.xy;
    vec2 q = uv * 2. - 1.;
    q.x *= resolution.x / resolution.y;

    float d = length(q) - 1.;

    myTime = mod(time, 122.);

    // Time Warp 
    if(myTime >= 110.)
    {
        float t = myTime - 110.; 
    
        float easeBack = 110.0 - 1.1 * t * t;
        float warpStrength = smoothstep(0., 3.0, t);

        myTime = easeBack - pow(length(q), 2.0) * 5.0 * warpStrength;
        myTime = max(0.0, myTime);
    }

    vec3 ro, ta;

    getCamera(ro, ta, q);

    vec3 f = normalize(ta - ro);
    vec3 r = normalize(cross(f, vec3(0, 1, 0)));
    vec3 u = normalize(cross(r, f));

    vec3 rd = normalize(q.x * r + q.y * u + f * 2.);

    vec3 col = vec3(0, 0, 0);

    col = scene(ro, rd);

    col=tonemap(col);
    col=pow(clamp(col,0.0,1.0),vec3(0.45));
    col = col*.2+.8*col * col*(3.-2.*col); // contrast

    col *= pow(16.0 * uv.x * uv.y * (1.0 - uv.x) * (1.0 - uv.y), 0.2);
    col *= fade;

    gl_FragColor = vec4(col, 1.0);
}