<jupyter_start><jupyter_text>Particle Distributionsauthor: <NAME>\Example showing how to you can work with particle distributions<jupyter_code>import numpy as np
import xarray as xr
import matplotlib.pylab as pl
import matplotlib.pyplot as plt

from pyrfu import mms, pyrf
from pyrfu.plot import plot_line, plot_spectr, plot_projection<jupyter_output><empty_output><jupyter_text>Define spacecraft index, time interval and data path<jupyter_code>mms_id = 3
tint = ["2015-12-02T01:14:15.000", "2015-12-02T01:15:13.000"]
data_path = "/Users/louisr/Documents/PhD/irfu-python-draft/pydata"<jupyter_output><empty_output><jupyter_text>Load Velocity Distribution Functions (VDFs)<jupyter_code>vdf_i = mms.get_data("PDi_fpi_brst_l2", tint, mms_id, data_path=data_path)
vdf_e = mms.get_data("PDe_fpi_brst_l2", tint, mms_id, data_path=data_path)<jupyter_output>Loading mms3_dis_dist_brst...
Loading mms3_des_dist_brst...<jupyter_text>Load supporting information<jupyter_code>b_xyz = mms.get_data("B_dmpa_fgm_brst_l2", tint, mms_id, data_path=data_path)
e_xyz = mms.get_data("e_dsl_edp_brst_l2", tint, mms_id, data_path=data_path)
sc_pot = mms.get_data("v_edp_brst_l2", tint, mms_id, data_path=data_path)<jupyter_output>Loading mms3_fgm_b_dmpa_brst_l2...
Loading mms3_edp_dce_dsl_brst_l2...
Loading mms3_edp_scpot_brst_l2...<jupyter_text>Example operations  Omnidirectional differential energy flux<jupyter_code>vdf_e_omni = mms.vdf_omni(vdf_e)<jupyter_output><empty_output><jupyter_text>Construt pitchangle distribution<jupyter_code>vdf_e_pad = mms.get_pitch_angle_dist(vdf_e, b_xyz, tint=tint, angles=24)<jupyter_output>notice : User defined number of pitch angles.<jupyter_text>Limit energy range<jupyter_code>vdf_e_lowen = mms.vdf_elim(vdf_e, [0, 200])<jupyter_output>Effective eint = [10.96, 191.15]<jupyter_text>Change units to differential energy flux<jupyter_code>vdf_e_deflux = mms.vdf_to_deflux(vdf_e)<jupyter_output><empty_output><jupyter_text>Change units to particle energy flux<jupyter_code>vdf_e_dpflux = mms.vdf_to_dpflux(vdf_e)<jupyter_output><empty_output><jupyter_text>Resample energy to 64 energy levels, reduces the time resolution<jupyter_code>vdf_e_e64 = mms.vdf_to_e64(vdf_e)
vdf_e_pa_lowen = mms.get_pitch_angle_dist(mms.vdf_elim(vdf_e_e64, [20, 200]), b_xyz, tint=tint, angles=18)
vdf_e_pa_lowen_spectr = xr.DataArray(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"])

vdf_e_pa_miden = mms.get_pitch_angle_dist(mms.vdf_elim(vdf_e_e64, [200, 2000]), b_xyz, tint=tint, angles=18)
vdf_e_pa_miden_spectr = xr.DataArray(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"])<jupyter_output>Effective eint = [20.40, 191.15]
notice : User defined number of pitch angles.<jupyter_text>Plot<jupyter_code>%matplotlib notebook

n_panels = 7;
f, axs = plt.subplots(n_panels, sharex="all", figsize=(8.5, 11))
f.subplots_adjust(hspace=0, left=.13, right=.87, bottom=.05, top=.95)

plot_line(axs[0], b_xyz)
plot_line(axs[0], pyrf.norm(b_xyz))
axs[0].legend(["$B_x$", "$B_y$", "$B_z$", "$|B|$"], bbox_to_anchor=(1, 1.05))
axs[0].set_ylabel("$B$ [nT]")
axs[0].set_title(f"MMS{mms_id:d}")


axs[1], caxs1 = plot_spectr(axs[1], mms.vdf_omni(vdf_e), yscale="log", 
                            cscale="log", cmap="jet")
axs[1].set_yticks(np.logspace(1, 4, 4))
caxs1.set_ylabel("PSD" + "\n" + "[s$^{3}$ cm$^{-6}$]")
axs[1].set_ylabel("$E_e$ [eV]")

e_lim = [20, 200]
vdf_e_pa_lowen = mms.get_pitch_angle_dist(mms.vdf_elim(vdf_e_e64, e_lim), b_xyz, tint=tint, angles=18)
vdf_e_pa_lowen_spectr = xr.DataArray(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"])

axs[2], caxs2 = plot_spectr(axs[2], vdf_e_pa_lowen_spectr, 
                            cscale="log", cmap="jet")
axs[2].set_yticks([0, 45, 90, 135])
caxs2.set_ylabel("PSD" + "\n" + "[s$^{3}$ cm$^{-6}$]")
axs[2].set_ylabel("$\\theta$ [deg.]")

axs[2].text(.03, .1, f"{e_lim[0]:d} < $E_e$ < {e_lim[1]:d} eV", transform=axs[2].transAxes,
            bbox=dict(boxstyle="square", ec=(1., 1., 1.), fc=(1., 1., 1.)))


e_lim = [200, 2000]
vdf_e_pa_miden = mms.get_pitch_angle_dist(mms.vdf_elim(vdf_e_e64, e_lim), b_xyz, tint=tint, angles=18)
vdf_e_pa_miden_spectr = xr.DataArray(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"])

axs[3], caxs3 = plot_spectr(axs[3], vdf_e_pa_miden_spectr, 
                            cscale="log", cmap="jet")
axs[3].set_yticks([0, 45, 90, 135])
caxs3.set_ylabel("PSD" + "\n" + "[s$^{3}$ cm$^{-6}$]")
axs[3].set_ylabel("$\\theta$ [deg.]")

axs[3].text(.03, .1, f"{e_lim[0]:d} < $E_e$ < {e_lim[1]:d} eV", transform=axs[3].transAxes,
            bbox=dict(boxstyle="square", ec=(1., 1., 1.), fc=(1., 1., 1.)))


pa_lim = [0, 15]
vdf_e_lowan = mms.get_pitch_angle_dist(mms.vdf_to_e64(vdf_e), b_xyz, tint=tint, angles=pa_lim)
vdf_e_lowan_spectr = xr.DataArray(np.nanmean(vdf_e_lowan.data, axis=2), 
                                  coords=[vdf_e_lowan.time.data, vdf_e_lowan.energy.data[0, :]], 
                                  dims=["time", "energy"])

axs[4], caxs4 = plot_spectr(axs[4], vdf_e_lowan_spectr, yscale="log", 
                            cscale="log", cmap="jet")
caxs4.set_ylabel("PSD" + "\n" + "[s$^{3}$ cm$^{-6}$]")
axs[4].set_ylabel("$E_e$ [eV]")

axs[4].text(.03, .1, f"{pa_lim[0]:d}$^\\circ$ < $\\theta$ < {pa_lim[1]:d}$^\\circ$", 
            transform=axs[4].transAxes,
            bbox=dict(boxstyle="square", ec=(1., 1., 1.), fc=(1., 1., 1.)))

pa_lim = [75, 105]
vdf_e_midan = mms.get_pitch_angle_dist(mms.vdf_to_e64(vdf_e), b_xyz, tint=tint, angles=pa_lim)
vdf_e_midan_spectr = xr.DataArray(np.nanmean(vdf_e_midan.data, axis=2), 
                                  coords=[vdf_e_midan.time.data, vdf_e_midan.energy.data[0, :]], 
                                  dims=["time", "energy"])

axs[5], caxs5 = plot_spectr(axs[5], vdf_e_midan_spectr, yscale="log", 
                            cscale="log", cmap="jet")
caxs5.set_ylabel("PSD" + "\n" + "[s$^{3}$ cm$^{-6}$]")
axs[5].set_ylabel("$E_e$ [eV]")

axs[5].text(.03, .1, f"{pa_lim[0]:d}$^\\circ$ < $\\theta$ < {pa_lim[1]:d}$^\\circ$", 
            transform=axs[5].transAxes,
            bbox=dict(boxstyle="square", ec=(1., 1., 1.), fc=(1., 1., 1.)))

pa_lim = [165, 180]
vdf_e_higan = mms.get_pitch_angle_dist(mms.vdf_to_e64(vdf_e), b_xyz, tint=tint, angles=pa_lim)
vdf_e_higan_spectr = xr.DataArray(np.nanmean(vdf_e_higan.data, axis=2), 
                                  coords=[vdf_e_higan.time.data, vdf_e_higan.energy.data[0, :]], 
                                  dims=["time", "energy"])

axs[6], caxs6 = plot_spectr(axs[6], vdf_e_higan_spectr, yscale="log", 
                            cscale="log", cmap="jet")
caxs6.set_ylabel("PSD" + "\n" + "[s$^{3}$ cm$^{-6}$]")
axs[6].set_ylabel("$E_e$ [eV]")

axs[6].text(.03, .1, f"{pa_lim[0]:d}$^\\circ$ < $\\theta$ < {pa_lim[1]:d}$^\\circ$", 
            transform=axs[6].transAxes,
            bbox=dict(boxstyle="square", ec=(1., 1., 1.), fc=(1., 1., 1.)))<jupyter_output><empty_output><jupyter_text>Project distribution onto the (E, ExB), (ExB, B), (B, E)  Compute pitchangle distribution with 17 angles<jupyter_code>vdf_e_pad = mms.get_pitch_angle_dist(vdf_e, b_xyz, tint, angles=17)<jupyter_output>notice : User defined number of pitch angles.<jupyter_text>Resample background magnetic field, electric field and ExB drift<jupyter_code>b_0 = pyrf.resample(b_xyz, vdf_e)
e_0 = pyrf.resample(e_xyz, vdf_e)
exb = pyrf.cross(e_0, b_0)<jupyter_output>/Users/louisr/opt/anaconda3/lib/python3.8/site-packages/pyrfu/pyrf/resample.py:223: UserWarning: Using averages in resample
  warnings.warn("Using averages in resample", UserWarning)<jupyter_text>Plot<jupyter_code>%matplotlib notebook

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=(8.5, 9))
gsp1 = f.add_gridspec(2, 3, hspace=0, bottom=.07, top=.99, left=.1, right=.9)

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

# Create 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=.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], cbar_pos="top")
axs10[0].set_xlabel("$V_{E}$ [Mm s$^{-1}$]")
axs10[0].set_ylabel("$V_{E\\times B}$ [Mm s$^{-1}$]")
caxs10.set_xlabel("log$_{10} f_e$ [s$^{3}$ 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], cbar_pos="top")
axs10[1].set_xlabel("$V_{E\\times B}$ [Mm s$^{-1}$]")
axs10[1].set_ylabel("$V_{B}$ [Mm s$^{-1}$]")
caxs11.set_xlabel("log$_{10} f_e$ [s$^{3}$ 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], cbar_pos="top")
axs10[2].set_xlabel("$V_{B}$ [Mm s$^{-1}$]")
axs10[2].set_ylabel("$V_{E}$ [Mm s$^{-1}$]")
caxs12.set_xlabel("log$_{10} f_e$ [s$^{3}$ m$^{-6}$]")


axs11[0].loglog(vdf_e_pad.energy.data[idx, :], vdf_e_pad.data.data[idx, :, 0], 
                label="$\\theta = 0$ deg.")
axs11[0].loglog(vdf_e_pad.energy.data[idx, :], vdf_e_pad.data.data[idx, :, 9], 
                label="$\\theta = 90$ deg.")
axs11[0].loglog(vdf_e_pad.energy.data[idx, :], vdf_e_pad.data.data[idx, :, -1], 
                label="$\\theta = 180$ deg.")

axs11[0].legend()
axs11[0].set_xlim([1e1, 1e3])
axs11[0].set_xlabel("$E_e$ [eV]")
axs11[0].set_ylim([1e-31, 2e-26])
axs11[0].set_ylabel("$f_e$ [s$^{3}$ cm$^{-6}$]")


colors = pl.cm.jet(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, :], 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.])
axs11[1].set_xlabel("$\\theta$ [deg.]")
axs11[1].set_ylim([1e-31, 2e-26])
axs11[1].set_ylabel("$f_e$ [s$^{3}$ cm$^{-6}$]")

axs11[1].set_xticks([0, 45, 90, 135, 180])

f.suptitle(time[0])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1944772458897266, 'fraction_numerical': 0.051049141611408846, 'mean_word_length': 4.036596523330283, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 51, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '486649', 'license': 'MIT', 'n_tokens_mistral': 5397, 'n_tokens_neox': 5009, 'n_words': 687, 'path': 'examples/01_mms/.ipynb_checkpoints/example_mms_particle_distributions-checkpoint.ipynb', 'repo_name': 'ablotekar/irfu-python', 'size': '11016'}