#!/usr/bin/env python3 """Chicago from the Michigan dunes: what air has to do, on a globe and on a flat plane. A ray is traced through a temperature profile over Lake Michigan from a camera on a dune to the Chicago skyline, on a globe of radius R and on a flat plane (R -> infinity). The profile is a surface inversion of amplitude DT (air warmer than the water by DT, decaying with e-folding depth d) on top of a standard 6.5 °C/km lapse. What comes out is the "cut": the true height at the city below which nothing can be seen, because the ray that would reach it grazes the water first. Three things are asked of the model. 1. the range of cuts that ordinary lake weather produces, from an airless globe (about 315 m from 59 m at 90.8 km) to full ducting (nothing hidden); 2. the cut on the evening of Levi Miller's chosen time-lapse (Joshua Nowicki, Grand Mere State Park, 30 April 2015), from the recorded weather; 3. the temperature gradient a FLAT plane would need to hide the same amount, and how that compares with the gradient at which air overturns. Nothing here is traced from a photograph; the frame measurements that the page compares against are made separately and stated as such. python3 chicago_refraction.py -> prints the table the page quotes python3 chicago_refraction.py --figure -> docs/long-path-cases/img/fig-chicago-cut.svg """ import math import os import sys HERE = os.path.dirname(os.path.abspath(__file__)) OUT = os.path.join(HERE, "..", "docs", "long-path-cases", "img", "fig-chicago-cut.svg") R_EARTH = 6_371_000.0 # m GAMMA_STD = 0.0065 # K/m, standard lapse AUTOCONVECTIVE = -0.0342 # K/m: the gradient at which k = 0 and air overturns # the 30 April 2015 vantage, from the 3DEP elevation model and the NOAA gauge # (Calumet Harbor 9087044, IGLD85): dune crest 234.5 m, lake 176.7 m, camera 1.5 m H_OBS = 234.5 - 176.7 + 1.5 # 59.3 m above the water D_WILLIS = 90_800.0 # m, Grand Mere crest to Willis Tower CHICAGO_GROUND = 4.0 # m, downtown street level above the lake TOWERS = { # roof heights above street, m (tip in brackets) "Willis": (442.1, 527.3), "Trump": (356.9, 423.2), "Aon": (346.3, 346.3), "Hancock": (343.7, 457.2), } def profile(dT, depth, T0=280.0, gamma=GAMMA_STD): """Temperature (K) as a function of height for a surface inversion of amplitude dT (K) with e-folding depth (m) laid over a lapse gamma (K/m); depth = 0 gives the bare lapse. Pressure is hydrostatic from 1013.25 hPa at the water.""" def T(h): inv = dT * (1.0 - math.exp(-h / depth)) if depth > 0 else 0.0 return T0 + inv - gamma * h return T GRID = 0.25 # m: the refractive-index table's spacing TOP = 800.0 def n_of_h(T): """Refractive index tabulated on a fine grid (hydrostatic pressure integrated on the same grid), with linear interpolation and a centred finite-difference gradient.""" kmg = 9.80665 * 0.0289644 / 8.314462 hs = [i * GRID for i in range(int(TOP / GRID) + 1)] p, ns = 1013.25, [] for h in hs: if h > 0: p *= math.exp(-kmg * GRID / T(h - GRID / 2)) ns.append(1.0 + 77.6e-6 * p / T(h)) ps = [1013.25] def n(h): x = min(max(h, 0.0), TOP - GRID) / GRID i = int(x); f = x - i return ns[i] * (1 - f) + ns[i + 1] * f def dn(h): i = min(max(int(h / GRID), 1), len(ns) - 2) return (ns[i + 1] - ns[i - 1]) / (2 * GRID) n.dn = dn n.T = T return n def trace(n, h0, elev0, dist, R=R_EARTH, ds=25.0): """Trace a ray from height h0 at elevation elev0 (rad) for ground distance dist. Returns the height at dist, or None if the ray meets the water first. dε/ds = cos ε / (R + h) [the ground curves away] + (dn/dh)/n · cos ε [the ray bends]""" h, e, s = h0, elev0, 0.0 while s < dist: step = min(ds, dist - s) c = math.cos(e) curv = (c / (R + h) if R else 0.0) + n.dn(h) / n(h) * c h += math.sin(e) * step e += curv * step s += step if h <= 0.0: return None return h def cut_height(n, h0=H_OBS, dist=D_WILLIS, R=R_EARTH): """True height at the target reached by the ray that just grazes the water: everything below it is hidden. Bisection on the launch elevation.""" lo, hi = -0.05, 0.05 for _ in range(50): mid = 0.5 * (lo + hi) if trace(n, h0, mid, dist, R) is None: lo = mid else: hi = mid h = trace(n, h0, hi, dist, R) return h if h is not None else 0.0 def k_of(n, h=2.0): """Refraction coefficient at height h: the ray's curvature as a fraction of the Earth's.""" return -R_EARTH * n.dn(h) / n(h) def airless_cut(h0=H_OBS, dist=D_WILLIS, R=R_EARTH): dh = math.sqrt(2 * R * h0) return (dist - dh) ** 2 / (2 * R) if dist > dh else 0.0 def flat_plane_gradient(target_cut, h0=H_OBS, dist=D_WILLIS): """The uniform temperature gradient (K/m) a flat plane needs for the same cut: rays must curve upward, so n must increase with height. Bisection on gamma.""" def cut_for(gamma): return cut_height(n_of_h(profile(0.0, 0.0, gamma=gamma)), h0, dist, R=None) lo, hi = 0.0, 0.5 # gamma in K/m, i.e. temperature FALLING with height at gamma for _ in range(40): mid = 0.5 * (lo + hi) if cut_for(mid) < target_cut: lo = mid else: hi = mid return -0.5 * (lo + hi) # sign convention: dT/dh FRAMES = os.path.join(HERE, "..", "docs", "long-path-cases", "data", "nowicki-2015-04-30-frames.json") LAKE_HEIGHTS = { # true height above the lake of each measured feature, m "willis_tip": 527.3 + CHICAGO_GROUND, "willis_roof": 442.1 + CHICAGO_GROUND, "willis_tier": 360.0 + CHICAGO_GROUND, "hancock_roof": 343.7 + CHICAGO_GROUND, "trump_roof": 356.9 + CHICAGO_GROUND, "aon_roof": 346.3 + CHICAGO_GROUND, } USABLE = range(10, 23) # seconds 10-22: caption gone, sky still bright enough for the masts AZ_SPAN_DEG = 1.296 # Willis to Hancock, from the Grand Mere crest PX_SPAN = 906 # the same two towers in the frame def fit_frames(path=FRAMES): """Per frame, least squares of pixels-above-horizon against true height: the x-intercept is the cut (the true height at the water horizon), the slope the vertical scale. Returns [(second, cut_m, m_per_px)] over the usable frames.""" import json out = [] for r in json.load(open(path))["frames"]: if r["second"] not in USABLE: continue pts = [(LAKE_HEIGHTS[k], r[k]) for k in LAKE_HEIGHTS if r.get(k) is not None] n = len(pts); sx = sum(p[0] for p in pts); sy = sum(p[1] for p in pts) sxx = sum(p[0] ** 2 for p in pts); sxy = sum(p[0] * p[1] for p in pts) slope = (n * sxy - sx * sy) / (n * sxx - sx ** 2) icpt = (sy - slope * sx) / n # one-sigma error on the cut from the residuals (ordinary least squares) res = sum((y - slope * x - icpt) ** 2 for x, y in pts) / (n - 2) dxx = sxx - sx ** 2 / n var_m, var_b, cov = res / dxx, res * sxx / (n * dxx), -res * (sx / n) / dxx dc_dm, dc_db = icpt / slope ** 2, -1.0 / slope err = math.sqrt(max(dc_dm ** 2 * var_m + dc_db ** 2 * var_b + 2 * dc_dm * dc_db * cov, 0.0)) out.append((r["second"], -icpt / slope, 1.0 / slope, err)) return out def horizontal_scale(dist=D_WILLIS): """Metres per pixel across the frame, from the towers' azimuths: refraction bends rays in the vertical, so this scale is the one it cannot touch.""" return dist * math.radians(AZ_SPAN_DEG) / PX_SPAN def sweep(): rows = [] for dT, depth in ((0.0, 0.0), (1.0, 50.0), (2.0, 50.0), (3.0, 50.0), (4.0, 50.0), (6.0, 50.0), (2.0, 20.0), (4.0, 20.0), (6.0, 20.0), (8.0, 20.0)): n = n_of_h(profile(dT, depth)) rows.append((dT, depth, k_of(n), cut_height(n))) return rows def main(): print(f"observer {H_OBS:.1f} m above the lake, {D_WILLIS/1000:.1f} km from Willis Tower") print(f"airless globe hides {airless_cut():.0f} m; standard air (k about 0.17) hides", end=" ") n = n_of_h(profile(0.0, 0.0)) print(f"{cut_height(n):.0f} m (k = {k_of(n):.3f} at 2 m)") print("\ninversion sweep (air warmer than water by dT, e-folding depth d):") print(" dT d k(2 m) cut") for dT, d, k, c in sweep(): print(f" {dT:4.1f} {d:4.0f} {k:5.2f} {c:5.0f} m") fits = fit_frames() cuts = [c for _, c, _, _ in fits]; scales = [m for _, _, m, _ in fits] print(f"\nthe frames, seconds 10-22: cut {min(cuts):.0f}-{max(cuts):.0f} m (median {sorted(cuts)[len(cuts)//2]:.0f}), " f"vertical scale {min(scales):.2f}-{max(scales):.2f} m/px against {horizontal_scale():.2f} m/px horizontally") for sec, c, m, e in fits: print(f" {sec:2d} s cut {c:5.0f} +/- {e:4.0f} m {m:.2f} m/px") print("\nthe flat plane's price for the frames' median cut and for standard air's:") for c in (160.0, 237.0): g = flat_plane_gradient(c) print(f" cut {c:.0f} m -> dT/dh = {g*1000:+.0f} °C/km ({g/AUTOCONVECTIVE:.1f} x the overturning gradient of {AUTOCONVECTIVE*1000:.1f} °C/km)") def figure(out=OUT): """Cut against inversion strength, with the airless and standard-air lines and the band the frames measured.""" W, H = 760, 420 L, R_, T, B = 70, 730, 40, 350 dts = [i * 0.25 for i in range(0, 25)] cuts = [cut_height(n_of_h(profile(dt, 50.0))) for dt in dts] cuts200 = [cut_height(n_of_h(profile(dt, 200.0))) for dt in dts] fits = fit_frames(); lo, hi = min(c for _, c, _, _ in fits), max(c for _, c, _, _ in fits) def X(dt): return L + dt / 6.0 * (R_ - L) def Y(c): return B - c / 340.0 * (B - T) o = [f'', f'', f'How much of Chicago the air hides, from the Grand Mere crest (59 m, 90.8 km)', f'', f'the evening of 30 April 2015, from the shore stations:', f'air 1–3 °C warmer than the open lake', f'', f'measured in the frames the same evening: {lo:.0f}–{hi:.0f} m hidden'] for c, lab, col in ((airless_cut(), f"airless globe: {airless_cut():.0f} m", "#a09a8c"), (cuts[0], f"standard air (no inversion), k ≈ 0.18: {cuts[0]:.0f} m", "#6E675B")): o.append(f'') o.append(f'{lab}') pts = " ".join(f"{X(dt):.1f},{Y(c):.1f}" for dt, c in zip(dts, cuts)) o.append(f'') pts2 = " ".join(f"{X(dt):.1f},{Y(c):.1f}" for dt, c in zip(dts, cuts200)) o.append(f'') o.append(f'the same warmth spread through 200 m: most of the city stays hidden') o.append(f'solid: the warmth held in the bottom 50 m. The gradient is what bends the light,') o.append(f'and at six degrees in fifty metres it reaches k ≈ 1 — the ray follows the lake') for c in range(0, 341, 100): o.append(f'') o.append(f'{c}') for dt in range(0, 7): o.append(f'') o.append(f'{dt}') o.append(f'') o.append(f'air warmer than the lake surface by (°C) — solid: within a 50 m layer; dashed: through 200 m') o.append(f'true height hidden at the city (m)') o.append(f'on a flat plane the same {lo:.0f}–{hi:.0f} m needs air cooling upward') o.append(f'at 125–150 °C per km, four times past overturning') o.append("") os.makedirs(os.path.dirname(out), exist_ok=True) open(out, "w").write("\n".join(o)) print("wrote", os.path.normpath(out)) if __name__ == "__main__": if "--figure" in sys.argv: figure() else: main()