// 'ringdans':
//   Inspired by the discontinued IKEA kids carpet 'ringdans', with animations but a very limited
//   selection of animals due to bad time planning capabilities.
// wungasaurus, Revision 2024, animated gif

// All SDFs from/based on https://iquilezles.org/articles/, MIT, iq
// [
float sdCircle(vec2 p, float r) {
  return length(p) - r;
}

float sdOrientedVesica(vec2 p, vec2 a, vec2 b, float w) {
  float r = 0.5 * length(b - a);
  float d = 0.5 * (r * r - w * w) / w;
  vec2  v = (b - a) / r;
  vec2  c = (b + a) * 0.5;
  vec2  q = 0.5 * abs(mat2(v.y, v.x, -v.x, v.y) * (p - c));
  vec3  h = (r * q.x < d * (q.y - r)) ? vec3(0.0, r, 0.0) : vec3(-d, 0.0, d + w);
  return length(q - h.xy) - h.z;
}

float sdEllipse(in vec2 p, in vec2 ab) {
  p = abs(p);
  if (p.x > p.y) {
    p  = p.yx;
    ab = ab.yx;
  }
  float l  = ab.y * ab.y - ab.x * ab.x;
  float m  = ab.x * p.x / l;
  float m2 = m * m;
  float n  = ab.y * p.y / l;
  float n2 = n * n;
  float c  = (m2 + n2 - 1.0) / 3.0;
  float c3 = c * c * c;
  float q  = c3 + m2 * n2 * 2.0;
  float d  = c3 + m2 * n2;
  float g  = m + m * n2;
  float co;
  if (d < 0.0) {
    float h  = acos(q / c3) / 3.0;
    float s  = cos(h);
    float t  = sin(h) * sqrt(3.0);
    float rx = sqrt(-c * (s + t + 2.0) + m2);
    float ry = sqrt(-c * (s - t + 2.0) + m2);
    co       = (ry + sign(l) * rx + abs(g) / (rx * ry) - m) / 2.0;
  } else {
    float h  = 2.0 * m * n * sqrt(d);
    float s  = sign(q + h) * pow(abs(q + h), 1.0 / 3.0);
    float u  = sign(q - h) * pow(abs(q - h), 1.0 / 3.0);
    float rx = -s - u - c * 4.0 + 2.0 * m2;
    float ry = (s - u) * sqrt(3.0);
    float rm = sqrt(rx * rx + ry * ry);
    co       = (ry / sqrt(rm - rx) + 2.0 * g / rm - m) / 2.0;
  }
  vec2 r = ab * vec2(co, sqrt(1.0 - co * co));
  return length(r - p) * sign(p.y - r.y);
}

vec3 sdJoint2DSphere(in vec2 p, in float l, in float a, float w) {
  if (abs(a) < 0.001) {
    float v  = p.y;
    p.y     -= clamp(p.y, 0.0, l);
    return vec3(length(p), p.x, v);
  }
  vec2  sc  = vec2(sin(a), cos(a));
  float ra  = 0.5 * l / a;
  p.x      -= ra;
  vec2  q   = p - 2.0 * sc * max(0.0, dot(sc, p));
  float u   = abs(ra) - length(q);
  float d   = (q.y < 0.0) ? length(q + vec2(ra, 0.0)) : abs(u);
  float s   = sign(a);
  float v   = ra * atan(s * p.y, -s * p.x);
  u         = u * s;
  if (v < 0.0) {
    if (s * p.x > 0.0) {
      v = abs(ra) * 6.283185 + v;
    } else {
      v = p.y;
      u = q.x + ra;
    }
  }
  return vec3(d - w, u, v);
}

float dot2(in vec2 v) {
  return dot(v, v);
}
float cro(in vec2 a, in vec2 b) {
  return a.x * b.y - a.y * b.x;
}
float sdBezier(in vec2 pos, in vec2 A, in vec2 B, in vec2 C) {
  vec2  a   = B - A;
  vec2  b   = A - 2.0 * B + C;
  vec2  c   = a * 2.0;
  vec2  d   = A - pos;
  float kk  = 1.0 / dot(b, b);
  float kx  = kk * dot(a, b);
  float ky  = kk * (2.0 * dot(a, a) + dot(d, b)) / 3.0;
  float kz  = kk * dot(d, a);
  float res = 0.0;
  float sgn = 0.0;
  float p   = ky - kx * kx;
  float q   = kx * (2.0 * kx * kx - 3.0 * ky) + kz;
  float p3  = p * p * p;
  float q2  = q * q;
  float h   = q2 + 4.0 * p3;
  if (h >= 0.0) { // 1 root
    h        = sqrt(h);
    vec2  x  = (vec2(h, -h) - q) / 2.0;
    vec2  uv = sign(x) * pow(abs(x), vec2(1.0 / 3.0));
    float t  = clamp(uv.x + uv.y - kx, 0.0, 1.0);
    vec2  q  = d + (c + b * t) * t;
    res      = dot2(q);
    sgn      = cro(c + 2.0 * b * t, q);
  } else { // 3 roots
    float z  = sqrt(-p);
    float v  = acos(q / (p * z * 2.0)) / 3.0;
    float m  = cos(v);
    float n  = sin(v) * 1.732050808;
    vec3  t  = clamp(vec3(m + m, -n - m, n - m) * z - kx, 0.0, 1.0);
    vec2  qx = d + (c + b * t.x) * t.x;
    float dx = dot2(qx), sx = cro(c + 2.0 * b * t.x, qx);
    vec2  qy = d + (c + b * t.y) * t.y;
    float dy = dot2(qy), sy = cro(c + 2.0 * b * t.y, qy);
    if (dx < dy) {
      res = dx;
      sgn = sx;
    } else {
      res = dy;
      sgn = sy;
    }
  }
  return sqrt(res) * sign(sgn);
}

float opUnion(float d1, float d2) {
  return min(d1, d2);
}
float opSubtraction(float d1, float d2) {
  return max(-d1, d2);
}
float opIntersection(float d1, float d2) {
  return max(d1, d2);
}
float opOnion(float sdf, float thickness) {
  return abs(sdf) - thickness;
}
float opSmoothExpUnion(float a, float b, float k) {
  float res = exp2(-k * a) + exp2(-k * b);
  return -log2(res) / k;
}

#define REPEATED(OUT_, SCENE_, REPEATED_P_, REPEATED_S_)                                                                                                       \
  do {                                                                                                                                                         \
    vec2  REPEATED_P = REPEATED_P_;                                                                                                                            \
    float REPEATED_S = REPEATED_S_;                                                                                                                            \
                                                                                                                                                               \
    vec2 REPEATED_ID = round(REPEATED_P / REPEATED_S);                                                                                                         \
    vec2 REPEATED_O  = sign(REPEATED_P - REPEATED_S * REPEATED_ID);                                                                                            \
                                                                                                                                                               \
    OUT_ = 1e20;                                                                                                                                               \
    for (int j = 0; j < 2; j++) {                                                                                                                              \
      for (int i = 0; i < 2; i++) {                                                                                                                            \
        vec2 REPEATED_RID = REPEATED_ID + vec2(i, j) * REPEATED_O;                                                                                             \
        vec2 REPEATED_R   = REPEATED_P - REPEATED_S * REPEATED_RID;                                                                                            \
        OUT_              = min(OUT_, SCENE_(REPEATED_R));                                                                                                     \
      }                                                                                                                                                        \
    }                                                                                                                                                          \
  } while (false)
// ]

// Noise functions (brown, grey)
// todo: unknown origin
// [
float noise3to1(in vec3 p) {
  mat3 m = mat3(1.0, 0.0, 0.0, 0.5, 1.2, 0.0, 0.0, 0.0, 1.0);
  vec3 s = m * p;
  return sin(s.x) * sin(s.y) * sin(s.z);
}
vec3 noise3to3(in vec3 p) {
  return vec3(noise3to1(p.xyz + vec3(1, 2, 3) * vec3(0.9, 0.7, 1.3)), noise3to1(p.zyx + vec3(7, 9, 8) * vec3(0.5, 1.2, 1.1)),
              noise3to1(p.yxz + vec3(3, 2, 5) * vec3(0.8, 0.3, 1.5)));
}

vec3 brown(in vec3 p) {
  return 0.5 * noise3to3(p) + 0.2 * noise3to3(p * 3.0) + 0.1 * noise3to3(p * 5.0);
}
vec3 grey(in vec3 p) {
  return (brown(brown(p * 0.1) * 5.0) - vec3(.5)) * 2.;
}
// ]

// Basic 2D math
// [
#define M_PI 3.1415926535897932384626433832795
vec2 rotate(vec2 uv, float th) {
  return mat2(cos(th), sin(th), -sin(th), cos(th)) * uv;
}
vec2 mirrorX(vec2 uv) {
  return vec2(-1., 1.) * uv;
}
vec2 mirrorY(vec2 uv) {
  return vec2(1., -1.) * uv;
}
// ]

// -------------------------------------------
float cocbpm   = 180.0;
float cocfac   = 1.0;
float looptime = 7.0;
float time() {
  return iTime * cocbpm / 30.0;
}
float csc() {
  return 0.015 * pow(sin(2.0 * time() * cocfac), 4.) + 0.95;
}
vec3 greyt(int idx, int var) {
  return grey(vec3(sin(2.0 * M_PI * time() / looptime), float(idx), float(var)));
}
float scrollOffs() {
  float stutter = 0.00 * pow(sin(time()), 4.0);
  return -1.5 - 1.0 / 1.0371 * time() / 3.5 + stutter;
}
float bgnoise(vec2 p) {
  //    p = sin(p/3.0);
  return 0.9 + 0.5 * (0.02358 * grey(p.xyy * 13.0).x + 0.0225 * grey(p.xyy * 71.0).x + 0.02394 * grey(p.xyy * 17.0).x);
}

// idx [0]
float sheepFoot(int idx, vec2 p, float a, float av, float b, float bv, float c, float cv) {
  float d = 1.;

  d = opUnion(sdJoint2DSphere(p, a + greyt(idx, 2).x * av, b + greyt(idx, 1).x * bv, c + greyt(idx, 0).x * cv).x, d);

  return d;
}

// idx [0…7]
float sheepCloud(int idx, vec2 p) {
  float d = 1.;

  d = opUnion(sdCircle(p + greyt(idx + 0, 1).xy * 0.02, 0.25 + greyt(idx + 0, 0).x * 0.01), d);
  d = opUnion(sdCircle(p - vec2(0.2, 0.2) + greyt(idx + 1, 1).xy * 0.02, 0.25 + greyt(idx + 1, 0).x * 0.01), d);
  d = opUnion(sdCircle(p - vec2(0.5, 0.1) + greyt(idx + 2, 1).xy * 0.03, 0.25 + greyt(idx + 2, 0).x * 0.02), d);
  d = opUnion(sdCircle(p - vec2(0.6, -0.1) + greyt(idx + 3, 1).xy * 0.01, 0.25 + greyt(idx + 3, 0).x * 0.01), d);
  d = opUnion(sdCircle(p - vec2(0.35, -0.25) + greyt(idx + 4, 1).xy * 0.03, 0.40 + greyt(idx + 4, 0).x * 0.03), d);
  d = opUnion(sdCircle(p - vec2(-0.1, -0.2) + greyt(idx + 5, 1).xy * 0.02, 0.25 + greyt(idx + 5, 0).x * 0.01), d);
  d = opUnion(sdCircle(p - vec2(0.05, -0.4) + greyt(idx + 6, 1).xy * 0.03, 0.25 + greyt(idx + 6, 0).x * 0.04), d);
  d = opUnion(sdCircle(p - vec2(0.25, -0.4) + greyt(idx + 7, 1).xy * 0.02, 0.25 + greyt(idx + 7, 0).x * 0.01), d);

  return d;
}

// idx [0…2]
float sheepHead(int idx, vec2 p) {
  float d = 1.;

  d = opUnion(sdOrientedVesica(p - vec2(-.5, -.15), vec2(0.25, 0.05) - greyt(idx + 0, 1).xy * 0.005, vec2(0.5, 0.), 0.06), d);
  d = opUnion(sdOrientedVesica(p - vec2(-.5, -.17), vec2(0.75, 0.05) - greyt(idx + 1, 1).xy * -0.01, vec2(0.5, 0.), 0.07), d);
  d = opUnion(sdEllipse(p - vec2(.0, .0), vec2(.14, .22) + greyt(idx + 2, 1).xy * .01), d);

  return d;
}

// idx [0…1]
float sheepEyes(int idx, vec2 p) {
  float d = 1.;

  d = opUnion(sdOrientedVesica(rotate(p, 1.24) - vec2(-.3, -.1), vec2(0.3, 0.05) - greyt(idx + 0, 1).xy * 0.005, vec2(0.35, 0.04), 0.013), d);
  d = opUnion(sdOrientedVesica(rotate(p, 1.25) - vec2(.02, -.03), vec2(0.01, 0.05) - greyt(idx + 1, 1).xy * -0.001, vec2(0.05, 0.03), 0.012), d);

  return d;
}

// idx [0…16]
float sheep(int idx, vec2 p, float sc) {
  float d = 1.;

  p += vec2(0.3, -0.1);
  p *= sc;

  d = opUnion(sheepCloud(idx + 0, p), d);
  d = opUnion(sheepFoot(idx + 8, mirrorY(p - vec2(-0.2, -0.05)), .45, .10, -1.3, 0.2, .05, 0.008), d);
  d = opUnion(sheepFoot(idx + 9, mirrorY(p - vec2(-0.2, -0.35)), .25, .05, -0.7, 0.2, .05, 0.009), d);
  d = opUnion(sheepFoot(idx + 10, rotate(mirrorY(p - vec2(0.65, -0.48)), 3.14 / 2.), .25, .10, -0.7, 0.2, .05, 0.01), d);
  d = opUnion(sheepFoot(idx + 11, rotate(mirrorY(p - vec2(0.85, -0.15)), 3.14 / 2.), .25, .10, -0.3, 0.2, .045, 0.01), d);
  d = opSubtraction(sheepHead(idx + 12, rotate(p - vec2(0.25, -0.05) - greyt(12, 1).xy * .01, 2.2 - sin(greyt(12, 0).x))), d);
  d = opUnion(sheepEyes(idx + 14, rotate(p - vec2(0.25, -0.05) - greyt(12, 1).xy * .01, 2.2 - sin(greyt(12, 0).x))), d);

  return d;
}

// idx [0]
float foxFoot(int idx, vec2 p) {
  float d = 1.;

  d = opUnion(opOnion(sdBezier(p, vec2(0.0, 0.0), vec2(0.0, -0.3), vec2(0.1, -0.25)), 0.008 * greyt(idx, 1).x + 0.03), d);

  return d;
}

// idx [0…5]
float fox(int idx, vec2 p, float sc) {
  float d = 1.;

  p *= sc;

  vec2 q = (p + vec2(.2, .0));

  if (p.x > 0.85) {
    p.y -= smoothstep(0., 1., 1.05 * cos(time() / 2.3) * 0.9 * (p.x - 0.85) * .7);
  }
  if (q.x > 0.65) {
    q.y += smoothstep(0., 1., 1.25 * sin(time() / 1.5) * 0.3 * (q.x - 0.65) * .7);
  }

  d = opUnion(sdOrientedVesica(p + vec2(.8, .1), vec2(-.3, -.2) - greyt(idx + 0, 1).xy * 0.005, vec2(0.5, 0.), 0.2), d);
  d = opSmoothExpUnion(sdOrientedVesica(p + vec2(.3, .2), vec2(-.2, .0) - greyt(idx + 0, 1).xy * 0.005, vec2(1.5, 0.1), 0.26), d, 128.);
  d = opSmoothExpUnion(sdOrientedVesica(q + vec2(.1, .2), vec2(-.2, .0) - greyt(idx + 0, 1).xy * 0.005, vec2(1.5, 0.1), 0.26), d, 128.);
  d = opUnion(sdOrientedVesica(p + vec2(-.3, .2), vec2(.3, .5), vec2(.4, 0.1), 0.08), d);
  d = opUnion(sdOrientedVesica(p + vec2(-.3, .2), vec2(.5, .47), vec2(.5, 0.05), 0.07), d);

  // feet
  d = opUnion(foxFoot(idx + 2, rotate(p + vec2(-.1, 0.35), 0.2 + 0.1 * (1. - sin(time())))), d);
  d = opUnion(foxFoot(idx + 3, rotate(p + vec2(.1, .35), 0.1 + 0.1 * (1. - sin(1.5 + time())))), d);
  d = opUnion(foxFoot(idx + 4, rotate(p + vec2(-.7, .3), 0.05 + 0.1 * (1. - sin(0.2 + time())))), d);
  d = opUnion(foxFoot(idx + 5, rotate(p + vec2(-.6, .35), 0.2 + 0.1 * (1. - sin(1.3 + time())))), d);

  // eyes
  vec2 pEyes = p + vec2(-0.2, .16);
  d          = opSubtraction(sdOrientedVesica(pEyes + vec2(0.05, .0) + greyt(idx + 0, 1).xy * 0.02, vec2(.5, .1), vec2(.5, 0.03), 0.02), d);
  d          = opSubtraction(sdOrientedVesica(pEyes + vec2(-.05, .01) + greyt((idx + 1), 1).xy * 0.02, vec2(.5, .1), vec2(.5, 0.023), 0.02), d);

  return d;
}

float snek(int idx, vec2 p, float sc) {
  float d = 1.0;

  p *= sc;

  float t = time();

  vec2 v0 = vec2(0.2, 0.4 + 0.11 * sin(5.1 + t));
  vec2 v1 = vec2(0.0, 0.4 + 0.18 * sin(5.1 + t));
  vec2 v2 = vec2(0.05, 0.0 + 0.07 * sin(0.0 + t));
  vec2 v3 = vec2(0.1, -0.8 + 0.03 * sin(1.2 + t));
  vec2 v4 = vec2(0.45, 0.05 + 0.03 * sin(2.1 + t));
  vec2 v5 = vec2(0.6, 0.4 + 0.04 * sin(3.4 + t));
  vec2 v6 = vec2(0.7, 0.1 + 0.04 * sin(3.4 + t));
  vec2 v7 = vec2(0.9, -0.74 + 0.04 * sin(5.1 + t));
  vec2 v8 = vec2(1.05, 0.3 + 0.08 * sin(5.1 + t));
  vec2 v9 = vec2(1.15, 0.8 + 0.05 * sin(5.1 + t));
  vec2 vA = vec2(1.3, 0.8 + 0.05 * sin(5.1 + t));
  vec2 vB = vec2(1.55, 0.75 + 0.05 * sin(5.1 + t));
  vec2 vC = vec2(1.55, 0.5 + 0.05 * sin(5.1 + t));

  p += vec2(1., 0.3);

  d = opUnion(opOnion(sdBezier(p, v0, v1, v2), 0.096), d);
  d = opUnion(opOnion(sdBezier(p, v2, v3, v4), 0.097), d);
  d = opUnion(opOnion(sdBezier(p, v4, v5, v6), 0.099), d);
  d = opUnion(opOnion(sdBezier(p, v6, v7, v8), 0.101), d);
  d = opUnion(opOnion(sdBezier(p, v8, v9, vA), 0.102), d);
  d = opUnion(opOnion(sdBezier(p, vA, vB, vC), 0.11), d);

  d = opSubtraction(sdOrientedVesica(p - vC - vec2(-0.5, 0.01 * sin(time())) + greyt(idx + 0, 1).xy * 0.02, vec2(.5, .1), vec2(.5, 0.03), 0.02), d);
  d = opSubtraction(
      sdOrientedVesica(p - vC - vec2(-0.54, -0.01 + 0.01 * sin(0.1 + time() / 2.0)) + greyt(idx + 0, 1).xy * 0.02, vec2(.5, .1), vec2(.5, 0.03), 0.02), d);

  return d;
}

float fw  = 1.7;
float shw = 1.8;
float snw = 1.7;

float scene(vec2 p) {
  float d = 1.;

  float sc = csc();

  d = opUnion(fox((0), p + vec2(0.3, 0.15), sc), d);
  d = opUnion(sheep((0 + 5), p - vec2(fw, 0.4), sc), d);
  d = opUnion(snek((0 + 5 + 16), p - vec2(fw + shw, 0.1), sc), d);
  d = opUnion(fox((0), p - vec2(-0.3 + fw + shw + snw, -0.15), sc), d);
  d = opUnion(sheep((0 + 5), p - vec2(fw + shw + snw + fw, 0.4), sc), d);
  d = opUnion(snek((0 + 5 + 16), p - vec2(fw + shw + snw + fw + shw, 0.1), sc), d);

  return d;
}

void mainImage(out vec4 fragColor, in vec2 fragCoord) {
  vec2 p = (2.0 * fragCoord - iResolution.xy) / iResolution.y;

// #define debug
#ifndef debug
  p.x += scrollOffs();
  p   *= vec2(1.2);

  float d;
  REPEATED(d, scene, p, fw + shw + snw);
#else
  float d = sheep(0, p, csc());
#endif

  // coloring
  vec3 col  = (d > 0.0) ? bgnoise(p) * vec3(57. / 255., 105. / 255., 38. / 255.) : vec3(181. / 255., 175. / 255., 159. / 255.);
  col      *= 1.0 - 0.1 * exp(-10.0 * abs(d));
  col      *= 1.8 + 0.01 * cos(150.0 * d);
  col       = mix(col, vec3(1.0), 1.0 - smoothstep(0.0, 0.009, abs(d)));

  fragColor = vec4(col, 1.0);
}
