Note to Deepnote users: to modify this notebook, duplicate it into your own workspace. Click the three dots (...) "Project settings" menu in the top right corner of the page then select "Duplicate project"

Making an Executable Paper with the Python in Heliophysics Community to Foster Open Science and Improve Reproducibility

Abstract

1 Introduction

2 Method

# pySPEDAS imports import pyspedas # PyTplot imports import pytplot from pytplot import tplot # SpacePy imports import matplotlib.pyplot as plt import numpy as np import spacepy.coordinates import spacepy.empiricals import spacepy.time import spacepy.toolbox # PlasmaPy imports import plasmapy from astropy import units as u # Kamodo imports from kamodo_ccmc.flythrough import SatelliteFlythrough as S import kamodo_ccmc.flythrough.model_wrapper as MW

kj/filesystem-disk-unix.c++:1690: warning: PWD environment variable doesn't match current directory; pwd = /

trange = ['2015-10-16/11:30', '2015-10-16/17:00']

def load_mms_data(time_range): """Loads MMS MEC, FGM, and FPI data.""" mec_vars = pyspedas.mms.mec(trange=time_range, time_clip=True, probe=1, level='l2', data_rate='srvy', varformat='*_r_gsm') # Cartesian GSM coords fgm_vars = pyspedas.mms.fgm(trange=time_range, time_clip=True, probe=1, level='l2', data_rate='srvy') ion_vars = pyspedas.mms.fpi(trange=time_range, datatype=['dis-moms', 'des-moms'], level='l2', data_rate='fast', time_clip=True, center_measurement=True) return mec_vars, fgm_vars, ion_vars

2.3 Compare MMS Data to an Empirical Model

def spacecraft_magnetopause_calculations(mms_mec_vars): """Returns epoch, pos, distance between spacecraft and magnetopause, magnetopause's distance from Earth, spacecraft's distance from Earth, and solar zenith angle. """ data = pytplot.get_data(mms_mec_vars[0]) pos_gsm = data.y ticks = spacepy.time.Ticktock(data.times, dtype='UNX') epoch = ticks.UTC c = spacepy.coordinates.Coords(pos_gsm, 'GSM', 'car', units='km', ticks=ticks) pos = c.convert('GSE', 'car').data # Get the Shue coefficients alpha = [] # Shue flaring angle standoff = spacepy.empiricals.getMPstandoff(ticks, alpha=alpha) # Shue subsolar standoff alpha = np.array(alpha) # Solar zenith angle of s/c position (angle with GSE +x) sza = np.arctan2(pos[:, 0], np.sqrt((pos[:, :2] ** 2).sum(axis=1))) # Radial distance to MP along Earth-SC line (application of Shue coefficients) mp_dist = standoff * (2. / (1 + np.cos(sza))) ** alpha # Radial distance to SC sc_dist = np.sqrt((pos ** 2).sum(axis=1)) / 6378 # How far is SC outside of MP? sc_to_mp = sc_dist - mp_dist return epoch, pos, sc_to_mp, mp_dist, sc_dist, sza

def find_magnetopause_crossings(sc_to_mp): """Returns the indices of elements before which zero crossings occur.""" return np.where(np.diff(np.sign(sc_to_mp)))[0]

mec_vars, fgm_vars, ion_vars = load_mms_data(trange) # loaded Tplot variables are printed as output epoch, pos, sc_to_mp, mp_dist, sc_dist, sza = spacecraft_magnetopause_calculations(mec_vars)

22-Nov-22 17:14:40: Loading pydata/mms1/mec/srvy/l2/epht89q/2015/10/mms1_mec_srvy_l2_epht89q_20151016_v2.0.0.cdf
Time clip was applied to: mms1_mec_r_gsm
Loaded variables:
mms1_mec_r_gsm
22-Nov-22 17:14:41: Loading pydata/mms1/fgm/srvy/l2/2015/10/mms1_fgm_srvy_l2_20151016_v4.18.0.cdf
The lengths of x and y do not match!
mms1_fgm_rdeltahalf_srvy_l2 is currently not in pytplot.
mms1_fgm_rdeltahalf_srvy_l2 is currently not in pytplot.
mms1_fgm_rdeltahalf_srvy_l2 is currently not in pytplot.
mms1_fgm_rdeltahalf_srvy_l2 is currently not in pytplot.
Time clip was applied to: mms1_fgm_b_gse_srvy_l2
Time clip was applied to: mms1_fgm_b_gsm_srvy_l2
Time clip was applied to: mms1_fgm_b_dmpa_srvy_l2
Time clip was applied to: mms1_fgm_b_bcs_srvy_l2
Time clip was applied to: mms1_fgm_flag_srvy_l2
Time clip was applied to: mms1_fgm_r_gse_srvy_l2
Time clip was applied to: mms1_fgm_r_gsm_srvy_l2
Time clip was applied to: mms1_fgm_hirange_srvy_l2
Time clip was applied to: mms1_fgm_bdeltahalf_srvy_l2
Time clip was applied to: mms1_fgm_stemp_srvy_l2
Time clip was applied to: mms1_fgm_etemp_srvy_l2
Time clip was applied to: mms1_fgm_mode_srvy_l2
Time clip error: No pytplot names were provided.
That name is currently not in pytplot
That name is currently not in pytplot
Loaded variables:
Epoch
mms1_fgm_b_gse_srvy_l2
mms1_fgm_b_gsm_srvy_l2
mms1_fgm_b_dmpa_srvy_l2
mms1_fgm_b_bcs_srvy_l2
mms1_fgm_flag_srvy_l2
Epoch_state
mms1_fgm_hirange_srvy_l2
mms1_fgm_bdeltahalf_srvy_l2
mms1_fgm_stemp_srvy_l2
mms1_fgm_etemp_srvy_l2
mms1_fgm_mode_srvy_l2
mms1_fgm_rdeltahalf_srvy_l2
mms1_fgm_b_dmpa_srvy_l2_bvec
mms1_fgm_b_dmpa_srvy_l2_btot
mms1_fgm_b_gse_srvy_l2_bvec
mms1_fgm_b_gse_srvy_l2_btot
mms1_fgm_b_gsm_srvy_l2_bvec
mms1_fgm_b_gsm_srvy_l2_btot
mms1_fgm_b_bcs_srvy_l2_bvec
mms1_fgm_b_bcs_srvy_l2_btot
22-Nov-22 17:14:53: Loading pydata/mms1/fpi/fast/l2/dis-moms/2015/10/mms1_fpi_fast_l2_dis-moms_20151016100000_v3.3.0.cdf
22-Nov-22 17:14:53: Loading pydata/mms1/fpi/fast/l2/dis-moms/2015/10/mms1_fpi_fast_l2_dis-moms_20151016120000_v3.3.0.cdf
22-Nov-22 17:14:53: Loading pydata/mms1/fpi/fast/l2/dis-moms/2015/10/mms1_fpi_fast_l2_dis-moms_20151016140000_v3.3.0.cdf
22-Nov-22 17:14:53: Loading pydata/mms1/fpi/fast/l2/dis-moms/2015/10/mms1_fpi_fast_l2_dis-moms_20151016160000_v3.3.0.cdf
22-Nov-22 17:14:53: Loading pydata/mms1/fpi/fast/l2/des-moms/2015/10/mms1_fpi_fast_l2_des-moms_20151016100000_v3.3.0.cdf
22-Nov-22 17:14:53: Loading pydata/mms1/fpi/fast/l2/des-moms/2015/10/mms1_fpi_fast_l2_des-moms_20151016120000_v3.3.0.cdf
22-Nov-22 17:14:53: Loading pydata/mms1/fpi/fast/l2/des-moms/2015/10/mms1_fpi_fast_l2_des-moms_20151016140000_v3.3.0.cdf
22-Nov-22 17:14:53: Loading pydata/mms1/fpi/fast/l2/des-moms/2015/10/mms1_fpi_fast_l2_des-moms_20151016160000_v3.3.0.cdf
22-Nov-22 17:14:53: /root/venv/lib/python3.7/site-packages/pytplot/importers/cdf_to_tplot.py:201: DeprecationWarning: `np.float` is a deprecated alias for the builtin `float`. To silence this warning, use `float` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.float64` here.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  delta_plus_var = delta_plus_var.astype(float)*np.float(si_conv.split('>')[0])

22-Nov-22 17:14:53: /root/venv/lib/python3.7/site-packages/pytplot/importers/cdf_to_tplot.py:213: DeprecationWarning: `np.float` is a deprecated alias for the builtin `float`. To silence this warning, use `float` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.float64` here.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  delta_minus_var = delta_minus_var.astype(float)*np.float(si_conv.split('>')[0])

Time clip was applied to: mms1_des_errorflags_fast
Time clip was applied to: mms1_des_startdelphi_count_fast
Time clip was applied to: mms1_des_startdelphi_angle_fast
Time clip was applied to: mms1_des_pitchangdist_lowen_fast
Time clip was applied to: mms1_des_pitchangdist_miden_fast
Time clip was applied to: mms1_des_pitchangdist_highen_fast
Time clip was applied to: mms1_des_energyspectr_px_fast
Time clip was applied to: mms1_des_energyspectr_mx_fast
Time clip was applied to: mms1_des_energyspectr_py_fast
Time clip was applied to: mms1_des_energyspectr_my_fast
Time clip was applied to: mms1_des_energyspectr_pz_fast
Time clip was applied to: mms1_des_energyspectr_mz_fast
Time clip was applied to: mms1_des_energyspectr_par_fast
Time clip was applied to: mms1_des_energyspectr_anti_fast
Time clip was applied to: mms1_des_energyspectr_perp_fast
Time clip was applied to: mms1_des_energyspectr_omni_fast
Time clip was applied to: mms1_des_numberdensity_fast
Time clip was applied to: mms1_des_numberdensity_err_fast
Time clip was applied to: mms1_des_densityextrapolation_low_fast
Time clip was applied to: mms1_des_densityextrapolation_high_fast
Time clip was applied to: mms1_des_bulkv_dbcs_fast
Time clip was applied to: mms1_des_bulkv_spintone_dbcs_fast
Time clip was applied to: mms1_des_bulkv_gse_fast
Time clip was applied to: mms1_des_bulkv_spintone_gse_fast
Time clip was applied to: mms1_des_bulkv_err_fast
Time clip was applied to: mms1_des_prestensor_dbcs_fast
Time clip was applied to: mms1_des_prestensor_gse_fast
Time clip was applied to: mms1_des_prestensor_err_fast
Time clip was applied to: mms1_des_temptensor_dbcs_fast
Time clip was applied to: mms1_des_temptensor_gse_fast
Time clip was applied to: mms1_des_temptensor_err_fast
Time clip was applied to: mms1_des_heatq_dbcs_fast
Time clip was applied to: mms1_des_heatq_gse_fast
Time clip was applied to: mms1_des_heatq_err_fast
Time clip was applied to: mms1_des_energy_fast
Time clip was applied to: mms1_des_energy_delta_fast
Time clip was applied to: mms1_des_temppara_fast
Time clip was applied to: mms1_des_tempperp_fast
Time clip was applied to: mms1_dis_errorflags_fast
Time clip was applied to: mms1_dis_startdelphi_count_fast
Time clip was applied to: mms1_dis_startdelphi_angle_fast
Time clip was applied to: mms1_dis_energyspectr_px_fast
Time clip was applied to: mms1_dis_energyspectr_mx_fast
Time clip was applied to: mms1_dis_energyspectr_py_fast
Time clip was applied to: mms1_dis_energyspectr_my_fast
Time clip was applied to: mms1_dis_energyspectr_pz_fast
Time clip was applied to: mms1_dis_energyspectr_mz_fast
Time clip was applied to: mms1_dis_energyspectr_omni_fast
Time clip was applied to: mms1_dis_spectr_bg_fast
Time clip was applied to: mms1_dis_numberdensity_bg_fast
Time clip was applied to: mms1_dis_numberdensity_fast
Time clip was applied to: mms1_dis_numberdensity_err_fast
Time clip was applied to: mms1_dis_densityextrapolation_low_fast
Time clip was applied to: mms1_dis_densityextrapolation_high_fast
Time clip was applied to: mms1_dis_bulkv_dbcs_fast
Time clip was applied to: mms1_dis_bulkv_spintone_dbcs_fast
Time clip was applied to: mms1_dis_bulkv_gse_fast
Time clip was applied to: mms1_dis_bulkv_spintone_gse_fast
Time clip was applied to: mms1_dis_bulkv_err_fast
Time clip was applied to: mms1_dis_prestensor_dbcs_fast
Time clip was applied to: mms1_dis_prestensor_gse_fast
Time clip was applied to: mms1_dis_prestensor_err_fast
Time clip was applied to: mms1_dis_pres_bg_fast
Time clip was applied to: mms1_dis_temptensor_dbcs_fast
Time clip was applied to: mms1_dis_temptensor_gse_fast
Time clip was applied to: mms1_dis_temptensor_err_fast
Time clip was applied to: mms1_dis_heatq_dbcs_fast
Time clip was applied to: mms1_dis_heatq_gse_fast
Time clip was applied to: mms1_dis_heatq_err_fast
Time clip was applied to: mms1_dis_energy_fast
Time clip was applied to: mms1_dis_energy_delta_fast
Time clip was applied to: mms1_dis_temppara_fast
Time clip was applied to: mms1_dis_tempperp_fast
That name is currently not in pytplot
That name is currently not in pytplot
That name is currently not in pytplot
22-Nov-22 17:14:58: /root/venv/lib/python3.7/site-packages/spacepy/coordinates.py:242: DeprecationWarning: No coordinate backend specified; using IRBEM. This default will change in the future.
  DeprecationWarning)

Loaded variables:
Epoch
Epoch_plus_var
Epoch_minus_var
mms1_des_errorflags_fast
mms1_des_compressionloss_fast
mms1_des_startdelphi_count_fast
mms1_des_startdelphi_angle_fast
mms1_des_pitchangdist_lowen_fast
mms1_des_pitchangdist_miden_fast
mms1_des_pitchangdist_highen_fast
mms1_des_energyspectr_px_fast
mms1_des_energyspectr_mx_fast
mms1_des_energyspectr_py_fast
mms1_des_energyspectr_my_fast
mms1_des_energyspectr_pz_fast
mms1_des_energyspectr_mz_fast
mms1_des_energyspectr_par_fast
mms1_des_energyspectr_anti_fast
mms1_des_energyspectr_perp_fast
mms1_des_energyspectr_omni_fast
mms1_des_numberdensity_fast
mms1_des_numberdensity_err_fast
mms1_des_densityextrapolation_low_fast
mms1_des_densityextrapolation_high_fast
mms1_des_bulkv_dbcs_fast
mms1_des_bulkv_spintone_dbcs_fast
mms1_des_bulkv_gse_fast
mms1_des_bulkv_spintone_gse_fast
mms1_des_bulkv_err_fast
mms1_des_prestensor_dbcs_fast
mms1_des_prestensor_gse_fast
mms1_des_prestensor_err_fast
mms1_des_temptensor_dbcs_fast
mms1_des_temptensor_gse_fast
mms1_des_temptensor_err_fast
mms1_des_heatq_dbcs_fast
mms1_des_heatq_gse_fast
mms1_des_heatq_err_fast
mms1_des_alpha_fast
mms1_des_alpha_delta_fast
mms1_des_energy_fast
mms1_des_energy_delta_fast
mms1_des_temppara_fast
mms1_des_tempperp_fast
mms1_dis_errorflags_fast
mms1_dis_compressionloss_fast
mms1_dis_startdelphi_count_fast
mms1_dis_startdelphi_angle_fast
mms1_dis_energyspectr_px_fast
mms1_dis_energyspectr_mx_fast
mms1_dis_energyspectr_py_fast
mms1_dis_energyspectr_my_fast
mms1_dis_energyspectr_pz_fast
mms1_dis_energyspectr_mz_fast
mms1_dis_energyspectr_omni_fast
mms1_dis_spectr_bg_fast
mms1_dis_numberdensity_bg_fast
mms1_dis_numberdensity_fast
mms1_dis_numberdensity_err_fast
mms1_dis_densityextrapolation_low_fast
mms1_dis_densityextrapolation_high_fast
mms1_dis_bulkv_dbcs_fast
mms1_dis_bulkv_spintone_dbcs_fast
mms1_dis_bulkv_gse_fast
mms1_dis_bulkv_spintone_gse_fast
mms1_dis_bulkv_err_fast
mms1_dis_prestensor_dbcs_fast
mms1_dis_prestensor_gse_fast
mms1_dis_prestensor_err_fast
mms1_dis_pres_bg_fast
mms1_dis_temptensor_dbcs_fast
mms1_dis_temptensor_gse_fast
mms1_dis_temptensor_err_fast
mms1_dis_heatq_dbcs_fast
mms1_dis_heatq_gse_fast
mms1_dis_heatq_err_fast
mms1_dis_energy_fast
mms1_dis_energy_delta_fast
mms1_dis_temppara_fast
mms1_dis_tempperp_fast
mms1_des_errorflags_fast_moms
mms1_des_compressionloss_fast_moms
mms1_dis_errorflags_fast_moms
mms1_dis_compressionloss_fast_moms
mms1_dis_errorflags_fast_moms_flagbars_full
mms1_dis_errorflags_fast_moms_flagbars_main
mms1_dis_errorflags_fast_moms_flagbars_others
mms1_dis_errorflags_fast_moms_flagbars_mini
mms1_des_pitchangdist_avg
mms1_des_errorflags_fast_moms_flagbars_full
mms1_des_errorflags_fast_moms_flagbars_main
mms1_des_errorflags_fast_moms_flagbars_others
mms1_des_errorflags_fast_moms_flagbars_mini

%matplotlib inline fig1, axs1 = plt.subplots(1, 1, figsize=(8, 4), sharey=True) fig2, axs2 = plt.subplots(1, 1, figsize=(8, 4), sharey=True) fig3, axs3 = plt.subplots(1, 1, figsize=(8, 4), sharey=True) fig4, axs4 = plt.subplots(1, 1, figsize=(8, 4), sharey=True) # Spacecraft distance outside magnetopause fig1.suptitle('(A) Distance Outside Magnetopause') axs1.plot(epoch, sc_to_mp) axs1.set(xlabel='Time', ylabel='Re') crossings = find_magnetopause_crossings(sc_to_mp) for x in crossings: axs1.plot(epoch[x], sc_to_mp[x], 'o', label='Crossing') axs1.legend() fig1.autofmt_xdate() # Spacecraft and magnetopause locations fig2.suptitle('(B) Distance From Earth') l_mp, = axs2.plot(epoch, mp_dist, label='Magnetopause') l_sc, = axs2.plot(epoch, sc_dist, label='Spacecraft') axs2.set(xlabel='Time', ylabel='Re') fig2.autofmt_xdate() fig2.legend([l_mp, l_sc], [l_mp.get_label(), l_sc.get_label()], loc='right') # Solar zenith angle fig3.suptitle('(C) Solar Zenith Angle') axs3.plot(epoch, np.degrees(sza)) axs3.set(xlabel='Time') fig3.autofmt_xdate() # Spacecraft coordinates fig4.suptitle('(D) Spacecraft Coordinates') l_x, = axs4.plot(epoch, pos[:, 0] / 6378, label='x') l_y, = axs4.plot(epoch, pos[:, 1] / 6378, label='y') l_z, = axs4.plot(epoch, pos[:, 2] / 6378, label='z') fig4.autofmt_xdate() fig4.legend([l_x, l_y, l_z], [l_x.get_label(), l_y.get_label(), l_z.get_label()], loc='right')

2.4 Augment Comparison with Plasma Parameters

from pyspedas import tinterpol tinterpol('mms1_fgm_b_gse_srvy_l2_btot', 'mms1_dis_numberdensity_fast') tinterpol('mms1_des_numberdensity_fast', 'mms1_dis_numberdensity_fast') tinterpol('mms1_des_temppara_fast', 'mms1_dis_numberdensity_fast') tinterpol('mms1_des_tempperp_fast', 'mms1_dis_numberdensity_fast')

tinterpol (linear) was applied to: mms1_fgm_b_gse_srvy_l2_btot-itrp
tinterpol (linear) was applied to: mms1_des_numberdensity_fast-itrp
tinterpol (linear) was applied to: mms1_des_temppara_fast-itrp
tinterpol (linear) was applied to: mms1_des_tempperp_fast-itrp

from pytplot import get_data fgm_b = get_data('mms1_fgm_b_gse_srvy_l2_btot-itrp') dis_n = get_data('mms1_dis_numberdensity_fast') dis_Tpara = get_data('mms1_dis_temppara_fast') dis_Tperp = get_data('mms1_dis_tempperp_fast') des_n = get_data('mms1_des_numberdensity_fast-itrp') des_Tpara = get_data('mms1_des_temppara_fast-itrp') des_Tperp = get_data('mms1_des_tempperp_fast-itrp')

dis_T = (dis_Tpara.y + 2*dis_Tperp.y) / 3.0 des_T = (des_Tpara.y + 2*des_Tperp.y) / 3.0

B = fgm_b.y * u.nT n_i = dis_n.y * u.cm**-3 n_e = des_n.y * u.cm**-3 T_i = dis_T * u.eV T_e = des_T * u.eV

Va = plasmapy.formulary.parameters.Alfven_speed(B, n_i, 'p') # convert to km / s Va = Va.to(u.km / u.s)

# ions beta_i = plasmapy.formulary.dimensionless.beta(T_i, n_i, B) # electrons beta_e = plasmapy.formulary.dimensionless.beta(T_e, n_e, B) # combined beta = beta_i + beta_e

n_i = np.float64(n_i.value) * u.cm ** -3 # convert values from float32 to float64 d_i = plasmapy.formulary.parameters.inertial_length(n_i, 'p+')

lamda_d = plasmapy.formulary.parameters.Debye_length(T_e, n_e)

omega_ci = plasmapy.formulary.parameters.gyrofrequency(B, 'p', to_hz=True)

r_i = plasmapy.formulary.parameters.gyroradius(B, 'p', T=T_i) # convert to km r_i = r_i.to(u.km)

DB = plasmapy.formulary.parameters.Bohm_diffusion(T_e, B)

omega_lh = plasmapy.formulary.parameters.lower_hybrid_frequency(B, n_i, 'p', to_hz=True)

omega_uh = plasmapy.formulary.parameters.upper_hybrid_frequency(B, n_e, to_hz=True)

from pytplot import store_data store_data('alfven_speed', data={'x': fgm_b.times, 'y': Va}) store_data('plasma_beta', data={'x': fgm_b.times, 'y': beta}) store_data('ion_inertial_length', data={'x': fgm_b.times, 'y': d_i}) store_data('debye_length', data={'x': fgm_b.times, 'y': lamda_d}) store_data('omega_ci', data={'x': fgm_b.times, 'y': omega_ci}) store_data('ion_gyroradius', data={'x': fgm_b.times, 'y': r_i}) store_data('bohm_diffusion_coeff', data={'x': fgm_b.times, 'y': DB}) store_data('omega_lh', data={'x': fgm_b.times, 'y': omega_lh}) store_data('omega_uh', data={'x': fgm_b.times, 'y': omega_uh})

from pytplot import options options('alfven_speed', 'ytitle', 'Va \\ (' + str(Va.unit) + ')') options('alfven_speed', 'legend_names', 'Alfvén speed') options('plasma_beta', 'ytitle', 'Beta') options('plasma_beta', 'legend_names', 'Plasma Beta') options('ion_inertial_length', 'ytitle', 'd_i \\ (' + str(d_i.unit) + ')') options('ion_inertial_length', 'legend_names', 'Ion inertial length') options('debye_length', 'ytitle', 'Lamda_d \\ (' + str(lamda_d.unit) + ')') options('debye_length', 'legend_names', 'Debye length') options('omega_ci', 'ytitle', 'omega_ci \\ (' + str(omega_ci.unit) + ')') options('omega_ci', 'legend_names', 'Ion gyrofrequency') options('ion_gyroradius', 'ytitle', 'r_i \\ (' + str(r_i.unit) + ')') options('ion_gyroradius', 'legend_names', 'Ion gyroradius') options('bohm_diffusion_coeff', 'ytitle', 'DB \\ (' + str(DB.unit) + ')') options('bohm_diffusion_coeff', 'legend_names', 'Bohm diffusion coefficient') options('omega_uh', 'ytitle', 'omega_uh \\ (' + str(omega_uh.unit) + ')') options('omega_uh', 'legend_names', 'Upper hybrid frequency') options('omega_lh', 'ytitle', 'omega_lh \\ (' + str(omega_lh.unit) + ')') options('omega_lh', 'legend_names', 'Lower hybrid frequency')

tplot(['alfven_speed', 'plasma_beta', 'ion_inertial_length', 'debye_length', 'omega_ci', 'ion_gyroradius', 'bohm_diffusion_coeff', 'omega_uh', 'omega_lh'])

tplot(['mms1_dis_energyspectr_omni_fast', 'mms1_dis_numberdensity_fast', 'mms1_fgm_b_gsm_srvy_l2_bvec', 'mms1_mec_r_gsm']) # Plot a multi-instrument figure

model = 'OpenGGCM_GM'

MW.Model_Variables(model)


The model accepts the standardized variable names listed below.
-----------------------------------------------------------------------------------
B_x : '['x component of magnetic field', 0, 'GSE', 'car', ['time', 'x', 'y', 'z'], 'nT']'
B_y : '['y component of magnetic field', 1, 'GSE', 'car', ['time', 'x', 'y', 'z'], 'nT']'
B_z : '['z component of magnetic field', 2, 'GSE', 'car', ['time', 'x', 'y', 'z'], 'nT']'
E_x : '['x component of electric field', 6, 'GSE', 'car', ['time', 'x', 'x', 'x'], 'mV/m']'
E_y : '['y component of electric field', 7, 'GSE', 'car', ['time', 'y', 'y', 'y'], 'mV/m']'
E_z : '['z component of electric field', 8, 'GSE', 'car', ['time', 'z', 'z', 'z'], 'mV/m']'
N_plasma : '['number density of plasma (hydrogen equivalent)', 12, 'GSE', 'car', ['time', 'x', 'y', 'z'], '1/cm**3']'
P_plasma : '['plasma pressure', 14, 'GSE', 'car', ['time', 'x', 'y', 'z'], 'pPa']'
eta : '['resistivity', 13, 'GSE', 'car', ['time', 'x', 'y', 'z'], 'm**2/s']'
j_x : '['current density, x component', 15, 'GSE', 'car', ['time', 'x', 'y', 'z'], 'muA/m**2']'
j_y : '['current density, y component', 16, 'GSE', 'car', ['time', 'x', 'y', 'z'], 'muA/m**2']'
j_z : '['current density, z component', 17, 'GSE', 'car', ['time', 'x', 'y', 'z'], 'muA/m**2']'
v_plasmax : '['x component of plasma velocity', 9, 'GSE', 'car', ['time', 'x', 'y', 'z'], 'km/s']'
v_plasmay : '['y component of plasma velocity', 10, 'GSE', 'car', ['time', 'x', 'y', 'z'], 'km/s']'
v_plasmaz : '['z component of plasma velocity', 11, 'GSE', 'car', ['time', 'x', 'y', 'z'], 'km/s']'

# Choose input values for ModelFlythrough function call file_dir = '/root/work/pydata/kamodo/' # full file path to where the model output data is stored variable_list = ['B_x','B_y', 'B_z'] # list of desired variable names coord_sys = "GSM-car" # GSM and cartesian because we're using x,y,z from above (see https://sscweb.gsfc.nasa.gov/users_guide/Appendix_C.shtml for a description of coordinate types) output_dir = '/root/work/output/ModelFlythrough/' # directory where the user wants the output stored file_name = 'OpenGGCMFlythrough' # what the user wants the output files to be named output_name = output_dir + file_name + ".csv" # output dir plus output file name without extension plot_coord = 'GSM' # coordinate system chosen for output plots # Convert to UTC timestamps from datetime import datetime, timezone sat_time = [time.replace(tzinfo=timezone.utc).timestamp() for time in epoch] # times range should be in the model data # Use MMS trajectory acquired from PySPEDAS/SpacePy sat_x = pos[:, 0] / 6378 sat_y = pos[:, 1] / 6378 sat_z = pos[:, 2] / 6378

# results = S.ModelFlythrough(model, file_dir, variable_list, sat_time, sat_x, sat_y, sat_z, # coord_sys, output_name=output_name, plot_coord=plot_coord)

results = S.WO.SF_read(output_name)

# Functionalize flythrough results kamodo_object = S.WO.Functionalize_SFResults(model,results) sat_fgm = pytplot.get_data("mms1_fgm_b_gsm_srvy_l2_bvec") sat_times = np.array(sat_fgm[0]) sat_B_x = np.array([b[0] for b in sat_fgm[1]]) sat_B_y = np.array([b[1] for b in sat_fgm[1]]) sat_B_z = np.array([b[2] for b in sat_fgm[1]]) # Functionalize MMS data kamodo_object = S.WO.Functionalize_TimeSeries(sat_times, 'B_xMMS', 'nT', sat_B_x, kamodo_object=kamodo_object) kamodo_object = S.WO.Functionalize_TimeSeries(sat_times, 'B_yMMS', 'nT', sat_B_y, kamodo_object=kamodo_object) kamodo_object = S.WO.Functionalize_TimeSeries(sat_times, 'B_zMMS', 'nT', sat_B_z, kamodo_object=kamodo_object)

Making an Executable Paper with the Python in Heliophysics Community to Foster Open Science and Improve Reproducibility

Abstract

1 Introduction

2 Method

2.3 Compare MMS Data to an Empirical Model

2.4 Augment Comparison with Plasma Parameters

3 Results

4 Discussion

Data Availability Statement

Conflict of Interest

Author Contributions

Funding

Acknowledgments

Footnotes

Appendix A. Team Member Descriptions

Appendix B. Deepnote Alternatives

Appendix C. PyHC Package Descriptions

References and Web Links