#!/usr/bin/env python3 """What the air had to be doing, case by case: the refraction coefficient each Chicago and Pontchartrain picture implies on a globe and on a flat plane, and the temperature gradient each coefficient is. The solver is the one at https://funwithscience.net/refraction-solver/ : for a camera at height h, a target at distance D and a hidden height H, globe (D - sqrt(2 R h / (1 - k)))^2 (1 - k) / (2 R) = H -> k flat the same hiding must come from rays curving AWAY from the surface, radius r with (D - sqrt(2 r h))^2 / (2 r) = H, k' = -R / r, which is the degeneracy k' = k - 1 exactly. dT/dh = k / (503 P / T^2) - 0.0342 (K per metre; Bislin's form) Regime bands on the axis are the solver's descriptors: inferior mirage k < 0, normal marine 0.10-0.20, minor looming 0.20-0.35, major 0.35-0.60, duct above. python3 long_path_k_bands.py -> prints the table python3 long_path_k_bands.py --figure -> docs/long-path-cases/img/fig-air-gradient.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") R = 6_371_000.0 AUTO = -0.0342 def hidden_globe(h, D, k): Re = R / (1.0 - k) dh = math.sqrt(2 * Re * h) return (D - dh) ** 2 / (2 * Re) if D > dh else 0.0 def k_globe(h, D, H): """Bisection on k for the hidden height H (H = 0 -> the largest k that just hides nothing).""" lo, hi = -3.0, 0.999 for _ in range(80): mid = 0.5 * (lo + hi) if hidden_globe(h, D, mid) > H: lo = mid # too much hidden: need more bending else: hi = mid return 0.5 * (lo + hi) def k_flat(h, D, H): return k_globe(h, D, H) - 1.0 def gradient(k, T=280.0, P=1013.25): """dT/dh in °C/km for a coefficient k at surface temperature T (K).""" return (k / (503.0 * P / T ** 2) + AUTO) * 1000.0 # (label, (camera height low, high) m, distance m, (hidden low, high) m, surface T K, note) CASES = { "Chicago": [ ("Nowicki time-lapse, 30 Apr 2015: measured", (59.3, 59.3), 90_800.0, (138.0, 203.0), 280.0, "cut 138-203 m from the frames"), ("p. 111 still, what the frame shows: a bound, k at least this", (2.0, 2.0), 62_355.0, (233.0, 263.0), 285.0, "at most 750-860 ft hidden from the beach"), ("p. 111 overlay's own number, 1,014 ft hidden", (2.0, 2.0), 62_355.0, (309.0, 309.0), 285.0, "8 in/mi^2 x 39^2 read as a hidden height"), ], "Pontchartrain": [ ("p. 112 photo A, 2017: the far towers gone; camera 2-5 m", (2.0, 5.0), 24_270.0, (23.0, 23.0), 303.0, "23 m towers fully hidden at 24 km"), ("p. 112 photo B and the linked clip: tops reach the shore, caps gone; camera 3-15 m", (3.0, 15.0), 24_270.0, (3.0, 23.0), 303.0, "constrains almost nothing"), ], } def solve_all(): rows = [] for place, cases in CASES.items(): for lab, (hlo, hhi), D, (Hlo, Hhi), T, note in cases: ks = [k_globe(h, D, H) for h in (hlo, hhi) for H in (Hlo, Hhi)] kg = (min(ks), max(ks)) kf = (kg[0] - 1.0, kg[1] - 1.0) rows.append((place, lab, (hlo, hhi), D, Hlo, Hhi, T, note, kg, kf)) return rows def main(): for place, lab, (hlo, hhi), D, Hlo, Hhi, T, note, kg, kf in solve_all(): g = sorted(kg); f = sorted(kf) print(f"{place}: {lab}") print(f" camera {hlo:.1f}-{hhi:.1f} m, {D/1000:.1f} km, hidden {Hlo:.0f}-{Hhi:.0f} m [{note}]") print(f" globe k {g[0]:+.2f} .. {g[1]:+.2f} dT/dh {gradient(g[0], T):+.0f} .. {gradient(g[1], T):+.0f} °C/km") print(f" flat k {f[0]:+.2f} .. {f[1]:+.2f} dT/dh {gradient(f[0], T):+.0f} .. {gradient(f[1], T):+.0f} °C/km" f" ({gradient(f[0], T)/(AUTO*1000):.1f}-{gradient(f[1], T)/(AUTO*1000):.1f} x overturning)") GLOBE, FLAT = "#e0b458", "#8b1e1e" # light: what a globe needs; dark: what a flat plane needs OUT_FIG = os.path.join(OUT, "fig-air-gradient.svg") def figure(out=OUT_FIG): """One picture: each case's implied temperature gradient on a globe (light) and on a flat plane (dark), on an axis of degrees C per kilometre, against the gradient at which air overturns and the band ordinary lake weather occupies.""" rows = [r for r in solve_all() if "photo B" not in r[1]] W, H = 900, 96 + len(rows) * 58 + 70 L, R_ = 60, 860 G0, G1 = -240.0, 60.0 def X(g): return L + (g - G0) / (G1 - G0) * (R_ - L) top = 96; rowh = 58; y0 = top; y1 = top + len(rows) * rowh short = {"Nowicki time-lapse, 30 Apr 2015: measured": "Chicago time-lapse, 30 April 2015 (measured)", "p. 111 still, what the frame shows: a bound, k at least this": "Chicago still, p. 111 (what the frame shows)", "p. 111 overlay's own number, 1,014 ft hidden": "Chicago still, p. 111 (the overlay’s own 1,014 ft)", "p. 112 photo A, 2017: the far towers gone; camera 2-5 m": "Pontchartrain photo A, p. 112 (camera 2–5 m)"} o = [f'', f'', f'How fast the air must cool with height for each picture to exist', f'on a globe', f'on a flat plane, for the same picture', # zones f'', f'', f'air cannot do this', f'cooling faster than −34 °C per km, the air is denser above than below and overturns', f'lake weather', f'standard air to a spring inversion over cold water', f''] y = top for _, lab, hh, D, Hlo, Hhi, T, note, kg, kf in rows: cy = y + 34 o.append(f'{short[lab]}') for (lo, hi), col, tag in ((sorted(kg), GLOBE, "globe"), (sorted(kf), FLAT, "flat")): glo, ghi = gradient(lo, T), gradient(hi, T) if ghi - glo < 4: o.append(f'') else: o.append(f'') txt = f"{glo:+.0f} to {ghi:+.0f}" if ghi - glo >= 4 else f"{0.5*(glo+ghi):+.0f}" o.append(f'{txt}') y += rowh for g in range(-240, 61, 40): o.append(f'') o.append(f'{g:+d}') o.append(f'') o.append(f'temperature change with height, °C per kilometre (− cooling upward, + an inversion)') o.append(f'on the measured evening the flat plane’s bar means air at −50 to −60 °C at the Willis roof above a 6 °C lakefront') o.append("") open(out, "w").write("\n".join(o)) print("wrote", os.path.normpath(out)) def airless(h, D): return hidden_globe(h, D, 0.0) if __name__ == "__main__": if "--figure" in sys.argv: figure() else: main()