// "Isolating Androids" by dr2 - 2020
// License: Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License

#define AA  1   // optional antialiasing

float PrBox2Df (vec2 p, vec2 b);
float PrSphDf (vec3 p, float r);
float PrRoundCylDf (vec3 p, float r, float rt, float h);
float PrTorusDf (vec3 p, float ri, float rc);
float Maxv3 (vec3 p);
mat3 StdVuMat (float el, float az);
vec2 Rot2D (vec2 q, float a);
float Fbm2 (vec2 p);

vec3 ltDir;
float tCur, dstFar, tRadB, tRadS;
int idObj;
const float pi = 3.14159;

#define VAR_ZERO min (iFrame, 0)

#define DMIN(id) if (d < dMin) { dMin = d;  idObj = id; }

float RobDf (vec3 p, float dMin)
{
  vec3 q;
  float szFac, rAngH, rAngA, d, s;
  s = sin (pi * tCur);
  rAngH = -0.7 * s;
  rAngA = 1.1 * s;
  szFac = 0.15;
  dMin /= szFac;
  p /= szFac;
  p.y -= 0.25;
  q = p;
  q.y -= 2.3;
  d = max (PrSphDf (q, 0.85), - q.y - 0.2);
  q = p;
  q.y -= 1.5;
  d = min (d, PrRoundCylDf (q.xzy, 0.9, 0.28, 0.75));
  q = p;
  q.x = abs (q.x);
  q.xy -= vec2 (1.05, 2.);
  q.yz = Rot2D (q.yz, rAngA * sign (p.x));
  q.y -= -0.5;
  d = min (d, PrRoundCylDf (q.xzy, 0.2, 0.15, 0.6));
  q = p;
  q.xz = Rot2D (q.xz, rAngH);
  q.x = abs (q.x);
  q.xy -= vec2 (0.3, 3.1);
  q.xy = Rot2D (q.xy, 0.2 * pi);
  q.y -= 0.25;
  d = min (d, PrRoundCylDf (q.xzy, 0.06, 0.04, 0.3));
  q = p;
  q.xy -= vec2 (0.4, 0.3);
  d = min (d, PrRoundCylDf (q.xzy, 0.25, 0.15, 0.65));
  q = p;
  q.x = - q.x;
  q.xy -= vec2 (0.4, 0.3);
  d = min (d, PrRoundCylDf (q.xzy, 0.25, 0.15, 0.65));
  DMIN (2);
  q = p;
  q.xz = Rot2D (q.xz, rAngH);
  q.x = abs (q.x);
  q -= vec3 (0.4, 2.7, 0.7);
  d = PrSphDf (q, 0.15);
  DMIN (3);
  return szFac * dMin;
}

float ObjDf (vec3 p)
{
  vec3 q;
  vec2 b;
  float dMin, d, a, t, ns;
  dMin = dstFar;
  p.yz = Rot2D (p.yz, 0.1 * pi);
  q = p;
  ns = 16.;
  t = 0.05 * tCur;
  d = length (vec2 (Rot2D (vec2 (length (q.xz) - tRadB, q.y),
     0.5 * atan (q.z, q.x)))) - tRadS;
  q.xz = Rot2D (q.xz, t);
  a = 2. * pi * (floor (ns * atan (q.z, - q.x) / (2. * pi)) + 0.5) / ns;
  q.xz = Rot2D (q.xz, a);
  q.x += tRadB;
  q.xy = Rot2D (q.xy, 0.5 * (a + t));
  b = vec2 (0.95 * tRadS, 0.65 * pi * tRadB / ns);
  d = max (d, - max (PrBox2Df (vec2 (abs (q.x) - tRadS, q.z), b),
     PrBox2Df (vec2 (abs (q.y) - tRadS, q.z), b)));
  DMIN (1);
  q.xy = Rot2D (q.xy, 2. * pi * (floor (4. * atan (q.y, - q.x) / (2. * pi)) / 4.)) -
     vec2 (-0.35, 0.06);
  q.xz = vec2 (q.z, - q.x);
  dMin = RobDf (q, dMin);
  return 0.7 * dMin;
}

float ObjRay (vec3 ro, vec3 rd)
{
  float dHit, d;
  dHit = 0.;
  for (int j = VAR_ZERO; j < 120; j ++) {
    d = ObjDf (ro + dHit * rd);
    if (d < 0.001 || dHit > dstFar) break;
    dHit += d;
  }
  return dHit;
}

vec3 ObjNf (vec3 p)
{
  vec4 v;
  vec2 e;
  e = vec2 (0.005, -0.005);
  for (int j = VAR_ZERO; j < 4; j ++) {
    v[j] = ObjDf (p + ((j < 2) ? ((j == 0) ? e.xxx : e.xyy) : ((j == 2) ? e.yxy : e.yyx)));
  }
  v.x = - v.x;
  return normalize (2. * v.yzw - dot (v, vec4 (1.)));
}

float TrObjDf (vec3 p)
{
  p.yz = Rot2D (p.yz, 0.1 * pi);
  return PrTorusDf (p.xzy, tRadS - 0.02, tRadB);
}

float TrObjRay (vec3 ro, vec3 rd)
{
  float dHit, d;
  dHit = 0.;
  for (int j = VAR_ZERO; j < 120; j ++) {
    d = TrObjDf (ro + dHit * rd);
    if (d < 0.0005 || dHit > dstFar) break;
    dHit += d;
  }
  return dHit;
}

vec3 TrObjNf (vec3 p)
{
  vec4 v;
  vec2 e;
  e = vec2 (0.005, -0.005);
  for (int j = VAR_ZERO; j < 4; j ++) {
    v[j] = TrObjDf (p + ((j < 2) ? ((j == 0) ? e.xxx : e.xyy) : ((j == 2) ? e.yxy : e.yyx)));
  }
  v.x = - v.x;
  return normalize (2. * v.yzw - dot (v, vec4 (1.)));
}

vec3 ErCol (vec3 rd)
{
  vec3 erDir, col, vn;
  float erRad, bs, ts;
  erDir = normalize (vec3 (0.02, -0.04, 1.));
  erRad = 0.04;
  col = vec3 (0.);
  bs = dot (rd, erDir);
  ts = bs * bs - 1. + erRad * erRad;
  if (ts > 0.) {
    ts = bs - sqrt (ts);
    if (ts > 0.) {
      vn = normalize ((ts * rd - erDir) / erRad);
      col = mix (vec3 (0.3, 0.4, 0.8), vec3 (1., 1., 0.95),
         smoothstep (0.2, 0.8, Fbm2 (6. * vn.xy + 7.1))) * (0.5 + 0.5 * max (- dot (vn, rd), 0.)) *
         (0.2 + 0.8 * max (dot (vn, ltDir), 0.));
    }
  }
  return col;
}

vec3 StarPat (vec3 rd, float scl)
{
  vec3 tm, qn, u;
  vec2 q;
  float f;
  tm = -1. / max (abs (rd), 0.0001);
  qn = - sign (rd) * step (tm.zxy, tm) * step (tm.yzx, tm);
  u = Maxv3 (tm) * rd;
  q = atan (vec2 (dot (u.zxy, qn), dot (u.yzx, qn)), vec2 (1.)) / pi;
  f = 0.57 * (Fbm2 (11. * dot (0.5 * (qn + 1.), vec3 (1., 2., 4.)) + 131.13 * scl * q) +
      Fbm2 (13. * dot (0.5 * (qn + 1.), vec3 (1., 2., 4.)) + 171.13 * scl * q.yx));
  return 8. * vec3 (1., 1., 0.8) * pow (f, 16.);
}

vec3 BgCol (vec3 rd)
{
  vec3 col;
  col = ErCol (rd);
  if (length (col) < 0.03) col += StarPat (rd, 16.);
  return col;
}

vec3 ShowScene (vec3 ro, vec3 rd)
{
  vec4 col4;
  vec3 col, colR, vn, roo;
  float dstObj, dstObjT;
  tRadB = 3.;
  tRadS = 0.8;
  roo = ro;
  dstObj = ObjRay (ro, rd);
  if (dstObj < dstFar) {
    ro += dstObj * rd;
    vn = ObjNf (ro);
    if (idObj == 1) col4 = vec4 (0.5, 0.5, 0.55, 0.2);
    else if (idObj == 2) col4 = vec4 (0., 1., 1., 0.2);
    else if (idObj == 3) col4 = vec4 (1., 0., 0., 0.2);
    col = col4.rgb * (0.2 + 0.2 * max (- dot (vn, ltDir), 0.) +
       0.8 * max (dot (vn, ltDir), 0.)) +
       col4.a * pow (max (dot (normalize (ltDir - rd), vn), 0.), 32.);
  } else {
    col = BgCol (rd);
  }
  ro = roo;
  dstObjT = TrObjRay (ro, rd);
  if (dstObjT < min (dstObj, dstFar)) {
    ro += dstObjT * rd;
    vn = TrObjNf (ro);
    col = mix (0.1 + 1.5 * BgCol (reflect (rd, vn)), col, 0.1 +
       0.9 * smoothstep (0.1, 0.7, - dot (rd, vn)));
  }
  return clamp (col, 0., 1.);
}

void mainImage (out vec4 fragColor, in vec2 fragCoord)
{
  mat3 vuMat;
  vec4 mPtr;
  vec3 ro, rd, col;
  vec2 canvas, uv;
  float el, az, zmFac, sr;
  canvas = iResolution.xy;
  uv = 2. * fragCoord.xy / canvas - 1.;
  uv.x *= canvas.x / canvas.y;
  tCur = iTime;
  mPtr = iMouse;
  mPtr.xy = mPtr.xy / canvas - 0.5;
  az = -0.07 * pi;
  el = -0.05 * pi;
  if (mPtr.z > 0.) {
    az += 2. * pi * mPtr.x;
    el += pi * mPtr.y;
  }
  vuMat = StdVuMat (el, az);
  ro = vuMat * vec3 (0., 0., -20.);
  zmFac = 6.;
  dstFar = 50.;
  ltDir = normalize (vec3 (1., 0.7, -1.));
#if ! AA
  const float naa = 1.;
#else
  const float naa = 3.;
#endif
  col = vec3 (0.);
  sr = 2. * mod (dot (mod (floor (0.5 * (uv + 1.) * canvas), 2.), vec2 (1.)), 2.) - 1.;
  for (float a = float (VAR_ZERO); a < naa; a ++) {
    rd = vuMat * normalize (vec3 (uv + step (1.5, naa) * Rot2D (vec2 (0.5 / canvas.y, 0.),
       sr * (0.667 * a + 0.5) * pi), zmFac));
    col += (1. / naa) * ShowScene (ro, rd);
  }
  fragColor = vec4 (pow (col, vec3 (0.8)), 1.);
}

float PrSphDf (vec3 p, float r)
{
  return length (p) - r;
}

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

float PrRoundCylDf (vec3 p, float r, float rt, float h)
{
  float dxy, dz;
  dxy = length (p.xy) - r;
  dz = abs (p.z) - h;
  return min (min (max (dxy + rt, dz), max (dxy, dz + rt)), length (vec2 (dxy, dz) + rt) - rt);
}

float PrTorusDf (vec3 p, float ri, float rc)
{
  return length (vec2 (length (p.xy) - rc, p.z)) - ri;
}

float Maxv3 (vec3 p)
{
  return max (p.x, max (p.y, p.z));
}

mat3 StdVuMat (float el, float az)
{
  vec2 ori, ca, sa;
  ori = vec2 (el, az);
  ca = cos (ori);
  sa = sin (ori);
  return mat3 (ca.y, 0., - sa.y, 0., 1., 0., sa.y, 0., ca.y) *
         mat3 (1., 0., 0., 0., ca.x, - sa.x, 0., sa.x, ca.x);
}

vec2 Rot2D (vec2 q, float a)
{
  vec2 cs;
  cs = sin (a + vec2 (0.5 * pi, 0.));
  return vec2 (dot (q, vec2 (cs.x, - cs.y)), dot (q.yx, cs));
}

const float cHashM = 43758.54;

vec2 Hashv2v2 (vec2 p)
{
  vec2 cHashVA2 = vec2 (37., 39.);
  return fract (sin (dot (p, cHashVA2) + vec2 (0., cHashVA2.x)) * cHashM);
}

float Noisefv2 (vec2 p)
{
  vec2 t, ip, fp;
  ip = floor (p);
  fp = fract (p);
  fp = fp * fp * (3. - 2. * fp);
  t = mix (Hashv2v2 (ip), Hashv2v2 (ip + vec2 (0., 1.)), fp.y);
  return mix (t.x, t.y, fp.x);
}

float Fbm2 (vec2 p)
{
  float f, a;
  f = 0.;
  a = 1.;
  for (int j = 0; j < 5; j ++) {
    f += a * Noisefv2 (p);
    a *= 0.5;
    p *= 2.;
  }
  return f * (1. / 1.9375);
}