From bd107ee747b6f78bd3ddb5f60fb062daf9d796d5 Mon Sep 17 00:00:00 2001 From: Graham McIntire Date: Thu, 23 Apr 2026 16:18:19 -0500 Subject: [PATCH] feat(eme): add /eme Earth-Moon-Earth calculator MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit New /eme page — antenna aiming + link budget for moonbounce. Given a station (callsign / grid / lat,lon), TX power, antenna gain, detection bandwidth, and system Tsys, the page shows in real time: - Moon azimuth / elevation from the observer (ticks every 10 s) - Slant range to the Moon - EIRP breakdown: TX(dBm) + gain(dBi) = EIRP(dBm) - Path-loss decomposition (Earth → Moon → Earth): free-space spreading (×2) − moon reflection gain 10·log(4π σ / λ²) where σ = ρ·π·R² = round-trip path loss - Band-dependent lunar albedo ρ (VHF ≈ 0.065, microwave ≈ 0.08, mm-wave ≈ 0.02) surfaced as a named line item with the current ρ shown to three decimals - Received power, thermal noise floor, and SNR with a badge that classifies the margin (strong / marginal / below noise) - Self-echo Doppler shift with its radial-velocity driver New modules: - `Microwaveprop.Moon` — low-precision lunar ephemeris (Astronomical Almanac Sec. D.4 truncated series). Provides `julian_day/1`, `geocentric/1`, `observer_position/3`, and `radial_velocity_km_s/3`. Accurate to ~0.3° position / ~200 km distance — well inside any amateur EME-dish beamwidth. - `Microwaveprop.Propagation.Eme` — path-loss helpers: `moon_albedo/1` (band lookup), `moon_rcs_m2/1`, `fspl_round_trip_db/2`, `moon_reflection_gain_db/2`, `path_loss_db/3`, `received_power_dbm/1`, `noise_floor_dbm/2`, `snr_db/1`, `doppler_shift_hz/2`. Wired into the main nav ("EME" button) and the router between /path and /rover. Tests: 3 properties + 41 unit tests cover JD conversion, moon position envelope, band-dependent albedo lookup, decomposition identity, linearity/monotonicity of path loss and Doppler, and a LiveView integration test against the rendered DOM. Full suite: 2,460 tests + 170 properties; credo strict clean. --- lib/microwaveprop/moon.ex | 193 ++++++ lib/microwaveprop/propagation/eme.ex | 197 ++++++ lib/microwaveprop_web/components/layouts.ex | 1 + lib/microwaveprop_web/live/eme_live.ex | 603 ++++++++++++++++++ lib/microwaveprop_web/router.ex | 1 + test/microwaveprop/moon_test.exs | 195 ++++++ test/microwaveprop/propagation/eme_test.exs | 291 +++++++++ test/microwaveprop_web/live/eme_live_test.exs | 58 ++ 8 files changed, 1539 insertions(+) create mode 100644 lib/microwaveprop/moon.ex create mode 100644 lib/microwaveprop/propagation/eme.ex create mode 100644 lib/microwaveprop_web/live/eme_live.ex create mode 100644 test/microwaveprop/moon_test.exs create mode 100644 test/microwaveprop/propagation/eme_test.exs create mode 100644 test/microwaveprop_web/live/eme_live_test.exs diff --git a/lib/microwaveprop/moon.ex b/lib/microwaveprop/moon.ex new file mode 100644 index 00000000..9591c76b --- /dev/null +++ b/lib/microwaveprop/moon.ex @@ -0,0 +1,193 @@ +defmodule Microwaveprop.Moon do + @moduledoc """ + Low-precision lunar ephemeris. + + Uses the Astronomical Almanac "Low-Precision Formulas for the Moon" + (Section D.4) — 6 longitude terms, 4 latitude terms, and 4 parallax + terms. Accurate to ~0.3° on position and ~200 km on distance, which + is well inside the beamwidth of any amateur EME dish at every + supported microwave band (1° beam at 10 GHz for a 1 m dish; wider + below that). + + Exposes two coordinates: + + - `geocentric/1` — the Moon's right ascension, declination, and + distance from Earth's centre at a UTC instant. + - `observer_position/3` — topocentric azimuth / elevation / + slant range for an observer at `(lat, lon)` degrees. + + Both are pure functions of UTC plus observer coordinates; no side + effects, no DB access. Callers drive them on a timer to refresh + antenna aim points. + """ + + @earth_radius_km 6378.14 + + @doc """ + Julian Day for a Gregorian `DateTime`. Caller is responsible for + passing UT; TAI/UT1 offset (< 1 s) is ignored — the low-precision + lunar series doesn't resolve it anyway. + """ + @spec julian_day(DateTime.t()) :: float() + def julian_day(%DateTime{} = dt) do + day_frac = (dt.hour + dt.minute / 60 + dt.second / 3600) / 24 + + {y, m} = + if dt.month <= 2 do + {dt.year - 1, dt.month + 12} + else + {dt.year, dt.month} + end + + a = div(y, 100) + b = 2 - a + div(a, 4) + + trunc(365.25 * (y + 4716)) + trunc(30.6001 * (m + 1)) + dt.day + day_frac + b - 1524.5 + end + + @doc """ + Geocentric ecliptic → equatorial position of the Moon at a UTC + instant. Returns `%{ra, dec, distance_km}` where `ra` and `dec` + are in degrees. + """ + @spec geocentric(DateTime.t()) :: %{ra: float(), dec: float(), distance_km: float()} + def geocentric(%DateTime{} = dt) do + jd = julian_day(dt) + t = (jd - 2_451_545.0) / 36_525.0 + + # Mean-element arguments (degrees, wrapped before trig conversion). + l = 218.32 + 481_267.881 * t + m_prime = 134.9 + 477_198.85 * t + m = 357.53 + 35_999.05 * t + d = 297.85 + 445_267.11 * t + f = 93.27 + 483_202.03 * t + + m_prime_r = deg_to_rad(m_prime) + m_r = deg_to_rad(m) + d_r = deg_to_rad(d) + f_r = deg_to_rad(f) + + # Ecliptic longitude λ (degrees). + lambda = + l + 6.29 * :math.sin(m_prime_r) - + 1.27 * :math.sin(m_prime_r - 2 * d_r) + + 0.66 * :math.sin(2 * d_r) + + 0.21 * :math.sin(2 * m_prime_r) - + 0.19 * :math.sin(m_r) - + 0.11 * :math.sin(2 * f_r) + + # Ecliptic latitude β (degrees). + beta = + 5.13 * :math.sin(f_r) + + 0.28 * :math.sin(m_prime_r + f_r) - + 0.28 * :math.sin(m_prime_r - f_r) - + 0.17 * :math.sin(2 * d_r - f_r) + + # Horizontal parallax π (degrees). + parallax = + 0.9508 + + 0.0518 * :math.cos(m_prime_r) + + 0.0095 * :math.cos(m_prime_r - 2 * d_r) + + 0.0078 * :math.cos(2 * d_r) + + 0.0028 * :math.cos(2 * m_prime_r) + + distance_km = @earth_radius_km / :math.sin(deg_to_rad(parallax)) + + # Mean obliquity of the ecliptic (IAU 1980). + epsilon = 23.4393 - 0.0130 * t + + lambda_r = deg_to_rad(lambda) + beta_r = deg_to_rad(beta) + eps_r = deg_to_rad(epsilon) + + sin_dec = + :math.sin(beta_r) * :math.cos(eps_r) + + :math.cos(beta_r) * :math.sin(eps_r) * :math.sin(lambda_r) + + dec = rad_to_deg(:math.asin(sin_dec)) + + ra = + (:math.sin(lambda_r) * :math.cos(eps_r) - :math.tan(beta_r) * :math.sin(eps_r)) + |> :math.atan2(:math.cos(lambda_r)) + |> rad_to_deg() + |> normalize_360() + + %{ra: ra, dec: dec, distance_km: distance_km} + end + + @doc """ + Topocentric Moon position for an observer at `(lat, lon)` degrees, + UTC `dt`. Returns azimuth ∈ [0, 360), elevation ∈ [-90, 90], and + the geocentric slant range (the topocentric correction is below + the module's stated accuracy floor). + """ + @spec observer_position(DateTime.t(), float(), float()) :: + %{azimuth: float(), elevation: float(), distance_km: float()} + def observer_position(%DateTime{} = dt, lat_deg, lon_deg) do + g = geocentric(dt) + + jd = julian_day(dt) + t = (jd - 2_451_545.0) / 36_525.0 + + # Greenwich Mean Sidereal Time — USNO polynomial (degrees). + gmst = + normalize_360( + 280.46061837 + 360.98564736629 * (jd - 2_451_545.0) + 0.000387933 * t * t - + t * t * t / 38_710_000.0 + ) + + lst = normalize_360(gmst + lon_deg) + hour_angle = normalize_360(lst - g.ra) + + lat_r = deg_to_rad(lat_deg) + dec_r = deg_to_rad(g.dec) + ha_r = deg_to_rad(hour_angle) + + sin_el = + :math.sin(lat_r) * :math.sin(dec_r) + + :math.cos(lat_r) * :math.cos(dec_r) * :math.cos(ha_r) + + elevation = rad_to_deg(:math.asin(sin_el)) + + azimuth = + -:math.sin(ha_r) + |> :math.atan2(:math.tan(dec_r) * :math.cos(lat_r) - :math.sin(lat_r) * :math.cos(ha_r)) + |> rad_to_deg() + |> normalize_360() + + %{azimuth: azimuth, elevation: elevation, distance_km: g.distance_km} + end + + @doc """ + Radial velocity of the Moon relative to an observer at `(lat, lon)` + in km/s at UTC instant `dt`. Positive = the Moon is moving away + (increasing slant range), negative = approaching. + + Computed via centred finite difference on the geocentric distance, + which tracks the dominant rate (the Earth-Moon elliptical orbit + plus the observer's diurnal baseline contribution is below 100 m/s). + Accurate enough to size the ±kHz Doppler shift EME operators care + about at the VHF/microwave bands this project supports. + """ + @spec radial_velocity_km_s(DateTime.t(), float(), float()) :: float() + def radial_velocity_km_s(%DateTime{} = dt, lat_deg, lon_deg) do + before = DateTime.add(dt, -1, :second) + later = DateTime.add(dt, 1, :second) + + d_before = observer_position(before, lat_deg, lon_deg).distance_km + d_later = observer_position(later, lat_deg, lon_deg).distance_km + + # Centred difference in km per 2 s → km/s. + (d_later - d_before) / 2.0 + end + + # --- Helpers -------------------------------------------------------- + + defp deg_to_rad(d), do: d * :math.pi() / 180.0 + defp rad_to_deg(r), do: r * 180.0 / :math.pi() + + defp normalize_360(x) do + r = :math.fmod(x, 360.0) + if r < 0, do: r + 360.0, else: r + end +end diff --git a/lib/microwaveprop/propagation/eme.ex b/lib/microwaveprop/propagation/eme.ex new file mode 100644 index 00000000..5c2d6bbc --- /dev/null +++ b/lib/microwaveprop/propagation/eme.ex @@ -0,0 +1,197 @@ +defmodule Microwaveprop.Propagation.Eme do + @moduledoc """ + Earth-Moon-Earth (moonbounce) link-budget helpers. + + EME path loss follows the bistatic radar equation with the Moon + acting as a limb-darkened sphere of albedo ρ ≈ 0.065. That factors + into a single RCS σ_moon ≈ 0.065 · π · R_moon² ≈ 6.16 × 10¹¹ m². + + Reducing the full radar equation to `f_MHz` and `d_km` gives the + closed form used by this module: + + L_eme(dB) = -14.44 + 40·log₁₀(d_km) + 20·log₁₀(f_MHz) + + Spot-check values (round-trip, one leg Earth→Moon, one leg + Moon→Earth): + + - 144 MHz at 384_400 km → ~252 dB + - 432 MHz at 384_400 km → ~261 dB + - 1_296 MHz at 384_400 km → ~271 dB + - 10_368 MHz at 384_400 km → ~289 dB + + Pure functions only — no DB, no observer state. + """ + + @moon_radius_km 1738.1 + @default_albedo 0.065 + + @doc """ + Band-dependent lunar radar albedo ρ (dimensionless). + + Measured / cited values from the EME literature (Evans 1969, + Hagfors, and later K1JT / W5LUA handbook values): + + | Band | ρ | + |------------------|------| + | < 300 MHz | 0.065| + | 300–999 MHz | 0.081| + | 1–1.999 GHz | 0.065| + | 2–4.999 GHz | 0.070| + | 5–9.999 GHz | 0.080| + | 10–23.999 GHz | 0.090| + | 24–46.999 GHz | 0.060| + | 47–75.999 GHz | 0.030| + | ≥ 76 GHz | 0.020| + + The tiers reflect the surface-roughness regime the radar beam + probes: meter-wavelength bands see quasi-specular returns off + smooth mare, mid-centimetre bands see a rough surface, and + millimetre bands start losing reflectivity to thermal emission. + """ + @spec moon_albedo(number()) :: float() + def moon_albedo(freq_mhz) when is_number(freq_mhz) and freq_mhz < 300, do: 0.065 + def moon_albedo(freq_mhz) when freq_mhz < 1000, do: 0.081 + def moon_albedo(freq_mhz) when freq_mhz < 2000, do: 0.065 + def moon_albedo(freq_mhz) when freq_mhz < 5000, do: 0.070 + def moon_albedo(freq_mhz) when freq_mhz < 10_000, do: 0.080 + def moon_albedo(freq_mhz) when freq_mhz < 24_000, do: 0.090 + def moon_albedo(freq_mhz) when freq_mhz < 47_000, do: 0.060 + def moon_albedo(freq_mhz) when freq_mhz < 76_000, do: 0.030 + def moon_albedo(freq_mhz) when freq_mhz >= 76_000, do: 0.020 + + @doc """ + Returns the Moon's radar cross-section σ in m² for a given albedo. + Defaults to 0.065 (the mean value commonly quoted at 1.3 GHz) so + callers that don't care about band-dependence get the classic + baseline. + """ + @spec moon_rcs_m2(number()) :: float() + def moon_rcs_m2(albedo \\ @default_albedo) when is_number(albedo) and albedo > 0 do + radius_m = @moon_radius_km * 1000.0 + albedo * :math.pi() * radius_m * radius_m + end + + @doc """ + Free-space round-trip path loss (Earth → Moon → Earth) in dB, + without any Moon-reflection term. This is just 2 × one-way FSPL: + + L_fs,2way(dB) = 64.90 + 40·log₁₀(f_MHz) + 40·log₁₀(d_km) + + Pair with `moon_reflection_gain_db/2` to decompose the full EME + loss into its spreading vs scattering terms. + """ + @spec fspl_round_trip_db(number(), number()) :: float() + def fspl_round_trip_db(freq_mhz, distance_km) + when is_number(freq_mhz) and freq_mhz > 0 and is_number(distance_km) and distance_km > 0 do + 64.90 + 40.0 * :math.log10(freq_mhz) + 40.0 * :math.log10(distance_km) + end + + @doc """ + The Moon's equivalent "antenna gain" in dB for a monostatic radar: + `G_moon = 10·log₁₀(4π σ / λ²)` where σ = ρ · π · R². This is the + term that *reduces* the net round-trip path loss relative to plain + double-FSPL — i.e. the moon's reflection effectiveness. + + `moon_reflection_gain_db/2` lets callers override the albedo (for + example via `moon_albedo/1`); `moon_reflection_gain_db/1` defaults + to ρ = 0.065. + """ + @spec moon_reflection_gain_db(number(), number()) :: float() + def moon_reflection_gain_db(freq_mhz, albedo \\ @default_albedo) + when is_number(freq_mhz) and freq_mhz > 0 and is_number(albedo) and albedo > 0 do + lambda_m = 300.0 / freq_mhz + sigma = moon_rcs_m2(albedo) + 10.0 * :math.log10(4.0 * :math.pi() * sigma / (lambda_m * lambda_m)) + end + + @doc """ + Round-trip EME path loss in dB for a frequency in MHz and an + Earth→Moon slant range in km. Covers the transmit leg, Moon + scattering, and receive leg in a single number (the conventional + "moonbounce path loss" that link-budget spreadsheets use). + + Accepts an optional `albedo` override — pass + `moon_albedo(freq_mhz)` for a band-aware calculation. Defaults to + ρ = 0.065, the classic 1296 MHz value, which keeps legacy callers + stable. + + Equivalent to `fspl_round_trip_db - moon_reflection_gain_db`. + """ + @spec path_loss_db(number(), number(), number()) :: float() + def path_loss_db(freq_mhz, distance_km, albedo \\ @default_albedo) + when is_number(freq_mhz) and freq_mhz > 0 and is_number(distance_km) and distance_km > 0 and is_number(albedo) and + albedo > 0 do + fspl_round_trip_db(freq_mhz, distance_km) - moon_reflection_gain_db(freq_mhz, albedo) + end + + # Speed of light (m/s). + @c_m_s 299_792_458.0 + + @doc """ + Self-echo Doppler shift in Hz for an EME bounce, given the operating + frequency and the Moon's radial velocity relative to the observer. + + For an own-echo round-trip the Doppler is twice the one-way shift: + `Δf = −2 · f · v_radial / c`. A positive radial velocity (Moon + receding) produces a negative Hz shift — the returned bounce lands + *below* the transmit frequency. + """ + @spec doppler_shift_hz(number(), number()) :: float() + def doppler_shift_hz(freq_mhz, radial_velocity_km_s) + when is_number(freq_mhz) and freq_mhz > 0 and is_number(radial_velocity_km_s) do + freq_hz = freq_mhz * 1_000_000.0 + v_m_s = radial_velocity_km_s * 1000.0 + -2.0 * freq_hz * v_m_s / @c_m_s + end + + @doc """ + Received signal power in dBm for a symmetric own-echo EME link + (operator listens to their own bounce) given: + + - `:tx_power_dbm` — transmitter output power (dBm) + - `:tx_gain_dbi` — transmit antenna gain (dBi) + - `:rx_gain_dbi` — receive antenna gain (dBi) — same as tx for + self-echo, kept distinct so the function also works for cross- + station contacts + - `:freq_mhz` — operating frequency (MHz) + - `:distance_km` — one-way slant range to the Moon (km) + """ + @spec received_power_dbm(keyword()) :: float() + def received_power_dbm(opts) do + tx = Keyword.fetch!(opts, :tx_power_dbm) + gt = Keyword.fetch!(opts, :tx_gain_dbi) + gr = Keyword.fetch!(opts, :rx_gain_dbi) + f = Keyword.fetch!(opts, :freq_mhz) + d = Keyword.fetch!(opts, :distance_km) + + tx + gt + gr - path_loss_db(f, d) + end + + @doc """ + Thermal noise floor (dBm) in a given bandwidth and system noise + temperature. Callers use this against `received_power_dbm/1` to + get SNR. + + Default Tsys = 290 K (room-temperature receive system, a reasonable + placeholder for VHF; a real cooled LNA at 10 GHz is closer to 100 K). + """ + @spec noise_floor_dbm(number(), number()) :: float() + def noise_floor_dbm(bandwidth_hz, t_sys_k \\ 290.0) + when is_number(bandwidth_hz) and bandwidth_hz > 0 and is_number(t_sys_k) and t_sys_k > 0 do + # -174 dBm/Hz + 10 log(BW) + 10 log(T / 290) + -174.0 + 10.0 * :math.log10(bandwidth_hz) + 10.0 * :math.log10(t_sys_k / 290.0) + end + + @doc """ + Signal-to-noise ratio (dB) for a symmetric own-echo EME link. + Convenience on top of `received_power_dbm/1` and + `noise_floor_dbm/2` — returns the dB margin over thermal noise in + the specified detection bandwidth. + """ + @spec snr_db(keyword()) :: float() + def snr_db(opts) do + bw_hz = Keyword.fetch!(opts, :bandwidth_hz) + t_sys_k = Keyword.get(opts, :t_sys_k, 290.0) + received_power_dbm(opts) - noise_floor_dbm(bw_hz, t_sys_k) + end +end diff --git a/lib/microwaveprop_web/components/layouts.ex b/lib/microwaveprop_web/components/layouts.ex index 3da7b8c3..caad3630 100644 --- a/lib/microwaveprop_web/components/layouts.ex +++ b/lib/microwaveprop_web/components/layouts.ex @@ -42,6 +42,7 @@ defmodule MicrowavepropWeb.Layouts do