// License CC0: The Monolith, take II
//  Result of messing around with refractions and gas giants when needing a break

#define PI                3.141592654
#define TAU               (2.0*PI)
#define TIME              iTime
#define RESOLUTION        iResolution
#define ROT(a)            mat2(cos(a), sin(a), -sin(a), cos(a))
#define L2(x)             dot(x, x)
#define PSIN(x)           (0.5 + 0.5*sin(x))
#define PCOS(x)           (0.5 + 0.5*cos(x))
#define TTIME             (TAU*TIME)

const float miss          = 1E6;
const float refrIndex     = 0.8;
const vec3  lightPos      = 200.0*vec3(1.5, 2.0, 1.0);

const vec3 skyCol1        = vec3(0.2, 0.4, 0.6);
const vec3 skyCol2        = vec3(0.4, 0.7, 1.0);
const vec3 sunCol         =  vec3(8.0,7.0,6.0)/8.0;

const vec3 boxDim         = vec3(1.0, 9.0, 4.0)/18.0;
const vec4 plane          = vec4(vec3(0.0, 1.0, 0.0), 0.5);

const vec4  s_giant     = vec4(vec3(-100.0, -45.0, -50.0), 40.0);
const vec4  s_miniGiant = vec4(vec3(0.0), 0.4);
const vec4  s_planet    = vec4(0.6*vec3(1.0, 0., 0.0), 0.05);
const vec3  p_normal    = normalize(vec3(0.4, 1.0, 0.0));
const vec4  p_rings     = vec4(p_normal, 0.0);
const vec3  innerSunDir = normalize(vec3(-0.0, 0.25, 0.5));

const mat2  rot0        = ROT(0.0);

float tanh_approx(float x) {
//  return tanh(x);
  float x2 = x*x;
  return clamp(x*(27.0 + x2)/(27.0+9.0*x2), -1.0, 1.0);
}

// https://iquilezles.org/articles/smin
float pmin(float a, float b, float k) {
  float h = clamp(0.5+0.5*(b-a)/k, 0.0, 1.0);

  return mix(b, a, h) - k*h*(1.0-h);
}

// https://stackoverflow.com/questions/15095909/from-rgb-to-hsv-in-opengl-glsl
vec3 hsv2rgb(vec3 c) {

  const vec4 K = vec4(1.0, 2.0 / 3.0, 1.0 / 3.0, 3.0);
  vec3 p = abs(fract(c.xxx + K.xyz) * 6.0 - K.www);
  return c.z * mix(K.xxx, clamp(p - K.xxx, 0.0, 1.0), c.y);
}

vec4 alphaBlend(vec4 back, vec4 front) {
  float w = front.w + back.w*(1.0-front.w);
  vec3 xyz = (front.xyz*front.w + back.xyz*back.w*(1.0-front.w))/w;
  return w > 0.0 ? vec4(xyz, w) : vec4(0.0);
}

vec4 alphaBlend(vec4 a, float adist, vec4 b, float bdist, out float dist) {
  bool t = adist < bdist;
  vec4 c = alphaBlend(t ? b : a, t ? a : b);
  dist = min(adist, bdist);
  return c;
}

// http://mercury.sexy/hg_sdf/
vec2 mod2(inout vec2 p, vec2 size) {
  vec2 c = floor((p + size*0.5)/size);
  p = mod(p + size*0.5,size) - size*0.5;
  return c;
}

float hash1(vec2 co) {
  return fract(sin(dot(co, vec2(12.9898,58.233))) * 13758.5453);
}

vec2 hash2(vec2 p) {
  p = vec2(dot (p, vec2 (127.1, 311.7)), dot (p, vec2 (269.5, 183.3)));
  return -1. + 2.*fract (sin (p)*43758.5453123);
}

// https://iquilezles.org/articles/intersectors
vec2 raySphere(vec3 ro, vec3 rd, vec4 s) {
    vec3 ce = s.xyz;
    float ra = s.w;
    vec3 oc = ro - ce;
    float b = dot( oc, rd );
    float c = dot( oc, oc ) - ra*ra;
    float h = b*b - c;
    if( h<0.0 ) return vec2(miss); // no intersection
    h = sqrt( h );
    return vec2( -b-h, -b+h );
}

// https://iquilezles.org/articles/intersectors
vec2 rayBox(vec3 ro, vec3 rd, vec3 boxSize, out vec3 outNormal )  {
    vec3 m = 1.0/rd; // can precompute if traversing a set of aligned boxes
    vec3 n = m*ro;   // can precompute if traversing a set of aligned boxes
    vec3 k = abs(m)*boxSize;
    vec3 t1 = -n - k;
    vec3 t2 = -n + k;
    float tN = max( max( t1.x, t1.y ), t1.z );
    float tF = min( min( t2.x, t2.y ), t2.z );
    if( tN>tF || tF<0.0) return vec2(miss); // no intersection
    outNormal = -sign(rd)*step(t1.yzx,t1.xyz)*step(t1.zxy,t1.xyz);
    return vec2( tN, tF );
}

// https://iquilezles.org/articles/intersectors
float rayPlane(vec3 ro, vec3 rd, vec4 p ) {
  return -(dot(ro,p.xyz)+p.w)/dot(rd,p.xyz);
}

float circle(vec2 p, float r) {
  return length(p) - r;
}

vec3 toSpherical(vec3 p) {
  float r   = length(p);
  float t   = acos(p.z/r);
  float ph  = atan(p.y, p.x);
  return vec3(r, t, ph);
}

float noise(vec2 p) {
  float a = sin(p.x);
  float b = sin(p.y);
  float c = 0.5 + 0.5*cos(p.x + p.y);
  float d = mix(a, b, c);
  return d;
}

// https://iquilezles.org/articles/fbm
float fbm(vec2 p) {
  const mat2 frot = mat2(0.80, 0.60, -0.60, 0.80);

  float f = 0.0;
  float a = 1.0;
  float s = 0.0;
  float m = 2.0;
  for (int x = 0; x < 2; ++x) {
    f += a*noise(p);
    p = frot*p*m;
    m += 0.01;
    s += a;
    a *= 0.5;
  }

  return f/s;
}

vec3 normal(vec2 p) {
  vec2 v;
  vec2 w;
  vec2 e = vec2(0.0001, 0);

  vec3 n;
  n.x = fbm(p + e.xy) - fbm(p - e.xy);
  n.y = 2.0*e.x;
  n.z = fbm(p + e.yx) - fbm(p - e.yx);

  return normalize(n);
}

// https://iquilezles.org/articles/warp
float warp(vec2 p, float e, float ttime, out vec2 v, out vec2 w) {
  vec2 vx = vec2(0.0, 0.5)*e;
  vec2 vy = vec2(3.2, 1.3)*e;

  vec2 wx = vec2(1.7, 9.2)*e;
  vec2 wy = vec2(8.3, 2.8)*e;

  vx *= ROT(ttime/1000.0);
  vy *= ROT(ttime/900.0);

  wx *= ROT(ttime/800.0);
  wy *= ROT(ttime/700.0);

  v = vec2(fbm(p + vx), fbm(p + vy));
  w = vec2(fbm(p + -3.0*v + wx), fbm(p + 3.0*v + wy));
  return fbm(p + vec2(2.25, 1.25)*w);
}

vec4 planetDim(float ttime) {
  vec4 s_p = s_planet;
  s_p.xz *= ROT(-(ttime-26.0*TAU)*0.01);
  return s_p;
}

vec3 postProcess(vec3 col, vec2 q)  {
  col=pow(clamp(col,0.0,1.0),vec3(1.0/2.2));
  col=col*0.6+0.4*col*col*(3.0-2.0*col);  // contrast
  col=mix(col, vec3(dot(col, vec3(0.33))), -0.4);  // satuation
  col*=0.5+0.5*pow(19.0*q.x*q.y*(1.0-q.x)*(1.0-q.y),0.7);  // vigneting
  return col;
}

vec4 stars(vec3 ro, vec3 rd, vec3 nrd) {
  vec3 col = vec3(0.0);
  vec3 srd = toSpherical(rd.xzy);

  const float m = 3.0;
  float a = 0.0;

  for (float i = 0.0; i < m; ++i) {
    vec2 pp = srd.yz+0.5*i;
    float s = i/(m-1.0);
    vec2 dim  = vec2(mix(0.05, 0.003, s)*PI);
    vec2 np = mod2(pp, dim);
    vec2 o  = hash2(np+127.0+i);
    float y = sin(srd.y);
    pp += o*dim*0.5;
    pp.y *= y;
    float l = length(pp);
    float l2 = L2(pp);

    float h1 = hash1(np+127.0+i);
    float h2 = fract(h1*113.0);
    float h3 = fract(h1*127.0);

    vec3 hsv = vec3(fract(0.025-0.4*h1*h1), mix(0.5, 0.125, s), 1.0);
    vec3 scol = mix(8.0*h2, 0.25*h2*h2, s)*hsv2rgb(hsv);

    if (h3 < y) {
      col += exp(-(4000.0/mix(1.0, 0.2, s))*max(l-0.001, 0.0))*scol;
    }
  }

  return vec4(col, 1.0);
}

float gasGiantShade(vec4 gasGiant, vec3 p) {
  vec2 rsi = raySphere(p, innerSunDir, gasGiant);

  // Fake, fake softshadows
  return rsi.x >= 0.0 && rsi.x < miss ? exp(-3.0*(rsi.y - rsi.x)/(rsi.x+0.1)) : 1.0;
}

float planetShade(vec3 p, float ttime) {
  vec4 s_p = planetDim(ttime);
  vec2 rsi = raySphere(p, innerSunDir, s_p);

  // Fake, fake softshadows
  return rsi.x >= 0.0 && rsi.x < miss ? exp(-3.0*(rsi.y - rsi.x)/(rsi.x+0.1)) : 1.0;
}

vec4 gasGiant(vec4 gasGiant, vec3 ro, vec3 rd, vec3 nrd, vec3 pmod, float ttime, out float dist) {
  vec2 rsi = raySphere(ro, rd, gasGiant);

  dist = rsi.x;

  if (rsi.x == miss || rsi.x <= 0.0) {
    return vec4(0.0);
  }

  vec3 pi = (ro + rd*rsi.x)-gasGiant.xyz;
  vec3 ni = normalize(pi);
  vec3 ppi = toSpherical(pi.xzy*pmod);

  float rings = 1.0;
  for (int i = 0; i < 4; ++i) {
    rings *= sin(-ppi.z*0.825+3.*pi.y/gasGiant.w*sqrt(float(i+3)));
  }

  rings = tanh_approx((rings)*1.5);
  float yf = smoothstep(.6, 1.0, abs(pi.y)/gasGiant.w);
  yf *= yf;

  vec2 v;
  vec2 w;
  float h = warp((ppi.zy+vec2(-0.005*ttime, 0.1))*2.0*vec2(1.0, 6.0), rings, ttime, v, w);
  rings *= rings;

  float ps = planetShade(pi, ttime);

  vec3 col1 = mix(vec3(0.5, 0.25, 0.2), vec3(0.5, 0.15, 0.2), rings);
  vec3 col2 = mix(vec3(0.3, 0.3, 0.15), vec3(0.3, 0.3, 0.25), rings);

  float diff = max(dot(ni, innerSunDir), 0.0);
  vec3 ccol = pow(diff, 0.75)*tanh_approx(pow(abs(h + 0.5), 1.5)) + (length(v)*col1 + length(w)*col2);

  vec3 col = vec3(0.0);
  col += ccol;
  col += vec3(0.7, 0.9, 0.8)*0.8*rings;
  col = max(col, 0.0);
  col = mix(col, vec3(0.85), yf);

  col *= ps;
  col *= mix(0.000, 1.0, diff);

  return vec4(col, smoothstep(0.0, 0.25, (rsi.y - rsi.x)/gasGiant.w));
}

vec4 planet(vec3 ro, vec3 rd, vec3 nrd, vec3 pmod, float ttime, out float dist) {
  vec4 s_p = planetDim(ttime);
  vec2 rsi = raySphere(ro, rd, s_p);

  dist = rsi.x;

  if (rsi.x == miss || rsi.x <= 0.0) {
    return vec4(0.0);
  }

  vec3 pi = ro + rd*rsi.x;
  vec3 ni = normalize(pi - s_p.xyz);
  vec3 ppi = toSpherical((pi-s_p.xyz).zxy*pmod);
  float h = fbm(ppi.yz*5.0);

  float yf = smoothstep(.6, 1.0, abs(pi.y)/s_planet.w);
  yf *= yf;

  float gs = gasGiantShade(s_miniGiant, pi);

  float diff = (pow(max(dot(ni, innerSunDir), 0.0), 1.0));
  vec3 col = mix(0.9*vec3(0.5, 0.5, 0.45), 0.1*vec3(0.5, 0.5, 0.8), pow(h*h, 2.0)*(1.0-yf));
  col *= mix(0.0125, 1.5, diff);
  col *= gs;

  return vec4(col, smoothstep(0.0, 0.5, (rsi.y - rsi.x)/s_planet.w));
}


vec4 rings(vec3 ro, vec3 rd, vec3 nrd, float ttime, out float dist) {
  float rpi = rayPlane(ro, rd, p_rings);

  dist = rpi;

  vec3 col = vec3(0.0);

  vec3 pi = ro + rd*rpi;
  vec3 npi = ro + nrd*rpi;
  float aa = length(npi-pi);

  vec2 pp = pi.xz;
  float r = length(pp);

  vec2 npp = npi.xz;
  float nr = length(npp);

  float d = circle(pp, 0.65);
  d = abs(d) - 0.06;

  float rings = 1.0;
  float nrings = 1.0;
  for (int i = 0; i < 5; ++i) {
    rings *= sin(100.0*r*sqrt(float(i+3)));
    nrings *= sin(100.0*nr*sqrt(float(i+3)));
  }

  float diff = max(dot(p_rings.xyz, innerSunDir), 0.0);

  rings = (max((rings + nrings), 0.0))/2.0;

  float gs = gasGiantShade(s_miniGiant, pi);
  float ps = planetShade(pi, ttime);

  float i = mix(0.0, 1.0, smoothstep(-aa, aa, -d))*rings;

  vec3 hsv = vec3(0.15+0.1*sin(-50.0*(r+nr)), 0.3*PSIN(200.0*(r+nr)), gs*ps*diff);
  col += hsv2rgb(hsv);

  return vec4(col, tanh_approx(i));
}

vec3 innerRender(vec3 ro, vec3 rd, vec3 nrd, vec3 enor, vec2 bd) {
  float ttime = TTIME;

  vec3 pmod = mix(vec3(-1.0), vec3(1.0), step(vec3(0.0), enor));
  ro += mix(0.2, 0.7, fract(enor.x*sqrt(3.0)+enor.y*sqrt(5.0)+enor.z*sqrt(7.0)))*enor/bd.x;
  float gdist;
  vec4 gcol = gasGiant(s_miniGiant, ro, rd, nrd, pmod, ttime, gdist);

  float pdist;
  vec4 pcol = planet(ro, rd, nrd, pmod, ttime, pdist);

  float rdist;
  vec4 rcol = rings(ro, rd, nrd, ttime, rdist);

  float dist = gdist;
  vec4 col = gcol;

  col = alphaBlend(col, dist, pcol, pdist, dist);
  col = alphaBlend(col, dist, rcol, rdist, dist);

  return col.xyz*col.w;
}

vec3 skyColor(vec3 ro, vec3 rd, vec3 nrd) {
  vec3 lightDir = normalize(lightPos - ro);
  vec3 back = vec3(0.0);

  float gdist;
  vec3 grd = rd;
  vec3 gnrd = nrd;
  grd.y = -grd.y;
  gnrd.y = -gnrd.y;
  if (rd.y > 0.0) {
    vec4 scol = stars(ro, grd, gnrd);
    vec4 gcol = gasGiant(s_giant, ro, grd, gnrd, vec3(1.0), 2000.0, gdist);
    vec4 col = alphaBlend(scol, gcol);
    return col.xyz*col.w;
  }

  float tp  = rayPlane(ro, rd, plane);

    vec3 pos  = ro + tp*rd;
    vec3 ld   = normalize(lightPos - pos);
    vec3 snor;
    vec2 rb   = rayBox(pos, ld, boxDim, snor);
    vec3 spos = pos + ld*rb.x;
    float it  = rb.y - rb.x;
    float sha = rb.x == miss ? 1.0 : (1.0-tanh_approx(it*6.0/(0.1+rb.x)));
    vec3 nor  = vec3(0.0, 1.0, 0.0);
    vec2 pos2 = pos.xz;
    float h   = fbm(pos2)+1.0;
    vec3 n    = normal(pos2);
    float diff= pow(max(dot(n, ld), 0.0), 1.0);
    vec3 icol = mix(0.5, 1.0, diff*h)*4.0*vec3(0.5, 0.5, 0.45)*sunCol*sha*dot(-rd, nor);

    vec3  pcol= vec3(0.5);

    vec3 col  = icol*pcol;
    col = clamp(col, 0.0, 1.0);
    float f   = exp(-0.05*(max(tp, 0.0)));
    return mix(back, col , f);
}


vec3 render(vec3 ro, vec3 rd, vec3 nrd) {
  vec3 skyCol = skyColor(ro, rd, nrd);
  vec3 col = vec3(0.0);

  float t   = 1E6;
  vec3 nor;
  vec2 rb = rayBox(ro, rd, boxDim, nor);

  if (rb.x < miss) {
    t = rb.x;
    float itd = rb.y - rb.x;
    vec3 pos  = ro + t*rd;
    vec3 anor = abs(nor);
    vec2 tp   = anor.x == 1.0 ? pos.yz : (anor.y == 1.0 ? pos.xz : pos.xy);
    vec2 bd   = anor.x == 1.0 ? boxDim.yz : (anor.y == 1.0 ? boxDim.xz : boxDim.xy);
    vec3 refr = refract(rd, nor, refrIndex);
    vec3 nrefr= refract(nrd, nor, refrIndex);
    vec3 refl = reflect(rd, nor);
    vec3 nrefl= reflect(nrd, nor);
    vec3 rcol = skyColor(pos, refl, nrefl);
    float fre = mix(0.0, 1.0, pow(1.0-dot(-rd, nor), 3.0));
    vec3 ld   = normalize(lightPos - pos);
    float dif = pow(max(dot(nor,ld),0.0), 3.0);
    float spe = pow(max(dot(reflect(-ld, nor), -rd), 0.), 50.);
    float lin = mix(0.0, 1.0, dif);
    vec3 lcol = 2.0*sqrt(sunCol);

    col = innerRender(pos, refr, nrefr, nor, bd);

    vec2 btp = (1.0*bd - abs(tp));
    float bdd = pmin(btp.x, btp.y, 0.0125);
    float bddd = exp(-10000.0*bdd*bdd);
    col *= 1.0 - bddd;

    vec3 diff = hsv2rgb(vec3(0.7, fre, 0.075*lin))*lcol;
    col += fre*rcol+diff+spe*lcol;
    if (refr == vec3(0.0)) {
      // Not expected to happen as the refraction index < 1.0
      col = vec3(1.0, 0.0, 0.0);
    }

    col = mix(col, skyCol, tanh_approx(1E-5/(itd*itd)));
  } else {
    // Ray intersected sky
    return skyCol;
  }

  return col;
}

void mainImage(out vec4 fragColor, in vec2 fragCoord) {
  vec2 q = fragCoord.xy/RESOLUTION.xy;
  vec2 p = -1.0 + 2.0*q;
  p.x *= RESOLUTION.x/RESOLUTION.y;
  vec2 np = p + 4.0/RESOLUTION.yy;

  vec3 ro = mix(0.5, 0.25, PCOS(TTIME/120.0*sqrt(3.0)))*vec3(2.0, 0, 0.2)+vec3(0.0, -0.125, 0.0);
  ro.yx *= ROT(-sin(TTIME/120.0*sqrt(0.5))*0.25);
  ro.xz *= ROT((TAU*(TIME-14.0)/120.0));
  vec3 la = vec3(0.0, 0.0, 0.0);

  vec3 ww = normalize(la - ro);
  vec3 uu = normalize(cross( vec3(0.0,1.0,0.0), ww));
  vec3 vv = normalize(cross(ww,uu));
  float rdd = 2.0;

  vec3 rd = normalize(p.x*uu + p.y*vv + rdd*ww);
  vec3 nrd = normalize(np.x*uu + np.y*vv + rdd*ww);

  vec3 col = render(ro, rd, nrd);
  col = mix(vec3(0.0), clamp(col, 0.0, 1.0), smoothstep(5.0, 12.0, TIME));
  fragColor = vec4(postProcess(col, q),1.0);
}