MUF, LUF and Reliability Optimization in Python
Anyone who has operated on shortwave knows the feeling: at noon the 20 m band roars with signals, while at midnight the same band sounds like static. HF sky-wave propagation depends on the ionosphere, and the ionosphere depends on the Sun. Choosing a frequency is therefore not a one-time decision. It is an optimization problem that changes every hour.
In this article we build a compact but physically meaningful model of a single-hop HF link. We then find the frequency that maximizes link reliability for every hour of the day.
The Example Problem
Consider a 2,500 km sky-wave link over a mid-latitude path (midpoint at 25°N) at the equinox, with a smoothed sunspot number $R_{12}=100$.
| Item | Value |
|---|---|
| Transmit power | 100 W (20 dBW) |
| Antenna gain (TX / RX) | 3 dBi / 3 dBi |
| Mode / bandwidth | SSB voice, 2.7 kHz |
| Required SNR | 10 dB |
| Reflecting layer | F2 layer, virtual height 300 km |
| Search range | 3 to 30 MHz |
Goal: for every hour $t$, find the operating frequency $f^*(t)$ that maximizes the probability that the link works.
The Mathematical Model
1. Geometry of a single hop
With Earth radius $R_E$, ground distance $D$ and reflection height $h$, the half central angle is
$$\alpha = \frac{D}{2R_E}$$
The take-off (elevation) angle $\Delta$ and the incidence angle $i$ at the layer are
$$\tan\Delta = \frac{\cos\alpha - \dfrac{R_E}{R_E+h}}{\sin\alpha}, \qquad \sin i = \frac{R_E\cos\Delta}{R_E+h}$$
The one-hop path length is
$$L = 2\sqrt{R_E^2 + (R_E+h)^2 - 2R_E(R_E+h)\cos\alpha}$$
The maximum usable frequency follows from the secant law:
$$\mathrm{MUF} = f_oF_2 \cdot \sec i$$
2. Day–night behavior of the F2 layer
The solar zenith angle $\chi$ at the path midpoint (latitude $\varphi$, declination $\delta$, hour angle $H = 15^\circ (t-12)$) satisfies
$$\cos\chi = \sin\varphi\sin\delta + \cos\varphi\cos\delta\cos H$$
The critical frequency blends a nighttime and a daytime value through a smooth switch:
$$S(t) = \frac{1}{2}\left[1+\tanh\left(3\cos\chi\right)\right], \qquad f_oF_2 = f_n + (f_d - f_n),S(t)$$
$$f_d = 3.5 + 0.05R_{12}, \qquad f_n = 2.0 + 0.015R_{12}$$
3. D-layer absorption
Lower frequencies are absorbed in the daytime D layer. An empirical form is
$$L_a(f) = \frac{677.2, I, \sec i_D}{(f+f_H)^2 + 10.2}, \qquad I = (1+0.0037R_{12})\left[\cos^{1.3}(0.881\chi) + 0.02\right]$$
Here $\chi$ is in degrees and capped at $102^\circ$, $f_H$ is an effective gyro-frequency, and $i_D$ is the incidence angle at the D layer (about 90 km).
4. Link budget
Free-space loss (with $f$ in MHz and $L$ in km) is
$$L_{fs} = 32.45 + 20\log_{10} f + 20\log_{10} L$$
The external noise figure decreases with frequency:
$$F_a = c - d\log_{10} f, \qquad N = F_a + 10\log_{10} B - 204 \ \ [\mathrm{dBW}]$$
Combining everything gives the signal-to-noise ratio:
$$\mathrm{SNR}(f,t) = P_t + G_t + G_r - L_{fs} - L_a - L_o - N$$
where $L_o$ lumps ground reflection, polarization and other losses.
5. Reliability as the objective function
Two independent things must go right for the link to work. First, the ionosphere must actually reflect the wave, meaning $f$ stays below the day-to-day MUF, which fluctuates by roughly $\sigma_M$ (a fraction of MUF). Second, the SNR must exceed the requirement, with fading spread $\sigma_S$. With $\Phi$ the standard normal CDF, the reliability is
$$R(f,t) = \Phi!\left(\frac{\mathrm{MUF}(t) - f}{\sigma_M,\mathrm{MUF}(t)}\right)\cdot \Phi!\left(\frac{\mathrm{SNR}(f,t) - \mathrm{SNR}_{\mathrm{req}}}{\sigma_S}\right)$$
The optimization problem is then
$$f^*(t) = \underset{3,\mathrm{MHz},\le, f,\le, 30,\mathrm{MHz}}{\arg\max}; R(f,t)$$
The first factor falls as $f$ approaches the MUF, while the second factor rises with $f$ because absorption and noise both drop. Their product has a clear peak.
Full Source Code
1 | import time |
Code Walkthrough
Section 1: Scenario parameters
All physical constants and design choices live at the top, so you can play with them easily. Change D_KM to test a different distance, R12 to move between solar minimum and maximum, or PT_DBW to see how much a linear amplifier really buys you. SIGMA_MUF and SIGMA_SNR control how conservative the optimizer will be: larger values mean less predictable propagation, which pushes the optimum further from the MUF.
Section 2: Geometry
hop_geometry implements the spherical-Earth formulas above. It returns three values: the total one-hop path length, the take-off angle, and $\sec i$ at the reflection layer. arctan2 is used instead of arctan so the quotient is handled safely. For a 2,500 km path the take-off angle is only about 7.5°, and $\sec i\approx 3.1$. This is why the MUF is roughly three times the critical frequency: a low take-off angle means a grazing incidence on the layer, so much higher frequencies can still be reflected.
SEC_D reuses the same take-off angle but evaluates the incidence angle at 90 km. That number tells us how obliquely the ray crosses the absorbing D layer.
Section 3: Ionosphere and link-budget model
Each formula from the mathematical section is one small function. The important design decision is that every function accepts scalars or NumPy arrays that broadcast against each other. There is no if statement and no explicit loop inside them. This is what makes the acceleration in Section 4 possible with zero code duplication.
cos_zenithcomputes $\cos\chi$ from latitude, declination and local time.fo_f2blends the day and night critical frequencies with the smooth $\tanh$ switch, so the transition at sunrise and sunset is gradual rather than a hard step.d_layer_absorptionclips $\chi$ at 102° (beyond which the D layer disappears), then applies the $1/[(f+f_H)^2+10.2]$ dependence. The small constant0.02keeps a minimal residual absorption at night.snr_dbassembles the link budget: transmit power plus antenna gains minus free-space loss, absorption and miscellaneous loss, minus the noise power.reliabilitymultiplies the two normal-CDF terms.ndtris SciPy’s fast, vectorized implementation of $\Phi$.
Section 4: Grid search, naive versus accelerated
Two versions of the same computation are provided.
grid_naiveuses a doubleforloop and callsreliabilityonce per (hour, frequency) pair. This is the most direct translation of the math, but Python-level loops are slow, and every scalar call pays NumPy’s function-call overhead.grid_vectorizedpassesfreqs[None, :](shape $1\times F$) andhours[:, None](shape $H\times 1$) to the same function. NumPy broadcasting expands them into a full $H\times F$ table in one shot, using compiled loops internally.
The script times both, prints the speed-up, and checks that the maximum difference between the two results is at floating-point rounding level. The vectorized version is the one used from here on. The naive version is kept only as a benchmark, and it becomes unusable as soon as you add another dimension such as the sunspot number in Section 6.
Section 5: Extracting the optimum
np.argmax(rel, axis=1) finds, for every hour, the column index of the highest reliability. From that index we read the optimal frequency f_star, the peak reliability r_star, and the SNR at the optimum.
The LUF is computed with a boolean mask: a frequency counts as usable when the SNR meets the requirement and the frequency lies below the MUF. argmax on the mask returns the first True, which is the lowest usable frequency. If no frequency qualifies, the hour gets NaN.
The grid has a 0.1 MHz resolution, so as a cross-check we refine the optimum at 03:00, 12:00 and 21:00 with scipy.optimize.minimize_scalar. It searches the continuous interval between 3 MHz and 1.3 times the MUF. The refined values should agree with the grid to within the grid spacing.
Section 6: Sensitivity to solar activity
To see how the answer changes over the solar cycle, we add a third axis. The arrays have shapes $H\times1\times1$, $1\times F\times1$ and $1\times1\times R$, so one call to reliability produces a full $H\times F\times R$ tensor of about 390,000 values. Taking argmax along the frequency axis yields f_star_map, the optimal frequency as a function of hour and $R_{12}$. Running this through the naive loop would be hopeless by comparison.
Section 7: Visualization
Everything is drawn in a single figure with six panels: four 2D plots and two 3D surfaces. plot_surface draws the 3D graphs, pcolormesh draws the heat map, and twinx gives the last panel two vertical axes. The finished image is also saved as a PNG so you can reuse it directly.
Execution Results
Console Output
=== Computation time ===
Evaluations : 13,279
Naive double loop : 1.6166 s
Vectorized (NumPy) : 0.0036 s
Speed-up : 449.3 x
Max abs difference : 6.66e-16
=== Model summary ===
Ground distance : 2500 km
Path length (1 hop) : 2623.6 km
Take-off angle : 7.53 deg
sec(i) at F2 layer : 3.107
=== Optimal frequency (R12 = 100) ===
Hour MUF f* f*/MUF SNR R*
0 10.94 8.30 0.759 18.51 0.849
2 11.01 8.30 0.754 18.51 0.850
4 11.83 8.90 0.752 18.84 0.860
6 18.64 14.40 0.772 18.92 0.858
8 25.45 20.40 0.802 17.77 0.815
10 26.27 21.80 0.830 15.96 0.738
12 26.34 22.10 0.839 15.21 0.703
14 26.27 21.80 0.830 15.96 0.738
16 25.45 20.40 0.802 17.77 0.815
18 18.64 14.40 0.772 18.92 0.858
20 11.83 8.90 0.752 18.84 0.860
22 11.01 8.30 0.754 18.51 0.850
=== Refinement with scipy.optimize.minimize_scalar ===
Hour Grid f* Refined f* Refined R*
3 8.50 8.470 0.8526
12 22.10 22.075 0.7027
21 8.50 8.470 0.8526
Graph Output
Reading the Results
Console output
The first block reports the computation time. Because the two implementations call exactly the same model function, the vectorized version reproduces the loop result to rounding error while running in a small fraction of the time.
The model summary confirms the geometry: a path of about 2,624 km, a take-off angle of roughly 7.5°, and $\sec i\approx 3.107$.
The optimal-frequency table shows the essence of the problem:
- At 00:00, the MUF is about 10.9 MHz and the best frequency is 8.3 MHz, only 76% of the MUF. The SNR is about 18.5 dB and the reliability about 0.85.
- At 12:00, the MUF rises to about 26.3 MHz and the best frequency is 22.1 MHz, or 84% of the MUF. The SNR is about 15.2 dB and the reliability about 0.70.
The ratio $f^*/\mathrm{MUF}$ stays between roughly 0.75 and 0.84 all day. This agrees with the classic operating rule of using about 85% of the MUF (the “frequency of optimum traffic”), but here the ratio comes out of the optimization instead of being assumed. The refinement table confirms that the continuous optimizer lands within 0.03 MHz of the grid result (for example 22.075 MHz against 22.1 MHz at noon).
Panel (1): MUF, LUF and the optimal frequency
The blue MUF curve is a flat plateau of about 11 MHz at night, rises steeply between 05:00 and 08:00 as the F2 layer ionizes, and peaks at about 26 MHz around noon. The orange LUF curve is the mirror image: at noon the D layer absorbs so strongly that frequencies below about 17.3 MHz cannot meet the 10 dB SNR requirement. At night the LUF sits at the 3 MHz floor of the search range, which only means that every frequency in the range meets the SNR target. The green usable window is therefore narrow at noon (roughly 17 to 26 MHz) and wide at night. The red optimal line runs inside this window and always stays below the 0.85 MUF dashed line, tracking the MUF at a safe distance.
Panel (2): 3D reliability surface
The surface shows reliability over time and frequency. It has a clear ridge that follows the MUF curve: high on the left of the ridge (frequency low enough to reflect and strong enough to be heard), and a sharp cliff on the right (frequency above the MUF, where the wave escapes into space). The red line traces the ridge crest. At night the ridge is low in frequency and broad along the frequency axis, and during the day it climbs to high frequencies. At noon its crest is visibly lower than at night, which is the fingerprint of D-layer absorption. The three-dimensional view makes it obvious that choosing a frequency without regard to the time of day would put you either over the cliff or in the valley.
Panel (3): Heat map
The same data viewed from above. The bright band is the region of high reliability, and the cyan dashed MUF marks the edge of the cliff. The white optimal line hugs the crest of the band, just under the cliff. Below the band the color fades gradually because of absorption and noise, while above it the color drops abruptly to zero. This asymmetry is the reason the optimum sits closer to the MUF than to the LUF but never touches it: the penalty for overshooting is far more severe than the penalty for undershooting.
Panel (4): Reliability versus frequency
Four cross sections at 00:00, 06:00, 12:00 and 18:00 make the trade-off concrete. Each curve rises slowly and falls quickly, and the dot marks its peak. The 00:00 curve peaks near 8 MHz, the 18:00 curve near 14 MHz, and the 12:00 curve near 22 MHz. Note how narrow the 00:00 curve is: it collapses beyond about 11 MHz. Picking 14 MHz at midnight (a common daytime choice) yields essentially zero reliability. The noon curve is broader but lower, peaking around 0.70.
Panel (5): 3D optimal frequency versus time and solar activity
This surface shows how the entire day-night pattern scales with the solar cycle. At solar minimum ($R_{12}=10$) the optimum is about 5.2 MHz at midnight and 12.0 MHz at noon. At $R_{12}=100$ it is 8.3 MHz and 22.1 MHz, and at $R_{12}=150$ it reaches 10.0 MHz and 27.5 MHz. The daytime plateau grows much faster than the nighttime floor, which reflects the stronger solar dependence of the daytime F2 layer. In practice this means that the same station needs a very different band plan at solar minimum than at solar maximum, with the higher bands opening only when the sunspot number is high.
Panel (6): Achievable performance at the optimum
The green line is the best reliability achievable at each hour. It is highest, at about 0.88, just before sunrise and after sunset, and lowest, at about 0.70, at noon. The red dashed line shows the SNR at the chosen frequency and follows the same pattern, staying well above the 10 dB requirement throughout. The dip at midday might look surprising, since the ionosphere is at its “best” then. The cause is that the optimum frequency is forced to be high (22 MHz), which costs free-space loss, and the residual absorption remains larger in daylight. The transitions at dawn and dusk give the best compromise: the MUF is high enough to allow a comfortable frequency, but the D layer is only weakly ionized.
Conclusion
We turned the qualitative rule of thumb “use a frequency somewhat below the MUF” into a quantitative optimization. By combining a geometric model, a simple ionosphere model, a link budget and a probabilistic reliability function, we found an optimal frequency for each hour of the day. The optimum lands naturally at 75 to 85 percent of the MUF, and it moves by more than a factor of two between night and day. Because every function was written to broadcast over NumPy arrays, adding another dimension such as solar activity cost nothing in code complexity, and a full three-dimensional sweep ran in a fraction of a second.
The same framework extends easily. You can add multi-hop paths, replace the toy ionosphere with real foF2 predictions, sweep the path distance, or compare antenna designs by changing GT_DBI and GR_DBI. Since the objective function is a plain Python function, any of these changes only requires editing the model in Section 3.
