prop/docs/algo_original.md
Graham McIntire a283ad9c66
Refactor HRRR fetch to batch points per hour
QSO enrichment now groups all path points by HRRR hour and creates
one batch job per hour instead of one job per point. The batch job
downloads the GRIB2 data once and extracts all needed points from
the same binary. Legacy single-point jobs are still supported for
backward compatibility.
2026-03-30 17:21:47 -05:00

16 KiB
Raw Permalink Blame History

Comprehensive Propagation Algorithm Documentation

Overview

The ntms-proptools suite consists of three interconnected applications for mm-wave (24/47 GHz) propagation analysis and prediction. All share common physics algorithms calibrated against confirmed radio contacts.

Calibration Contacts:

  • 47 GHz Jun 01 2024 98.8 km (09:00 CDT)
  • 47 GHz Nov 01 2025 116.0 km (07:44 CST)
  • 24 GHz Nov 01 2025 101.0 km (07:54 CST)
  • 24 GHz Sep 07 2002 542.1 km (longest, 07:35 CDT)

PART 1: PROP-CAST (Real-Time Station Scoring)

Purpose

Live multi-station propagation scoring across ~500-mile DFW radius using current ASOS observations.

Data Inputs

Surface Observations (ASOS/NWS):

  • Temperature (°F → °C conversion: (tF - 32) * 5/9)
  • Dew point (°F → °C)
  • Relative humidity (%)
  • Wind speed (knots → m/s: knots * 1.15078)
  • Wind direction (degrees)
  • Sea-level pressure (hPa)
  • Altimeter setting (inches Hg → hPa: inches * 33.8639)
  • Cloud cover (categorical: CLR, FEW, SCT, BKN, OVC, VV)

Frequency Configurations:

// 47 GHz Configuration
47: {
  label: "47 GHz",
  accentHex: "#00ffa3",
  attenPerKm: (rho) => 0.006 * rho + 0.040,
  humPenalty: 1.0,
  seasonAdj: { 6:-0, 7:-0, 8:-0, 5:-0, 9:-0 },
  // Atmospheric window: between 22.235 GHz H₂O peak and 60 GHz O₂ peak
}

// 24 GHz Configuration
24: {
  label: "24 GHz",
  accentHex: "#44ddff",
  attenPerKm: (rho) => 0.0067 * rho + 0.006,
  humPenalty: 1.6,
  seasonAdj: { 6:-8, 7:-10, 8:-10, 5:-4, 9:-4 },
  // Adjacent to 22.235 GHz H₂O absorption line shoulder
}

Physics Calculations

1. Absolute Humidity Calculation

function absHumidity(rhPct, tempC) {
  // Magnus formula for saturation vapor pressure
  const es = 6.112 * Math.exp((17.67 * tempC) / (tempC + 243.5));
  // Absolute humidity in g/m³
  return (217 * (rhPct/100) * es) / (tempC + 273.15);
}

Formula: ρ = 217 * (RH/100) * e_s / T

  • RH: relative humidity (%)
  • e_s: saturation vapor pressure (hPa) via Magnus formula
  • T: absolute temperature (K)
  • Result: ρ in g/m³

2. Gaseous Absorption by Frequency

attenPerKm(rho) = O₂_coefficient + H₂O_coefficient * rho

47 GHz: 0.040 + 0.006·ρ dB/km

  • O₂ component: ~0.040 dB/km (from 60 GHz O₂ line wing)
  • H₂O component: ~0.006·ρ dB/km (22.235 GHz H₂O line wing)

24 GHz: 0.006 + 0.0067·ρ dB/km

  • O₂ component: ~0.006 dB/km
  • H₂O component: ~0.0067·ρ dB/km (stronger H₂O sensitivity on 22.235 GHz shoulder)

Scoring Functions

All return 0-100 scale representing propagation favorability.

Score: Humidity

function scoreHumidity(rho, freq) {
  const r = rho * FREQ_CONFIG[freq].humPenalty;
  if (r <= 6)   return 100;
  if (r <= 9)   return 95  - ((r-6)/3)  * 20;
  if (r <= 13)  return 75  - ((r-9)/4)  * 30;
  if (r <= 18)  return 45  - ((r-13)/5) * 35;
  return Math.max(0, 10 - (r-18)*2);
}

Penalty Factor per Frequency:

  • 47 GHz: multiplier = 1.0
  • 24 GHz: multiplier = 1.6 (more sensitive due to H₂O line proximity)

Thresholds:

  • r ≤ 6: Score 100 (dry/excellent)
  • r ≤ 9: Linear decay from 95 → 75 (excellent → good)
  • r ≤ 13: Linear decay from 75 → 45 (good → marginal)
  • r ≤ 18: Linear decay from 45 → 10 (marginal → poor)
  • r > 18: Score 0 (very humid/negligible)

Score: Wind Speed

function scoreWind(mph) {
  if (mph < 2)  return 100;
  if (mph < 5)  return 90;
  if (mph < 9)  return 72;
  if (mph < 14) return 42;
  if (mph < 20) return 16;
  return 0;
}

Rationale: Wind causes mechanical mixing of boundary layer, destroying temperature inversion. Still air promotes stable inversion.

Score: Sky Cover

function scoreSky(c) {
  return ({
    CLR: 100, SKC: 100,
    FEW: 88,
    SCT: 60,
    BKN: 25,
    OVC: 5, VV: 5
  })[c] ?? 50;
}

METAR Cloud Cover Legend:

  • CLR/SKC: Clear (100)
  • FEW (1/82/8): 88
  • SCT (3/84/8): 60 (scattering)
  • BKN (5/87/8): 25 (broken)
  • OVC (8/8): 5 (overcast)
  • VV (vertical visibility): 5

Score: Season (Month-Based)

function scoreSeason(month, freq) {
  const b = ({
    1:88, 2:84, 3:72, 4:62, 5:55,
    6:42, 7:28, 8:28, 9:52, 10:68,
    11:96, 12:88
  })[month] ?? 50;
  return Math.max(0, Math.min(100, b + FREQ_CONFIG[freq].seasonAdj[month]??0));
}

Base Monthly Scores:

  • Peak: November (96), December (88), January (88)
  • Summer worst: July/August (28)
  • 24 GHz additional penalties: June (-8), July (-10), August (-10), May (-4), September (-4)

Rationale: Seasonal humidity variation; summer monsoon moisture increases H₂O absorption.

Score: Time of Day (Inversion-Based)

function scoreTimeOfDay(utcH, utcM, month) {
  const isCDT = month >= 3 && month <= 10;
  const local = (utcH + utcM/60 + (isCDT ? -5 : -6) + 24) % 24;
  const sr = [7.4,7.3,7.0,6.7,6.35,6.25,6.35,6.65,6.9,7.1,7.35,7.45][month-1];
  const d = local - sr;  // hours relative to sunrise

  if (d >= -0.5 && d <= 1.5) return {s:100, lbl:"Peak — inversion maximum"};
  if (d > 1.5 && d <= 3.0)   return {s:78,  lbl:"Good — inversion eroding"};
  if (d > -1.5 && d < -0.5)  return {s:82,  lbl:"Pre-sunrise — building"};
  if (d > 3.0 && d <= 6.0)   return {s:38,  lbl:"Marginal — BL mixing"};
  if (local >= 21 || local <= 3) return {s:68,lbl:"Overnight — cooling"};
  if (d > 6.0) return {s:18,lbl:"Afternoon — duct unlikely"};
  return {s:48,lbl:"Early morning"};
}

Window Tiers:

  1. Peak Window (sunrise -30min to +1.5h): Score 100 — Inversion maximum strength
  2. Good Window (+1.5h to +3h after sunrise): Score 78 — Inversion still present but eroding
  3. Pre-Sunrise (-1.5h to -30 min): Score 82 — Inversion building under clear skies
  4. Marginal Window (+3h to +6h after sunrise): Score 38 — Convective boundary layer growth
  5. Overnight Cooling (21:00 to 03:00 local): Score 68 — Weak inversion reformation
  6. Afternoon (+6h after sunrise): Score 18 — Full convective mixing, inversion destroyed

Sunrise Lookup Table (Monthly, local decimal hours):

Jan:7.4  Feb:7.3  Mar:7.0  Apr:6.7  May:6.35  Jun:6.25
Jul:6.35 Aug:6.65 Sep:6.9  Oct:7.1  Nov:7.35  Dec:7.45

Score: Temperature-Dew Point Depression

function scoreTdDep(tF, tdF) {
  const d = tF - tdF;
  if (d > 22) return 96;
  if (d > 14) return 80;
  if (d > 8)  return 60;
  if (d > 4)  return 38;
  return 18;
}

Rationale: Large Td depression indicates dry air above surface, favorable for stable stratification.

Thresholds (°F):

  • 22°F: 96 (very dry aloft)

  • 14°F: 80 (good stability)

  • 8°F: 60 (moderate)

  • 4°F: 38 (marginal)

  • ≤4°F: 18 (humid aloft)

Score: Pressure Tendency

function scorePressure(mb, prev) {
  if (prev == null) return 58;
  const d = mb - prev;
  if (d > 2.5) return 98;  // Rising pressure → anticyclone → stable
  if (d > 0.8) return 82;
  if (d > -0.5) return 60;
  if (d > -2.0) return 32;
  return 10;              // Falling pressure → cyclone → unstable
}

Composite Propagation Score

function composite(f) {
  return Math.round(
    f.humidity  * 0.26 +
    f.wind      * 0.18 +
    f.sky       * 0.15 +
    f.timeOfDay * 0.14 +
    f.tdDep     * 0.11 +
    f.season    * 0.09 +
    f.pressure  * 0.07
  );
}

Weights:

Factor Weight Importance
Humidity 26% Most critical (absorption)
Wind 18% Destroys inversion
Sky 15% Less impact than humidity
Time of Day 14% Inversion strength varies
Td Depression 11% Aloft stability
Season 9% Long-term trend
Pressure 7% Weak indicator

Total: 100%

Scoring Tiers

Score Band Label Hex
80-100 EXCELLENT Green/cyan #00ffa3
65-79 GOOD Light cyan #7dffd4
50-64 MARGINAL Yellow #ffe566
33-49 POOR Orange #ff9044
0-32 NEGLIGIBLE Red #ff4f4f

PART 2: PROP-ANALYZER (Historical Contact Analysis)

Purpose

Analyze atmospheric conditions (ASOS + Skew-T sounding) during historical radio contacts to identify propagation mechanisms.

Additional Physics Calculations

Refractivity Profile (Modified N-Units)

// Saturation vapor pressure (Buck equation)
function satVapPres(tC) {
  return 6.1121 * Math.exp((18.678 - tC/234.5) * (tC/(257.14 + tC)));
}

// N-unit profile (ITU-R P.453-14)
const N = 77.6 * P/T + 3.73e5 * e / (T * T)

// M-unit profile (accounts for earth curvature)
const M = N + 0.157 * (hght - sfc_hght)

Where:

  • P: pressure (hPa)
  • T: absolute temperature (K)
  • e: water vapor pressure (hPa) = satVapPres(Td)
  • h: height MSL (m)
  • sfc_hght: surface height MSL (m)

Duct Detection via M-Profile

Definition: A duct exists where dM/dh < 0 (M decreases with height)

const ducts = [];
let inDuct = false, ductBase = null, ductBaseM = null;

for (let i = 1; i < refractProfile.length; i++) {
  const dM = refractProfile[i].M - refractProfile[i-1].M;

  if (dM < 0 && !inDuct) {
    inDuct = true;
    ductBase = refractProfile[i-1].hAGL;
    ductBaseM = refractProfile[i-1].M;
  } else if (dM >= 0 && inDuct) {
    inDuct = false;
    const top = refractProfile[i-1].hAGL;
    const strength = ductBaseM - refractProfile[i-1].M;
    if (strength > 2) {  // Filter: only significant ducts
      ducts.push({ base: ductBase, top, strength: strength.toFixed(1) });
    }
  }
}

M-Unit Gradient Interpretation:

  • dM/dh < -157 /km: Super-refractive (ducting probable)
  • dM/dh < -79 /km: Enhanced refraction (extended range)
  • dM/dh > 0: Sub-refractive (reduced range)
  • dM/dh ≈ -40 /km: Standard atmosphere

Inversion Detection

const inversions = [];

for (let i = 1; i < sorted.length; i++) {
  if (sorted[i].tmpc > sorted[i-1].tmpc) {  // Temperature increases with height
    inversions.push({
      base: sorted[i-1].hght - sfc.hght,
      top: sorted[i].hght - sfc.hght,
      strength: sorted[i].tmpc - sorted[i-1].tmpc,  // °C
      rate: strength / ((top - base) / 100)  // °C per 100m
    });
  }
}

// Merge adjacent inversions within 200m gap
// Keep only meaningful: strength >= 0.5°C and base < 5000m

Stability Indices

K-Index:

const kIndex = (T850 - T500) + Td850 - (T700 - Td700)

Lifted Index (Simplified):

const li = T500 - (Tsfc - (p500_hght - sfc_hght) * 0.00976)
  • LI < 0: Unstable (convection likely)
  • LI > 0: Stable (inversion maintained)

Precipitable Water:

PW = Σ [(MR_i + MR_{i-1})/2 * ΔP / (9.81 * 10)]
// MR = mixing ratio = 622 * e / (P - e)

Boundary Layer Depth

// Potential temperature
theta = T + (9.8 * h / 1000)

// Find height where theta > theta_sfc + 2°C

Purpose

Complete point-to-point RF link analysis including path clearance, Fresnel zone obstruction, and margin calculations.

Path Geometry

Maidenhead Grid Decode: Field (20°×10°) → Square (2°×1°) → Subsquare (5'×2.5') → Extended

Distance: Great-circle haversine

Path Clearance Physics

Fresnel Zone Radius

function fresnelRadius(d1m, d2m, lambdaM) {
  return Math.sqrt(lambdaM * d1m * d2m / (d1m + d2m));
}
// λ = 0.3 / f_GHz

Earth Bulge Correction

function earthBulge(frac, distKm, K = 4/3) {
  const d1 = frac * distKm * 1000;
  const d2 = (1 - frac) * distKm * 1000;
  return (d1 * d2) / (2 * K * 6371000);
}

Effective K-factor from Surface Refractivity:

function effectiveKFactor(N_surface) {
  const N_std = 315;
  const dN_std = -40;
  const dN_est = dN_std - (N_surface - N_std) * 0.25;
  const Ke = 1 / (1 + 6371 * dN_est * 1e-6);
  return Math.max(0.5, Math.min(5.0, Ke));
}

function surfaceRefractivity(tC, rhPct, mslpHpa) {
  const T = tC + 273.15;
  const P = mslpHpa ?? 1013.25;
  const es = 6.112 * Math.exp(17.67 * tC / (tC + 243.5));
  const e = (rhPct / 100) * es;
  return (77.6 * P / T) + (3.73e5 * e / (T * T));
}

Knife-Edge Diffraction

function knifeEdgeLoss(v) {
  if (v <= -0.7787) return 0;      // Clear
  if (v <= 0) return -20 * Math.log10(0.5 - 0.62 * v);
  if (v <= 1) return -20 * Math.log10(0.5 * Math.exp(-0.95 * v));
  if (v <= 2.4) {
    const inner = Math.max(0, 0.1184 - (0.38 - 0.1 * v) ** 2);
    return -20 * Math.log10(0.4 - Math.sqrt(inner));
  }
  return 20 * Math.log10(v) + 13.0;  // Asymptotic
}

Diffraction Verdicts:

  • v < -0.7787: CLEAR (0 dB loss)
  • -0.7787 < v < 0: FRESNEL_MINOR (< 3 dB)
  • 0 < v < 2.4: FRESNEL_PARTIAL (3-10 dB)
  • v > 2.4: BLOCKED (>10 dB)

Free-Space Path Loss (ITU-R P.525)

function fspl(distKm, fGhz) {
  return 20 * Math.log10(distKm) + 20 * Math.log10(fGhz) + 92.45;
}

Gaseous Absorption (ITU-R P.676)

function gasLoss(distKm, rhoGm3, freq) {
  return (cfg.o2dBkm + cfg.h2oCoeff * rhoGm3) * distKm;
}
Frequency O₂ (dB/km) H₂O Coeff Notes
47 GHz 0.040 0.006 Clear window
24 GHz 0.006 0.0067 H₂O shoulder
10 GHz 0.003 0.0005 Clear window

EIRP

eirp = pDbm + gtxDbi - feedDb

Receiver Sensitivity

sensitivity = -174 + nfDb + 10 * Math.log10(bw)
// CW: bw=500 Hz, SSB: bw=2700 Hz

Duct Enhancement Estimate

function ductEnhancement(propScore) {
  if (propScore >= 80) return -14;  // 14 dB improvement
  if (propScore >= 65) return -10;
  if (propScore >= 50) return -6;
  if (propScore >= 33) return -2;
  return 0;
}

Calibrated against confirmed contacts.

Success Probability

function marginToSuccess(marginDb, propScore) {
  // Margin factor (piecewise linear)
  let marginPct;
  if (marginDb <= 0)       marginPct = 0;
  else if (marginDb <= 10) marginPct = (marginDb / 10) * 20;
  else if (marginDb <= 15) marginPct = 20 + ((marginDb - 10) / 5) * 20;
  else if (marginDb <= 20) marginPct = 40 + ((marginDb - 15) / 5) * 20;
  else if (marginDb <= 25) marginPct = 60 + ((marginDb - 20) / 5) * 20;
  else if (marginDb <= 30) marginPct = 80 + ((marginDb - 25) / 5) * 20;
  else marginPct = 100;

  // Propagation factor (±30% modulation)
  const propFactor = 0.70 + (propScore / 100) * 0.60;

  return Math.max(0, Math.min(99, Math.round(marginPct * propFactor)));
}

Margin → Success Calibration:

  • 0 dB → 0%, 10 dB → 20%, 15 dB → 40%, 20 dB → 60%, 25 dB → 80%, 30+ dB → 100%

Propagation Modulation:

  • Score 100 (excellent): ×1.30
  • Score 50 (neutral): ×1.00
  • Score 0 (poor): ×0.70

Integration: Data Flow

Weather Data (ASOS/NWS)
    ↓
[PropCast] → Composite Score 0-100
    ↓
    ├→ Real-time station dashboard
    └→ Input to LinkBudget (propScore parameter)

Sounding Data (IEM RAOB)
    ↓
[PropAnalyzer] → Refractivity profiles, ducts, inversions
    ↓
    ├→ Skew-T visualization
    ├→ dN/dh gradient analysis
    └→ Mechanism identification

Path Geometry + Elevation
    ↓
[LinkBudget] → Profile analysis, Fresnel zones, clearance
    ↓
    ├→ Diffraction loss calculation
    ├→ Link margin = RX Power - Sensitivity
    └→ Success % = marginToSuccess(margin, propScore)

Constants & Thresholds Reference

Constant Value Source
Humidity penalty 47 GHz 1.0 Minimal H₂O line sensitivity
Humidity penalty 24 GHz 1.6 Adjacent 22.235 GHz H₂O peak
Duct M-unit min depth 2 Noise filter
Inversion min strength 0.5°C Below is noise
Inversion height limit 5000m AGL Above irrelevant
Earth radius 6371 km WGS-84 mean
Standard K-factor 4/3 Standard atmosphere
Standard surface N 315 ITU-R P.453
Standard dN/dh -40 /km ITU-R P.453

ITU-R References

  • P.453-14: Radio refractivity
  • P.525: Free-space path loss
  • P.676: Gaseous absorption
  • P.526: Diffraction