Snell's Window: The Ocean Seen from Below

When you look upward from underwater, the entire 180° sky is compressed into a cone with half-angle $\arcsin(n_a/n_w)$. Outside this cone relative to the local normal, light coming from the air cannot reach the eye, and the interface becomes a mirror through total internal reflection. A rippling water surface bends and distorts that boundary, but does not change the local critical-angle law itself.

Physical Model: Measured Constants and Explicit Model Coefficients

This implementation is intentionally configured as a clear-water hybrid: it uses the specified seawater conditions for the interface refractive index, but uses measured pure-water absorption and molecular scattering for spectral attenuation. Dissolved substances, phytoplankton, and suspended-particle scattering from any particular ocean are not included.

The initial AI implementation used $a_B=0.015\,\mathrm{m^{-1}}$ for blue absorption and also a scalar B_BACK=0.0025. The first value is much closer to cyan wavelengths longer than 440 nm, and the second was named as a backscattering coefficient but was actually used as a total extinction value. Both have been corrected. Beam extinction now follows the explicit relation $\mathbf c=\mathbf a+\mathbf b$.

Mathematical Notes: Exact Checks and Range-Limited Approximations

The entire hemisphere fits exactly inside the cone. Sampling sky zenith angles from 0–90° and refracting each into water, the maximum reachable angle matches the critical angle to within machine precision.

The bright rim is not drawn. The last 5° of sky just before the horizon occupies 0.51% of the angular radius. The full 360° horizon panorama is squeezed into a thin ring. The radius in a perspective view depends on the projection, so 0.51% should not be presented as a universal pixel-width ratio.

Fresnel. The numerically computed $R_0$ agrees with the analytic formula $[(n_w-n_a)/(n_w+n_a)]^2=0.02105$, and $R\to1$ holds at the critical angle.

Interface radiance. The Fresnel power ratio satisfies $R+T=1$, but that alone does not account for the full radiance conversion. Transmitted skylight also obeys the $L/n^2$ invariant, so air-side radiance is multiplied by $n_w^2$ when it passes into the water.

The caustic focal distance checks out. Surface ripples with amplitude a and wavenumber k focus refracted sunlight at f = 1/(0.25339·a·k²). The swell used here (λ = 0.8–1.8 m, a = 8–30 mm) focuses at 7.7–10.8 m, consistent with the distance between the surface and the seafloor. Without these ripples, swell alone would focus beyond 100 m, producing no caustics at all.

Caustic estimator. The gain uses a finite-difference Jacobian of the refracted-sunlight landing map. On a flat surface the raw value is 1, and curvature focuses light. The displayed value adds fold clamping and empirical normalization on top, so it need not equal 1 in flat regions. It is noted explicitly that FLOOR_NORM and BEAM_NORM are empirical correction factors, not proofs of global energy conservation.

Shortcuts Tried and Discarded

In the initial AI implementation, the caustic Jacobian was implemented as a monolithic block. The caustic Jacobian is expensive because it evaluates the surface normal three times per sample. An analytic version using the Hessian of the wave field and the normal-incidence factor $1-1/n$ was cheaper, but it diverged significantly from the finite-difference forward map because the base refracted sunlight enters at an oblique angle.

The current estimator starts with an initial guess inverted from the mean plane, evaluates only a single local branch, clamps fold singularities, and applies offline normalization. It is not an exact inverse caustic solution, nor does it sum all preimages at folds.

The ray estimator intentionally uses a wide stencil and two octaves of surface. Octaves smaller than that stencil were excluded as an explicit band-limiting approximation. Without reproducible sampling artifacts to compare against, this choice is not claimed to be a verified physical equivalent.

From Equations to Production Code

The equations below are paired with the code currently running. The visual model is owned by shaders/snells-window.frag, while browser-side auditing, loading, interaction, and audio state are owned by main.js.

Each runtime file has exactly one owner.

• shaders/vertex.vert — full-screen WebGL vertex stage

• shaders/snells-window.frag — all optics, surface, caustics, and water GLSL

• main.js — explicit shader URLs, loading, compilation, controls, and auditing

• underwatersound.mp3 — optional original audio recording

There is no inline shader fallback path. If the GLSL file is missing, an error is raised visibly rather than silently falling back to an old JavaScript copy as a second source of truth.

const SHADER_URLS = Object.freeze({
    vertex: new URL('./shaders/vertex.vert', import.meta.url),
    fragment: new URL('./shaders/snells-window.frag', import.meta.url)
});

1. Snell's Cone and the Critical Angle

Let $n_a$ be the refractive index of air and $n_w$ that of water; then

The air-side horizon is at $\theta_a=90^\circ$, so the image seen from underwater is

Implementation — main.js, runAudit()

const thetaC = Math.asin(1 / N_W);
let maxSeen = 0;
for (let d = 0; d <= 90; d += 0.5) {
    maxSeen = Math.max(
        maxSeen,
        Math.asin(Math.sin(d / DEG) / N_W)
    );
}

The shader uses a backward camera ray, i.e., the water-to-air direction. Therefore refract(I,N,eta) in GLSL takes $\eta=n_w/n_a=N\_W$.

vec3 rr0 = refract(rd, -n, N_W);

refract() returns a zero vector when the value inside the square root is negative. This is the total-internal-reflection branch outside the local cone.

2. Exact Unpolarized Dielectric Fresnel

Let $c_w=\cos\theta_w$ be the cosine of the angle of incidence on the water side; then

At normal incidence,

Implementation — shaders/snells-window.frag, fresnelWA()

float fresnelWA(float cosW) {
    float s = sqrt(max(0.0, 1.0 - cosW * cosW)) * N_W;
    if (s >= 1.0) return 1.0;
    float ca = sqrt(1.0 - s * s);
    float rs = (N_W * cosW - ca) / (N_W * cosW + ca);
    float rp = (N_W * ca - cosW) / (N_W * ca + cosW);
    return clamp(0.5 * (rs * rs + rp * rp), 0.0, 1.0);
}

The sign of the rp amplitude in the code is opposite to the convention shown above. However, since it is squared in the power term, $R_p$ itself is the same.

3. Fresnel Coverage and the $n^2$ Radiance Law

At a lossless refractive interface,

Applying Fresnel transmittance,

The shader approximates sub-pixel coverage with three surface normal samples. It averages the transmittance and refraction direction together, then evaluates the sky once in the weighted direction.

Implementation — shaders/snells-window.frag, main()

float T0 = 1.0 - R0, T1 = 1.0 - R1, T2 = 1.0 - R2;
vec3 rr0 = refract(rd, -n, N_W);
vec3 rr1 = refract(rd, -nA, N_W);
vec3 rr2 = refract(rd, -nB, N_W);
vec3 transmittedDirection = rr0 * T0 + rr1 * T1 + rr2 * T2;
float transmission = (T0 + T1 + T2) / 3.0;
float R = 1.0 - transmission;

if (dot(transmittedDirection, transmittedDirection) > 1e-6) {
    above = skyColor(normalize(transmittedDirection), sd, 1.0)
          * (N_W * N_W);
}
col = above * (1.0 - R) + mirrored * R;

This is coverage anti-aliasing, not full supersampling. The three refracted sky rays are combined into one weighted direction, and the reflection branch still uses only a single coarse-scale normal.

4. Water Absorption, Scattering, and Beam Extinction

The abbreviated water-column model separates absorption $\mathbf a$, total scattering $\mathbf b$, and beam extinction $\mathbf c$.

The 8-step single-scattering approximation is as follows.

Implementation — shaders/snells-window.frag, constants and lookDown(), main()

const vec3 A_W = vec3(0.4524, 0.0654, 0.0083);
const vec3 B_W = vec3(0.0006, 0.0015, 0.0038);
const vec3 C_W = A_W + B_W;

vec3 sunHere = SUN_I * sunT
             * exp(-C_W * (-p.y) / max(-dMean.y, 0.2));
inSc += B_W * (sunHere * 0.13 * gg + vec3(0.10))
       * exp(-C_W * tt) * ds;
col = col * exp(-C_W * min(pathLen, 60.0)) + inSc;

0.13 is an effective phase/source weighting, not a measured volume scattering function. This is therefore a consistent abbreviated single-scattering closure model, not a spectral oceanic radiative transfer equation (RTE) solution.

5. Mean-Centered Surface and the First-Intersection Equation

The procedural octave function returns only positive values. Without mean correction, the previous surface had a mean of roughly $+0.4\,\mathrm{m}$, causing the depth control labeled "depth below mean surface" to be off by that amount. The corrected field is as follows.

Here $\mu_N$ is the spatially measured mean for each LOD. The first root of the equation below gives where the camera ray ${\mathbf r}(s)=\mathbf o+s\mathbf d$ meets the surface.

Implementation — shaders/snells-window.frag, seaHeight3() and traceSurface()

const float SEA_MEAN_3 = 0.4107;

float seaHeight3(vec2 xz, float t) {
    float freq = 0.22, amp = 0.32, choppy = 3.0;
    vec2 uv = vec2(xz.x * 0.75, xz.y);
    float h = 0.0;
    for (int i = 0; i < 3; i++) {
        h += (seaOct((uv + t) * freq, choppy)
            + seaOct((uv - t) * freq, choppy)) * amp;
        uv = OCTM * uv;
        freq *= 1.9;
        amp *= 0.22;
        choppy = mix(choppy, 1.0, 0.2);
    }
    return h - SEA_MEAN_3 + chop(xz, t) * uChop;
}

for (int i = 1; i <= 24; i++) {
    float ti = mix(tLo, tHi, float(i) / 24.0);
    vec3 p = ro + rd * ti;
    float h = p.y - seaHeight3(p.xz, t);
    if (h * hPrev <= 0.0) {
        float a = tPrev, b = ti, ha = hPrev;
        for (int j = 0; j < 8; j++) {
            float m = 0.5 * (a + b);
            vec3 pm = ro + rd * m;
            float hm = pm.y - seaHeight3(pm.xz, t);
            if (hm * ha <= 0.0) {
                b = m;
            } else {
                a = m;
                ha = hm;
            }
        }
        return 0.5 * (a + b);
    }
    tPrev = ti;
    hPrev = h;
}

Unlike the previous silent fallback path, the function returns -1.0 when no sign change is found and no intersection can be located. It no longer fabricates an arbitrary intersection at the top search plane.

6. Caustic Landing Map and Finite-Difference Jacobian

Let $\mathbf s=(s_x,s_z)$ be the surface source coordinate, $\mathbf d(\mathbf s)$ the refracted sun direction, $h(\mathbf s)$ the surface height, and $y_f$ the floor height; then the landing map is as follows.

The geometric-optics raw local gain is approximated as follows.

The value actually displayed on screen is as follows.

Here $N$, $\epsilon_J$, and $G_{max}$ are correction and normalization terms, not optical constants. The raw map therefore returns 1 on a flat surface, but the corrected return value need not equal 1.

Implementation — shaders/snells-window.frag, causticGainFloor()

float surfaceY = seaHeight4(sp, t);
vec3 d = refract(-sd, seaNormal4(sp, 0.05, t), 1.0 / N_W);
vec2 land = sp + d.xz * ((surfaceY - p.y) / (-d.y));

vec2 jx = (Lx - L0) / eps;
vec2 jz = (Lz - L0) / eps;
return min(
    1.0 / max(abs(jx.x * jz.y - jz.x * jx.y), 0.02),
    FLOOR_GMAX
) * FLOOR_NORM;

The previous code used -p.y as the travel height, implicitly placing all source points at $y=0$. It now incorporates the computed wave height. The estimator still uses a single initial guess inverted from the mean plane and does not find or sum all preimages of caustic folds. FLOOR_GMAX and FLOOR_NORM are explicit normalization and correction constants.

7. Focal Length Under the Small-Slope Approximation

For near-normal incidence with $h(x)=a\cos(kx)$, the deflection of small-slope refracted rays produces the following focal length.

Implementation — shaders/snells-window.frag, chop()

return 0.030 * sin(3.5 * (xz.x * 0.92 + xz.y * 0.39) + t * 1.7)
     + 0.014 * sin(5.5 * (xz.x * 0.31 - xz.y * 0.95) - t * 2.1)
     + 0.008 * sin(8.0 * (xz.x * 0.71 + xz.y * 0.70) + t * 2.7);

At the default uChop=1, the focal lengths for these modes are approximately 10.7, 9.3, and 7.8 m. This relationship is an estimate valid under small-slope, near-normal-incidence conditions and is not directly used in the runtime caustic solution.

Problems That Arose While Working with AI

Near the critical angle, the Fresnel reflectance rises sharply from 0.2 to 1.0 within about 2°. A single sample per pixel therefore undersamples this real glint, making it appear as speckling.

This problem produced wrong output roughly three times before it was measured. The first three fixes produced pixel-identical output. The issue became clear only after rendering intermediate values as debug images: the surface normals were perfectly smooth, and $R$ was a pure binary 0/1 value at the fractal boundary. The problem was not the normals; the Fresnel term at that point was effectively behaving like a step function.

An earlier attempt to average only $R$ also failed. That stencil was derived from the pixel footprint, ts * pixAng / max(rd.y, 0.12) ≈ 3.6 cm, while the normals vary at a scale of several meters. The fix that actually worked was widening the stencil to the scale at which the normals exist.

float nEps = clamp(0.35 + ts * 0.03, 0.35, 1.6);   // metres, was clamp(ts*0.012, 0.05, 0.9)

Seafloor caustics and rays suffered the same class of error separately. It was corrected by applying pixel-footprint LOD (lodC, lodA) so that caustic gain fades toward 1 as the footprint grows larger.

Residual speckling remains at the edge of the window. Fixing it properly would require supersampling, which has not yet been applied. The water column still uses only single scattering; multiple scattering is not included.

The current 3-tap interface filter averages transmittance coverage and refraction direction together.

The AI-written code averaged only $R$. When the center tap is in total internal reflection and a neighboring tap transmits, the center refract() returns a zero vector, but after averaging, $(1-R)$ was still nonzero. As a result, interface energy was silently discarded. (This fix was my own personal correction.)

No white balance is applied anywhere. Red genuinely disappears within a few meters. The absorption-only $1/e$ lengths computed with the corrected coefficients are 2.21 / 15.29 / 120.48 m for R/G/B, and the direct beam extinction lengths $1/(a+b)$ are 2.21 / 14.95 / 82.64 m. This is why the result looks blue, and it is also why divers use red filters.

A personal note: working with AI, I find that if you extract and record the values from the reference papers underlying the implementation, it works much faster than I could.

Playing around with GLSL like this is something I always wanted to do since I was young, once I became a programmer. But because I genuinely lacked the talent, I was mostly limited to ray-marching techniques and tweaking the colors and shapes of things others had written. It feels like the world has changed radically.

It's enjoyable to be able to do things I couldn't do before.

References

• E. O. Hulburt, Optics of distilled and natural water, JOSA 35 (1945) 698.

• X. Quan, E. S. Fry, Empirical equation for the index of refraction of seawater, Applied Optics 34 (1995) 3477–3480. DOI: 10.1364/AO.34.003477

• R. M. Pope, E. S. Fry, Absorption spectrum (380–700 nm) of pure water, Applied Optics 36 (1997) 8710–8723. DOI: 10.1364/AO.36.008710

• A. Morel, Optical properties of pure water and pure sea water, in

Optical Aspects of Oceanography (1974), 1–24.

• NASA, Ocean Optics Protocols for Satellite Ocean Color Sensor Validation, Rev. 4, Vol. IV, Table 1.1 (Pope–Fry absorption and Morel scattering values). NASA Ocean Optics Protocols PDF

• R. W. Preisendorfer, Hydrologic Optics, Vol. II: Foundations, §2.6,

  radiance $n^2$ law. NOAA Repository

• T. Nishita et al., Display of the Earth taking into account atmospheric scattering, SIGGRAPH 1993. DOI: 10.1145/166117.166140

• C. D. Mobley, Light and Water: Radiative Transfer in Natural Waters, Academic Press, 1994.