Particle Distributions: Units, Rebinning, Projections

author: Louis Richard

Shows basic operations on FPI burst-mode electron velocity distribution functions (VDFs): omnidirectional phase-space density (PSD, mms.vdf_omni), pitch-angle distributions (mms.get_pitch_angle_dist), energy limits (mms.vdf_elim), conversions to differential energy flux (mms.psd2def) and differential particle flux (mms.psd2dpf), and rebinning to 64 energy channels (mms.vdf_to_e64). It plots pitch-angle spectrograms in two energy ranges and energy spectrograms in three pitch-angle ranges, and projections of the VDF onto the planes defined by E, E\(\times\)B and B with mms.vdf_projection, using the FGM magnetic field, EDP electric field and spacecraft potential.

%matplotlib inline
import numpy as np
import xarray as xr
import matplotlib as mpl
import matplotlib.pyplot as plt

from pyrfu import mms, pyrf
from pyrfu.plot import (
    use_pyrfu_style,
    plot_line,
    plot_spectr,
    plot_projection,
    make_labels,
)

use_pyrfu_style(usetex=False)

Define spacecraft index, time interval and data path

mms.db_init(default="local", local="/data/mms")
tint = ["2015-12-02T01:14:15.000", "2015-12-02T01:15:13.000"]
mms_id = 3

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

Load the velocity distribution functions (VDFs)

vdf_i = mms.get_data("pdi_fpi_brst_l2", tint, mms_id)
vdf_e = mms.get_data("pde_fpi_brst_l2", tint, mms_id)

[05-Oct-26 17:59:33] INFO: Loading mms3_dis_dist_brst...
[05-Oct-26 17:59:34] INFO: Loading mms3_des_dist_brst...

Load supporting data

b_xyz = mms.get_data("b_dmpa_fgm_brst_l2", tint, mms_id)
e_xyz = mms.get_data("e_dsl_edp_brst_l2", tint, mms_id)
sc_pot = mms.get_data("v_edp_brst_l2", tint, mms_id)

[05-Oct-26 17:59:35] INFO: Loading mms3_fgm_b_dmpa_brst_l2...
[05-Oct-26 17:59:35] INFO: Loading mms3_edp_dce_dsl_brst_l2...
[05-Oct-26 17:59:35] INFO: Loading mms3_edp_scpot_brst_l2...

Example operations

Omnidirectional PSD

vdf_e_omni = mms.vdf_omni(vdf_e)

Compute the PAD

vdf_e_pad = mms.get_pitch_angle_dist(vdf_e, b_xyz, tint=tint, angles=24)

[05-Oct-26 17:59:35] INFO: User defined number of pitch angles.

Limit the energy range

vdf_e_lowen = mms.vdf_elim(vdf_e, [0, 200])

[05-Oct-26 17:59:42] INFO: Effective eint = [10.96, 191.15]

Convert to differential energy flux (DEF)

vdf_e_deflux = mms.psd2def(vdf_e)

Convert to differential particle flux

vdf_e_dpflux = mms.psd2dpf(vdf_e)

Rebin to 64 energy channels (halves the time resolution)

vdf_e_e64 = mms.vdf_to_e64(vdf_e)

e_lim_low = [2e1, 2e2]
vdf_e_pa_lowen = mms.get_pitch_angle_dist(
    mms.vdf_elim(vdf_e_e64, e_lim_low), b_xyz, tint=tint, angles=18
)
vdf_e_pa_lowen_spectr = xr.DataArray(
    1e12 * np.nanmean(vdf_e_pa_lowen.data, axis=1),
    coords=[vdf_e_pa_lowen.time.data, vdf_e_pa_lowen.theta.data[0, :]],
    dims=["time", "theta"],
)

e_lim_mid = [2e2, 2e3]
vdf_e_pa_miden = mms.get_pitch_angle_dist(
    mms.vdf_elim(vdf_e_e64, e_lim_mid), b_xyz, tint=tint, angles=18
)

vdf_e_pa_miden_spectr = xr.DataArray(
    1e12 * np.nanmean(vdf_e_pa_miden.data, axis=1),
    coords=[vdf_e_pa_miden.time.data, vdf_e_pa_miden.theta.data[0, :]],
    dims=["time", "theta"],
)

pa_lim_low = [0.0, 15.0]
vdf_e_lowan = mms.get_pitch_angle_dist(vdf_e_e64, b_xyz, tint=tint, angles=pa_lim_low)
vdf_e_lowan_spectr = xr.DataArray(
    1e12 * np.nanmean(vdf_e_lowan.data, axis=2),
    coords=[vdf_e_lowan.time.data, vdf_e_lowan.energy.data[0, :]],
    dims=["time", "energy"],
)

pa_lim_mid = [75.0, 105.0]
vdf_e_midan = mms.get_pitch_angle_dist(vdf_e_e64, b_xyz, tint=tint, angles=pa_lim_mid)
vdf_e_midan_spectr = xr.DataArray(
    1e12 * np.nanmean(vdf_e_midan.data, axis=2),
    coords=[vdf_e_midan.time.data, vdf_e_midan.energy.data[0, :]],
    dims=["time", "energy"],
)

pa_lim_hig = [165.0, 180.0]
vdf_e_higan = mms.get_pitch_angle_dist(vdf_e_e64, b_xyz, tint=tint, angles=pa_lim_hig)
vdf_e_higan_spectr = xr.DataArray(
    1e12 * np.nanmean(vdf_e_higan.data, axis=2),
    coords=[vdf_e_higan.time.data, vdf_e_higan.energy.data[0, :]],
    dims=["time", "energy"],
)

[05-Oct-26 17:59:43] INFO: Effective eint = [20.40, 191.15]
[05-Oct-26 17:59:43] INFO: User defined number of pitch angles.
[05-Oct-26 17:59:44] INFO: Effective eint = [216.45, 1790.88]
[05-Oct-26 17:59:44] INFO: User defined number of pitch angles.
[05-Oct-26 17:59:46] INFO: User defined pitch angle limits.
[05-Oct-26 17:59:47] INFO: User defined pitch angle limits.
[05-Oct-26 17:59:48] INFO: User defined pitch angle limits.

Plot

f, axs = plt.subplots(7, sharex="all", figsize=(8, 10))
f.subplots_adjust(hspace=0, left=0.1, right=0.82, bottom=0.05, top=0.95)

plot_line(axs[0], b_xyz)
plot_line(axs[0], pyrf.norm(b_xyz), color="k")
axs[0].legend(["$B_x$", "$B_y$", "$B_z$", "$|B|$"], ncols=4, frameon=False)
axs[0].set_ylim([-30, 90])
axs[0].set_ylabel("$B$ [nT]")
axs[0].set_title(f"MMS{mms_id:d}")


axs[1], caxs1 = plot_spectr(
    axs[1],
    1e12 * mms.vdf_omni(vdf_e),
    yscale="log",
    cscale="log",
    clim=[1e-18, 1e-13],
    cmap="Spectral_r",
)
axs[1].set_yticks(np.logspace(1, 4, 4))
caxs1.set_ylabel(r"$f_e~[\mathrm{s}^{3}~\mathrm{m}^{-6}]$")
axs[1].set_ylabel(r"$E_e~[\mathrm{eV}]$")


axs[2] = plot_spectr(
    axs[2],
    vdf_e_pa_lowen_spectr,
    cscale="log",
    clim=caxs1.get_ylim(),
    colorbar="none",
    cmap="Spectral_r",
)
axs[2].set_yticks([0, 45, 90, 135])
axs[2].set_ylabel(r"$\theta~[\mathrm{{deg.}}]$")

axs[2].text(
    0.05,
    0.8,
    f"{e_lim_low[0]:3.0f}"
    + r"$~\mathrm{{eV}} < E_e < $"
    + f"{e_lim_low[1]:3.0f}"
    + r"$~\mathrm{{eV}}$",
    transform=axs[2].transAxes,
    bbox=dict(boxstyle="square", ec=(1.0, 1.0, 1.0), fc=(1.0, 1.0, 1.0)),
)


axs[3] = plot_spectr(
    axs[3],
    vdf_e_pa_miden_spectr,
    cscale="log",
    clim=[1e-18, 1e-13],
    colorbar="none",
    cmap="Spectral_r",
)
axs[3].set_yticks([0, 45, 90, 135])
axs[3].set_ylabel(r"$\theta~[\mathrm{{deg.}}]$")

axs[3].text(
    0.05,
    0.8,
    f"{e_lim_mid[0]:3.0f}"
    + r"$~\mathrm{{eV}} < E_e < $"
    + f"{e_lim_mid[1]:3.0f}"
    + r"$~\mathrm{{eV}}$",
    transform=axs[3].transAxes,
    bbox=dict(boxstyle="square", ec=(1.0, 1.0, 1.0), fc=(1.0, 1.0, 1.0)),
)


axs[4] = plot_spectr(
    axs[4],
    vdf_e_lowan_spectr,
    yscale="log",
    cscale="log",
    clim=[1e-18, 1e-13],
    colorbar="none",
    cmap="Spectral_r",
)
axs[4].set_ylabel(r"$E_e~[\mathrm{eV}]$")

axs[4].text(
    0.05,
    0.8,
    f"{pa_lim_low[0]:3.0f}"
    + r"$~\mathrm{deg.} < \theta < $"
    + f"{pa_lim_low[1]:3.0f}"
    + r"$~\mathrm{deg.}$",
    transform=axs[4].transAxes,
)


axs[5] = plot_spectr(
    axs[5],
    vdf_e_midan_spectr,
    yscale="log",
    cscale="log",
    clim=[1e-18, 1e-13],
    colorbar="none",
    cmap="Spectral_r",
)
axs[5].set_ylabel(r"$E_e~[\mathrm{eV}]$")

axs[5].text(
    0.05,
    0.8,
    f"{pa_lim_mid[0]:3.0f}"
    + r"$~\mathrm{deg.} < \theta < $"
    + f"{pa_lim_mid[1]:3.0f}"
    + r"$~\mathrm{deg.}$",
    transform=axs[5].transAxes,
)


axs[6] = plot_spectr(
    axs[6],
    vdf_e_higan_spectr,
    yscale="log",
    cscale="log",
    clim=[1e-18, 1e-13],
    colorbar="none",
    cmap="Spectral_r",
)
axs[6].set_ylabel(r"$E_e~[\mathrm{eV}]$")

axs[6].text(
    0.05,
    0.8,
    f"{pa_lim_hig[0]:3.0f}"
    + r"$~\mathrm{deg.} < \theta < $"
    + f"{pa_lim_hig[1]:3.0f}"
    + r"$~\mathrm{deg.}$",
    transform=axs[6].transAxes,
)

pos_axs6 = axs[6].get_position()
pos_cax1 = caxs1.get_position()
x0 = pos_cax1.x0
y0 = pos_axs6.y0
width = pos_cax1.width
height = pos_cax1.y0 + pos_cax1.height - y0
caxs1.set_position([x0, y0, width, height])

Project the distribution onto the (E, E\(\times\)B), (E\(\times\)B, B) and (B, E) planes

Compute the PAD with 17 pitch-angle bins

vdf_e_pad = mms.get_pitch_angle_dist(vdf_e, b_xyz, tint, angles=17)

[05-Oct-26 17:59:51] INFO: User defined number of pitch angles.

Resample the magnetic and electric fields and compute E\(\times\)B

b_0 = pyrf.resample(b_xyz, vdf_e)
e_0 = pyrf.resample(e_xyz, vdf_e)
exb = pyrf.cross(e_0, b_0)

Plot the projected distributions

idx = 1339
x = e_0.data[idx, :]
y = exb.data[idx, :]
z = b_0.data[idx, :]
time = list(pyrf.datetime642iso8601(vdf_e.time.data[idx]))

f = plt.figure(figsize=(9, 7.5))
gsp1 = f.add_gridspec(2, 3, hspace=0, bottom=0.07, top=0.99, left=0.1, right=0.9)

gsp10 = gsp1[0, :].subgridspec(1, 3, hspace=0)
gsp11 = gsp1[1, :].subgridspec(1, 2, hspace=0)

# Create the axes in the grid spec
axs10 = [f.add_subplot(gsp10[i]) for i in range(3)]
axs11 = [f.add_subplot(gsp11[i]) for i in range(2)]

f.subplots_adjust(wspace=0.4)
v_x, v_y, f_mat = mms.vdf_projection(
    vdf_e, time, np.vstack([x, y, -z]), sc_pot, e_lim=15
)
axs10[0], caxs10 = plot_projection(
    axs10[0], v_x, v_y, f_mat * 1e12, vlim=12e3, clim=[-18, -13], colorbar="top"
)
axs10[0].set_xlabel(r"$V_{E}~[\mathrm{Mm}~\mathrm{s}^{-1}]$")
axs10[0].set_ylabel(r"$V_{E\times B}~[\mathrm{Mm}~\mathrm{s}^{-1}]$")
caxs10.set_xlabel(r"$\mathrm{log}_{10}f_e~[\mathrm{s}^{3}~\mathrm{m}^{-6}]$")

v_x, v_y, f_mat = mms.vdf_projection(
    vdf_e, time, np.vstack([y, z, -x]), sc_pot, e_lim=15
)
axs10[1], caxs11 = plot_projection(
    axs10[1], v_x, v_y, f_mat * 1e12, vlim=12e3, clim=[-18, -13], colorbar="top"
)
axs10[1].set_xlabel(r"$V_{E\times B}~[\mathrm{Mm}~\mathrm{s}^{-1}]$")
axs10[1].set_ylabel(r"$V_{B}~[\mathrm{Mm}~\mathrm{s}^{-1}]$")
caxs11.set_xlabel(r"$\mathrm{log}_{10}f_e~[\mathrm{s}^{3}~\mathrm{m}^{-6}]$")

v_x, v_y, f_mat = mms.vdf_projection(
    vdf_e, time, np.vstack([z, x, -y]), sc_pot, e_lim=15
)
axs10[2], caxs12 = plot_projection(
    axs10[2], v_x, v_y, f_mat * 1e12, vlim=12e3, clim=[-18, -13], colorbar="top"
)
axs10[2].set_xlabel(r"$V_{B}~[\mathrm{Mm}~\mathrm{s}^{-1}]$")
axs10[2].set_ylabel(r"$V_{E}~[\mathrm{Mm}~\mathrm{s}^{-1}]$")
caxs12.set_xlabel(r"$\mathrm{log}_{10}f_e~[\mathrm{s}^{3}~\mathrm{m}^{-6}]$")


axs11[0].loglog(
    vdf_e_pad.energy.data[idx, :],
    1e12 * vdf_e_pad.data.data[idx, :, 0],
    label=r"$\theta = 0~\mathrm{deg.}$",
)
axs11[0].loglog(
    vdf_e_pad.energy.data[idx, :],
    1e12 * vdf_e_pad.data.data[idx, :, 9],
    label=r"$\theta = 90~\mathrm{deg.}$",
)
axs11[0].loglog(
    vdf_e_pad.energy.data[idx, :],
    1e12 * vdf_e_pad.data.data[idx, :, -1],
    label=r"$\theta = 180~\mathrm{deg.}$",
)

axs11[0].legend(loc="lower left")
axs11[0].set_xlim([1e1, 1e3])
axs11[0].set_xlabel(r"$E_e~[\mathrm{eV}]$")
axs11[0].set_ylim([1e-18, 1e-13])
axs11[0].set_ylabel(r"$f_e~[\mathrm{s}^{3}~\mathrm{m}^{-6}]$")

colormap = mpl.colormaps.get_cmap("hot_desaturated")
colors = colormap(np.linspace(0, 1, len(vdf_e_pad.energy[idx, :])))
for i_en in range(len(vdf_e_pad.energy[idx, :])):
    axs11[1].semilogy(
        vdf_e_pad.theta.data[idx, :],
        1e12 * vdf_e_pad.data.data[idx, i_en, :],
        color=colors[i_en],
        label=f"{vdf_e_pad.energy.data[idx, i_en]:5.2f} eV",
    )

axs11[1].set_xlim([0, 180.0])
axs11[1].set_xlabel(r"$\theta~[\mathrm{deg.}]$")
axs11[1].set_ylim([1e-18, 1e-13])
axs11[1].set_ylabel(r"$f_e~[\mathrm{s}^{3}~\mathrm{m}^{-6}]$")

axs11[1].set_xticks([0, 45, 90, 135, 180])
make_labels(axs10, (0.03, 0.90), pad=0, color="w")
make_labels(axs11, (0.03, 0.94), pad=3, color="k")
f.suptitle(time[0])

[05-Oct-26 17:59:57] WARNING: In making xyz a right handed orthogonal coordinate system, z (in-plane 2) was changed from [-10.36005306  -1.19106516 -34.19923401] to [-0.34608795  0.10882445 -0.93186929]. Please verify that this is according to your intentions.
[05-Oct-26 17:59:57] WARNING: In making xyz a right handed orthogonal coordinate system, z (in-plane 2) was changed from [-2.18016844  5.88453277  1.49689565] to [-0.3873171   0.91798646  0.08535992]. Please verify that this is according to your intentions.
[05-Oct-26 17:59:57] WARNING: In making xyz a right handed orthogonal coordinate system, y (in-plane 1) was changed from [ 2.18016844 -5.88453277 -1.49689565] to [ 0.3873171  -0.91798646 -0.08535992]. Please verify that this is according to your intentions.

Text(0.5, 0.98, '2015-12-02T01:14:55.191087000')