v2.3.27

pyrfu/mms/__init__.py

@@ -11,6 +11,7 @@
 from .lh_wave_analysis import lh_wave_analysis
 from .whistler_b2e import whistler_b2e
 from .get_feeps_omni import get_feeps_omni
+from .remove_imoms_background import remove_imoms_background
 from .remove_idist_background import remove_idist_background
 from .remove_edist_background import remove_edist_background
 from .psd_moments import psd_moments

pyrfu/mms/eis_skymap_combine_sc.py

@@ -103,7 +103,10 @@ def eis_skymap_combine_sc(skymaps, method: str = "mean"):
 
     out = xr.Dataset(out_dict)
 
-    out.attrs["energy_dminus"] = common_minu
-    out.attrs["energy_dplus"] = common_plus
+    out.attrs["delta_energy_minus"] = np.tile(common_minu, (len(ref_probe.time), 1))
+    out.attrs["delta_energy_plus"] = np.tile(common_plus, (len(ref_probe.time), 1))
+    out.attrs["energy0"] = common_energy[0, :]
+    out.attrs["energy1"] = common_energy[1, :]
+    out.attrs["species"] = "ions"
 
     return out

pyrfu/mms/fk_power_spectrum_4sc.py

@@ -3,6 +3,7 @@
 
 # Built-in imports
 import bisect
+import pdb
 
 # 3rd party imports
 import numpy as np
@@ -116,16 +117,20 @@ def fk_power_spectrum_4sc(
     times = e[0].time
     use_linear = df is not None
 
-    idx = time_clip(e[0].time, tints)
+    idx = time_clip(e[0].time, list(tints))
 
     # If odd, remove last data point (as is done in irf_wavelet)
     if len(idx) % 2:
         idx = idx[:-1]
 
     if use_linear:
-        cwt_options = dict(linear=df, returnpower=False, wavelet_width=5.36 * w_width)
+        cwt_options = dict(
+            linear=df, return_power=False, wavelet_width=5.36 * w_width, cut_edge=False
+        )
     else:
-        cwt_options = dict(nf=num_f, returnpower=False, wavelet_width=5.36 * w_width)
+        cwt_options = dict(
+            nf=num_f, return_power=False, wavelet_width=5.36 * w_width, cut_edge=False
+        )
 
     w = [wavelet(e[i], **cwt_options) for i in range(4)]
 
@@ -138,23 +143,23 @@ def fk_power_spectrum_4sc(
 
     fk_power = 0
     for i in range(4):
-        fk_power += w[i].data * np.conj(w[i].data) / 4
+        fk_power += (w[i].data * np.conj(w[i].data) / 4).astype(np.float64)
 
     n = int(np.floor(nt / cav) - 1)
-    pos_av = cav / 2 + np.arange(n) * cav
+    pos_av = cav / 2 + np.arange(n + 1) * cav
     av_times = times[pos_av.astype(np.int64)]
 
     b_avg = resample(b_avg, av_times)
 
     r = [resample(r[i], av_times) for i in range(4)]
 
-    cx12, cx13, cx14 = [np.zeros((n, num_f), dtype="complex128") for _ in range(3)]
-    cx23, cx24, cx34 = [np.zeros((n, num_f), dtype="complex128") for _ in range(3)]
+    cx12, cx13, cx14 = [np.zeros((n + 1, num_f), dtype="complex128") for _ in range(3)]
+    cx23, cx24, cx34 = [np.zeros((n + 1, num_f), dtype="complex128") for _ in range(3)]
 
-    power_avg = np.zeros((n, num_f), dtype="complex128")
+    power_avg = np.zeros((n + 1, num_f), dtype="complex128")
 
     for m, pos_avm in enumerate(pos_av):
-        lb, ub = [int(pos_avm - cav / 2 + 1), int(pos_avm + cav / 2)]
+        lb, ub = [int(pos_avm - cav / 2), int(pos_avm + cav / 2)]
 
         cx12[m, :] = np.nanmean(
             w[0].data[lb:ub, :] * np.conj(w[1].data[lb:ub, :]), axis=0
@@ -185,7 +190,7 @@ def fk_power_spectrum_4sc(
     th24 = np.arctan2(np.imag(cx24), np.real(cx24))
     th34 = np.arctan2(np.imag(cx34), np.real(cx34))
 
-    w_mat = 2 * np.pi * np.tile(w[0].frequency.data, (n, 1))
+    w_mat = 2 * np.pi * np.tile(w[0].frequency.data, (n + 1, 1))
 
     # Convert phase difference to time delay
     dt12, dt13, dt14, dt23, dt24, dt34 = [
@@ -215,24 +220,21 @@ def fk_power_spectrum_4sc(
     # Compute phase speeds
     r = [r[i].data for i in range(4)]
 
-    k_x, k_y, k_z = [np.zeros((n, num_f)) for _ in range(3)]
-
-    # Volumetric tensor with SC1 as center.
-    dr = np.reshape(np.hstack(r[1:]), (n, 3, 3)) - np.reshape(
-        np.tile(r[0], (1, 3)), (n, 3, 3)
-    )
-    dr = np.transpose(dr, [0, 2, 1])
-
-    # Delay tensor with SC1 as center.
-    # dT = np.reshape(np.hstack([dt2,dt3,dt4]),(N,num_f,3))
-    tau = np.dstack([dt2, dt3, dt4])
-
-    for ii in range(num_f):
-        m = np.linalg.solve(dr, np.squeeze(tau[:, ii, :]))
+    k_x, k_y, k_z = [np.zeros((n + 1, num_f)) for _ in range(3)]
 
-        k_x[:, ii] = 2 * np.pi * w[0].frequency[ii].data * m[:, 0]
-        k_y[:, ii] = 2 * np.pi * w[0].frequency[ii].data * m[:, 1]
-        k_z[:, ii] = 2 * np.pi * w[0].frequency[ii].data * m[:, 2]
+    for ii in range(n + 1):
+        dr = np.array(
+            [
+                r[1][ii, :] - r[0][ii, :],
+                r[2][ii, :] - r[0][ii, :],
+                r[3][ii, :] - r[0][ii, :],
+            ]
+        )
+        for jj in range(num_f):
+            m = np.linalg.solve(dr, [dt2[ii, jj], dt3[ii, jj], dt4[ii, jj]])
+            k_x[ii, jj] = 2 * np.pi * w[0].frequency[jj].data * m[0]
+            k_y[ii, jj] = 2 * np.pi * w[0].frequency[jj].data * m[1]
+            k_z[ii, jj] = 2 * np.pi * w[0].frequency[jj].data * m[2]
 
     k_x, k_y, k_z = [k / 1e3 for k in [k_x, k_y, k_z]]
 
@@ -273,21 +275,11 @@ def fk_power_spectrum_4sc(
 
         power_k_mag_f[nn, k_number] += np.real(power_avg[:, nn])
 
-    # power_k_x_f[power_k_x_f == 0] 	= np.nan
-    # power_k_y_f[power_k_y_f == 0] 	= np.nan
-    # power_k_z_f[power_k_z_f == 0] 	= np.nan
-    # power_k_mag_f[power_k_mag_f == 0] = np.nan
-
     power_k_x_f /= np.max(power_k_x_f)
     power_k_y_f /= np.max(power_k_y_f)
     power_k_z_f /= np.max(power_k_z_f)
     power_k_mag_f /= np.max(power_k_mag_f)
 
-    # power_k_x_f[power_k_x_f < 1.0e-6] 		= 1e-6
-    # power_k_y_f[power_k_y_f < 1.0e-6] 		= 1e-6
-    # power_k_z_f[power_k_z_f < 1.0e-6] 		= 1e-6
-    # power_k_mag_f[power_k_mag_f < 1.0e-6] 	= 1e-6
-
     frequencies = w[0].frequency.data
     idx_f = np.arange(num_f)
 
@@ -316,21 +308,11 @@ def fk_power_spectrum_4sc(
 
         power_k_perp_k_par[k_par_number, k_perp_number] += np.real(power_avg[:, nn])
 
-    # power_k_x_k_y[power_k_x_k_y == 0] = np.nan
-    # power_k_x_k_z[power_k_x_k_z == 0] = np.nan
-    # power_k_y_k_z[power_k_y_k_z == 0] = np.nan
-    # power_k_perp_k_par[power_k_perp_k_par == 0] = np.nan
-
     power_k_x_k_y /= np.max(power_k_x_k_y)
     power_k_x_k_z /= np.max(power_k_x_k_z)
     power_k_y_k_z /= np.max(power_k_y_k_z)
     power_k_perp_k_par /= np.max(power_k_perp_k_par)
 
-    # power_k_x_k_y(power_k_x_k_y < 1.0e-6) 				= 1e-6
-    # power_k_x_k_z(power_k_x_k_z < 1.0e-6) 				= 1e-6
-    # power_k_y_k_z(power_k_y_k_z < 1.0e-6) 				= 1e-6
-    # power_k_perp_k_par[power_k_perp_k_par < 1.0e-6] 	= 1e-6
-
     out_dict = {
         "k_x_f": (["k_x", "f"], power_k_x_f.T),
         "k_y_f": (["k_x", "f"], power_k_y_f.T),