Problem Set 3

Code

dir = "docs/courses/epss261/homework"
if isdir(dir)
    cd(dir)
    Pkg.activate(".")
    Pkg.resolve()
    Pkg.instantiate()
end

1 Minimum and Maximum Variance Analysis

Minimize the variance $σ^2$ Eq. (8.4) of Ch. 8 by writing out the derivatives with respect to the vector components. Show that the solution is given by Equation (8.7) which is an eigenvalue problem with orthogonal eigenvectors.

\[ \sigma^2=\frac{1}{M} \sum_{m=1}^M\left|\left(\boldsymbol{B}^{(m)}-\langle\boldsymbol{B}\rangle\right) \cdot \hat{\boldsymbol{n}}\right|^2 \tag{1}\]

\[ \sum_{\nu=1}^3 M_{\mu \nu}^B n_\nu=\lambda n_\mu \tag{2}\]

Expand the dot product in components:

\[ \left(\boldsymbol{B}^{(m)} - \langle \boldsymbol{B} \rangle\right) \cdot \hat{\boldsymbol{n}} = \sum_{\mu=1}^3 \left(B_\mu^{(m)} - \langle B_\mu \rangle\right) n_\mu, \]

so that

\[ \sigma^2 = \frac{1}{M} \sum_{m=1}^M \left(\sum_{\mu=1}^3 \left(B_\mu^{(m)} - \langle B_\mu \rangle\right) n_\mu \right)^2 = \sum_{\mu,\nu=1}^3 n_\mu n_\nu \left[ \frac{1}{M} \sum_{m=1}^M \left(B_\mu^{(m)} - \langle B_\mu \rangle\right) \left(B_\nu^{(m)} - \langle B_\nu \rangle\right) \right] \]

Define the magnetic variance matrix

\[ M_{\mu\nu}^B \equiv \left\langle B_\mu B_\nu \right\rangle - \langle B_\mu \rangle \langle B_\nu \rangle \]

so that

\[ \sigma^2 = \sum_{\mu,\nu=1}^3 n_\mu M_{\mu\nu}^B n_\nu. \]

The Lagrangian is given by:

\[ L(\hat{\boldsymbol{n}}, \lambda) = \sigma^2 - \lambda\left(\sum_{\mu=1}^3 n_\mu^2 - 1\right) = \sum_{\mu,\nu=1}^3 n_\mu M_{\mu\nu}^B n_\nu - \lambda\left(\sum_{\mu=1}^3 n_\mu^2 - 1\right). \]

Differentiate with Respect to the Components $n_\alpha$ equal to zero for each component

\[ \frac{\partial}{\partial n_\alpha} \left(\sum_{\mu,\nu} n_\mu M_{\mu\nu}^B n_\nu\right) - \frac{\partial}{\partial n_\alpha} \left[\lambda\left(\sum_{\mu=1}^3 n_\mu^2 - 1\right)\right] = 2\sum_{\nu=1}^3 M_{\alpha\nu}^B n_\nu - 2\lambda n_\alpha. \]

we get:

\[ 2\sum_{\nu=1}^3 M_{\alpha\nu}^B n_\nu - 2\lambda n_\alpha = 0. \]

This is exactly Equation (8.7): $\sum_{\nu=1}^3 M_{\mu \nu}^B n_\nu = \lambda n_\mu$. The equation is an eigenvalue problem for the symmetric $3\times 3$ matrix $M_{\mu\nu}^B$. The eigenvalues $\lambda$ are real, and the corresponding eigenvectors $\hat{\boldsymbol{n}}$ are orthogonal.

2 Electromagnetic ion cyclotron waves

Electromagnetic ion cyclotron waves are left-hand polarized waves excited by velocity space anisotropy (perpendicular gradients in velocity distribution function). They are typically below the ion cyclotron frequency for each species (H+, He+ or O+). MMS observations on 2015 Dec. 14, show that such waves were excited after the passage of a shock.

Code

using Speasy
spz = speasy

using Dates
using DimensionalData
using CairoMakie, SpacePhysicsMakie
using SPEDAS
using Unitful
using DSP

2.1 Spectrogram with cyclotron frequencies

[a] Obtain MMS 1 (or 2,3,4) data between 13:27 and 13:29 UT, from the FGM instrument in GSM coordinates, and plot the spectrogram of the X or Y component between 0-2Hz. Overlay the 3 cyclotron frequencies, as shown in Fig. 5 of the paper (for subset of time). [Note: you should create a low-pass filtered version of the average field at 5-10s window to compute the cyclotron frequencies with; for that use tsmooth2 (IDL) or pyspedas.interpol (PySPEDAS) after degap/clip/deflaging the data first].

https://cdaweb.gsfc.nasa.gov/misc/NotesM.html#MMS1_FGM_SRVY_L2

Code

t0 = "2015-12-14T13:10:00"
t1 = "2015-12-14T13:50:00"
tvars = "cda/MMS1_FGM_SRVY_L2/mms1_fgm_b_gsm_srvy_l2_clean"
da = get_data(tvars, t0, t1) |> DimArray
da = SPEDAS.setmeta(da, "LABLAXIS" => "B")

Code

using SignalAnalysis
using PlasmaFormulary

p = SPEDAS.pspectrum(da[:, 1])
begin
    f = tplot(p; colorrange=(1e-7, 1), alpha=0.7)
    ax = current_axis()
    ylims!(ax, 0, 2)

    B_ts = smooth(da, 5u"s")[:, 4] .* 1u"nT"

    ic_H = gyrofrequency.(B_ts, :"H+"; to_hz = true) ./ 1u"Hz"
    ic_He = gyrofrequency.(B_ts, :"He+"; to_hz = true) ./ 1u"Hz"
    ic_He2 = gyrofrequency.(B_ts, :"He2+"; to_hz = true) ./ 1u"Hz"
    ic_O = gyrofrequency.(B_ts, :"O+"; to_hz = true) ./ 1u"Hz"
    tplot_panel!.(current_axis(), [ic_H, ic_He, ic_He2, ic_O])
    f
end

2.2 Band-pass filtered data

[b] Plot the GSM time series of the band-pass filtered data between 0.2 – 2 Hz (again using tsmooth2/pyspedas.interpol to first low-pass filter the data at 0.2Hz, subtract it from the original to create the high-pass residual >0.2Hz, then low-pass filter that with a 2Hz window).

Code

time = DateTime("2015-12-14T13:20:00") .. DateTime("2015-12-14T13:30:00")
sda = da[time, 1:3]
# Band-pass filter the data between 0.2-2 Hz
# First step: Low-pass filter at 0.2 Hz
lp_02Hz = smooth(sda, 5u"s")  # 0.2 Hz ≈ 5 second window
# High-pass by subtracting low-pass from original
hp_02Hz = amap(-, sda, lp_02Hz)
# Second step: Low-pass filter at 2 Hz (0.5 second window)
bp_filtered = smooth(hp_02Hz, 0.5u"s")

tplot([lp_02Hz, hp_02Hz, bp_filtered])

Actually, the frequency response of an $M$ point moving average filter is not a perfect low-pass filter.

\[ H[f]=\frac{\sin (\pi f M)}{M \sin (\pi f)} \]

Code

using TimeseriesUtilities: tfilter
bp_filtered2 = tfilter(sda, 0.2u"Hz", 2u"Hz")
tplot([hp_02Hz, bp_filtered, bp_filtered2])

2.3 Minimum Variance Analysis

[c] Use minvar_matrix_make (IDL or pyspedas) applied on the band-pass filtered data to obtain the minimum variance matrix, on a sliding window. You can control the width and shift of the window (see keywords). Rotate the data in i,j,k (max,int,min) coord’s using tvector_rotate (in IDL / Py - SPEDAS; see how to use it from class examples, or from cribs / readthedocs). Plot the data in that coordinate system.

Code

time = DateTime("2015-12-14T13:20:00") .. DateTime("2015-12-14T13:30:00")
B_ts = bp_filtered[time, 1:3]
B_mva_ts = setmeta(mva(B_ts, B_ts), "LABL_PTR_1" => ["Bl", "Bm", "Bn"])

B_ts_2 = bp_filtered2[time, 1:3]
B_mva_ts_2 = setmeta(mva(B_ts_2, B_ts_2), "LABL_PTR_1" => ["Bl", "Bm", "Bn"])

tplot([B_ts, B_mva_ts, B_ts_2, B_mva_ts_2])

2.4 Minvar direction

[d] The code minvar_matrix_make allows you to output the eigenvectors and eigenvalues. Plot the eigenvalues and evaluate the confidence in the minvar direction. Plot the angle between the minvar direction and the direction of the ambient magnetic field (the 5-10s average above).

Code

using LinearAlgebra
using Statistics
using Dates

# Load magnetic field data around the shock crossing
t0 = DateTime("2008-09-05T15:30:00")
t1 = DateTime("2008-09-05T15:40:00")
b_data = get_data("cda/THA_L2_FGM/tha_fgs_gse", t0, t1) |> DimArray
B = b_data.data  # Get the magnetic field components

# Get eigenvalues and eigenvectors
F = mva_eigen(B)

# The minimum variance direction (shock normal)
n̂ = F.vectors[:, 3]  # Third eigenvector corresponds to minimum variance

# Calculate θBn (angle between B-field and shock normal)
B_avg = mean(B, dims=1)[1, :]
B_avg_norm = B_avg / norm(B_avg)
θBn = acosd(abs(dot(B_avg_norm, n̂)))

# Calculate Bn (shock-normal magnetic field component)
Bn = dot(B_avg, n̂)

println("Results of Minimum Variance Analysis:")
println("Shock normal (GSE coordinates): ", round.(n̂, digits=3))
println("θBn: ", round(θBn, digits=1), "°")
println("Bn: ", Bn)
check_mva_eigen(F; verbose=true)

2.5 Hodogram

[e] There are several bursts of EMIC wave power in the interval. Pick one clean burst with 5-10 cycles and plot the B-field hodogram in the plane of polarization (max,int). (Bandpass filtering first will help). Is it left-hand polarized as expected? [Note that the minvar code ensures a right-hand orthogonal system, but you may have to switch vectors around if the angle to B-field is not less than 90deg].

Code

interval = DateTime("2015-12-14T13:40:20") .. DateTime("2015-12-14T13:40:40")
let
    sda = da[interval, 1:3]
    bp_filtered = tfilter(sda, 0.8u"Hz", 2u"Hz")
    B_mva_ts = mva(bp_filtered, bp_filtered)
    f = Figure(; size=(1200, 800))
    Bx = B_mva_ts[:, 1]
    By = B_mva_ts[:, 2]
    Bz = B_mva_ts[:, 3]
    tplot(f[1, 1], [sda, bp_filtered, B_mva_ts])
    tplot(f[1, 2], [Bx, By, Bz]; link_yaxes=true)
    time = dims(Bx, 1).val
    color = (time .- time[1]) ./ (time[end] - time[1])
    lines(f[1, 3], Bx.data, By.data, color=color)
    f
end

Hodogram analysis indicated left handed polarization during the interval.

3 Shock crossings

Shock crossings were observed by ARTEMIS P1 (TH-B) in the solar wind on 2013-07-09 20:40UT and on 2014-06-09 16:58UT, both captured in Fast Survey and published in Fig. 3 (analysis results in Table 1 of Zhou et al., 202022). Pick one of the two to study.

Code

using LinearAlgebra
using Statistics
using DataFrames, TimeSeries

3.1 Shock normal

[a] Determine, using minimum variance analysis, the shock normal, the angle between the B-field and the shock normal, and the shock-normal magnetic field component, Bn.

The shock angle, $θ_{Bn}$, the angle between the directions of the upstream magnetic field and the shock normal;

Code

# Load magnetic field data around the shock crossing
t0 = DateTime("2013-07-09T20:37:00")  # Start a few minutes before the crossing
t1 = DateTime("2013-07-09T20:45:00")  # End a few minutes after
b_data = get_data("cda/THB_L2_FGM/thb_fgs_gse", t0, t1) |> DimArray
tplot(b_data; add_title=true)

Code

# Perform MVA analysis using the built-in function
F = mva_eigen(b_data)

# The minimum variance direction (shock normal)
n̂ = F.vectors[:, 3]  # Third eigenvector corresponds to minimum variance

# Calculate θBn (angle between B-field and shock normal)
B_up = vec(mean(b_data[1:10, :].data; dims=1))
θBn = SPEDAS.angle(B_up, n̂)

# Calculate Bn (shock-normal magnetic field component)
Bn = dot(B_up, n̂)

println("Results of Minimum Variance Analysis for ARTEMIS P1 shock crossing (2013-07-09):")
println("Shock normal (GSE coordinates): ", round.(n̂, digits=3))
println("θ_Bn: ", round(θBn, digits=1), "°")
println("Bn: ", Bn)
check_mva_eigen(F, verbose=true)

3.2 Minimum variance direction

[b] Note that there are significant waves downstream but not upstream, which is common for weak (subcritical) shocks. Starting from a small interval barely encompassing the shock (<1min in length), increase the interval by ~0.5 min and recalculate. Then increase again and recalculate. Do so for ~10-20 times until you reach ~10min encompassing the shock. The solution will initially fluctuate from one to the next choice of intervals, then stabilize, then start being jittery again because you are including the waves. Plot the minimum variance direction (two angles, elevation and azimuth, in GSE) and Bn, all as a function of the interval chosen. Pick and report the best solution, based on the shortest time interval choice when the solution becomes stable.

Code

azimuth(x) = rad2deg(atan(x[2], x[1]))
elevation(x) = rad2deg(asin(x[3]))

# Central time of the shock crossing for ARTEMIS P1 (2013-07-09)
t_start = DateTime("2013-07-09T20:39:00")
time_shock = DateTime("2013-07-09T20:39:30")

n = 20
t_ends = [time_shock + Second(30 * (i - 1)) for i in 1:n]
b_data = get_data("cda/THB_L2_FGM/thb_fgs_gse", t_start, t_ends[end]) |> DimArray
b_data.metadata[:units] = "nT"

n̂s = map(t_ends) do t_end
    b_data_subset = b_data[t_start..t_end, :]
    F = mva_eigen(b_data_subset)
    F.vectors[:, 3]
end

Bn = dot.(Ref(B_up), n̂s)
direction = DimArray(
    [azimuth.(n̂s) elevation.(n̂s)],
    (Ti(t_ends), Y([:azimuth, :elevation]));
    name=:direction
)
Bns = abs.(DimArray(Bn, Ti(t_ends); name=:Bn))
ΔBn = abs.(diff(Bns))
ΔBn = rebuild(ΔBn; name=:ΔBn)

tvars = [b_data, direction, [Bns, ΔBn]]
f, axes = tplot(tvars; add_title=true)

The solution becomes stable around “2013-07-09T20:43:00”, about 4 mins after the shock crossing. after “2013-07-09T20:45:00” the solution becomes jittery because the waves are included.

3.3 Error estimates

[c] The noise in Bn and the minimum variance directions (angular noise) are tabulated in equations 8.24 and 8.23 of Chapter 8. Plot the $ΔB_n$ and $Δφ$ separately. Check if they support your choice of interval.

The angular error estimates (in radians) is

\[ \begin{aligned} \left|\Delta \varphi_{i j}\right|=\left|\Delta \varphi_{j i}\right| & =\left\langle\left\langle\left(\Delta x_{i j}\right)^2\right\rangle\right\rangle^{1 / 2}=\left\langle\left\langle\left(\Delta x_{j i}\right)^2\right\rangle\right\rangle^{1 / 2} \\ & =\sqrt{\frac{\lambda_3}{(M-1)} \frac{\left(\lambda_i+\lambda_j-\lambda_3\right)}{\left(\lambda_i-\lambda_j\right)^2}}, \quad i \neq j \end{aligned} \]

The composite statistical error estimate for $\left\langle\mathbf{B} \cdot \mathbf{x}_3\right\rangle$ is

\[ \left|\Delta\left\langle\mathbf{B} \cdot \mathbf{x₃}\right\rangle\right|=\sqrt{\frac{\lambda_3}{M-1}+\left(\Delta \varphi_{32}\langle\mathbf{B}\rangle \cdot \mathbf{x}_2\right)^2+\left(\Delta \varphi_{31}\langle\mathbf{B}\rangle \cdot \mathbf{x}_1\right)^2} \]

See MinimumVarianceAnalysis.jl for more details.

Code

using SPEDAS.MinimumVarianceAnalysis: Δφij, B_x3_error

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

Convert a NamedTuple array to a DimStack of DimArrays.
"""
function nt2ds(nt_arr, dim; fields=propertynames(first(nt_arr)))
    DimStack([
        DimArray(getfield.(nt_arr, field), dim; name=field)
        for field in fields
    ])
end

function nt2ds(nt_arr; sym=:time)
    dim = Dim{sym}(getfield.(nt_arr, sym))
    # filter the time dimension
    fields = propertynames(first(nt_arr))
    fields = filter(field -> field != sym, fields)
    nt2ds(nt_arr, dim; fields)
end

res = map(t_ends) do t_end
    b_data_subset = b_data[t_start..t_end, :]
    F = mva_eigen(b_data_subset)
    M = size(b_data_subset, 1)
    B = mean(b_data_subset.data; dims=1)

    λ₁ = F.values[1]
    λ₂ = F.values[2]
    λ₃ = F.values[3]
    Δφ12 = Δφij(λ₁, λ₂, λ₃, M)
    Δφ13 = Δφij(λ₁, λ₃, λ₃, M)
    Δφ23 = Δφij(λ₂, λ₃, λ₃, M)
    B_x3_err = B_x3_error(F, M, B)
    (; time=t_end, Δφ12, Δφ13, Δφ23, B_x3_err)
end

errors = nt2ds(res)
tvars = [
    b_data,
    [errors.Δφ12, errors.Δφ13, errors.Δφ23],
    errors.B_x3_err,
]
tplot(tvars)

The error estimates support our choice of interval.

4 Reproducibility

This document was built using these direct dependencies

Code

using Pkg
Pkg.status()

--- title: Problem Set 3 number-sections: true engine: julia execute: eval: false --- ```{julia} dir = "docs/courses/epss261/homework" if isdir(dir) cd(dir) Pkg.activate(".") Pkg.resolve() Pkg.instantiate() end ``` ## Minimum and Maximum Variance Analysis > Minimize the variance $σ^2$ Eq. (8.4) of Ch. 8 by writing out the derivatives with respect to the vector components. Show that the solution is given by Equation (8.7) which is an eigenvalue problem with orthogonal eigenvectors. $$ \sigma^2=\frac{1}{M} \sum_{m=1}^M\left|\left(\boldsymbol{B}^{(m)}-\langle\boldsymbol{B}\rangle\right) \cdot \hat{\boldsymbol{n}}\right|^2 $$ {#eq-Chen84} $$ \sum_{\nu=1}^3 M_{\mu \nu}^B n_\nu=\lambda n_\mu $$ {#eq-Chen87} Expand the dot product in components: $$ \left(\boldsymbol{B}^{(m)} - \langle \boldsymbol{B} \rangle\right) \cdot \hat{\boldsymbol{n}} = \sum_{\mu=1}^3 \left(B_\mu^{(m)} - \langle B_\mu \rangle\right) n_\mu, $$ so that $$ \sigma^2 = \frac{1}{M} \sum_{m=1}^M \left(\sum_{\mu=1}^3 \left(B_\mu^{(m)} - \langle B_\mu \rangle\right) n_\mu \right)^2 = \sum_{\mu,\nu=1}^3 n_\mu n_\nu \left[ \frac{1}{M} \sum_{m=1}^M \left(B_\mu^{(m)} - \langle B_\mu \rangle\right) \left(B_\nu^{(m)} - \langle B_\nu \rangle\right) \right] $$ Define the magnetic variance matrix $$ M_{\mu\nu}^B \equiv \left\langle B_\mu B_\nu \right\rangle - \langle B_\mu \rangle \langle B_\nu \rangle $$ so that $$ \sigma^2 = \sum_{\mu,\nu=1}^3 n_\mu M_{\mu\nu}^B n_\nu. $$ The Lagrangian is given by: $$ L(\hat{\boldsymbol{n}}, \lambda) = \sigma^2 - \lambda\left(\sum_{\mu=1}^3 n_\mu^2 - 1\right) = \sum_{\mu,\nu=1}^3 n_\mu M_{\mu\nu}^B n_\nu - \lambda\left(\sum_{\mu=1}^3 n_\mu^2 - 1\right). $$ Differentiate with Respect to the Components $n_\alpha$ equal to zero for each component $$ \frac{\partial}{\partial n_\alpha} \left(\sum_{\mu,\nu} n_\mu M_{\mu\nu}^B n_\nu\right) - \frac{\partial}{\partial n_\alpha} \left[\lambda\left(\sum_{\mu=1}^3 n_\mu^2 - 1\right)\right] = 2\sum_{\nu=1}^3 M_{\alpha\nu}^B n_\nu - 2\lambda n_\alpha. $$ we get: $$ 2\sum_{\nu=1}^3 M_{\alpha\nu}^B n_\nu - 2\lambda n_\alpha = 0. $$ This is exactly Equation (8.7): $\sum_{\nu=1}^3 M_{\mu \nu}^B n_\nu = \lambda n_\mu$. The equation is an eigenvalue problem for the symmetric $3\times 3$ matrix $M_{\mu\nu}^B$. The eigenvalues $\lambda$ are real, and the corresponding eigenvectors $\hat{\boldsymbol{n}}$ are orthogonal.  ## Electromagnetic ion cyclotron waves > Electromagnetic ion cyclotron waves are left-hand polarized waves excited by velocity space anisotropy (perpendicular gradients in velocity distribution function). They are typically below the ion cyclotron frequency for each species (H+, He+ or O+). MMS observations on 2015 Dec. 14, show that such waves were excited after the passage of a shock. ```{julia} using Speasy spz = speasy using Dates using DimensionalData using CairoMakie, SpacePhysicsMakie using SPEDAS using Unitful using DSP ``` ### Spectrogram with cyclotron frequencies > [a] Obtain MMS 1 (or 2,3,4) data between 13:27 and 13:29 UT, from the FGM instrument in GSM coordinates, and plot the spectrogram of the X or Y component between 0-2Hz. Overlay the 3 cyclotron frequencies, as shown in Fig. 5 of the paper (for subset of time). [Note: you should create a low-pass filtered version of the average field at 5-10s window to compute the cyclotron frequencies with; for that use tsmooth2 (IDL) or pyspedas.interpol (PySPEDAS) after degap/clip/deflaging the data first]. https://cdaweb.gsfc.nasa.gov/misc/NotesM.html#MMS1_FGM_SRVY_L2 ```{julia} t0 = "2015-12-14T13:10:00" t1 = "2015-12-14T13:50:00" tvars = "cda/MMS1_FGM_SRVY_L2/mms1_fgm_b_gsm_srvy_l2_clean" da = get_data(tvars, t0, t1) |> DimArray da = SPEDAS.setmeta(da, "LABLAXIS" => "B") ``` ```{julia} using SignalAnalysis using PlasmaFormulary p = SPEDAS.pspectrum(da[:, 1]) begin f = tplot(p; colorrange=(1e-7, 1), alpha=0.7) ax = current_axis() ylims!(ax, 0, 2) B_ts = smooth(da, 5u"s")[:, 4] .* 1u"nT" ic_H = gyrofrequency.(B_ts, :"H+"; to_hz = true) ./ 1u"Hz" ic_He = gyrofrequency.(B_ts, :"He+"; to_hz = true) ./ 1u"Hz" ic_He2 = gyrofrequency.(B_ts, :"He2+"; to_hz = true) ./ 1u"Hz" ic_O = gyrofrequency.(B_ts, :"O+"; to_hz = true) ./ 1u"Hz" tplot_panel!.(current_axis(), [ic_H, ic_He, ic_He2, ic_O]) f end ``` ### Band-pass filtered data > [b] Plot the GSM time series of the band-pass filtered data between 0.2 – 2 Hz (again using tsmooth2/pyspedas.interpol to first low-pass filter the data at 0.2Hz, subtract it from the original to create the high-pass residual >0.2Hz, then low-pass filter that with a 2Hz window). ```{julia} time = DateTime("2015-12-14T13:20:00") .. DateTime("2015-12-14T13:30:00") sda = da[time, 1:3] # Band-pass filter the data between 0.2-2 Hz # First step: Low-pass filter at 0.2 Hz lp_02Hz = smooth(sda, 5u"s") # 0.2 Hz ≈ 5 second window # High-pass by subtracting low-pass from original hp_02Hz = amap(-, sda, lp_02Hz) # Second step: Low-pass filter at 2 Hz (0.5 second window) bp_filtered = smooth(hp_02Hz, 0.5u"s") tplot([lp_02Hz, hp_02Hz, bp_filtered]) ``` Actually, the frequency response of an $M$ point moving average filter is not a perfect low-pass filter. $$ H[f]=\frac{\sin (\pi f M)}{M \sin (\pi f)} $$ ```{julia} using TimeseriesUtilities: tfilter bp_filtered2 = tfilter(sda, 0.2u"Hz", 2u"Hz") tplot([hp_02Hz, bp_filtered, bp_filtered2]) ``` ### Minimum Variance Analysis > [c] Use minvar_matrix_make (IDL or pyspedas) applied on the band-pass filtered data to obtain the minimum variance matrix, on a sliding window. You can control the width and shift of the window (see keywords). Rotate the data in i,j,k (max,int,min) coord’s using tvector_rotate (in IDL / Py - SPEDAS; see how to use it from class examples, or from cribs / readthedocs). Plot the data in that coordinate system. ```{julia} time = DateTime("2015-12-14T13:20:00") .. DateTime("2015-12-14T13:30:00") B_ts = bp_filtered[time, 1:3] B_mva_ts = setmeta(mva(B_ts, B_ts), "LABL_PTR_1" => ["Bl", "Bm", "Bn"]) B_ts_2 = bp_filtered2[time, 1:3] B_mva_ts_2 = setmeta(mva(B_ts_2, B_ts_2), "LABL_PTR_1" => ["Bl", "Bm", "Bn"]) tplot([B_ts, B_mva_ts, B_ts_2, B_mva_ts_2]) ``` ### Minvar direction > [d] The code minvar_matrix_make allows you to output the eigenvectors and eigenvalues. Plot the eigenvalues and evaluate the confidence in the minvar direction. Plot the angle between the minvar direction and the direction of the ambient magnetic field (the 5-10s average above). ```{julia} using LinearAlgebra using Statistics using Dates # Load magnetic field data around the shock crossing t0 = DateTime("2008-09-05T15:30:00") t1 = DateTime("2008-09-05T15:40:00") b_data = get_data("cda/THA_L2_FGM/tha_fgs_gse", t0, t1) |> DimArray B = b_data.data # Get the magnetic field components # Get eigenvalues and eigenvectors F = mva_eigen(B) # The minimum variance direction (shock normal) n̂ = F.vectors[:, 3] # Third eigenvector corresponds to minimum variance # Calculate θBn (angle between B-field and shock normal) B_avg = mean(B, dims=1)[1, :] B_avg_norm = B_avg / norm(B_avg) θBn = acosd(abs(dot(B_avg_norm, n̂))) # Calculate Bn (shock-normal magnetic field component) Bn = dot(B_avg, n̂) println("Results of Minimum Variance Analysis:") println("Shock normal (GSE coordinates): ", round.(n̂, digits=3)) println("θBn: ", round(θBn, digits=1), "°") println("Bn: ", Bn) check_mva_eigen(F; verbose=true) ``` ### Hodogram > [e] There are several bursts of EMIC wave power in the interval. Pick one clean burst with 5-10 cycles and plot the B-field hodogram in the plane of polarization (max,int). (Bandpass filtering first will help). Is it left-hand polarized as expected? [Note that the minvar code ensures a right-hand orthogonal system, but you may have to switch vectors around if the angle to B-field is not less than 90deg].  ```{julia} interval = DateTime("2015-12-14T13:40:20") .. DateTime("2015-12-14T13:40:40") let sda = da[interval, 1:3] bp_filtered = tfilter(sda, 0.8u"Hz", 2u"Hz") B_mva_ts = mva(bp_filtered, bp_filtered) f = Figure(; size=(1200, 800)) Bx = B_mva_ts[:, 1] By = B_mva_ts[:, 2] Bz = B_mva_ts[:, 3] tplot(f[1, 1], [sda, bp_filtered, B_mva_ts]) tplot(f[1, 2], [Bx, By, Bz]; link_yaxes=true) time = dims(Bx, 1).val color = (time .- time[1]) ./ (time[end] - time[1]) lines(f[1, 3], Bx.data, By.data, color=color) f end ``` Hodogram analysis indicated left handed polarization during the interval. ## Shock crossings > Shock crossings were observed by ARTEMIS P1 (TH-B) in the solar wind on 2013-07-09 20:40UT and on 2014-06-09 16:58UT, both captured in Fast Survey and published in Fig. 3 (analysis results in Table 1 of Zhou et al., 202022). Pick one of the two to study. ```{julia} using LinearAlgebra using Statistics using DataFrames, TimeSeries ``` ### Shock normal > [a] Determine, using minimum variance analysis, the shock normal, the angle between the B-field and the shock normal, and the shock-normal magnetic field component, Bn. The shock angle, $θ_{Bn}$, the angle between the directions of the upstream magnetic field and the shock normal; ```{julia} # Load magnetic field data around the shock crossing t0 = DateTime("2013-07-09T20:37:00") # Start a few minutes before the crossing t1 = DateTime("2013-07-09T20:45:00") # End a few minutes after b_data = get_data("cda/THB_L2_FGM/thb_fgs_gse", t0, t1) |> DimArray tplot(b_data; add_title=true) ``` ```{julia} # Perform MVA analysis using the built-in function F = mva_eigen(b_data) # The minimum variance direction (shock normal) n̂ = F.vectors[:, 3] # Third eigenvector corresponds to minimum variance # Calculate θBn (angle between B-field and shock normal) B_up = vec(mean(b_data[1:10, :].data; dims=1)) θBn = SPEDAS.angle(B_up, n̂) # Calculate Bn (shock-normal magnetic field component) Bn = dot(B_up, n̂) println("Results of Minimum Variance Analysis for ARTEMIS P1 shock crossing (2013-07-09):") println("Shock normal (GSE coordinates): ", round.(n̂, digits=3)) println("θ_Bn: ", round(θBn, digits=1), "°") println("Bn: ", Bn) check_mva_eigen(F, verbose=true) ``` ### Minimum variance direction > [b] Note that there are significant waves downstream but not upstream, which is common for weak (subcritical) shocks. Starting from a small interval barely encompassing the shock (<1min in length), increase the interval by ~0.5 min and recalculate. Then increase again and recalculate. Do so for ~10-20 times until you reach ~10min encompassing the shock. The solution will initially fluctuate from one to the next choice of intervals, then stabilize, then start being jittery again because you are including the waves. Plot the minimum variance direction (two angles, elevation and azimuth, in GSE) and Bn, all as a function of the interval chosen. Pick and report the best solution, based on the shortest time interval choice when the solution becomes stable. ```{julia} azimuth(x) = rad2deg(atan(x[2], x[1])) elevation(x) = rad2deg(asin(x[3])) # Central time of the shock crossing for ARTEMIS P1 (2013-07-09) t_start = DateTime("2013-07-09T20:39:00") time_shock = DateTime("2013-07-09T20:39:30") n = 20 t_ends = [time_shock + Second(30 * (i - 1)) for i in 1:n] b_data = get_data("cda/THB_L2_FGM/thb_fgs_gse", t_start, t_ends[end]) |> DimArray b_data.metadata[:units] = "nT" n̂s = map(t_ends) do t_end b_data_subset = b_data[t_start..t_end, :] F = mva_eigen(b_data_subset) F.vectors[:, 3] end Bn = dot.(Ref(B_up), n̂s) direction = DimArray( [azimuth.(n̂s) elevation.(n̂s)], (Ti(t_ends), Y([:azimuth, :elevation])); name=:direction ) Bns = abs.(DimArray(Bn, Ti(t_ends); name=:Bn)) ΔBn = abs.(diff(Bns)) ΔBn = rebuild(ΔBn; name=:ΔBn) tvars = [b_data, direction, [Bns, ΔBn]] f, axes = tplot(tvars; add_title=true) ``` The solution becomes stable around "2013-07-09T20:43:00", about 4 mins after the shock crossing. after "2013-07-09T20:45:00" the solution becomes jittery because the waves are included. ### Error estimates > [c] The noise in Bn and the minimum variance directions (angular noise) are tabulated in equations 8.24 and 8.23 of Chapter 8. Plot the $ΔB_n$ and $Δφ$ separately. Check if they support your choice of interval. The angular error estimates (in radians) is $$ \begin{aligned} \left|\Delta \varphi_{i j}\right|=\left|\Delta \varphi_{j i}\right| & =\left\langle\left\langle\left(\Delta x_{i j}\right)^2\right\rangle\right\rangle^{1 / 2}=\left\langle\left\langle\left(\Delta x_{j i}\right)^2\right\rangle\right\rangle^{1 / 2} \\ & =\sqrt{\frac{\lambda_3}{(M-1)} \frac{\left(\lambda_i+\lambda_j-\lambda_3\right)}{\left(\lambda_i-\lambda_j\right)^2}}, \quad i \neq j \end{aligned} $$ The composite statistical error estimate for $\left\langle\mathbf{B} \cdot \mathbf{x}_3\right\rangle$ is $$ \left|\Delta\left\langle\mathbf{B} \cdot \mathbf{x₃}\right\rangle\right|=\sqrt{\frac{\lambda_3}{M-1}+\left(\Delta \varphi_{32}\langle\mathbf{B}\rangle \cdot \mathbf{x}_2\right)^2+\left(\Delta \varphi_{31}\langle\mathbf{B}\rangle \cdot \mathbf{x}_1\right)^2} $$ See [MinimumVarianceAnalysis.jl](https://juliaspacephysics.github.io/MinimumVarianceAnalysis.jl/dev/) for more details. ```{julia} using SPEDAS.MinimumVarianceAnalysis: Δφij, B_x3_error """ nt2ds(nt_arr, dim; fields=propertynames(first(nt_arr))) Convert a NamedTuple array to a DimStack of DimArrays. """ function nt2ds(nt_arr, dim; fields=propertynames(first(nt_arr))) DimStack([ DimArray(getfield.(nt_arr, field), dim; name=field) for field in fields ]) end function nt2ds(nt_arr; sym=:time) dim = Dim{sym}(getfield.(nt_arr, sym)) # filter the time dimension fields = propertynames(first(nt_arr)) fields = filter(field -> field != sym, fields) nt2ds(nt_arr, dim; fields) end res = map(t_ends) do t_end b_data_subset = b_data[t_start..t_end, :] F = mva_eigen(b_data_subset) M = size(b_data_subset, 1) B = mean(b_data_subset.data; dims=1) λ₁ = F.values[1] λ₂ = F.values[2] λ₃ = F.values[3] Δφ12 = Δφij(λ₁, λ₂, λ₃, M) Δφ13 = Δφij(λ₁, λ₃, λ₃, M) Δφ23 = Δφij(λ₂, λ₃, λ₃, M) B_x3_err = B_x3_error(F, M, B) (; time=t_end, Δφ12, Δφ13, Δφ23, B_x3_err) end errors = nt2ds(res) tvars = [ b_data, [errors.Δφ12, errors.Δφ13, errors.Δφ23], errors.B_x3_err, ] tplot(tvars) ``` The error estimates support our choice of interval. ## Reproducibility ::: {.callout-note collapse="true"} ## This document was built using these direct dependencies ```{julia} using Pkg Pkg.status() ``` :::