// Kestrel, a stand-in world from zaisma.com/developers/backgrounds. MIT licensed: free to start from.

float h1(float n) { return fract(sin(n * 127.1) * 43758.5453); }
float h2(vec2 p) { return fract(sin(dot(p, vec2(127.1, 311.7))) * 43758.5453); }
vec2 h22(vec2 p) { p = vec2(dot(p, vec2(127.1, 311.7)), dot(p, vec2(269.5, 183.3))); return fract(sin(p) * 43758.5453); }
float noise(vec2 p) { vec2 i = floor(p), f = fract(p); vec2 u = f * f * (3.0 - 2.0 * f);
  return mix(mix(h2(i), h2(i + vec2(1.0, 0.0)), u.x), mix(h2(i + vec2(0.0, 1.0)), h2(i + vec2(1.0, 1.0)), u.x), u.y); }
float fbm(vec2 p) { float v = 0.0, a = 0.5; for (int i = 0; i < 5; i++) { v += a * noise(p); p = p * 2.03 + 1.7; a *= 0.5; } return v; }
// Soft particles in a grid of cells: each cell holds at most one, at a random spot, drifting with the grid.
float motes(vec2 p, float size, float chance, float seed) {
  vec2 id = floor(p), f = fract(p) - 0.5; vec2 r = h22(id + seed) - 0.5;
  float on = step(1.0 - chance, h2(id + seed * 1.3));
  float d = length(f - r * 0.7);
  return on * smoothstep(size, 0.0, d);
}

float edge(vec2 p, vec2 a, vec2 b) { vec2 e = b - a; return (e.y * (p.x - a.x) - e.x * (p.y - a.y)) / length(e); }
float quad(vec2 p, vec2 a, vec2 b, vec2 c, vec2 d) { return max(max(edge(p, a, b), edge(p, b, c)), max(edge(p, c, d), edge(p, d, a))); }
// The fighter from behind and above. x: distance, y: 1 on the tails.
vec2 jet(vec2 q) {
  vec2 m = vec2(abs(q.x), q.y);
  float body = length(vec2(q.x, q.y - clamp(q.y, -0.45, 0.5))) - 0.05 * (1.0 - smoothstep(0.03, 0.5, q.y) * 0.92);
  float wing = quad(m, vec2(0.05, 0.14), vec2(0.05, -0.2), vec2(0.46, -0.22), vec2(0.46, -0.16));
  float tail = quad(m, vec2(0.04, -0.36), vec2(0.04, -0.48), vec2(0.2, -0.5), vec2(0.2, -0.45));
  float fin = quad(m, vec2(0.06, -0.3), vec2(0.07, -0.47), vec2(0.105, -0.44), vec2(0.14, -0.3));
  return vec2(min(min(body, wing), min(tail, fin)), step(fin, min(body, min(wing, tail)) + 0.004));
}
// One cloud dome from cell g + o, soft-blended into acc.
void dome(vec2 g, vec2 f, vec2 o, inout vec4 acc) {
  vec2 r = h22(g + o), c = f - o - 0.2 - 0.6 * r;
  float R = 0.6 + 0.4 * fract(r.x * 7.3 + r.y), up = sqrt(max(R * R - dot(c, c), 0.0)), k = up * up * up * (1.0 + 2.0 * fract(r.y * 5.1));
  acc += vec4(k * c * (0.6 + up), k * (up + 0.1), k * up);
}
// Value noise with its gradient.
vec3 noised(vec2 p) {
  vec2 i = floor(p), f = fract(p), u = f * f * (3.0 - 2.0 * f), du = 6.0 * f * (1.0 - f);
  float a = h2(i), b = h2(i + vec2(1.0, 0.0)), c = h2(i + vec2(0.0, 1.0)), d = h2(i + 1.0), k = a - b - c + d;
  return vec3(a + (b - a) * u.x + (c - a) * u.y + k * u.x * u.y, du * (vec2(b - a, c - a) + k * u.yx));
}
// Thunderhead density at q (> 0 inside): an outline with a crown and turret, cauliflower bulges of many sizes (fine: the
// small ones too), the foot thinning to wisps.
vec3 towerA, towerB; float horizon;
float cloud(vec2 q, float fine) {
  vec2 ya = (q.y - horizon) / vec2(towerA.z, towerB.z);
  vec2 xa = (q.x - vec2(towerA.x, towerB.x)) / (vec2(towerA.y, towerB.y) * (1.0 - 0.15 * clamp(ya, 0.0, 1.0)));
  vec2 top = max(ya - 0.62, 0.0) / 0.38, e = 1.0 - xa * xa - top * top - smoothstep(0.12, -0.04, ya) * 0.9 - step(ya, vec2(-0.04)) * 9.0;
  vec2 tu = (q - vec2(towerA.x + towerA.y * 0.7, horizon + towerA.z * 0.62)) / (towerA.y * vec2(0.6, 1.05));
  float env = max(max(e.x, e.y), 1.0 - dot(tu, tu));
  vec2 b = q * 13.0 + vec2(0.0, -iTime * 0.015);
  float n = abs(noise(b) * 2.0 - 1.0) * 0.45 + abs(noise(b * 2.3 + 3.1) * 2.0 - 1.0) * 0.28
          + (abs(noise(b * 5.1 + 1.7) * 2.0 - 1.0) * 0.15 + abs(noise(b * 10.7 + 5.3) * 2.0 - 1.0) * 0.06) * fine;
  return env + (n - 0.32 - 0.1 * (1.0 - fine)) * 1.5 - smoothstep(0.1, -0.4, env);
}
void mainImage(out vec4 color, in vec2 pixel) {
  vec2 uv = pixel / iResolution.xy; float ar = iResolution.x / iResolution.y; vec2 p = vec2(uv.x * ar, uv.y);
  float px = 1.5 / iResolution.y, t = iTime, wide = clamp((ar - 0.6) / 0.9, 0.0, 1.0);
  // Bank: w is the rolled point.
  float bank = 0.06 + 0.05 * sin(t * 0.13);
  mat2 roll = mat2(cos(bank), -sin(bank), sin(bank), cos(bank));
  vec2 cen = vec2(ar * 0.5, 0.5), w = roll * (p - cen) + cen;
  vec2 sunP = vec2(ar * mix(0.6, 0.64, wide), mix(0.84, 0.8, wide)), sun = roll * (sunP - cen) + cen;
  float xn = (w.x - cen.x) / max(ar, 1.0), hz = mix(0.54, 0.46, wide) - 0.06 * xn * xn;
  float r = length(w - sun), glow = exp(-abs(w.x - sun.x) * 1.3), sy = w.y - hz;
  // Sky: pale horizon, cerulean, cobalt.
  vec3 col = mix(mix(vec3(0.84, 0.91, 0.98), vec3(0.36, 0.6, 0.96), smoothstep(0.0, 0.2, sy)), vec3(0.1, 0.3, 0.78), smoothstep(0.24, 0.6, sy));
  vec3 haze = mix(vec3(0.84, 0.9, 0.98), vec3(1.0, 0.95, 0.86), glow * 0.7);
  // The deck: a plane below rushing toward us.
  if (sy < 0.0) {
    float v = -sy, z = 0.1 / (v + 0.004), fp = 0.8 / ((v + 0.004) * (v + 0.004)) * px;
    vec2 wp = vec2((w.x - cen.x) * z, z) * 8.0 + vec2(0.0, t * 1.8);
    vec3 m1 = noised(wp * 0.1 + 3.0), m2 = noised(wp * 0.26 + 7.0);  // big masses
    float mass = m1.x * 0.65 + m2.x * 0.35;
    vec2 mg = (m1.yz * 0.065 + m2.yz * 0.091) * 7.0;
    wp += vec2(noise(wp * 0.3), noise(wp * 0.3 + 5.2)) * 1.4;  // warp
    vec2 g = floor(wp), f = fract(wp); vec4 a = vec4(0.0);
    dome(g, f, vec2(-1.0), a); dome(g, f, vec2(0.0, -1.0), a); dome(g, f, vec2(1.0, -1.0), a);
    dome(g, f, vec2(-1.0, 0.0), a); dome(g, f, vec2(0.0), a); dome(g, f, vec2(1.0, 0.0), a);
    dome(g, f, vec2(-1.0, 1.0), a); dome(g, f, vec2(0.0, 1.0), a); dome(g, f, vec2(1.0), a);
    vec4 dm = vec4(a.w / (a.z + 0.001) * 1.6, normalize(a.xyz + vec3(0.0, 0.0, 0.001)));
    vec3 n = normalize(vec3(-mg * 1.3, 1.0) + (dm.yzw - vec3(0.0, 0.0, 1.0)) * 0.95);
    float diff = dot(n, vec3(0.3, 0.81, 0.51)) * 0.5 + 0.5, e = 1.0 - n.z, rim = e * e * max(dot(n.xy, vec2(0.35, 0.94)), 0.0);
    vec3 cc = vec3(0.52, 0.62, 0.84) * (0.7 + 0.3 * n.z) + vec3(1.0, 0.95, 0.86) * diff * diff * 0.45 + vec3(1.0, 0.86, 0.62) * rim * (0.3 + 0.4 * glow);
    cc *= (0.78 + 0.22 * smoothstep(0.2, 1.1, dm.x)) * (1.0 + (noise(wp * 3.1) - 0.5) * 0.14 * smoothstep(0.3, 0.05, fp));
    cc = mix(cc, vec3(0.3, 0.42, 0.64), smoothstep(0.35, 0.05, dm.x) * 0.5);
    cc *= 0.84 + 0.26 * smoothstep(0.3, 0.75, mass);
    // Gaps onto a hazier layer below where the masses thin.
    vec3 low = mix(vec3(0.3, 0.4, 0.6), haze, 0.3) * (0.92 + 0.16 * noise(wp * 0.7));
    cc = mix(low, cc, smoothstep(0.41, 0.5, mass + (dm.x - 0.5) * 0.15));
    cc = mix(cc, vec3(0.84, 0.89, 0.97), smoothstep(0.12, 0.45, fp));
    col = mix(cc, haze, 1.0 - exp(-z * 0.2));
    col = mix(col, haze, smoothstep(0.015, 0.0, v));
  }
  // Thunderheads, lit by how little cloud lies toward the sun; feet in haze.
  towerA = vec3(cen.x + 0.36 * ar, mix(0.11, 0.12, wide), 0.35); towerB = vec3(cen.x - 0.24 * ar, 0.06, 0.18); horizon = hz;
  if (sy > -0.05 && ((abs(w.x - towerA.x) < 0.3 && sy < 0.5) || (abs(w.x - towerB.x) < 0.14 && sy < 0.26))) {
    float D = cloud(w, 1.0), yy = sy / 0.4;
    vec2 L = normalize(sun - w);
    float D1 = cloud(w + L * 0.015, 0.0), T = exp(-max(D1, 0.0) * 1.3 - max(cloud(w + L * 0.06, 0.0), 0.0) * 0.9);
    float lobe = clamp((cloud(w, 0.0) - D1) * 3.0 + 0.4, 0.0, 1.0);
    vec3 tc = mix(mix(vec3(0.42, 0.49, 0.68), vec3(0.6, 0.66, 0.84), smoothstep(0.0, 0.7, yy)), vec3(1.0, 0.97, 0.92), max(T, lobe * 0.5));
    tc += vec3(1.0, 0.95, 0.85) * smoothstep(0.6, 1.0, yy) * T * 0.15;
    tc += vec3(1.0, 0.88, 0.66) * T * smoothstep(0.25, 0.0, D) * 0.35;
    tc = mix(tc, haze, smoothstep(0.3, 0.0, yy) * 0.75 + 0.1 + 0.3 * step(w.x, cen.x));
    col = mix(col, tc, smoothstep(0.0, mix(0.45, 0.07, smoothstep(0.0, 0.35, yy)), D));
  }
  // Sun: disc, glare, ice halo, streak.
  vec2 ds = w - sun; float rr = (r - 0.24) / 0.014;
  vec3 glare = vec3(1.0, 0.92, 0.74) * (exp(-r * 2.6) * 0.32 + exp(-r * 10.0) * 0.55 + exp(-r * 36.0) * 0.8);
  col += glare + mix(vec3(0.9, 0.55, 0.35), vec3(0.55, 0.75, 1.0), smoothstep(0.23, 0.255, r)) * exp(-rr * rr) * 0.12;
  col += vec3(1.0, 0.9, 0.75) * exp(-abs(ds.y) * 90.0) * exp(-abs(ds.x) * 1.8) * 0.28;
  col = mix(col, vec3(1.0, 0.99, 0.95), smoothstep(0.028 + px, 0.028 - px, r));
  // The fighter and its contrails.
  vec2 jc = vec2(ar * mix(0.5, 0.61, wide), mix(0.66, 0.62, wide) + 0.004 * sin(t * 0.8));
  float js = mix(0.13, 0.15, wide);
  vec2 q = (p - jc) / js; q.y /= 0.55;
  float back = max(-q.y - 0.47, 0.0), trail = abs(abs(q.x) - 0.035 - back * 0.11) / (px / js + back * 0.012);
  float lane = step(0.0, -q.y - 0.47) * exp(-back * 0.3);
  col = mix(col * (1.0 - 0.12 * exp(-trail * trail * 0.25) * lane), vec3(1.0), exp(-trail * trail) * lane * 0.95);
  if (abs(q.x) < 0.6 && abs(q.y) < 0.6) {
    vec2 jd = jet(q); float d = jd.x, aa = px / js * 1.4;
    vec3 hull = mix(vec3(0.2, 0.23, 0.28), vec3(0.4, 0.44, 0.5), smoothstep(-0.3, 0.4, q.x + q.y * 0.3));
    hull = mix(hull, mix(vec3(0.14, 0.16, 0.2), vec3(0.62, 0.64, 0.66), step(0.0, q.x)), jd.y);
    hull += vec3(1.0, 0.9, 0.7) * smoothstep(-0.02, 0.0, d) * 0.45 * step(0.0, q.y + 0.1);
    float can = length((q - vec2(0.0, 0.2)) * vec2(1.9, 0.75)) - 0.06;
    hull = mix(hull, vec3(0.08, 0.1, 0.16), smoothstep(0.01, -0.01, can));
    hull += vec3(1.0, 0.95, 0.85) * smoothstep(0.035, 0.0, length(q - vec2(0.015, 0.23))) * 0.9;
    col = mix(col, hull, smoothstep(aa, -aa, d));
  }
  color = vec4(col + glare * 0.2, 1.0);
}