Analysis Tools

Wavelet Analysis Tools

There are several tools in PySPEDAS for performing wavelet analysis on tplot variables or data arrays.

The wavelet98 routine implements wavelet transforms as described in Torrence & Compo (1998). It works with bare data arrays.

Reference: Torrence, C. and G. P. Compo, 1998: A Practical Guide to Wavelet Analysis. <I>Bull. Amer. Meteor. Soc.</I>, 79, 61-78.

pyspedas.wavelet98(y1, dt, s0=None, dj=None, j_scales=None, pad=0, mother='MORLET', param=-1, lag1=0.0, siglvl=0.95, recon=False, fft_theor=None, do_wave=True, do_daughter=False)[source]

Compute the wavelet transform of a 1D time series using the Torrence & Compo (1998) algorithm.

Parameters:

  • y1 (array_like) – Time series to be analyzed.

  • dt (float) – Time step of the time series (sampling time).

  • s0 (float, optional) – Smallest scale of the wavelet transform. Default is 2 * dt.

  • dj (float, optional) – Spacing between discrete scales of the wavelet transform. Default is 0.125. A smaller value will give better scale resolution, but it will be slower.

  • j_scales (int, optional) – Number of scales minus one to use in the wavelet transform. Scales range from s0 to s0 * 2^(j_scales * dj). Default is int(np.floor(np.log2(n * dt / s0) / dj)).

  • pad (int {0, 1, 2}, optional::) – pad = 0, no padding (default) pad = 1, pad the time series to the next power of 2 pad = 2, pad the time series to the next power of 4

  • mother (string {'MORLET', 'PAUL', 'DOG'}, optional) – A string giving the mother wavelet to use. Default is ‘MORLET’.

  • param (float, optional) –

    Parameter for the mother wavelet. Default is -1:

    For 'Morlet' this is k0 (wavenumber), default is 6.
For 'Paul' this is m (order), default is 4.
For 'DOG' this is m (m-th derivative), default is 2.

  • lag1 (float, optional) – Lag-1 autocorrelation of the time series used for SIGNIF levels. Default is 0.0.

  • siglvl (float, optional) – Significance level to use for the wavelet transform. Default is 0.95.

  • recon (bool, optional) – If True, reconstruct the time series from the wavelet transform. Default is False.

  • fft_theor (array_like, optional) – Theoretical background spectrum as a function of Fourier frequency. This will be smoothed by the wavelet function.

  • do_wave (bool, optional) – If True, return the wavelet transform. Default is True.

  • do_daughter (bool, optional) – If True, return the wavelet transform of the mother wavelet. Default is False.

Returns:

  • wave (array_like) – Wavelet transform of the time series. This is a complex array of dimensions (N,j_scales+1).

  • scale (array_like) – Scales of the wavelet transform.

  • period (array_like) – Periods of the wavelet transform.

  • signif (array_like) – Significance levels of the wavelet transform.

  • daughter (array_like) – Wavelet wavelets, if do_daughter is True. Else, an empty array.

  • y2 (array_like) – Reconstructed time series, if recon is True. Else, -1.

  • fft_theor_out (array_like) – Output theoretical background spectrum (smoothed by the wavelet function) that was used in the wavelet transform.

References

Torrence, C. and G. P. Compo, 1998: A Practical Guide to Wavelet Analysis. <I>Bull. Amer. Meteor. Soc.</I>, 79, 61-78.

The wavelet2 routine is a wrapper for wavelet98. It also works with bare data arrays.

pyspedas.wavelet2(y, dt, prange=None, frange=None, pad=False, period=None, dj=None, param=None)[source]

Wrapper for the Torrence & Compo wavelet98.py function. Uses Morlet mother wavelet.

Parameters:

  • y (ndarray) – Input data, 1D or 2D (n x d2)

  • dt (float) – Time step

  • prange (list or tuple, optional) – Period range [min, max]

  • frange (list or tuple, optional) – Frequency range [min, max] (overrides prange if set)

  • pad (bool, optional) – If True, pad the time series for FFT

  • period (array, optional) – Output: periods used (filled by wavelet)

  • dj (float, optional) – Scale spacing

  • param (float, optional) – Morlet parameter (w0)

Returns:

  • wave (ndarray) – Wavelet transform (n, jv+1, d2) (n = number of time points, jv = number of scales, d2 = 1 or 2)

  • period (ndarray) – Periods corresponding to each scale (jv+1)

The wavelet routine provides an interface to the external pywavelets package, and works with tplot variables. If the wavelet scales are unspecified, they are derived using an algorithm similar to that used by the IDL SPEDAS wav_data routine.

pyspedas.wavelet(names, newname=None, suffix='_pow', wavename='morl', scales=None, method='fft', sampling_period=1.0)[source]

Find the wavelet transformation of a tplot variable.

Parameters:

  • names (str/list of str) – List of pytplot names.

  • newname (str/list of str, optional) – List of new names for tplot variables. Default: If not given, then a suffix is applied.

  • suffix (str, optional) – A suffix to apply. Default: ‘_pow’.

  • wavename (str, optional) – The name of the continuous wavelet function to apply. Examples: ‘gaus1’, ‘morl’, ‘cmorlB-C’. Default: ‘morl’

  • scales (list of float, optional) – The wavelet scales to use. Default: None

  • method (str, optional) – Either ‘fft’ for frequency domain convolution, or ‘conv’ for numpy.convolve. Default: ‘fft’

  • sampling_period (float, optional) – The sampling period for the frequencies output. Default: 1.0

Return type:

A list of tplot variables that contain the wavelet power.

Example

>>> import numpy as np
>>> import pyspedas
>>> from pyspedas import time_float
>>> from pyspedas.analysis.wavelet import wavelet
>>> # Create a tplot variable that contains a wave.
>>> t = np.arange(4000.)
>>> y = np.sin(2*np.pi*t/32.)
>>> y2 = np.sin(2*np.pi*t/64.)
>>> y[1000:3000] = y2[1000:3000]
>>> var = 'sin_wav'
>>> time = time_float('2010-01-01') + 10*t
>>> pyspedas.store_data(var, data={'x':time, 'y':y})
>>> # Gaussian Derivative wavelets transformation.
>>> powervar = wavelet(var, wavename='gaus1')
>>> pvar = powervar[0]
>>> # Define plotting parameters and plot.
>>> pyspedas.options(pvar, 'colormap', 'jet')
>>> pyspedas.ylim(pvar, 0.001, 0.1)
>>> pyspedas.options(pvar, 'ylog', True)
>>> pyspedas.options(pvar, 'ytitle', pvar)
>>> pyspedas.tplot([var, pvar])

The wav_data routine is a Python implementation of the IDL SPEDAS wav_data routine. This is the routine to use if you want to reproduce results from IDL SPEDAS.

pyspedas.wav_data(varname, period=None, prange=None, frange=None, param=None, avg_period=6.0, nmorlet=None, tplot_prefix=None, magrat=False, per_axis=False, kolom=None, normval=None, normconst=1.0, get_components=False, tint_pow=None, fraction=False, rotmats=None, get_rotmats=False, wid=None, dimennum=None, rotate_pow=False, maxpoints=32768, rbin=None, cross1=False, cross2=False, trange=None, resolution=0)[source]

Python implementation of wav_data.pro. Computes wavelet transform of tplot variable, using the Torrence & Compo (1998) algorithm.

Stores results in pytplot variables with normalization, masking, smoothing, and options.

This function calls wavelet2.pro which sets various parameters for the wavelet transform in wavelet98.py.

Parameters:

  • varname (str) – Name of tplot variable to analyze

  • period (array-like, optional) – Array of periods to use for wavelet analysis

  • prange (array-like, optional) – Period range [min_period, max_period] in same units as time

  • frange (array-like, optional) – Frequency range [min_freq, max_freq] - converted to prange

  • param (float, optional) – Wavelet parameter (default=6 for Morlet wavelet)

  • avg_period (float, optional) – Averaging period for smoothing (default=6.0)

  • nmorlet (float, optional) – Number of Morlet wavelets, converted to param

  • tplot_prefix (str, optional) – Prefix for output tplot variable names

  • magrat (bool, optional) – Compute magnitude ratio analysis

  • per_axis (bool, optional) – Use period (True) or frequency (False) for y-axis

  • kolom (array-like, optional) – Kolmogorov spectrum for normalization

  • normval (array-like or float, optional) – Normalization values

  • normconst (float, optional) – Normalization constant (default=1.0)

  • get_components (bool, optional) – Store individual component power spectra

  • tint_pow (float, optional) – Power law for frequency normalization

  • fraction (bool, optional) – Normalize by total power fraction

  • rotmats (array-like, optional) – Precomputed rotation matrices for rotate_pow

  • get_rotmats (bool, optional) – Return rotation matrices

  • wid (array-like, optional) – Smoothing width array

  • dimennum (int, optional) – Select specific component of input variable

  • rotate_pow (bool, optional) – Perform field-aligned coordinate analysis

  • maxpoints (int, optional) – Maximum number of time points (default=32768)

  • rbin (int, optional) – Rebinning factor

  • cross1 (bool, optional) – Cross-correlation analysis type 1

  • cross2 (bool, optional) – Cross-correlation analysis type 2

  • trange (array-like, optional) – Time range for analysis

  • resolution (int, optional) – Time resolution reduction factor. If > 0, reduces time resolution of output data using reduce_tres function (default=0, no reduction)

Returns:

Dictionary with the following keys:

’wave’: ndarray, wavelet transform data ‘power’: ndarray, power spectral data ‘period’: ndarray, period/frequency scale values ‘yax’: ndarray, y-axis values (period or frequency) ‘time’: ndarray, time values ‘returned_rotmats’: ndarray or None, rotation matrices if get_rotmats=True and rotate_pow=True

Return type:

dict

The wave_signif routine computes significance levels for a wavelet transform.

pyspedas.wave_signif(Y, dt, scale, sigtest, **kwargs)[source]

Compute significance levels for wavelet transforms.

Parameters:

Yfloat or array_like

The time series, or, the VARIANCE of the time series. (If this is a single number, it is assumed to be the variance…)

dtfloat

Amount of time between each Y value, i.e. the sampling time.

scalearray_like

The vector of scale indices, from previous call to WAVELET.

sigtestint {0, 1, 2}

Type of significance test:

If 0 (the default), then just do a regular chi-square test,
   i.e. Eqn (18) from Torrence & Compo.
If 1, then do a "time-average" test, i.e. Eqn (23).
         In this case, DOF should be set to NA, the number
         of local wavelet spectra that were averaged together.
         For the Global Wavelet Spectrum, this would be NA=N,
         where N is the number of points in your time series.
If 2, then do a "scale-average" test, i.e. Eqns (25)-(28).
         In this case, DOF should be set to a
         two-element vector [S1,S2], which gives the scale
         range that was averaged together.
         e.g. if one scale-averaged scales between 2 and 8,
             then DOF=[2,8].
kwargsdict

Dictionary, containing optional additional parameters:

lag1: Lag-1 autocorrelation (default: 0.0)
siglvl: Significance level (default: 0.95)
mother: Mother wavelet ('MORLET', 'PAUL', 'DOG')
param: Mother wavelet parameter
dof: Degrees of freedom for significance tests
gws: Global wavelet spectrum (optional)
confidence: If True, return confidence intervals.

Returns:

signifarray_like

Significance levels for the wavelet power spectrum.

outputsdict

Dictionary containing optional output parameters:

'Cdelta': Reconstruction factor
'period': Vector of Fourier periods corresponding to scales
'fft_theor': Theoretical red-noise spectrum as function of period
'Savg': Scale average (for sigtest=2 only)
'Smid': Scale midpoint (for sigtest=2 only)

Available wavelet definitions

The following routines are the predefined wavelet functions used by the wavelet98 routines.

Morlet wavelet:

pyspedas.analysis.wavelet98.morlet(k0, scale, k)[source]

Morlet wavelet in the Fourier domain.

Parameters:

  • k0 (float) – Non-dimensional wavenumber of the Morlet wavelet. If set to -1, the function will use the default value of 6.0.

  • scale (float) – Wavelet scale.

  • k (array_like) – Array of angular wavenumbers (rad / time unit). Typically created as in wavelet98: k = [0, k_pos…, k_neg…].

Returns:

  • vmorlet (ndarray) – Complex-valued Morlet wavelet evaluated in the Fourier domain at the frequencies given by k. Shape matches k.

  • period (float) – Equivalent Fourier period corresponding to the provided scale.

  • coi (float) – Cone-of-influence value for this scale (in same units as period).

  • dofmin (int) – Minimum degrees of freedom for significance testing (no smoothing).

  • cdelta (float) – Reconstruction factor used for inverse wavelet reconstruction. May be -1 if undefined for the chosen k0.

  • psi0 (float) – Time-domain normalization constant of the mother wavelet.

Notes

Implementation follows Torrence & Compo (1998). The returned vmorlet is normalized such that total energy equals the signal length N (see inline code comment “total energy=N”).

Paul wavelet:

pyspedas.analysis.wavelet98.paul(m, scale, k)[source]

Paul wavelet in the Fourier domain.

Parameters:

  • m (float) – Order of the Paul wavelet. If set to -1, the function will use the default value of 4.0.

  • scale (float) – Wavelet scale.

  • k (array_like) – Array of angular wavenumbers (rad / time unit). Typically created as in wavelet98: k = [0, k_pos…, k_neg…].

Returns:

  • vpaul (ndarray) – Complex-valued Paul wavelet evaluated in the Fourier domain at the frequencies given by k. Shape matches k.

  • period (float) – Equivalent Fourier period corresponding to the provided scale.

  • coi (float) – Cone-of-influence value for this scale (in same units as period).

  • dofmin (int) – Minimum degrees of freedom for significance testing (no smoothing).

  • cdelta (float) – Reconstruction factor used for inverse wavelet reconstruction. May be -1 if undefined for the chosen m.

  • psi0 (float) – Time-domain normalization constant of the mother wavelet.

Notes

Implementation follows Torrence & Compo (1998). The returned vpaul is normalized such that total energy equals the signal length N (see inline code comment “total energy=N”).

Difference of gaussians wavelet:

pyspedas.analysis.wavelet98.dog(m, scale, k)[source]

Derivative of Gaussian (DOG) wavelet in the Fourier domain.

Parameters:

  • m (int or float) – Order of the DOG wavelet (m-th derivative). If set to -1, the function will use the default value of 2.

  • scale (float) – Wavelet scale.

  • k (array_like) – Array of angular wavenumbers (rad / time unit). Typically created as in this module: k = [0, k_pos…, k_neg…].

Returns:

  • gauss (ndarray) – Complex-valued DOG wavelet evaluated in the Fourier domain at the frequencies given by k. Shape matches k.

  • period (float) – Equivalent Fourier period corresponding to the provided scale.

  • coi (float) – Cone-of-influence value for this scale (in same units as period).

  • dofmin (int) – Minimum degrees of freedom for significance testing (no smoothing).

  • cdelta (float) – Reconstruction factor used for inverse wavelet reconstruction. May be -1 if undefined for the chosen m.

  • psi0 (float) – Time-domain normalization constant of the mother wavelet.

Notes

Implementation follows Torrence & Compo (1998). Returned gauss is normalized such that total energy equals the signal length when used consistently with the rest of the transform implementation in this file.

Generalized 3-D Particle Distribution Tools

The tools documented in this section are not intended to be called directly by PySPEDAS users; rather, they are provided as building blocks for mission-specific 3-D particle distribution tools. Mission-specific wrappers will generally be needed to load the particle data to be operated on, perform any calibration, sanitization, or other preliminary steps, then populate the data structures used by the general-purpose particle tools.

For documentation of mission-specific particle tools, see the “Mission Specific Tools” page.

Plasma Moments

This group of routines calculates plasma moments (density, velocity, fluxes, pressure tensors, etc.) from 3-D particle distributions (with two dimensions being azimuthal and elevation angles, and the third dimension representing energy bins).

Zijin Zhang has written a useful guide to the mathematical definitions of the various moment quantities, and how the calculations are discretized in order to apply them to observed particle distributions:

https://juliaspacephysics.github.io/VelocityDistributionFunctions.jl/dev/moment_note/

moments_3d

This routine takes a data structure containing the particle distribution function, and other information like angle and energy bin definitions and sizes, and returns a dictionary containing plasma moments generated from the particle distributions.

pyspedas.moments_3d(data_in, sc_pot=0, no_unit_conversion=False)[source]

Calculates plasma moments from 3D data structure

Parameters:

  • data_in (dict) –

    Particle data structure with entries:

    'charge'
'mass'
'energy'
'denergy'
'theta'
'dtheta'
'phi'
'dphi'
'bins'
'data'

  • sc_pot (float) – Spacecraft potential

  • no_unit_conversion (bool) – Flag indicating that datta is already in eflux and no unit conversion is required

Note

The calculations were mostly heisted from Davin Larson’s IDL SPEDAS version

Returns:

Dictionary of plasma moments with entries:

'density'
'flux'
'eflux'
'qflux'
'mftens'
'velocity'
'ptens'
'ttens'
'vthermal'
'avgtemp'
'magt3'
't3'
'symm'
'symm_theta'
'symm_phi'
'symm_ang'

Return type:

dict

Examples

spd_pgs_moments

Basically a wrapper around moments_3d

pyspedas.spd_pgs_moments(data_in, sc_pot=0)[source]

Calculates moments from a simplified particle data structure. Simply a wrapper for moments_3d right now

Parameters:

  • data_in (dict) – Particle data structure

  • sc_pot (float) – Spacecraft potential

Returns:

Dictionary containing moments

Return type:

dict

spd_pgs_moments_tplot

Converts a dictionary (as returned by moments_3d) to tplot variables

pyspedas.spd_pgs_moments_tplot(moments, x=None, prefix='', suffix='', coords='DSL', use_mms_sdc_units=False)[source]

Creates tplot variables from moments dictionaries

Parameters:

  • moments (dict) – Dictionary containing moments values returned by moments_3d

  • x (numpy.ndarray) – The x-axis (time) values

  • prefix (str) – Name prefix for the output variables

  • suffix (str) – Name suffix for the output variables

  • coords (str) – Coordinate system to set for non-field-aligned moments that have coordinates Default: ‘DSL’

  • use_mms_sdc_units (bool) – If True, convert pressure tensor and heat flux values and units to nPa and mW/m^2 respectively, for compatibility with MMS SDC products. Default: False

Returns:

List of tplot variables returned

Return type:

list of str

Other quantities derived from 3-D particle distributions

spd_pgs_do_fac

pyspedas.particles.spd_part_products.spd_pgs_do_fac.spd_pgs_do_fac(data_in, mat)[source]

Applies field aligned coordinate transformation to input data

Parameters:

  • data_in (dict) – Particle data structure to be rotated

  • mat (numpy.ndarray) – The 3x3 field-aligned rotation matrix to apply to the data

Returns:

Rotated particle data structure

Return type:

dict

spd_pgs_regrid

pyspedas.particles.spd_part_products.spd_pgs_regrid(data, regrid_dimen)[source]

Regrid a data dictionary

Parameters:

  • data (dict) – The dictionary containing the input data

  • regrid_dimen (list of int) – The lengths of the new angle and energy arrays to compute

Returns:

A dictionary containing the regridded data

Return type:

dict

Slices of 3-D particle distributions

This set of routines creates 1-D and 2-D slices through 3-D particle distributions.

slice1d_plot

This routine plots the values along the x or y axis of a 2-D slice.

pyspedas.slice1d_plot(the_slice, direction, value, xrange=None, yrange=None)[source]

Create 1D plot from a 2D particle slice.

If the ‘value’ argument is a scalar, this plots a cut through the distribution at the nearest point in the specified direction

If the ‘value’ argument is a two-element list, this sums over the values between the min and max of the list

Parameters:

  • the_slice (dict) – 2D slice returned by slice2d

  • direction (str) – Axis to plot - ‘x’ or ‘y’

  • value (float or list of float) – If direction is ‘x’, this is the y-value to create a 1D plot at (and vice versa) can also be a range of values, e.g., [-1000, 1000] to sum over the y-values from -1000 to +1000

slice2d

pyspedas.slice2d(dists, time=None, samples=None, window=None, center_time=False, interpolation='geometric', trange=None, resolution=None, rotation='xy', custom_rotation=None, subtract_bulk=False, energy=False, log=False, erange=None, smooth=None, thetarange=None, zdirrange=None, average_angle=None, sum_angle=None, mag_data=None, vel_data=None, sun_data=None, slice_x=None, slice_z=None)[source]

Returns an interpolated 2D slice of 3D particle data for plotting

Interpolation methods:

2D Interpolation:

Data points within the specified theta or z-axis range are projected onto the slice plane and linearly interpolated onto a regular 2D grid.

Geometric:

Each point on the plot is given the value of the bin it intersects. This allows bin boundaries to be drawn at high resolutions.

Input

dists: list of dicts

List of 3D particle data structures

Basic Keywords

trange: list of str or list of float

Two-element time range over which data will be averaged (optional)

time: str

Time at which the slice will be computed (optional)

samples: int

Numer of samples nearest to TIME to average (default 1)

window: int or float

Length in seconds from TIME over which data will be averaged.

center_time: bool

Flag denoting that TIME should be midpoint for window instead of beginning.

interpolation: str

Interpolation method to use; options: ‘geometric’ for geometric interpolation and ‘2d’ for 2D interpolation (described above)

Orientation Keywords

custom_rotation: str or np.ndarray

Applies a custom rotation matrix to the data. Input may be a 3x3 rotation matrix or a tplot variable containing matrices. If the time window covers multiple matrices they will be averaged.

rotation: str

Aligns the data relative to the magnetic field and/or bulk velocity. This is applied after the CUSTOM_ROTATION. (BV and BE are invariant between coordinate systems)

Use MAG_DATA keyword to specify magnetic field vector. Use VEL_DATA keyword to specify bulk velocity (optional).

‘BV’: The x-axis is parallel to B field; the bulk velocity defines the x-y plane ‘BE’: The x-axis is parallel to B field; the B x V(bulk) vector defines the x-y plane ‘xy’: (default) The x-axis is along the data’s x-axis and y is along the data’s y axis ‘xz’: The x-axis is along the data’s x-axis and y is along the data’s z axis ‘yz’: The x-axis is along the data’s y-axis and y is along the data’s z axis ‘xvel’: The x-axis is along the data’s x-axis; the x-y plane is defined by the bulk velocity ‘perp’: The x-axis is the bulk velocity projected onto the plane normal to the B field; y is B x V(bulk) ‘perp_xy’: The data’s x & y axes are projected onto the plane normal to the B field ‘perp_xz’: The data’s x & z axes are projected onto the plane normal to the B field ‘perp_yz’: The data’s y & z axes are projected onto the plane normal to the B field

mag_data: str

Name of tplot variable containing magnetic field data or 3-vector. This will be used for slice plane alignment and must be in the same coordinates as the particle data.

vel_data: str

Name of tplot variable containing the bulk velocity data or 3-vector. This will be used for slice plane alignment and must be in the same coordinates as the particle data. If not set the bulk velocity will be automatically calculated from the distribution (when needed).

Other Keywords

resolution: int

Integer specifying the resolution along each dimension of the slice (defaults: 2D interpolation: 150, geometric: 500)

smooth: int

An odd integer >=3 specifying the width of a smoothing window in # of points. Smoothing is applied to the final plot using a gaussian convolution. Even entries will be incremented, 0 and 1 are ignored.

energy: bool

Flag to plot data against energy (in eV) instead of velocity.

log: bool

Flag to apply logarithmic scaling to the radial measure (i.e. energy/velocity). (on by default if energy=True)

erange: list of float

Two element list specifying the energy range to be used in eV

thetarange: list of float

(2D interpolation only): angle range, in degrees (-90, +90), used to calculate the slice default = [-20, 20]; will override zdirrange

zdirrange: list of float

(2D interpolation only): Z-Axis range, in km/s, used to calculate slice. Ignored if called with THETARANGE.

subtract_bulk: bool

Flag to subtract the bulk velocity vector

rtype:

Dictionary containing 2D slice of 3D particle data

slice2d_plot

pyspedas.slice2d_plot(the_slice, xrange=None, yrange=None, zrange=None, colormap=None, olines=8, contours=False, plotsize=10, title=None, save_png=None, save_jpeg=None, save_svg=None, save_pdf=None, save_eps=None, dpi=None, display=True)[source]

Creates plots of 2D particle slices

Magnetic Null Finding

For missions such as MMS or Cluster, with at least four spacecraft in a relatively close tetrahedron-like configuration, measuring the magnetic field simultaneously at four distinct locations allows the calculation of field gradients in each field component along the X, Y, and Z directions (in other words, a Jacobian matrix). This information is sufficient to find the location of magnetic null points (where all three field components are zero), and infer the topology of the magnetic field at the null point.

pyspedas.find_magnetic_nulls_fote(positions=None, fields=None, smooth_fields=True, smooth_npts=10, smooth_median=True, scale_factor=1.0)[source]

Find magnetic null points, using the First Order Taylor Expansion (FOTE) method, from a set of four-point magnetic field observations.

Parameters:

  • positions (list of str) – A 4-element list of tplot variable names representing the probe positions

  • fields (list of str) – A 4-element list of tplot variable names representing the magnetic field measurements

  • smooth_fields (bool) – If True, perform boxcar averaging on the fields Default: True

  • smooth_npts (int) – Number of points to use in the boxcar smoothing, if smoothing is enabled.

  • smooth_median (bool) – If True and smoothing is enabled, use median filtering Default: True

  • scale_factor (float) – The scale factor passed to the lingradest routine to scale some of the distances Default: 1.0

Returns:

A list of tplot variables describing the nulls found:

'null_pos': Position of the null point, in the same coordinate system as imput positions
'null_bary_dist': Distance between the null point and the barycenter of the tetrahedron
'null_bary_dist_types': A composite variable more suitable for plotting, with the null to barycenter distances,
      superimposed with symbols representing the type of each null found
'null_sc_distances': The distances from the null to each of the four spacecraft
'null_fom': Figures of merit 'eta' and 'xi', roughly representing the confidence in the null location and null type.
      Lower is better, with values less than 0.4 denoting fairly reliable detection and classification
'null_typecode': The type of each null point found, with values from 0-6. See classify_null_type() for interpretation.
'max_reconstruction_error': The maximum error out of the four s/c, when using the calculated Jacobian and field
      at the barycenter to reconstruct the field vectors at each spacecraft.  Should be extremely close to zero.

Return type:

list of str

Notes

This routine uses tetrahedral interpolation to estimate the magnetic field (B0) and field gradients at the tetrahedron barycenter. From the field gradients we can construct a Jacobian matrix (J). The field in that region can be expressed, to a first order approximation, as

B = B0 + JV

where V is the position vector relative to the barycenter. At a magnetic null V_null, all components of B are zero, so we have JV_null = -B0, which can be solved to get V_null.

As long as the four field measurements all differ, a null point will always be found. For it to be credible, it should be in the neighborhood where the linear gradient approximation is expected to be valid, i.e. some smallish multiple of the tetrahedron size.

The topology of the field around the null can be inferred from the eigenvalues and eigenvectors of the estimated Jacobian. See the classify_null_type() function for details.

References

Fu, H. S., A. Vaivads, Y. V. Khotyaintsev, V. Olshevsky, M. André, J. B. Cao, S. Y. Huang, A. Retinò, and G. Lapenta (2015), How to find magnetic nulls and reconstruct field topology with MMS data?. J. Geophys. Res. Space Physics, 120, 3758–3782. doi: 10.1002/2015JA021082.

Paschmann, G., Daly, P. (1998), Analysis Methods for Multi-Spacecraft Data, ISSR

Example

>>> import pyspedas
>>> from pyspedas import tplot
>>> data = pyspedas.projects.mms.fgm(probe=[1, 2, 3, 4], trange=['2015-09-19/07:40', '2015-09-19/07:45'], data_rate='srvy', time_clip=True, varformat='*_gse_*', get_fgm_ephemeris=True)
>>> fields = ['mms'+prb+'_fgm_b_gse_srvy_l2' for prb in ['1', '2', '3', '4']]
>>> positions = ['mms'+prb+'_fgm_r_gse_srvy_l2' for prb in ['1', '2', '3', '4']]
>>> null_vars = pyspedas.find_magnetic_nulls_fote(fields=fields, positions=positions, smooth_fields=True,smooth_npts=10,smooth_median=True)
>>> tplot(null_vars)
pyspedas.classify_null_type(lambdas_in)[source]

Determine the topological type of a magnetic null, given the eigenvalues of the Jacobian matrix.

Parameters:

lambdas_in (An array of 3 complex-valued eigenvalues of the Jacobian matrix,) – in no particular order.

Returns:

  • int –

    An integer representing the null type:

    0: Unknown/unable to characterize
1: X-type null
2: O-type null (island or plasmoid)
3: A-type (radial) null
4: B-type (radial) null
5: A_s-type (spiral) null
6: B_s-type (spiral) null
7: X-type, degenerated from A type
8: X-type, degenerated from B type
9: O-type, degenerated from A_s type
10: O-type, degenerated from B_s type

  • See Table 1, Fu et al 2015 for details

References

Fu, H. S., A. Vaivads, Y. V. Khotyaintsev, V. Olshevsky, M. André, J. B. Cao, S. Y. Huang, A. Retinò, and G. Lapenta (2015), How to find magnetic nulls and reconstruct field topology with MMS data?. J. Geophys. Res. Space Physics, 120, 3758–3782. doi: 10.1002/2015JA021082.

Paschmann, G., Daly, P. (1998), Analysis Methods for Multi-Spacecraft Data, ISSR

pyspedas.lingradest(Bx1, Bx2, Bx3, Bx4, By1, By2, By3, By4, Bz1, Bz2, Bz3, Bz4, R1, R2, R3, R4, scale_factor=1000.0)[source]

Calculate magnetic field gradients, divergence, curl, and field line curvature from 4-point observations

Parameters:

  • Bx1, Bx2, Bx3, Bx4 (ndarray)

  • By1, By2, By3, By4 (ndarray)

  • Bz1, Bz2, Bz3, Bz4 (ndarray) – Components of the B-field from the four probes versus time

  • R1, R2, R3, R4 (ndarray) – Position vectors for the four probes versus time

  • scale_factor (float) – Scaling divisor to apply to internal distance calculations. Default: 1000.0

References

Linear Gradient/Curl Estimator technique: see Chanteur, ISSI, 1998, Ch. 11

Notes

Based on the IDL version (lingradest.pro), which was originally designed for Cluster by A. Runov (2003)

Returns:

A dictionary containing:

Rbary, # position of barycenter
dR1, dR2, dR3, dR4, # Distance of barycenter from each probe
Bxbc, Bybc, Bzbc, Bbc
LGBx, LGBy, LGBz,
LCxB, LCyB, LCzB, LD,
curv_x_B, curv_y_B, curv_z_B, RcurvB

Return type:

dict