// Starward, 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 n3(vec2 p) { return (0.5 * noise(p) + 0.25 * noise(p * 2.03 + 1.7) + 0.125 * noise(p * 4.07 + 3.3)) / 0.875; }
// Stars: one per cell at most, soft, denser in the galaxy's band.
vec3 starLayer(vec2 p, float cell, float chance, float px, float seed, vec2 band) {
  vec2 g = p / cell, id = floor(g), r = h22(id + seed) - 0.5;
  float inBand = exp(-pow(dot((id + 0.5) * cell, band) - 1.0, 2.0) / 0.025);
  float on = step(1.0 - chance * (1.0 + 2.0 * inBand), h2(id + seed * 1.7)), b = h2(id + seed + 9.1);
  float d = length((fract(g) - 0.5 - r * 0.7) * cell), s = px * (0.85 + 1.3 * b * b);
  vec3 tint = mix(vec3(1.0, 0.8, 0.6), vec3(0.68, 0.85, 1.0), h2(id + seed + 3.7));
  return tint * on * (0.2 + 0.85 * b * b) * exp(-d * d / (s * s));
}
// Ring density by radius (planet radii): bands, a gap, soft edges; w is a pixel.
float ringDensity(float rr, float w) {
  float k = smoothstep(1.32 - w, 1.4 + w, rr) * smoothstep(2.45 + w, 2.32 - w, rr);
  float gap = 1.0 - 0.85 * smoothstep(0.04 + w, 0.015, abs(rr - 1.98));
  return k * gap * (0.55 + 0.28 * sin(rr * 29.0) + 0.17 * sin(rr * 67.0 + 2.0)) * mix(0.55, 1.0, smoothstep(1.4, 1.75, rr));
}
// A diffraction spike along a.y, a pixel wide at the core.
float spike(vec2 a, float px, float fall) { return exp(-abs(a.x) / (px + 0.012 * abs(a.y))) * exp(-abs(a.y) * fall); }
void mainImage(out vec4 color, in vec2 pixel) {
  vec2 uv = pixel / iResolution.xy; float ar = iResolution.x / iResolution.y, px = 1.5 / iResolution.y, t = iTime * 0.02;
  vec2 p = vec2(uv.x * ar, uv.y); float wide = smoothstep(0.7, 1.4, ar);
  // The giant: right of centre when wide, high on a phone; its sun off the upper-left limb.
  vec2 pc = mix(vec2(ar * 0.62, 0.64), vec2(ar * 0.7, 0.4), wide); float R = mix(0.25, 0.38, wide);
  vec3 L = normalize(vec3(-0.6, 0.55, -0.28)), N = normalize(vec3(-0.2, 1.0, 0.24));
  vec2 sunDir = normalize(L.xy), sun = pc + sunDir * R * mix(1.2, 1.32, wide);
  vec2 sd = p - sun; float sr = length(sd), glow = 1.0 / (1.0 + sr * 3.5);
  // Nebula: teal, rose, gold toward the sun; filaments, dark dust, the galaxy's band.
  vec2 band = normalize(vec2(0.32, 1.0)); band /= dot(vec2(ar * 0.5, mix(0.88, 0.8, wide)), band);
  vec2 q = p * 1.25 + vec2(t, -0.3 * t);
  float w = n3(q * 1.6 + 5.0), neb = fbm(q + w * 1.5 - t * 0.3), hg = noise(q * 0.8 + 1.0) + 0.4 * glow;
  vec3 col = vec3(0.004, 0.012, 0.03) * (1.0 + 1.5 * uv.y);
  vec3 nc = mix(vec3(0.0, 0.42, 0.56), vec3(0.85, 0.16, 0.5), smoothstep(0.36, 0.5, hg));
  nc = mix(nc, vec3(1.0, 0.56, 0.16), smoothstep(0.58, 0.75, hg));
  col += nc * (pow(neb, 2.6) * (1.0 + 1.4 * glow) + pow(1.0 - abs(neb * 2.0 - 1.1), 8.0) * 0.5);
  float bandD = dot(p, band) - 1.0;
  col += vec3(0.5, 0.47, 0.46) * exp(-bandD * bandD / 0.01) * (0.25 + 0.75 * w) * (1.0 - 0.6 * exp(-pow(bandD / 0.025 + 2.0 * w - 1.0, 2.0)));
  col *= 1.0 - 0.7 * smoothstep(0.58, 0.72, w) * (1.0 - smoothstep(0.75, 0.9, neb));
  col += vec3(1.0, 0.75, 0.45) * pow(neb, 5.0) * smoothstep(0.55, 0.72, w) * 3.0 * glow;
  col += starLayer(p, 0.021, 0.2, px, 1.0, band) * 0.75 + starLayer(p, 0.07, 0.2, px * 1.3, 5.0, band);
  // Rays from the sun through the dust, slowly turning.
  float ang = atan(sd.y, sd.x), rays = noise(vec2(ang * 9.0 + t * 2.0, 0.0)) * noise(vec2(ang * 23.0 - t * 3.0, 3.0));
  col += vec3(1.0, 0.8, 0.55) * rays * (0.08 + 0.3 * neb) * exp(-sr * 2.6) * smoothstep(0.0, 0.1, sr);
  // Planet and ring plane in planet radii.
  vec2 d = (p - pc) / R; float r2 = dot(d, d), r = sqrt(r2), pxR = px / R;
  float zs = sqrt(max(0.0, 1.0 - r2)), zq = -(N.x * d.x + N.y * d.y) / N.z, rq = length(vec3(d, zq));
  float toward = max(dot(d / max(r, 0.001), sunDir), 0.0);
  float atm = exp(-max(r - 1.0, 3.0 - 3.0 * r) * R * mix(80.0, 35.0, toward)) * (0.1 + 1.4 * pow(toward, 6.0) + 0.4 * toward * toward);
  if (r < 1.0 + pxR) {
    vec3 P = vec3(d, zs), T = normalize(vec3(-N.y, N.x, 0.0));
    float lat = dot(P, N), lon = atan(dot(P, T), dot(P, cross(N, T))) + iTime * 0.006;
    // Bands warped by slow turbulence, calm at the poles.
    float tur = n3(vec2(lon * 3.0, lat * 14.0)), la = lat + 0.05 * (tur - 0.5) * smoothstep(0.97, 0.85, abs(lat));
    float bv = 0.5 + 0.3 * sin(la * 23.0) + 0.14 * sin(la * 57.0 + 1.3) + 0.08 * sin(la * 9.0 + 2.0);
    vec3 pcol = mix(vec3(0.33, 0.11, 0.05), vec3(0.98, 0.8, 0.55), clamp(bv, 0.0, 1.0));
    pcol = mix(pcol, vec3(0.78, 0.4, 0.18), smoothstep(0.62, 0.8, tur) * 0.6);
    pcol = mix(pcol, vec3(0.45, 0.58, 0.72), smoothstep(0.6, 0.92, abs(lat)) * 0.6);
    // Light: soft terminator with a red edge, the ring's shadow, a hot rim toward the sun.
    float dl = dot(P, L), lit = smoothstep(-0.12, 0.5, dl), s = -dot(P, N) / dot(L, N);
    float shade = 1.0 - 0.75 * step(0.0, s) * ringDensity(length(P + s * L), 0.02);
    vec3 pl = pcol * (lit * shade * 1.25 + 0.012) + vec3(0.13, 0.03, 0.01) * smoothstep(0.14, 0.0, abs(dl + 0.02));
    pl += vec3(0.6, 0.3, 0.15) * pow(1.0 - zs, 3.0) * pow(toward, 3.0);
    col = mix(col, pl, smoothstep(1.0 + pxR, 1.0 - pxR, r));
  }
  col += vec3(1.0, 0.72, 0.45) * atm + vec3(0.2, 0.33, 0.5) * atm * (1.0 - toward);
  // Rings in front of or beside the planet: glowing toward the sun, cut by its shadow.
  if (r > 1.0 - pxR || zq > zs) {
    vec3 Q = vec3(d, zq); float b = dot(Q, L), disc = b * b - dot(Q, Q) + 1.0;
    float shadow = 1.0 - 0.8 * smoothstep(-0.02, 0.06, disc) * step(b, 0.0);
    float fwd = 0.3 + 1.2 * pow(dot(normalize(Q.xy), sunDir) * 0.5 + 0.5, 3.0);
    vec3 rc = mix(vec3(0.95, 0.72, 0.5), vec3(1.0, 0.92, 0.8), 0.5 + 0.5 * sin(rq * 13.0));
    float front = (zq > zs) ? 1.0 : smoothstep(1.0 - pxR, 1.0 + pxR, r);
    col = mix(col, col * 0.35 + rc * fwd * shadow * 0.8, clamp(ringDensity(rq, pxR / N.z) * 0.85, 0.0, 1.0) * front);
  }
  // Two moons in crescent.
  for (int i = 0; i < 2; i++) {
    float fi = float(i), mr = R * mix(0.12, 0.05, fi);
    vec2 md = (p - pc - R * mix(vec2(-2.3, -0.75), vec2(1.15, 1.25), fi) * mix(vec2(0.42, 1.25), vec2(1.0), wide)) / mr;
    float m2 = dot(md, md);
    if (m2 < 1.5) {
      vec3 M = vec3(md, sqrt(max(0.0, 1.0 - m2)));
      vec3 mm = vec3(0.68, 0.58, 0.5) * (0.7 + 0.4 * noise(md * 3.0 + fi * 7.0));
      mm = mm * (smoothstep(-0.05, 0.45, dot(M, L)) * 1.2 + 0.015) + vec3(0.6, 0.42, 0.27) * pow(1.0 - M.z, 3.0) * max(dot(md, sunDir), 0.0);
      col = mix(col, mm, smoothstep(1.0 + px / mr, 1.0 - px / mr, sqrt(m2)));
    }
  }
  // The sun: a hot core, spikes, a blue lens streak, a warm bloom.
  vec2 sx = vec2(sd.x + sd.y, sd.y - sd.x) * 0.7071;
  float spikes = spike(sd, px, 7.0) + spike(sd.yx, px, 7.0) + 0.5 * (spike(sx, px, 12.0) + spike(sx.yx, px, 12.0));
  col += vec3(1.0, 0.95, 0.85) * (exp(-sr * sr / 0.00025) * 2.2 + spikes * 0.5) + vec3(0.035, 0.021, 0.01) / (sr * 6.0 + 0.2);
  col += vec3(0.33, 0.45, 0.6) * exp(-abs(sd.y) / (px * 1.2 + 0.0015)) * exp(-abs(sd.x) * 1.8);
  // Grade: soft shoulder, vignette, lift toward deep blue.
  col = col / (1.0 + col * 0.4);
  col *= 1.0 - 0.35 * pow(length(uv - vec2(0.55, 0.5)) * 1.1, 2.5);
  color = vec4(pow(col, vec3(0.9)) + vec3(0.004, 0.008, 0.016), 1.0);
}