GitHub

Solar Energetic Particle Onset Analysis

Setup

First, let's load the required packages and set up the analysis environment:

using SolarEnergeticParticle
using Dates
using SpacePhysicsMakie, CairoMakie

The time range for the analysis is set as follows:

# Analysis time range (October 28-29, 2021 SEP event)
tmin = DateTime(2021, 10, 28)
tmax = DateTime(2021, 10, 29)

# Define background time range (before the event)
background_range = (DateTime(2021, 10, 28, 10, 0, 0), DateTime(2021, 10, 28, 12, 0, 0))
(Dates.DateTime("2021-10-28T10:00:00"), Dates.DateTime("2021-10-28T12:00:00"))

Load SEP Data

Load data from the SOHO ERNE-HED instrument:

data = get_data("SOHO_ERNE-HED_L2-1MIN", tmin, tmax)
tplot(data.PH)
Example block output

Perform Onset Analysis

Now we can perform the onset determination analysis using find_onset and find_peak functions.

# Select a specific channel
channel = 4
I = select_channel(data.PH, channel)
peak = find_peak(I)
onset = find_onset(I, background_range)
(onset_time = Dates.DateTime("2021-10-28T16:57:17.639"), bg_timerange = (Dates.DateTime("2021-10-28T10:00:00"), Dates.DateTime("2021-10-28T12:00:00")), bg_mean = 0.00010347501f0, bg_std = 0.000102109436f0)

Visualization

Create a comprehensive plot showing the onset analysis results:

f = tplot(I)

_tvspan!(ax, tmin, tmax; kwargs...) = vspan!(ax, tmin, Dates.value.(tmax); kwargs...)

function onsetplot!(ax, onset_time, bg_timerange, bg_mean, bg_std, peak_time = nothing, n = 2)
    # Add background region if specified
    onset_plot = tlines!(ax, onset_time, color = :red, label = "Onset $(round(onset_time, Second))")
    bg_plot = _tvspan!(ax, bg_timerange...; alpha = 0.2, color = :gray, label = "Background")
    bg_stat_plot1 = hlines!(bg_mean, linestyle = :dash, color = :cyan, label = "Mean of Background")
    bg_stat_plot2 = hlines!(bg_mean + n*bg_std, linestyle = :dash, color = :cyan, label = "Mean + $n Std of Background")
    plots = [onset_plot, bg_plot, bg_stat_plot1, bg_stat_plot2]

    if peak_time !== nothing
        peak_plot = tlines!(ax, peak_time, color = :orange, linewidth = 2, label = "Peak $(round(peak_time, Second))")
        push!(plots, peak_plot)
    end

    return plots
end

onsetplot!(f.axes[1], onset..., peak[1])
ylims!(1e-4, 5e-1)
axislegend(f.axes[1]; position = :lt)
f
Example block output

References