Ocean

A near gale at dusk on an open sea. Four 512² FFT cascades of a JONSWAP wind sea and swell are evolved, inverse-transformed and packed into mipmapped textures every frame. A clipmap of eleven levels displaces the surface, and the fragment shader filters slopes and their moments, so distant waves turn into roughness instead of sparkle. The water's body is shaded after the Atlas water model as Ethan Truong's GodotOceanWaves (MIT) implements it: a dark albedo lit by the sun and sky, and the sun scattered through the backs of raised waves. Foam forms where the horizontal deformation compresses the crests, and a persistent buffer keeps it as lace behind them. The sky comes from Hillaire's atmosphere lookup tables and a volumetric cloud deck, which an environment bakes into the water's reflections, and bloom and the final pass finish the frame. The eye stands still above the swell; the sun runs from below the horizon to 60°.

Files, in reading order

  1. scene.tsCreates the scene: its resources, its settings and what each frame records.
  2. atmosphere.tsA planet's atmosphere after Hillaire's 2020 lookup tables: transmittance, multiple scattering and the sky seen from the eye.
  3. clouds.tsA stratocumulus deck over a 40 km box round the eye, tiled so it drifts without a new bake.
  4. foam.tsWhitecap lace: a tile of torn foam that layer cuts to the share of the water the breaking waves call for, so a fresh cap closes over and an older one erodes into lace along the wind.
  5. spectrum.tsThe sea's spectrum under a wind of speed metres per second, blowing towards direction.
  6. surface.tsThe sea's surface: a clipmap displaced by the cascades, shaded with filtered slope covariance, the water body, sun glint, foam and mist, then the sky where the sea leaves the depth clear.
  7. waves.tsFour FFT cascades of an ocean spectrum, mipmapped as displacement, slope and moment maps.
  8. shading.ts
  9. wave-mips.tsBox-filtered mips of the wave maps, three levels a dispatch.

Kit modules bloom, cameras, clock, environment, exposure, finish, history, memo, mipmaps, owner, pipelines, targets, tone. Also @typegpu/noise, typegpu, wgpu-matrix.

scene.ts

import {  bloom,  cameras,  clock,  environment,  exposure,  finish,  owner,  pipelines,  targets,  tone,  type Frame,  type Surface,} from '@typegpu/kit';import { d, type TgpuRoot } from 'typegpu';import { mat4 } from 'wgpu-matrix';import type { Settings } from '../settings.ts';import { createAtmosphere } from './atmosphere.ts';import { createClouds } from './clouds.ts';import { createFoam } from './foam.ts';import { foamingAt, spectrumAt } from './spectrum.ts';import { createSurface } from './surface.ts';import { createWaves } from './waves.ts';const HDR = 'rgba16float';const DEPTH = 'depth32float';const UP = [0, 1, 0] as const;/** The waves travel towards +z, at the camera, out of the low sun at −z. */const WIND_DIRECTION = Math.PI / 2;/** The sun sits 15° right of the view's centre, so its glitter runs across the sea. */const SUN_AZIMUTH = 0.26;/** * The eye stands `height` metres above the still sea, or `clear` above the highest crest beneath * it. Over the sea's loop that crest stays below 0.78 of the significant wave height at any wind, * and is `crest` of it at the default. Angles are in radians. */const EYE = {  height: 4,  clear: 1.5,  crest: 0.7312573,  pitch: 0.028,  pitchRange: [-0.35, 0.6] as const,  across: (55 * Math.PI) / 180,  tallest: 1.1,  turning: 14,  near: 0.5,  far: 20_000,};const ENVIRONMENT_SIZE = 128;const ENVIRONMENT_TEXEL = Math.PI / (2 * ENVIRONMENT_SIZE);/** The deck's drift in metres per second, and the seconds between the sky's rebakes. */const DRIFT = { speed: 14, sky: 2 };const EXPOSURE = { key: 0.15, lowest: 1, highest: 3000 };/** The key darkens to `share` as the sun sets through `fade` degrees, so dusk stays dark. */const DUSK = { share: 0.25, fade: [-4, 10] as const };/** The sea opens `START` seconds into its loop, after `HISTORY.steps` steps of foam history. */const START = 195.4;const HISTORY = { steps: 24, step: 0.125 };/** Seconds a frame may advance the sea, so a stalled tab resumes calmly. */const LONGEST_STEP = 0.1;/** * A near gale at dusk: FFT waves with foam where the crests compress, and a low sun under a cloud * deck lighting the backs of the raised waves and laying a glint path to the eye. */export async function createScene(root: TgpuRoot, surface: Surface) {  const own = owner();  const values = { wind: 20, sun: 0.5, clouds: 1 };  const screen = own(    targets.create(root, surface, {      hdr: { format: HDR, sampled: true },      depth: { format: DEPTH },    }),  );  const view = own(cameras.custom(root));  // Hosts turn the goal as an orbit's; the eye stays put and looks the other way, at the target  const goal = { yaw: 0, pitch: EYE.pitch, distance: 1 };  const look = { ...goal };  const waves = own(createWaves(root));  const foam = own(createFoam(root, { wind: WIND_DIRECTION }));  const atmosphere = own(createAtmosphere(root, { haze: 3.5 }));  const clouds = own(createClouds(root, atmosphere, view, surface, { coverage: 1 }));  // Reflections and ambient light see the sky and deck from the still sea, metres below the eye  const light = own(environment.create(root, { size: ENVIRONMENT_SIZE }));  const skySource = light.add((direction: d.v3f) => {    'use gpu';    const deck = clouds.march(d.vec3f(), direction, ENVIRONMENT_TEXEL);    return atmosphere.skyRadiance(direction) * deck.w + deck.xyz;  });  const sea = own(    createSurface(root, {      view,      waves,      atmosphere,      clouds,      foam,      light,      format: HDR,      depth: DEPTH,    }),  );  const glow = own(    bloom.create(root, {      source: () => screen.textures.hdr,      settings: { threshold: 3, knee: 0.6 },      levels: 3,    }),  );  const metering = own(    exposure.create(root, { source: () => screen.textures.hdr, settings: EXPOSURE }),  );  const final = own(    finish.create(root, {      format: surface.format,      source: () => screen.textures.hdr,      glow,      exposure: metering,      tone: tone.agx,      settings: { exposure: 1, vignette: 0.3, glowIntensity: 0.08 },    }),  );  await Promise.all([    pipelines.compile([      ...waves.pipelines,      ...foam.pipelines,      ...sea.pipelines,      ...atmosphere.pipelines,      ...clouds.pipelines,      glow,      metering,      final,    ]),    light.initAsync(),  ]);  let rebuild = true;  let rebake = true;  let skyBaked = -Infinity;  const eye = new Float32Array(3);  const ahead = new Float32Array(3);  const viewMatrix = mat4.create();  const projection = mat4.create();  const write = (key?: keyof typeof values) => {    if (key === undefined || key === 'wind') {      waves.foaming.write(foamingAt(values.wind));      const spectrum = spectrumAt(values.wind, WIND_DIRECTION);      waves.spectrum.write(spectrum);      eye[1] = Math.max(EYE.height, EYE.crest * spectrum.waveHeight + EYE.clear);      // Changes between frames coalesce: the next frame rebuilds the sea once      rebuild = true;    }    if (key === undefined || key === 'sun') {      const elevation = (values.sun * Math.PI) / 180;      const sun = d.vec3f(        Math.cos(elevation) * Math.sin(SUN_AZIMUTH),        Math.sin(elevation),        -Math.cos(elevation) * Math.cos(SUN_AZIMUTH),      );      atmosphere.place(sun, EYE.height);      const [from, to] = DUSK.fade;      const day = Math.min(Math.max((values.sun - from) / (to - from), 0), 1);      metering.settings.patch({        key: EXPOSURE.key * (DUSK.share + (1 - DUSK.share) * day * day),      });      rebake = true;    }    if (key === undefined || key === 'clouds') {      clouds.cover(values.clouds);      rebake = true;    }  };  write();  sea.place(eye[0]!, eye[2]!);  const settings: Settings<typeof values> = {    values,    knobs: {      wind: { label: 'wind (m/s)', min: 4, max: 24, step: 0.5 },      sun: { label: 'sun height (°)', min: -4, max: 60, step: 0.5 },      clouds: { label: 'cloud cover', min: 0, max: 1, step: 0.05 },    },    write,  };  // The clouds' first bake also fills their noise volume, longer than a frame, so it runs at load  const warm = root['~unstable'].createCommandEncoder();  const warming = warm.beginComputePass();  clouds.bake(warming);  foam.bake(warming);  waves.rebuild(warming);  warming.end();  warm.submit();  for (let step = 0; step < HISTORY.steps; step++) {    waves.step.patch({ time: START - (HISTORY.steps - step) * HISTORY.step, dt: HISTORY.step });    const encoder = root['~unstable'].createCommandEncoder();    const history = encoder.beginComputePass();    waves.record(history);    history.end();    encoder.submit();  }  rebuild = false;  const aim = () => {    goal.pitch = Math.min(Math.max(goal.pitch, EYE.pitchRange[0]), EYE.pitchRange[1]);    const level = Math.cos(look.pitch);    ahead[0] = -level * Math.sin(look.yaw);    ahead[1] = -Math.sin(look.pitch);    ahead[2] = -level * Math.cos(look.yaw);    mat4.lookAt(      eye,      [eye[0]! + ahead[0]!, eye[1]! + ahead[1]!, eye[2]! + ahead[2]!],      UP,      viewMatrix,    );    const aspect = surface.width / surface.height;    const tall = Math.min(2 * Math.atan(Math.tan(EYE.across / 2) / aspect), EYE.tallest);    mat4.perspective(tall, aspect, EYE.near, EYE.far, projection);    view.write(viewMatrix, projection);  };  return {    view: {      goal,      /** Turns the eye to its goal at once. */      jump() {        Object.assign(look, goal);      },    },    settings,    frame(frame: Frame) {      const turning = clock.ease(frame.dt, EYE.turning);      look.yaw += (goal.yaw - look.yaw) * turning;      look.pitch += (goal.pitch - look.pitch) * turning;      aim();      screen.fit();      const dt = Math.min(frame.dt, LONGEST_STEP);      waves.step.patch({ time: START + frame.time, dt });      const drifted = frame.time * DRIFT.speed;      clouds.move([Math.cos(WIND_DIRECTION) * drifted, Math.sin(WIND_DIRECTION) * drifted]);      const bakeSky = rebake || frame.time - skyBaked >= DRIFT.sky;      const compute = frame.encoder.beginComputePass();      if (rebake) {        atmosphere.bake(compute);        clouds.bake(compute);      }      clouds.record(compute);      if (rebuild) waves.rebuild(compute);      waves.record(compute);      compute.end();      rebuild = false;      if (bakeSky) {        light.bake(frame.encoder, skySource);        skyBaked = frame.time;      }      rebake = false;      const pass = frame.encoder.beginRenderPass({        colorAttachments: [{ view: screen.views.hdr, loadOp: 'clear', storeOp: 'store' }],        depthStencilAttachment: {          view: screen.views.depth,          depthClearValue: 1,          depthLoadOp: 'clear',          depthStoreOp: 'discard',        },      });      sea.draw(pass);      pass.end();      glow.record(frame);      metering.record(frame);      final.record(frame, surface.view(frame));    },    destroy: own.destroy,  };}

atmosphere.ts

import { owner } from '@typegpu/kit';import { d, std, tgpu, type TgpuComputePass, type TgpuRoot } from 'typegpu';/** Planet and atmosphere radii in kilometres. */const GROUND_RADIUS = 6360;const TOP_RADIUS = 6460;const ATMOSPHERE_HEIGHT = TOP_RADIUS - GROUND_RADIUS;/** Scattering and extinction per kilometre at sea level, and scale heights in kilometres. */const RAYLEIGH_SCATTERING = d.vec3f(5.802e-3, 13.558e-3, 33.1e-3);const RAYLEIGH_HEIGHT = 8;const MIE_SCATTERING = 3.996e-3;const MIE_EXTINCTION = 4.44e-3;const MIE_HEIGHT = 1.2;const MIE_G = 0.8;const OZONE_EXTINCTION = d.vec3f(0.65e-3, 1.881e-3, 0.085e-3);const OZONE_CENTER = 25;const OZONE_WIDTH = 15;const TRANSMITTANCE_SIZE = d.vec2f(256, 64);const MULTI_SCATTER_SIZE = 32;const SKY_VIEW_SIZE = d.vec2f(192, 108);const TRANSMITTANCE_STEPS = 40;const MULTI_SCATTER_STEPS = 20;const SPHERE_SAMPLES = 8;const SKY_VIEW_STEPS = 32;const TILE = 8;const SUN_ANGULAR_RADIUS = 0.004675;const SUN_SOLID_ANGLE = Math.PI * SUN_ANGULAR_RADIUS ** 2;const SUN_COS_EDGE = Math.cos(SUN_ANGULAR_RADIUS);const Settings = d.struct({  sunDirection: d.vec3f,  haze: d.f32,  /** The eye's distance from the planet's centre in kilometres. */  eyeRadius: d.f32,});const Medium = d.struct({ rayleigh: d.vec3f, mie: d.vec3f, extinction: d.vec3f });const mediumAt = (r: number, haze: number) => {  'use gpu';  const height = std.max(r - GROUND_RADIUS, 0);  const rayleigh = RAYLEIGH_SCATTERING * std.exp(-height / RAYLEIGH_HEIGHT);  const mieDensity = std.exp(-height / MIE_HEIGHT) * haze;  const ozone = OZONE_EXTINCTION * std.max(1 - std.abs(height - OZONE_CENTER) / OZONE_WIDTH, 0);  return Medium({    rayleigh,    mie: d.vec3f(MIE_SCATTERING * mieDensity),    extinction: rayleigh + MIE_EXTINCTION * mieDensity + ozone,  });};/** Distance to the nearest intersection in front of the origin, negative when the ray misses. */const intersectSphere = (origin: d.v3f, dir: d.v3f, radius: number) => {  'use gpu';  const b = std.dot(origin, dir);  const c = std.dot(origin, origin) - d.f32(radius * radius);  const disc = b * b - c;  if (disc < 0) {    return d.f32(-1);  }  const near = -b - std.sqrt(disc);  if (near > 0) {    return near;  }  return -b + std.sqrt(disc);};const rayLength = (origin: d.v3f, dir: d.v3f) => {  'use gpu';  const ground = intersectSphere(origin, dir, GROUND_RADIUS);  if (ground > 0) {    return ground;  }  return std.max(intersectSphere(origin, dir, TOP_RADIUS), 0);};const subUv = (uv: d.v2f, size: d.v2f) => {  'use gpu';  return (uv + 0.5 / size) * (size / (size + 1));};const fromSubUv = (uv: d.v2f, size: d.v2f) => {  'use gpu';  return (uv - 0.5 / size) * (size / (size - 1));};const transmittanceUv = (r: number, mu: number) => {  'use gpu';  const h = std.sqrt(TOP_RADIUS * TOP_RADIUS - GROUND_RADIUS * GROUND_RADIUS);  const rho = std.sqrt(std.max(r * r - GROUND_RADIUS * GROUND_RADIUS, 0));  const disc = r * r * (mu * mu - 1) + TOP_RADIUS * TOP_RADIUS;  const dist = std.max(-r * mu + std.sqrt(std.max(disc, 0)), 0);  const dMin = TOP_RADIUS - r;  const dMax = rho + h;  return d.vec2f((dist - dMin) / (dMax - dMin), rho / h);};/** Angle between the nadir and the horizon seen from radius r. */const horizonBeta = (r: number) => {  'use gpu';  return std.acos(std.sqrt(std.max(r * r - GROUND_RADIUS * GROUND_RADIUS, 0)) / r);};const rayleighPhase = (cosTheta: number) => {  'use gpu';  return (3 / (16 * Math.PI)) * (1 + cosTheta * cosTheta);};const miePhase = (cosTheta: number) => {  'use gpu';  const g2 = MIE_G * MIE_G;  return (    ((3 / (8 * Math.PI)) * ((1 - g2) * (1 + cosTheta * cosTheta))) /    ((2 + g2) * std.pow(1 + g2 - 2 * MIE_G * cosTheta, 1.5))  );};export interface AtmosphereOptions {  /** Multiplies the aerosol density: 1 is a clear day, 3 a hazy one. Default 1. */  readonly haze?: number;}export type Atmosphere = ReturnType<typeof createAtmosphere>;/** * A planet's atmosphere after Hillaire's 2020 lookup tables: transmittance, multiple scattering * and the sky seen from the eye. Radiance is per unit of sunlight above the air. */export function createAtmosphere(root: TgpuRoot, options: AtmosphereOptions = {}) {  const own = owner();  const settings = root.createUniform(Settings, {    sunDirection: d.vec3f(0, 1, 0),    haze: options.haze ?? 1,    eyeRadius: GROUND_RADIUS,  });  const lutSampler = root.createSampler({    magFilter: 'linear',    minFilter: 'linear',    addressModeU: 'clamp-to-edge',    addressModeV: 'clamp-to-edge',  });  const createLut = (width: number, height: number) => {    const texture = own(      root        .createTexture({ size: [width, height], format: 'rgba16float' })        .$usage('sampled', 'storage'),    );    return {      target: texture.createView(d.textureStorage2d('rgba16float', 'write-only')),      view: texture.createView(d.texture2d(d.f32)),    };  };  const transmittanceLut = createLut(TRANSMITTANCE_SIZE.x, TRANSMITTANCE_SIZE.y);  const multiScatterLut = createLut(MULTI_SCATTER_SIZE, MULTI_SCATTER_SIZE);  const skyViewLut = createLut(SKY_VIEW_SIZE.x, SKY_VIEW_SIZE.y);  const sunlight = root.createMutable(d.vec3f);  const transmittance = (r: number, mu: number) => {    'use gpu';    const uv = transmittanceUv(r, mu);    return std.textureSampleLevel(transmittanceLut.view.$, lutSampler.$, uv, 0).xyz;  };  const sunlightAt = (pos: d.v3f, sun: d.v3f) => {    'use gpu';    let light = d.vec3f();    if (intersectSphere(pos, sun, GROUND_RADIUS) < 0) {      light = transmittance(std.length(pos), std.dot(pos, sun) / std.length(pos));    }    return light;  };  const bakeTransmittance = tgpu.computeFn({    workgroupSize: [TILE, TILE],    in: { gid: d.builtin.globalInvocationId },  })(({ gid }) => {    'use gpu';    if (gid.x >= TRANSMITTANCE_SIZE.x || gid.y >= TRANSMITTANCE_SIZE.y) return;    const uv = (d.vec2f(gid.xy) + 0.5) / TRANSMITTANCE_SIZE;    const h = std.sqrt(TOP_RADIUS * TOP_RADIUS - GROUND_RADIUS * GROUND_RADIUS);    const rho = h * uv.y;    const r = std.sqrt(rho * rho + GROUND_RADIUS * GROUND_RADIUS);    const dMin = TOP_RADIUS - r;    const dist = dMin + uv.x * (rho + h - dMin);    let mu = d.f32(1);    if (dist > 0) {      mu = std.clamp((h * h - rho * rho - dist * dist) / (2 * r * dist), -1, 1);    }    const origin = d.vec3f(0, r, 0);    const dir = d.vec3f(std.sqrt(1 - mu * mu), mu, 0);    const step = rayLength(origin, dir) / TRANSMITTANCE_STEPS;    let depth = d.vec3f();    for (let i = 0; i < TRANSMITTANCE_STEPS; i++) {      const pos = origin + dir * ((d.f32(i) + 0.5) * step);      depth += mediumAt(std.length(pos), settings.$.haze).extinction * step;    }    std.textureStore(transmittanceLut.target.$, gid.xy, d.vec4f(std.exp(depth * -1), 1));  });  /** Second and higher scattering orders summed as a geometric series. */  const bakeMultiScatter = tgpu.computeFn({    workgroupSize: [TILE, TILE],    in: { gid: d.builtin.globalInvocationId },  })(({ gid }) => {    'use gpu';    if (gid.x >= MULTI_SCATTER_SIZE || gid.y >= MULTI_SCATTER_SIZE) return;    const uv = (d.vec2f(gid.xy) + 0.5) / MULTI_SCATTER_SIZE;    const sunCos = uv.x * 2 - 1;    const sun = d.vec3f(std.sqrt(1 - sunCos * sunCos), sunCos, 0);    const origin = d.vec3f(0, GROUND_RADIUS + uv.y * ATMOSPHERE_HEIGHT, 0);    let luminance = d.vec3f();    let scattered = d.vec3f();    for (let a = 0; a < SPHERE_SAMPLES; a++) {      for (let b = 0; b < SPHERE_SAMPLES; b++) {        const theta = (Math.PI * (d.f32(a) + 0.5)) / SPHERE_SAMPLES;        const phi = (2 * Math.PI * (d.f32(b) + 0.5)) / SPHERE_SAMPLES;        const dir = d.vec3f(          std.sin(theta) * std.cos(phi),          std.cos(theta),          std.sin(theta) * std.sin(phi),        );        const step = rayLength(origin, dir) / MULTI_SCATTER_STEPS;        let throughput = d.vec3f(1);        for (let j = 0; j < MULTI_SCATTER_STEPS; j++) {          const pos = origin + dir * ((d.f32(j) + 0.5) * step);          const medium = mediumAt(std.length(pos), settings.$.haze);          const segment = std.exp(medium.extinction * -step);          const scattering = medium.rayleigh + medium.mie;          const integrated = (1 - segment) / std.max(medium.extinction, d.vec3f(1e-7));          luminance +=            throughput * scattering * sunlightAt(pos, sun) * (integrated / (4 * Math.PI));          scattered += throughput * scattering * integrated;          throughput *= segment;        }      }    }    const samples = SPHERE_SAMPLES * SPHERE_SAMPLES;    const psi = luminance / samples / (1 - scattered / samples);    std.textureStore(multiScatterLut.target.$, gid.xy, d.vec4f(psi, 1));  });  /** Sky radiance seen from the eye, by view zenith and azimuth from the sun. */  const bakeSkyView = tgpu.computeFn({    workgroupSize: [TILE, TILE],    in: { gid: d.builtin.globalInvocationId },  })(({ gid }) => {    'use gpu';    if (gid.x >= SKY_VIEW_SIZE.x || gid.y >= SKY_VIEW_SIZE.y) return;    const uv = fromSubUv((d.vec2f(gid.xy) + 0.5) / SKY_VIEW_SIZE, SKY_VIEW_SIZE);    const r = settings.$.eyeRadius;    const beta = horizonBeta(r);    const horizonAngle = Math.PI - beta;    let zenithAngle = d.f32(0);    if (uv.y < 0.5) {      const coord = 1 - 2 * uv.y;      zenithAngle = horizonAngle * (1 - coord * coord);    } else {      const coord = uv.y * 2 - 1;      zenithAngle = horizonAngle + beta * coord * coord;    }    const azimuthCos = -(uv.x * uv.x * 2 - 1);    const azimuthSin = std.sqrt(std.max(1 - azimuthCos * azimuthCos, 0));    const dir = d.vec3f(      std.sin(zenithAngle) * azimuthCos,      std.cos(zenithAngle),      std.sin(zenithAngle) * azimuthSin,    );    const sunCos = settings.$.sunDirection.y;    const sun = d.vec3f(std.sqrt(std.max(1 - sunCos * sunCos, 0)), sunCos, 0);    const origin = d.vec3f(0, r, 0);    const cosTheta = std.dot(dir, sun);    const phaseR = rayleighPhase(cosTheta);    const phaseM = miePhase(cosTheta);    const step = rayLength(origin, dir) / SKY_VIEW_STEPS;    let luminance = d.vec3f();    let throughput = d.vec3f(1);    for (let i = 0; i < SKY_VIEW_STEPS; i++) {      const pos = origin + dir * ((d.f32(i) + 0.5) * step);      const height = std.length(pos);      const medium = mediumAt(height, settings.$.haze);      const segment = std.exp(medium.extinction * -step);      const multiUv = d.vec2f(        (std.dot(pos, sun) / height) * 0.5 + 0.5,        (height - GROUND_RADIUS) / ATMOSPHERE_HEIGHT,      );      const psi = std.textureSampleLevel(multiScatterLut.view.$, lutSampler.$, multiUv, 0).xyz;      const source =        sunlightAt(pos, sun) * (medium.rayleigh * phaseR + medium.mie * phaseM) +        psi * (medium.rayleigh + medium.mie);      const integrated = (1 - segment) / std.max(medium.extinction, d.vec3f(1e-7));      luminance += throughput * source * integrated;      throughput *= segment;    }    std.textureStore(skyViewLut.target.$, gid.xy, d.vec4f(luminance, 1));  });  /** Sky radiance per unit of sunlight in a world direction, from the eye. */  const skyRadiance = (dir: d.v3f) => {    'use gpu';    const r = settings.$.eyeRadius;    const beta = horizonBeta(r);    const horizonAngle = Math.PI - beta;    const zenithAngle = std.acos(std.clamp(dir.y, -1, 1));    let v = d.f32(0);    if (zenithAngle < horizonAngle) {      v = (1 - std.sqrt(std.saturate(1 - zenithAngle / horizonAngle))) * 0.5;    } else {      v = std.sqrt(std.saturate((zenithAngle - horizonAngle) / beta)) * 0.5 + 0.5;    }    const flat = dir.xz / std.max(std.length(dir.xz), 1e-4);    const sun = settings.$.sunDirection;    const sunFlat = sun.xz / std.max(std.length(sun.xz), 1e-4);    const u = std.sqrt(std.saturate(-std.dot(flat, sunFlat) * 0.5 + 0.5));    const uv = subUv(d.vec2f(u, v), SKY_VIEW_SIZE);    return std.textureSampleLevel(skyViewLut.view.$, lutSampler.$, uv, 0).xyz;  };  const eyeTransmittance = (dir: d.v3f) => {    'use gpu';    const origin = d.vec3f(0, settings.$.eyeRadius, 0);    if (intersectSphere(origin, dir, GROUND_RADIUS) > 0) {      return d.vec3f();    }    return transmittance(settings.$.eyeRadius, dir.y);  };  const bakeLighting = tgpu.computeFn({ workgroupSize: [1] })(() => {    'use gpu';    sunlight.$ = eyeTransmittance(settings.$.sunDirection);  });  const transmittancePipeline = root.createComputePipeline({ compute: bakeTransmittance });  const multiScatterPipeline = root.createComputePipeline({ compute: bakeMultiScatter });  const skyViewPipeline = root.createComputePipeline({ compute: bakeSkyView });  const lightingPipeline = root.createComputePipeline({ compute: bakeLighting });  /** Like `aerialPerspective`, for cloud hiding `opacity` of the sky: haze adds only there. */  const veiled = (color: d.v3f, opacity: number, dir: d.v3f, distance: number) => {    'use gpu';    const extinction = mediumAt(d.f32(GROUND_RADIUS), settings.$.haze).extinction;    const fog = std.exp(extinction * (-distance / 1000));    const horizon = skyRadiance(std.normalize(d.vec3f(dir.x, 0.002, dir.z)));    return color * fog + horizon * ((1 - fog) * opacity);  };  return {    pipelines: [transmittancePipeline, multiScatterPipeline, skyViewPipeline, lightingPipeline],    skyRadiance,    /** A unit vector towards the sun. */    sunDirection: () => {      'use gpu';      return d.vec3f(settings.$.sunDirection);    },    /** The share of sunlight reaching the eye's altitude. */    sunTransmittance: () => {      'use gpu';      return d.vec3f(sunlight.$);    },    /** The share of sunlight reaching `metres` above the sea, lit after the sun sets at the eye. */    sunTransmittanceAt: (metres: number) => {      'use gpu';      const position = d.vec3f(0, GROUND_RADIUS + metres / 1000, 0);      return sunlightAt(position, settings.$.sunDirection);    },    /** The sun's disc, per unit of sunlight. */    sunRadiance: (dir: d.v3f) => {      'use gpu';      const cosine = std.dot(dir, settings.$.sunDirection);      const disc = std.smoothstep(SUN_COS_EDGE - 3e-6, SUN_COS_EDGE + 3e-6, cosine);      return eyeTransmittance(dir) * (disc / SUN_SOLID_ANGLE);    },    /** Extinction and in-scattering along a near-level ray of `distance` metres. */    aerialPerspective: (color: d.v3f, dir: d.v3f, distance: number) => {      'use gpu';      return veiled(color, d.f32(1), dir, distance);    },    veiled,    /** Moves the sun and the eye, in metres above the sea. Record `bake` after it. */    place(sunDirection: d.v3f, eyeHeight: number) {      settings.patch({        sunDirection,        eyeRadius: GROUND_RADIUS + Math.max(eyeHeight, 0.5) / 1000,      });    },    /** Records the four bakes; record them again after `place`. */    bake(pass: TgpuComputePass) {      transmittancePipeline        .with(pass)        .dispatchWorkgroups(TRANSMITTANCE_SIZE.x / TILE, TRANSMITTANCE_SIZE.y / TILE);      multiScatterPipeline        .with(pass)        .dispatchWorkgroups(MULTI_SCATTER_SIZE / TILE, MULTI_SCATTER_SIZE / TILE);      skyViewPipeline        .with(pass)        .dispatchWorkgroups(Math.ceil(SKY_VIEW_SIZE.x / TILE), Math.ceil(SKY_VIEW_SIZE.y / TILE));      lightingPipeline.with(pass).dispatchWorkgroups(1);    },    destroy: own.destroy,  };}

clouds.ts

import {  cameras,  history,  memo,  owner,  targets,  tone,  type CameraBase,  type SizeSource,} from '@typegpu/kit';import { perlin3d, randf } from '@typegpu/noise';import { d, std, tgpu, type TgpuComputePass, type TgpuRoot } from 'typegpu';import type { Atmosphere } from './atmosphere.ts';const NOISE_SIZE = 96;/** The volumes' texels, and the box they fill: metres across, base height and depth. */const VOLUME = d.vec3u(192, 48, 192);const EXTENT = 40_000;const BASE = 800;const HEIGHT = 1600;/** Extinction per metre at density 1. */const SIGMA = 0.005;const SHADOW_SIZE = VOLUME.x * 3;const VIEW_STEPS = 48;const BAKE_STEPS = 40;/** Below this transmittance the march is behind the surface the eye sees. */const SURFACE = 0.3;/** The share of each new march the view keeps, so the steps' offsets average out over frames. */const BLEND = 0.1;const CHECKER = 2;const GOLDEN = 0.618034;/** View heights over which the march fades into the atmosphere below the deck. */const HORIZON = { from: 0.008, to: 0.02 };const BILLOWS = { metres: 1600, octave: 2.7, erosion: 0.45 };/** A low sun's horizon glow on the deck's base; `cone` is a cosine, about 14° round the sun. */const UNDERLIGHT = {  gain: 2.2,  reach: 4000,  cone: 0.97,  fadeFrom: 0.1,  fadeTo: 0.3,  /** Sampled this far round the horizon from the sun, clear of its aureole. */  cosine: Math.cos(0.6),  sine: Math.sin(0.6),};/** The glow averages the bright band under the deck's edge at these angles from the sun. */const BAND = [0.6, -0.6, 1.4, -1.4].map((angle) => ({  cosine: Math.cos(angle),  sine: Math.sin(angle),}));/** The horizon bank and the overcast beyond the box; `bank` and `veils` are view heights. */const FAR = {  bank: 0.0055,  veils: 0.15,  veil: 0.97,  layer: BASE + HEIGHT / 2,  levels: [0.1, 0.35, 0.6],  hazed: 150_000,  resolved: [20_000, 50_000] as const,  pile: 0.02,  forward: 0.85,  lit: 0.4,  under: 0.25,  beneath: [0.02, 0.15] as const,  metres: 25_000,};/** Zenith light on the deck per unit of escaping light, and the share of its colour kept. */const SKY_FILL = { gain: 0.8, tint: 0.5 };/** The deck's weather. Lengths are noise periods in metres that divide the box, so it tiles. */const WEATHER = {  large: EXTENT / 3,  small: EXTENT / 6,  strength: EXTENT,  edge: 0.3,  spread: 0.5,  cell: 0.7,  overcast: [0.6, 1] as const,  floor: 0.6,};/** The clouds' shape. `metres` and `detail` are noise periods that divide the box. */const SHAPE = {  metres: EXTENT / 8,  warp: 0.15,  range: [0.5, 0.9] as const,  detail: EXTENT / 13,  erosion: 0.4,};/** The cells' vertical profile, in shares of the box's height. */const PROFILE = { base: 0.03, rise: 0.06, rolls: 0.15, lowest: 0.3 };/** Rays climbing to the deck leave the low haze early: they cross this share of its depth. */const DECK_HAZE = 0.01;/** How far a point looks for the billows that shade it, and how dense their shade is. */const SELF_SHADOW = { metres: 180, depth: 3 };/** Hillaire's scattering octaves, a colour channel each; the third holds diffused light. */const ORDERS = [  { channel: 0, attenuation: 1, weight: 1, anisotropy: 0.65 },  { channel: 1, attenuation: 0.5, weight: 0.55, anisotropy: 0.3 },] as const;const DIFFUSE = { channel: 2 };const Shape = d.struct({  coverage: d.f32,  drift: d.vec2f,  /** The frame count, and the share of the new march kept: 1 when the whole view is marched. */  turn: d.f32,  blend: d.f32,});const viewWriteLayout = tgpu.bindGroupLayout({  target: { storageTexture: d.textureStorage2d('rgba16float', 'write-only') },});const cloudsLayout = tgpu.bindGroupLayout({ seen: { texture: d.texture2d(d.f32) } });const resolveLayout = tgpu.bindGroupLayout({  fresh: { texture: d.texture2d(d.f32) },  kept: { texture: d.texture2d(d.f32) },  target: { storageTexture: d.textureStorage2d('rgba16float', 'write-only') },});/** Where a direction meets the screen, as clip x, y and w. */const vanishing = (direction: d.v3f) => {  'use gpu';  const clip = cameras.current.$.viewProjection * d.vec4f(direction, 0);  return d.vec3f(clip.x, clip.y, clip.w);};const vanishingBefore = history.previously(vanishing);const henyeyGreenstein = (cosTheta: number, g: number) => {  'use gpu';  const denom = 1 + g * g - 2 * g * cosTheta;  return (1 - g * g) / (4 * Math.PI * denom * std.sqrt(denom));};const segmentTransmission = (density: number, distance: number) => {  'use gpu';  return std.exp(-std.max(density, 0) * SIGMA * std.max(distance, 0));};/** Diffuse transmission through a slab in the diffusion limit: one over its depth. */const diffuseTransmission = (opticalDepth: number) => {  'use gpu';  return 1 / (1 + 0.75 * (1 - ORDERS[0].anisotropy) * std.max(opticalDepth, 0));};const worley = (point: d.v3f, cells: number) => {  'use gpu';  const p = point * cells;  const base = std.floor(p);  let nearest = d.f32(2);  for (const z of std.range(-1, 2)) {    for (const y of std.range(-1, 2)) {      for (const x of std.range(-1, 2)) {        const cell = base + d.vec3f(d.f32(x), d.f32(y), d.f32(z));        randf.seed3(std.mod(std.mod(cell, d.f32(cells)) + cells, d.f32(cells)));        const feature = cell + d.vec3f(randf.sample(), randf.sample(), randf.sample());        nearest = std.min(nearest, std.distance(p, feature));      }    }  }  return std.saturate(1 - nearest);};const dome = (cone: number) => {  'use gpu';  return 1 - (1 - cone) * (1 - cone);};const volumeUv = (point: d.v3f) => {  'use gpu';  return d.vec3f(point.x / EXTENT + 0.5, (point.y - BASE) / HEIGHT, point.z / EXTENT + 0.5);};/** How far the deck has thinned out towards the box's sides, from 0 to 1. */const rimOf = (uv: d.v2f) => {  'use gpu';  return std.smoothstep(0.38, 0.5, std.max(std.abs(uv.x - 0.5), std.abs(uv.y - 0.5)));};export interface CloudsOptions {  /** Sky the deck covers, from 0 to 1. Default 0.7. */  readonly coverage?: number;}export type Clouds = ReturnType<typeof createClouds>;/** A stratocumulus deck over a 40 km box round the eye, tiled so it drifts without a new bake. */export function createClouds(  root: TgpuRoot,  atmosphere: Atmosphere,  view: CameraBase,  screen: SizeSource,  options: CloudsOptions = {},) {  const own = owner();  const shape = root.createUniform(Shape, {    coverage: options.coverage ?? 0.7,    drift: d.vec2f(),    turn: 0,    blend: 1,  });  const repeat = root.createSampler({    minFilter: 'linear',    magFilter: 'linear',    addressModeU: 'repeat',    addressModeV: 'repeat',    addressModeW: 'repeat',  });  const volumeSampler = root.createSampler({    minFilter: 'linear',    magFilter: 'linear',  });  /** The deck tiles its box across but not up. */  const deckSampler = root.createSampler({    minFilter: 'linear',    magFilter: 'linear',    addressModeU: 'repeat',    addressModeW: 'repeat',  });  const volume = () =>    own(      root        .createTexture({          size: [VOLUME.x, VOLUME.y, VOLUME.z],          dimension: '3d',          format: 'rgba16float',        })        .$usage('sampled', 'storage'),    );  const noise = own(    root      .createTexture({        size: [NOISE_SIZE, NOISE_SIZE, NOISE_SIZE],        dimension: '3d',        format: 'rgba8unorm',      })      .$usage('sampled', 'storage'),  );  const noiseWrite = noise.createView(d.textureStorage3d('rgba8unorm', 'write-only'));  const noiseRead = noise.createView(d.texture3d(d.f32));  const density = volume();  const densityRead = density.createView(d.texture3d(d.f32));  const densityWrite = density.createView(d.textureStorage3d('rgba16float', 'write-only'));  const lighting = volume();  const lightRead = lighting.createView(d.texture3d(d.f32));  const lightWrite = lighting.createView(d.textureStorage3d('rgba16float', 'write-only'));  const shadowMap = own(    root      .createTexture({ size: [SHADOW_SIZE, SHADOW_SIZE], format: 'rgba8unorm' })      .$usage('sampled', 'storage'),  );  const shadowRead = shadowMap.createView(d.texture2d(d.f32));  const shadowWrite = shadowMap.createView(d.textureStorage2d('rgba8unorm', 'write-only'));  const half = own(    targets.create(root, screen, {      fresh: { format: 'rgba16float', sampled: true, storage: true, scale: 0.5 },      even: { format: 'rgba16float', sampled: true, storage: true, scale: 0.5 },      odd: { format: 'rgba16float', sampled: true, storage: true, scale: 0.5 },    }),  );  type Half = typeof half.textures.fresh;  const resolveGroupOf = memo((fresh: Half, kept: Half, target: Half) =>    root.createBindGroup(resolveLayout, {      fresh: fresh.createView(d.texture2d(d.f32)),      kept: kept.createView(d.texture2d(d.f32)),      target: target.createView(d.textureStorage2d('rgba16float', 'write-only')),    }),  );  const writeGroupOf = memo((texture: Half) =>    root.createBindGroup(viewWriteLayout, {      target: texture.createView(d.textureStorage2d('rgba16float', 'write-only')),    }),  );  const readGroupOf = memo((texture: Half) =>    root.createBindGroup(cloudsLayout, { seen: texture.createView(d.texture2d(d.f32)) }),  );  // Billowy Perlin-Worley in red, fine Worley in green, Perlin in blue, coarse Worley in alpha  const gradients = own(perlin3d.staticCache({ root, size: d.vec3u(16) }));  const bakeNoise = root.pipe(gradients.inject()).createComputePipeline({    compute: tgpu.computeFn({      workgroupSize: [4, 4, 4],      in: { gid: d.builtin.globalInvocationId },    })(({ gid }) => {      'use gpu';      if (std.any(std.ge(gid, d.vec3u(NOISE_SIZE)))) return;      const uv = (d.vec3f(gid) + 0.5) / NOISE_SIZE;      let perlin = d.f32(0);      for (const octave of tgpu.unroll(std.range(3))) {        perlin += perlin3d.sample(uv * (16 * 2 ** octave)) * 0.5 ** (octave + 1);      }      const base = std.saturate(perlin * 0.5 + 0.5);      const lumps = worley(uv, d.f32(4));      const billows = std.saturate((base + 1 - lumps) / (2 - lumps));      std.textureStore(noiseWrite.$, gid, d.vec4f(billows, worley(uv, d.f32(8)), base, lumps));    }),  });  const slice = (ground: d.v2f, metres: number, at: number) => {    'use gpu';    return std.textureSampleLevel(      noiseRead.$,      repeat.$,      d.vec3f(ground.x / metres, at, ground.y / metres),      0,    );  };  /** The weather over a ground point: its cells' share of the air, their strength, overcast. */  const weatherAt = (ground: d.v2f) => {    'use gpu';    const large = slice(ground, WEATHER.large, 0.31).w;    const small = slice(ground, WEATHER.small, 0.67).y;    const strength = slice(ground, WEATHER.strength, 0.13).z;    const cells = std.max(large, small * 0.9) * std.mix(0.7, 1.25, strength);    const cover = shape.$.coverage;    const overcast = std.smoothstep(WEATHER.overcast[0], WEATHER.overcast[1], cover);    const scattered =      std.saturate((cells - (1 - std.pow(cover, WEATHER.spread))) / WEATHER.edge) * WEATHER.cell;    const local =      std.mix(        scattered,        std.max(scattered, WEATHER.floor + (1 - WEATHER.floor) * std.smoothstep(0.2, 0.8, small)),        overcast,      ) * std.smoothstep(0, 0.05, cover);    return d.vec3f(local, large, overcast);  };  const bakeDensity = root.createComputePipeline({    compute: tgpu.computeFn({      workgroupSize: [4, 4, 4],      in: { gid: d.builtin.globalInvocationId },    })(({ gid }) => {      'use gpu';      if (std.any(std.ge(gid, VOLUME))) return;      const uv = (d.vec3f(gid) + 0.5) / d.vec3f(VOLUME);      // The volume holds the deck in its own frame, which drifts      const ground = (uv.xz - 0.5) * EXTENT;      const weather = weatherAt(ground);      const local = weather.x;      const large = weather.y;      const overcast = weather.z;      const base = PROFILE.base + PROFILE.rolls * overcast * slice(ground, SHAPE.metres, 0.83).z;      const ceiling = std.mix(PROFILE.lowest, 1, local * std.mix(0.5, 1, large));      const profile =        std.smoothstep(base, base + PROFILE.rise, uv.y) *        (1 - std.smoothstep(ceiling * 0.5, ceiling, uv.y));      const point = d.vec3f(ground.x, uv.y * HEIGHT, ground.y) / SHAPE.metres;      const twist = std.textureSampleLevel(noiseRead.$, repeat.$, point * 0.5 + 0.21, 0);      const warp = (twist.yz - 0.5) * SHAPE.warp;      const billows = std.textureSampleLevel(        noiseRead.$,        repeat.$,        point + d.vec3f(warp.x, 0, warp.y),        0,      ).x;      // An overcast deck keeps its billows as thicker and thinner parts rather than holes      const form =        std.mix(1, 0.5, overcast) * std.smoothstep(SHAPE.range[0], SHAPE.range[1], billows) +        0.5 * overcast;      // Cells keep only the shape's strongest parts where the weather is weak      const value = std.saturate((form * profile - (1 - local)) / std.max(local, 0.15));      const detail = std.textureSampleLevel(        noiseRead.$,        repeat.$,        (point * SHAPE.metres) / SHAPE.detail,        0,      ).y;      const carved = std.saturate(value - detail * SHAPE.erosion * (1 - value));      std.textureStore(densityWrite.$, gid, d.vec4f(carved, 0, 0, 1));    }),  });  const densityAt = (uv: d.v3f) => {    'use gpu';    if (uv.y <= 0 || uv.y >= 1) return d.f32(0);    return std.textureSampleLevel(densityRead.$, deckSampler.$, uv, 0).x;  };  /** Sunlight through the deck towards the sun, and sky light escaping upwards. */  const bakeLighting = root.createComputePipeline({    compute: tgpu.computeFn({      workgroupSize: [4, 4, 4],      in: { gid: d.builtin.globalInvocationId },    })(({ gid }) => {      'use gpu';      if (std.any(std.ge(gid, VOLUME))) return;      const uv = (d.vec3f(gid) + 0.5) / d.vec3f(VOLUME);      const amount = std.textureLoad(densityRead.$, gid, 0).x;      let orders = d.vec3f();      let skyEscape = d.f32(0);      if (amount > 0) {        const sun = atmosphere.sunDirection();        const length = std.min(((1 - uv.y) * HEIGHT) / std.max(sun.y, 0.08), 8000);        let opticalDepth = d.f32(0);        for (const tap of std.range(16)) {          const start = length * (d.f32(tap) / 16) ** 2;          const end = length * ((d.f32(tap) + 1) / 16) ** 2;          const at = uv + (sun * ((start + end) * 0.5)) / d.vec3f(EXTENT, HEIGHT, EXTENT);          opticalDepth += densityAt(at) * (end - start) * SIGMA;        }        const path = (1 - uv.y) * HEIGHT;        let skyDepth = d.f32(0);        for (const tap of std.range(8)) {          const start = path * (d.f32(tap) / 8) ** 2;          const end = path * ((d.f32(tap) + 1) / 8) ** 2;          const at = uv + d.vec3f(0, (start + end) * 0.5, 0) / HEIGHT;          skyDepth += densityAt(at) * (end - start) * SIGMA;        }        skyEscape = diffuseTransmission(skyDepth);        for (const order of tgpu.unroll(ORDERS)) {          orders[order.channel] = order.weight * std.exp(-opticalDepth * order.attenuation);        }        // Scattered isotropically, at 1/4π, it leaves as μ T / π, a diffuser's radiance        orders[DIFFUSE.channel] = 4 * std.max(sun.y, 0.005) * diffuseTransmission(opticalDepth);      }      std.textureStore(lightWrite.$, gid, d.vec4f(orders * amount, skyEscape * amount));    }),  });  const bakeShadow = root.createComputePipeline({    compute: tgpu.computeFn({ workgroupSize: [8, 8], in: { gid: d.builtin.globalInvocationId } })(      ({ gid }) => {        'use gpu';        if (gid.x >= SHADOW_SIZE || gid.y >= SHADOW_SIZE) return;        const xz = (d.vec2f(gid.xy) + 0.5) / SHADOW_SIZE;        const sun = atmosphere.sunDirection();        const step = HEIGHT / std.max(sun.y, 0.08) / 32;        let transmittance = d.f32(1);        for (const tap of std.range(32)) {          const at =            d.vec3f(xz.x, 0, xz.y) +            (sun * ((d.f32(tap) + 0.5) * step)) / d.vec3f(EXTENT, HEIGHT, EXTENT);          transmittance *= segmentTransmission(densityAt(at), step);          if (transmittance < 0.005) break;        }        std.textureStore(shadowWrite.$, gid.xy, d.vec4f(transmittance, 0, 0, 1));      },    ),  });  /** Erosion of the deck's edges, fading where it would alias; `detailed` adds a finer octave. */  const billowsOf = (detailed: boolean) => (point: d.v3f, footprint: number) => {    'use gpu';    const at = point / BILLOWS.metres;    const cells = std.textureSampleLevel(noiseRead.$, repeat.$, at, 0);    const fine = 1 - std.smoothstep(1, 2.5, (footprint * NOISE_SIZE) / BILLOWS.metres);    const smaller =      1 - std.smoothstep(1, 2.5, (footprint * NOISE_SIZE * BILLOWS.octave) / BILLOWS.metres);    let small = d.f32(0.5);    if (detailed) {      const finer = std.textureSampleLevel(noiseRead.$, repeat.$, at * BILLOWS.octave + 0.37, 0).y;      small = std.mix(0.5, dome(finer), smaller);    }    const billow = 0.5 * dome(cells.y) + 0.25 * (1 - cells.z) + 0.25 * small;    return billow * BILLOWS.erosion * fine;  };  const erodedOf = (detailed: boolean) => {    const billowsAt = billowsOf(detailed);    return (amount: number, height: number, point: d.v3f, footprint: number) => {      'use gpu';      const billow = billowsAt(point, footprint);      // Billows bite deepest into the tops; the bases stay flat where the air condenses      const bite = billow * std.mix(0.6, 1.3, std.smoothstep(0.05, 0.45, height));      const eaten = std.saturate((amount * (1 - 0.5 * bite) - bite) / std.max(1 - bite, 0.05));      // Cumulus ends at a rim rather than fading out, so faint density rises quickly to a body      return std.max(eaten, 0.3 * std.smoothstep(0.02, 0.1, eaten));    };  };  const eroded = erodedOf(true);  const shading = erodedOf(false);  const erodedAt = (point: d.v3f, footprint: number) => {    'use gpu';    const uv = volumeUv(point);    return shading(densityAt(uv), uv.y, point, footprint);  };  /** The far bank and veils behind the deck: their light, and the share of sky showing through. */  const farAt = (ray: d.v3f, sunlight: d.v3f, skylight: d.v3f, glow: d.v3f) => {    'use gpu';    const height = ray.y / std.length(ray);    if (height > FAR.veils || height < -0.01) return d.vec4f(0, 0, 0, 1);    const around = std.atan2(ray.x, -ray.z) / (2 * Math.PI);    const outline = std.textureSampleLevel(      noiseRead.$,      repeat.$,      d.vec3f(around * 3, 0.5, 0.21),      0,    );    const heavy = std.smoothstep(0, 0.6, shape.$.coverage);    const crest = FAR.bank * (0.35 + outline.z * 0.9 + outline.w * 0.4) * heavy;    const bank = 1 - std.smoothstep(crest - 0.0008, crest + 0.0008, height);    const sun = atmosphere.sunDirection();    const lowest = std.max(height, FAR.layer / FAR.metres / 8);    const distance = FAR.layer / lowest;    const heading = ray.xz / std.length(ray.xz);    // Only an overcast goes on past the box: scattered cells that far are smaller than a pixel    const overcast = std.smoothstep(WEATHER.overcast[0], WEATHER.overcast[1], shape.$.coverage);    let cells = d.f32(0);    let thickness = d.f32(0);    if (overcast > 0) {      // The ray crosses the deck at several heights, and only strong cells reach the upper ones      for (const level of tgpu.unroll(FAR.levels)) {        const reach = (BASE + level * HEIGHT) / lowest;        const across = heading * reach;        const beyond = rimOf(across / EXTENT + 0.5);        const weather = weatherAt(across - shape.$.drift);        const averaged = std.mix(          std.saturate(weather.x / WEATHER.cell),          shape.$.coverage,          std.smoothstep(FAR.resolved[0], FAR.resolved[1], reach),        );        const strong = std.saturate((averaged - level * 0.5) / (1 - level * 0.5));        cells = std.max(cells, strong * beyond);        thickness +=          std.saturate((weather.x - WEATHER.floor) / (1 - WEATHER.floor)) / FAR.levels.length;      }    }    const piled = 1 - std.pow(1 - cells, 1 + FAR.pile / std.max(height, 1e-3));    const veil = FAR.veil * overcast * piled * std.smoothstep(FAR.bank * 0.3, FAR.bank, height);    const through = (1 - bank) * (1 - veil);    // A low sun lights an overcast's base from beyond the deck's edge; a high one diffuses through    const forward = henyeyGreenstein(std.dot(ray, sun) / std.length(ray), FAR.forward);    const high = std.smoothstep(FAR.beneath[0], FAR.beneath[1], sun.y);    const base = std.mix(1, FAR.under * std.mix(1.4, 0.6, thickness), overcast);    const ambient = sunlight * (FAR.lit * std.max(sun.y, 0)) + skylight * 0.4;    const diffuse = ambient * base;    const lit = (sunlight * forward + glow * 0.15) * (1 - overcast * high) + diffuse;    // Under an overcast the haze takes only the light the deck lets through    const hazed = std.mix(DECK_HAZE, 1, std.smoothstep(EXTENT / 2, FAR.hazed, distance));    const dimmed = std.mix(1, FAR.under, overcast);    const hazy = std.mix(1, dimmed, high);    const veils = atmosphere.veiled(lit * veil, veil * hazy, ray, distance * hazed);    // Under a high sun the bank is the same deck as the layer, further out    const shade = std.mix(skylight * 0.35, ambient * dimmed, high * overcast);    const banked = atmosphere.veiled(shade * bank, bank * hazy, ray, FAR.metres);    return d.vec4f(veils * (1 - bank) + banked, through);  };  /**   * Marches a ray from `eye` through the deck: the light per unit of sunlight and the transmittance   * left. Billows finer than `pixel` fade; only a `billowed` march erodes and self-shades.   */  const marcher =    (steps: number, billowed: boolean) =>    (eye: d.v3f, ray: d.v3f, pixel: number, jitter: number) => {      'use gpu';      const up = std.max(ray.y, 0.001);      const near = std.max((BASE - eye.y) / up, 0);      const side =        (d.vec2f(std.sign(ray.xz)) * (EXTENT / 2) - eye.xz) /        std.max(std.abs(ray.xz), d.vec2f(1e-4));      const far = std.min((BASE + HEIGHT - eye.y) / up, std.min(std.abs(side.x), std.abs(side.y)));      const range = std.max(far - near, 0);      const sun = atmosphere.sunDirection();      const mu = std.clamp(std.dot(ray, sun), -1, 1);      let phases = d.vec3f();      for (const order of tgpu.unroll(ORDERS)) {        phases[order.channel] = henyeyGreenstein(mu, order.anisotropy);      }      phases[DIFFUSE.channel] = 1 / (4 * Math.PI);      // The deck's own altitude keeps the sun a little after it sets at the eye      const sunlight = atmosphere.sunTransmittanceAt(BASE + HEIGHT * 0.25);      const zenith = atmosphere.skyRadiance(d.vec3f(0, 1, 0));      const skylight =        std.mix(d.vec3f(tone.luminance(zenith)), zenith, SKY_FILL.tint) * SKY_FILL.gain;      const aside = d.vec2f(        sun.x * UNDERLIGHT.cosine - sun.z * UNDERLIGHT.sine,        sun.x * UNDERLIGHT.sine + sun.z * UNDERLIGHT.cosine,      );      let band = d.vec3f();      for (const turn of tgpu.unroll(BAND)) {        const along = d.vec2f(          sun.x * turn.cosine - sun.z * turn.sine,          sun.x * turn.sine + sun.z * turn.cosine,        );        band += atmosphere.skyRadiance(std.normalize(d.vec3f(along.x, 0.03, along.y)));      }      const glow =        band *        ((UNDERLIGHT.gain / BAND.length) *          (1 - std.smoothstep(UNDERLIGHT.fadeFrom, UNDERLIGHT.fadeTo, sun.y)));      const toward = std.normalize(d.vec2f(sun.x, sun.z) + d.vec2f(1e-5, 0));      const reach = (EXTENT * 0.44) / std.max(std.abs(toward.x), std.abs(toward.y));      const below = std.normalize(d.vec3f(aside.x, -1, aside.y)) * SELF_SHADOW.metres;      const overhead = std.normalize(d.vec3f(sun.x, 1, sun.z)) * SELF_SHADOW.metres;      const drift = d.vec3f(shape.$.drift.x, 0, shape.$.drift.y);      const ember = std.smoothstep(UNDERLIGHT.cone, 0.999, mu);      let throughput = d.f32(1);      let radiance = d.vec3f();      let hazeDepth = d.f32(0);      if (far > near && ray.y > HORIZON.from) {        for (const tap of std.range(steps)) {          // Steps start short where the ray enters the deck, the part the eye sees          const start = (d.f32(tap) / steps) ** 2;          const end = ((d.f32(tap) + 1) / steps) ** 2;          const step = range * (end - start);          const distance = near + range * std.mix(start, end, jitter);          const world = eye + ray * distance;          const point = world - drift;          const uv = volumeUv(point);          const footprint = distance * pixel;          const amount = std.textureSampleLevel(densityRead.$, deckSampler.$, uv, 0).x;          if (amount < 1e-3) continue;          const light = std.textureSampleLevel(lightRead.$, deckSampler.$, uv, 0);          let local = amount * (1 - rimOf(volumeUv(world).xz));          if (billowed) {            // The finer billows show only on the surface the eye sees, not behind it            if (throughput > SURFACE) {              local = eroded(local, uv.y, point, footprint);            } else {              local = shading(local, uv.y, point, footprint);            }          }          const segment = segmentTransmission(local, step);          const inverse = 1 / std.max(amount, 1e-6);          // Folds above a point dim the light from above, folds below it the horizon's glow          const above = billowed ? erodedAt(point + overhead, footprint) : d.f32(0);          const shade = segmentTransmission(above, d.f32(SELF_SHADOW.metres) * SELF_SHADOW.depth);          let underglow = d.vec3f();          if (ember > 0 && uv.y < 0.35) {            const under = billowed              ? erodedAt(point + below, footprint)              : densityAt(volumeUv(point + below));            const inward = reach - std.dot(world.xz, toward);            const gap = std.exp(-std.max(inward, 0) / UNDERLIGHT.reach);            underglow =              glow *              (segmentTransmission(under, d.f32(SELF_SHADOW.metres) * SELF_SHADOW.depth) *                (1 - std.smoothstep(0.05, 0.35, uv.y)) *                gap *                ember);          }          const lit =            (sunlight * (std.dot(light.xyz, phases) * inverse) +              skylight * std.saturate(light.w * inverse)) *              shade +            underglow;          radiance += lit * (throughput * (1 - segment));          hazeDepth += distance * throughput * (1 - segment);          throughput *= segment;          if (throughput < 0.01) break;        }      }      const fade = std.smoothstep(HORIZON.from, HORIZON.to, ray.y);      const cover = (1 - throughput) * fade;      const deck = atmosphere.veiled(        radiance * fade,        cover,        ray,        (hazeDepth / std.max(1 - throughput, 1e-4)) * DECK_HAZE,      );      const beyond = farAt(ray, sunlight, skylight, glow);      const through = 1 - cover;      return d.vec4f(deck + beyond.xyz * through, through * beyond.w);    };  const viewMarch = marcher(VIEW_STEPS, true);  const bakeMarch = marcher(BAKE_STEPS, false);  const viewKernel = root.pipe(view.inject()).createComputePipeline({    compute: tgpu.computeFn({ workgroupSize: [8, 8], in: { gid: d.builtin.globalInvocationId } })(      ({ gid }) => {        'use gpu';        const size = std.textureDimensions(viewWriteLayout.$.target);        let at = gid.xy;        if (shape.$.blend < 1) {          at = d.vec2u(gid.x * CHECKER + ((gid.y + d.u32(shape.$.turn)) % CHECKER), gid.y);        }        if (at.x >= size.x || at.y >= size.y) return;        const camera = cameras.current.$;        const lens = cameras.lensOf(camera);        const pixel = std.max(2 / (lens.x * d.f32(size.x)), 2 / (lens.y * d.f32(size.y)));        const uv = (d.vec2f(at) + 0.5) / d.vec2f(size);        const ray = cameras.rayOf(camera, d.vec2f(uv.x * 2 - 1, 1 - uv.y * 2)).direction;        // Each pixel's step offset moves on each frame, so the slices average out in the history        randf.seed2(d.vec2f(at));        const jitter = std.fract(randf.sample() + shape.$.turn * GOLDEN);        std.textureStore(          viewWriteLayout.$.target,          at,          viewMarch(camera.position, ray, pixel, jitter),        );      },    ),  });  const resolve = root.pipe(view.inject()).createComputePipeline({    compute: tgpu.computeFn({ workgroupSize: [8, 8], in: { gid: d.builtin.globalInvocationId } })(      ({ gid }) => {        'use gpu';        const size = std.textureDimensions(resolveLayout.$.target);        if (gid.x >= size.x || gid.y >= size.y) return;        const fresh = std.textureLoad(resolveLayout.$.fresh, gid.xy, 0);        // History outside the new march's neighbourhood is stale: the sun, cover or view moved        let low = d.vec4f(fresh);        let high = d.vec4f(fresh);        for (const y of std.range(-1, 2)) {          for (const x of std.range(-1, 2)) {            const at = std.clamp(d.vec2i(gid.xy) + d.vec2i(x, y), d.vec2i(0), d.vec2i(size) - 1);            const near = std.textureLoad(resolveLayout.$.fresh, at, 0);            low = std.min(low, near);            high = std.max(high, near);          }        }        const uv = (d.vec2f(gid.xy) + 0.5) / d.vec2f(size);        const ray = cameras.rayOf(cameras.current.$, d.vec2f(uv.x * 2 - 1, 1 - uv.y * 2)).direction;        // The deck is far enough that last frame saw each direction where its camera points it        const before = vanishingBefore(ray);        let result = d.vec4f(fresh);        if (before.z > 0) {          const ndc = before.xy / before.z;          const then = d.vec2f(ndc.x * 0.5 + 0.5, 0.5 - ndc.y * 0.5);          if (std.all(std.ge(then, d.vec2f(0))) && std.all(std.le(then, d.vec2f(1)))) {            const kept = std.textureSampleLevel(resolveLayout.$.kept, volumeSampler.$, then, 0);            // A pixel the march skipped this frame holds last frame's march, already in the history            const marched =              shape.$.blend >= 1 || gid.x % CHECKER === (gid.y + d.u32(shape.$.turn)) % CHECKER;            const share = std.select(0, shape.$.blend, marched);            result = std.mix(std.clamp(kept, low, high), fresh, share);          }        }        std.textureStore(resolveLayout.$.target, gid.xy, result);      },    ),  });  let turn = 0;  let resolved = 0;  let started = false;  let noiseBaked = false;  return {    pipelines: [bakeNoise, bakeDensity, bakeLighting, bakeShadow, viewKernel, resolve],    /** The light and transmittance along any ray, for an environment's bake; see `marcher`. */    march: (eye: d.v3f, ray: d.v3f, pixel: number) => {      'use gpu';      return bakeMarch(eye, ray, pixel, 0.5);    },    /** The eye's view at screen `uv`, top left first, from the half-size march. Bind `group`. */    seen: (uv: d.v2f) => {      'use gpu';      return std.textureSampleLevel(cloudsLayout.$.seen, volumeSampler.$, uv, 0);    },    /** The share of sunlight the deck lets through to a point on the sea. */    shadow: (xz: d.v2f) => {      'use gpu';      const sun = atmosphere.sunDirection();      // Where sunlight enters the deck's base; under the deck's thinned edge it reaches the sea      const entry = (xz + sun.xz * (BASE / std.max(sun.y, 0.01))) / EXTENT + 0.5;      const drifted = entry - shape.$.drift / EXTENT;      const shade = std.textureSampleLevel(shadowRead.$, repeat.$, drifted, 0).x;      return std.mix(shade, 1, rimOf(std.clamp(entry, d.vec2f(0), d.vec2f(1))));    },    /** The read group of the eye's view, for drawings that call `seen`. */    get group() {      return readGroupOf(resolved === 0 ? half.textures.even : half.textures.odd);    },    /** Moves the deck to `drift` metres from where it started. */    move(drift: readonly [number, number]) {      shape.patch({ drift: d.vec2f(drift[0], drift[1]) });    },    /** Sets the share of the sky the deck covers, from 0 to 1. Record `bake` after it. */    cover(coverage: number) {      shape.patch({ coverage });    },    /** Records the volume's bake: after creation, and after the sun or the cover changes. */    bake(pass: TgpuComputePass) {      if (!noiseBaked) {        const n = NOISE_SIZE / 4;        bakeNoise.with(pass).dispatchWorkgroups(n, n, n);        noiseBaked = true;      }      const [x, y, z] = [VOLUME.x / 4, VOLUME.y / 4, VOLUME.z / 4];      bakeDensity.with(pass).dispatchWorkgroups(x, y, z);      bakeLighting.with(pass).dispatchWorkgroups(x, y, z);      bakeShadow.with(pass).dispatchWorkgroups(SHADOW_SIZE / 8, SHADOW_SIZE / 8);    },    /** Records the march of half the eye's view, at half size, and its blend into the history. */    record(pass: TgpuComputePass) {      const resized = half.fit();      const whole = !started || resized;      shape.patch({ turn, blend: whole ? 1 : BLEND });      turn = (turn + 1) % 4096;      started = true;      const { fresh, even, odd } = half.textures;      const [width, height] = fresh.props.size;      const groups = [Math.ceil(width / 8), Math.ceil(height / 8)] as const;      const across = whole ? width : Math.ceil(width / CHECKER);      viewKernel        .with(pass)        .with(writeGroupOf(fresh))        .dispatchWorkgroups(Math.ceil(across / 8), groups[1]);      resolved = 1 - resolved;      const [kept, target] = resolved === 0 ? [odd, even] : [even, odd];      resolve        .with(pass)        .with(resolveGroupOf(fresh, kept, target))        .dispatchWorkgroups(...groups);    },    destroy: own.destroy,  };}

foam.ts

import { mipmaps, owner } from '@typegpu/kit';import { perlin2d } from '@typegpu/noise';import { d, std, tgpu, type TgpuComputePass, type TgpuRoot } from 'typegpu';const TILE = 16;/** The foam tile; a level `t` keeps about 1 − t of its red field. Lengths are in metres. */const PATTERN = {  texels: 1024,  mips: 11,  cells: 8,  size: 9,  streaks: 2.5,  stretch: 16,  drawn: 0.5,  detail: 4.3,  holes: 0.35,  blur: d.vec2f(0.03, 0.25),};/** Foam within `FILM` of its level is a film that hides `FILM_OPACITY` of the water. */const FILM = 0.25;const FILM_OPACITY = 0.45;/** The share of the water below which foam fades out instead of breaking into specks. */const SPARSE = 0.2;/** Octaves at whole multiples of the lattice, which the cache wraps, keep the tile seamless. */const fbm = (at: d.v2f, octaves: number, roughness: number) => {  'use gpu';  let p = d.vec2f(at);  let sum = d.f32(0);  let weight = d.f32(0.5);  for (let octave = 0; octave < octaves; octave++) {    sum += perlin2d.sample(p) * weight;    p = p * 2 + 0.37;    weight *= roughness;  }  return sum;};/** The pixel's share a `field` covers above `level`: crisp up close, averaged when `blurred`. */const above = (field: number, level: number, blurred: number) => {  'use gpu';  const margin = field - level;  const pixel = std.max(std.fwidth(margin), 1e-4);  const crisp = std.smoothstep(-pixel, pixel, margin);  const share = std.saturate(1 - level);  return std.mix(crisp, share, blurred) * std.smoothstep(0, SPARSE, share);};export type Foam = ReturnType<typeof createFoam>;/** * Whitecap lace: a tile of torn foam that `layer` cuts to the share of the water the breaking * waves call for, so a fresh cap closes over and an older one erodes into lace along the wind. */export function createFoam(root: TgpuRoot, options: { wind: number }) {  const own = owner();  const wind = d.vec2f(Math.cos(options.wind), Math.sin(options.wind));  const tile = own(    root      .createTexture({        size: [PATTERN.texels, PATTERN.texels],        format: 'rgba16float',        mipLevelCount: PATTERN.mips,      })      .$usage('sampled', 'storage'),  );  const tileWrite = tile.createView(d.textureStorage2d('rgba16float', 'write-only'), {    baseMipLevel: 0,    mipLevelCount: 1,  });  const tileRead = tile.createView(d.texture2d(d.f32));  const levels = mipmaps.create(root, tile, { combine: mipmaps.mean });  const tiling = root.createSampler({    magFilter: 'linear',    minFilter: 'linear',    mipmapFilter: 'linear',    addressModeU: 'repeat',    addressModeV: 'repeat',    maxAnisotropy: 4,  });  const lattice = own(perlin2d.staticCache({ root, size: d.vec2u(PATTERN.cells, PATTERN.cells) }));  const bake = root.pipe(lattice.inject()).createComputePipeline({    compute: tgpu.computeFn({      workgroupSize: [TILE, TILE],      in: { gid: d.builtin.globalInvocationId },    })(({ gid }) => {      'use gpu';      if (gid.x >= PATTERN.texels || gid.y >= PATTERN.texels) return;      const at = ((d.vec2f(gid.xy) + 0.5) / PATTERN.texels) * PATTERN.cells;      // Warped by its own low octaves, the clumps and holes come out irregular and torn      const warp = d.vec2f(fbm(at + 5.2, 2, 0.5), fbm(at + 11.7, 2, 0.5));      const field = 0.5 + 0.5 * std.tanh(fbm(at + warp * 1.6, 7, 0.62) * 3.4);      const bubbles = 0.5 + 0.5 * std.tanh(fbm(at * 16 + 7.3, 3, 0.6) * 5);      std.textureStore(tileWrite.$, gid.xy, d.vec4f(field, bubbles, 0, 1));    }),  });  return {    pipelines: [bake, levels],    /**     * Foam at an undisplaced point: the share of the pixel it hides, and the bubbles' brightness.     * `age` runs from 0 for a fresh cap to 1 for spent foam; `footprint` is the pixel in metres.     */    layer: (xz: d.v2f, share: number, age: number, footprint: number) => {      'use gpu';      const caps = std.textureSample(tileRead.$, tiling.$, xz / PATTERN.size);      const along = d.vec2f(        std.dot(xz, d.vec2f(-wind.y, wind.x)),        std.dot(xz, wind) / PATTERN.stretch,      );      const streaks = std.textureSample(tileRead.$, tiling.$, along / PATTERN.streaks + 0.61);      const blurred = std.smoothstep(PATTERN.blur.x, PATTERN.blur.y, footprint);      const fine = std.textureSample(        tileRead.$,        tiling.$,        (xz * PATTERN.detail) / PATTERN.size + 0.37,      ).x;      const field = std.mix(caps.x, streaks.x, PATTERN.drawn * age);      const dense = std.mix(std.smoothstep(0, FILM, field - (1 - share)), 1, blurred);      // Inside, the foam thins over the finer octave into lace and bright clumps      const body = std.mix(FILM_OPACITY, 1, dense) * (1 - PATTERN.holes * (1 - fine));      return d.vec2f(above(field, 1 - share, blurred) * body, fine);    },    /** Records the tile and its mips, once, before the first frame. */    bake(pass: TgpuComputePass) {      bake.with(pass).dispatchWorkgroups(PATTERN.texels / TILE, PATTERN.texels / TILE);      levels.record(pass);    },    destroy: own.destroy,  };}

spectrum.ts

import { d, std } from 'typegpu';export const N = 512;export const LOG2N = 9;export const CASCADES = 4;/** The wind waves' cascade, whose crests break and pass light. */export const WIND = 1;/** The ripples' cascade only shades: its waves are shorter than the mesh resolves. */export const RIPPLES = 3;/** Four packed complex signals carry eight real fields per cascade. */export const SIGNALS = 4;export const MIP_LEVELS = LOG2N + 1;/** Rotated, differently sized domains keep the shorter cascades from showing repeated tiles. */const PATCH_SIZES = [2048, 613, 101, 12];const PATCH_ROTATIONS = [13, 47, 101, 71];/** Wavelength boundaries in metres between the swell, wind wave, chop and ripple cascades. */const BAND_EDGES = [640, 35, 5, 0.9, 0.06];const BAND_TAPER = 0.15;const CASCADE_GAINS = [  { height: 1.2, chop: 1 },  { height: 1.15, chop: 1 },  { height: 1, chop: 1 },  { height: 1, chop: 0.3 },];const GRAVITY = 9.81;/** Frequencies snap to multiples of 2π over this many seconds, so the sea repeats exactly. */export const LOOP_PERIOD = 200;const SWELL_PERIOD_RATIO = 1.7;const SWELL_PEAKEDNESS = 7;const SWELL_NARROWING = 1.6;const PEAK_BROADENING = 0.45;const PEAK_BROADENING_WIDTH = 0.3;const ISOTROPIC_FLOOR = 0.35;export const SpectrumParams = d.struct({  windDirection: d.f32,  peakFrequency: d.f32,  energyScale: d.f32,  peakedness: d.f32,  directionality: d.f32,  swellDirection: d.f32,  swellFrequency: d.f32,  swellScale: d.f32,  waveHeight: d.f32,  seed: d.f32,});/** `rotation` holds the cosine and sine of the angle that takes patch space to world space. */export const Cascade = d.struct({  patchSize: d.f32,  kLow: d.f32,  kHigh: d.f32,  maxWavelength: d.f32,  rotation: d.vec2f,  heightGain: d.f32,  chop: d.f32,});/** * A cascade's foam rule, after GodotOceanWaves: foam grows at `grow` per unit its Jacobian falls * below `whitecap`, and decays by `decay`, each second. Padded to a uniform array's 16-byte stride. */export const Foaming = d.struct({ whitecap: d.f32, grow: d.f32, decay: d.size(8, d.f32) });export const cascadeValues: d.Infer<typeof Cascade>[] = PATCH_SIZES.map((patchSize, c) => {  const angle = (PATCH_ROTATIONS[c]! * Math.PI) / 180;  return {    patchSize,    kLow: (2 * Math.PI) / BAND_EDGES[c]!,    kHigh: (2 * Math.PI) / BAND_EDGES[c + 1]!,    maxWavelength: BAND_EDGES[c]!,    rotation: d.vec2f(Math.cos(angle), Math.sin(angle)),    heightGain: CASCADE_GAINS[c]!.height,    chop: CASCADE_GAINS[c]!.chop,  };});/** Each cascade's waves lie within this many cells of its lattice's centre; FFTs skip the rest. */export const REACH = cascadeValues.map(({ kHigh, patchSize }) =>  Math.min(N / 2 - 1, Math.floor(((1 + BAND_TAPER) * kHigh * patchSize) / (2 * Math.PI)) + 1),);export const patchToWorld = (v: d.v2f, rotation: d.v2f) => {  'use gpu';  return d.vec2f(rotation.x * v.x - rotation.y * v.y, rotation.y * v.x + rotation.x * v.y);};export const worldToPatch = (v: d.v2f, rotation: d.v2f) => {  'use gpu';  return d.vec2f(rotation.x * v.x + rotation.y * v.y, rotation.x * v.y - rotation.y * v.x);};/** Rotates a symmetric tensor packed as (xx, yy, xy) from patch space to world space. */export const tensorToWorld = (t: d.v3f, rotation: d.v2f) => {  'use gpu';  const c = rotation.x;  const s = rotation.y;  return d.vec3f(    c * c * t.x - 2 * c * s * t.z + s * s * t.y,    s * s * t.x + 2 * c * s * t.z + c * c * t.y,    c * s * t.x + (c * c - s * s) * t.z - c * s * t.y,  );};/** Determinant of the horizontal deformation, packed as (xx, zz, xz). */export const surfaceJacobian = (deformation: d.v3f) => {  'use gpu';  return (1 + deformation.x) * (1 + deformation.y) - deformation.z * deformation.z;};/** The height gradient through the horizontal deformation, bounded where cells fold. */export const surfaceSlope = (gradient: d.v2f, deformation: d.v3f) => {  'use gpu';  const determinant = std.max(surfaceJacobian(deformation), 0.05);  const x = ((1 + deformation.y) * gradient.x - deformation.z * gradient.y) / determinant;  const z = ((1 + deformation.x) * gradient.y - deformation.z * gradient.x) / determinant;  const scale = std.max(std.sqrt(x * x + z * z) / 8, 1);  return d.vec2f(x, z) / scale;};/** First-order slope covariance after horizontal deformation: M⁻¹ C M⁻ᵀ. */export const surfaceCovariance = (cov: d.v3f, deformation: d.v3f) => {  'use gpu';  const determinant = std.max(surfaceJacobian(deformation), 0.05);  const a = (1 + deformation.y) / determinant;  const b = -deformation.z / determinant;  const c = (1 + deformation.x) / determinant;  return d.vec3f(    a * a * cov.x + 2 * a * b * cov.z + b * b * cov.y,    b * b * cov.x + 2 * b * c * cov.z + c * c * cov.y,    a * b * cov.x + (a * c + b * b) * cov.z + b * c * cov.y,  );};/** Deep water dispersion: omega and d omega / dk. */export const dispersion = (k: number) => {  'use gpu';  const omega = std.sqrt(GRAVITY * k);  return d.vec2f(omega, (0.5 * GRAVITY) / omega);};function jonswapEnergyScale(waveHeight: number, peakPeriod: number, gamma: number) {  const omegaP = (2 * Math.PI) / peakPeriod;  const steps = 4000;  const from = 0.3 * omegaP;  const dw = (12 * omegaP - from) / steps;  let m0 = 0;  for (let i = 0; i < steps; i++) m0 += jonswap(from + (i + 0.5) * dw, omegaP, gamma) * dw;  return (waveHeight * waveHeight) / 16 / m0;}/** A young sea, shorter and steeper than a developed one, spread wide so its crests cross. */const SEA = {  height: [0.25, 0.012],  period: [2, 0.27],  peakedness: 3.3,  directionality: 0.6,  swell: 0.25,  angle: 0.35,};/** The sea's spectrum under a wind of `speed` metres per second, blowing towards `direction`. */export function spectrumAt(speed: number, direction: number): d.Infer<typeof SpectrumParams> {  const waveHeight = SEA.height[0] + SEA.height[1] * speed * speed * Math.sqrt(speed / 18);  const peakPeriod = SEA.period[0] + SEA.period[1] * speed;  const swellPeriod = peakPeriod * SWELL_PERIOD_RATIO;  return {    windDirection: direction,    peakFrequency: (2 * Math.PI) / peakPeriod,    energyScale: jonswapEnergyScale(      Math.sqrt(1 - SEA.swell * SEA.swell) * waveHeight,      peakPeriod,      SEA.peakedness,    ),    peakedness: SEA.peakedness,    directionality: SEA.directionality,    swellDirection: direction + SEA.angle,    swellFrequency: (2 * Math.PI) / swellPeriod,    swellScale: jonswapEnergyScale(SEA.swell * waveHeight, swellPeriod, SWELL_PEAKEDNESS),    waveHeight,    seed: 3.7,  };}/** Only the wind waves foam: their whitecap rises from `calm` to `gale` with the wind. */const FOAMING = { calm: 0.45, gale: 0.9, grow: 67.5, decay: 1.15 };/** Each cascade's foam rule under a wind of `speed` metres per second. */export function foamingAt(speed: number): d.Infer<typeof Foaming>[] {  const gale = Math.min(Math.max((speed - 6) / 16, 0), 1);  const none = { whitecap: 0, grow: 0, decay: 0 };  const whitecap = FOAMING.calm + (FOAMING.gale - FOAMING.calm) * gale;  const foaming = { whitecap, grow: FOAMING.grow, decay: FOAMING.decay };  return Array.from({ length: CASCADES }, (_, c) => (c === WIND ? foaming : none));}const jonswap = (omega: number, omegaP: number, gamma: number) => {  'use gpu';  const sigma = std.select(0.09, 0.07, omega <= omegaP);  const dw = omega - omegaP;  const r = std.exp(-(dw * dw) / (2 * sigma * sigma * omegaP * omegaP));  const ratio = omegaP / omega;  const ratio4 = ratio * ratio * ratio * ratio;  const omega5 = omega * omega * omega * omega * omega;  return ((GRAVITY * GRAVITY) / omega5) * std.exp(-1.25 * ratio4) * std.pow(gamma, r);};const wrapAngle = (theta: number) => {  'use gpu';  return std.atan2(std.sin(theta), std.cos(theta));};/** Donelan-Banner spreading, broadened round the peak so crests stay short. */const donelanBanner = (r: number, theta: number, narrowing: number) => {  'use gpu';  let beta = 2.61 * std.pow(r, 1.3);  if (r >= 1.6) {    beta = std.pow(10, -0.4 + 0.8393 * std.exp(-0.567 * std.log(r * r)));  } else if (r >= 0.95) {    beta = 2.28 * std.pow(r, -1.3);  }  const lnr = std.log(r);  const width = PEAK_BROADENING_WIDTH * PEAK_BROADENING_WIDTH;  beta = beta * narrowing * (1 - PEAK_BROADENING * std.exp(-(lnr * lnr) / (2 * width)));  const sech = 1 / std.cosh(beta * wrapAngle(theta));  return (beta / (2 * std.tanh(beta * Math.PI))) * sech * sech;};const spreading = (r: number, theta: number, params: d.Infer<typeof SpectrumParams>) => {  'use gpu';  const lobe = donelanBanner(r, theta, std.mix(0.35, 1.1, params.directionality));  return std.mix(lobe, 0.5 / Math.PI, (1 - params.directionality) * ISOTROPIC_FLOOR);};const bandEdge = (k: number, edge: number) => {  'use gpu';  return std.smoothstep(edge * (1 - BAND_TAPER), edge * (1 + BAND_TAPER), k);};/** Expected |h0(k)|² of a lattice cell: a wind sea plus a swell from another direction. */export const waveVariance = (  kVec: d.v2f,  cascade: d.Infer<typeof Cascade>,  params: d.Infer<typeof SpectrumParams>,) => {  'use gpu';  const k = std.length(kVec);  if (k < 1e-4) {    return d.f32(0);  }  const disp = dispersion(k);  const kWorld = patchToWorld(kVec, cascade.rotation);  const angle = std.atan2(kWorld.y, kWorld.x);  const windSea =    jonswap(disp.x, params.peakFrequency, params.peakedness) *    params.energyScale *    spreading(disp.x / params.peakFrequency, angle - params.windDirection, params);  const swell =    jonswap(disp.x, params.swellFrequency, d.f32(SWELL_PEAKEDNESS)) *    params.swellScale *    donelanBanner(disp.x / params.swellFrequency, angle - params.swellDirection, SWELL_NARROWING);  const dk = (2 * Math.PI) / cascade.patchSize;  const band = bandEdge(k, cascade.kLow) * (1 - bandEdge(k, cascade.kHigh));  const gain = cascade.heightGain * cascade.heightGain;  return ((0.5 * (windSea + swell) * disp.y * dk * dk) / k) * band * gain;};/** Conjugate cell in a centred power-of-two FFT lattice. Integer keys avoid -0 on the axes. */export const mirrorCell = (cell: d.v2u) => {  'use gpu';  return d.vec2u((N - cell.x) & (N - 1), (N - cell.y) & (N - 1));};

surface.ts

import { cameras, environment, owner, type CameraBase, type Environment } from '@typegpu/kit';import { d, std, tgpu, type TgpuRenderCommands, type TgpuRoot } from 'typegpu';import { fresnelWater, meanFresnel, roughFresnel, smithLambda, sunGlint } from './shading.ts';import type { Waves } from './waves.ts';import type { Atmosphere } from './atmosphere.ts';import type { Clouds } from './clouds.ts';import type { Foam } from './foam.ts';import {  CASCADES,  N,  patchToWorld,  RIPPLES,  surfaceCovariance,  surfaceJacobian,  surfaceSlope,  tensorToWorld,  WIND,} from './spectrum.ts';const GRID = 128;const LEVELS = 11;const BASE_CELL = 0.25;const HOLE = GRID / 2 - 2;const MORPH = { start: GRID / 2 - 16, end: GRID / 2 - 4 };/** * After GodotOceanWaves (MIT, Ethan Truong) and the water BSDF of Atlas (GDC 2019): water and foam * albedos, and the green a wave scatters forward. */const WATER = d.vec3f(0.0100228, 0.0196066, 0.0272118);const FOAM = d.vec3f(0.75, 0.73, 0.7);const SCATTER = d.vec3f(0.9, 1.15, 0.85);/** The slopes' spread the body's Fresnel and masking see. */const SLOPES = 0.65;/** The body takes the sun's light by what enters the water, in full from a sun at 10°. */const SUN_FULL = Math.sin((10 * Math.PI) / 180);const ENTERING_FULL = SUN_FULL * (1 - fresnelWater(SUN_FULL));/** A wave lights through from this many metres below the still sea. */const SCATTER_DEPTH = 2.5;const FOAM_FADE = 0.0075;/** The share of the water foam hides, at most `most`, so the water shows through the lace. */const LACE = {  held: 0.45,  gathered: 0.4,  converging: 1,  gather: 3,  whitecap: 0.65,  fresh: 4,  most: 0.8,};/** Foam streaks upwind of a crest: `taps` samples `step` metres apart, each fading by `fade`. */const TRAIL = { taps: 4, step: 2, fade: 0.7 };/** Mist from `begin` to `end` metres, the light on the sea blended `aerial` to the sky. */const MIST = {  begin: 300,  end: 3000,  curve: 0.5,  density: 0.7,  tint: d.vec3f(0.029, 0.069, 0.092),  aerial: 0.626,  scatter: 0.05,};/** The share of the light's colour the eye adapts away; little at dusk, so dusk stays warm. */const ADAPTATION = { dusk: 0.7, day: 1, from: 0.035, to: 0.2 };/** The disc's peak per unit of sunlight, dim enough that the stratus over it shows. */const SUN_PEAK = 20;/** The deck hides the disc below this share of sky let through, and shows it from `clear`. */const DISC = { hidden: 0.35, clear: 0.9 };const LUMA = d.vec3f(0.2126, 0.7152, 0.0722);/** The ripples' share of the visible slope; the rest is roughness, so they don't shimmer. */const SURFACE_DETAIL = 0.7;const GLINT_CAP = 4;/** The share of the slopes' width the sky's reflection blurs over. */const REFLECTION_BLUR = 0.08;const RIPPLE_BIAS = 2;/** Seconds of swell motion the ripples' footprint spans: two frames, past which they alias. */const CARRIED = 2 / 60;/** Beyond this distance, wave detail filters into roughness. */const DETAIL_RANGE = 250;const ROUGHNESS = { near: 0.03, far: 0.1, ramp: 400 };/** Below the horizon the next wave, tilted skywards, mirrors `share` of the sky from `rise` up. */const NEXT_WAVE = { share: 0.25, rise: 0.35 };const Level = d.struct({ origin: d.vec2f, cellSize: d.f32 });const Slopes = d.struct({  slope: d.vec2f,  deformation: d.vec3f,  covariance: d.vec3f,  height: d.f32,});const gridLayout = tgpu.vertexLayout(d.arrayOf(d.vec2f));function clipmapMesh() {  const vertices: d.v2f[] = [];  for (let j = 0; j <= GRID; j++) {    for (let i = 0; i <= GRID; i++) vertices.push(d.vec2f(i - GRID / 2, j - GRID / 2));  }  const plate: number[] = [];  const ring: number[] = [];  const holeStart = (GRID - HOLE) / 2;  const holeEnd = (GRID + HOLE) / 2;  for (let j = 0; j < GRID; j++) {    for (let i = 0; i < GRID; i++) {      const a = j * (GRID + 1) + i;      const quad = [a, a + GRID + 1, a + 1, a + 1, a + GRID + 1, a + GRID + 2];      plate.push(...quad);      if (i < holeStart || i >= holeEnd || j < holeStart || j >= holeEnd) ring.push(...quad);    }  }  return { vertices, plate, ring };}/** Blend towards the next coarser level near a ring's outer edge. */const morphFactor = (grid: d.v2f) => {  'use gpu';  return std.smoothstep(MORPH.start, MORPH.end, std.max(std.abs(grid.x), std.abs(grid.y)));};/** Snaps odd vertices onto the coarser lattice so neighbouring levels share edges. */const morphGrid = (grid: d.v2f, morph: number) => {  'use gpu';  return grid - std.fract(grid * 0.5) * 2 * morph;};const roughnessFloor = (distance: number) => {  'use gpu';  const alpha = std.mix(ROUGHNESS.near, ROUGHNESS.far, std.smoothstep(0, ROUGHNESS.ramp, distance));  return 0.5 * alpha * alpha;};const skyAlong = (direction: d.v3f, roughness: number) => {  'use gpu';  const below = std.smoothstep(0, 0.05, -direction.y);  const height = std.mix(std.max(direction.y, 0.01), std.max(-direction.y, NEXT_WAVE.rise), below);  const seen = std.normalize(d.vec3f(direction.x, height, direction.z));  return environment.reflection(seen, roughness) * std.mix(1, NEXT_WAVE.share, below);};/** * The sea's surface: a clipmap displaced by the cascades, shaded with filtered slope covariance, * the water body, sun glint, foam and mist, then the sky where the sea leaves the depth clear. */export function createSurface(  root: TgpuRoot,  options: {    view: CameraBase;    waves: Waves;    atmosphere: Atmosphere;    clouds: Clouds;    foam: Foam;    light: Environment;    format: GPUTextureFormat;    depth: GPUTextureFormat;  },) {  const { view, waves, atmosphere, clouds, foam: lace, light, format, depth } = options;  const own = owner();  const mesh = clipmapMesh();  const grid = own(    root      .createBuffer(gridLayout.schemaForCount(mesh.vertices.length), mesh.vertices)      .$usage('vertex'),  );  const plate = own(    root.createBuffer(d.arrayOf(d.u32, mesh.plate.length), mesh.plate).$usage('index'),  );  const ring = own(    root.createBuffer(d.arrayOf(d.u32, mesh.ring.length), mesh.ring).$usage('index'),  );  const levels = root.createUniform(d.arrayOf(Level, LEVELS));  const placed = Array.from({ length: LEVELS }, (_, level) => ({    origin: d.vec2f(),    cellSize: BASE_CELL * 2 ** level,  }));  const vertex = tgpu.vertexFn({    in: { grid: d.vec2f, instance: d.builtin.instanceIndex },    out: { position: d.builtin.position, world: d.vec3f, rest: d.vec2f },  })((input) => {    'use gpu';    const level = levels.$[input.instance]!;    const morph = morphFactor(input.grid);    const rest = level.origin + morphGrid(input.grid, morph) * level.cellSize;    const footprint = level.cellSize * (1 + morph);    let offset = d.vec3f();    for (const c of tgpu.unroll(std.range(RIPPLES))) {      const cascade = waves.cascades.$[c]!;      // Cascades whose waves the cells cannot resolve fade out instead of aliasing      const weight =        1 - std.smoothstep(cascade.maxWavelength * 0.5, cascade.maxWavelength * 2, footprint);      const lod = std.max(std.log2((footprint * N) / cascade.patchSize), 0);      const sample = waves.displacementAt(rest, c, lod);      const horizontal = patchToWorld(sample.xz, cascade.rotation);      offset += d.vec3f(horizontal.x, sample.y, horizontal.y) * weight;    }    const world = d.vec3f(rest.x, 0, rest.y) + offset;    return {      position: cameras.current.$.viewProjection * d.vec4f(world, 1),      world,      rest,    };  });  /** Partly adapts the eye to the light falling on the sea from the sun and the sky above it. */  const balance = () => {    'use gpu';    const sun = atmosphere.sunDirection();    const falling =      atmosphere.sunTransmittance() * std.max(sun.y, 0.02) +      environment.irradiance(d.vec3f(0, 1, 0)) +      1e-6;    const adaptation = std.mix(      ADAPTATION.dusk,      ADAPTATION.day,      std.smoothstep(ADAPTATION.from, ADAPTATION.to, sun.y),    );    return std.pow(d.vec3f(std.dot(falling, LUMA)) / falling, d.vec3f(adaptation));  };  /**   * The cascades at a rest point, filtered over the pixel's world footprint, `dx` by `dy` and   * `wide` metres across: the resolved slope, the horizontal deformation, the slopes' covariance   * and the height of the waves the mesh resolves.   */  const slopesAt = (rest: d.v2f, dx: d.v2f, dy: d.v2f, wide: number, distance: number) => {    'use gpu';    let gradient = d.vec2f();    let deformation = d.vec3f();    let cov = d.vec3f();    let height = d.f32(0);    const carried = waves.spectrum.$.peakFrequency * waves.spectrum.$.waveHeight * (0.5 * CARRIED);    for (const c of tgpu.unroll(std.range(CASCADES))) {      const rotation = waves.cascades.$[c]!.rotation;      // Ripples finer than twice the pixel, or than the distance the swell carries them in a frame,      // move into the slope variance: the glint lobe widens steadily instead of flickering      const bias =        c === RIPPLES ? std.max(d.f32(RIPPLE_BIAS), carried / std.max(wide, 1e-6)) : d.f32(1);      const slopes =        c > WIND          ? waves.smoothDerivatives(rest, c, dx * bias, dy * bias)          : waves.derivatives(rest, c, dx, dy);      const displacement = waves.displacement(rest, c, dx, dy);      const moment =        c > WIND          ? waves.smoothMoments(rest, c, dx * bias, dy * bias)          : waves.moments(rest, c, dx, dy);      const local = d.vec3f(slopes.zw, displacement.w);      const detail = c === RIPPLES ? d.f32(SURFACE_DETAIL) : d.f32(1);      const visible = slopes.xy * detail;      gradient += patchToWorld(visible, rotation);      deformation += tensorToWorld(local, rotation) * detail;      cov += tensorToWorld(        d.vec3f(          moment.x - visible.x * visible.x,          moment.y - visible.y * visible.y,          moment.z - visible.x * visible.y,        ),        rotation,      );      if (c < RIPPLES) {        height += displacement.y;      }    }    // Filter linear gradients and moments first, then apply the nonlinear deformation    cov = surfaceCovariance(cov, deformation);    const unresolved = roughnessFloor(distance);    const covX = std.max(cov.x, 0) + unresolved;    const covY = std.max(cov.y, 0) + unresolved;    const limit = 0.99 * std.sqrt(covX * covY);    return Slopes({      slope: surfaceSlope(gradient, deformation),      deformation,      covariance: d.vec3f(covX, covY, std.clamp(cov.z, -limit, limit)),      height,    });  };  /**   * The foam's share of the pixel and its bubbles' brightness. The wind waves' coarse foam history   * marks broken crests. Within it the sea's sharpest compression makes fresh caps, and the rest   * gathers where the short waves converge.   */  const foamAt = (rest: d.v2f, dx: d.v2f, dy: d.v2f, footprint: number, deformation: d.v3f) => {    'use gpu';    const windAngle = waves.spectrum.$.windDirection;    const downwind = d.vec2f(std.cos(windAngle), std.sin(windAngle)) * TRAIL.step;    let breaking = waves.smoothMoments(rest, WIND, dx, dy).w;    for (let tap = 1; tap <= TRAIL.taps; tap++) {      const behind = waves.smoothMoments(rest + downwind * d.f32(tap), WIND, dx, dy).w;      breaking = std.max(breaking, behind * std.pow(TRAIL.fade, d.f32(tap)));    }    const jacobian = surfaceJacobian(deformation);    const fresh = std.saturate((LACE.whitecap - jacobian) * LACE.fresh);    const gathered = std.saturate((LACE.converging - jacobian) * LACE.gather);    const share =      std.saturate(breaking * (LACE.held + LACE.gathered * gathered) + fresh) * LACE.most;    return lace.layer(rest, share, 1 - breaking, footprint);  };  const fragment = tgpu.fragmentFn({    in: { world: d.vec3f, rest: d.vec2f },    out: d.vec4f,  })((input) => {    'use gpu';    const eye = cameras.current.$.position;    const distance = std.length(input.world - eye);    const ray = (input.world - eye) / distance;    const lodScale = std.exp2(std.max(std.log2(distance / DETAIL_RANGE), 0));    const dx = std.dpdx(input.rest);    const dy = std.dpdy(input.rest);    const footprint = std.max(std.length(dx), std.length(dy));    const lodDx = dx * lodScale;    const lodDy = dy * lodScale;    const sea = slopesAt(input.rest, lodDx, lodDy, footprint * lodScale, distance);    const cov = sea.covariance;    let normal = std.normalize(d.vec3f(-sea.slope.x, 1, -sea.slope.y));    const v = ray * -1;    const facing = std.dot(normal, v);    if (facing < 0.02) {      normal = std.normalize(normal - ray * (facing - 0.02));    }    const sun = atmosphere.sunDirection();    const sunColor = atmosphere.sunTransmittance() * clouds.shadow(input.world.xz);    const layer = foamAt(input.rest, lodDx, lodDy, footprint, sea.deformation);    const foam = layer.x * std.exp(-distance * FOAM_FADE);    const nv = std.max(std.dot(normal, v), 2e-5);    const nl = std.max(std.dot(normal, sun), 2e-5);    const fresnel = roughFresnel(nv, SLOPES);    // The sun scattered forward through raised waves, back towards the eye, and onto lit faces    const viewMask = smithLambda(SLOPES, nv);    const through =      std.max(sea.height + SCATTER_DEPTH, 0) *      std.pow(std.max(std.dot(sun, ray), 0), 4) *      std.pow(0.5 - 0.5 * std.dot(sun, normal), 3);    const back = 0.5 * nv * nv;    const lit = SCATTER * ((through + back) / (1 + viewMask)) + d.vec3f(0.5 * nl);    const sunIn = std.max(sun.y, 0);    const entering = std.saturate((sunIn * (1 - fresnelWater(sunIn))) / ENTERING_FULL);    // Sky low on the horizon glances off the surface, so the body sees mostly the sky above it    const ambient = environment.irradiance(std.normalize(normal + d.vec3f(0, 1, 0))) / Math.PI;    const water = sunColor * (WATER * lit * ((1 - fresnel) * entering)) + ambient * WATER;    // Foam lies on the surface, a rough heap of bubbles that shade each other from a glancing sun    const bubbles = FOAM * std.mix(0.55, 1.05, layer.y);    const froth = (sunColor * std.saturate(std.dot(normal, sun)) + ambient) * bubbles;    const body = std.mix(water, froth, foam);    const bitangent = std.normalize(std.cross(d.vec3f(1, 0, 0), normal));    const tangent = std.cross(normal, bitangent);    const vT = d.vec3f(std.dot(v, tangent), nv, std.dot(v, bitangent));    const lT = d.vec3f(std.dot(sun, tangent), std.dot(sun, normal), std.dot(sun, bitangent));    const width = std.sqrt(cov.x + cov.y);    // Sea slopes are Gaussian, and the cube's blur has long tails that would gather the bright    // band into every rough reflection: a narrower blur keeps the water mirroring what is above it    const reflection = skyAlong(      std.reflect(ray, normal),      std.sqrt(std.min(width * REFLECTION_BLUR, 1)),    );    const surface =      (reflection * meanFresnel(vT, cov) +        std.min(sunGlint(vT, lT, cov, sunColor), sunColor * GLINT_CAP)) *      (1 - foam);    const falling =      atmosphere.sunTransmittance() * std.max(sun.y, 0) + environment.irradiance(d.vec3f(0, 1, 0));    const mistColor =      std.mix(MIST.tint * falling, skyAlong(ray, 0.9), MIST.aerial) +      sunColor * (std.pow(std.max(std.dot(ray, sun), 0), 8) * MIST.scatter);    const mist = std.pow(std.smoothstep(MIST.begin, MIST.end, distance), MIST.curve) * MIST.density;    const color = std.mix(body + surface, mistColor, mist);    return d.vec4f(atmosphere.aerialPerspective(color, ray, distance) * balance(), 1);  });  // At the far plane, unlike `common.fullScreenTriangle`, so the sky fills only what the sea left  const fullscreen = tgpu.vertexFn({    in: { index: d.builtin.vertexIndex },    out: { position: d.builtin.position, ndc: d.vec2f },  })(({ index }) => {    'use gpu';    const ndc = d.vec2f(d.f32((index << 1) & 2), d.f32(index & 2)) * 2 - 1;    return { position: d.vec4f(ndc, 1, 1), ndc };  });  const background = tgpu.fragmentFn({ in: { ndc: d.vec2f }, out: d.vec4f })(({ ndc }) => {    'use gpu';    const ray = cameras.rayOf(cameras.current.$, ndc).direction;    const deck = clouds.seen(d.vec2f(ndc.x, -ndc.y) * 0.5 + 0.5);    // The disc is capped so its bloom keeps one size as the sun climbs and brightens    const disc = atmosphere.sunRadiance(ray);    const peak = std.max(disc.x, std.max(disc.y, disc.z));    const shown = std.smoothstep(DISC.hidden, DISC.clear, deck.w);    const sky =      atmosphere.skyRadiance(ray) + disc * (std.min(1, SUN_PEAK / std.max(peak, 1e-6)) * shown);    return d.vec4f((sky * deck.w + deck.xyz) * balance(), 1);  });  const configured = root.pipe(view.inject()).pipe(light.inject());  const ocean = configured.createRenderPipeline({    attribs: { grid: gridLayout.attrib },    vertex,    fragment,    targets: { format },    primitive: { cullMode: 'none' },    depthStencil: { format: depth, depthWriteEnabled: true, depthCompare: 'less' },  });  const behind = configured.createRenderPipeline({    vertex: fullscreen,    fragment: background,    targets: { format },    depthStencil: { format: depth, depthWriteEnabled: false, depthCompare: 'less-equal' },  });  return {    pipelines: [ocean, behind],    /** Centres the clipmap's levels under the eye; each snaps to twice its cell. */    place(x: number, z: number) {      for (const level of placed) {        const snap = 2 * level.cellSize;        level.origin = d.vec2f(          Math.floor(x / snap + 0.5) * snap,          Math.floor(z / snap + 0.5) * snap,        );      }      levels.write(placed);    },    /** Draws the sea, then the sky where the sea left the depth clear. */    draw(pass: TgpuRenderCommands) {      const drawn = ocean.with(pass).with(gridLayout, grid);      drawn.withIndexBuffer(plate).drawIndexed(mesh.plate.length, 1);      drawn.withIndexBuffer(ring).drawIndexed(mesh.ring.length, LEVELS - 1, 0, 0, 1);      behind.with(pass).with(clouds.group).draw(3);    },    destroy: own.destroy,  };}

waves.ts

import { owner } from '@typegpu/kit';import { randf } from '@typegpu/noise';import { d, std, tgpu, type TgpuComputePass, type TgpuRoot } from 'typegpu';import { createWaveMips } from './wave-mips.ts';import {  CASCADES,  Cascade,  cascadeValues,  dispersion,  Foaming,  LOG2N,  LOOP_PERIOD,  MIP_LEVELS,  mirrorCell,  N,  REACH,  SIGNALS,  SpectrumParams,  surfaceJacobian,  waveVariance,  worldToPatch,} from './spectrum.ts';const TILE = 16;const EVOLVE_GROUP = 256;const HALF_N = N / 2;const OMEGA_0 = (2 * Math.PI) / LOOP_PERIOD;/** How far the waves' crests lean and sharpen, as a share of the linear waves' motion. */const CHOPPINESS = 1.3;const STRIDES = Array.from({ length: LOG2N }, (_, stage) => 1 << stage);const SIGNAL_INDICES = Array.from({ length: SIGNALS }, (_, i) => i);/** Kernels visit only the cells within `REACH`, listed across cascades from these starts. */const SIDES = REACH.map((reach) => 2 * reach + 1);const startsOf = (counts: number[]) =>  counts.map((_, c) => counts.slice(0, c).reduce((a, b) => a + b, 0));const ROW_STARTS = startsOf(SIDES);const ROWS = ROW_STARTS[CASCADES - 1]! + SIDES[CASCADES - 1]!;const CELL_STARTS = startsOf(SIDES.map((side) => side * side));const CELLS = CELL_STARTS[CASCADES - 1]! + SIDES[CASCADES - 1]! ** 2;const reachOf = tgpu.const(d.arrayOf(d.u32, CASCADES), REACH);const rowStarts = tgpu.const(d.arrayOf(d.u32, CASCADES), ROW_STARTS);const cellStarts = tgpu.const(d.arrayOf(d.u32, CASCADES), CELL_STARTS);const cascadeAt = (starts: number[]) => {  const later = starts.slice(1);  return (index: number) => {    'use gpu';    let c = d.u32(0);    for (const start of tgpu.unroll(later)) {      if (index >= start) {        c += 1;      }    }    return c;  };};const cascadeOfRow = cascadeAt(ROW_STARTS);const cascadeOfCell = cascadeAt(CELL_STARTS);const Step = d.struct({ time: d.f32, dt: d.f32 });const rowsLayout = tgpu.bindGroupLayout({  data: { storage: d.arrayOf(d.vec2f), access: 'mutable' },});const field = d.textureStorage2dArray('rgba16float', 'write-only');const columnsLayout = tgpu.bindGroupLayout({  rows: { storage: d.arrayOf(d.vec2f), access: 'readonly' },  displacement: { storageTexture: field },  derivatives: { storageTexture: field },  moments: { storageTexture: field },});const cmul = (a: d.v2f, b: d.v2f) => {  'use gpu';  return d.vec2f(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);};/** a + i·b for two spectra whose inverse transforms are real. */const pack = (a: d.v2f, b: d.v2f) => {  'use gpu';  return d.vec2f(a.x - b.y, a.y + b.x);};/** Where thread `t` loads its two samples, `t` and `t` + N/2, in bit-reversed order. */const reversedOf = (t: number) => {  'use gpu';  return std.reverseBits(t) >>> (32 - LOG2N);};/** Inverse-transforms `count` signals of `lines` in place: bit-reversed in, natural order out. */const butterfliesOf = (lines: { $: d.v2f[] }, count: number) => {  const signals = SIGNAL_INDICES.slice(0, count);  return (t: number) => {    'use gpu';    std.workgroupBarrier();    for (const stride of tgpu.unroll(STRIDES)) {      const j = t & (stride - 1);      const first = (t - j) * 2 + j;      const angle = (Math.PI / stride) * d.f32(j);      const twiddle = d.vec2f(std.cos(angle), std.sin(angle));      for (const signal of tgpu.unroll(signals)) {        const at = signal * N + first;        const a = d.vec2f(lines.$[at]!);        const b = cmul(lines.$[at + stride]!, twiddle);        lines.$[at] = a + b;        lines.$[at + stride] = a - b;      }      std.workgroupBarrier();    }  };};const row = tgpu.workgroupVar(d.arrayOf(d.vec2f, N));const rowButterflies = butterfliesOf(row, 1);/** One workgroup inverse-transforms one row within reach; samples beyond the reach are zero. */const transformRows = tgpu.computeFn({  workgroupSize: [HALF_N],  in: { lid: d.builtin.localInvocationId, wid: d.builtin.workgroupId },})(({ lid, wid }) => {  'use gpu';  const t = lid.x;  const cascade = cascadeOfRow(wid.x);  const reach = reachOf.$[cascade]!;  const base =    ((cascade * SIGNALS + wid.y) * N + HALF_N - reach + wid.x - rowStarts.$[cascade]!) * N;  const reversed = reversedOf(t);  row.$[reversed] = d.vec2f();  row.$[reversed + 1] = d.vec2f();  if (t + reach >= HALF_N) {    row.$[reversed] = d.vec2f(rowsLayout.$.data[base + t]!);  }  if (t <= reach) {    row.$[reversed + 1] = d.vec2f(rowsLayout.$.data[base + t + HALF_N]!);  }  rowButterflies(t);  rowsLayout.$.data[base + t] = d.vec2f(row.$[t]!);  rowsLayout.$.data[base + t + HALF_N] = d.vec2f(row.$[t + HALF_N]!);});const columns = tgpu.workgroupVar(d.arrayOf(d.vec2f, SIGNALS * N));const columnButterflies = butterfliesOf(columns, SIGNALS);export type Waves = ReturnType<typeof createWaves>;/** Four FFT cascades of an ocean spectrum, mipmapped as displacement, slope and moment maps. */export function createWaves(root: TgpuRoot) {  const own = owner();  const spectrum = root.createUniform(SpectrumParams);  const cascades = root.createUniform(d.arrayOf(Cascade, CASCADES), cascadeValues);  const step = root.createUniform(Step, { time: 0, dt: 0 });  const foaming = root.createUniform(d.arrayOf(Foaming, CASCADES));  const foam = own(root.createMutable(d.arrayOf(d.f32, CASCADES * N * N)));  const amplitudes = own(root.createMutable(d.arrayOf(d.vec4f, CASCADES * N * N)));  const omegas = own(root.createMutable(d.arrayOf(d.f32, CASCADES * N * N)));  const lattice = own(    root.createBuffer(d.arrayOf(d.vec2f, CASCADES * SIGNALS * N * N)).$usage('storage'),  );  const signals = lattice.as('mutable');  const sampler = root.createSampler({    magFilter: 'linear',    minFilter: 'linear',    mipmapFilter: 'linear',    addressModeU: 'repeat',    addressModeV: 'repeat',    maxAnisotropy: 8,  });  const createMap = () =>    own(      root        .createTexture({ size: [N, N, CASCADES], format: 'rgba16float', mipLevelCount: MIP_LEVELS })        .$usage('sampled', 'storage'),    );  const maps = [createMap(), createMap(), createMap()] as const;  const displacement = maps[0].createView(d.texture2dArray(d.f32));  const derivatives = maps[1].createView(d.texture2dArray(d.f32));  const moments = maps[2].createView(d.texture2dArray(d.f32));  const mips = createWaveMips(root, maps);  const top = { baseMipLevel: 0, mipLevelCount: 1 };  /** h0(k), conj(h0(-k)) and the looped frequency of every lattice cell of every cascade. */  const initialSpectrum = tgpu.computeFn({    workgroupSize: [TILE, TILE],    in: { gid: d.builtin.globalInvocationId },  })(({ gid }) => {    'use gpu';    if (gid.x >= N || gid.y >= N) return;    const cascade = Cascade(cascades.$[gid.z]!);    const params = SpectrumParams(spectrum.$);    const dk = (2 * Math.PI) / cascade.patchSize;    const kVec = (d.vec2f(gid.xy) - N / 2) * dk;    const seed = params.seed + d.f32(gid.z) * 17;    const valid = std.select(d.f32(1), d.f32(0), gid.x === 0 || gid.y === 0);    randf.seed3(d.vec3f(d.vec2f(gid.xy), seed));    const h0 =      d.vec2f(randf.normal(0, 1), randf.normal(0, 1)) *      std.sqrt(0.5 * waveVariance(kVec, cascade, params) * valid);    randf.seed3(d.vec3f(d.vec2f(mirrorCell(gid.xy)), seed));    const h0Mirror =      d.vec2f(randf.normal(0, 1), randf.normal(0, 1)) *      std.sqrt(0.5 * waveVariance(kVec * -1, cascade, params) * valid);    const idx = (gid.z * N + gid.y) * N + gid.x;    amplitudes.$[idx] = d.vec4f(h0, h0Mirror.x, -h0Mirror.y);    const omega = dispersion(std.max(std.length(kVec), 1e-6)).x;    const looped = std.floor(omega / OMEGA_0) * OMEGA_0;    omegas.$[idx] = looped;  });  const evolve = tgpu.computeFn({    workgroupSize: [EVOLVE_GROUP],    in: { gid: d.builtin.globalInvocationId },  })(({ gid }) => {    'use gpu';    if (gid.x >= CELLS) return;    const c = cascadeOfCell(gid.x);    const reach = reachOf.$[c]!;    const side = reach * 2 + 1;    const local = gid.x - cellStarts.$[c]!;    const cell = d.vec2u(local % side, std.intdiv(local, side)) + (HALF_N - reach);    const idx = (c * N + cell.y) * N + cell.x;    const h0 = amplitudes.$[idx]!;    const phase = omegas.$[idx]! * step.$.time;    const rotation = d.vec2f(std.cos(phase), std.sin(phase));    const h = cmul(h0.xy, rotation) + cmul(h0.zw, d.vec2f(rotation.x, -rotation.y));    const dk = (2 * Math.PI) / cascades.$[c]!.patchSize;    const k = (d.vec2f(cell) - N / 2) * dk;    const kLen = std.max(std.length(k), 1e-6);    const ih = d.vec2f(-h.y, h.x);    const base = c * SIGNALS * N * N + cell.y * N + cell.x;    // Displacement x, height; displacement z, slope x; slope z, ∂Dx/∂x; ∂Dz/∂z, ∂Dx/∂z    signals.$[base] = pack(ih * (-k.x / kLen), h);    signals.$[base + N * N] = pack(ih * (-k.y / kLen), ih * k.x);    signals.$[base + 2 * N * N] = pack(ih * k.y, h * ((k.x * k.x) / kLen));    signals.$[base + 3 * N * N] = pack(h * ((k.y * k.y) / kLen), h * ((k.x * k.y) / kLen));  });  /** Packs a transformed texel's displacement, slopes and moments, and carries its foam. */  const packTexel = (texel: d.v2u, c: number) => {    'use gpu';    const sign = std.select(d.f32(1), d.f32(-1), ((texel.x + texel.y) & 1) === 1);    const s0 = d.vec2f(columns.$[texel.y]!) * sign;    const s1 = d.vec2f(columns.$[N + texel.y]!) * sign;    const s2 = d.vec2f(columns.$[2 * N + texel.y]!) * sign;    const s3 = d.vec2f(columns.$[3 * N + texel.y]!) * sign;    // With -i·k/|k|, a negative gain compresses the crests; a positive one sharpens the troughs    const lambda = -CHOPPINESS * cascades.$[c]!.chop;    const gradient = d.vec2f(s1.y, s2.x);    const xx = lambda * s2.y;    const zz = lambda * s3.x;    const xz = lambda * s3.y;    const jacobian = surfaceJacobian(d.vec3f(xx, zz, xz));    const rule = Foaming(foaming.$[c]!);    const cell = (c * N + texel.x) * N + texel.y;    const foamed = std.saturate(      foam.$[cell]! * std.exp(-rule.decay * step.$.dt) +        std.max(rule.whitecap - jacobian, 0) * rule.grow * step.$.dt,    );    foam.$[cell] = foamed;    std.textureStore(      columnsLayout.$.displacement,      texel,      c,      d.vec4f(lambda * s0.x, s0.y, lambda * s1.x, xz),    );    std.textureStore(columnsLayout.$.derivatives, texel, c, d.vec4f(gradient, xx, zz));    std.textureStore(      columnsLayout.$.moments,      texel,      c,      d.vec4f(gradient.x * gradient.x, gradient.y * gradient.y, gradient.x * gradient.y, foamed),    );  };  /** One workgroup inverse-transforms one column of a cascade's signals and packs its texels. */  const transformColumns = tgpu.computeFn({    workgroupSize: [HALF_N],    in: { lid: d.builtin.localInvocationId, wid: d.builtin.workgroupId },  })(({ lid, wid }) => {    'use gpu';    const t = lid.x;    const x = wid.x;    const c = wid.y;    const reach = reachOf.$[c]!;    const reversed = reversedOf(t);    for (const signal of tgpu.unroll(SIGNAL_INDICES)) {      const column = (c * SIGNALS + signal) * N * N + x;      let low = d.vec2f();      let high = d.vec2f();      if (t + reach >= HALF_N) {        low = d.vec2f(columnsLayout.$.rows[column + t * N]!);      }      if (t <= reach) {        high = d.vec2f(columnsLayout.$.rows[column + (t + HALF_N) * N]!);      }      columns.$[signal * N + reversed] = d.vec2f(low);      columns.$[signal * N + reversed + 1] = d.vec2f(high);    }    columnButterflies(t);    for (const half of tgpu.unroll([0, HALF_N])) {      packTexel(d.vec2u(x, t + half), c);    }  });  const initialPipeline = root.createComputePipeline({ compute: initialSpectrum });  const evolvePipeline = root.createComputePipeline({ compute: evolve });  const rowsPipeline = root.createComputePipeline({ compute: transformRows });  const columnsPipeline = root.createComputePipeline({ compute: transformColumns });  const rowsGroup = root.createBindGroup(rowsLayout, { data: lattice });  const columnsGroup = root.createBindGroup(columnsLayout, {    rows: lattice,    displacement: maps[0].createView(field, top),    derivatives: maps[1].createView(field, top),    moments: maps[2].createView(field, top),  });  const patchUv = (xz: d.v2f, c: number) => {    'use gpu';    return worldToPatch(xz, cascades.$[c]!.rotation) / cascades.$[c]!.patchSize;  };  /** Samples a wave field filtered over a pixel's world footprint, `dx` by `dy`. */  const filtered = (map: d.texture2dArray<d.F32>, xz: d.v2f, c: number, dx: d.v2f, dy: d.v2f) => {    'use gpu';    const uv = patchUv(xz, c);    const cascade = cascades.$[c]!;    const uvDx = worldToPatch(dx, cascade.rotation) / cascade.patchSize;    const uvDy = worldToPatch(dy, cascade.rotation) / cascade.patchSize;    return std.textureSampleGrad(map, sampler.$, uv, c, uvDx, uvDy);  };  /** Cubic B-spline through four bilinear taps, for fields magnified past their texels. */  const smoothField = (map: d.texture2dArray<d.F32>, xz: d.v2f, c: number) => {    'use gpu';    const texel = patchUv(xz, c) * N - 0.5;    const base = std.floor(texel);    const f = texel - base;    const f2 = f * f;    const f3 = f2 * f;    const w0 = (f2 * 3 - f3 - f * 3 + 1) / 6;    const w1 = (f3 * 3 - f2 * 6 + 4) / 6;    const w2 = (f2 * 3 - f3 * 3 + f * 3 + 1) / 6;    const w3 = f3 / 6;    const g0 = w0 + w1;    const g1 = w2 + w3;    const at0 = (base + 0.5 + w1 / g0 - 1) / N;    const at1 = (base + 0.5 + w3 / g1 + 1) / N;    const a = std.textureSampleLevel(map, sampler.$, d.vec2f(at0.x, at0.y), c, 0);    const b = std.textureSampleLevel(map, sampler.$, d.vec2f(at1.x, at0.y), c, 0);    const e = std.textureSampleLevel(map, sampler.$, d.vec2f(at0.x, at1.y), c, 0);    const h = std.textureSampleLevel(map, sampler.$, d.vec2f(at1.x, at1.y), c, 0);    return (a * g0.x + b * g1.x) * g0.y + (e * g0.x + h * g1.x) * g1.y;  };  const smoothed = (map: d.texture2dArray<d.F32>, xz: d.v2f, c: number, dx: d.v2f, dy: d.v2f) => {    'use gpu';    if (std.max(std.length(dx), std.length(dy)) < cascades.$[c]!.patchSize / N) {      return smoothField(map, xz, c);    }    return filtered(map, xz, c, dx, dy);  };  return {    pipelines: [initialPipeline, evolvePipeline, rowsPipeline, columnsPipeline, mips.pipeline],    spectrum,    cascades,    step,    foaming,    /** The displacement at mip level `lod`, for kernels and vertices. */    displacementAt: (xz: d.v2f, c: number, lod: number) => {      'use gpu';      return std.textureSampleLevel(displacement.$, sampler.$, patchUv(xz, c), c, lod);    },    /** The fields filtered over a pixel's world footprint, `dx` by `dy`, for fragments. */    displacement: (xz: d.v2f, c: number, dx: d.v2f, dy: d.v2f) => {      'use gpu';      return filtered(displacement.$, xz, c, dx, dy);    },    derivatives: (xz: d.v2f, c: number, dx: d.v2f, dy: d.v2f) => {      'use gpu';      return filtered(derivatives.$, xz, c, dx, dy);    },    moments: (xz: d.v2f, c: number, dx: d.v2f, dy: d.v2f) => {      'use gpu';      return filtered(moments.$, xz, c, dx, dy);    },    /**     * Like `derivatives` and `moments`, but smooth where a pixel is smaller than the cascade's     * texel. Use them together: the same positive weights keep the slopes' variance non-negative.     */    smoothDerivatives: (xz: d.v2f, c: number, dx: d.v2f, dy: d.v2f) => {      'use gpu';      return smoothed(derivatives.$, xz, c, dx, dy);    },    smoothMoments: (xz: d.v2f, c: number, dx: d.v2f, dy: d.v2f) => {      'use gpu';      return smoothed(moments.$, xz, c, dx, dy);    },    /** Records the initial amplitudes of the spectrum `spectrum` holds; write it first. */    rebuild(pass: TgpuComputePass) {      initialPipeline.with(pass).dispatchWorkgroups(N / TILE, N / TILE, CASCADES);    },    record(pass: TgpuComputePass) {      evolvePipeline.with(pass).dispatchWorkgroups(Math.ceil(CELLS / EVOLVE_GROUP));      rowsPipeline.with(pass).with(rowsGroup).dispatchWorkgroups(ROWS, SIGNALS);      columnsPipeline.with(pass).with(columnsGroup).dispatchWorkgroups(N, CASCADES);      mips.record(pass);    },    destroy: own.destroy,  };}

shading.ts

import { d, std } from 'typegpu';const WATER_IOR = 1.333;/** Below this cosine the lobe is subpixel; keeps intermediate maths in f32 range. */const MIN_COSINE = 1e-6;/** Unpolarized air-to-water Fresnel reflectance. */export const fresnelWater = (cosI: number) => {  'use gpu';  const sinT2 = (1 - cosI * cosI) / (WATER_IOR * WATER_IOR);  const cosT = std.sqrt(std.max(1 - sinT2, 0));  const rs = (cosI - WATER_IOR * cosT) / (cosI + WATER_IOR * cosT);  const rp = (WATER_IOR * cosI - cosT) / (WATER_IOR * cosI + cosT);  return std.saturate(0.5 * (rs * rs + rp * rp));};/** Smith masking for a Beckmann slope distribution, Walter's rational fit. */const beckmannLambda = (a: number) => {  'use gpu';  if (a >= 1.6) {    return d.f32(0);  }  return (1 - 1.259 * a + 0.396 * a * a) / (3.535 * a + 2.181 * a * a);};/** Slope variance of a covariance packed as (xx, yy, xy) along the azimuth of v, times sin². */const projectedVariance = (v: d.v3f, cov: d.v3f) => {  'use gpu';  return cov.x * v.x * v.x + 2 * cov.z * v.x * v.z + cov.y * v.z * v.z;};const smithParam = (v: d.v3f, cov: d.v3f) => {  'use gpu';  return std.max(v.y, MIN_COSINE) / std.sqrt(2 * std.max(projectedVariance(v, cov), 0) + 1e-8);};const slopeDensity = (zeta: d.v2f, cov: d.v3f) => {  'use gpu';  const det = std.max(cov.x * cov.y - cov.z * cov.z, 1e-8);  const quad =    (cov.y * zeta.x * zeta.x - 2 * cov.z * zeta.x * zeta.y + cov.x * zeta.y * zeta.y) / det;  // Rounded moments can make this slightly negative; negligible tails skip subnormal maths  if (quad > 160) {    return d.f32(0);  }  return std.exp(-0.5 * std.max(quad, 0)) / (2 * Math.PI * std.sqrt(det));};/** Sun glint over Gaussian slopes, with vectors in the resolved normal's frame: x, up, z. */export const sunGlint = (v: d.v3f, l: d.v3f, cov: d.v3f, sun: d.v3f) => {  'use gpu';  if (v.y <= 0 || l.y <= 0) {    return d.vec3f();  }  const halfway = v + l;  const h = halfway / std.max(std.length(halfway), MIN_COSINE);  const cosH = std.max(h.y, MIN_COSINE);  const p = slopeDensity(d.vec2f(-h.x / cosH, -h.z / cosH), cov);  const h4 = cosH * cosH * cosH * cosH;  const masking = 1 + beckmannLambda(smithParam(v, cov)) + beckmannLambda(smithParam(l, cov));  return sun * ((p * fresnelWater(std.dot(v, h))) / (4 * h4 * std.max(v.y, MIN_COSINE) * masking));};/** Fresnel over slopes of deviation `sigma` at `cosine`, Bruneton's fit, as GodotOceanWaves. */export const roughFresnel = (cosine: number, sigma: number) => {  'use gpu';  const oneMinus = 1 - std.saturate(cosine);  return (    0.02 + 0.98 * (std.pow(oneMinus, 5 * std.exp(-2.69 * sigma)) / (1 + 22.7 * std.pow(sigma, 1.5)))  );};/** Smith's Λ for slopes of width `alpha` seen at `cosine` from the mean normal. */export const smithLambda = (cosine: number, alpha: number) => {  'use gpu';  return beckmannLambda(cosine / (alpha * std.sqrt(std.max(1 - cosine * cosine, 1e-6))));};/** Fresnel averaged over the slope distribution seen from v, Bruneton's fit. */export const meanFresnel = (v: d.v3f, cov: d.v3f) => {  'use gpu';  const horizontal = v.x * v.x + v.z * v.z + 1e-8;  const sigmaV = std.sqrt(std.max(projectedVariance(v, cov), 0) / horizontal);  const oneMinus = 1 - v.y;  const fifth = oneMinus * oneMinus * oneMinus * oneMinus * oneMinus;  return 0.02 + 0.98 * ((fifth * std.exp(-2.69 * sigmaV)) / (1 + 22.7 * std.pow(sigmaV, 1.5)));};

wave-mips.ts

import {  d,  std,  tgpu,  type SampledFlag,  type StorageFlag,  type TgpuComputePass,  type TgpuRoot,  type TgpuTexture,} from 'typegpu';const GROUP = 16;const storage = d.textureStorage2dArray('rgba16float', 'write-only');const layout = tgpu.bindGroupLayout({  src: { texture: d.texture2dArray(d.f32) },  first: { storageTexture: storage },  second: { storageTexture: storage },  third: { storageTexture: storage },});const held = tgpu.workgroupVar(d.arrayOf(d.vec4f, GROUP * GROUP));/** Rounds half floats towards zero, below the smallest normal in steps of 2⁻²⁴. */const towardZero = (value: number) => {  'use gpu';  if (std.abs(value) < 2 ** -14) {    return std.trunc(value * 2 ** 24) / 2 ** 24;  }  return std.bitcastU32toF32(std.bitcastF32toU32(value) & 0xffffe000);};/** Rounds as `rgba16float` stores do on Apple GPUs, so held and stored averages match. */const asStored = (value: d.v4f) => {  'use gpu';  return d.vec4f(    towardZero(value.x),    towardZero(value.y),    towardZero(value.z),    towardZero(value.w),  );};const averageHeld = (index: number, stride: number) => {  'use gpu';  const sum =    held.$[index]! +    held.$[index + stride]! +    held.$[index + stride * GROUP]! +    held.$[index + stride * GROUP + stride]!;  return asStored(sum * 0.25);};/** Each thread averages a texel of the first level, and the workgroup those into two more. */const fillThree = tgpu.computeFn({  workgroupSize: [GROUP, GROUP],  in: { gid: d.builtin.globalInvocationId, lid: d.builtin.localInvocationId },})(({ gid, lid }) => {  'use gpu';  const index = lid.y * GROUP + lid.x;  let value = d.vec4f();  if (std.all(std.lt(gid.xy, std.textureDimensions(layout.$.first)))) {    const at = gid.xy * 2;    const sum =      std.textureLoad(layout.$.src, at, gid.z, 0) +      std.textureLoad(layout.$.src, at + d.vec2u(1, 0), gid.z, 0) +      std.textureLoad(layout.$.src, at + d.vec2u(0, 1), gid.z, 0) +      std.textureLoad(layout.$.src, at + d.vec2u(1, 1), gid.z, 0);    value = asStored(sum * 0.25);    std.textureStore(layout.$.first, gid.xy, gid.z, value);  }  held.$[index] = d.vec4f(value);  std.workgroupBarrier();  const half = d.vec2u(gid.x >>> 1, gid.y >>> 1);  const second = (lid.x & 1) === 0 && (lid.y & 1) === 0;  if (second && std.all(std.lt(half, std.textureDimensions(layout.$.second)))) {    const average = averageHeld(index, d.u32(1));    held.$[index] = d.vec4f(average);    std.textureStore(layout.$.second, half, gid.z, average);  }  std.workgroupBarrier();  const quarter = d.vec2u(gid.x >>> 2, gid.y >>> 2);  const third = (lid.x & 3) === 0 && (lid.y & 3) === 0;  if (third && std.all(std.lt(quarter, std.textureDimensions(layout.$.third)))) {    std.textureStore(layout.$.third, quarter, gid.z, averageHeld(index, d.u32(2)));  }});type WaveMap = TgpuTexture<{  size: readonly number[];  format: 'rgba16float';  mipLevelCount?: number | undefined;}> &  SampledFlag &  StorageFlag;const only = (baseMipLevel: number) => ({ baseMipLevel, mipLevelCount: 1 });/** * Box-filtered mips of the wave maps, three levels a dispatch. The maps refill every frame, and the * kit's `mipmaps`, a level a dispatch, waits between three times as many. */export function createWaveMips(root: TgpuRoot, maps: readonly WaveMap[]) {  const pipeline = root.createComputePipeline({ compute: fillThree });  const dispatches = maps.flatMap((map) => {    const { size, mipLevelCount = 1 } = map.props;    return Array.from({ length: Math.floor((mipLevelCount - 1) / 3) }, (_, batch) => {      const top = batch * 3;      const group = root.createBindGroup(layout, {        src: map.createView(d.texture2dArray(d.f32), only(top)),        first: map.createView(storage, only(top + 1)),        second: map.createView(storage, only(top + 2)),        third: map.createView(storage, only(top + 3)),      });      const x = Math.ceil(Math.max(size[0]! >> (top + 1), 1) / GROUP);      const y = Math.ceil(Math.max(size[1]! >> (top + 1), 1) / GROUP);      return (pass: TgpuComputePass) =>        pipeline.with(pass).with(group).dispatchWorkgroups(x, y, size[2]);    });  });  return {    pipeline,    record(pass: TgpuComputePass) {      for (const dispatch of dispatches) dispatch(pass);    },  };}