Reduced Ion Distributions in User-Defined Coordinates#

author: Louis Richard

Computes reduced ion velocity distribution functions (VDFs) from FPI burst-mode skymaps with Monte-Carlo integration (mms.reduce), in the user-defined coordinates \(\mathbf{\hat{n}}\), \(\mathbf{\hat{t}_1}\), \(\mathbf{\hat{t}_2}\) built from a shock normal and the upstream magnetic field. It plots the 1D reduced distributions along each axis as time series, and the 2D reduced distributions in the three planes averaged over 2 s.

[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)

Load data#

Define time interval and spacecraft index#

[2]:

mms.db_init(default="local", local="/data/mms")
tint = ["2015-12-28T03:57:10", "2015-12-28T03:59:00"]
mms_id = 1

[05-Oct-26 18:04:49] INFO: Updating MMS data access configuration in /homelocal/louisr/.config/pyrfu/mms_config.json...
[05-Oct-26 18:04:49] WARNING: No system keyring available: MMS SDC credentials saved in plain text in /homelocal/louisr/.local/share/python_keyring/keyring_pass.cfg

[3]:

b_dmpa = mms.get_data("b_dmpa_fgm_brst_l2", tint, mms_id)

[05-Oct-26 18:04:49] INFO: Loading mms1_fgm_b_dmpa_brst_l2...

Load the spacecraft attitude#

Used to convert between Geocentric Solar Ecliptic (GSE) and spacecraft coordinates.

[4]:

defatt = mms.load_ancillary("defatt", tint, mms_id)

[05-Oct-26 18:04:50] INFO: Loading ancillary defatt files...

Load the ion VDF skymaps and their uncertainties#

The uncertainty is \(\delta f / f = 1/\sqrt{n}\), with \(n\) the number of counts.

[5]:

vdf_i = mms.get_data("pdi_fpi_brst_l2", tint, mms_id)
vdf_i_err = mms.get_data("pderri_fpi_brst_l2", tint, mms_id)

[05-Oct-26 18:04:53] INFO: Loading mms1_dis_dist_brst...
[05-Oct-26 18:04:53] INFO: Loading mms1_dis_disterr_brst...

Remove the one-count level#

Phase-space density values below 1.1 times their uncertainty (about one count) are set to zero. This also makes mms.reduce faster.

[6]:

vdf_i.data.data[vdf_i.data.data < 1.1 * vdf_i_err.data.data] = 0.0

Define the coordinate system of the projection#

Shock normal vector in GSE#

Get it from irf_shock_normal or irf_shock_gui.

[7]:

n_vec = np.array([0.9580, -0.2708, -0.0938])
n_vec /= np.linalg.norm(n_vec)

Upstream magnetic field in GSE#

[8]:

b_u = [-1.0948, -2.6270, 1.6478]

\(\mathbf{\hat{t}_2}\) vector in GSE#

Same vectors as in Johlander et al., 2016, PRL.

[9]:

t2_vec = np.cross(n_vec, b_u) / np.linalg.norm(np.cross(n_vec, b_u))

\(\mathbf{\hat{t}_1}\) vector in GSE#

[10]:

t1_vec = np.cross(t2_vec, n_vec)

Build the time series of \(\mathbf{\hat{n}}\), \(\mathbf{\hat{t}_1}\) and \(\mathbf{\hat{t}_2}\)#

[11]:

n_t = len(b_dmpa.time.data)
n_gse = pyrf.ts_vec_xyz(b_dmpa.time.data, np.tile(n_vec[np.newaxis, :], [n_t, 1]))
t1_gse = pyrf.ts_vec_xyz(b_dmpa.time.data, np.tile(t1_vec[np.newaxis, :], [n_t, 1]))
t2_gse = pyrf.ts_vec_xyz(b_dmpa.time.data, np.tile(t2_vec[np.newaxis, :], [n_t, 1]))

Transform the vectors to spacecraft coordinates with the spacecraft attitude#

[12]:

n_dmpa = mms.dsl2gse(n_gse, defatt, -1)
t1_dmpa = mms.dsl2gse(t1_gse, defatt, -1)
t2_dmpa = mms.dsl2gse(t2_gse, defatt, -1)

Build the transformation matrices from spacecraft coordinates#

[13]:

nt1t2 = np.transpose(np.stack([n_dmpa.data, t1_dmpa.data, t2_dmpa.data]), [1, 2, 0])
nt1t2 = pyrf.ts_tensor_xyz(b_dmpa.time.data, nt1t2)

t1t2n = np.transpose(np.stack([t1_dmpa.data, t2_dmpa.data, n_dmpa.data]), [1, 2, 0])
t1t2n = pyrf.ts_tensor_xyz(b_dmpa.time.data, t1t2n)

t2nt1 = np.transpose(np.stack([t2_dmpa.data, n_dmpa.data, t1_dmpa.data]), [1, 2, 0])
t2nt1 = pyrf.ts_tensor_xyz(b_dmpa.time.data, t2nt1)

Define the velocity grids along the three vectors#

[14]:

vn_lim = np.array([-800.0, 800.0], dtype=np.float64)
vt1_lim = vn_lim
vt2_lim = vn_lim + 300.0

vg_1d_n = 1e3 * np.linspace(vn_lim[0], vn_lim[1], 100)
vg_1d_t1 = 1e3 * np.linspace(vt1_lim[0], vt1_lim[1], 100)
vg_1d_t2 = 1e3 * np.linspace(vt2_lim[0], vt2_lim[1], 100)

Reduce the distribution along the three vectors#

[15]:

n_mc = 200

f1dn = mms.reduce(vdf_i, projection_dim="1d", xyz=nt1t2, n_mc=n_mc, vg=vg_1d_n)
f1dt1 = mms.reduce(vdf_i, projection_dim="1d", xyz=t1t2n, n_mc=n_mc, vg=vg_1d_t1)
f1dt2 = mms.reduce(vdf_i, projection_dim="1d", xyz=t2nt1, n_mc=n_mc, vg=vg_1d_t2)

  0%|                               | 0/733 [00:00<?, ?it/s]100%|████████████████████| 733/733 [00:06<00:00, 119.55it/s]
100%|████████████████████| 733/733 [00:04<00:00, 149.82it/s]
100%|████████████████████| 733/733 [00:04<00:00, 152.24it/s]

Plot#

[16]:

f, axs = plt.subplots(4, sharex="all", figsize=(9, 6.5))
f.subplots_adjust(hspace=0.0, left=0.15, right=0.85, bottom=0.08, top=0.95)

plot_line(axs[0], b_dmpa)
plot_line(axs[0], pyrf.norm(b_dmpa), color="k")
axs[0].legend([r"$B_{x}$", r"$B_{y}$", r"$B_{z}$", r"$|\mathbf{B}|$"], ncol=4)
axs[0].set_ylabel("$B_{dmpa}~[\mathrm{nT}]$")

axs[1], cax1 = plot_spectr(
    axs[1], f1dn, cscale="log", clim=[1e-4, 1e3], cmap="Spectral_r"
)
axs[1].set_ylim(vn_lim)
axs[1].set_ylabel("$v_n~[\\mathrm{km}~\\mathrm{s}^{-1}]$")
axs[1].set_yticks([-500.0, 0.0, 500.0])

axs[2] = plot_spectr(
    axs[2], f1dt1, cscale="log", clim=[1e-4, 1e3], colorbar="none", cmap="Spectral_r"
)
axs[2].set_ylim(vt1_lim)
axs[2].set_ylabel("$v_{t1}~[\\mathrm{km}~\\mathrm{s}^{-1}]$")
axs[2].set_yticks([-500.0, 0.0, 500.0])

axs[3] = plot_spectr(
    axs[3], f1dt2, cscale="log", clim=[1e-4, 1e3], colorbar="none", cmap="Spectral_r"
)
axs[3].set_ylim(vt2_lim)
axs[3].set_ylabel("$v_{t2}~[\\mathrm{km}~\\mathrm{s}^{-1}]$")
axs[3].set_yticks([0.0, 500.0, 1000.0])

pos_axs3 = axs[3].get_position()
pos_cax1 = cax1.get_position()
x0 = pos_cax1.x0
y0 = pos_axs3.y0
width = pos_cax1.width
height = pos_cax1.y0 + pos_cax1.height - y0
cax1.set_position([x0, y0, width, height])

cax1.set_ylabel("$F_i~[\\mathrm{s}~\\mathrm{m}^{-4}]$")
f.align_ylabels(axs)

<>:7: SyntaxWarning: invalid escape sequence '\m'
<>:7: SyntaxWarning: invalid escape sequence '\m'
/tmp/ipykernel_517108/2743122119.py:7: SyntaxWarning: invalid escape sequence '\m'
  axs[0].set_ylabel("$B_{dmpa}~[\mathrm{nT}]$")

../../_images/examples_01_mms_example_mms_reduced_ion_dist_33_1.png

2D projection of the ion VDF#

Define the averaging interval: center time \(\pm 1~\mathrm{s}\)#

[17]:

t_2d = np.datetime64("2015-12-28T03:57:40.300")
t_2d = pyrf.extend_tint([t_2d, t_2d], [-1, 1])

Define the velocity grid for the 2D projections#

[18]:

vg_2d = np.linspace(-1500, 1500, 200) * 1e3

Reduce the ion distributions in the \((n, t_1)\), \((t_2, n)\) and \((t_1, t_2)\) planes#

[19]:

f2Dnt1 = mms.reduce(
    pyrf.time_clip(vdf_i, t_2d),
    xyz=nt1t2,
    dim="2d",
    base="cart",
    n_mc=n_mc * 5,
    vg=vg_2d,
)
f2Dt2n = mms.reduce(
    pyrf.time_clip(vdf_i, t_2d),
    xyz=t2nt1,
    dim="2d",
    base="cart",
    n_mc=n_mc * 5,
    vg=vg_2d,
)
f2Dt1t2 = mms.reduce(
    pyrf.time_clip(vdf_i, t_2d),
    xyz=t1t2n,
    dim="2d",
    base="cart",
    n_mc=n_mc * 5,
    vg=vg_2d,
)

  0%|                                | 0/14 [00:00<?, ?it/s]100%|███████████████████████| 14/14 [00:00<00:00, 62.62it/s]
100%|███████████████████████| 14/14 [00:00<00:00, 63.90it/s]
100%|███████████████████████| 14/14 [00:00<00:00, 66.18it/s]

Plot the 2D reduced distributions#

[20]:

f, axs = plt.subplots(1, 3, figsize=(7, 3.0))
f.subplots_adjust(wspace=0.7, left=0.11, right=0.97, bottom=0.2, top=0.78)
axs[0], cax0 = plot_spectr(
    axs[0], f2Dnt1.mean(axis=0), cscale="log", colorbar="top", cmap="Spectral_r"
)
axs[0].set_xlim([-1e3, 1e3])
axs[0].set_ylim([-1e3, 1e3])
axs[0].set_aspect("equal")
axs[0].set_xlabel("$v_{n}~[\\mathrm{km}~\\mathrm{s}^{-1}]$")
axs[0].set_ylabel("$v_{t1}~[\\mathrm{km}~\\mathrm{s}^{-1}]$")

axs[1] = plot_spectr(
    axs[1], f2Dt2n.mean(axis=0), cscale="log", colorbar="none", cmap="Spectral_r"
)
axs[1].set_xlim([-0.5e3, 1.5e3])
axs[1].set_ylim([-1e3, 1e3])
axs[1].set_aspect("equal")
axs[1].set_xlabel("$v_{t2}~[\\mathrm{km}~\\mathrm{s}^{-1}]$")
axs[1].set_ylabel("$v_{n}~[\\mathrm{km}~\\mathrm{s}^{-1}]$")

axs[2] = plot_spectr(
    axs[2], f2Dt1t2.mean(axis=0), cscale="log", colorbar="none", cmap="Spectral_r"
)
axs[2].set_xlim([-1e3, 1e3])
axs[2].set_ylim([-0.5e3, 1.5e3])
axs[2].set_aspect("equal")
axs[2].set_xlabel("$v_{t1}~[\\mathrm{km}~\\mathrm{s}^{-1}]$")
axs[2].set_ylabel("$v_{t2}~[\\mathrm{km}~\\mathrm{s}^{-1}]$")

pos_axs2 = axs[2].get_position()
pos_cax0 = cax0.get_position()
x0 = pos_cax0.x0
y0 = pos_cax0.y0 - 0.01
width = pos_axs2.x0 + pos_axs2.width - pos_cax0.x0
height = 0.02
cax0.set_position([x0, y0, width, height])

cax0.set_xlabel("$F_i~[\\mathrm{s}^2~\\mathrm{m}^{-5}]$")

[20]:

Text(0.5, 0, '$F_i~[\\mathrm{s}^2~\\mathrm{m}^{-5}]$')

../../_images/examples_01_mms_example_mms_reduced_ion_dist_43_1.png

[ ]:

[ ]: