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
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"
# 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
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)
%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')
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)
# 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)
kamodo_object.plot('B_xMMS', 'B_x')
kamodo_object.plot('B_yMMS', 'B_y')
kamodo_object.plot('B_zMMS', 'B_z')