Added wave data processing based on RBSP data

el_paso/processing/__init__.py

@@ -15,6 +15,9 @@
 from el_paso.processing.fold_pitch_angles_and_flux import fold_pitch_angles_and_flux
 from el_paso.processing.get_real_time_tipsod import get_real_time_tipsod
 from el_paso.processing.magnetic_field_utils import MagFieldVarTypes
+from el_paso.processing.clean_magnetometer_data import clean_magnetometer_variables
+from el_paso.processing.calculate_L_mlt_mlat_fce_eq import add_derived_params
+
 
 __all__ = [
     "MagFieldVarTypes",
@@ -31,4 +34,6 @@
     "fold_pitch_angles_and_flux",
     "get_real_time_tipsod",
     "magnetic_field_utils",
+    "clean_magnetometer_variables",
+    "add_derived_params",
 ]

el_paso/processing/clean_magnetometer_data.py

@@ -0,0 +1,51 @@
+import numpy as np
+
+
+def clean_magnetometer_variables(variables):
+    """
+    Clean magnetometer data using Bt only.
+
+    Steps:
+    1. Keep only samples where Bt > 0.
+    2. Remove spikes where |ΔBt| >= 500 nT between consecutive samples.
+    3. Apply the resulting mask consistently to all input variables.
+
+    Parameters
+    ----------
+    variables : dict
+        Dictionary of variable objects with 'Bt' key included.
+        Each variable must implement .get_data() -> np.ndarray and
+        .set_data(filtered_data, mode).
+
+    Returns
+    -------
+    dict
+        Updated dictionary with filtered data arrays applied to all variables.
+    """
+
+    # Extract the magnetic field magnitude
+    Bt_var = variables["Bt"]
+    Bt = np.asarray(Bt_var.get_data())
+
+    # Step 1: Filter for Bt > 0
+    mask_positive = Bt > 0
+
+    # Step 2: Compute deltaB (prepend 0 to match original length)
+    deltaB = np.diff(Bt, prepend=Bt[0])
+
+    # Step 3: Filter out spikes (absolute change less than 500 nT)
+    mask_spike = np.abs(deltaB) < 500
+
+    # Combined mask
+    mask = mask_positive & mask_spike
+
+    # Step 4: Apply consistent mask to all variables
+    for var in variables.values():
+        data = np.asarray(var.get_data())
+        if data.shape[0] != mask.shape[0]:
+            raise ValueError(
+                f"Data length mismatch for variable. Expected {mask.shape[0]}, got {data.shape[0]}"
+            )
+        var.set_data(data[mask], unit="same")
+
+    return variables

el_paso/processing/compute_L_mlt_mlat_fce_eq.py

@@ -0,0 +1,64 @@
+import numpy as np
+from astropy import units as u
+from el_paso.variable import Variable
+
+
+def add_derived_params(variables, Re=6371.0):
+    """
+    Compute derived magnetospheric parameters (L, MLT, MLAT, fce, fce_eq)
+    and add them to a dictionary of EL_PASO Variable objects.
+
+    Parameters
+    ----------
+    variables : dict
+        Dictionary of EL_PASO Variable objects. Must contain:
+        - "Bt" : magnetic field strength (nT)
+        - "Coordinates" : Cartesian coordinates [x, y, z] in km
+    Re : float, optional
+        Earth radius in km (default: 6371).
+
+    Returns
+    -------
+    dict
+        Updated dictionary with new Variable entries:
+        - "L" : L-shell parameter (dimensionless)
+        - "mlat" : magnetic latitude (degrees)
+        - "mlt" : magnetic local time (hours)
+        - "fce" : electron gyrofrequency (Hz)
+        - "fce_eq" : equatorial electron gyrofrequency (Hz)
+    """
+
+    # Extract data arrays
+    Bt = np.asarray(variables["Bt"].get_data())
+    coords = np.asarray(variables["Coordinates"].get_data())
+
+    # Normalize coordinates by Earth radius
+    x, y, z = coords[:, 0] / Re, coords[:, 1] / Re, coords[:, 2] / Re
+
+    # Compute derived spatial parameters
+    r_xy = np.hypot(x, y)
+    r = np.sqrt(x**2 + y**2 + z**2)
+    mlat_rads = np.arctan2(z, r_xy)
+    mlat = np.degrees(mlat_rads)
+
+    # Dipole geometry formulas
+    L = r / np.cos(mlat_rads) ** 2
+    mlt = np.degrees(np.arctan2(y, x)) / 15.0 + 12.0  # convert deg→hours, center at noon
+    mlt = np.mod(mlt, 24.0)  # wrap to 0–24 h range
+
+    # Physical constants
+    q = 1.602e-19  # electron charge (C)
+    me = 9.109e-31  # electron mass (kg)
+
+    # Electron gyrofrequency and equatorial version
+    fce = (q * Bt * 1e-9) / (2 * np.pi * me)
+    fce_eq = fce * (np.cos(mlat_rads) ** 6) / np.sqrt(1 + 3 * np.sin(mlat_rads) ** 2)
+
+    # Create EL_PASO Variable objects with proper units
+    variables["L"] = Variable(u.dimensionless_unscaled, data=L)
+    variables["mlat"] = Variable(u.deg, data=mlat)
+    variables["mlt"] = Variable(u.hourangle, data=mlt)
+    variables["fce"] = Variable(u.Hz, data=fce)
+    variables["fce_eq"] = Variable(u.Hz, data=fce_eq)
+
+    return variables

scripts/inspect_cdf_file.py

@@ -34,6 +34,9 @@ def inspect_cdf_file(file_path: str) -> None:
     for var in variable_names:
         var_attrs_full = cdf_file.varattsget(var)
         vdr_info = cdf_file.varinq(var)
+        if vdr_info.Last_Rec < 0:
+            print(f"{var}: no records")
+            continue
         var_data = cdf_file.varget(var)
 
         var_shape = var_data.shape  # type: ignore[reportAttributeAccessIssue]