Four-Spacecraft Wave Dispersion Relation#
author: Louis Richard
Computes the frequency–wavenumber (k–ω) power spectrum of a perpendicular component of the magnetic field fluctuations (SCM burst) in field-aligned coordinates (FAC) with
mms.fk_power_spectrum_4sc, using the positions (MEC) and the background magnetic field (FGM burst) of the four spacecraft. Plots the power as a function of frequency and of \(|\mathbf{k}|\), \(k_x\), \(k_y\) and \(k_z\), and in four wave vector planes, with the mean wave vector and frequency and their
standard deviations.[1]:
%matplotlib inline
import matplotlib.pyplot as plt
from pyrfu import mms, pyrf
from pyrfu.plot import use_pyrfu_style, plot_spectr
use_pyrfu_style(usetex=False)
Define time interval and data path#
[2]:
mms.db_init(default="local", local="/data/mms")
tint = ["2015-10-16T13:05:24.000", "2015-10-16T13:05:50.000"]
[05-Oct-26 17:38:57] INFO: Updating MMS data access configuration in /homelocal/louisr/.config/pyrfu/mms_config.json...
[05-Oct-26 17:38:57] WARNING: No system keyring available: MMS SDC credentials saved in plain text in /homelocal/louisr/.local/share/python_keyring/keyring_pass.cfg
Load data#
Load FGM and SCM magnetic fields#
[3]:
b_xyz_mms = [mms.get_data("b_gse_fgm_brst_l2", tint, ic) for ic in range(1, 5)]
b_scm_mms = [mms.get_data("b_gse_scm_brst_l2", tint, ic) for ic in range(1, 5)]
[05-Oct-26 17:38:57] INFO: Loading mms1_fgm_b_gse_brst_l2...
[05-Oct-26 17:38:57] INFO: Loading mms2_fgm_b_gse_brst_l2...
[05-Oct-26 17:38:57] INFO: Loading mms3_fgm_b_gse_brst_l2...
[05-Oct-26 17:38:57] INFO: Loading mms4_fgm_b_gse_brst_l2...
[05-Oct-26 17:38:57] INFO: Loading mms1_scm_acb_gse_scb_brst_l2...
[05-Oct-26 17:38:57] INFO: Loading mms2_scm_acb_gse_scb_brst_l2...
[05-Oct-26 17:38:57] INFO: Loading mms3_scm_acb_gse_scb_brst_l2...
[05-Oct-26 17:38:57] INFO: Loading mms4_scm_acb_gse_scb_brst_l2...
Load spacecraft positions#
[4]:
tint_long = pyrf.extend_tint(tint, [-60, 60])
r_xyz_mms = [mms.get_data("r_gse_mec_srvy_l2", tint_long, ic) for ic in range(1, 5)]
[05-Oct-26 17:38:57] INFO: Loading mms1_mec_r_gse...
[05-Oct-26 17:38:57] INFO: Loading mms2_mec_r_gse...
[05-Oct-26 17:38:57] INFO: Loading mms3_mec_r_gse...
[05-Oct-26 17:38:57] INFO: Loading mms4_mec_r_gse...
Convert SCM magnetic field to field-aligned coordinates#
[5]:
b_xyz_avg = pyrf.avg_4sc(b_xyz_mms)
b_fac_mms = [pyrf.convert_fac(b_scm, b_xyz_avg, [1, 0, 0]) for b_scm in b_scm_mms]
Compute four-spacecraft dispersion relation#
Define time interval of the wave burst#
[6]:
tints = ["2015-10-16T13:05:26.500", "2015-10-16T13:05:27.000"]
Select a perpendicular component of the magnetic field fluctuations#
[7]:
b_para_mms = [b_fac[:, 0] for b_fac in b_fac_mms]
Compute the k–ω power spectrum#
[8]:
f_range = [300.0, 1000.0]
power = mms.fk_power_spectrum_4sc(
b_para_mms,
r_xyz_mms,
b_xyz_mms,
tints,
df=10.0,
num_k=500,
cav=4,
w_width=2.0,
f_range=f_range,
)
/homelocal/louisr/irfu-python/.venv/lib/python3.12/site-packages/pyrfu/mms/fk_power_spectrum_4sc.py:302: ComplexWarning: Casting complex values to real discards the imaginary part
fk_power += (w[i].data * np.conj(w[i].data) / 4).astype(np.float64)
[05-Oct-26 17:39:15] INFO: Computing power versus (kx,f); (ky,f), (kz,f), (k,f)
[05-Oct-26 17:39:15] INFO: Computing power versus (kx,ky); (kx,kz); (ky,kz)
[05-Oct-26 17:39:15] INFO: Computing power versus kperp,kpara
Change units of wave vector and frequency for plotting#
[9]:
# Convert wave vector and frequency of the power spectrum to km^-1 and kHz
power = power.assign_coords(
k_mag=power.k_mag.data * 1e3,
k_x=power.k_x.data * 1e3,
k_y=power.k_y.data * 1e3,
k_z=power.k_z.data * 1e3,
k_perp=power.k_perp.data * 1e3,
k_para=power.k_para.data * 1e3,
)
power = power.assign_coords(f=power.f.data * 1e-3)
# Convert wave vector and frequency of the peak power
power.attrs = {
key: value * 1e3 if "k_" in key else value * 1e-3
for key, value in power.attrs.items()
}
Plot#
[10]:
f, axs = plt.subplots(2, 4, figsize=(11, 5))
f.subplots_adjust(wspace=0.6, hspace=0.9)
f_lim = [
0,
]
axs[0, 0], cax00 = plot_spectr(
axs[0, 0], power.k_mag_f, cscale="log", cmap="jet", clim=[1e-6, 1e0], colorbar="top"
)
axs[0, 0].axhline(f_range[0], color="gray", linestyle="--")
axs[0, 0].axhline(f_range[1], color="gray", linestyle="--")
axs[0, 0].errorbar(
power.attrs["k_mag_avg"],
power.attrs["f_avg"],
xerr=power.attrs["k_mag_std"],
yerr=power.attrs["f_std"],
fmt="o",
color="black",
markersize=4,
)
axs[0, 0].set_xlim([0, 0.3])
axs[0, 0].set_ylim([0, 2])
axs[0, 0].set_xlabel(r"$k~[\mathrm{km}^{-1}]$")
axs[0, 0].set_ylabel(r"$f~[\mathrm{kHz}]$")
cax00.set_xlabel(r"$P(k, f)/P_{max}$")
axs[0, 1], cax01 = plot_spectr(
axs[0, 1], power.k_x_f, cscale="log", cmap="jet", clim=[1e-6, 1e0], colorbar="top"
)
axs[0, 1].axhline(f_range[0], color="gray", linestyle="--")
axs[0, 1].axhline(f_range[1], color="gray", linestyle="--")
axs[0, 1].errorbar(
power.attrs["k_x_avg"],
power.attrs["f_avg"],
xerr=power.attrs["k_x_std"],
yerr=power.attrs["f_std"],
fmt="o",
color="black",
markersize=4,
)
axs[0, 1].set_xlim([-0.3, 0.3])
axs[0, 1].set_ylim([0, 2])
axs[0, 1].set_xlabel(r"$k_x~[\mathrm{km}^{-1}]$")
axs[0, 1].set_ylabel(r"$f~[\mathrm{kHz}]$")
cax01.set_xlabel(r"$P(k_x, f)/P_{max}$")
axs[0, 2], cax02 = plot_spectr(
axs[0, 2], power.k_y_f, cscale="log", cmap="jet", clim=[1e-6, 1e0], colorbar="top"
)
axs[0, 2].axhline(f_range[0], color="gray", linestyle="--")
axs[0, 2].axhline(f_range[1], color="gray", linestyle="--")
axs[0, 2].errorbar(
power.attrs["k_y_avg"],
power.attrs["f_avg"],
xerr=power.attrs["k_y_std"],
yerr=power.attrs["f_std"],
fmt="o",
color="black",
markersize=4,
)
axs[0, 2].set_xlim([-0.3, 0.3])
axs[0, 2].set_ylim([0, 2])
axs[0, 2].set_xlabel(r"$k_y~[\mathrm{km}^{-1}]$")
axs[0, 2].set_ylabel(r"$f~[\mathrm{kHz}]$")
cax02.set_xlabel(r"$P(k_y, f)/P_{max}$")
axs[0, 3], cax03 = plot_spectr(
axs[0, 3], power.k_z_f, cscale="log", cmap="jet", clim=[1e-6, 1e0], colorbar="top"
)
axs[0, 3].axhline(f_range[0], color="gray", linestyle="--")
axs[0, 3].axhline(f_range[1], color="gray", linestyle="--")
axs[0, 3].errorbar(
power.attrs["k_z_avg"],
power.attrs["f_avg"],
xerr=power.attrs["k_z_std"],
yerr=power.attrs["f_std"],
fmt="o",
color="black",
markersize=4,
)
axs[0, 3].set_xlim([-0.3, 0.3])
axs[0, 3].set_ylim([0, 2])
axs[0, 3].set_xlabel(r"$k_z~[\mathrm{km}^{-1}]$")
axs[0, 3].set_ylabel(r"$f~[\mathrm{kHz}]$")
cax03.set_xlabel(r"$P(k_z, f)/P_{max}$")
axs[1, 0], cax10 = plot_spectr(
axs[1, 0],
power.k_perp_k_para,
cscale="log",
cmap="jet",
clim=[1e-6, 1e0],
colorbar="top",
)
axs[1, 0].errorbar(
power.attrs["k_perp_avg"],
power.attrs["k_para_avg"],
xerr=power.attrs["k_perp_std"],
yerr=power.attrs["k_para_std"],
fmt="o",
color="black",
markersize=4,
)
axs[1, 0].set_xlim([0, 0.3])
axs[1, 0].set_ylim([-0.3, 0.3])
axs[1, 0].set_xlabel(r"$k_\bot~[\mathrm{km}^{-1}]$")
axs[1, 0].set_ylabel(r"$k_\parallel~[\mathrm{km}^{-1}]$")
cax10.set_xlabel(r"$P(k_\bot, k_\parallel)/P_{max}$")
axs[1, 1], cax11 = plot_spectr(
axs[1, 1], power.k_x_k_y, cscale="log", cmap="jet", clim=[1e-6, 1e0], colorbar="top"
)
axs[1, 1].errorbar(
power.attrs["k_x_avg"],
power.attrs["k_y_avg"],
xerr=power.attrs["k_x_std"],
yerr=power.attrs["k_y_std"],
fmt="o",
color="black",
markersize=4,
)
axs[1, 1].set_xlim([-0.3, 0.3])
axs[1, 1].set_ylim([-0.3, 0.3])
axs[1, 1].set_xlabel(r"$k_x~[\mathrm{km}^{-1}]$")
axs[1, 1].set_ylabel(r"$k_y~[\mathrm{km}^{-1}]$")
cax11.set_xlabel(r"$P(k_x, k_y)/P_{max}$")
axs[1, 2], cax12 = plot_spectr(
axs[1, 2], power.k_x_k_z, cscale="log", cmap="jet", clim=[1e-6, 1e0], colorbar="top"
)
axs[1, 2].errorbar(
power.attrs["k_x_avg"],
power.attrs["k_z_avg"],
xerr=power.attrs["k_x_std"],
yerr=power.attrs["k_z_std"],
fmt="o",
color="black",
markersize=4,
)
axs[1, 2].set_xlim([-0.3, 0.3])
axs[1, 2].set_ylim([-0.3, 0.3])
axs[1, 2].set_xlabel(r"$k_x~[\mathrm{km}^{-1}]$")
axs[1, 2].set_ylabel(r"$k_z~[\mathrm{km}^{-1}]$")
cax12.set_xlabel(r"$P(k_x, k_z)/P_{max}$")
axs[1, 3], cax13 = plot_spectr(
axs[1, 3], power.k_y_k_z, cscale="log", cmap="jet", clim=[1e-6, 1e0], colorbar="top"
)
axs[1, 3].errorbar(
power.attrs["k_y_avg"],
power.attrs["k_z_avg"],
xerr=power.attrs["k_y_std"],
yerr=power.attrs["k_z_std"],
fmt="o",
color="black",
markersize=4,
)
axs[1, 3].set_xlim([-0.3, 0.3])
axs[1, 3].set_ylim([-0.3, 0.3])
axs[1, 3].set_xlabel(r"$k_y~[\mathrm{km}^{-1}]$")
axs[1, 3].set_ylabel(r"$k_z~[\mathrm{km}^{-1}]$")
cax13.set_xlabel(r"$P(k_y, k_z)/P_{max}$")
[10]:
Text(0.5, 0, '$P(k_y, k_z)/P_{max}$')
[ ]: