Interplanetary Shocks in the Normal Incidence Frame#
author: Louis Richard
Computes the normal and parameters of an interplanetary (IP) shock from upstream and downstream averages of MMS burst-mode FGM and FPI ion data. It then transforms the magnetic field, the ion bulk velocity and the ion velocity distribution reduced along the shock normal to the normal incidence frame (NIF). The same procedure can be used at the bow shock.
[1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from pyrfu import mms, pyrf
from pyrfu.plot import use_pyrfu_style, plot_spectr, plot_line
use_pyrfu_style(usetex=False)
Set the path to the MMS data.
[2]:
mms.db_init(default="local", local="/data/mms")
[05-Oct-26 17:58:30] INFO: Updating MMS data access configuration in /homelocal/louisr/.config/pyrfu/mms_config.json...
[05-Oct-26 17:58:30] WARNING: No system keyring available: MMS SDC credentials saved in plain text in /homelocal/louisr/.local/share/python_keyring/keyring_pass.cfg
Non-exhaustive list of IP shocks observed by MMS in burst mode#
[3]:
ipshock_events = [
{
"tint": ["2017-10-24T08:25:45.000000", "2017-10-24T08:27:45.000000"],
"tintu": ["2017-10-24T08:26:36.000000", "2017-10-24T08:26:42.000000"],
"tintd": ["2017-10-24T08:26:54.000000", "2017-10-24T08:27:04.000000"],
},
{
"tint": ["2018-01-08T06:40:45.000000", "2018-01-08T06:41:30.000000"],
"tintu": ["2018-01-08T06:40:45.000000", "2018-01-08T06:41:00.000000"],
"tintd": ["2018-01-08T06:41:20.000000", "2018-01-08T06:41:30.000000"],
},
{
"tint": ["2022-02-01T22:17:00.000000", "2022-02-01T22:20:00.000000"],
"tintu": ["2022-02-01T22:18:08.000000", "2022-02-01T22:18:26.000000"],
"tintd": ["2022-02-01T22:18:43.000000", "2022-02-01T22:18:56.000000"],
},
{
"tint": ["2022-02-11T10:23:20.000000", "2022-02-11T10:25:10.000000"],
"tintu": ["2022-02-11T10:24:05.000000", "2022-02-11T10:24:07.000000"],
"tintd": ["2022-02-11T10:24:16.000000", "2022-02-11T10:24:22.000000"],
},
]
Select an IP shock#
[4]:
event_num = 0
reference_frame = "nif"
mms_id = 2
[5]:
tint = ipshock_events[event_num]["tint"]
tintu = ipshock_events[event_num]["tintu"]
tintd = ipshock_events[event_num]["tintd"]
Load data#
[6]:
# Spacecraft location
r_gse = mms.get_data("r_gse_mec_srvy_l2", tint, mms_id)
# Magnetic field
b_gse = mms.get_data("b_gse_fgm_brst_l2", tint, mms_id)
# Ion number density and bulk velocity
n_i = mms.get_data("ni_fpi_brst_l2", tint, mms_id)
v_gse_i = mms.get_data("vi_gse_fpi_brst_l2", tint, mms_id)
# Electron number density
n_e = mms.get_data("ne_fpi_brst_l2", tint, mms_id)
vdf_i = mms.get_data("pdi_fpi_brst_l2", tint, mms_id)
vdf_err_i = mms.get_data("pderri_fpi_brst_l2", tint, mms_id)
vdf_i.data.data[vdf_i.data.data < 1.1 * vdf_err_i.data.data] = 0.0
[05-Oct-26 17:58:30] INFO: Loading mms2_mec_r_gse...
[05-Oct-26 17:58:30] INFO: Loading mms2_fgm_b_gse_brst_l2...
[05-Oct-26 17:58:30] INFO: Loading mms2_dis_numberdensity_brst...
[05-Oct-26 17:58:31] INFO: Loading mms2_dis_bulkv_gse_brst...
[05-Oct-26 17:58:31] INFO: Loading mms2_des_numberdensity_brst...
[05-Oct-26 17:58:31] INFO: Loading mms2_dis_dist_brst...
[05-Oct-26 17:58:31] INFO: Loading mms2_dis_disterr_brst...
Compute the shock normal and shock parameters (NIF only)#
[7]:
if reference_frame.lower() == "nif":
plp = {}
plp["b_u"] = np.mean(pyrf.time_clip(b_gse, tintu).data, axis=0)
plp["n_u"] = np.mean(pyrf.time_clip(n_i, tintu).data, axis=0)
plp["v_u"] = np.mean(pyrf.time_clip(v_gse_i, tintu).data, axis=0)
plp["b_d"] = np.mean(pyrf.time_clip(b_gse, tintd).data, axis=0)
plp["n_d"] = np.mean(pyrf.time_clip(n_i, tintd).data, axis=0)
plp["v_d"] = np.mean(pyrf.time_clip(v_gse_i, tintd).data, axis=0)
plp["r_xyz"] = np.mean(r_gse.data, axis=0)
# Shock normal and shock speed
nst = pyrf.shock_normal(plp)
# Use the mixed method 3 normal
nvec = nst["n"]["mx_3"]
# and the SB shock speed
v_sh = nst["v_sh"]["sb"]["mx_3"]
# Add them to plp
plp["nvec"] = nvec
plp["v_sh"] = v_sh
plp["ref_sys"] = reference_frame.lower()
# Shock coordinate system (n, t1, t2)
t2vec = np.cross(nvec, plp["b_u"])
t2vec /= np.linalg.norm(t2vec)
t1vec = np.cross(t2vec, nvec)
t1vec /= np.linalg.norm(t1vec)
t1vec = np.tile(t1vec, (len(vdf_i.time), 1))
t2vec = np.tile(t2vec, (len(vdf_i.time), 1))
nvec = np.tile(nvec, (len(vdf_i.time), 1))
# Rotation matrix to the shock coordinate system
r_mat = np.transpose(np.stack([nvec, t1vec, t2vec]), [1, 2, 0])
r_mat = pyrf.ts_tensor_xyz(vdf_i.time.data, r_mat)
# Other shock parameters
shp = pyrf.shock_parameters(plp)
# Velocity of the NIF in the spacecraft frame
v_nif = shp["v_nif_u"]
elif reference_frame.lower() == "s/c":
r_mat = np.tile(np.eye(3), (len(vdf_i.time), 1, 1))
r_mat = pyrf.ts_tensor_xyz(vdf_i.time.data, r_mat)
v_nif = np.zeros(3)
nvec = r_mat.data[:, :, 0]
Reduce the ion VDF along the shock normal (NIF) or x GSE#
[8]:
# Monte Carlo iterations per bin: decrease for speed, increase for a smoother plot.
n_mc = 2e2
if reference_frame.lower() == "nif":
# Velocity grid
v_lim = [0.0, 1000.0] # km/s
v_1d = np.linspace(v_lim[0], v_lim[1], 100) * 1e3
# Reduce the ion distribution
f1d = mms.reduce(vdf_i, projection_dim="1d", xyz=r_mat, n_mc=n_mc, vg=v_1d)
f1d = f1d.assign_coords(vx=f1d.vx.data - np.dot(v_nif, nvec[0, :]) / 1e3)
elif reference_frame.lower() == "s/c":
# Velocity grid
v_lim = [-1000.0, 0.0] # km/s
v_1d = np.linspace(v_lim[0], v_lim[1], 100) * 1e3
# Reduce the ion distribution
f1d = mms.reduce(vdf_i, projection_dim="1d", xyz=r_mat, n_mc=n_mc, vg=v_1d)
f1d = f1d.assign_coords(vx=f1d.vx.data - np.dot(v_nif, nvec[0, :]) / 1e3)
0%| | 0/800 [00:00<?, ?it/s]100%|████████████████████| 800/800 [00:05<00:00, 139.93it/s]
Plot#
[9]:
legend_options = dict(
handlelength=1, ncol=1, frameon=False, loc="upper left", bbox_to_anchor=(1.0, 1.0)
)
[10]:
f, axs = plt.subplots(5, sharex="all", figsize=(6, 8))
f.subplots_adjust(hspace=0, left=0.1, right=0.82, bottom=0.05, top=0.95)
plot_line(axs[0], pyrf.norm(b_gse), color="k")
axs[0].set_ylabel("$|\\mathbf{B}|~[\\mathrm{nT}]$")
plot_line(axs[1], pyrf.new_xyz(b_gse, r_mat.data[0, ...]))
axs[1].axhline(0.0, color="k", linestyle="--")
axs[1].legend(["$B_{n}$", "$B_{t1}$", "$B_{t2}$"], **legend_options)
axs[1].set_ylabel("$\\mathbf{B}~[\\mathrm{nT}]$")
plot_line(axs[2], n_i, label="$n_i$")
plot_line(axs[2], n_e, label="$n_e$")
axs[2].legend(**legend_options)
axs[2].set_ylabel("$n~[\\mathrm{cm}^{-3}]$")
plot_line(axs[3], pyrf.new_xyz(v_gse_i - v_nif / 1e3, r_mat.data[0, ...]))
axs[3].legend(["$V_{n}$", "$V_{t2}$", "$V_{t2}$"], **legend_options)
axs[3].set_ylabel("$\\mathbf{V}_i~[\\mathrm{km} \\mathrm{s}^{-1}]$")
axs[4], caxs4 = plot_spectr(axs[4], f1d, cscale="log", cmap="Spectral_r")
axs[4].set_ylabel("$v_n~[\\mathrm{km} \\mathrm{s}^{-1}]$")
caxs4.set_ylabel("$F_i~[\\mathrm{s}~\\mathrm{m}^{-4}]$")
f.align_ylabels(axs)
f.suptitle(f"MMS {mms_id} - IP shock in the {reference_frame.upper()} frame")
[10]:
Text(0.5, 0.98, 'MMS 2 - IP shock in the NIF frame')
