Epidaurus

An acoustics lab for the ancient theatre at Epidaurus, in the browser. It solves the wave equation across a cross-section of the theatre (stage, orchestra and all 55 seat rows) with a finite-difference time-domain (FDTD) solver in WebGPU compute shaders. It then tests the best-known explanation of why the theatre sounds the way it does: that the stepped seat rows act as an acoustic filter. You can watch the wave climb the cavea, compare stepped rows against a smooth slope and a seated audience, and listen to a voice and a murmuring crowd through each.

Open the lab (needs WebGPU) Validation suite

Central claim

Declercq & Dekeyser (2007) modelled Epidaurus with diffraction theory. They argued that the periodic seat rows act as a filter: they suppress low-frequency background noise, such as a murmuring audience or wind, and pass the higher frequencies of the voice. They also argued that waves scattered back from the rows help the listener. Popular accounts usually give the crossover as about 500 Hz.

Others are more cautious. Psarras et al. (2013) measured the theatre and found the response dominated by a dip near 180 Hz, with a lift between 500 and 1500 Hz. Lokki et al. (2013) simulated it with 3D FDTD and beam tracing, matched those measurements, and found that the measurements do not fully support the simulated filtering. Both papers come up again below.

This project rebuilds the mechanism from first principles and measures it three ways:

• the dispersion of waves guided along an infinite staircase of rows;
• the transmission of a plane wave across a flight of rows;
• the full theatre, with a voice in the orchestra and crowd noise in the seats.

What the simulation shows

All results are at Δx = 2 cm, resolved to 2.1 kHz, using Rindel's dimensions. The voice is at the orchestra centre, 1.5 m up, and the crowd noise comes from row 12. Every seat is probed at ear height, 0.8 m above the tread. Levels are relative to the free field at the same distance, so the 2D spreading divides out. The full tables are in docs/results.md.

1. The rows carry a slow surface wave whose band ends at 180 Hz. A single row with Bloch-periodic sides (see Method) carries a guided wave just below the sound line. Near the zone edge it bends away from the sound line and flattens at 180 Hz, below the empty-lattice value c/2d = 207 Hz for the 0.35 × 0.75 m rows. Above that edge, up to the next modes at about 250–290 Hz, nothing can travel along the steps as a guided wave. It is backscattered instead. This is the frequency of the dip Psarras et al. measured in the real theatre.

2. A flight of rows is a band-stop filter at grazing incidence. Twenty-four rows under a plane wave arriving 5° above the slope behave as follows:

• Below 170 Hz, stepped and smooth are within 3 dB of each other.
• At 180 Hz the stepped rows' level drops by 8–14 dB, and stays down to about 700 Hz.
• It recovers by 900 Hz.
• From 1.1 kHz up, the steps give 2–3 dB more than a smooth slope.

3. In the theatre, the steps separate the crowd from the voice between 200 and 700 Hz. Sound from the orchestra reaches each row from below at a steeper angle than noise travelling along the rows from other seats. Averaged over all 55 rows, the steps change the level relative to a smooth slope by:

Path	125 Hz	250 Hz	500 Hz	1 kHz
Voice from the orchestra	+0.6 dB	−6.6 dB	+0.8 dB	−0.8 dB
Crowd noise from row 12	+2.6 dB	−11.8 dB	−8.3 dB	+0.3 dB

Against a smooth slope, the steps improve the voice-to-crowd ratio by about 5 dB in the 250 Hz octave and 9 dB in the 500 Hz octave. At 125 Hz (wind, the lowest male voice) and at 1 kHz they do nothing.

4. An audience matters more than the steps. The seated audience is modelled as a porous layer fitted to Rindel's absorption data. Against bare stepped rows, it takes another 6–25 dB out of the voice at 250 Hz–1 kHz, more for the higher rows and frequencies. That is grazing attenuation over the crowd, and the crowd also fills the corners of the steps. In a full theatre the step filter is largely buried under the audience's own absorption.

5. Backscatter adds late energy rather than clarity. Declercq & Dekeyser credit the backscattered waves with helping the listener. In this model they arrive after the direct sound, and on stepped rows the clarity C50 is about 10 dB lower than on a smooth slope (for example, 29 against 38 dB at row 5). With no reverberation, C50 stays at about 20 dB or more everywhere (lowest: 19.7 dB at row 33). The backscatter is audible as a short tail, not as extra clarity.

Verdict. The mechanism is real, and the frequencies are right where the measurements put them. The steps do work as a filter, but not as a high-pass with a crossover near 500 Hz. They are a band-stop: it starts at the guided-wave band edge, 180 Hz, and is selective for sound that grazes the rows, which is mostly crowd noise. Wind and the lowest frequencies pass. Once an audience is present, its absorption dominates. Lokki et al. reached a similar conclusion: the measured response is explained by the seat rows, but not as a simple filter.

Listening

The Listen panel plays a synthesized voice and a synthesized crowd (murmur plus wind) through the listening responses at the chosen seat, and switches between stepped, smooth, audience and dry. Switching is a crossfade, so differences can be heard within a phrase. A live spectrogram shows the output.

• Voice. A Klatt-style formant synthesizer speaks the Persians' battle cry from Aeschylus, ō paides Hellēnōn, ite. It is generated in code, so there is nothing to license.
• Crowd. The murmur is seven synthetic voices of random syllables, plus gusting low-passed brown noise for wind.
• 2D → 3D correction. A 2D simulation is a line source. Its response falls off as 1/√r and carries a 1/√(jω) tilt. Each response is filtered by √(jω), multiplied by 1/√(2πc·r) with r = c·t, and integrated once to turn the mass-rate source into a pressure source. This is exact for the direct sound and specular reflections in the far field, and approximate for diffraction. It is checked to within 8% of the analytic 3D waveform in the validation suite.
• Bandwidth. At 2 cm the responses stop at 2.1 kHz, so the theatre sounds band-limited. The 1 cm grid doubles that.

The theatre

The cross-section runs along the theatre's axis. Every dimension comes from the literature; assumptions are marked.

Dimension	Value	Source
Orchestra radius	10.15 m (9.77 m earth floor + 0.38 m limestone threshold)	Dörpfeld & Reisch 1896, p. 124
Passage to the first step	2.10 m at the nine central stairways (2.84 m at the wings); one step below the first row	Dörpfeld & Reisch 1896, p. 123; Rindel's first-row diameter of 24.6 m agrees
Lower cavea	34 rows, riser 0.35 m, tread 0.75 m (25.0°)	Rindel 2025, Table I, phase 1
Upper cavea	21 rows, riser 0.43 m (derived), tread 0.75 m	derived so that the 55-row mean riser is Rindel's phase-2 0.38 m
Diazoma	1.82 m	unverified: Wikipedia only; not in Dörpfeld's text, Rindel or Lokki
Proskenion	3.53 m high (optional, off by default)	Dörpfeld & Reisch 1896, p. 127: reconstructed from the colonnade, matching Vitruvius' 10–12 ft
Proskenion depth	2 m	assumed; the skene behind it is not modelled
Alternative rows	riser 0.367 m, tread 0.746 m	Lokki et al. 2013 survey of the lower cavea (toggle Dimensions)
Stone and audience absorption	octave bands 63 Hz–8 kHz	Rindel 2025, Table III
Ear, voice height	0.8 m above the seat; 1.5 m above the orchestra floor	assumed; Lokki et al. use 1.5 m for their orchestra source

On the grid, risers and treads are rounded to whole cells so the rows stay exactly periodic: at 2 cm this gives 0.36 × 0.76 m. The smooth variant is the plane through every tread centre, so a listener's ear is the same height above the surface in both. The audience variant fills each tread to 0.9 m with the porous layer. No images or drawings of the theatre are used; the profile is generated from the numbers above.

Method

Linear acoustics on a staggered grid. Pressure p lives at cell centres and velocity v on cell faces (the acoustic Yee scheme):

$$ \rho,\frac{\partial \mathbf{v}}{\partial t} = -\nabla p, \qquad \frac{\partial p}{\partial t} = -\rho c^2,\nabla\cdot\mathbf{v} $$

$$ v_x^{n+\frac12}\big|_{i+\frac12,j} = v_x^{n-\frac12} - \frac{\Delta t}{\rho,\Delta x}\left(p^n_{i+1,j} - p^n_{i,j}\right), \qquad p^{n+1}_{i,j} = p^n_{i,j} - \frac{\rho c^2 \Delta t}{\Delta x}\left(v_x\big|_{i+\frac12} - v_x\big|_{i-\frac12} + v_y\big|_{j+\frac12} - v_y\big|_{j-\frac12}\right) $$

The scheme is second order in space and time, with c = 343 m/s, ρ = 1.204 kg/m³ and fp32 on the GPU. The Courant number cΔt/Δx = 0.6 is under the 2D limit of 1/√2. One run of the full cavea at 2 cm is 4.8 M cells for 10,290 steps (0.36 s).

Absorbing edges. A split-field PML (Berenger 1994) splits the pressure as p = p_x + p_y, and each part is damped by its own axis' σ, graded cubically over 32 cells for a nominal 10⁻⁵ reflection. The measured reflection is −99 dB.

Stone. Every face between air and stone is a locally reacting impedance Z(s) = R + sM + K/s (Botteldooren 1995 for impedance boundaries in FDTD). The half-cell momentum balance is

$$ \rho_f,\frac{w^{+} - w^{-}}{\Delta t} = \frac{2}{\Delta x}\left(p - R,\bar w - M,\frac{w^{+} - w^{-}}{\Delta t} - K,\bar u\right), \qquad \bar u = u^{-} + \tfrac{\Delta t}{2},\bar w $$

Here w is the velocity into the wall, u its displacement and ρ_f the density of the fluid half-cell. It is discretized with the trapezoidal rule, which maps a passive impedance to a passive discrete one, so the update is stable for any R, M, K ≥ 0. Limestone is fitted to Rindel's octave-band absorption through the closed-form random-incidence (Paris) integral: R = 6.0·10⁴ Pa·s/m and K = 8.5·10⁷ Pa/m, rms error 0.014.

Audience. The audience is an equivalent fluid, the Rayleigh model with flow resistivity σ and porosity φ:

$$ \frac{\rho}{\phi},\frac{\partial \mathbf{v}}{\partial t} + \sigma,\mathbf{v} = -\nabla p, \qquad \frac{\phi}{\rho c^2},\frac{\partial p}{\partial t} = -\nabla\cdot\mathbf{v} $$

A one-parameter fit to the audience absorption wants σ ≈ 180 kPa·s/m². No practical grid can resolve that: the wave dies within 1.5 cm at 500 Hz. With porosity as a second parameter, and the effective drag σφ held to 10 kPa·s/m², the fit is as good: σ = 44 kPa·s/m², φ = 0.23, rms 0.09 against the table. The layer is then resolved. The drag uses Crank–Nicolson, so it is stable however stiff it gets. Faces between air and the porous layer take the larger of the two densities, so the layer's stiffness ρc²/φ cannot outrun the CFL limit. Impedance faces use ρ/φ when the fluid beside the stone is porous. The validation suite and two stability tests check both choices.

Level relative to free field. Each run injects a zero-mass pulse, the derivative of a Gaussian, 40 dB down at the top of the band, into one cell, and records every seat. The transfer function H = P/S comes from regularized spectral division. It is compared with the exact 2D free field at the same distance:

$$ |H_{\text{free}}| = \frac{\Delta x^2}{\Delta t},\frac{\omega,|H_0^{(2)}(kr)|}{4c^2} $$

The result is smoothed in sixth-octaves, and band levels are energy means over each octave. Geometric spreading divides out, so this excess attenuation carries over from 2D to the real theatre, up to the out-of-plane effects listed under Limitations.

Bloch cell. One row, with sheared periodic sides: the right neighbour of the last column is the first column, one riser higher. The complex field obeys f(r + D) = f(r)·e^iφ, with D = (tread, riser). It is solved as coupled real and imaginary parts, with rigid stone and a PML above. A broadband kick is followed by spectral peaks at 24 Bloch phases from 0 to π. The result is compared with the empty-lattice sound line folded into the zone, f = c|k + 2πm/d|/2π, with d = |D| = 0.828 m. This runs on the CPU in a Web Worker.

Strip. A flight of 24 rows preceded by a smooth lead-in, under a plane wave from a tapered line source launched over the lead-in. It is compared with the grating frequencies for a wave at grazing angle θ:

• c/(d(1 + cos θ)), where the backward diffraction order starts to radiate;
• c/(4h), the quarter-wave resonance of a riser-deep step.

Validation

validation.html runs every case on the GPU and compares it with a closed-form answer. npm test runs the same cases on the CPU reference solver. Results on the GPU at Δx = 2 cm (all cases are defined in cells, so the errors scale with the grid):

Case	Reference	Error
Shaders against the CPU reference (every branch, including the sheared wrap)	CpuFdtd in double precision	6·10⁻⁶ and 1·10⁻⁵ of peak
Free-field line source, 30/60/120 cells along an axis	exact 2D Green's function convolved with the pulse	1.2 / 2.4 / 4.9%
Same, along the diagonal	same	0.8 / 0.9 / 1.3%
Peak decay	√(r₁/r₂)	0.5%
2D → 3D correction, 120 cells	3D point source, 1/r	8.2%
Rigid floor	direct + image source	1.7 / 2.6 / 7.9% (out to 200 cells, nearly grazing)
Impedance wall: matched, R + K, R + M	(Z − ρc)/(Z + ρc)	max |ΔR| 0.006 / 0.004 / 0.004
Porous layer: σ = 5 kPa·s/m², and the audience fit	transfer matrix, Zs = −jZc·cot(kd)	max |ΔR| 0.016 / 0.040
Rectangular-groove grating, 0.75 m period	energy of all Floquet orders	balance within 0.5%; first order opens 2.1% from c/d
PML	large-domain reference	−99 dB

The waveform error along the axes grows linearly with distance. That is the numerical dispersion of a second-order scheme, about 0.04% per cell at 8 cells per wavelength; diagonals are four times better at this Courant number. The reflection measurements correct the phase with the scheme's own wavenumber, sin(ωΔt/2)/(cΔt) = sin(k̃Δx/2)/Δx, so that dispersion in the duct is not counted against the boundary.

Controls

Input	Action
Click the field	Move the listener to the nearest row
Shift-click	Move the crowd noise
Wheel / drag / double-click	Zoom / pan / fit
Space, R	Pause, fire the voice again
1 – 3	Stepped rows, smooth slope, with audience
A	Compute every variant, for the voice and for the crowd
L	Listen
↑ / ↓	Listener row

The URL keeps the scene: ?variant=smooth&row=40&noise=12&dx=0.01&cavea=lower&proskenion=on&source=logeion&dims=lokki. The 1 cm grid needs about 0.6 GB of GPU memory for the full cavea, and lower only roughly halves that. Each run takes about four times as many cells and twice as many steps.

Code map

Path	Contents
src/shaders/fdtd.wgsl	Velocity (impedance, porous, PML), pressure and probe kernels
src/physics/fdtd.ts	Grid, CFL, PML profiles, impedance coefficients and CpuFdtd, which mirrors the shader
src/physics/materials.ts	Impedance model, random-incidence absorption, Rayleigh layer, fits to Rindel's table
src/physics/analytic.ts	2D and 3D Green's functions for the pulse, grating frequencies
src/physics/spectra.ts	FFT, deconvolution, smoothing, spectrogram, 2D → 3D correction
src/physics/acoustics.ts	Bessel functions, level re free field, clarity, listening responses, resampling
src/physics/theatre.ts	Sourced dimensions, the profile, the three variants, positions
src/physics/experiments.ts	Theatre runs and their analysis
src/physics/rows.ts	Bloch cell, strip, grating marks
src/physics/cases.ts	Validation cases, shared by the tests and the validation page
src/gpu/	WebGPU solver and field renderer
src/audio/	Formant voice, crowd and wind synthesis; Web Audio auralization
src/analysis/	Result summary, SVG plots
scripts/figures.ts	README figures and docs/results.md from docs/results.json

Tests

npm test runs 40 checks with Node's built-in test runner:

• Solver. CFL, exact discrete energy conservation in a closed box, blow-up above 1/√2, impedance walls that only absorb, a porous layer on stone steps staying stable, the sheared periodic wrap, numerical phase velocity against the dispersion relation, and the rigid material matching the outer wall.
• Materials. The Paris integral (maximum 0.951 near ζ = 1.6, and against numerical integration), stone and audience fits, and penetration depth.
• Signal processing. FFT, deconvolution, band centres, and the 2D → 3D correction of the analytic Green's function.
• Theatre. Dimensions against their sources, snapping, the smooth slope through tread centres, the variants, the proskenion, Bessel functions, and a free-field run at 0 dB re free field.
• Rows and listening. Grating frequencies, the guided wave below the sound line and below c/2d, the strip, the zero-mass source, the listening response as a unit impulse at 10 m, and resampling.
• Validation. Every case on the validation page, on the CPU at Δx = 5 cm (2.5 cm for the grating).

Limitations

• 2D. The cross-section ignores the cavea's curvature in plan (it focuses towards the orchestra), the side walls of each wedge, the stairways and the paths around the side. Each seat is in the plane of the source. Lokki et al.'s 3D results are the reference for those effects. A 3D mode at low frequency is the natural next step.
• Bandwidth. 2.1 kHz at 2 cm (4.3 kHz at 1 cm). The 2 kHz octave is not reported.
• Idealized materials. Uniform limestone fitted to one table: no weathering, cracks, missing or restored blocks, vegetation or soil in the joints. The orchestra floor is the same stone (Rindel lists "beaten earth or stone" with one set of coefficients).
• The audience is an equivalent fluid, not people. It cannot reproduce the sharp rise between 250 and 500 Hz in Rindel's table (fit rms 0.09), and listeners' ears sit inside it.
• Unverified dimensions. The diazoma width (1.82 m) has no primary source here, and the proskenion depth is assumed. Both are away from the lower rows where most of the analysis happens. The definitive survey, von Gerkan & Müller-Wiener's Das Theater von Epidauros (1961), was not available.
• The synthesized voice is intelligible as speech, but it is not a recording of a trained actor's voice.

References

• Declercq, N. F. & Dekeyser, C. S. A. Acoustic diffraction effects at the Hellenistic amphitheater of Epidaurus: Seat rows responsible for the marvelous acoustics. J. Acoust. Soc. Am. 121(4), 2011–2022 (2007).
• Lokki, T., Southern, A., Siltanen, S. & Savioja, L. Acoustics of Epidaurus – studies with room acoustics modelling methods. Acta Acustica united with Acustica 99(1), 40–47 (2013).
• Psarras, S., Hatziantoniou, P., Kountouras, M., Tatlas, N.-A., Mourjopoulos, J. N. & Skarlatos, D. Measurements and analysis of the Epidaurus ancient theatre acoustics. Acta Acustica united with Acustica 99(1), 30–39 (2013).
• Rindel, J. H. Acoustical aspects of the development of Greek theaters in the 4th century B.C.E. J. Acoust. Soc. Am. 157(3), 2042–2066 (2025).
• Dörpfeld, W. & Reisch, E. Das griechische Theater, ch. II.6, "Das Theater in Epidauros", pp. 120–133 (Barth & von Hirst, Athens, 1896).
• Shankland, R. S. Acoustics of Greek theatres. Physics Today 26(10), 30–35 (1973).
• Yee, K. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media. IEEE Trans. Antennas Propag. 14(3), 302–307 (1966).
• Berenger, J.-P. A perfectly matched layer for the absorption of electromagnetic waves. J. Comput. Phys. 114(2), 185–200 (1994).
• Botteldooren, D. Finite-difference time-domain simulation of low-frequency room acoustic problems. J. Acoust. Soc. Am. 98(6), 3302–3308 (1995).
• Kuttruff, H. Room Acoustics, 6th ed., §2.6 (CRC Press, 2016): random-incidence absorption of locally reacting surfaces.
• Klatt, D. H. Software for a cascade/parallel formant synthesizer. J. Acoust. Soc. Am. 67(3), 971–995 (1980).
• Press, W. H. et al. Numerical Recipes, 2nd ed., §6.5 (Cambridge, 1992): Bessel functions J₀ and Y₀.

The design notes are in docs/PLAN.md.

Building

Requirements:

• Node.js 22.18 or newer (includes npm)
• A browser with WebGPU to run it: a recent Chrome or Edge, Safari 26+, or a recent Firefox

Clone and install dependencies:

git clone https://github.com/ylfeid/epidaurus.git
cd epidaurus
npm install

Start the dev server, then open the URL it prints (usually http://localhost:5173):

npm run dev

Run the tests:

npm test

Regenerate the README figures and tables. First, in the lab, run Compute A/B, Compute dispersion and Compute transmission, then save JSON.stringify(epidaurus.exportResults()) from the console as docs/results.json. Then:

npm run figures

Make a production build. This type-checks the code and writes the site to dist/:

npm run build

Serve that build locally to check it (usually at http://localhost:4173):

npm run preview

dist/ is plain static files with relative paths, so it can be uploaded to any static host, including under a subpath.

Deploying

.github/workflows/deploy.yml runs the tests, builds and publishes to GitHub Pages on pushes to main. Enable it under Settings → Pages → Source: GitHub Actions.

License

MIT