Chapter 26 — Meteor Scatter & Passive Radar (optional)¶
!!! info "Before you start" Prerequisites: Ch 5 (Hands-on SDR) · Maths Lab: Lab B (Matched Filtering) · ~40 min · Intermediate
This chapter is OPTIONAL and hardware-flavoured, in the spirit of Chapter 5 (Hands-on SDR) and Chapter 6 (Hydrogen Line). It points at one of the cheapest and friendliest real experiments in amateur radio astronomy: counting meteors by radio. You can do it for real with a $30 RTL-SDR and a scrap of wire — but you need no hardware and no network to work through this chapter. Every cell runs end-to-end on the plain
janskybase environment using the simulation injansky.meteor. The one real-hardware recipe is clearly marked and never required.
When a meteoroid the size of a sand grain hits the upper atmosphere at tens of kilometres per second, it ablates: it vaporises and ionises the air along its path, leaving a thin column of free electrons ~80–120 km up. That column lasts a fraction of a second to a few seconds — and while it lasts, it reflects radio waves. Point a VHF receiver at a powerful transmitter that is below your horizon (so you normally hear nothing), and each meteor trail briefly bounces that transmitter into your antenna: a short "ping." This is forward meteor scatter, and counting those pings over a night is a genuine, contributable measurement — the meteor rate, which amateurs report to networks like the IMO and RMOB.
What you'll learn¶
- What forward meteor scatter is, and the two cheap ways amateurs do it: the GRAVES radar at 143.050 MHz (Europe) and reflecting a distant FM carrier (worldwide).
- The two classic echo shapes — underdense (instant rise, exponential decay) and overdense (longer plateau) — and the physics that sets them.
- How to read pings in a power time series and a spectrogram, and how to count them with a threshold (the sensitivity-vs-false-alarm trade-off).
- Why the meteor rate over a night or a shower is the real scientific product.
- The neighbouring idea of passive (bistatic) radar — using a transmitter you
don't control to detect aircraft or meteors — and where the open-source
blah2receiver fits in.
This chapter connects directly to the SDR machinery of
Chapter 5 (the spectrogram is the same scipy.signal
tool) and to the hardware field guide in
docs/projects.md and
docs/field-notes.md.
The setup (and the real-hardware path)¶
Forward scatter is bistatic: the transmitter (Tx) and your receiver (Rx) are in different places, hundreds to ~2000 km apart, and a meteor trail at the midpoint acts like a momentary mirror.
meteor trail (ionised column, ~100 km up)
\ * /
Tx ───╳────── Rx Tx is BELOW Rx's horizon:
(GRAVES / (you) no direct signal, only meteor "pings"
FM tower)
Real-hardware recipe (you do NOT need this to run the chapter). Tune an RTL-SDR in narrow-band mode to a quiet frequency just off a strong distant carrier, point a small directional VHF antenna (a 3-element Yagi is plenty) up at the sky, and record the audio/power. Two proven targets:
- GRAVES (Europe), 143.050 MHz. A French space-surveillance radar that floodlights the sky. Tune ~143.049 MHz (a slight offset so a ping appears as an audible tone) and watch for pings. See RTL-SDR.com: GRAVES meteors.
- A distant FM carrier (worldwide). No GRAVES nearby? Tune an unused ~88 MHz FM frequency whose nearest user is hundreds of km away and point the antenna up; meteors reflect that far transmitter into your receiver. See BAA: Meteor Reflections.
It is always legal to listen — radio astronomy is receive-only (Field Notes). The catch from Chapter 5 applies: the RTL-SDR's 8-bit ADC has tiny dynamic range, so a band-pass filter and pointing away from strong local FM help a lot.
Everything below uses the jansky.meteor simulation, so the whole chapter is
reproducible offline. We write the single-ping echo shapes ourselves in Section 2, then
reach for the packaged jansky.meteor copies -- simulate_meteor_timeseries,
detect_pings and its MeteorDetection result, and eventually underdense_echo /
overdense_echo themselves -- for the bulk, repetitive work of building and scanning a
whole night of pings.
Back of the envelope¶
Before any code: the diagram above already has enough geometry in it to answer a real question about why meteor trails have to be so high up. Forward-scatter links (GRAVES, distant FM carriers) routinely span hundreds to ~2000 km between transmitter and receiver, as noted above -- and the ionised trail, sitting near the midpoint, has to be tall enough that both Tx and Rx can see it over the curve of the Earth. Meteors ablate at roughly $h \approx 100$ km altitude.
The question: for a trail at $h = 100$ km, what is the longest Tx-Rx ground separation $L_\mathrm{max}$ that simple horizon geometry allows?
Rules of the napkin: decades matter, factors of two don't. Treat the trail as a mirror sitting at the midpoint, and treat "visible" as "above the geometric horizon" -- the classic horizon distance from a ground station to a point of height $h$ is $d = \sqrt{2 R_\oplus h}$, and you need two such sightlines (one from Tx, one from Rx).
R_earth = 6371.0 # km -- Earth's mean radius
h_trail = 100.0 # km -- typical meteor ablation height (see the intro above)
# Your one line: horizon distance from one ground station to a point at height h_trail is
# d = sqrt(2 * R_earth * h_trail); Tx and Rx each need their own sightline to the trail, so
# the max separation is twice that.
L_max_guess_km = None # <- replace None with an expression in R_earth and h_trail
Grade it with the course's envelope checker, jansky.envelope.check, introduced in
Chapter 1: within half a decade is envelope-grade,
within one decade is the right order of magnitude, further out usually means a units slip.
from jansky.envelope import check
# The answer is worked in the reveal below -- commit to your own number first.
# The check stores only its base-10 logarithm, so a stray glance can't spoil it.
check(L_max_guess_km, expected_log10=3.3537, name="max Tx-Rx separation via a 100 km trail", units="km")
[max Tx-Rx separation via a 100 km trail] no guess yet — fill in the cell above, then re-run this one.
False
Reveal — the envelope answer
One ground station's horizon distance to a point at height $h$:
$$ d = \sqrt{2 R_\oplus h} = \sqrt{2 \times 6371~\mathrm{km} \times 100~\mathrm{km}} \approx 1129~\mathrm{km}. $$
Both Tx and Rx need their own sightline to the trail, so the longest usable separation is about twice that:
$$ L_\mathrm{max} \;\approx\; 2d \;\approx\; 2258~\mathrm{km}. $$
import math
d = math.sqrt(2 * R_earth * h_trail)
L_max_guess_km = 2 * d # ~2258 km
That lands right on the chapter's own numbers -- "hundreds to ~2000 km apart" for real GRAVES/FM links, and the ~600 km baseline used later in Section 5.
Where the envelope leaks (and where this chapter patches it):
- We used a flat, unrefracted horizon. Real radio waves bend gently toward the Earth (roughly the "4/3 Earth radius" rule), stretching the true horizon a few percent further -- inside our factor-of-two tolerance, but it's why real link budgets use an effective Earth radius rather than the geographic one.
- We put the trail exactly at the geometric midpoint and ignored antenna gain and pointing. Section 5's bistatic-range ellipse is the fuller picture: a single delay measurement only constrains you to an ellipse, not a point.
- We said nothing about whether the trail reflects strongly enough to be heard at that range -- that depends on whether it's under- or overdense and how fast it decays (Section 2), and on the threshold you set to call it a detection at all (Section 4).
- Real trails ablate anywhere from ~80-120 km, not a fixed 100 km -- turn the knob below to see how much that range matters.
Turn the knob¶
Invert the question: how tall would a trail have to be to close a hypothetical transatlantic, 4000 km forward-scatter hop? Reuse your expression, solved for $h$ instead of $L$.
L_hop_km = 4000.0 # km -- a hypothetical transatlantic forward-scatter path
# Solve d = sqrt(2 * R_earth * h) for h, with d = L_hop_km / 2.
h_needed_guess_km = None # <- your line, solved for h
check(h_needed_guess_km, expected_log10=2.4969, name="trail height for a 4000 km hop", units="km")
[trail height for a 4000 km hop] no guess yet — fill in the cell above, then re-run this one.
False
Reveal — the scaling answer
$$ h \;=\; \frac{(L/2)^2}{2 R_\oplus} \;=\; \frac{2000^2}{2\times6371}~\mathrm{km} \;\approx\; 314~\mathrm{km}. $$
h_needed_guess_km = (L_hop_km / 2) ** 2 / (2 * R_earth) # ~314 km
Three hundred kilometres is well above where meteors actually ablate (~80-120 km, Section 2) -- that high up there's essentially no air left to ionise a trail in the first place. That's the honest geometric reason forward meteor scatter tops out around ~2000-2500 km: past that, no physically achievable trail height clears the horizon for both ends of the link, no matter how good your receiver is.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal as sig
from jansky import signals, plotting, meteor
plotting.use_jansky_style()
# Seed everything for reproducible figures. 143 -> the GRAVES band, in MHz.
SEED = 143
rng = signals.rng(SEED)
# The meteor helpers this chapter reaches for once the physics is in the open --
# for the bulk/repetitive work of building and scanning a whole night of pings.
print("jansky.meteor helpers used for the bulk work below:")
for name in (
"simulate_meteor_timeseries",
"detect_pings",
"MeteorDetection",
):
print(f" meteor.{name}")
# Is real SDR capture available? (Purely informational -- nothing below needs it.)
try:
from rtlsdr import RtlSdr # provided by `pyrtlsdr`, needs a real dongle
HAVE_RTLSDR = True
except Exception as exc:
RtlSdr = None
HAVE_RTLSDR = False
_why = exc
if HAVE_RTLSDR:
print("\npyrtlsdr importable -- real GRAVES/FM capture *may* be possible.")
print("(We still use the simulation so results are reproducible.)")
else:
print(f"\npyrtlsdr not available -> using the SIMULATION (base env, no hardware).")
print(f" reason: {type(_why).__name__}: {_why}")
jansky.meteor helpers used for the bulk work below: meteor.simulate_meteor_timeseries meteor.detect_pings meteor.MeteorDetection pyrtlsdr not available -> using the SIMULATION (base env, no hardware). reason: ImportError: Error loading librtlsdr. Make sure librtlsdr (and all of its dependencies) are in your path
2. Echo shapes — underdense vs overdense¶
What a single ping looks like depends on how much ionisation the meteor laid down. The dividing line is the trail's electron line density $q$ (electrons per metre); the natural break is near $q \approx 2\times10^{14}\ \mathrm{e^-/m}$.
Underdense trails — the exponential decay¶
In a faint (underdense) trail the electrons are sparse enough that the radio wave passes through and each electron scatters more or less independently. The echo rises essentially instantly when the trail forms, then decays exponentially as the column spreads out and the free electrons recombine and diffuse away. The decay is set by ambipolar diffusion — electrons and ions drift apart, smearing the sharp column until it no longer reflects coherently. The received amplitude follows
$$ A(t) \;=\; A_0\,\exp\!\left(-\frac{t-t_0}{\tau}\right), \qquad \tau \;=\; \frac{\lambda^2}{16\pi^2 D}, $$
where $\lambda$ is the radio wavelength and $D$ is the ambipolar diffusion coefficient (which grows with altitude as the air thins). At VHF, $\tau$ ranges from milliseconds to about a second. Let's write that formula out ourselves below.
Overdense trails — the plateau¶
In a dense (overdense) trail there are so many electrons that the column behaves like a solid metallic cylinder: the wave can't penetrate, it reflects off the surface, and it keeps reflecting until the expanding column's critical density drops below the radio frequency. These echoes last longer (up to many seconds for bright fireballs) and have a ragged, roughly flat-topped (plateau) envelope rather than a clean exponential. Overdense trails are rarer (bright meteors are rare), so in a typical record most pings are underdense.
Let's write both formulas out and plot them.
# A fine time axis for one ping, in seconds.
t = np.linspace(0, 3.0, 3000)
t0 = 0.3 # s -- the instant the trail forms
amplitude = 10.0 # arbitrary echo-amplitude units
# Underdense echo: instantaneous rise at t0, then exponential decay
# A(t) = A0 * exp(-(t - t0) / tau)
# -- the formula derived above, written out in flat NumPy (a few decay times shown).
after = t >= t0
u_fast = np.zeros_like(t)
u_fast[after] = amplitude * np.exp(-(t[after] - t0) / 0.10) # tau = 0.10 s (high, thin air)
u_slow = np.zeros_like(t)
u_slow[after] = amplitude * np.exp(-(t[after] - t0) / 0.40) # tau = 0.40 s (lower, denser air)
# Overdense echo: a longer, ragged plateau of a given duration, smooth sin^0.5 rise/fall.
duration = 1.8 # s
window = (t >= t0) & (t <= t0 + duration)
phase = (t[window] - t0) / duration
o_echo = np.zeros_like(t)
o_echo[window] = amplitude * np.sin(np.pi * phase) ** 0.5
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4.5), sharey=True)
ax1.plot(t, u_fast, color="#1f77b4", lw=1.8, label=r"$\tau = 0.10$ s (high, thin air)")
ax1.plot(t, u_slow, color="#2ca02c", lw=1.8, label=r"$\tau = 0.40$ s (lower, denser)")
ax1.axvline(0.3, color="0.5", ls=":", lw=1)
ax1.set_title("Underdense echo: instant rise, exp. decay")
ax1.set_xlabel("time [s]")
ax1.set_ylabel("received amplitude")
ax1.legend(fontsize=9)
ax2.plot(t, o_echo, color="#d62728", lw=1.8, label="duration = 1.8 s")
ax2.axvline(0.3, color="0.5", ls=":", lw=1)
ax2.set_title("Overdense echo: longer ragged plateau")
ax2.set_xlabel("time [s]")
ax2.legend(fontsize=9)
fig.tight_layout()
plt.show()
The left panel shows two underdense echoes with different decay times $\tau$.
Both jump up the instant the trail forms, then fall off as $e^{-(t-t_0)/\tau}$.
Higher trails sit in thinner air where the diffusion coefficient $D$ is larger,
so they smear out faster (smaller $\tau$, blue); lower trails linger (green).
The right panel shows an overdense echo: it switches on, holds a rough plateau
for its full duration, then switches off — a qualitatively different shape that
betrays a much brighter meteor. In the next section we let
simulate_meteor_timeseries sprinkle a random mix of these into a noisy record.
3. A night of meteors — the power record¶
A real meteor monitor logs received power vs time, all night long. Meteors
arrive at random, so the number in any interval is Poisson-distributed about
a mean rate; the background between pings is receiver/sky noise (the radiometer
noise of Chapter 3).
We just wrote the underdense and overdense shapes by hand for a single ping; a whole
night means stamping down many of them at random (Poisson) times, mixed 80%/20%
underdense/overdense, on top of noise. Retyping that placement loop every time would just
invite bugs, so here we reach for meteor.simulate_meteor_timeseries -- the packaged
version of exactly the two formulas above, scaled up with random timing and a mean
rate per minute. It returns the time axis, the power, and the true ping times (our
ground truth for the counting section).
We simulate a 2-minute stretch at a brisk ~30 meteors/minute so several pings are visible at once.
# Simulate a 2-minute record. The helper seeds its own RNG for reproducibility.
DURATION_S = 120.0 # seconds (2 minutes)
RATE_PER_MIN = 30.0 # mean meteors per minute (Poisson)
SAMPLE_RATE = 100.0 # samples per second
NOISE = 1.0 # background noise sigma
times, power, true_pings = meteor.simulate_meteor_timeseries(
duration_s=DURATION_S,
rate_per_min=RATE_PER_MIN,
sample_rate=SAMPLE_RATE,
noise=NOISE,
seed=SEED,
)
print(f"record length : {DURATION_S:.0f} s ({times.size:,} samples @ {SAMPLE_RATE:.0f} Hz)")
print(f"mean rate set : {RATE_PER_MIN:.0f} meteors/min")
print(f"TRUE pings : {true_pings.size} (this realisation; Poisson scatter)")
print(f"expected mean : {RATE_PER_MIN * DURATION_S / 60:.0f} (rate x minutes)")
record length : 120 s (12,000 samples @ 100 Hz) mean rate set : 30 meteors/min TRUE pings : 58 (this realisation; Poisson scatter) expected mean : 60 (rate x minutes)
fig, ax = plt.subplots(figsize=(12, 4.5))
ax.plot(times, power, color="#1f77b4", lw=0.6)
# Mark the TRUE ping arrival times along the top.
ax.plot(
true_pings,
np.full_like(true_pings, power.max() * 1.05),
"v",
color="#d62728",
ms=6,
label="true ping times",
)
ax.set_xlabel("time [s]")
ax.set_ylabel("received power")
ax.set_title(f"A night of meteors: power record ({RATE_PER_MIN:.0f}/min, {DURATION_S:.0f} s)")
ax.legend(loc="upper right", fontsize=9)
fig.tight_layout()
plt.show()
Each red marker is a real meteor. Most pings are the underdense spikes — a sharp jump followed by a fast decay back into the noise. A handful are taller and wider: those are the overdense plateaus from brighter meteors. The flat band near the bottom is the noise floor; a ping has to poke clearly above it to be counted, which is the whole problem of the next section.
The spectrogram view¶
A real receiver records not just power but a band, so the natural display is the
same waterfall / spectrogram from Chapter 5: time on
one axis, frequency on the other, brightness = power. To show what an operator
actually sees, we promote our power record to a narrow-band signal: model each
ping as the off-tuned carrier tone (a meteor reflecting GRAVES appears as a
short tone a few hundred Hz from your tuned frequency), amplitude-modulated by the
echo envelope, in complex noise. Then we hand it to scipy.signal.spectrogram.
# Promote the power envelope to a narrow-band IQ signal so we can make a
# spectrogram exactly like Chapter 5. Each ping shows up as a short tone offset
# from the tuned frequency (here +350 Hz, a typical GRAVES audio offset).
fs_audio = SAMPLE_RATE
t_audio = times
carrier_offset = 35.0 # Hz offset of the reflected carrier from tuned freq
# Envelope = the (clipped) echo part of the power record above the noise.
envelope = np.clip(power, 0, None)
iq_audio = envelope * np.exp(1j * 2 * np.pi * carrier_offset * t_audio)
# Add complex baseband noise.
g = signals.rng(SEED)
iq_audio += (g.normal(0, NOISE, t_audio.size) + 1j * g.normal(0, NOISE, t_audio.size)) / np.sqrt(2)
f_spec, t_spec, Sxx = sig.spectrogram(
iq_audio,
fs=fs_audio,
window="hann",
nperseg=128,
noverlap=96,
return_onesided=False,
scaling="density",
)
f_spec = np.fft.fftshift(f_spec)
Sxx = np.fft.fftshift(Sxx, axes=0)
Sxx_db = 10 * np.log10(Sxx + 1e-9)
fig, ax = plt.subplots(figsize=(12, 4.5))
extent = [t_spec[0], t_spec[-1], f_spec[0], f_spec[-1]]
im = ax.imshow(Sxx_db, aspect="auto", extent=extent, origin="lower", cmap="inferno")
ax.axhline(carrier_offset, color="cyan", ls=":", lw=1, alpha=0.6, label="reflected carrier offset")
ax.set_xlabel("time [s]")
ax.set_ylabel("frequency offset [Hz]")
ax.set_title("Spectrogram (waterfall): each ping is a brief tone at the carrier offset")
ax.legend(loc="upper right", fontsize=9)
fig.colorbar(im, ax=ax, label="power [dB]", fraction=0.030, pad=0.02)
fig.tight_layout()
plt.show()
Read it like the Chapter 5 waterfall: the meteor pings are the bright vertical dashes sitting on the dotted carrier-offset line — short in time, narrow in frequency. This is precisely what GRAVES monitoring software (e.g. Spectrum Lab, HROFFT) shows; operators (and automatic counters) tally these streaks. The power-vs-time plot and this waterfall are two views of the same events. Now we count them automatically.
4. Counting them — detection and the threshold trade-off¶
To turn a wiggly power record into a number, we threshold it: a ping is any
excursion that rises more than $k$ noise-sigmas above the baseline. We estimate the noise
level with the median absolute deviation (MAD) (so a few huge pings don't inflate the
noise estimate), find samples above threshold $\times \sigma$, and merge samples within
min_separation_s into a single event at its peak -- three named steps, written out below
in flat NumPy.
The only knob that matters is threshold. Set it too low and noise spikes get
counted as meteors (false alarms); set it too high and you miss the faint
underdense pings (low sensitivity). Let's run it at a sensible value and
compare to the truth.
THRESHOLD = 6.0 # in units of noise sigma (MAD estimate)
MIN_SEP_S = 0.3 # merge crossings within this many seconds into one event
# 1. A robust noise sigma from the median absolute deviation (MAD) -- a few huge
# pings can't inflate this the way they would a plain standard deviation.
med = np.median(power)
mad = np.median(np.abs(power - med))
sigma = 1.4826 * mad # MAD -> Gaussian-equivalent sigma
# 2. Flag every sample that rises more than THRESHOLD sigma above the baseline.
above = np.where(power - med > THRESHOLD * sigma)[0]
# 3. Merge consecutive/nearby crossings into one event per meteor, keeping the
# time of each event's power maximum.
groups = np.split(above, np.where(np.diff(times[above]) > MIN_SEP_S)[0] + 1)
det_times = np.array([times[g[np.argmax(power[g])]] for g in groups])
det_count = det_times.size
print(f"threshold : {THRESHOLD:.1f} sigma")
print(f"DETECTED pings: {det_count}")
print(f"TRUE pings : {true_pings.size}")
print(f"difference : {det_count - true_pings.size:+d}")
# Match detections to truth to estimate completeness & false-alarm rate.
TOL = 0.5 # s: a detection within TOL of a true ping counts as a real match
matched_true = 0
for tp in true_pings:
if det_times.size and np.min(np.abs(det_times - tp)) < TOL:
matched_true += 1
false_alarms = 0
for dt in det_times:
if true_pings.size == 0 or np.min(np.abs(true_pings - dt)) >= TOL:
false_alarms += 1
completeness = matched_true / true_pings.size if true_pings.size else 0.0
print(
f"\ncompleteness : {completeness:5.1%} ({matched_true}/{true_pings.size} true pings recovered)"
)
print(f"false alarms : {false_alarms} (detections with no true ping nearby)")
threshold : 6.0 sigma DETECTED pings: 48 TRUE pings : 58 difference : -10 completeness : 89.7% (52/58 true pings recovered) false alarms : 2 (detections with no true ping nearby)
# Show the detections on the power record.
fig, ax = plt.subplots(figsize=(12, 4.5))
ax.plot(times, power, color="#1f77b4", lw=0.6, zorder=1)
# The robust threshold line that detect_pings effectively uses (MAD-based).
med = np.median(power)
mad = np.median(np.abs(power - med))
sigma = 1.4826 * mad
ax.axhline(
med + THRESHOLD * sigma,
color="0.4",
ls="--",
lw=1.2,
label=f"threshold = {THRESHOLD:.0f}$\\sigma$",
)
ax.plot(
true_pings,
np.full_like(true_pings, power.max() * 1.10),
"v",
color="#d62728",
ms=6,
label=f"true ({true_pings.size})",
)
ax.plot(
det_times,
np.full_like(det_times, power.max() * 0.97),
"^",
color="#2ca02c",
ms=6,
label=f"detected ({det_count})",
)
ax.set_xlabel("time [s]")
ax.set_ylabel("received power")
ax.set_title("Detected pings vs truth")
ax.legend(loc="upper right", fontsize=9)
fig.tight_layout()
plt.show()
Sweeping the threshold¶
We just wrote the MAD-threshold-and-merge logic once, by hand. Sweeping it over nineteen
threshold values means running that same three-step logic nineteen times, so here we reach
for meteor.detect_pings -- the packaged, unit-tested copy of exactly the cell above,
returning a MeteorDetection(times, count) -- rather than retyping it inside a loop. The
trade-off becomes obvious once we sweep threshold and watch how the detected count, the
completeness, and the false alarms move. We expect: low thresholds
over-count (false alarms from noise), high thresholds under-count (missed
faint pings), and a broad sweet spot in between where the detected count tracks
the truth.
thresholds = np.arange(3.0, 12.1, 0.5)
counts, comps, fas = [], [], []
for th in thresholds:
d = meteor.detect_pings(times, power, threshold=th, min_separation_s=0.3)
counts.append(d.count)
mt = sum(1 for tp in true_pings if d.times.size and np.min(np.abs(d.times - tp)) < TOL)
fa = sum(1 for dt in d.times if true_pings.size == 0 or np.min(np.abs(true_pings - dt)) >= TOL)
comps.append(mt / true_pings.size if true_pings.size else 0.0)
fas.append(fa)
counts, comps, fas = map(np.asarray, (counts, comps, fas))
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4.5))
ax1.plot(thresholds, counts, "o-", color="#1f77b4", label="detected count")
ax1.axhline(
true_pings.size, color="#d62728", ls="--", lw=1.5, label=f"true count = {true_pings.size}"
)
ax1.axvline(THRESHOLD, color="0.5", ls=":", lw=1)
ax1.set_xlabel(r"threshold [$\sigma$]")
ax1.set_ylabel("number of pings")
ax1.set_title("Detected count vs threshold")
ax1.legend(fontsize=9)
ax2.plot(thresholds, comps, "o-", color="#2ca02c", label="completeness")
ax2.plot(thresholds, fas / max(fas.max(), 1), "s--", color="#ff7f0e", label="false alarms (scaled)")
ax2.axvline(THRESHOLD, color="0.5", ls=":", lw=1)
ax2.set_xlabel(r"threshold [$\sigma$]")
ax2.set_ylabel("fraction")
ax2.set_title("Sensitivity vs false alarms")
ax2.legend(fontsize=9)
fig.tight_layout()
plt.show()
The picture is the classic detection trade-off. At low thresholds (left of each plot) the detected count balloons above the true count: the detector is firing on noise, so false alarms (orange) are high even though completeness (green) is near 100%. At high thresholds the noise false alarms vanish but completeness falls — the faintest underdense pings never clear the bar, so the count drops below truth. In between, around 6$\sigma$, the detected count sits closest to the true count with few false alarms. There is no universally "correct" threshold; you pick it for your noise level and how much you fear false positives versus missed meteors — and then keep it fixed so your rates are comparable night to night.
The real product: meteor rate¶
A single count isn't the science — the rate over time is. Amateurs bin pings per hour all through a night (and especially across a shower like the Perseids or Geminids) and the rate curve, corrected for sensitivity, is what gets contributed to networks like the IMO Video Meteor Network and RMOB. Radio has a superpower here that optical lacks: it works in daylight and through cloud, so radio is how the daytime meteor showers (the Arietids, the $\zeta$-Perseids) are monitored at all. Let's turn our detections into a rate.
# Bin the detected pings into a rate (per minute) over the record.
bin_s = 20.0
edges = np.arange(0, DURATION_S + bin_s, bin_s)
det_hist, _ = np.histogram(det_times, bins=edges)
true_hist, _ = np.histogram(true_pings, bins=edges)
centers = 0.5 * (edges[:-1] + edges[1:])
to_per_min = 60.0 / bin_s
fig, ax = plt.subplots(figsize=(11, 4.5))
ax.step(centers, true_hist * to_per_min, where="mid", color="#d62728", lw=1.8, label="true rate")
ax.step(
centers,
det_hist * to_per_min,
where="mid",
color="#2ca02c",
lw=1.8,
ls="--",
label="detected rate",
)
ax.axhline(RATE_PER_MIN, color="0.5", ls=":", lw=1.2, label=f"mean set = {RATE_PER_MIN:.0f}/min")
ax.set_xlabel("time [s]")
ax.set_ylabel("meteor rate [per minute]")
ax.set_title("Meteor rate over the record (the contributable product)")
ax.legend(fontsize=9)
fig.tight_layout()
plt.show()
The detected rate (green) tracks the true rate (red) and both scatter around the mean we set (dotted) — that scatter is just Poisson counting statistics, the fundamental noise of any rate measurement. Over a real night you would see this curve climb during a shower's peak; the shape of that climb, the zenithal hourly rate, is the headline number meteor astronomers care about.
5. Passive radar — a brief detour¶
Forward meteor scatter is one member of a bigger family: passive (bistatic) radar, where you exploit a powerful transmitter you don't own or control — an FM tower, a digital-TV transmitter, GRAVES — as your illuminator, and only ever receive. A conventional radar transmits a pulse and times the echo; passive radar instead listens with two channels:
- a reference channel, an antenna pointed at the illuminator to record a clean copy of the transmitted signal;
- a surveillance channel, an antenna pointed at the sky/airspace to catch the faint reflections off targets (aircraft, the Moon, ionised meteor trails).
You then cross-correlate the two channels. The correlation peaks at the time
delay of the echo (which gives the bistatic range — the extra path length via
the target) and at its Doppler shift (which gives the target's line-of-sight
velocity). Sweep delay vs Doppler and you get the range–Doppler map that
reveals each moving target as a blob. This needs two coherent receivers — the
same phase-coherence problem as amateur interferometry in
Chapter 17 and
Field Notes — which is why people use a dual-channel
dongle or a shared clock. The open-source blah2
receiver (see docs/field-notes.md) implements exactly
this with two RTL-SDRs and an FM station, and will draw you a live range–Doppler
map of the aircraft overhead.
Here is the tiny piece of geometry that turns a measured time delay into a
bistatic range — no hardware, just arithmetic with astropy constants.
from astropy import constants as const, units as u
# Bistatic geometry: signal travels Tx -> target -> Rx. The "bistatic range" is
# that total path length; the EXTRA path over the direct Tx->Rx line is what the
# cross-correlation delay measures.
tx_rx_baseline = 600.0 * u.km # distance from transmitter to receiver
measured_delay = 1.7 * u.ms # extra delay of the echo vs the direct path
c = const.c
extra_path = (c * measured_delay).to(u.km) # bistatic range SUM minus baseline
bistatic_range_sum = (tx_rx_baseline + extra_path).to(u.km) # R_tx + R_rx
print(f"Tx-Rx baseline : {tx_rx_baseline:.0f}")
print(f"measured echo delay : {measured_delay:.2f}")
print(f"extra path (c x delay): {extra_path:.0f}")
print(f"bistatic range sum : {bistatic_range_sum:.0f} (= R_tx + R_rx)")
print("\nThe target lies on an ELLIPSE with the Tx and Rx at its foci")
print(f"(the locus of points whose path sum = {bistatic_range_sum:.0f}).")
print("A second receiver, or the target's Doppler, pins it down further.")
Tx-Rx baseline : 600 km measured echo delay : 1.70 ms extra path (c x delay): 510 km bistatic range sum : 1110 km (= R_tx + R_rx) The target lies on an ELLIPSE with the Tx and Rx at its foci (the locus of points whose path sum = 1110 km). A second receiver, or the target's Doppler, pins it down further.
That ellipse is the heart of bistatic radar: a single delay measurement does
not give a point, it gives the iso-range ellipse of all locations whose
total path $R_\mathrm{tx}+R_\mathrm{rx}$ equals the measured value, with the
transmitter and receiver at the two foci. Add the Doppler shift (constrains
the velocity) and a second receiver (a second ellipse — they intersect at the
target) and you localise the aircraft or meteor. blah2 does all of this in real
time; this one-liner is just the first rung. For meteors specifically, the same
correlation lets you measure a trail's range and radial speed, turning the
ping-counter of this chapter into a proper meteor radar.
From napkin to package¶
Both echo shapes were typed out in the open in Section 2, because their functional form
is this chapter's physics. But everything downstream -- the night-long record, the shower
burst in Exercise 2, the decay-time fit in Exercise 3 -- needs many of these echoes stamped
down at once, so the course keeps one tested copy of each formula in jansky.meteor, and
calls it by name from here on. Two asserts prove the packaged versions compute exactly the
lines you wrote above:
t_check = np.linspace(0, 3.0, 500)
# Underdense: instantaneous rise then exponential decay -- the Section 2 formula.
after_check = t_check >= 0.3
inline_under = np.zeros_like(t_check)
inline_under[after_check] = 10.0 * np.exp(-(t_check[after_check] - 0.3) / 0.10)
packaged_under = meteor.underdense_echo(t_check, t0=0.3, amplitude=10.0, decay_time=0.10)
# Overdense: the ragged plateau -- the Section 2 formula.
window_check = (t_check >= 0.3) & (t_check <= 0.3 + 1.8)
phase_check = (t_check[window_check] - 0.3) / 1.8
inline_over = np.zeros_like(t_check)
inline_over[window_check] = 10.0 * np.sin(np.pi * phase_check) ** 0.5
packaged_over = meteor.overdense_echo(t_check, t0=0.3, amplitude=10.0, duration=1.8)
assert np.allclose(inline_under, packaged_under)
assert np.allclose(inline_over, packaged_over)
print("inline formulas == jansky.meteor.underdense_echo / overdense_echo -- promoted.")
inline formulas == jansky.meteor.underdense_echo / overdense_echo -- promoted.
Try it yourself¶
Exercise 1 — Change the rate, compare detected vs true¶
simulate_meteor_timeseries takes rate_per_min. Re-run it at a low rate
(say 5/min) and a high rate (say 60/min), detect with the same threshold, and
compare detected vs true counts. At high rates pings start to overlap within
min_separation_s and get merged — does the detector under-count? A working
starter is filled in.
# Exercise 1: sweep the meteor rate and compare detected vs true counts.
for rate in (5.0, 15.0, 60.0):
tms, pwr, tp = meteor.simulate_meteor_timeseries(
duration_s=120.0,
rate_per_min=rate,
sample_rate=100.0,
noise=1.0,
seed=SEED,
)
d = meteor.detect_pings(tms, pwr, threshold=6.0, min_separation_s=0.3)
print(
f"rate {rate:5.0f}/min -> true {tp.size:3d} detected {d.count:3d} "
f"(diff {d.count - tp.size:+d})"
)
# TODO: at 60/min do detections fall *below* truth as pings overlap and merge?
# Try shrinking min_separation_s and see whether the count recovers.
rate 5/min -> true 9 detected 9 (diff +0) rate 15/min -> true 29 detected 25 (diff -4) rate 60/min -> true 124 detected 33 (diff -91)
Solution
The starter already runs the sweep; the # TODO asks whether shrinking min_separation_s rescues the count at high rates. It helps a little but cannot close the gap: once pings arrive faster than they decay they physically overlap, so many genuine meteors merge into one detected event no matter how small the merge window is. The detector under-counts because the record is saturated, not because the threshold is wrong.
tms, pwr, tp = meteor.simulate_meteor_timeseries(
duration_s=120.0, rate_per_min=60.0, sample_rate=100.0, noise=1.0, seed=SEED,
)
print(f"true pings at 60/min: {tp.size}")
for ms in (0.3, 0.15, 0.05):
d = meteor.detect_pings(tms, pwr, threshold=6.0, min_separation_s=ms)
print(f" min_separation_s={ms:.2f} -> detected {d.count}")
Expected: ~124 true pings, but only ~33 detected at min_separation_s=0.3, rising to ~37 (0.15) and ~45 (0.05) — recovering some, yet still far below truth. Takeaway: at high meteor rates the count flattens out (saturates) because overlapping trails blend together; the honest fix is a faster sample rate / shorter decay times or correcting the rate for this dead-time, not just a smaller merge window.
Exercise 2 — Simulate a shower (a burst of pings)¶
A meteor shower is a burst of activity. Build a record where the rate is low
in the first and last thirds but high in the middle, and watch the rate curve
spike. (Hint: concatenate three simulate_meteor_timeseries calls with different
rate_per_min, offsetting the time axes and ping times. A working version is
provided.)
# Exercise 2: a shower = a burst of high rate sandwiched by quiet background.
seg = 60.0 # seconds per segment
rates = (8.0, 50.0, 8.0) # quiet / SHOWER PEAK / quiet
all_times, all_power, all_pings, offset = [], [], [], 0.0
for i, r in enumerate(rates):
tms, pwr, tp = meteor.simulate_meteor_timeseries(
duration_s=seg,
rate_per_min=r,
sample_rate=100.0,
noise=1.0,
seed=SEED + i,
)
all_times.append(tms + offset)
all_power.append(pwr)
all_pings.append(tp + offset)
offset += seg
times_s = np.concatenate(all_times)
power_s = np.concatenate(all_power)
pings_s = np.concatenate(all_pings)
det_s = meteor.detect_pings(times_s, power_s, threshold=6.0, min_separation_s=0.3)
edges = np.arange(0, offset + 10, 10.0)
hist, _ = np.histogram(det_s.times, bins=edges)
ctr = 0.5 * (edges[:-1] + edges[1:])
fig, ax = plt.subplots(figsize=(11, 4))
ax.step(ctr, hist * (60.0 / 10.0), where="mid", color="#2ca02c", lw=1.8)
ax.axvspan(seg, 2 * seg, color="#d62728", alpha=0.10, label="shower peak window")
ax.set_xlabel("time [s]")
ax.set_ylabel("detected rate [per minute]")
ax.set_title("Exercise 2: a simulated meteor shower (burst of pings)")
ax.legend(fontsize=9)
fig.tight_layout()
plt.show()
print(f"total detected over the 'night': {det_s.count} (true {pings_s.size})")
total detected over the 'night': 40 (true 64)
Solution
The starter is the worked answer: concatenate three simulate_meteor_timeseries segments (quiet / peak / quiet), shifting each segment's time axis and its ping times by a running offset so the three records join seamlessly, then bin the detections into a rate curve. Using a different seed per segment (SEED + i) avoids stitching three identical realisations together.
seg = 60.0
rates = (8.0, 50.0, 8.0) # quiet / SHOWER PEAK / quiet
all_times, all_power, all_pings, offset = [], [], [], 0.0
for i, r in enumerate(rates):
tms, pwr, tp = meteor.simulate_meteor_timeseries(
duration_s=seg, rate_per_min=r, sample_rate=100.0, noise=1.0, seed=SEED + i)
all_times.append(tms + offset); all_power.append(pwr); all_pings.append(tp + offset)
offset += seg
times_s, power_s, pings_s = map(np.concatenate, (all_times, all_power, all_pings))
det_s = meteor.detect_pings(times_s, power_s, threshold=6.0, min_separation_s=0.3)
Expected: the detected-rate step curve is low (~8/min) in the outer thirds and spikes sharply inside the shaded middle minute (~50/min) — exactly the shape of a real shower peak. Takeaway: the contributable product is this rate-vs-time curve, and a shower shows up as a clear burst above the quiet background.
Exercise 3 — Estimate the decay time of an echo¶
The underdense decay time $\tau$ encodes the trail's altitude (via the diffusion
coefficient $D$). Make a single clean underdense echo with underdense_echo, then
recover $\tau$ by fitting a straight line to $\ln A$ vs $t$ on the decaying
part — the slope is $-1/\tau$. Confirm you get back the value you put in.
# Exercise 3: recover tau from an underdense echo by log-linear fitting.
t_e = np.linspace(0, 2.0, 2000)
TRUE_TAU = 0.25
echo = meteor.underdense_echo(t_e, t0=0.2, amplitude=12.0, decay_time=TRUE_TAU)
# Fit only the decaying part, well above the (here zero) floor.
mask = (t_e > 0.2) & (echo > 0.05 * echo.max())
slope, intercept = np.polyfit(t_e[mask] - 0.2, np.log(echo[mask]), 1)
tau_fit = -1.0 / slope
print(f"true tau : {TRUE_TAU:.3f} s")
print(f"recovered tau: {tau_fit:.3f} s")
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.semilogy(t_e, np.clip(echo, 1e-3, None), color="#1f77b4", lw=1.6, label="echo")
ax.semilogy(
t_e[mask],
np.exp(intercept + slope * (t_e[mask] - 0.2)),
color="#d62728",
ls="--",
lw=1.6,
label=rf"fit: $\tau$={tau_fit:.3f}s",
)
ax.set_xlabel("time [s]")
ax.set_ylabel("amplitude (log)")
ax.set_title("Exercise 3: recovering the decay time")
ax.legend(fontsize=9)
fig.tight_layout()
plt.show()
# TODO: add noise to the echo (echo += rng.normal(0, 0.3, echo.size)) and see how
# the recovered tau degrades -- real fits need the noise floor handled.
true tau : 0.250 s recovered tau: 0.250 s
Solution
The clean fit in the starter recovers $\tau$ exactly because $\ln A$ is perfectly linear in $t$ for an underdense echo. The # TODO asks what noise does: add Gaussian noise, and the recovered $\tau$ is biased high, because near the tail the noise floor stops the curve from falling, flattening the apparent slope ($-1/\tau$) — the classic reason real fits must mask out or model the noise floor before taking the log.
rng = signals.rng(SEED)
for noise_sigma in (0.0, 0.1, 0.3):
echo_n = meteor.underdense_echo(t_e, t0=0.2, amplitude=12.0, decay_time=TRUE_TAU)
echo_n = echo_n + rng.normal(0, noise_sigma, echo_n.size)
mask = (t_e > 0.2) & (echo_n > 0.05 * echo_n.max())
slope, _ = np.polyfit(t_e[mask] - 0.2, np.log(np.clip(echo_n[mask], 1e-6, None)), 1)
print(f"noise={noise_sigma:.1f} -> tau_fit={-1.0/slope:.3f} s (true {TRUE_TAU:.3f})")
Expected: tau_fit ≈ 0.250 at zero noise, ~0.255 at $\sigma=0.1$, and ~0.34 at $\sigma=0.3$ — the bias grows with noise and always over-estimates $\tau$. Takeaway: log-linear fitting recovers the decay time cleanly only above the noise floor; with real, noisy pings you must restrict the fit to high-SNR samples (or fit the amplitude directly) to avoid biasing $\tau$ upward.
!!! tip "Contribute & compare — meteor counts" Radio meteor detection works in daylight and through cloud, and amateurs feed real counts to the American Meteor Society. An RTL-SDR on a distant VHF transmitter (or the GRAVES radar) is enough to start; the Society of Amateur Radio Astronomers coordinates amateur meteor work and offers student/teacher grants to fund a build.
Recap¶
- Forward meteor scatter reflects a below-horizon transmitter off a meteor's ionised trail, producing a short ping in your receiver. It's the cheapest real radio astronomy there is — an RTL-SDR plus a wire, aimed off GRAVES (143.050 MHz) in Europe or a distant FM carrier worldwide.
- Pings come in two flavours: underdense (instant rise, exponential decay set
by ambipolar diffusion —
underdense_echo) and overdense (a longer ragged plateau from a dense, metallic-like trail —overdense_echo). simulate_meteor_timeseriesbuilds a realistic noisy record with Poisson-random pings; the same events appear as bright dashes in the Chapter 5 spectrogram.detect_pings(returning aMeteorDetection) counts pings by MAD-robust thresholding. The threshold is a sensitivity vs false-alarm trade-off; the scientific product is not one count but the meteor rate over a night or shower, which amateurs contribute to the IMO/RMOB networks.- Passive radar generalises the idea: cross-correlate a reference and a
surveillance channel to get a target's range and Doppler from a
transmitter you don't control — the realm of
blah2(Field Notes).
What's next¶
This was an optional hardware side-quest off the main interferometry arc. The
coherence problem behind passive radar's two channels is the same one that makes
real interferometry hard — picked up in
Chapter 17 (Coherent Interferometry with KrakenSDR)
and the Field Notes. For more buildable projects — IBT,
hydrogen-line horns, VLF, and the meteor and passive-radar hacks above — see
docs/projects.md. And it is always legal to just
listen: tune a WebSDR right now and, if you're in
Europe, see if you can catch a GRAVES ping yourself.