Fortran Physics Demos

Small, self-contained Fortran programs exploring the math behind everyday physical scenarios. Each program is a single .f90 file with no dependencies beyond a Fortran 2008 compiler.

• milk_drop_ripple.f90 — ripple eigenmodes and rebound-drip mixing when a drop of milk lands in a glass of water.
• speed_of_sound.f90 — how the speed of sound differs across dry air, humid air, fresh water, and salt water, and how each varies with temperature.

General build pattern (per-program compile commands are below):

gfortran -O2 -std=f2008 -Wall -o <program_name> <program_name>.f90
./<program_name>

---

Milk Drop in a Glass of Water

A single drop of milk falls into the center of a 12 cm diameter glass of water. This program models two things about that event:

1. Ripple wave frequencies — the standing-wave (eigenmode) patterns of the ripples in the circular glass, and the frequency of each mode.
2. Milk fraction of the "rebound drip" — when a drop hits a pool it digs a crater, the crater collapses, and a thin jet (a "Worthington jet") shoots back up, sometimes pinching off a secondary droplet (the "rebound drip" / crown-splash droplet familiar from Harold Edgerton's high-speed milk-drop photographs). The program estimates what fraction of that rebound droplet is milk versus entrained water.

The physics — what's real and what's a toy model

The ripple frequencies rest on real physics. Surface waves in water follow the capillary-gravity dispersion relation:

omega(k)^2 = ( g*k + (sigma/rho)*k^3 ) * tanh(k*h)

where k is the wavenumber, g gravity, sigma surface tension, rho water density, and h the water depth. The glass wall imposes a rigid (no-radial-flow) boundary condition, which restricts the allowed wavenumbers to k_mn = root_n(J_m') / a, i.e. the zeros of the derivative of the Bessel function J_m, where a is the glass radius — the same math used for vibrating drumheads and Chladni plate patterns. The program finds these Bessel roots numerically and reports frequency, wavenumber, and wavelength for azimuthal modes m = 0..4 and radial modes n = 1..3.

The milk-fraction number is a simplified, illustrative model, not a measured fact. How much of the rebound droplet is actually the original milk (versus water pulled in from the crater) depends on turbulent mixing during crater collapse, contact-line effects, and viscosity contrasts between milk and water — this is genuinely an open, empirical fluid- mechanics question that in practice requires two-phase (VOF / level-set) Navier–Stokes simulation with species tracking to answer precisely. There is no simple textbook formula for it.

What the program actually does: it estimates crater depth from the impact Weber number using a rough depth ~ d * We^(1/4) scaling, assumes the milk drop mixes uniformly into that estimated crater volume, and reports the resulting concentration as a stand-in for the rebound droplet's milk fraction. The proportionality constants (c_crater_depth, c_crater_radius in the source) are illustrative guesses, not values calibrated against experiment — treat that output as a toy estimate to play with, not a scientific result.

Compiling and running

gfortran -O2 -std=f2008 -Wall -o milk_drop_ripple milk_drop_ripple.f90
./milk_drop_ripple

This prints:

• Impact velocity, Weber/Reynolds/Froude numbers for the falling drop.
• Estimated crater depth/radius and the resulting milk-fraction estimate.
• A table of ripple eigenmodes (m, n, wavenumber, wavelength, frequency) for the 12 cm glass.

Changing the scenario

All scenario parameters live as parameter declarations near the top of program milk_drop_ripple in the source file:

• glass_diameter — diameter of the glass (default 0.12 m / 12 cm)
• drop_diameter — size of the milk drop (default 4 mm, a typical eyedropper-sized drop)
• drop_height — release height above the water surface (default 2 cm)
• water_depth, fluid densities, viscosity, and surface tension are set in the physics_constants module at the top of the file.

Edit these and recompile to explore other drop sizes, glass sizes, or release heights.

---

Speed of Sound: Air vs Fresh Water vs Salt Water

How fast sound travels through dry air, humid air, fresh water, and seawater, and how each of those speeds changes with temperature. Unlike the milk-drop rebound-mixing model above, every formula used here is an established result from acoustics/thermodynamics references, not a toy estimate.

The physics

Air — ideal-gas result. Sound in a gas is an adiabatic pressure disturbance:

c_air = sqrt( gamma * R_specific * T )

with T the absolute temperature (K), gamma = Cp/Cv = 1.400 for dry air, and R_specific = R_universal / M_air = 287.05 J/(kg K). Notably this does not depend on pressure — raising pressure increases both a gas's stiffness and its density by the same factor, so they cancel. It does depend on molar mass, which is why humid air (avg. molar mass lowered by lighter H2O molecules, M=18 g/mol vs dry air's M=29 g/mol) carries sound very slightly faster than dry air at the same temperature. The program includes this humidity correction via the Magnus/Tetens saturation-vapor-pressure formula.

Fresh water — Marczak (1997) polynomial. In a liquid, c = sqrt(K_s / rho) where K_s is the adiabatic bulk modulus, but K_s(T) and rho(T) for water don't have simple closed forms, so pure-water sound speed is fit empirically instead:

c(T) = 1402.385 + 5.038813 T - 5.799136e-2 T^2 + 3.287156e-4 T^3
       - 1.398845e-6 T^4 + 2.787860e-9 T^5        (T in deg C)

Valid 0–95°C at 1 atm.

Salt water — Mackenzie (1981) equation. The standard oceanographic sound-speed formula, valid T 0–30°C, salinity S 25–40 ppt (parts per thousand), depth D 0–8000 m:

c(T,S,D) = 1448.96 + 4.591 T - 5.304e-2 T^2 + 2.374e-4 T^3
           + 1.340 (S - 35) + 1.630e-2 D + 1.675e-7 D^2
           - 1.025e-2 T (S - 35) - 7.139e-13 T D^3

The program uses standard open-ocean surface salinity (S = 35 ppt) at D = 0 (surface), so the comparison isolates the temperature and salinity effects rather than depth/pressure effects.

Sample results (0–40°C)

 T(C)   c_air,dry   c_air,100%RH   c_freshwater   c_seawater
  0.0     331.3         331.7         1402.4        1449.0
 10.0     337.3         338.1         1447.3        1489.8
 20.0     343.2         344.7         1482.4        1521.5
 30.0     349.0         351.8         1509.2        1545.4
 40.0     354.8         359.7         1528.9        1562.9

A few takeaways:

• Water carries sound roughly 4.3–4.4x faster than air — water's much higher bulk modulus outweighs its higher density.
• Salinity adds about 1.3 m/s per ppt: at 20°C, seawater (35 ppt) is about 39 m/s faster than fresh water.
• Air's speed rises almost linearly with temperature (~0.6 m/s per °C) because sqrt(T) is nearly linear over a modest temperature span. Water's isn't as clean — pure water's sound speed actually peaks around 74°C and decreases above that, a quirk tied to water's density anomaly, visible in the higher-order terms of the Marczak polynomial (not shown in the 0–40°C range printed by the program).

Compiling and running

gfortran -O2 -std=f2008 -Wall -o speed_of_sound speed_of_sound.f90
./speed_of_sound

This prints a table of all four sound speeds from 0–40°C in 5°C steps, plus ratios and the salinity effect at 20°C.

Changing the scenario

Edit the parameter declarations in program speed_of_sound_demo:

• salinity_ocean — seawater salinity in ppt (default 35, typical open ocean)
• depth_surface — depth in meters (default 0)
• Tmin, Tmax, dT — temperature sweep range and step size

The air_speed, freshwater_speed, and seawater_speed functions in the sound_speed module can also be called directly with other temperatures, humidity levels, salinities, or depths within each formula's valid range (noted above).