Have you ever marvelled at the dancing light patterns at the bottom of a swimming pool? I have, and it sparked a question: could an AI come up with “realistic” physics to simulate it?
I prompted Gemini with the concept, it produced code in Python, which it then translated into the Wolfram Language. The entire simulation was coded by the LLM, essentially “prompted into life.” I was stunned by the depth of physics it managed to implement from just a description:
- “Realistic” Water Surface: The AI didn’t just make simple waves. It implemented a technique called spectral synthesis, using Fast Fourier Transforms (FFTs) to create a complex, natural-looking water surface.
- Accurate Optics: It correctly modelled light rays from a sun source and applied Snell’s Law to calculate how
each ray bends (refracts) as it enters the water.
each ray bends (refracts) as it enters the water.
- The Caustic Patterns: Most impressively, it knew to accumulate the energy of these refracted rays on the pool floor, which is what creates the bright, focused lines that form the “dancing” caustic patterns.
To me, this is a mind-blowing glimpse into the future. We are moving beyond simple code generation to a world where AI can be a partner in scientific exploration and creative coding.
Geometric Optics Simulation of Light Patterns on Swimming Pool Bottoms
This notebook implements a physically-based simulation of caustic patterns -
the beautiful "dancing light" formations seen at the bottom of swimming pools.
The simulation combines:
- Spectral synthesis (FFT-based Tessendorf waves) for realistic water surface
- Geometric optics (Snell's Law) for light refraction
- Ray tracing and accumulation for caustic intensity computation
the beautiful "dancing light" formations seen at the bottom of swimming pools.
The simulation combines:
- Spectral synthesis (FFT-based Tessendorf waves) for realistic water surface
- Geometric optics (Snell's Law) for light refraction
- Ray tracing and accumulation for caustic intensity computation
1. Physical Constants and Helper Functions
1. Physical Constants and Helper Functions
We begin by defining the core physics functions that govern light behaviour.
These functions are vectorized to work efficiently on 2D arrays representing
the water surface.
These functions are vectorized to work efficiently on 2D arrays representing
the water surface.
1.1 Sun Direction Vector
1.1 Sun Direction Vector
The sun direction is specified in spherical coordinates:
- Elevation: angle above the horizon (0° = horizon, 90° = directly overhead)
- Azimuth: compass direction in the horizontal plane
We convert this to a 3D Cartesian unit vector pointing DOWNWARD (from sun to surface).
- Elevation: angle above the horizon (0° = horizon, 90° = directly overhead)
- Azimuth: compass direction in the horizontal plane
We convert this to a 3D Cartesian unit vector pointing DOWNWARD (from sun to surface).
In[]:=
sunDirection[elevDeg_,azDeg_]:=Module[{elev,az,sx,sy,sz},elev=elevDeg*Degree;az=azDeg*Degree;sx=Cos[elev]*Cos[az];sy=Cos[elev]*Sin[az];sz=-Sin[elev];Normalize[{sx,sy,sz}]];
1.2 Surface Gradient Computation
1.2 Surface Gradient Computation
To calculate how light refracts at the water surface, we need the surface normals,
which are derived from the height field gradients. We use central finite differences
with periodic boundary conditions to ensure seamless tiling.
Mathematical formulation:
∂h/∂x ≈ (h[i+1,j] - h[i-1,j]) / (2Δx)
∂h/∂y ≈ (h[i,j+1] - h[i,j-1]) / (2Δy)
which are derived from the height field gradients. We use central finite differences
with periodic boundary conditions to ensure seamless tiling.
Mathematical formulation:
∂h/∂x ≈ (h[i+1,j] - h[i-1,j]) / (2Δx)
∂h/∂y ≈ (h[i,j+1] - h[i,j-1]) / (2Δy)
In[]:=
gradient2D[h_,dx_,dy_]:=Module[{ny,nx,hPadX,hPadY,gradX,gradY},{ny,nx}=Dimensions[h];(*ProperperiodicpaddingusingArrayPad*)hPadX=ArrayPad[h,{{0,0},{1,1}},"Periodic"];(*Padcolumns*)hPadY=ArrayPad[h,{{1,1},{0,0}},"Periodic"];(*Padrows*)(*Centraldifferences*)gradX=Table[(hPadX[[j,i+1]]-hPadX[[j,i-1]])/(2*dx),{j,ny},{i,2,nx+1}];gradY=Table[(hPadY[[j+1,i]]-hPadY[[j-1,i]])/(2*dy),{j,2,ny+1},{i,nx}];{gradX,gradY}];
2. Spectral Wave Surface Generation
2. Spectral Wave Surface Generation
The water surface is generated using Fourier synthesis, following the method
popularized by Tessendorf (2001) for ocean rendering. This approach is preferred
over PDE-based simulations because it:
1. Produces more aesthetically realistic pool-scale ripples
2. Allows precise control over wavelength spectrum
3. Is computationally efficient via FFT
4. Guarantees periodic boundaries (seamless tiling)
popularized by Tessendorf (2001) for ocean rendering. This approach is preferred
over PDE-based simulations because it:
1. Produces more aesthetically realistic pool-scale ripples
2. Allows precise control over wavelength spectrum
3. Is computationally efficient via FFT
4. Guarantees periodic boundaries (seamless tiling)
2.1 Physics Background: Dispersion Relation
2.1 Physics Background: Dispersion Relation
Water waves exhibit DISPERSION: their phase velocity depends on wavelength.
For water of finite depth H, the dispersion relation is:
ω² = g·k·tanh(k·H)
where:
ω = angular frequency (rad/s)
g = gravitational acceleration (9.81 m/s²)
k = wavenumber = 2π/λ (rad/m)
H = water depth (m)
This reduces to:
- Deep water (kH >> 1): ω² ≈ g·k (short wavelengths)
- Shallow water (kH << 1): ω² ≈ g·k²·H (long wavelengths)
For water of finite depth H, the dispersion relation is:
ω² = g·k·tanh(k·H)
where:
ω = angular frequency (rad/s)
g = gravitational acceleration (9.81 m/s²)
k = wavenumber = 2π/λ (rad/m)
H = water depth (m)
This reduces to:
- Deep water (kH >> 1): ω² ≈ g·k (short wavelengths)
- Shallow water (kH << 1): ω² ≈ g·k²·H (long wavelengths)
2.2 Create Spectral Surface Function
2.2 Create Spectral Surface Function
createSurface[Nx_,Ny_,Lx_,Ly_,depth_,rms_,lambda0_,seed_]:=Module[{dx,dy,g,kx,ky,KX,KY,k,omega,k0,sigma,S,h0,h0conj,scale,h0test,x,y,X,Y},dx=Lx/Nx;dy=Ly/Ny;g=9.81;SeedRandom[seed];(*Frequencygrid*)kx=2*Pi*Join[Range[0,Floor[Nx/2]-1],Range[-Ceiling[Nx/2],-1]]/Lx;ky=2*Pi*Join[Range[0,Floor[Ny/2]-1],Range[-Ceiling[Ny/2],-1]]/Ly;KX=ConstantArray[kx,Ny];KY=Transpose[ConstantArray[ky,Nx]];k=Sqrt[KX^2+KY^2];(*Dispersion*)omega=Sqrt[g*k*Tanh[k*depth]];omega=omega/.{Indeterminate->0,ComplexInfinity->0};(*Spectrum*)k0=2*Pi/lambda0;sigma=0.16*k0;S=Exp[-0.5*((k-k0)/sigma)^2]*Exp[-(k/(3*k0))^4];S=S/.{Indeterminate->0};(*Randomamplitudes*)h0=(RandomVariate[NormalDistribution[],{Ny,Nx}]+I*RandomVariate[NormalDistribution[],{Ny,Nx}])*Sqrt[S/2];(*Hermitiansymmetry*)h0conj=Conjugate[h0[[Table[Mod[-i,Ny]+1,{i,0,Ny-1}],Table[Mod[-i,Nx]+1,{i,0,Nx-1}]]]];(*ScaletoRMS*)h0test=Re[InverseFourier[h0+h0conj,FourierParameters->{1,-1}]];scale=rms/StandardDeviation[Flatten[h0test]];h0=h0*scale;h0conj=h0conj*scale;(*Spatialgrid*)x=Range[0,Lx-dx,dx];y=Range[0,Ly-dy,dy];X=ConstantArray[x,Ny];Y=Transpose[ConstantArray[y,Nx]];<|"omega"->omega,"h0"->h0,"h0conj"->h0conj,"X"->X,"Y"->Y,"Lx"->Lx,"Ly"->Ly,"dx"->dx,"dy"->dy,"depth"->depth|>];
2.3 Time Evolution of Surface
2.3 Time Evolution of Surface
The height field at time t is computed via inverse FFT of the time-evolved spectrum:
H(k,t) = h₀(k)·exp(iωt) + h₀*(-k)·exp(-iωt)
h(x,y,t) = IFFT[H(k,t)]
This analytical time evolution in Fourier space ensures the field remains real-valued
and obeys the correct dispersion relation.
H(k,t) = h₀(k)·exp(iωt) + h₀*(-k)·exp(-iωt)
h(x,y,t) = IFFT[H(k,t)]
This analytical time evolution in Fourier space ensures the field remains real-valued
and obeys the correct dispersion relation.
In[]:=
getHeight[surf_,t_]:=Re[InverseFourier[surf["h0"]*Exp[I*surf["omega"]*t]+surf["h0conj"]*Exp[-I*surf["omega"]*t],FourierParameters->{1,-1}]];
3. Caustics Rendering via Ray Tracing
3. Caustics Rendering via Ray Tracing
This is the core rendering function that simulates light refraction through
the water surface and accumulates ray hits on the pool bottom.
Physical process:
1. Sun rays approach the surface from direction I (incident)
2. Rays refract at the air-water interface following Snell's Law
3. Refracted rays propagate through water until hitting the bottom
4. Each ray contributes intensity weighted by Fresnel transmission
5. Ray hits are accumulated into a 2D histogram (the caustic image)
the water surface and accumulates ray hits on the pool bottom.
Physical process:
1. Sun rays approach the surface from direction I (incident)
2. Rays refract at the air-water interface following Snell's Law
3. Refracted rays propagate through water until hitting the bottom
4. Each ray contributes intensity weighted by Fresnel transmission
5. Ray hits are accumulated into a 2D histogram (the caustic image)
3.1 Main Rendering Function
3.1 Main Rendering Function
In[]:=
renderFrame[surf_,t_,res_]:=Module[{h,X,Y,dx,dy,Lx,Ly,depth,gradX,gradY,denom,nx,ny,nz,sunDir,sx,sy,sz,cosi,eta,k,sqrtk,c2,tx,ty,tz,s,bx,by,r0,R,T,w,good,xv,yv,wv,ix,iy,img,mean,maxval},(*Surface*)h=getHeight[surf,t];X=surf["X"];Y=surf["Y"];dx=surf["dx"];dy=surf["dy"];Lx=surf["Lx"];Ly=surf["Ly"];depth=surf["depth"];(*Gradients*){gradX,gradY}=gradient2D[h,dx,dy];(*Normals*)denom=Sqrt[1+gradX^2+gradY^2];nx=-gradX/denom;ny=-gradY/denom;nz=1/denom;(*Refraction*)sunDir=sunDirection[62,25];{sx,sy,sz}=sunDir;cosi=Clip[-(nx*sx+ny*sy+nz*sz),{0,1}];eta=1/1.333;k=1-eta^2*(1-cosi^2);sqrtk=Sqrt[Clip[k,{0,Infinity}]];c2=eta*cosi-sqrtk;tx=eta*sx+c2*nx;ty=eta*sy+c2*ny;tz=eta*sz+c2*nz;(*Intersectbottom*)s=(-depth-h)/(tz-1*^-20);bx=Mod[X+s*tx,Lx];by=Mod[Y+s*ty,Ly];(*Weights*)r0=((1-1.333)/(1+1.333))^2;R=r0+(1-r0)*(1-cosi)^5;T=1-R;w=T*cosi/Clip[nz,{1*^-8,Infinity}];(*Filtervalidrays(downward)*)good=Flatten[Position[Flatten[tz],_?((#<0)&)]];xv=Flatten[bx][[good]];yv=Flatten[by][[good]];wv=Flatten[w][[good]];(*Accumulate*)img=ConstantArray[0.0,{res,res}];ix=Clip[Floor[(xv/Lx)*res]+1,{1,res}];iy=Clip[Floor[(yv/Ly)*res]+1,{1,res}];Do[img[[iy[[i]],ix[[i]]]]+=wv[[i]],{i,Length[wv]}];(*Simpleblur*)img=(img+RotateRight[img,{0,1}]+RotateLeft[img,{0,1}]+RotateRight[img,{1,0}]+RotateLeft[img,{1,0}])/5.0;(*Tonemap*)mean=Mean[Flatten[img]];img=img/(mean+1*^-12);img=Log[1+70*img];maxval=Max[Flatten[img]];img=img/(maxval+1*^-12);img^0.8];
4. Visualization and Color Mapping
4. Visualization and Color Mapping
We use a custom color map that mimics the appearance of water:
- Dark blue in shadows (low intensity)
- Turquoise in mid-tones
- Bright cyan to white in caustic highlights
- Dark blue in shadows (low intensity)
- Turquoise in mid-tones
- Bright cyan to white in caustic highlights
5. Simulation Parameters
5. Simulation Parameters
Physical parameters chosen to match realistic swimming pool conditions:
5.1 Geometry
5.1 Geometry
◼
Domain: 2m * 2m pool section
◼
Grid: 128 * 128 points (sufficient for smooth caustics)
◼
Pool depth: 1.2 m (typical shallow end)
5.2 Wave Characteristics
5.2 Wave Characteristics
◼
Dominant wavelength λ₀ = 0.45 m (45 cm) - realistic pool ripples
◼
Wave height RMS = 8 mm - gentle wind-driven ripples
◼
Spectral width = 16% - natural wave variability
5.3 Lighting
5.3 Lighting
◼
Sun elevation: 62° (mid-afternoon)
◼
Sun azimuth: 25° (sun from slightly east of north)
◼
Refractive indices: n_air = 1.0, n_water = 1.333
6. Run Simulation
6. Run Simulation
Now we execute the simulation in three steps:
1. Create the spectral wave surface
2. Render a test frame to verify correctness
3. Generate full animation sequence
1. Create the spectral wave surface
2. Render a test frame to verify correctness
3. Generate full animation sequence
7. Generate Animation
7. Generate Animation
Generate 60 frames spanning 2 seconds (30 fps) to visualize the time evolution
of caustic patterns as waves propagate across the surface.
of caustic patterns as waves propagate across the surface.
8. Generate Frames and Export
8. Generate Frames and Export
Conclusion
Conclusion
Except from the very first section and these conclusions, the entire post is LLM generated. During a discussion about some ray tracing experiments about black holes, I asked the LLM while driving to generate code for these caustic waves. It produced python code which I then had it translate to Wolfram Language code, which it then annotated.
My “native” programming language is the Wolfram Language. Sometimes it seems to work better for me to have LLMs produce code in Python and then translate it to the Wolfram Language, with which I then keep working. Of course, I could also use the ExternalEvaluator to run external code, but I often prefer having everything in the Wolfram Language.
My “native” programming language is the Wolfram Language. Sometimes it seems to work better for me to have LLMs produce code in Python and then translate it to the Wolfram Language, with which I then keep working. Of course, I could also use the ExternalEvaluator to run external code, but I often prefer having everything in the Wolfram Language.
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
Dancing light at the bottom of a pool (AI-generated caustic waves)
by Marco Thiel
Wolfram Community, STAFF PICKS, December 11, 2025
https://community.wolfram.com/groups/-/m/t/3589458
by Marco Thiel
Wolfram Community, STAFF PICKS, December 11, 2025
https://community.wolfram.com/groups/-/m/t/3589458