API Reference

A variable v of a type derived from AbstractDataVariable should at least implement:

• Base.parent(v): the parent array of the variable

Optional:

• times(v): the timestamps of the variable
• units(v): the units of the variable
• meta(v): the metadata of the variable
• name(v): the name of the variable

DataSet <: AbstractDataSet

A concrete dataset with a name, data (parameters), and metadata.

Construct a DataSet from a name and data, with optional metadata.

DataSet(ld::LDataSet; kwargs...)

Create a concrete DataSet from a Dataset template with specified data.

See also: LDataSet

Instrument <: AbstractInstrument

Fields

• name: The name of the instrument
• metadata: Additional metadata
• datasets: Collection of datasets

Construct an Instrument from a dictionary.

keyword-based constructor

LDataSet <: AbstractDataSet

A template for generating datasets with parameterized naming patterns.

Fields

• format: Format string pattern for the dataset name
• data: Dictionary of variable patterns
• metadata: Additional metadata

Examples

using SPEDAS.MMS

# Access FPI dataset specification
lds = mms.datasets.fpi_moms

# Create a concrete dataset with specific parameters
ds = DataSet(lds; probe=1, data_rate="fast", data_type="des")

The format string and variable patterns use placeholders like {probe}, {data_rate}, which are replaced with actual values when creating a concrete DataSet.

Construct a LDataSet from a dictionary.

NoMetadata

Indicates an object has no metadata. But unlike using nothing, get, keys and haskey will still work on it, get always returning the fallback argument. keys returns () while haskey always returns false.

Project <: AbstractProject

A representation of a project or mission containing instruments and datasets.

Fields

• name: The name of the project
• metadata: Additional metadata
• instruments: Collection of instruments
• datasets: Collection of datasets

keyword-based constructor

Create a new product with the composed function

_getfield(v, name, default)

Return the field from a composite v for the given name, or the given default if no field is present.

See also: getfield.

abbr(p)

Get the abbreviation (abbr) of p.

Check if a string is in Day of Year format (YYYY-DDD).

@get(collection, key, default)

Short-circuiting version of get. See also @something.

Short-circuiting version of _getfield.

cotrans(A, in, out)
cotrans(A, out; in=get_coord(A))

Transform the data to the out coordinate system from the in coordinate system.

This function automatically choose between Julia's GeoCotrans (if available) and Fortran's IRBEM implementation.

References:

• IRBEM-LIB: compute magnetic coordinates and perform coordinate conversions (Documentation, IRBEM.jl)
• SPEDAS Cotrans

Coordinate systems, transformations, and geomagnetic field models.

Coordinate systems

• GEO: Geocentric geographic cartesian coordinate system (x [𝐋], y [𝐋], z [𝐋]).
• GSM: Geocentric Solar Magnetic (GSM) coordinate system.
• GSE: Geocentric Solar Ecliptic (GSE) coordinate system.
• GEI: Geocentric Equatorial Inertial (GEI) coordinate system.
• GDZ: Geodetic (GDZ) coordinate system (altitude [km], latitude [deg], longitude [deg]).
• MAG: Magnetic (MAG) coordinate system.
• SPH: Geocentric geographic spherical (SPH) coordinate system (r [𝐋], θ [deg], φ [deg]).
• J2000: J2000 (J2000) coordinate system.

Coordinate transformations

• geo2gei, gei2geo: Transform between GEO and GEI coordinate systems.
• geo2gsm, gsm2geo: Transform between GEO and GSM coordinate systems.
• gei2gsm, gsm2gei: Transform between GEI and GSM coordinate systems.
• gse2gsm, gsm2gse: Transform between GSE and GSM coordinate systems.

References

• Coordinate transformations between geocentric systems

International Geomagnetic Reference Field (IGRF)

The International Geomagnetic Reference Field (IGRF) is a standard mathematical description of the Earth's main magnetic field. It is used widely in studies of the Earth's deep interior, crust, ionosphere, and magnetosphere.

Functions

• igrf_B: Compute the geomagnetic field (IGRF-14, dipole model)

Examples

r, θ, φ = 6500., 30., 4.
t = Date(2021, 3, 28)
Br, Bθ, Bφ = igrf_B(r, θ, φ, t)

# Input position in geodetic coordinates, output magnetic field in East-North-Up (ENU) coordinates
Be, Bn, Bu = igrf_B(GDZ(0, 60.39299, 5.32415), t)

References

• IAGA - NOAA/NCEI
• IGRF-14 Evaluation

Elsewhere

• SatelliteToolboxGeomagneticField.jl: Models to compute the geomagnetic field (IGRF-13, dipole model)
• ppigrf: Pure Python code to calculate IGRF model predictions.
• geopack: Python code to calculate IGRF model predictions.

Geocentric Solar Magnetospheric (GSM)

X points sunward from Earth's center. The X-Z plane is defined to contain Earth's dipole axis (positive North).

Schmidt_normalization(l, m)

Compute Schmidt normalization factor for degree l and order m.

Reference: Geomagnetism and Schmidt quasi-normalization

calc_dipole_angle(g10, g11, h11)

Calculate dipole angle (θ, φ) and dipole strength (b0) from spherical harmonic coefficients g10, g11, h11.

θ: dipole tilt angle (radians) φ: dipole longitude/phase (radians) b0: dipole strength (nT)

calc_dipole_geo(time)

Compute dipole direction in GEO coordinates. [IRBEM]

calc_sun_gei(ra, dec)

Calculate sun direction vector in GEI given right ascension ra and declination dec.

calculate_gst(time)

Calculate Greenwich sidereal time (GST) in radians from the given time.

Reference: https://aa.usno.navy.mil/faq/GAST, https://github.com/JuliaAstro/AstroLib.jl/blob/main/src/ct2lst.jl

calculate_gst_alt(time)

Alternative implementation of Greenwich sidereal time calculation based on the reference algorithm from pyspedas.cotrans_tools.csundir_vect.

car2sph(x, y, z)

Convert (x, y, z) in Cartesian coordinate to (r, colat [deg], lon [deg]) in spherical coordinate.

cdipdir(time)

Compute dipole direction in GEO coordinates. [PySPEDAS]

References

• https://pyspedas.readthedocs.io/en/latest/coords.html#pyspedas.cotranstools.cotranslib.cdipdir

csundir(time)

Calculate the direction of the sun.

Returns

• gst: Greenwich mean sidereal time (radians)
• slong: Longitude along ecliptic (radians)
• sra: Right ascension (radians)
• sdec: Declination of the sun (radians)
• obliq: Inclination of Earth's axis (radians)

Notes

• This function calculates various parameters related to the sun's position needed for coordinate transformations.

gdz2sph(alt, lat, lon)

Convert (alt [km], lat [deg], lon [deg]) in Geodetic coordinate to Spherical geocentric coordinate.

gei2geo(x, t)

Transforms x vector from Geocentric Equatorial Inertial (GEI) to Geographic (GEO) coordinates at time t.

gei2gsm(x, t)

Transforms x vector from Geocentric Equatorial Inertial (GEI) to Geocentric Solar Magnetic (GSM) coordinates at time t.

gei2gsm_mat(time)

Compute the GEI to GSM transformation matrix.

First axis: sun direction (xS, yS, zS) Second axis: y-axis (cross product of dipole and sun, normalized) Third axis: z-axis (cross product of sun and y-axis, normalized)

geo2gei(x, t)

Transforms x vector from Geographic (GEO) to Geocentric Equatorial Inertial (GEI) coordinates at time t.

geo2gsm(x, t)

Transforms x vector from Geographic (GEO) to Geocentric Solar Magnetic (GSM) coordinates at time t.

get_dipole_terms(g, h)

Compute dipole parameters (θ, φ, x0, y0, z0, b0) from IGRF coefficients.

Returns a named tuple: (θ, φ, x0, y0, z0, b0)

get_igrf_coeffs(time)

Get IGRF-14 coefficients for a given time.

Similar to IRBEM implementation, but with higher precision (IRBEM uses year as the time unit).

gse2gsm(x, t)

Transforms x vector from Geocentric Solar Ecliptic (GSE) to Geocentric Solar Magnetic (GSM) coordinates at time t.

See also: gse2gsm_mat

gse2gsm(x, t)

Transforms x vector from Geocentric Solar Ecliptic (GSE) to Geocentric Solar Magnetic (GSM) coordinates at time t.

gse2gsm_mat(time)

Compute the GSE to GSM transformation matrix T3 = <- psi,X>, where the rotation angle psi is the GSE-GSM angle.

This transformation is a rotation in the GSE YZ plane from the GSE Z axis to the GSM Z axis.

gsm2gei(x, t)

Transforms x vector from Geocentric Solar Magnetic (GSM) to Geocentric Equatorial Inertial (GEI) coordinates at time t.

gsm2geo(x, t)

Transforms x vector from Geocentric Solar Magnetic (GSM) to Geographic (GEO) coordinates at time t.

gsm2gse(x, t)

Transforms x vector from Geocentric Solar Magnetic (GSM) to Geocentric Solar Ecliptic (GSE) coordinates at time t.

igrf_B(r, θ, φ, t; max_degree=IGRF_degree) -> (Br, Bθ, Bφ)

Calculate IGRF model components in geocentric coordinates (r [km], θ [rad], φ [rad]) at time t.

Parameters

• r: radius [km]
• θ: colatitude [rad]
• φ: longitude [rad], positive east
• maxdegree: highest degree of expansion (1 <= maxdegree <= 13)

igrf_B(𝐫::CoordinateVector{GDZ}, t; max_degree=IGRF_degree) -> (Be, Bn, Bu)

Calculate IGRF model components in east, north, up (ENU) coordinates for geodetic coordinates 𝐫 at time t.

igrf_Bd(r, θ, φ, t; max_degree=IGRF_degree) -> (Br, Bθ, Bφ)

Calculate IGRF model components in geocentric coordinates (r [km], θ [deg], φ [deg]) at time t.

Parameters

• r: radius [km]
• θ: colatitude [deg]
• φ: longitude [deg]
• maxdegree: highest degree of expansion (1 <= maxdegree <= 13)

Format into matrix form (l, m)

Debouncer

A struct that creates a debounced version of a function f.

The debounced function will only execute after a specified delay (in seconds) of inactivity. If called again before the delay expires, the timer resets.

DualAxisData(data1, data2; meta=nothing)

A type for handling dual-axis data where each field represents data for a different axis. The first field is plotted against the left y-axis and the second field against the right y-axis.

Fields

• data1: Data for the left y-axis
• data2: Data for the right y-axis
• metadata: Metadata for the data (e.g., title)

DualPlot is the plot type associated with plotting function dualplot. Check the docstring for dualplot for further information.

FunctionPlot is the plot type associated with plotting function functionplot. Check the docstring for functionplot for further information.

LinesPlot is the plot type associated with plotting function linesplot. Check the docstring for linesplot for further information.

MultiPlot is the plot type associated with plotting function multiplot. Check the docstring for multiplot for further information.

PanelPlot is the plot type associated with plotting function panelplot. Check the docstring for panelplot for further information.

SpecPlot is the plot type associated with plotting function specplot. Check the docstring for specplot for further information.

Convert the vector of vectors into a single vector of curves

Add labels to a figure, automatically searching for blocks to label.

Notes

• https://github.com/brendanjohnharris/Foresight.jl/blob/main/src/Layouts.jl#L2

Add labels to a grid of layouts

Notes

• See tag_facet in egg for reference

Only add legend when the axis contains multiple labels

Only add legend when the axis contains multiple labels

auto_figure_size(tas)

Calculate an appropriate figure size based on the number of plots in the list. Returns a tuple of (width, height) in pixels.

Get the default orientation for a legend based on the position

No docstring defined.

Plot type

The plot type alias for the dualplot function is DualPlot.

Attributes

color2 = (Makie.wong_colors())[end - 1] — No docs available.

linestyle2 = :dash — No docs available.

marker2 = :rect — No docs available.

plotfunc = scatterlines! — No docs available.

dualplot! is the mutating variant of plotting function dualplot. Check the docstring for dualplot for further information.

Setup the panel with both primary and secondary y-axes

Format datetime ticks with time on top and date on bottom.

No docstring defined.

Plot type

The plot type alias for the functionplot function is FunctionPlot.

Attributes

plotfunc = tplot_panel! — No docs available.

functionplot! is the mutating variant of plotting function functionplot. Check the docstring for functionplot for further information.

functionplot!(ax, f, tmin, tmax; kwargs...)

Interactive plot of a function f on ax for a time range from tmin to tmax

functionplot(gp, f, tmin, tmax; kwargs...)

Interactively plot a function over a time range on a grid position

No docstring defined.

Plot type

The plot type alias for the linesplot function is LinesPlot.

Attributes

labels = nothing — No docs available.

linesplot! is the mutating variant of plotting function linesplot. Check the docstring for linesplot for further information.

linesplot(gp, ta)

Plot a multivariate time series on a panel

Create and configure a secondary y-axis

No docstring defined.

Plot type

The plot type alias for the multiplot function is MultiPlot.

Attributes

plotfunc = plot! — No docs available.

multiplot! is the mutating variant of plotting function multiplot. Check the docstring for multiplot for further information.

multiplot!(ax, fs, t0, t1; plotfunc=plot2spec, kwargs...)

Specialized multiplot function for functionplot. Merge specs before plotting so as to cycle through them.

multiplot(gp, tas::MultiPlottable, args...; axis=(;), kwargs...)

Setup the panel on a position and plot multiple time series on it

No docstring defined.

Plot type

The plot type alias for the panelplot function is PanelPlot.

Attributes

color = @inherit linecolor — No docs available.

cycle = [:color] — No docs available.

plotfunc = plot! — No docs available.

panelplot! is the mutating variant of plotting function panelplot. Check the docstring for panelplot for further information.

Plot attributes for a time array (axis + labels)

Determine the plotting function for a given data type. Extend this for custom data types to integrate with the plotting system.

Determine the plotting function for a given data type. Extend this for custom data types to integrate with the plotting system.

Plot attributes for a time array (labels)

Set an attribute if all values are equal and non-empty

No docstring defined.

Plot type

The plot type alias for the specplot function is SpecPlot.

Attributes

specplot! is the mutating variant of plotting function specplot. Check the docstring for specplot for further information.

Plot heatmap of a time series on the same axis

specplot(gp, ta)

Plot a spectrogram on a panel

spectrogram_y_values(ta; check=false, center=false, transform=identity)

Get y-axis values from a spectrogram array. Can return either bin centers or edges. By default, return bin edges for better compatibility.

Arguments

• check: If true, check if values are constant along time
• center: If true, return bin centers instead of edges
• transform: Optional transform function for edge calculation (e.g., log for logarithmic bins)

Reference: Makie.edges

Add vertical lines to a plot

tplot(f, tas; legend=(; position=Right()), link_xaxes=true, link_yaxes=false, rowgap=5, transform=transform_pipeline, kwargs...)

Lay out multiple time series across different panels (rows) on one Figure / GridPosition f

If legend is nothing, no legend will be added to the plot. Otherwise, legend can be a NamedTuple containing options for legend placement and styling. By default, the time series are transformed via transform_pipeline, which is extensible via transform.

See also: tplot_panel, transform_pipeline, transform

tplot_panel!(ax, args...; kwargs...)

Generic entry point for adding plots to an existing axis ax.

Transforms the arguments to appropriate types and calls the plotting function. Dispatches to appropriate implementation based on the plotting trait of the transformed arguments.

tplot_panel(gp, args...; kwargs...)

Generic entry point for plotting different types of data on a grid position gp.

Transforms the arguments to appropriate types and calls the plotting function. Dispatches to appropriate implementation based on the plotting trait of the transformed arguments.

tplot_panel_s!(ax::Axis, data; kwargs...)

Plot data on an axis.

transform(args...; kwargs...)

Transform data into plottable format (e.g., DimArray).

Extend with transform(x::MyType) for custom types.

transform_pipeline(x)

Transform data for plotting with the following pipeline:

1. Custom transformations (transform(x))
2. String -> SpeasyProduct

See also: transform

Convert x to DateTime

Reference:

• https://docs.makie.org/dev/explanations/dim-converts#Makie.DateTimeConversion
• https://github.com/MakieOrg/Makie.jl/issues/442
• https://github.com/MakieOrg/Makie.jl/blob/master/src/dim-converts/dates-integration.jl

Julia-based Space Physics Environment Data Analysis Software

See the Documentation for more information.

SPEDAS.DEFAULTS

A global constant that holds default parameters:

• add_title::Bool defaults to false.
• add_colorbar::Bool defaults to true.
• delay : in seconds, the time interval between updates. Default is 0.25.
• resample::Int : the number of points to resample to. Default is 6070.

DiffQ(v, t; dims=1)

Difference quotient of v with respect to t.

Calculate the composite statistical error estimate for ⟨B·x₃⟩: |Δ⟨B·x₃⟩| = √(λ₃/(M-1) + (Δφ₃₂⟨B⟩·x₂)² + (Δφ₃₁⟨B⟩·x₁)²)

Parameters:

• λ₁, λ₂, λ₃: eigenvalues in descending order
• M: number of samples
• B: mean magnetic field vector
• x₁, x₂, x₃: eigenvectors

Constant Thickness Approach (CTA) for determining boundary normal and velocity. Based on the method described in Haaland et al., Annales Geophysicae, 2004.

ConstantVelocityApproach(positions, times, durations)

Given durations of the boundary crossings, calculate the thickness of the boundary

CVA(positions, times)

Constant Velocity Approach (CVA) for determining boundary normal and velocity. Solve timing equation: $D * m = Δts$

Parameters:

• positions: Positions of 4 spacecraft (4×3 array)
• times: Times of boundary crossing for each spacecraft

Discontinuity Analyzer (DA) for analyzing properties of discontinuities using multi-spacecraft measurements.

_linear_binedges(centers)

Calculate bin edges assuming linear spacing.

amap(f, a, b)

Apply a function f to the intersection of a and b.

https://github.com/rafaqz/DimensionalData.jl/issues/914

Return the angle between two vectors.

binedges(centers; transform=identity)

Calculate bin edges from bin centers.

• For linear spacing, edges are placed halfway between centers.
• For transformed spacing, edges are placed halfway between transformed centers.

Arguments

• transform: Function to transform the space (e.g., log for logarithmic spacing)

Example

centers = [1.0, 2.0, 3.0]
edges = binedges(centers)               # Returns [0.5, 1.5, 2.5, 3.5]
edges = binedges(centers, transform=log)  # Returns edges in log space

check_mva_eigen(F; r=5, verbose=false)

Check the quality of the MVA result.

If λ₁ ≥ λ₂ ≥ λ₃ are 3 eigenvalues of the constructed matrix M, then a good indicator of nice fitting LMN coordinate system should have λ₂ / λ₃ > r.

common_timerange(arrays; query=timeDimType)

Get the common time range (intersection) across multiple arrays. If there is no overlap, returns nothing.

cotrans(A, in, out)
cotrans(A, out; in=get_coord(A))

Transform the data to the out coordinate system from the in coordinate system.

This function automatically choose between Julia's GeoCotrans (if available) and Fortran's IRBEM implementation.

References:

• IRBEM-LIB: compute magnetic coordinates and perform coordinate conversions (Documentation, IRBEM.jl)
• SPEDAS Cotrans

current_density(B, V)

Calculate the current density time series from the magnetic field (B) and plasma velocity (V) time series.

Assume 1-D structure along the z-direction. Remember to transform the coordinates of B and V first (e.g. using mva

dimarrayify(x)

Convert x or values of x to DimArray(s).

dropna(da::DimArray, query)

Remove slices containing NaN values along dimensions other than query.

Transform matrix-like A to n×m shape

fac_mat(vec::AbstractVector; xref=[1.0, 0.0, 0.0])

Generates a field-aligned coordinate (FAC) transformation matrix for a vector.

Arguments

• vec: A 3-element vector representing the magnetic field

fill_gaps(times, data; resolution, margin)

Given a sorted vector of time stamps times and corresponding data values, this function inserts missing time stamps with a value of NaN if the gap between consecutive time stamps is larger than resolution + margin.

• If the gap is only slightly larger (within margin of the resolution), no gap is inserted.
• The function supports numeric times or DateTime (with appropriate resolution types).

Arguments

• times: Sorted vector of time stamps.
• resolution: The expected time difference between consecutive time stamps.
• margin: Allowed deviation from resolution before inserting missing time stamps.

Returns

A tuple (full_times, full_values) where:

• full_times is a vector containing all time stamps (original and inserted).
• full_values is a vector of data values with NaN for inserted gaps.

References

• https://pyspedas.readthedocs.io/en/latest/modules/pytplot/tplotmath/degap.html

find_continuous_timeranges(x, max_dt)

Find continuous time ranges for x (e.g. times or DimArray). max_dt is the maximum time gap between consecutive times.

find_outliers(A, [method, window]; dim = 1, kw...)

Find outliers in data A along the specified dim dimension.

Returns a Boolean array whose elements are true when an outlier is detected in the corresponding element of A.

The default method is :median (other option is :mean), which uses the median absolute deviation (MAD) to detect outliers. When the length of A is greater than 256, it uses a moving window of size 16.

See also: find_outliers_median, find_outliers_mean, isoutlier - MATLAB

find_outliers_mean(x::AbstractVector, window; threshold = 3)

Find outliers that are defined as elements more than three standard deviations from the mean.

This method is faster but less robust than find_outliers_median.

find_outliers_median(x, window; threshold=3)

Find outliers that are defined as elements more than threshold=3 times the scaled median absolute deviation (MAD) from the median.

When window is set to a integer, a moving window of that size is used to compute local MAD. Otherwise, global statistics are used.

References

• Median absolute deviation - Wikipedia

Group x into windows based on every and period.

interpolate_nans(da; interp=LinearInterpolation)

Interpolate only the NaN values in da along the specified dimension dims. Non-NaN values are preserved exactly as they are.

The default interpolation method interp is LinearInterpolation.

interpolate_outliers!(x, t, outliers)

Interpolate outliers in x using interpolation of neighboring non-outlier values.

lingradest(
    Bx1, Bx2, Bx3, Bx4,
    By1, By2, By3, By4,
    Bz1, Bz2, Bz3, Bz4,
    R1, R2, R3, R4
)

SPEDAS-argument-compatible version of lingradest.

lingradest(B1, B2, B3, B4, R1, R2, R3, R4)

Compute spatial derivatives such as grad, div, curl and curvature using reciprocal vector technique (linear interpolation).

Arguments

• B1, B2, B3, B4: 3-element vectors giving magnetic field measurements at each probe
• R1, R2, R3, R4: 3-element vectors giving the probe positions

Returns

A named tuple containing: • Rbary: Barycenter position • Bbc: Magnetic field at the barycenter • Bmag: Magnetic field magnitude at the barycenter • LGBx, LGBy, LGBz: Linear gradient estimators for each component • LD: Linear divergence estimator • LCB: Linear curl estimator • curvature: Field‐line curvature vector • R_c: Field‐line curvature radius

References

Based on the method of Chanteur (ISSI, 1998, Ch. 11).

• lingradest.pro
• lingradest.py

lingradest(B1::AbstractDimArray, args...)

Method for handling dimensional arrays. Takes AbstractDimArray inputs with a time dimension and returns a DimStack containing all computed quantities.

lingradest(B1::MatrixLike, args...)

Vectorized method for simplified usage. Returns a StructArray containing the results.

Convert slowness vector $𝐦 = 𝐧/V$ to normal vector and velocity

mva(V, B=V; kwargs...)

Rotate a timeseries V into the LMN coordinates based on the reference field B.

Arguments

• V: The timeseries data to be transformed, where each column represents a component
• B: The reference field used to determine the minimum variance directions, where each column represents a component

See also: mva_eigen, rotate

mva_eigen(B::AbstractMatrix; sort=(;), check=false) -> F::Eigen

Perform minimum variance analysis, returning Eigen factorization object F which contains the eigenvalues in F.values and the eigenvectors in the columns of the matrix F.vectors.

Set check=true to check the reliability of the result.

The kth eigenvector can be obtained from the slice F.vectors[:, k].

nt2ds(nt_arr, dim; fields=propertynames(first(nt_arr)))

Convert a NamedTuple array to a DimStack of DimArrays.

Vector rejection

Phase factor exp (i φ) satisfies the following equation

$\exp (4 i φ) = \exp (-2 i γ)$

where

$γ = \arctan (2 Re(𝐮)^𝐓 Im(𝐮) /(Re(𝐮)^2-Im(𝐮)^2))$

polarization(S0, S1, S2, S3)
polarization(S::StokesVector)

Compute the degree of polarization (p) from Stoke parameters or a Stokes vector.

Reference

• Wikipedia
• Stokes parameters

polarization(S)

Compute the degree of polarization (DOP) p^2 from spectral matrix S.

\[\begin{aligned} p^2 &= 1-\frac{(tr 𝐒)^2-(tr 𝐒^2)}{(tr 𝐒)^2-n^{-1}(tr 𝐒)^2} \\ &= \frac{n(tr 𝐒^2)-(tr 𝐒)^2}{(n-1)(tr 𝐒)^2} \end{aligned}\]

\[𝐑 = ∑_α (𝐫_α-𝐫_b) (𝐫_α-𝐫_b)' = ∑_α 𝐫_α 𝐫_α'-𝐫_b 𝐫_b'\]

with $𝐫_b = ∑_α 𝐫_α / N$ and N is the number of positions.

References

• Paschmann and Daly [2] Paschmann & Daly, 2008. Section 4.7

proj(a, b)

Vector projection of a vector a on (or onto) a nonzero vector b.

References: Wikipedia

See also: sproj, oproj

pspectrum(x::AbstractDimArray, spec::Spectrogram)
pspectrum(x::AbstractDimArray; nfft=256, noverlap=128, window=hamming)

Compute the power spectrum (time-frequency representation) of a time series using the short-time Fourier transform.

Returns a DimArray with frequency and original time dimensions.

See also: DSP.Spectrogram, DSP.stft

Reference

• Matlab

reciprocal_vector(rα, rβ, rγ, rλ)

Compute the reciprocal vector $𝒌_α$ for a vertex of a tetrahedron given the position vectors of all vertices.

The vertices (α, β, γ, λ) must form a cyclic permutation of (1, 2, 3, 4).

reciprocal_vector(rα, r0s::AbstractVector{<:AbstractVector})

Generalised reciprocal vector for N != 4

\[𝐪_α = 𝐑^{-1} 𝐫_α\]

See also: reciprocal_vector, position_tensor

reciprocal_vector(r_βα, r_βγ, r_βλ)

Compute the reciprocal vector $𝒌_α$ for a vertex of a tetrahedron given the relative position vectors.

\[𝒌_α = \frac{𝐫_{βγ} × 𝐫_{βλ}}{𝐫_{βα} ⋅ (𝐫_{βγ} × 𝐫_{βλ})}\]

where $𝐫_{αβ} = r_β - r_α$ are relative position vectors.

References

• Multi-spacecraft analysis methods revisited : 4.3 Properties of reciprocal vectors

Compute the set of reciprocal vectors {$𝒌_α$}, which is also called the reciprocal base of the tetrahedron.

See also: reciprocal_vector

Rectify the time step of a DimArray to be uniform.

replace_outliers!(A, s, [find_method, window]; kwargs...)

Finds outliers in A and replaces them with s (by default: NaN).

See also: find_outliers, filloutliers - MATLAB

replace_outliers!(A, method, [find_method, window]; kwargs...)
replace_outliers!(A, method, outliers; kwargs...)

Replaces outliers in A with values determined by the specified method.

Outliers can be detected using find_outliers with optional find_method and window parameters or specified directly as a Boolean array outliers.

method can be one of the following:

• :linear: Linear interpolation of neighboring, nonoutlier values
• :previous: Previous nonoutlier value
• :next: Next nonoutlier value
• :nearest: Nearest nonoutlier value

See also: filloutliers - MATLAB

replace_outliers(A; args...; kw...)

Non-mutable version of replace_outliers!.

resample(arr, n=DEFAULTS.resample; dim=1, verbose=false)

Resample an array along the dimension dim to n points. If the original length is less than or equal to n, the original array is returned unchanged.

rotate(ts::AbstractMatrix, mat::AbstractMatrix)

Coordinate-aware transformation of vector/matrix by rotation matrix(s) mat(s). Assume ts is a matrix of shape (n, 3).

Set the coordinate system.

Updates legend names and axis labels if they include the coordinate system. Also updates the dimension name if it contains the coordinate system.

Reference:

• https://pyspedas.readthedocs.io/en/latest/modules/pytplot/dataattgetterssetters.html#set_coords

smooth(da::AbstractDimArray, window; dim=Ti, suffix="_smoothed", kwargs...)

Smooths a time series by computing a moving average over a sliding window.

The size of the sliding window can be either:

• Quantity: A time duration that will be converted to number of samples based on data resolution
• Integer: Number of samples directly

Arguments

• dims=Ti: Dimension along which to perform smoothing (default: time dimension)
• suffix="_smoothed": Suffix to append to the variable name in output
• kwargs...: Additional arguments passed to RollingWindowArrays.rolling

smooth_spectral_matrix!(S_smooth, S, aa)

In-place version of smooth_spectral_matrix that writes results to a pre-allocated array.

smooth_spectral_matrix(S, aa)

Smooth the spectral matrix $S(f)$ by applying a weighted average over frequency. The smoothing uses a symmetric window aa (for example, a Hamming window) of length M.

Arguments

• S: Spectral matrix array of size $N_{freq}, n, n$ where n is the number of components.
• aa: Weighting vector of length M.

spectral_matrix(Xf)

Compute the spectral matrix $S$ defined by

\[S_{ij}(f) = X_i(f) X_j^*(f),\]

where $X_i(f)$=Xf[f, i] is the FFT of the i-th component and * denotes complex conjugation.

spectral_matrix(X, window)

Compute the spectral matrix $S(f)$ given the time series data X.

Returns a 3-D array of size $N_{freq}, n, n$, where $N_{freq} = \lfloor N/2 \rfloor$ and n is the dimensionality (number of components).

Arguments

• X: Matrix where each column is a component of the multivariate time series, or a vector of vectors.
• window: A window function (optional). If not provided, a rectangular window (no windowing) is used.

Scalar projection

tclip(d, t0, t1)

Clip a dimension or DimArray to a time range [t0, t1].

tclips(xs...; trange=common_timerange(xs...))

Clip multiple arrays to a common time range trange.

If trange is not provided, automatically finds the common time range across all input arrays.

tcross(x, y; dims=TimeDim, stack=nothing)

Compute the cross product of two (arrays of) vectors along the dims dimension.

References:

• https://docs.xarray.dev/en/stable/generated/xarray.cross.html

tderiv(data; dims = Ti)

Compute the time derivative of data.

See also: deriv_data - PySPEDAS

tderiv(data, times; dims = 1)

Compute the time derivative of data with respect to times.

tdot(x, y; dims=TimeDim)

Dot product of two arrays x and y along the dims dimension.

Calculate tetrahedron quality factors

tfilter(da, Wn1, Wn2=samplingrate(da) / 2; designmethod=Butterworth(2))

By default, the max frequency corresponding to the Nyquist frequency is used.

References

• https://docs.juliadsp.org/stable/filters/
• https://www.mathworks.com/help/signal/ref/filtfilt.html
• https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.filtfilt.html

Issues

• DSP.jl and Unitful.jl: https://github.com/JuliaDSP/DSP.jl/issues/431

tgroupby(x::AbstractDimArray, every, period = every, start_by = :window)

Returns a vector of AbstractDimArrays grouped by time intervals defined by every and period.

tinterp(A, B, interp=LinearInterpolation)

Interpolate A to times in B

tinterp(A, t; interp=LinearInterpolation)

Interpolate time series A at time point(s) t. Returns interpolated value for single time point or DimArray for multiple time points.

tinterp_nans(da::AbstractDimArray; query=timeDimType, kwargs...)

Interpolate only the NaN values in da along the specified dimensions query. Non-NaN values are preserved exactly as they are.

See also interpolate_nans

tmask!(da, t0, t1)
tmask!(da, it::Interval)
tmask!(da, its)

Mask all data values within the specified time range(s) (t0, t1) / it / its with NaN.

tmask(da, args...; kwargs...)

Non-mutable version of tmask!. See also tmask!.

https://github.com/JuliaLang/julia/issues/54542

tmean(x; query=TimeDim)

Calculate the mean of x along the query dimension (default: TimeDim).

It returns a value if x is a vector along the query dimension, otherwise returns a DimArray with the specified dimension dropped.

tnorm_combine(x; dims=timedim(x), name=:magnitude)

Calculate the norm of each slice along query dimension and combine it with the original components.

tresample(da::DimArray, n=DEFAULTS.resample; dimtype=Ti)

Resample a DimArray specifically along its dimension of type dimtype to n points. Throws an error if no dimension of type dimtype is found in the array.

tshift(x; dim=TimeDim, t0=nothing, new_dim=nothing)

Shift the dim of x by t0.

tsplit(da::AbstractDimArray, dim=Ti)

Splits up data along dimension dim.

tstack(vectors::AbstractVector{<:AbstractVector{T}})

Stack a time series of vectors into a matrix.

By default, each row in the output matrix represents a time point from the input vector of vectors.

tsubtract(x, f=nanmedian; dims=timedim(x))

Subtract a statistic (default function f: nanmedian) along dimensions (default: time dimension) from x.

tview(d, t0, t1)

View a dimension or DimArray in time range [t0, t1].

twavpol(x)

A convenience wrapper around wavpol that works with DimensionalData arrays.

It automatically extracts the time dimension and returns the results as a DimStack with properly labeled dimensions.

Wave normal angle is the angle between (wnx, wny) and the vertical |wnz| Use the imaginary parts of off-diagonals. Define:$A = Im(S₁₂), B = Im(S₁₃), C = Im(S₂₃)$

wavpol(X, fs=1; nfft=256, noverlap=div(nfft, 2), smooth_t=_smooth_t(nfft), smooth_f=hamming(3), nbuffers=Threads.nthreads())

Perform polarization analysis of n-component time series data.

For each FFT window (with specified overlap), the routine:

1. Applies a time-domain window function and computes the FFT to construct the spectral matrix $S(f)$
2. Applies frequency smoothing using a window function
3. Computes wave parameters: power, degree of polarization, wave normal angle, ellipticity, and helicity

The analysis assumes the data are in a right-handed, field-aligned coordinate system (with Z along the ambient magnetic field).

Arguments

• X: Matrix where each column is a component of the multivariate time series
• fs: Sampling frequency (default: 1)

Keywords

• nfft: Number of points for FFT (default: 256)
• noverlap: Number of overlapping points between windows (default: nfft÷2)
• smooth_t: Time domain window function (default: Hann window)
• smooth_f: Frequency domain smoothing window (default: 3-point Hamming window)
• nbuffers: Number of pre-allocated buffers for parallel processing (default: number of threads)

Returns

A named tuple containing:

• indices: Time indices for each FFT window
• freqs: Frequency array
• power: Power spectral density, normalized by frequency bin width and window function
• degpol: Degree of polarization [0,1]
• waveangle: Wave normal angle [0,π/2]
• ellipticity: Wave ellipticity [-1,1], negative for left-hand polarized
• helicity: Wave helicity

See also: polarization, wave_normal_angle, wpol_helicity

window_bf_sizes(window)

Converts a window specification to backward and forward window sizes.

When window is a positive integer scalar, the window is centered about the current element and contains window-1 neighboring elements. If window is even, then the window is centered about the current and previous elements.

wpol_helicity(S, waveangle)

Compute helicity and ellipticity for a single frequency.

Arguments

• S: Spectral matrix for a single frequency, size (3,3)
• waveangle: Wave normal angle for this frequency

Returns

• helicity: Average helicity across the three components
• ellipticity: Average ellipticity across the three components

permutedims is needed for series in Makie

Δφij(λᵢ, λⱼ, λ₃, M)

Calculate the phase error between components i and j according to: |Δφᵢⱼ| = |Δφⱼᵢ| = √(λ₃/(M-1) * (λᵢ + λⱼ - λ₃)/(λᵢ - λⱼ)²)

Parameters:

• λᵢ: eigenvalue i
• λⱼ: eigenvalue j
• λ₃: smallest eigenvalue (λ₃)
• M: number of samples

Convert angular frequency to frequency

Reference: https://www.wikiwand.com/en/articles/Angular_frequency

@load_project_config(file)

Load configuration from a file and export all key-value pairs as constants. The macro evaluates in the calling module's context.