Extracting the magnetopause from the MHD emissivity

Extracting the magnetopause from the MHD emissivity#

Instead of trying to smooth the LaTeP surface and deal with the discrete nature of the particle simulation, it would suffice to show that the maximum emissivity calculated from the MHD input coincides with the mean maximum curve of the particle model. This is to be expected, since the particles follow the motion dictated by the input fields and the position of the magnetopause should be able to be defined both by the field and the emissivity interchangeably. We load the cube and convert to spherical coordinates:

cube_MHD = readcube.Cube("MHD_OpenGGCM", step=0.05)
Q_MHD_GSE = cube_MHD.Q_GSE()
R, Theta, Phi, Q_MHD_sph = cube_MHD.cartesian_to_spherical(rmax = 15, dtheta = 0.5, dphi = 0.5)

The XY and XZ slices of the MHD emissivity cube are plotted in spherical coordinates:

_images/b2e2efa585ed0d5626ed833489b8169e0128a529842fda50a06e29f90a42b981.png

To better understand the radial profile of the emissivity, we plot it versus the radius for different polar \(\theta\) directions denoted above. We notice similar peaks as in the LaTeP case but with much smoother behavior as well as moderate contribution from the cusps, particularly due to the mask that is applied to remove the artificial inner-magnetosphere emission at \(R \leq 4 \,RE\).

_images/3a0bed59ffd1c792fca0e78b9eb0a2bed02db0373d8299dd968edf7b53d7914f.png

By following the previous methods, we can extract the maximum intensity per radial direction and plot it as a surface. Comparing with the LaTeP case this returns a smoother surface and further points away from the cusps. We can also try to limit the radial extend of the cusps to above \(5\,RE\), and thus gain more information on the directions that are obscured by them.

Here we try to keep both the global maximum peak, and the second maximum. In the MHD case, the second peak is indeed detected in the outer boundary of the magnetopause, and can be utilized to construct a surface without the problem of obscuring some points by the curvature of the cusps. In the LaTeP case, the second peak jumps from one boundary to another and oftentimes overshoots for the directions that we do not have two physical maxima. This is caused by the “bumpy” behavior of the radial profiles, meaning that small disturbances along the first maximum will be counted as a higher peak than the second maximum. This also causes the overshoot of the subsolar distance when fitting empirical models to the LaTeP surface. By limiting the radial maxima as well as the magnitude, we have essentially forced it to jump to these bumps for certain directions, resulting in an even higher variability in the surface, than if we took the total maximum. This could potentially be mitigated by smoothing the radial profiles before extracting the maxima.

Moving back to the analysis of the MHD cube, we can utilize this smooth behavior to construct a surface where we keep both the first and second maxima, keeping both the information about the strong indentation and the outer boundary:

# We do not discard something in particular here - some masking on the second peak will be done after we find the indexes
threshold = 10
r_skip = 0
print(f"Discard radius below {R[0, r_skip, 0]} Re and Emissivity above {threshold}")

# Prepare arrays to store up to two peak indices for each (theta, phi)
ind_peaks_below_thresh = np.full((Q_MHD_sph.shape[0], Q_MHD_sph.shape[2], 2), -1, dtype=int)

# Loop over theta and phi
for i_theta in range(Q_MHD_sph.shape[0]):
    for i_phi in range(Q_MHD_sph.shape[2]):
        profile = Q_MHD_sph[i_theta, r_skip:, i_phi]  # profile along r

        # Find local maxima
        peaks, _ = find_peaks(profile)

        # Filter by threshold
        peaks_below = peaks[profile[peaks] < threshold]

        # Sort peaks by their value (descending)
        if len(peaks_below) > 0:
            sorted_peaks = peaks_below[np.argsort(profile[peaks_below])[::-1]]
            # Keep up to two peaks
            for j, peak in enumerate(sorted_peaks[:2]):
                ind_peaks_below_thresh[i_theta, i_phi, j] = peak + r_skip

# Calculate the rmax surface for all (theta, phi) and peaks, and discards the ones that clip to the MHD mask
rmax_surface = np.full_like(R[:, 0, :], np.nan)
for i_theta in range(ind_peaks_below_thresh.shape[0]):
    for i_phi in range(ind_peaks_below_thresh.shape[1]):
        for j in range(2):
            idx = ind_peaks_below_thresh[i_theta, i_phi,j]
            if idx>20:
                rmax_surface[i_theta, i_phi] = R[i_theta, idx, i_phi]

Discard radius below 1.0 Re and Emissivity above 10