Chapter 23 — Solar & Jupiter Radio Astronomy¶
!!! info "Before you start" Prerequisites: Ch 2 (The Physics of Radio Emission) · Maths Lab: none · ~40 min · Intermediate
The Sun is, by a wide margin, the brightest radio source in the sky — and Jupiter is the only planet you can hear with a wire antenna in your backyard. Both shout loudest in the low-frequency, decametric band around 20 MHz, which is exactly why they are the twin targets of NASA's Radio JOVE project: a build-it-yourself dual-dipole antenna and an SDR that lets a student capture solar radio bursts and Jovian storms with real, archivable data.
This chapter is the science behind Radio JOVE. We will:
- meet the Sun as an intense, variable radio source (Hey & Southworth, 1946);
- see why a low observing frequency probes high in the corona, through the plasma frequency $f_p \propto \sqrt{n_e}$ and a coronal density model;
- simulate a type II burst — the radio fingerprint of a coronal-mass-ejection (CME) shock — and derive the shock speed from its frequency drift;
- simulate Jupiter's decametric (DAM) emission, the Io-controlled storms that Radio JOVE was built to hear.
Everything runs fully offline on the base environment using simulation: the
physics lives in jansky.solar, and we add realistic noise with
jansky.signals.rng. No telescope, no network.
Learning goals¶
By the end of this chapter you will be able to
- explain why the Sun and Jupiter are observed near 20 MHz, and what part of the corona a given frequency samples;
- relate electron density, plasma frequency, and height in the corona with the Newkirk model;
- read a dynamic spectrum (frequency vs time waterfall) and recognise a frequency-drifting burst;
- turn a type II burst's drift rate into a CME shock speed, and understand how the assumed density model biases that number;
- describe Jovian L-bursts and S-bursts and the Io–Jupiter interaction.
This chapter connects back to Chapter 16 (data formats): the spectrograms we
build here are exactly the Radio-Sky Spectrograph (.sps) products that the
Radio JOVE community archives and shares.
1. Orientation — the Sun (and a planet) burst onto the radio sky¶
Karl Jansky heard the Milky Way, not the Sun: in the early 1930s the Sun happened to be near solar minimum and quiet. The Sun's radio debut came with World War II radar. In February 1942, British radar operators were jammed by a powerful noise source that turned out to be the Sun during an intense sunspot group — work led by J. S. Hey, published once secrecy lifted as "Solar Radiations in the 4–6 Metre Radio Wave-Length Band" (Hey 1946, Nature 157, 47). In parallel, G. C. Southworth had already detected the quiet Sun's thermal microwave emission during the war (Southworth 1945, J. Franklin Inst. 239, 285). Together they revealed two faces of the radio Sun: a steady thermal glow, and violent, sporadic bursts tied to flares and active regions. See the papers timeline for both, alongside Ryle & Vonberg's 1946 first interferometric solar observations.
A decade later, Burke & Franklin (1955, JGR 60, 213) stumbled on intense bursts at 22.2 MHz that tracked with the rotation of Jupiter — the first detection of radio emission from a planet, and the birth of planetary radio astronomy.
Why observe near 20 MHz?¶
A low frequency is not an accident — it is forced by the physics. A plasma cannot transmit an electromagnetic wave below its plasma frequency $f_p$; below $f_p$ the wave is reflected or absorbed. The corona's density (and so its $f_p$) falls with height, so each observing frequency "sees down" only to the layer where the local $f_p$ matches. Observe at 20 MHz and you are probing the corona a few solar radii out — and you are below Earth's ionospheric cutoff only on good nights, which is its own amateur sport. This is Section 2.
The solar burst zoo (types I–V), in one paragraph¶
Solar radio bursts are classified I–V by their appearance on a dynamic spectrum (frequency vs time): Type I — long-lived noise storms of many short bursts over active regions; Type II — slow frequency drift (minutes), the signature of a shock (often a CME) climbing through the corona; Type III — fast drift (seconds), beams of electrons streaming outward along open field lines; Type IV — broadband continuum from a moving or stationary plasmoid after a flare; Type V — short continuum sometimes following type III. The star of this chapter is the type II: because its drift is slow and tied to a physical shock moving through a known density profile, you can measure the shock's speed — "low-hanging real science" with a ~1 m antenna (see the field notes).
Setup¶
We import the plasma / coronal-density toolkit from jansky.solar and the
course plotting style. We do not re-implement the plasma-frequency relation,
the Newkirk density model, or the burst-drift forward/inverse calculation —
every numerical step below calls a helper from that module. We seed every random
draw with jansky.signals.rng so the chapter is reproducible.
The helpers we use:
plasma_frequency(n_e)— electron plasma frequency $f_p$ [MHz] from density [cm$^{-3}$];density_from_plasma_frequency(f)— its inverse;newkirk_density(r, fold)— Newkirk (1961) coronal density at heliocentric radius $r$ [solar radii];newkirk_radius(n_e, fold)— inverse Newkirk: radius for a density;type_ii_track(speed, ...)— forward-model a type II's frequency-time drift;shock_speed_from_track(times, freqs, ...)— recover the shock speed from a drift;R_SUN_KM— the solar radius in km.
import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
from jansky.plotting import use_jansky_style, show_image
from jansky.signals import rng
from jansky.solar import (
R_SUN_KM,
plasma_frequency,
density_from_plasma_frequency,
newkirk_density,
newkirk_radius,
type_ii_track,
shock_speed_from_track,
)
use_jansky_style()
SEED = 23
generator = rng(SEED) # one seeded generator for the whole chapter
print(f"solar radius R_sun = {R_SUN_KM:,.0f} km = {(R_SUN_KM * u.km).to(u.Rsun):.3f}")
print(f"plasma frequency at n_e = 1e8 cm^-3 : {plasma_frequency(1e8):.1f} MHz")
solar radius R_sun = 695,700 km = 1.000 solRad plasma frequency at n_e = 1e8 cm^-3 : 89.8 MHz
2. Plasma frequency & the corona¶
The key relation for everything in this chapter is the electron plasma frequency: the natural oscillation frequency of a cloud of free electrons,
$$ f_p \;=\; \frac{1}{2\pi}\sqrt{\frac{n_e\,e^2}{\varepsilon_0\,m_e}} \;\approx\; 8.98\times10^{-3}\,\sqrt{n_e\,[\mathrm{cm^{-3}}]}\ \ \mathrm{MHz}. $$
The headline is the square-root: $f_p \propto \sqrt{n_e}$. A plasma is opaque to waves below $f_p$, so the radio emission we receive at frequency $\nu$ was generated where the local plasma frequency is $\nu$ (the fundamental) or $\nu/2$ (the harmonic, $2 f_p$). To turn frequency into height, we need how density falls with height — the Newkirk (1961) coronal model,
$$ n_e(r) \;=\; f \cdot 4.2\times10^{4}\;\,10^{\,4.32/r}\ \ \mathrm{cm^{-3}}, $$
with $r$ in solar radii from Sun centre. The fold factor $f$ (1–4) scales the model up for denser streamers and active regions. Putting the two together gives $f_p(r)$: a curve that falls steeply with height. Low frequency ⇒ high in the corona. Let us plot all three.
# Heliocentric radius grid, from just above the surface to ~10 solar radii.
r = np.linspace(1.05, 10.0, 400) # [solar radii]
n_e = newkirk_density(r, fold=1.0) # [cm^-3]
fp = plasma_frequency(n_e) # fundamental plasma frequency [MHz]
# Common observing references. The Sun emits at the FUNDAMENTAL (f_p) and the
# HARMONIC (2 f_p); the height a band probes depends on which you assume.
jove_band = (18.0, 28.0) # Radio JOVE listens over ~8 MHz around ~20-24 MHz
f_20 = 20.0 # a round "metre-wave" reference
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# (a) density vs height
axes[0].semilogy(r, n_e, color="#d62728", lw=2)
axes[0].set_xlabel("heliocentric radius $r$ [$R_\\odot$]")
axes[0].set_ylabel("electron density $n_e$ [cm$^{-3}$]")
axes[0].set_title("Newkirk (1961) coronal density")
# (b) plasma frequency vs height, with observing bands marked
axes[1].semilogy(r, fp, color="#1f77b4", lw=2, label="fundamental $f_p$")
axes[1].semilogy(r, 2 * fp, color="#1f77b4", lw=1.5, ls="--", label="harmonic $2 f_p$")
axes[1].axhspan(*jove_band, color="#2ca02c", alpha=0.18, label="Radio JOVE band (~18-28 MHz)")
axes[1].axhline(f_20, color="#7f7f7f", ls=":", lw=1.2, label="20 MHz")
axes[1].set_xlabel("heliocentric radius $r$ [$R_\\odot$]")
axes[1].set_ylabel("plasma frequency [MHz]")
axes[1].set_title("Plasma frequency vs height")
axes[1].legend(fontsize=9)
plt.tight_layout()
plt.show()
Where does 20 MHz come from? Invert the model: find the height at which the
fundamental (and the harmonic) plasma frequency equals our observing references.
This is newkirk_radius composed with density_from_plasma_frequency.
def height_of_frequency(f_mhz, harmonic=1, fold=1.0):
"Heliocentric radius (R_sun) at which emission at f_mhz originates."
f_p = np.asarray(f_mhz, dtype=float) / harmonic
n_e = density_from_plasma_frequency(f_p)
return newkirk_radius(n_e, fold=fold)
for f in (20.0, 24.0, 80.0):
r_fund = height_of_frequency(f, harmonic=1)
r_harm = height_of_frequency(f, harmonic=2)
print(
f"{f:5.0f} MHz -> fundamental from r = {r_fund:.2f} R_sun, "
f"harmonic (2 f_p) from r = {r_harm:.2f} R_sun"
)
print()
print(
"So 20 MHz fundamental emission is born around "
f"{height_of_frequency(20.0):.1f} solar radii out -- well up in the corona,"
)
print("which is why low frequency = high altitude. A 200 MHz burst, by contrast,")
print(f"comes from r = {height_of_frequency(200.0):.2f} R_sun -- right down near the surface.")
20 MHz -> fundamental from r = 2.08 R_sun, harmonic (2 f_p) from r = 2.94 R_sun 24 MHz -> fundamental from r = 1.94 R_sun, harmonic (2 f_p) from r = 2.65 R_sun 80 MHz -> fundamental from r = 1.32 R_sun, harmonic (2 f_p) from r = 1.62 R_sun So 20 MHz fundamental emission is born around 2.1 solar radii out -- well up in the corona, which is why low frequency = high altitude. A 200 MHz burst, by contrast, comes from r = 1.06 R_sun -- right down near the surface.
Notice the asymmetry the square-root creates: density spans orders of magnitude from the low corona to a few $R_\odot$, but plasma frequency, going as $\sqrt{n_e}$, falls more gently. A factor-of-100 drop in density is only a factor-of-10 drop in $f_p$. That is what lets a single ~8 MHz-wide receiver, parked near 20 MHz, watch a disturbance climb across a real swath of the corona — the engine of the type II measurement in Section 3.
3. A type II burst → a CME shock speed¶
Now the payoff. A coronal mass ejection drives a shock outward through the corona. At the shock front, electrons are accelerated and the local plasma rings at its plasma frequency, radiating at $f_p$ and $2 f_p$. As the shock climbs, the density (and so $f_p$) drops, so the emission drifts down in frequency. That slow downward drift on the dynamic spectrum is the type II burst.
type_ii_track(speed_km_s, r_start, harmonic, fold) forward-models exactly this:
a shock moving radially at constant speed from r_start, emitting at
harmonic$\times f_p$ at each height, using the Newkirk density. We pick a CME
shock at 700 km/s — a fairly typical value — starting at $1.5\,R_\odot$ and
emitting at the harmonic ($2 f_p$), and read off its frequency-time track.
V_TRUE = 700.0 # true CME shock speed [km/s]
R_START = 1.5 # shock launches at 1.5 solar radii
HARMONIC = 2 # observe the harmonic band (2 f_p)
DURATION = 600.0 # 10 minutes of drift
times, freqs = type_ii_track(
V_TRUE,
r_start=R_START,
duration_s=DURATION,
harmonic=HARMONIC,
fold=1.0,
)
print(f"Type II track for a {V_TRUE:.0f} km/s shock from r = {R_START} R_sun:")
print(
f" starts at {freqs[0]:.1f} MHz, drifts down to {freqs[-1]:.1f} MHz "
f"over {DURATION / 60:.0f} min"
)
drift = (freqs[0] - freqs[-1]) / DURATION
print(f" mean drift rate: {-drift * 1e3:.2f} kHz/s (slow -- the hallmark of a type II)")
Type II track for a 700 km/s shock from r = 1.5 R_sun: starts at 101.3 MHz, drifts down to 39.1 MHz over 10 min mean drift rate: -103.68 kHz/s (slow -- the hallmark of a type II)
Make it look like real data: a dynamic spectrum with noise¶
A receiver does not hand you a clean curve; it hands you a waterfall — power
in each (time, frequency) pixel, mostly noise, with the burst a bright ridge
tracing the track. We paint the track onto a noisy time–frequency grid with a
Gaussian profile in frequency, exactly as signals.synthetic_spectrum adds
Gaussian noise elsewhere in the course.
# Build a (frequency, time) dynamic spectrum around the burst's band.
f_axis = np.linspace(20.0, 110.0, 240) # MHz (receiver passband)
t_axis = np.linspace(0.0, DURATION, 300) # s
# Interpolate the burst's frequency at each time sample.
f_burst = np.interp(t_axis, times, freqs)
# Start from thermal noise, then add a Gaussian-in-frequency ridge at f_burst(t).
dynspec = generator.normal(0.0, 1.0, size=(f_axis.size, t_axis.size))
ridge_width = 2.5 # MHz, the burst's instantaneous bandwidth
amplitude = 8.0 # peak in units of the noise sigma
FF, _ = np.meshgrid(f_axis, t_axis, indexing="ij")
dynspec += amplitude * np.exp(-0.5 * ((FF - f_burst[None, :]) / ridge_width) ** 2)
fig, ax = plt.subplots(figsize=(10, 5.5))
extent = [0.0, DURATION / 60.0, f_axis[0], f_axis[-1]] # x=min, y=MHz
show_image(
dynspec,
ax=ax,
aspect="auto",
extent=extent,
title="Simulated type II burst — a slow drift to lower frequency",
)
ax.plot(times / 60.0, freqs, color="cyan", lw=1.3, ls="--", label="true shock track")
ax.set_xlabel("time [min]")
ax.set_ylabel("frequency [MHz]")
ax.grid(False)
ax.legend(loc="upper right")
plt.tight_layout()
plt.show()
Recover the shock speed¶
This is the science. The drift rate encodes the shock's speed through the density
gradient. shock_speed_from_track inverts the whole chain: each observed
frequency → plasma frequency (÷ harmonic) → density → height (Newkirk) → a
straight-line fit of height vs time, whose slope is the radial speed. We feed
it the burst track as if we had measured it from the waterfall and compare to
the true input.
v_recovered = shock_speed_from_track(times, freqs, harmonic=HARMONIC, fold=1.0)
print(f"true shock speed : {V_TRUE:.0f} km/s")
print(f"recovered shock speed : {v_recovered:.0f} km/s")
print(f"error : {100 * (v_recovered - V_TRUE) / V_TRUE:+.2f} %")
print()
print("With the SAME density model and harmonic we used to make the burst, the")
print("recovery is essentially exact -- the method is self-consistent. Real life")
print("is messier: you must ASSUME a coronal model, and that assumption bites.")
true shock speed : 700 km/s recovered shock speed : 700 km/s error : -0.00 % With the SAME density model and harmonic we used to make the burst, the recovery is essentially exact -- the method is self-consistent. Real life is messier: you must ASSUME a coronal model, and that assumption bites.
Let us see the inversion geometrically: convert the track to height vs time and confirm it is a straight line whose slope is the speed.
# Reproduce the inversion by hand to visualise it.
f_p = freqs / HARMONIC
n_e_track = density_from_plasma_frequency(f_p)
r_track = newkirk_radius(n_e_track) # solar radii
slope, intercept = np.polyfit(times, r_track, 1) # R_sun / s
v_fit = slope * R_SUN_KM # km/s
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(
times / 60.0, r_track, "o", ms=3, color="#1f77b4", label="height inferred from each frequency"
)
ax.plot(
times / 60.0,
(intercept + slope * times),
color="#d62728",
lw=2,
label=f"linear fit -> {v_fit:.0f} km/s",
)
ax.set_xlabel("time [min]")
ax.set_ylabel("heliocentric radius [$R_\\odot$]")
ax.set_title("Type II inversion: a straight climb through the corona")
ax.legend()
plt.tight_layout()
plt.show()
Typical CME speeds run from a leisurely ~250 km/s up to the fastest Earth-directed events at 2000+ km/s; the fastest can reach Earth in under a day and drive geomagnetic storms. The type II's drift rate is one of the few ways to estimate the shock speed while it is still in the corona, from the ground, at amateur frequencies. Let us confirm faster shocks drift faster.
print(
f"{'V_in [km/s]':>12} {'start MHz':>10} {'end MHz':>9} {'drift kHz/s':>12} {'V_out [km/s]':>13}"
)
for v_in in (300.0, 700.0, 1200.0, 2000.0):
t_k, f_k = type_ii_track(v_in, r_start=R_START, duration_s=DURATION, harmonic=HARMONIC)
drift_rate = -(f_k[-1] - f_k[0]) / DURATION * 1e3 # kHz/s
v_out = shock_speed_from_track(t_k, f_k, harmonic=HARMONIC)
print(f"{v_in:12.0f} {f_k[0]:10.1f} {f_k[-1]:9.1f} {drift_rate:12.2f} {v_out:13.0f}")
V_in [km/s] start MHz end MHz drift kHz/s V_out [km/s]
300 101.3 62.2 65.20 300
700 101.3 39.1 103.68 700
1200 101.3 26.2 125.28 1200
2000 101.3 17.2 140.23 2000
The density model is an assumption — and it biases the speed¶
The catch the field notes warned about: you never know the coronal density;
you assume a model. If the real corona is denser than you think (a streamer,
an active region — a larger fold), but you invert with fold=1, your inferred
heights are wrong and so is your speed. Below we generate a burst in a corona
that is 3× denser (fold=3, e.g. a dense streamer) but analyse it with
the default fold=1, and watch the recovered speed go astray.
# Truth: a 700 km/s shock in a DENSE (fold=3) corona.
t_dense, f_dense = type_ii_track(
V_TRUE, r_start=R_START, duration_s=DURATION, harmonic=HARMONIC, fold=3.0
)
# Analyst's choices: the correct fold vs the (wrong) default fold=1.
v_correct = shock_speed_from_track(t_dense, f_dense, harmonic=HARMONIC, fold=3.0)
v_wrongfold = shock_speed_from_track(t_dense, f_dense, harmonic=HARMONIC, fold=1.0)
# And: what if you mistake the harmonic band for the fundamental?
v_wrongharm = shock_speed_from_track(t_dense, f_dense, harmonic=1, fold=3.0)
print(f"true speed : {V_TRUE:.0f} km/s")
print(f"recovered with CORRECT model (fold=3, h=2) : {v_correct:.0f} km/s")
print(
f"recovered assuming WRONG density (fold=1) : {v_wrongfold:.0f} km/s "
f"({100 * (v_wrongfold - V_TRUE) / V_TRUE:+.0f}%)"
)
print(
f"recovered assuming WRONG harmonic (h=1) : {v_wrongharm:.0f} km/s "
f"({100 * (v_wrongharm - V_TRUE) / V_TRUE:+.0f}%)"
)
print()
print("Density fold and harmonic both rescale the inferred heights; a wrong")
print("choice systematically biases the speed. Real type II analyses quote the")
print("coronal model used precisely because the answer depends on it.")
true speed : 700 km/s recovered with CORRECT model (fold=3, h=2) : 700 km/s recovered assuming WRONG density (fold=1) : 487 km/s (-30%) recovered assuming WRONG harmonic (h=1) : 448 km/s (-36%) Density fold and harmonic both rescale the inferred heights; a wrong choice systematically biases the speed. Real type II analyses quote the coronal model used precisely because the answer depends on it.
!!! tip "Contribute & compare — ESA Solar Radio Burst Tracker" The drifting dynamic-spectrum bursts in this section are exactly what ESA's Solar Radio Burst Tracker citizen-science project asks you to find: you annotate Type III solar radio bursts in real Solar Orbiter / RPW dynamic spectra to build the mission's first burst catalogue. No hardware or expertise needed — a direct encounter with real spacecraft solar-radio data.
4. Jupiter's decametric (DAM) emission¶
Switch targets. Jupiter is a strong, sporadic source of decametric (DAM) radio emission — wavelengths of tens of metres, i.e. frequencies up to ~40 MHz with most activity below ~30 MHz. It is generated by the cyclotron maser instability: electrons spiralling along Jupiter's enormous magnetic field near its poles radiate coherently. Crucially, much of the emission is controlled by Io: the volcanic moon ploughs through Jupiter's magnetosphere, and the Io–Jupiter flux tube of current it sets up modulates the storms. Whether you hear Jupiter depends on the System III longitude (which way the magnetic field points) and on Io's orbital phase — the famous "Io-A/B/C" storm regions.
On a dynamic spectrum the storms show two characteristic features:
- L-bursts ("long") — smooth, second-to-minute-scale rises and falls; on an SDR they sound like ocean waves washing in and out;
- S-bursts ("short") — millisecond crackles that drift rapidly in frequency; they sound like popcorn or a crackling fire.
Both live in the same ~18–28 MHz band Radio JOVE listens to. Let us simulate a Jovian storm spectrogram with pure NumPy — no new helper needed — by laying L-burst blobs and fast-drifting S-burst streaks onto a noisy time–frequency grid.
# Jovian DAM spectrogram: 6 minutes over the Radio JOVE band.
f_jove = np.linspace(16.0, 30.0, 280) # MHz
t_jove = np.linspace(0.0, 360.0, 1400) # s (fine sampling for S-bursts)
JJ, TT = np.meshgrid(f_jove, t_jove, indexing="ij") # (freq, time)
jov = generator.normal(0.0, 1.0, size=(f_jove.size, t_jove.size)) # noise floor
# --- L-bursts: broadband, smooth, tens-of-seconds envelopes -----------------
def add_lburst(grid, t0, dur, f_lo, f_hi, amp):
"Smooth broadband swell centred at time t0."
t_env = np.exp(-0.5 * ((TT - t0) / (dur / 2.355)) ** 2)
f_env = ((JJ > f_lo) & (JJ < f_hi)).astype(float)
# soften the band edges
f_env = np.clip(f_env + 0.4 * np.exp(-0.5 * ((JJ - 0.5 * (f_lo + f_hi)) / 4.0) ** 2), 0, 1)
grid += amp * t_env * f_env
for t0, dur, flo, fhi, amp in [
(70.0, 40.0, 18.0, 27.0, 4.0),
(150.0, 55.0, 17.0, 24.0, 3.0),
(250.0, 35.0, 20.0, 29.0, 4.5),
]:
add_lburst(jov, t0, dur, flo, fhi, amp)
# --- S-bursts: fast NEGATIVE frequency drift, millisecond crackles ----------
s_generator = rng(SEED + 1)
n_sbursts = 240
drift_mhz_per_s = -25.0 # S-bursts drift down ~ -10 to -40 MHz/s
for _ in range(n_sbursts):
ts = s_generator.uniform(20.0, 340.0) # onset time [s]
fs = s_generator.uniform(22.0, 29.0) # onset frequency [MHz]
life = s_generator.uniform(0.05, 0.4) # duration [s]
amp = s_generator.uniform(5.0, 9.0)
f_of_t = fs + drift_mhz_per_s * (TT - ts)
mask = (TT >= ts) & (TT <= ts + life)
jov += amp * mask * np.exp(-0.5 * ((JJ - f_of_t) / 0.35) ** 2)
print(f"Jovian spectrogram: {jov.shape[0]} channels x {jov.shape[1]} time samples")
print(
f" {3} L-burst swells + {n_sbursts} S-burst crackles over "
f"{t_jove[-1] / 60:.0f} min, all in the {f_jove[0]:.0f}-{f_jove[-1]:.0f} MHz JOVE band"
)
Jovian spectrogram: 280 channels x 1400 time samples 3 L-burst swells + 240 S-burst crackles over 6 min, all in the 16-30 MHz JOVE band
fig, ax = plt.subplots(figsize=(11, 5.5))
extent = [0.0, t_jove[-1], f_jove[0], f_jove[-1]]
show_image(
jov,
ax=ax,
aspect="auto",
extent=extent,
vmax=6.0,
title="Simulated Jovian decametric storm (Radio JOVE band)",
)
ax.set_xlabel("time [s]")
ax.set_ylabel("frequency [MHz]")
ax.grid(False)
# Annotate the two burst types.
ax.annotate(
"L-burst\n(smooth swell)",
xy=(70, 22.5),
xytext=(15, 17.0),
color="white",
fontsize=9,
arrowprops=dict(arrowstyle="->", color="white"),
)
ax.annotate(
"S-bursts\n(fast crackles)",
xy=(300, 27.0),
xytext=(250, 16.7),
color="white",
fontsize=9,
arrowprops=dict(arrowstyle="->", color="white"),
)
plt.tight_layout()
plt.show()
Zoom into one second to see the S-bursts resolve into individual fast-drifting streaks — the millisecond fine structure that needs a fast, wide-bandwidth receiver to catch (and whose physics is still not fully settled).
# Zoom: a 1.5 s window rich in S-bursts.
tmask = (t_jove >= 298.0) & (t_jove <= 304.0)
fig, ax = plt.subplots(figsize=(10, 4.5))
extent_zoom = [t_jove[tmask][0], t_jove[tmask][-1], f_jove[0], f_jove[-1]]
show_image(
jov[:, tmask],
ax=ax,
aspect="auto",
extent=extent_zoom,
vmax=7.0,
title="Zoom: S-bursts resolve into fast-drifting streaks",
)
ax.set_xlabel("time [s]")
ax.set_ylabel("frequency [MHz]")
ax.grid(False)
plt.tight_layout()
plt.show()
What a Radio JOVE station actually detects. A Radio JOVE kit is a pair of
half-wave dipoles strung between poles (a phased dual-dipole aimed at the
meridian) feeding an SDRplay RSP1B SDR, recording an ~8 MHz band around
20–24 MHz with Radio-Sky Spectrograph software. Point it on a night when
Jupiter is up and Io is in an active phase, and the waterfall you record looks
much like the figures above: L-burst swells and S-burst crackles riding on the
galactic background. The same station catches solar type III and type II bursts
by day. See Projects, Kits & Hacks for the kit, and
Chapter 16 for the .sps file format these
spectrograms are saved and shared in.
Now the real thing: an actual .sps recording¶
Everything above was simulated. Here we load a real Radio-Sky Spectrograph
.sps file — a dual-polarisation 16–32 MHz recording from the AJ4CO Radio
JOVE station (Dave Typinski, High Springs FL) — with jansky.formats.read_sps,
and plot its dynamic spectrum. After subtracting each channel's time-median (to
flatten the receiver bandpass), the time-varying decametric sky emerges: the
slowly drifting galactic background, terrestrial interference, and — on an
active night — the solar and Jovian bursts whose physics we modelled above. This
is exactly the file format the Radio JOVE community archives and shares
(Chapter 16); read_sps returns it as a
Spectrogram with both polarisation feeds in meta["feeds"].
The file is fetched once (~12 MB) and cached; if you are offline the cell prints a note and the simulated figures above stand in.
# --- A REAL Radio JOVE observation, read with jansky.formats.read_sps --------
from jansky import data, formats
try:
sps_path = data.fetch("radiojove-sps") # ~12 MB, cached after the first run
obs = formats.read_sps(sps_path)
except Exception as exc: # offline / unavailable -> the simulated figures stand
obs = None
print(f"(Real .sps unavailable: {exc})")
print("This needs a one-time ~12 MB download; the simulated figures above stand.")
if obs is not None:
m = obs.meta
feed = m["feed_names"][0]
print(f"{m['obsname']} - {m['author']}, {m['obsloc']}")
print(
f"{m['nchannels']} channels, "
f"{m['fmax_hz'] / 1e6:.0f}-{m['fmin_hz'] / 1e6:.0f} MHz, "
f"{obs.power.shape[0]:,} sweeps over {obs.times[-1] / 60:.0f} min "
f"(feeds: {m['feed_names']})"
)
pw = obs.power # (n_sweep, n_chan) for the first feed
# Subtract each channel's time-median to flatten the receiver bandpass and
# reveal the time-varying sky (the standard Radio JOVE background subtraction).
disp = pw - np.median(pw, axis=0, keepdims=True)
disp = disp[:, ::-1] # reorder channels ascending: 16 -> 32 MHz
lo, hi = np.percentile(disp, [2.5, 99.5])
fig, ax = plt.subplots(figsize=(11, 5))
show_image(
disp.T,
ax=ax,
aspect="auto",
extent=[0.0, obs.times[-1] / 60.0, m["fmin_hz"] / 1e6, m["fmax_hz"] / 1e6],
vmin=lo,
vmax=hi,
title=f"Real Radio JOVE .sps observation - {m['author']} ({feed} feed)",
)
ax.set_xlabel("time [min]")
ax.set_ylabel("frequency [MHz]")
ax.grid(False)
plt.tight_layout()
plt.show()
AJ4CO DPS - Dave Typinski, High Springs, FL 300 channels, 32-16 MHz, 10,412 sweeps over 27 min (feeds: ['RR', 'LL'])
Predicting DAM storm windows with JPL Horizons¶
Jupiter's decametric (DAM) storms are not random: their probability is controlled by two
geometric clocks, both of which JPL Horizons
will compute for any night and site (via astroquery.jplhorizons):
- Central-meridian longitude (CML, System III) — which Jovian magnetic longitude faces Earth.
Horizons reports it directly as the sub-observer planetodetic longitude (
PDObsLon); it sweeps $360^\circ$ every $9.925$ h (Jupiter's radio rotation period). - Io phase — where Io sits in its $42.46$ h orbit, measured from superior geocentric conjunction. Io's flux tube drives the strongest sources, so DAM clusters into the classic Io-A / Io-B / Io-C regions in the $(\mathrm{CML},\,\text{Io phase})$ plane (Bigg 1964; Carr, Desch & Alexander 1983; ranges after Louis et al. 2017).
We query Horizons for Jupiter (id 599) and Io (id 501) over a night (here Jupiter's 2026-01-11 opposition, from Goldstone), derive Io phase from the two bodies' line-of-sight and transverse offsets, and flag the hours when Jupiter is **up *and*** the geometry lands in a DAM source region — exactly the prediction a Radio JOVE or GAVRT observer makes before a session. As everywhere in this course, the cell falls back to bundled values so it runs fully offline.
import numpy as np
from jansky import plotting
plotting.use_jansky_style()
AU_KM = 1.495978707e8
SITE = {"lon": -116.805, "lat": 35.247, "elevation": 1.0} # Goldstone, CA (deg E, deg N, km)
NIGHT = {
"start": "2026-01-11 00:00",
"stop": "2026-01-12 00:00",
"step": "20m",
} # Jupiter opposition
def io_phase_deg(ra_j, dec_j, delta_j, ra_i, delta_i):
"""Io phase in degrees: 0 deg at superior geocentric conjunction, increasing with Io's
(prograde) orbital motion. Built from Io's offset from Jupiter resolved into a transverse
(on-sky, east) component and a line-of-sight (range) component."""
dra = ((ra_i - ra_j) + 180.0) % 360.0 - 180.0
x_km = np.radians(dra * np.cos(np.radians(dec_j))) * (delta_j * AU_KM) # +east (transverse)
z_km = (delta_i - delta_j) * AU_KM # +behind (superior conj.)
return np.degrees(np.arctan2(x_km, z_km)) % 360.0
try:
from astroquery.jplhorizons import Horizons
jup = Horizons(id="599", location=SITE, epochs=NIGHT, id_type=None).ephemerides()
io = Horizons(id="501", location=SITE, epochs=NIGHT, id_type=None).ephemerides()
jd = np.array(jup["datetime_jd"], dtype=float)
hours = (jd - jd[0]) * 24.0
elev = np.array(jup["EL"], dtype=float) # altitude above horizon, deg
cml = np.array(jup["PDObsLon"], dtype=float) # System III CML, deg
iph = io_phase_deg(
np.array(jup["RA"], dtype=float),
np.array(jup["DEC"], dtype=float),
np.array(jup["delta"], dtype=float),
np.array(io["RA"], dtype=float),
np.array(io["delta"], dtype=float),
)
source = "live JPL Horizons"
except Exception as exc: # offline: bundled Horizons output for 2026-01-11, Goldstone
print(f"Horizons unavailable ({type(exc).__name__}); using the bundled fallback.")
hours = np.array(
[
0.0,
0.333,
0.667,
1.0,
1.333,
1.667,
2.0,
2.333,
2.667,
3.0,
3.333,
3.667,
4.0,
4.333,
4.667,
5.0,
5.333,
5.667,
6.0,
6.333,
6.667,
7.0,
7.333,
7.667,
8.0,
8.333,
8.667,
9.0,
9.333,
9.667,
10.0,
10.333,
10.667,
11.0,
11.333,
11.667,
12.0,
12.333,
12.667,
13.0,
13.333,
13.667,
14.0,
14.333,
14.667,
15.0,
15.333,
15.667,
16.0,
16.333,
16.667,
17.0,
17.333,
17.667,
18.0,
18.333,
18.667,
19.0,
19.333,
19.667,
20.0,
20.333,
20.667,
21.0,
21.333,
21.667,
22.0,
22.333,
22.667,
23.0,
23.333,
23.667,
24.0,
]
)
elev = np.array(
[
-7.9,
-4.4,
-0.9,
2.8,
6.5,
10.4,
14.2,
18.2,
22.2,
26.2,
30.2,
34.3,
38.4,
42.5,
46.6,
50.7,
54.7,
58.7,
62.6,
66.3,
69.8,
72.9,
75.4,
76.8,
76.8,
75.5,
73.0,
70.0,
66.5,
62.8,
58.9,
54.9,
50.8,
46.8,
42.7,
38.6,
34.5,
30.4,
26.4,
22.3,
18.4,
14.4,
10.5,
6.7,
3.0,
-0.7,
-4.2,
-7.7,
-11.0,
-14.2,
-17.2,
-20.0,
-22.6,
-24.9,
-27.0,
-28.8,
-30.3,
-31.4,
-32.1,
-32.5,
-32.5,
-32.0,
-31.2,
-30.0,
-28.5,
-26.7,
-24.5,
-22.1,
-19.5,
-16.7,
-13.6,
-10.4,
-7.1,
]
)
cml = np.array(
[
43.4,
55.5,
67.6,
79.7,
91.8,
103.9,
116.0,
128.1,
140.2,
152.3,
164.3,
176.4,
188.5,
200.6,
212.7,
224.8,
236.9,
249.0,
261.1,
273.2,
285.3,
297.4,
309.5,
321.6,
333.6,
345.7,
357.8,
9.9,
22.0,
34.1,
46.2,
58.3,
70.4,
82.5,
94.6,
106.7,
118.8,
130.8,
142.9,
155.0,
167.1,
179.2,
191.3,
203.4,
215.5,
227.6,
239.7,
251.8,
263.9,
276.0,
288.1,
300.1,
312.2,
324.3,
336.4,
348.5,
0.6,
12.7,
24.8,
36.9,
49.0,
61.1,
73.2,
85.3,
97.3,
109.4,
121.5,
133.6,
145.7,
157.8,
169.9,
182.0,
194.1,
]
)
iph = np.array(
[
91.4,
94.2,
97.1,
100.0,
102.9,
105.8,
108.7,
111.6,
114.5,
117.4,
120.3,
123.1,
126.0,
128.9,
131.8,
134.6,
137.5,
140.3,
143.2,
146.0,
148.9,
151.7,
154.5,
157.4,
160.2,
163.0,
165.8,
168.6,
171.4,
174.3,
177.0,
179.8,
182.6,
185.4,
188.2,
191.0,
193.8,
196.6,
199.4,
202.2,
205.0,
207.8,
210.6,
213.4,
216.2,
219.0,
221.8,
224.6,
227.4,
230.2,
233.0,
235.8,
238.6,
241.4,
244.3,
247.1,
249.9,
252.8,
255.6,
258.4,
261.3,
264.1,
267.0,
269.8,
272.7,
275.5,
278.4,
281.2,
284.0,
286.9,
289.7,
292.6,
295.4,
]
)
source = "offline fallback (JPL Horizons — Jupiter & Io, 2026-01-11 UT, Goldstone)"
up = elev > 0.0
print(f"Data source: {source}")
print(f"Jupiter above the horizon for {np.count_nonzero(up) * 20 / 60:.1f} h of this night.")
print(
f"CML sweeps {cml.min():.0f}–{cml.max():.0f}°; Io phase {iph[up].min():.0f}"
f"–{iph[up].max():.0f}° while Jupiter is up."
)
# Classic Io-controlled DAM source regions in (CML_III, Io phase) -- approximate occurrence-
# probability boxes (Bigg 1964; Carr, Desch & Alexander 1983; ranges after the 26-yr Nancay
# statistics of Louis et al. 2017). Io-C wraps past 360 deg in CML.
IO_SOURCES = {
"Io-B": {"cml": (90, 180), "io": (80, 140), "color": "#0072B2"},
"Io-A": {"cml": (195, 290), "io": (195, 265), "color": "#D55E00"},
"Io-C": {"cml": (280, 420), "io": (260, 340), "color": "#009E73"},
}
def dam_region(cml_deg, io_deg):
"""Return the Io-DAM source label for a (CML, Io phase) point, or None."""
for name, b in IO_SOURCES.items():
lo, hi = b["cml"]
c = cml_deg + 360.0 if (hi > 360.0 and cml_deg < lo) else cml_deg
if lo <= c <= hi and b["io"][0] <= io_deg <= b["io"][1]:
return name
return None
labels = np.array([dam_region(c, p) if u else None for c, p, u in zip(cml, iph, up)])
hit = labels != None # noqa: E711 (object array; compare to None elementwise)
def _hhmm(h):
m = int(round(h * 60))
return f"{m // 60:02d}:{m % 60:02d}"
# Group consecutive in-region samples into separate windows (CML cycles every ~10 h,
# so a given Io source can light up more than once a night).
windows, i, n = [], 0, len(hit)
while i < n:
if hit[i]:
j = i
while j + 1 < n and hit[j + 1] and labels[j + 1] == labels[i]:
j += 1
windows.append((hours[i], hours[j], labels[i]))
i = j + 1
else:
i += 1
if windows:
print()
for h0, h1, src in windows:
print(f"Predicted DAM window: {_hhmm(h0)}–{_hhmm(h1)} UT ({src} source, Jupiter up)")
else:
print("\nNo Io-DAM window while Jupiter is up this night.")
Data source: live JPL Horizons Jupiter above the horizon for 14.0 h of this night. CML sweeps 1–358°; Io phase 100–216° while Jupiter is up. Predicted DAM window: 01:20–03:40 UT (Io-B source, Jupiter up) Predicted DAM window: 14:20–14:40 UT (Io-A source, Jupiter up)
Figure 5 — left: Jupiter's altitude through the night (grey = below the horizon). Right: the night's track in the $(\mathrm{CML},\,\text{Io phase})$ plane over the Io-DAM source boxes; filled points are observable (Jupiter up), open points are below the horizon. Where an observable point lands in a box, a storm is likely.
fig, (axL, axR) = plt.subplots(1, 2, figsize=(13, 5.2))
# Left: altitude vs UT, with the below-horizon band shaded.
axL.axhspan(-90, 0, color="0.92", zorder=0)
axL.plot(hours, elev, color="tab:purple", lw=2.2)
axL.scatter(hours[hit], elev[hit], s=42, color="tab:red", zorder=5, label="DAM likely")
axL.axhline(0, color="0.4", lw=1)
axL.set_xlim(0, 24)
axL.set_ylim(-35, 90)
axL.set_xlabel("UT hour (2026-01-11)")
axL.set_ylabel("Jupiter altitude [deg]")
axL.set_title("Is Jupiter up?")
if hit.any():
axL.legend(loc="upper right", fontsize=9)
# Right: (CML, Io phase) plane with the Io-DAM source boxes.
from matplotlib.patches import Rectangle
for name, b in IO_SOURCES.items():
lo, hi = b["cml"]
for x0 in [lo, lo - 360] if hi > 360 else [lo]:
axR.add_patch(
Rectangle(
(x0, b["io"][0]),
(hi - lo),
b["io"][1] - b["io"][0],
facecolor=b["color"],
alpha=0.18,
edgecolor=b["color"],
lw=1.4,
)
)
axR.text(
np.clip((lo + hi) / 2, 15, 345) - (360 if hi > 360 else 0),
(b["io"][0] + b["io"][1]) / 2,
name,
color=b["color"],
fontsize=10,
ha="center",
va="center",
fontweight="bold",
)
axR.scatter(cml[~up], iph[~up], s=26, facecolors="none", edgecolors="0.6", label="below horizon")
axR.scatter(cml[up], iph[up], s=34, color="tab:purple", label="observable")
axR.scatter(cml[hit], iph[hit], s=70, color="tab:red", zorder=5, label="DAM likely")
axR.set_xlim(0, 360)
axR.set_ylim(0, 360)
axR.set_xticks(range(0, 361, 60))
axR.set_yticks(range(0, 361, 60))
axR.set_xlabel("CML (System III) [deg]")
axR.set_ylabel("Io phase [deg]")
axR.set_title("Where the storms live")
axR.legend(loc="upper left", fontsize=8, framealpha=0.9)
fig.tight_layout()
plt.show()
5. Try it yourself¶
These all reuse jansky.solar; copy a cell above and change the marked values.
Exercise 1 — A different CME¶
Re-run the type II analysis for a fast CME, say 1500 km/s, and a slow
one, 350 km/s. For each, (a) make the track with type_ii_track, (b) recover
the speed with shock_speed_from_track, and (c) report the mean drift rate in
kHz/s. Which drifts out of the band first? (Hint: the cell below is a starting
point.)
# TODO: try other speeds, and add your own. Then plot drift rate vs speed.
for v in (350.0, 1500.0):
t_x, f_x = type_ii_track(v, r_start=1.5, duration_s=600.0, harmonic=2)
v_back = shock_speed_from_track(t_x, f_x, harmonic=2)
rate = -(f_x[-1] - f_x[0]) / 600.0 * 1e3
print(
f"V_in {v:6.0f} km/s -> recovered {v_back:6.0f} km/s, "
f"drift {rate:5.2f} kHz/s, band {f_x[0]:.1f}->{f_x[-1]:.1f} MHz"
)
V_in 350 km/s -> recovered 350 km/s, drift 71.98 kHz/s, band 101.3->58.2 MHz V_in 1500 km/s -> recovered 1500 km/s, drift 132.53 kHz/s, band 101.3->21.8 MHz
Solution
The starter cell already does parts (a)-(c); the table from it reads:
| $V_\text{in}$ [km/s] | recovered [km/s] | mean drift [kHz/s] | band [MHz] |
|---|---|---|---|
| 350 | 350 | 71.98 | 101.3 -> 58.2 |
| 1500 | 1500 | 132.53 | 101.3 -> 21.8 |
Two things to read off:
Recovery is exact. Because we analyse with the same density model and harmonic used to build the track (Newkirk
fold=1,harmonic=2),shock_speed_from_trackinverts the forward model self-consistently -- so the recovered speed equals the input to the digit, for both shocks.Faster shock => faster (steeper) drift, and it leaves the band first. Both tracks start at the same frequency (same
r_start=1.5$R_\odot$ => same launch density => same $2f_p \approx 101$ MHz), but the 1500 km/s shock climbs through the density gradient faster, so it sweeps down to ~22 MHz in 10 min while the 350 km/s shock only reaches ~58 MHz. The drift rate is not simply proportional to speed: 1500 km/s is ~4.3x faster than 350 km/s, yet its drift is only ~1.8x steeper -- because the same $\Delta r$ spans a smaller $\Delta f$ once the shock is higher up (the $\sqrt{n_e}$ flattening of Section 2).
Which drifts out of the band first? Take the receiver passband floor used in the dynamic-spectrum cell (20 MHz). Extending the track, the fast shock crosses 20 MHz at $\approx 11$ min; the slow shock never gets below ~58 MHz in 30 min -- it stays in band the whole time. The fast CME exits the bottom first.
import numpy as np
for v in (1500.0, 350.0):
t, f = type_ii_track(v, r_start=1.5, duration_s=1800.0, harmonic=2, n_points=2000)
below = np.where(f < 20.0)[0]
when = f"{t[below[0]]/60:.1f} min" if below.size else "never (stays in band)"
print(f"V={v:6.0f} km/s -> reaches 20 MHz at {when}")
# V= 1500 km/s -> reaches 20 MHz at 11.1 min
# V= 350 km/s -> never (stays in band)
Exercise 2 — Where did the shock start, and how dense?¶
Vary r_start (e.g. 1.2, 1.5, 2.0 $R_\odot$): a shock launched lower starts at
a higher frequency. Then vary the density fold used to generate the burst
while keeping the analysis at fold=1, and quantify the resulting speed bias as a
percentage. How big a fold error gives a 25% speed error?
# TODO: scan r_start and the generate/analyse fold mismatch.
print("r_start scan (V=700, fold=1, harmonic=2):")
for rs in (1.2, 1.5, 2.0):
_, f_rs = type_ii_track(700.0, r_start=rs, harmonic=2)
print(f" r_start={rs} R_sun -> burst spans {f_rs[0]:.1f} -> {f_rs[-1]:.1f} MHz")
print("\\nfold-mismatch scan (true fold varies; analysed with fold=1):")
for true_fold in (1.0, 2.0, 3.0, 4.0):
t_m, f_m = type_ii_track(700.0, r_start=1.5, harmonic=2, fold=true_fold)
v_m = shock_speed_from_track(t_m, f_m, harmonic=2, fold=1.0)
print(f" true fold={true_fold} -> recovered {v_m:6.0f} km/s ({100 * (v_m - 700) / 700:+.0f}%)")
r_start scan (V=700, fold=1, harmonic=2): r_start=1.2 R_sun -> burst spans 232.2 -> 58.0 MHz r_start=1.5 R_sun -> burst spans 101.3 -> 39.1 MHz r_start=2.0 R_sun -> burst spans 44.2 -> 24.9 MHz \nfold-mismatch scan (true fold varies; analysed with fold=1): true fold=1.0 -> recovered 700 km/s (-0%) true fold=2.0 -> recovered 553 km/s (-21%) true fold=3.0 -> recovered 487 km/s (-30%) true fold=4.0 -> recovered 448 km/s (-36%)
Solution
(a) Launch height sets the starting frequency. A shock launched lower starts
in denser plasma (higher $f_p$), so the burst begins at a higher frequency. From
the starter cell ($V=700$ km/s, fold=1, harmonic=2):
| $r_\text{start}$ [$R_\odot$] | start -> end [MHz] |
|---|---|
| 1.2 | 232.2 -> 58.0 |
| 1.5 | 101.3 -> 39.1 |
| 2.0 | 44.2 -> 24.9 |
Dropping the launch from 2.0 to 1.2 $R_\odot$ raises the onset frequency from ~44 MHz to ~232 MHz -- a real type II launched deep in the corona starts well above the metre-wave band and only enters a 20 MHz receiver later in its life.
(b) Density-fold mismatch biases the speed. Generate the burst in a denser
corona (fold $>1$) but analyse with the default fold=1:
true fold |
recovered [km/s] | error |
|---|---|---|
| 1.0 | 700 | 0% |
| 1.5 | 608 | -13% |
| 2.0 | 553 | -21% |
| 3.0 | 487 | -30% |
| 4.0 | 448 | -36% |
The bias is always negative: assuming a thinner corona than reality places
each frequency at too small a height and compresses the inferred climb, so the
fitted slope (the speed) comes out too low. The effect is roughly multiplicative --
for true fold=3 the recovered speed is ~0.70x the truth across a wide range of
input speeds.
How big a fold error gives a 25% speed error? Scanning finely, the -25% point
falls at true fold $\approx 2.4$ (-24% at 2.3, -26% at 2.5):
import numpy as np
for tf in np.arange(2.0, 2.7, 0.1):
t, f = type_ii_track(700.0, r_start=1.5, harmonic=2, fold=tf)
v = shock_speed_from_track(t, f, harmonic=2, fold=1.0)
print(f"true fold={tf:.1f} -> {v:.0f} km/s ({100*(v-700)/700:+.0f}%)")
# true fold=2.4 -> 521 km/s (-26%) (approx -25% near fold ~ 2.4)
So mistaking a streamer ~2.4x denser than the quiet corona for a fold=1 corona
costs you roughly a quarter of the shock speed -- which is exactly why type II
analyses always quote the coronal density model they assumed.
Exercise 3 — Design a Jupiter storm¶
Build your own Jovian spectrogram (copy the Section 4 cells): change the
S-burst drift rate (drift_mhz_per_s, try $-10$ to $-40$ MHz/s), the number
of S-bursts, and the L-burst timings, to mimic an Io-B storm (intense,
S-burst-rich) versus a quieter non-Io event. Add a slowly varying galactic
background by replacing the flat noise floor with one that ramps with time.
# TODO: your storm. Starter: an S-burst-rich "Io-B"-like burst train.
demo = rng(SEED + 7).normal(0.0, 1.0, size=(120, 600))
fdemo = np.linspace(18.0, 28.0, 120)
tdemo = np.linspace(0.0, 60.0, 600)
Fd, Td = np.meshgrid(fdemo, tdemo, indexing="ij")
g = rng(SEED + 8)
for _ in range(80):
ts, fs, life = g.uniform(2, 55), g.uniform(24, 27.5), g.uniform(0.05, 0.3)
f_of_t = fs - 30.0 * (Td - ts)
demo += 7.0 * ((Td >= ts) & (Td <= ts + life)) * np.exp(-0.5 * ((Fd - f_of_t) / 0.3) ** 2)
fig, ax = plt.subplots(figsize=(9, 4))
show_image(
demo, ax=ax, aspect="auto", vmax=6.0, extent=[0, 60, 18, 28], title="Your storm (edit me)"
)
ax.set_xlabel("time [s]")
ax.set_ylabel("frequency [MHz]")
ax.grid(False)
plt.tight_layout()
plt.show()
Solution
This is open-ended; the goal is to feel how the parameters map onto the two
storm archetypes. Reuse the Section 4 idioms (a noisy (freq, time) grid, Gaussian
S-burst streaks, smooth L-burst swells) and change only the marked knobs.
Io-B vs non-Io. An Io-B storm is intense and S-burst-rich: many bursts, high amplitude. A quieter non-Io event has few, weaker S-bursts and leans on smooth L-burst swells. In practice that is two numbers -- the S-burst count and amplitude:
storms = {"Io-B (active)": dict(n_sbursts=300, amp_hi=9.0),
"non-Io (quiet)": dict(n_sbursts=30, amp_hi=4.5)}
S-burst drift rate. Real S-bursts drift down at roughly $-10$ to $-40$ MHz/s;
set drift_mhz_per_s anywhere in that range. Steeper (more negative) values make
the streaks more nearly vertical in a fixed time window -- they sweep the whole band
in a few tens of milliseconds.
A galactic background that ramps with time. Replace the flat noise floor with one whose level rises across the pass (the galactic source drifting through the beam): scale each time column by a ramp, instead of using a single $\sigma$.
import numpy as np
from jansky.signals import rng
f_jove = np.linspace(16.0, 30.0, 280)
t_jove = np.linspace(0.0, 360.0, 1400)
JJ, TT = np.meshgrid(f_jove, t_jove, indexing="ij")
gen = rng(SEED + 11)
bg_level = 1.0 + 1.5 * (t_jove / t_jove[-1]) # rises ~2.5x across the pass
jov = gen.normal(0.0, 1.0, size=JJ.shape) * bg_level[None, :]
# --- Io-B: many strong, fast-drifting S-bursts -----------------------------
g = rng(SEED + 12)
drift_mhz_per_s = -30.0 # steeper = more vertical streaks
for _ in range(300): # non-Io: drop this to ~30
ts = g.uniform(20.0, 340.0)
fs = g.uniform(22.0, 29.0)
life = g.uniform(0.05, 0.4)
amp = g.uniform(6.0, 9.0) # non-Io: lower, e.g. uniform(3.0, 4.5)
f_of_t = fs + drift_mhz_per_s * (TT - ts)
mask = (TT >= ts) & (TT <= ts + life)
jov += amp * mask * np.exp(-0.5 * ((JJ - f_of_t) / 0.35) ** 2)
Plot it exactly as in Section 4 (show_image(jov, ..., vmax=6.0)). You should see:
the Io-B grid is a dense field of crackles brightening toward the right as the
background ramps up; the non-Io grid is mostly the rising galactic floor with a
sparse scatter of faint streaks. Both stay in the same ~18-28 MHz Radio JOVE band --
the difference is entirely in how much S-burst activity Io's flux tube switches
on, which is the physical point of the Io control.
!!! tip "Contribute & compare — Radio JOVE"
The DAM storms and the real .sps file in this chapter come from NASA/SARA Radio JOVE,
the classic ~20 MHz Jupiter-and-Sun project. With a kit (or any 15–30 MHz receiver) you can
record your own storms and submit them to the Radio JOVE archive —
two decades of observer data to compare against, keyed to Io phase and time. Project home:
radiojove.gsfc.nasa.gov; community and grants via the
Society of Amateur Radio Astronomers.
!!! tip "Observe it for real — GAVRT" Jupiter's GHz synchrotron radiation belts (and total solar flux) are exactly what NASA's GAVRT program has students measure on a real 34 m Deep Space Network antenna at Goldstone, operated remotely over the web. The Jupiter Quest campaign compares your radio-belt measurements against Juno; Solar Patrol tracks the solar flux modelled above. Free after a short online training — the no-hardware route to research-grade data.
Recap¶
- The Sun is the sky's brightest radio source — a steady thermal glow (Southworth 1945) plus violent bursts (Hey 1946) — and Jupiter is the only planet you can hear with a wire (Burke & Franklin 1955). Both peak in the decametric band near 20 MHz.
- A plasma is opaque below its plasma frequency, $f_p \propto \sqrt{n_e}$. The corona thins with height, so low frequency ⇒ high in the corona — 20 MHz emission is born a few solar radii out. The Newkirk model ties density, $f_p$, and height together.
- Type II bursts are the radio signature of a CME shock climbing the corona; their slow downward frequency drift plus a density model gives the shock speed. With a self-consistent model the recovery is exact; in practice the assumed density fold and harmonic systematically bias the answer — which is why papers always state the coronal model used. Typical CMEs run from ~250 to 2000+ km/s.
- Jupiter's DAM emission (cyclotron maser, Io-controlled) appears as smooth L-bursts and fast-crackling S-bursts in the same ~18–28 MHz band — the signal a Radio JOVE dual-dipole + SDR records.
- This is real, low-cost science: see the field notes on solar bursts as "low-hanging" amateur work, Projects for the Radio JOVE kit, and the papers timeline for the founding observations.
What's next¶
You now have the science of Radio JOVE — the physics of every waterfall it
records. The natural next step is the data path: the Radio-Sky Spectrograph
(.sps) format the community archives, covered in
Chapter 16: Data Formats & the Ecosystem.
Combine the two and you can take a real archived JOVE spectrogram, find a type II,
and measure a CME's speed yourself.