API Reference

This page contains the auto-generated documentation for the rxn_location package's core modules and utility functions.

Core Physics & Models

rxn_location.rx_model_funcs

get_bis(b_vec_1, b_vec_2)

Calculate the bisector field components of two magnetic field lines.

Parameters:

Name	Type	Description	Default
b_vec_1	array-like of shape (3,)	First input magnetic field vector.	required
b_vec_2	array-like of shape (3,)	Second input magnetic field vector.	required

Returns:

Name	Type	Description
bis_field_1	float	The magnitude of the bisected field line corresponding to the first input magnetic field vector.
bis_field_2	float	The magnitude of the bisected field line corresponding to the second input magnetic field vector.

Source code in src/rxn_location/rx_model_funcs.py

def get_bis(b_vec_1, b_vec_2):
    r"""
    Calculate the bisector field components of two magnetic field lines.

    Parameters
    ----------
    b_vec_1 : array-like of shape (3,)
        First input magnetic field vector.
    b_vec_2 : array-like of shape (3,)
        Second input magnetic field vector.

    Returns
    -------
    bis_field_1 : float
        The magnitude of the bisected field line corresponding to the first input magnetic field vector.
    bis_field_2 : float
        The magnitude of the bisected field line corresponding to the second input magnetic field vector.
    """
    b_vec_1 = np.array(b_vec_1)
    b_vec_2 = np.array(b_vec_2)

    mag_vec_1 = np.linalg.norm(b_vec_1)
    mag_vec_2 = np.linalg.norm(b_vec_2)

    unit_vec_1 = b_vec_1 / mag_vec_1
    unit_vec_2 = b_vec_2 / mag_vec_2

    unit_vec_bisect = (unit_vec_1 + unit_vec_2) / np.linalg.norm(unit_vec_1 + unit_vec_2)

    # Cross product of the two vectors with the bisector
    b_vec_1_cross = np.cross(b_vec_1, unit_vec_bisect)
    b_vec_2_cross = np.cross(b_vec_2, unit_vec_bisect)

    # Magnitude of the cross product vectors (should be the same for both for a symmetric
    # reconnection). These magnitude are the reconnecting components of the magnetic field.
    bis_field_1 = np.linalg.norm(b_vec_1_cross)
    bis_field_2 = np.linalg.norm(b_vec_2_cross)
    return bis_field_1, bis_field_2

get_ca(b_vec, angle_unit='radians')

Calculate the cone angle of a magnetic field vector.

Parameters:

Name	Type	Description	Default
b_vec	array-like of shape (3,)	Input magnetic field vector.	required
angle_unit	str	Preferred unit of angle returned by the code. Default is "radians".	'radians'

Raises:

Type	Description
KeyError	If angle_unit is not set to "radians" or "degrees".

Returns:

Name	Type	Description
angle	float	Returns the arctangent of the y- and z-components of the magnetic field vector.

Source code in src/rxn_location/rx_model_funcs.py

def get_ca(b_vec, angle_unit="radians"):
    r"""
    Calculate the cone angle of a magnetic field vector.

    Parameters
    ----------
    b_vec : array-like of shape (3,)
        Input magnetic field vector.
    angle_unit : str, optional
        Preferred unit of angle returned by the code. Default is "radians".

    Raises
    ------
    KeyError
        If `angle_unit` is not set to "radians" or "degrees".

    Returns
    -------
    angle : float
        Returns the arctangent of the y- and z-components of the magnetic field vector.
    """
    angle = np.arctan(b_vec[1] / b_vec[2])

    if angle_unit == "radians":
        return angle
    elif angle_unit == "degrees":
        return angle * 180 / np.pi
    else:
        raise KeyError("angle_unit must be radians or degrees")

get_rxben(b_vec_1, b_vec_2)

Calculate the reconnection energy between two magnetic field lines.

It has the following mathematical expression:

.. math:: rexben = 0.5 (|\vec{B_1}| + |\vec{B_2}|) (1 - \hat{B_1} \cdot \hat{B_2})

Parameters:

Name	Type	Description	Default
b_vec_1	array-like of shape (3,)	First input magnetic field vector.	required
b_vec_2	array-like of shape (3,)	Second input magnetic field vector.	required

Returns:

Name	Type	Description
rxben	float	Reconnection field energy density in nPa.

Source code in src/rxn_location/rx_model_funcs.py

def get_rxben(b_vec_1, b_vec_2):
    r"""
    Calculate the reconnection energy between two magnetic field lines.

    It has the following mathematical expression:

    .. math:: rexben = 0.5 (|\vec{B_1}| + |\vec{B_2}|) (1 - \hat{B_1} \cdot \hat{B_2})

    Parameters
    ----------
    b_vec_1 : array-like of shape (3,)
        First input magnetic field vector.
    b_vec_2 : array-like of shape (3,)
        Second input magnetic field vector.

    Returns
    -------
    rxben : float
        Reconnection field energy density in nPa.
    """

    # alpha = - 14.87 * np.pi / 180  # radians (From Hesse2013)
    b_vec_1 = np.array(b_vec_1)
    b_vec_2 = np.array(b_vec_2)
    mag_vec_1 = np.linalg.norm(b_vec_1)
    mag_vec_2 = np.linalg.norm(b_vec_2)

    unit_vec_1 = b_vec_1 / mag_vec_1
    unit_vec_2 = b_vec_2 / mag_vec_2

    # The bisector vector of thwo input vectors
    unit_vec_bisec = (unit_vec_1 + unit_vec_2) / np.linalg.norm(unit_vec_1 + unit_vec_2)

    # Cross product of the two input vectors with the bisector vector to get the reconnection
    # component of the magnetic field.
    rx_b_1 = np.cross(b_vec_1, unit_vec_bisec)
    rx_b_2 = np.cross(b_vec_2, unit_vec_bisec)

    rx_b_mag_1 = np.linalg.norm(rx_b_1)
    rx_b_mag_2 = np.linalg.norm(rx_b_2)

    # Reconnection energy (from Hesse2013)
    rx_en = rx_b_mag_1**2 * rx_b_mag_2**2

    return rx_en

get_shear(b_vec_1, b_vec_2, angle_unit='radians')

Calculate the magnetic shear angle between two magnetic field vectors.

Parameters:

Name	Type	Description	Default
b_vec_1	array-like of shape (3,)	First input magnetic field vector.	required
b_vec_2	array-like of shape (3,)	Second input magnetic field vector.	required
angle_unit	str	Preferred unit of angle returned by the code. Default is "radians".	'radians'

Raises:

Type	Description
KeyError	If angle_unit is not set to "radians" or "degrees".

Returns:

Name	Type	Description
angle	float	Angle between the two vectors, in radians by default.

Source code in src/rxn_location/rx_model_funcs.py

def get_shear(b_vec_1, b_vec_2, angle_unit="radians"):
    r"""
    Calculate the magnetic shear angle between two magnetic field vectors.

    Parameters
    ----------
    b_vec_1 : array-like of shape (3,)
        First input magnetic field vector.
    b_vec_2 : array-like of shape (3,)
        Second input magnetic field vector.
    angle_unit : str, optional
        Preferred unit of angle returned by the code. Default is "radians".

    Raises
    ------
    KeyError
        If `angle_unit` is not set to "radians" or "degrees".

    Returns
    -------
    angle : float
        Angle between the two vectors, in radians by default.
    """
    unit_vec_1 = b_vec_1 / np.linalg.norm(b_vec_1)
    unit_vec_2 = b_vec_2 / np.linalg.norm(b_vec_2)
    angle = np.arccos(np.dot(unit_vec_1, unit_vec_2))

    if angle_unit == "radians":
        return angle
    elif angle_unit == "degrees":
        return angle * 180 / np.pi
    else:
        raise KeyError("angle_unit must be radians or degrees")

get_sw_params(probe=None, omni_level='hro', time_clip=True, trange=None, mms_probe_num=None, latest_version=False, verbose=False)

Retrieve and process solar wind parameters and spacecraft location from the OMNI and MMS databases.

Parameters:

Name	Type	Description	Default
probe	str	The probe to use. Default is "None".	None
omni_level	str	The omni data level to use. Options are "hro" and "hro2". Default is "hro".	'hro'
time_clip	bool	If True, the data will be clipped to the time range specified by trange. Default is True.	True
trange	list or an array of length 2	The time range to use. Should in the format [start, end], where start and end times should be a string in the format "YYYY-MM-DD HH:MM:SS".	None
mms_probe_num	str	The MMS probe to use. Options are "1", "2", "3" and "4". Default is None.	None
verbose	bool	If True, print out a few messages and the solar wind parameters. Default is False.	False

Raises:

Type	Description
ValueError: If the probe is not one of the options.
ValueError: If the trange is not in the correct format.

Returns:

Name	Type	Description
sw_params	dict	The solar wind parameters.

Source code in src/rxn_location/rx_model_funcs.py

def get_sw_params(
    probe=None,
    omni_level="hro",
    time_clip=True,
    trange=None,
    mms_probe_num=None,
    latest_version=False,
    verbose=False,
):
    r"""
    Retrieve and process solar wind parameters and spacecraft location from the OMNI and MMS databases.

    Parameters
    ----------
    probe : str
        The probe to use. Default is "None".
    omni_level : str
        The omni data level to use. Options are "hro" and "hro2". Default is "hro".
    time_clip : bool
        If True, the data will be clipped to the time range specified by trange. Default is True.
    trange : list or an array of length 2
        The time range to use. Should in the format [start, end], where start and end times should
        be a string in the format "YYYY-MM-DD HH:MM:SS".
    mms_probe_num : str
        The MMS probe to use. Options are "1", "2", "3" and "4". Default is None.
    verbose : bool
        If True, print out a few messages and the solar wind parameters. Default is False.

    Raises
    ------
    ValueError: If the probe is not one of the options.
    ValueError: If the trange is not in the correct format.

    Returns
    -------
    sw_params : dict
        The solar wind parameters.
    """

    if trange is None:
        raise ValueError(
            "trange must be specified as a list of start and end times in the format"
            "'YYYY-MM-DD HH:MM:SS'."
        )

    # Check if trange is either a list or an array of length 2
    if not isinstance(trange, (list, np.ndarray)) or len(trange) != 2:
        raise ValueError(
            "trange must be specified as a list or array of length 2 in the format"
            + "'YYYY-MM-DD HH:MM:SS'"
        )

    # Download the OMNI data (default level of "hro_1min") for the specified timerange.
    omni_varnames = ["BX_GSE", "BY_GSM", "BZ_GSM", "proton_density", "Vx", "Vy", "Vz", "SYM_H", "T"]
    omni_vars = spd.omni.data(
        trange=trange, varnames=omni_varnames, level=omni_level, time_clip=time_clip
    )

    omni_time = spd.get_data(omni_vars[0])[0]
    # vprint(2, f"omni_time: {omni_time}")
    omni_bx_gse = spd.get_data(omni_vars[0])[1]
    omni_by_gsm = spd.get_data(omni_vars[1])[1]
    omni_bz_gsm = spd.get_data(omni_vars[2])[1]
    omni_np = spd.get_data(omni_vars[3])[1]
    omni_vx = spd.get_data(omni_vars[4])[1]
    omni_vy = spd.get_data(omni_vars[5])[1]
    omni_vz = spd.get_data(omni_vars[6])[1]
    omni_sym_h = spd.get_data(omni_vars[7])[1]
    omni_t_p = spd.get_data(omni_vars[8])[1]

    # Convert omni_time to datetime objects from unix time
    omni_time_datetime = [datetime.datetime.fromtimestamp(t, datetime.timezone.utc) for t in omni_time]
    # Get trange in datetime format
    omni_trange_time_object = [pd.to_datetime(trange[0]), pd.to_datetime(trange[1])]
    # Add utc as timezone to omni_trange_time_object
    omni_trange_time_object = [
        t.replace(tzinfo=datetime.timezone.utc) for t in omni_trange_time_object
    ]

    # Get mms postion in GSM coordinates for the specified time range

    if mms_probe_num is not None:
        # Define mms time as the center of the time range
        mms_time_t0 = pd.to_datetime(trange[0])
        mms_time_dt = pd.to_datetime(trange[1]) - mms_time_t0
        mms_time = mms_time_t0 + mms_time_dt / 2

        # Define mms time range as +/- 10 seconds from the mms time
        mms_trange = [mms_time - pd.Timedelta("5 second"), mms_time + pd.Timedelta("5 second")]
        mms_mec_trange = [
            mms_time - pd.Timedelta("30 second"),
            mms_time + pd.Timedelta("30 second"),
        ]
        # Convert mms time range to string
        mms_trange = [
            mms_trange[0].strftime("%Y-%m-%d %H:%M:%S"),
            mms_trange[1].strftime("%Y-%m-%d %H:%M:%S"),
        ]
        mms_mec_trange = [
            mms_mec_trange[0].strftime("%Y-%m-%d %H:%M:%S"),
            mms_mec_trange[1].strftime("%Y-%m-%d %H:%M:%S"),
        ]
        # Convert mms time range to datetime and localize to UTC
        mms_trange_time_object = [
            pd.to_datetime(mms_trange[0]).replace(tzinfo=datetime.timezone.utc),
            pd.to_datetime(mms_trange[1]).replace(tzinfo=datetime.timezone.utc),
        ]
        mms_mec_trange_time_object = [
            pd.to_datetime(mms_mec_trange[0]).replace(tzinfo=datetime.timezone.utc),
            pd.to_datetime(mms_mec_trange[1]).replace(tzinfo=datetime.timezone.utc),
        ]

        mms_mec_varnames = [f"mms{mms_probe_num}_mec_r_gsm"]
        _ = mms.mec(
            trange=mms_mec_trange,
            varnames=mms_mec_varnames,
            probe=mms_probe_num,
            data_rate="srvy",
            level="l2",
            time_clip=time_clip,
            latest_version=latest_version,
        )
        mms_mec_time = spd.get_data(f"mms{mms_probe_num}_mec_r_gsm")[0]
        # Convert mms fgm time to datetime
        mms_mec_time = np.array([datetime.datetime.fromtimestamp(t, datetime.timezone.utc) for t in mms_mec_time])

        # Position of MMS in GSM coordinates in earth radii (r_e) units
        r_e = 6378.137  # Earth radius in km
        mms_sc_pos = spd.get_data(f"mms{mms_probe_num}_mec_r_gsm")[1:3][0] / r_e

        # TODO: Find out why adding "mms_fgm_varnames" as a variable causes the code to give out no
        # data.
        mms_fgm_varnames = [f"mms{mms_probe_num}_fgm_b_gsm_srvy_l2"]
        _ = mms.fgm(
            trange=mms_trange,
            probe=mms_probe_num,
            time_clip=time_clip,
            latest_version=latest_version,
        )
        # mms_fgm_time = spd.get_data(mms_fgm_varnames[0])[0]

        mms_fgm_b_gsm = spd.get_data(f"mms{mms_probe_num}_fgm_b_gsm_srvy_l2_bvec")[1:4][0]
        mms_fgm_time = spd.get_data(f"mms{mms_probe_num}_fgm_b_gsm_srvy_l2_bvec")[0]
        # Convert mms fgm time to datetime
        mms_fgm_time = np.array([datetime.datetime.fromtimestamp(t, datetime.timezone.utc) for t in mms_fgm_time])

        try:
            data_rate = "fast"
            mms_fpi_varnames = [f"mms{mms_probe_num}_dis_bulkv_gse_{data_rate}"]
            _ = mms.fpi(
                trange=mms_trange,
                probe=mms_probe_num,
                data_rate=data_rate,
                level="l2",
                datatype="dis-moms",
                varnames=mms_fpi_varnames,
                time_clip=time_clip,
                latest_version=latest_version,
            )
            _ = spd.cotrans(
                name_in=f"mms{mms_probe_num}_dis_bulkv_gse_{data_rate}",
                name_out=f"mms{mms_probe_num}_dis_bulkv_gsm_{data_rate}",
                coord_in="gse",
                coord_out="gsm",
            )
            mms_fpi_time = spd.get_data(f"mms{mms_probe_num}_dis_bulkv_gsm_{data_rate}")[0]
            # Convert mms_fpi_time to datetime from unix time
            mms_fpi_time = np.array([datetime.datetime.fromtimestamp(x, datetime.timezone.utc) for x in mms_fpi_time])
            mms_fpi_bulkv_gsm = spd.get_data(f"mms{mms_probe_num}_dis_bulkv_gsm_{data_rate}")[1:4][
                0
            ]
            if verbose:
                vprint(2, 
                    f"\n \033[1;31m {data_rate} mode data found for MMS{mms_probe_num} \033[0m \n"
                )
        except:
            # Get the data from the FPI
            data_rate = "brst"
            mms_fpi_varnames = [f"mms{mms_probe_num}_dis_bulkv_gse_{data_rate}"]
            _ = mms.fpi(
                trange=mms_trange,
                probe=mms_probe_num,
                data_rate=data_rate,
                level="l2",
                datatype="dis-moms",
                varnames=mms_fpi_varnames,
                time_clip=time_clip,
                latest_version=latest_version,
            )
            _ = spd.cotrans(
                name_in=f"mms{mms_probe_num}_dis_bulkv_gse_{data_rate}",
                name_out=f"mms{mms_probe_num}_dis_bulkv_gsm_{data_rate}",
                coord_in="gse",
                coord_out="gsm",
            )
            mms_fpi_time = spd.get_data(f"mms{mms_probe_num}_dis_bulkv_gsm_{data_rate}")[0]
            # Convert mms_fpi_time to datetime from unix time
            mms_fpi_time = [datetime.datetime.fromtimestamp(x, datetime.timezone.utc) for x in mms_fpi_time]
            mms_fpi_bulkv_gsm = spd.get_data(f"mms{mms_probe_num}_dis_bulkv_gsm_{data_rate}")[1:4][
                0
            ]
            if verbose:
                vprint(2, 
                    f"\n \033[1;32m {data_rate} mode data found for MMS{mms_probe_num} \033[0m \n"
                )
    else:
        mms_time = None
        mms_sc_pos = None
        mms_fgm_b_gsm = None
        mms_fpi_bulkv_gsm = None
        pass

    # Create the dataframe for OMNI data using omni_time_datetime as the index
    omni_df = pd.DataFrame(
        {
            "time": omni_time,
            "bx_gsm": omni_bx_gse,
            "by_gsm": omni_by_gsm,
            "bz_gsm": omni_bz_gsm,
            "vx": omni_vx,
            "vy": omni_vy,
            "vz": omni_vz,
            "np": omni_np,
            "sym_h": omni_sym_h,
            "t_p": omni_t_p,
        },
        index=omni_time_datetime,
    )

    # Add UTC as time zone to the index of omni_df, handling both naive and aware datetimes safely
    omni_df.index = omni_df.index.tz_localize(None).tz_localize("UTC")

    # Get the mean values of the parameters from OMNI data for the time range betwwen
    time_imf = np.nanmean(
        omni_df["time"].loc[omni_trange_time_object[0] : omni_trange_time_object[1]]
    )
    b_imf_x = np.nanmean(
        omni_df["bx_gsm"].loc[omni_trange_time_object[0] : omni_trange_time_object[1]]
    )
    b_imf_y = np.nanmean(
        omni_df["by_gsm"].loc[omni_trange_time_object[0] : omni_trange_time_object[1]]
    )
    b_imf_z = np.nanmean(
        omni_df["bz_gsm"].loc[omni_trange_time_object[0] : omni_trange_time_object[1]]
    )
    vx_imf = np.nanmean(omni_df["vx"].loc[omni_trange_time_object[0] : omni_trange_time_object[1]])
    vy_imf = np.nanmean(omni_df["vy"].loc[omni_trange_time_object[0] : omni_trange_time_object[1]])
    vz_imf = np.nanmean(omni_df["vz"].loc[omni_trange_time_object[0] : omni_trange_time_object[1]])
    np_imf = np.nanmean(omni_df["np"].loc[omni_trange_time_object[0] : omni_trange_time_object[1]])
    sym_h_imf = np.nanmean(
        omni_df["sym_h"].loc[omni_trange_time_object[0] : omni_trange_time_object[1]]
    )
    tp_imf = np.nanmean(omni_df["t_p"].loc[omni_trange_time_object[0] : omni_trange_time_object[1]])

    if b_imf_z > 15 or b_imf_z < -18:
        warnings.warn(
            f"The given parameters produced the z-component of IMF field (b_imf_z) {b_imf_z} nT,"
            f"which is out of range in which model is valid (-18 nT < b_imf_z < 15 nT)"
        )

    time_imf_hrf = datetime.datetime.fromtimestamp(time_imf, datetime.timezone.utc)

    v_imf = [vx_imf, vy_imf, vz_imf]
    b_imf = [b_imf_x, b_imf_y, b_imf_z]

    imf_clock_angle = np.arctan2(b_imf[1], b_imf[2]) * 180 / np.pi
    if imf_clock_angle < 0:
        imf_clock_angle += 360
    if mms_probe_num is not None:
        mean_mms_sc_pos = np.round(np.nanmean(mms_sc_pos, axis=0), decimals=2)
    else:
        mean_mms_sc_pos = None

    vprint(1, "IMF parameters found:", color="green")
    if verbose:
        vprint(2, 
            tabulate(
                [
                    ["Time of observation (UTC)", f"{time_imf_hrf}"],
                    [
                        "IMF Magnetic field [GSM] (nT)",
                        f"[{b_imf[0]:.2f}, {b_imf[1]:.2f}, {b_imf[2]:.2f}]",
                    ],
                    ["IMF Proton density (1/cm^-3)", f"{np_imf:.2f}"],
                    [
                        "IMF Plasma velocity (km/sec)",
                        f"[{v_imf[0]:.2f}, {v_imf[1]:.2f}, {v_imf[2]:.2f}]",
                    ],
                    ["IMF clock angle (degrees)", f"{imf_clock_angle:.2f}"],
                    ["IMF Sym H", f"{sym_h_imf:.2f}"],
                    [
                        "MMS position (GSM) (R_E)",
                        f"[{mean_mms_sc_pos[0]:.2f}, {mean_mms_sc_pos[1]:.2f}, "
                        f"{mean_mms_sc_pos[2]:.2f}]",
                    ],
                ],
                headers=["Parameter", "Value"],
                tablefmt="fancy_grid",
                floatfmt=".2f",
                numalign="center",
            )
        )

    # Check if the values are finite, if not then assign a default value to each of them
    if ~(np.isfinite(np_imf)):
        np_imf = 5
    if ~(np.isfinite(vx_imf)):
        vx_imf = -500
    if ~(np.isfinite(vy_imf)):
        vy_imf = 0
    if ~(np.isfinite(vz_imf)):
        vz_imf = 0
    if ~(np.isfinite(sym_h_imf)):
        sym_h_imf = -1

    m_proton = 1.672e-27  # Mass of proton in SI unit

    rho = np_imf * m_proton * 1.15  # NOTE to self: Unit is fine, do not worry about it

    #  Solar wind ram pressure in nPa, including roughly 4% Helium++ contribution
    p_dyn = 1.6726e-6 * 1.15 * np_imf * (vx_imf**2 + vy_imf**2 + vz_imf**2)

    if p_dyn > 8.5 or p_dyn < 0.5:
        warnings.warn(
            f"The given parameters produced a dynamic pressure of {p_dyn} nPa which is out of"
            f" range in which model is valid (0.5 nPa < p_dyn < 8.5 nPa)",
        )
    param = [p_dyn, sym_h_imf, b_imf_y, b_imf_z, 0, 0, 0, 0, 0, 0]

    # Compute the dipole tilt angle
    ps = gp.recalc(time_imf)

    # Create a dataframe for mms fgm data
    df_fgm = pd.DataFrame(data=mms_fgm_b_gsm, columns=["Bx", "By", "Bz"], index=mms_fgm_time)

    # Get the mean of the fgm data for time range between mms_trange_time_object[0] and
    bx_mean = np.nanmean(df_fgm["Bx"].loc[mms_trange_time_object[0] : mms_trange_time_object[1]])
    by_mean = np.nanmean(df_fgm["By"].loc[mms_trange_time_object[0] : mms_trange_time_object[1]])
    bz_mean = np.nanmean(df_fgm["Bz"].loc[mms_trange_time_object[0] : mms_trange_time_object[1]])

    # Create a dataframe for mms fpi data
    df_fpi = pd.DataFrame(data=mms_fpi_bulkv_gsm, columns=["Vx", "Vy", "Vz"], index=mms_fpi_time)

    # Get the mean of the fpi data for time range between mms_trange_time_object[0] and
    vx_mean = np.nanmean(df_fpi["Vx"].loc[mms_trange_time_object[0] : mms_trange_time_object[1]])
    vy_mean = np.nanmean(df_fpi["Vy"].loc[mms_trange_time_object[0] : mms_trange_time_object[1]])
    vz_mean = np.nanmean(df_fpi["Vz"].loc[mms_trange_time_object[0] : mms_trange_time_object[1]])

    # Create a dataframe for mms sc position data
    df_mec = pd.DataFrame(data=mms_sc_pos, columns=["X", "Y", "Z"], index=mms_mec_time)

    # Get the mean of the sc position data for time range between mms_trange_time_object[0] and
    x_mean = np.nanmean(
        df_mec["X"].loc[mms_mec_trange_time_object[0] : mms_mec_trange_time_object[1]]
    )
    y_mean = np.nanmean(
        df_mec["Y"].loc[mms_mec_trange_time_object[0] : mms_mec_trange_time_object[1]]
    )
    z_mean = np.nanmean(
        df_mec["Z"].loc[mms_mec_trange_time_object[0] : mms_mec_trange_time_object[1]]
    )

    # Make a dictionary of all the solar wind parameters
    sw_dict = {}
    sw_dict["time"] = time_imf
    sw_dict["b_imf"] = b_imf
    sw_dict["v_imf"] = v_imf
    sw_dict["np"] = np_imf
    sw_dict["rho"] = rho
    sw_dict["ps"] = ps
    sw_dict["p_dyn"] = p_dyn
    sw_dict["sym_h"] = sym_h_imf
    sw_dict["t_p"] = tp_imf
    sw_dict["imf_clock_angle"] = imf_clock_angle
    sw_dict["param"] = param
    sw_dict["mms_time"] = mms_time
    sw_dict["mms_sc_pos"] = [x_mean, y_mean, z_mean]
    sw_dict["mms_b_gsm"] = [bx_mean, by_mean, bz_mean]
    sw_dict["mms_v_gsm"] = [vx_mean, vy_mean, vz_mean]

    return sw_dict

get_vcs(b_vec_1, b_vec_2, n_1, n_2)

Calculate the Cassak-Shay asymmetric reconnection exhaust velocity.

Parameters:

Name	Type	Description	Default
b_vec_1	array-like of shape (3,)	First input magnetic field vector.	required
b_vec_2	array-like of shape (3,)	Second input magnetic field vector.	required
n_1	float	Plasma number density of the first component.	required
n_2	float	Plasma number density of the second component.	required

Returns:

Name	Type	Description
vcs	float	The exhaust velocity in km/s.

Source code in src/rxn_location/rx_model_funcs.py

def get_vcs(b_vec_1, b_vec_2, n_1, n_2):
    r"""
    Calculate the Cassak-Shay asymmetric reconnection exhaust velocity.

    Parameters
    ----------
    b_vec_1 : array-like of shape (3,)
        First input magnetic field vector.
    b_vec_2 : array-like of shape (3,)
        Second input magnetic field vector.
    n_1 : float
        Plasma number density of the first component.
    n_2 : float
        Plasma number density of the second component.

    Returns
    -------
    vcs : float
        The exhaust velocity in km/s.
    """
    va_p1 = 21.812  # conv. nT, m_P/cm ^ 3 product to km/s cassak-shay

    mag_vec_1 = np.linalg.norm(b_vec_1)
    mag_vec_2 = np.linalg.norm(b_vec_2)

    unit_vec_1 = b_vec_1 / mag_vec_1
    unit_vec_2 = b_vec_2 / mag_vec_2

    # The bisector vector of thwo input vectors
    unit_vec_bisec = (unit_vec_1 + unit_vec_2) / np.linalg.norm(unit_vec_1 + unit_vec_2)

    # Cross product of the two input vectors with the bisector vector to get the reconnection
    # component of the magnetic field.
    rx_b_1 = np.cross(b_vec_1, unit_vec_bisec)
    rx_b_2 = np.cross(b_vec_2, unit_vec_bisec)

    # Magnitude of the reconnection component of the magnetic fields
    rx_mag_1 = np.linalg.norm(rx_b_1)
    rx_mag_2 = np.linalg.norm(rx_b_2)

    vcs = va_p1 * np.sqrt(
        rx_mag_1 * rx_mag_2 * (rx_mag_1 + rx_mag_2) / (rx_mag_1 * n_2 + rx_mag_2 * n_1)
    )

    return vcs

line_fnc(r0=np.array([0, 0, 0]), b_msh=np.array([0, 0, 0]), r=0)

Function to compute the equation of a line along the direction of the magnetic field.

Parameters:

Name	Type	Description	Default
r0	array	Starting position of the line	array([0, 0, 0])
b_msh	array	Magnetic field from magnetosheath at the magnetopause	array([0, 0, 0])
r	float	Position of the point where the line is computed	0

Source code in src/rxn_location/rx_model_funcs.py

def line_fnc(
    r0=np.array([0, 0, 0]),
    b_msh=np.array([0, 0, 0]),
    r=0,
):
    """
    Function to compute the equation of a line along the direction of the magnetic field.

    Parameters
    ----------
    r0 : array
        Starting position of the line
    b_msh : array, optional
        Magnetic field from magnetosheath at the magnetopause
    r : float, optional
        Position of the point where the line is computed
    """
    # Compute the direction of the magnetic field line
    b_msh_norm = np.linalg.norm(b_msh)
    b_msh_dir = b_msh / b_msh_norm

    # return the line function
    return r0 + r * b_msh_dir

line_fnc_der(x, y)

Function to give a line interpolation function

Parameters:

Name	Type	Description	Default
x	array	x coordinates of the points	required
y	array	y coordinates of the points	required

Returns:

Name	Type	Description
line_intrp	function	Function to interpolate the line

Source code in src/rxn_location/rx_model_funcs.py

def line_fnc_der(x, y):
    """
    Function to give a line interpolation function

    Parameters
    ----------
    x : array
        x coordinates of the points
    y : array
        y coordinates of the points

    Returns
    -------
    line_intrp : function
        Function to interpolate the line
    """

    nans, x_nan = nan_helper(x)
    x[nans] = np.interp(x_nan(nans), x_nan(~nans), x[~nans])

    nans, y_nan = nan_helper(y)
    y[nans] = np.interp(y_nan(nans), y_nan(~nans), y[~nans])

    line_intrp = sp.interpolate.CubicSpline(x, y)
    return line_intrp

model_run(*args)

Calculate the magnetic field and reconnection parameters at a specific grid point.

This function is designed to be mapped concurrently across a grid of points on the magnetopause. It returns the magnetic field values and associated reconnection metrics based on the given model.

Parameters:

Name	Type	Description	Default
args	tuple	A single tuple containing all arguments packed together for multiprocessing. The elements are: (j, k, y_max, z_max, dr, m_p, ro, alpha, rmp, sw_params, model_type).	()

Returns:

Type	Description
tuple	A tuple containing the following elements: (j, k, bx, by, bz, shear, rx_en, va_cs, bisec_msp, bisec_msh, x_shu, y_shu, z_shu, b_msx, b_msy, b_msz)

Source code in src/rxn_location/rx_model_funcs.py

def model_run(*args):
    r"""
    Calculate the magnetic field and reconnection parameters at a specific grid point.

    This function is designed to be mapped concurrently across a grid of points on the magnetopause.
    It returns the magnetic field values and associated reconnection metrics based on the given model.

    Parameters
    ----------
    args : tuple
        A single tuple containing all arguments packed together for multiprocessing.
        The elements are: (j, k, y_max, z_max, dr, m_p, ro, alpha, rmp, sw_params, model_type).

    Returns
    -------
    tuple
        A tuple containing the following elements:
        (j, k, bx, by, bz, shear, rx_en, va_cs, bisec_msp, bisec_msh, x_shu, y_shu, z_shu, b_msx, b_msy, b_msz)
    """
    j = args[0][0]
    k = args[0][1]
    y_max = args[0][2]
    z_max = args[0][3]
    dr = args[0][4]
    m_p = args[0][5]
    ro = args[0][6]
    alpha = args[0][7]
    rmp = args[0][8]
    sw_params = args[0][9]
    model_type = args[0][-1]

    y0 = int(j * dr) - y_max
    z0 = int(k * dr) - z_max
    rp = np.sqrt(y0**2 + z0**2)  # Projection of r into yz-plane

    d_theta = np.pi / 100

    for index in range(0, 100):

        theta = index * d_theta
        r = ro * (2 / (1 + np.cos(theta))) ** alpha
        zp = r * np.sin(theta)  # not really in z direction, but a distance in yz plane
        x0 = r * np.cos(theta)

        if y0 == 0:
            signy = 1.0
        else:
            signy = np.sign(y0)

        if z0 == 0:
            signz = 1.0
        else:
            signz = np.sign(z0)

        if rp <= zp:

            x_shu = (r - m_p) * np.cos(theta)
            phi = np.arctan2(z0, y0)

            if abs(y0) == 0 or abs(z0) == 0:
                if abs(y0) == 0:
                    y_shu = 0
                    z_shu = (r - m_p) * np.sin(theta)
                elif abs(z0) == 0:
                    z_shu = 0
                    y_shu = (r - m_p) * np.sin(theta)
            else:
                z_shu = np.sqrt((rp - 1.0) ** 2 / (1 + np.tan(phi) ** (-2)))
                y_shu = z_shu / np.tan(phi)

            m_proton = 1.672e-27  # Mass of proton in SI unit
            n_sh = sw_params["rho"] * (1.509 * np.exp(x_shu / rmp) + 0.1285) / m_proton

            y_shu = abs(y_shu) * signy
            z_shu = abs(z_shu) * signz

            # Cooling JGR 2001 Model, equation 9 to 12
            # the distance from the focus to the magnetopause surface
            A = 2
            ll = 3 * rmp / 2 - x0
            b_msx = -A * (
                -sw_params["b_imf"][0] * (1 - rmp / (2 * ll))
                + sw_params["b_imf"][1] * (y0 / ll)
                + sw_params["b_imf"][2] * (z0 / ll)
            )
            b_msy = A * (
                -sw_params["b_imf"][0] * (y0 / (2 * ll))
                + sw_params["b_imf"][1] * (2 - y0**2 / (ll * rmp))
                - sw_params["b_imf"][2] * (y0 * z0 / (ll * rmp))
            )
            b_msz = A * (
                -sw_params["b_imf"][0] * (z0 / (2 * ll))
                - sw_params["b_imf"][1] * (y0 * z0 / (ll * rmp))
                + sw_params["b_imf"][2] * (2 - z0**2 / (ll * rmp))
            )
            try:
                if model_type == "t96":
                    bx_ext, by_ext, bz_ext = gp.t96.t96(
                        sw_params["param"], sw_params["ps"], x_shu, y_shu, z_shu
                    )
                elif model_type == "t01":
                    bx_ext, by_ext, bz_ext = gp.t01.t01(
                        sw_params["param"], sw_params["ps"], x_shu, y_shu, z_shu
                    )
            except Exception:
                vprint(2, f"Skipped for {x_shu, y_shu, z_shu}")
                pass

            bx_igrf, by_igrf, bz_igrf = gp.igrf_gsm(x_shu, y_shu, z_shu)

            # vprint(2, j, k, bx_ext, bx_igrf)
            bx = bx_ext + bx_igrf
            by = by_ext + by_igrf
            bz = bz_ext + bz_igrf

            shear = get_shear([bx, by, bz], [b_msx, b_msy, b_msz], angle_unit="degrees")

            rx_en = get_rxben([bx, by, bz], [b_msx, b_msy, b_msz])
            va_cs = get_vcs([bx, by, bz], [b_msx, b_msy, b_msz], n_sh, 0.1)
            bisec_msp, bisec_msh = get_bis([bx, by, bz], [b_msx, b_msy, b_msz])
            break

    return (
        j,
        k,
        bx,
        by,
        bz,
        shear,
        rx_en,
        va_cs,
        bisec_msp,
        bisec_msh,
        x_shu,
        y_shu,
        z_shu,
        b_msx,
        b_msy,
        b_msz,
    )

nan_helper(y)

Helper to handle indices and logical indices of NaNs.

Input: - y, 1d numpy array with possible NaNs Output: - nans, logical indices of NaNs - index, a function, with signature indices= index(logical_indices), to convert logical indices of NaNs to "equivalent" indices Example: >>> # linear interpolation of NaNs >>> nans, x= nan_helper(y) >>> y[nans]= np.interp(x(nans), x(~nans), y[~nans])

Source code in src/rxn_location/rx_model_funcs.py

def nan_helper(y):
    """Helper to handle indices and logical indices of NaNs.

    Input:
        - y, 1d numpy array with possible NaNs
    Output:
        - nans, logical indices of NaNs
        - index, a function, with signature indices= index(logical_indices),
          to convert logical indices of NaNs to "equivalent" indices
    Example:
        >>> # linear interpolation of NaNs
        >>> nans, x= nan_helper(y)
        >>> y[nans]= np.interp(x(nans), x(~nans), y[~nans])
    """

    return np.isnan(y), lambda z: z.nonzero()[0]

ridge_finder_multiple(image=[None, None, None, None], convolution_order=[1, 1, 1, 1], t_range=['2016-12-24 15:08:00', '2016-12-24 15:12:00'], dt=5, b_imf=[-5, 0, 0], b_msh=[-5, 0, 0], v_msh=[-200, 50, 50], xrange=[-15.1, 15], yrange=[-15.1, 15], mms_probe_num='1', mms_sc_pos=[0, 0], dr=0.5, dipole_tilt_angle=None, p_dyn=None, imf_clock_angle=None, np_imf=None, v_imf=[None, None, None], sym_h=None, sigma=[2.2, 2.2, 2.2, 2.2], mode='nearest', alpha=1.0, vmin=[None, None, None, None], vmax=[None, None, None, None], cmap_list=['viridis', 'viridis', 'viridis', 'viridis'], draw_patch=[True, True, True, True], draw_ridge=[True, True, True, True], save_fig=True, fig_name='new', fig_format='png', c_label=[None, None, None, None], wspace=0.1, hspace=0.1, fig_size=(10, 10), box_style=None, title_y_pos=0.95, interpolation='nearest', tsy_model='t96', dark_mode=True, rc_file_name='rc_file.csv', rc_folder='data', save_rc_file=False, fig_version='v001', df_jet_reversal=None, **kwargs)

Finds ridges in an image and plot the points with maximum ridge value on the given image.

Parameters:

Name	Type	Description	Default
image	list of numpy arrays	List of images to be plotted.	[None, None, None, None]
convolution_order	list of ints	List of the order of the convolution to be used while smoothing each image. Values must be non-negative integers. Default is [1, 1, 1, 1].	[1, 1, 1, 1]
t_range	list of str	The time range to find the ridge in. Default is ["2016-12-24 15:08:00", "2016-12-24 15:12:00"].	['2016-12-24 15:08:00', '2016-12-24 15:12:00']
dt	float	The time differential, in minutes, for observation if "t_range" has only one element. Default is 5 minutes.	5
b_imf	list of floats	The IMF vector. Default is [-5, 0, 0].	[-5, 0, 0]
b_msh	list of floats	The MSH vector. Default is [-5, 0, 0].	[-5, 0, 0]
v_msh	list of floats	The MSH velocity vector. Default is [-200, 50, 50].	[-200, 50, 50]
xrange	list of floats	The range of x-values for image. Default is [-15.1, 15].	[-15.1, 15]
yrange	list of floats	The range of y-values for image. Default is [-15.1, 15].	[-15.1, 15]
mms_probe_num	str	The probe number of the MMS spacecraft. Default is 1.	'1'
mms_sc_pos	list of floats	The position of the spacecraft in the image. Default is [0, 0].	[0, 0]
dr	float	The step size for the grid. Default is 0.5.	0.5
dipole_tilt_angle	float	The dipole tilt angle. Default is None.	None
p_dyn	float	The dynamic pressure. Default is None.	None
imf_clock_angle	float	The IMF clock angle. Default is None.	None
sigma	list of floats	List of sigmas to be used for the ridge plot. Default is [2.2, 2.2, 2.2, 2.2].	[2.2, 2.2, 2.2, 2.2]
mode	str	The mode of the filter. Can be "nearest", "reflect", "constant", "mirror", "wrap" or "linear". Default is "nearest".	'nearest'
alpha	float	The alpha value for the filter. Default is 0.5.	1.0
vmin	list of floats	List of vmin values for the ridge plot. Default is [None, None, None, None].	[None, None, None, None]
vmax	list of floats	List of vmax values for the ridge plot. Default is [None, None, None, None].	[None, None, None, None]
cmap_list	list of str	List of colormaps to be used for the ridge plot. Default is ["viridis", "viridis", "viridis", "viridis"].	['viridis', 'viridis', 'viridis', 'viridis']
draw_patch	list of bool	Whether to draw the circular patch. Default is [True, True, True, True].	[True, True, True, True]
draw_ridge	list of bool	Whether to draw the ridge line. Default is [True, True, True, True].	[True, True, True, True]
save_fig	bool	Whether to save the figure. Default is True.	True
fig_name	str	The name of the figure. Default is "new".	'new'
fig_format	str	The format of the figure. Default is "png".	'png'
c_label	list of str	List of colorbar labels. Default is [None, None, None, None].	[None, None, None, None]
wspace	float	The width space between subplots. Default is 0.1.	0.1
hspace	float	The height space between subplots. Default is 0.1.	0.1
fig_size	tuple of floats	The size of the figure. Default is (10, 10).	(10, 10)
box_style	dict	The style of the box. Default is None.	None
title_y_pos	float	The y-position of the title. Default is 0.95.	0.95
interpolation	str	The interpolation method for imshow. Default is "nearest". Options are "nearest", "bilinear", "bicubic", "spline16", "spline36", "hanning", "hamming", "hermite", "kaiser", "quadric", "catrom", "gaussian", "bessel", "mitchell"	'nearest'
tsy_model	str	The Tsyganenko model to be used. Default is "t96".	't96'
dark_mode	bool	Sets the dark mode for the plot and adjusts the color of labels and tickmarks accordingly. Default is True.	True
rc_file_name	str	The name of the file to save the reconnection line data. Default is "rc_file.csv".	'rc_file.csv'
rc_folder	str	The folder to save the reconnection line data. Default is "data".	'data'
save_rc_file	bool	Whether to save the reconnection line data. Default is False.	False
fig_version	str	The version of the figure. Default is "v001".	'v001'
df_jet_reversal	pandas DataFrame	The dataframe containing the jet reversal times. Default is None.	None

Raises:

Type	Description
ValueError: If the image is not a numpy array.

Returns:

Name	Type	Description
ridge_points	ndarray

Source code in src/rxn_location/rx_model_funcs.py

def ridge_finder_multiple(
    image=[None, None, None, None],
    convolution_order=[1, 1, 1, 1],
    t_range=["2016-12-24 15:08:00", "2016-12-24 15:12:00"],
    dt=5,
    b_imf=[-5, 0, 0],
    b_msh=[-5, 0, 0],
    v_msh=[-200, 50, 50],
    xrange=[-15.1, 15],
    yrange=[-15.1, 15],
    mms_probe_num="1",
    mms_sc_pos=[0, 0],
    dr=0.5,
    dipole_tilt_angle=None,
    p_dyn=None,
    imf_clock_angle=None,
    np_imf=None,
    v_imf=[None, None, None],
    sym_h=None,
    sigma=[2.2, 2.2, 2.2, 2.2],
    mode="nearest",
    alpha=1.0,
    vmin=[None, None, None, None],
    vmax=[None, None, None, None],
    cmap_list=["viridis", "viridis", "viridis", "viridis"],
    draw_patch=[True, True, True, True],
    draw_ridge=[True, True, True, True],
    save_fig=True,
    fig_name="new",
    fig_format="png",
    c_label=[None, None, None, None],
    wspace=0.1,
    hspace=0.1,
    fig_size=(10, 10),
    box_style=None,
    title_y_pos=0.95,
    interpolation="nearest",
    tsy_model="t96",
    dark_mode=True,
    rc_file_name="rc_file.csv",
    rc_folder="data",
    save_rc_file=False,
    fig_version="v001",
    df_jet_reversal=None,
    **kwargs,
):
    r"""
    Finds ridges in an image and plot the points with maximum ridge value on the given image.

    Parameters
    ----------
    image : list of numpy arrays
        List of images to be plotted.
    convolution_order : list of ints
        List of the order of the convolution to be used while smoothing each image. Values must be
        non-negative integers. Default is [1, 1, 1, 1].
    t_range : list of str
            The time range to find the ridge in. Default is ["2016-12-24 15:08:00",
            "2016-12-24 15:12:00"].
    dt : float, optional
        The time differential, in minutes, for observation if "t_range" has only one element.
        Default is 5 minutes.
    b_imf : list of floats, optional
            The IMF vector. Default is [-5, 0, 0].
    b_msh : list of floats, optional
            The MSH vector. Default is [-5, 0, 0].
    v_msh : list of floats, optional
            The MSH velocity vector. Default is [-200, 50, 50].
    xrange : list of floats, optional
            The range of x-values for image. Default is [-15.1, 15].
    yrange : list of floats, optional
            The range of y-values for image. Default is [-15.1, 15].
    mms_probe_num : str, optional
            The probe number of the MMS spacecraft. Default is 1.
    mms_sc_pos : list of floats, optional
            The position of the spacecraft in the image. Default is [0, 0].
    dr : float, optional
            The step size for the grid. Default is 0.5.
    dipole_tilt_angle : float, optional
            The dipole tilt angle. Default is None.
    p_dyn : float, optional
            The dynamic pressure. Default is None.
    imf_clock_angle : float, optional
            The IMF clock angle. Default is None.
    sigma : list of floats, optional
            List of sigmas to be used for the ridge plot. Default is [2.2, 2.2, 2.2, 2.2].
    mode : str
            The mode of the filter. Can be "nearest", "reflect", "constant", "mirror", "wrap" or
            "linear". Default is "nearest".
    alpha : float
            The alpha value for the filter. Default is 0.5.
    vmin : list of floats, optional
            List of vmin values for the ridge plot. Default is [None, None, None, None].
    vmax : list of floats, optional
            List of vmax values for the ridge plot. Default is [None, None, None, None].
    cmap_list : list of str, optional
            List of colormaps to be used for the ridge plot. Default is ["viridis", "viridis",
            "viridis", "viridis"].
    draw_patch : list of bool, optional
            Whether to draw the circular patch. Default is [True, True, True, True].
    draw_ridge : list of bool, optional
            Whether to draw the ridge line. Default is [True, True, True, True].
    save_fig : bool, optional
            Whether to save the figure. Default is True.
    fig_name : str, optional
            The name of the figure. Default is "new".
    fig_format : str, optional
            The format of the figure. Default is "png".
    c_label : list of str, optional
            List of colorbar labels. Default is [None, None, None, None].
    wspace : float, optional
            The width space between subplots. Default is 0.1.
    hspace : float, optional
            The height space between subplots. Default is 0.1.
    fig_size : tuple of floats, optional
            The size of the figure. Default is (10, 10).
    box_style : dict, optional
            The style of the box. Default is None.
    title_y_pos : float, optional
            The y-position of the title. Default is 0.95.
    interpolation : str, optional
            The interpolation method for imshow. Default is "nearest".
            Options are "nearest", "bilinear", "bicubic", "spline16", "spline36", "hanning",
            "hamming", "hermite", "kaiser", "quadric", "catrom", "gaussian", "bessel", "mitchell"
    tsy_model : str, optional
            The Tsyganenko model to be used. Default is "t96".
    dark_mode : bool, optional
        Sets the dark mode for the plot and adjusts the color of labels and tickmarks accordingly.
        Default is True.
    rc_file_name : str, optional
        The name of the file to save the reconnection line data. Default is "rc_file.csv".
    rc_folder : str, optional
        The folder to save the reconnection line data. Default is "data".
    save_rc_file : bool, optional
        Whether to save the reconnection line data. Default is False.
    fig_version : str, optional
        The version of the figure. Default is "v001".
    df_jet_reversal : pandas DataFrame, optional
        The dataframe containing the jet reversal times. Default is None.

    Raises
    ------
    ValueError: If the image is not a numpy array.

    Returns
    -------
    ridge_points : ndarray
    """
    if image is None:
        raise ValueError("No image given")

    if len(t_range) == 1:
        # Check if t_range is a datetime object
        if isinstance(t_range[0], datetime.datetime):
            t_range_date = t_range[0]
        else:
            t_range_date = datetime.datetime.strptime(t_range[0], "%Y-%m-%d %H:%M:%S")
        t_range_date_min = t_range_date - datetime.timedelta(minutes=dt)
        t_range_date_max = t_range_date + datetime.timedelta(minutes=dt)
        t_range = [
            t_range_date_min.strftime("%Y-%m-%d %H:%M:%S"),
            t_range_date_max.strftime("%Y-%m-%d %H:%M:%S"),
        ]

    if dark_mode:
        plt.style.use("dark_background")
        # tick_color = "w"  # color of the tick lines
        mtick_color = "w"  # color of the minor tick lines
        label_color = "w"  # color of the tick labels
        clabel_color = "w"  # color of the colorbar label
    else:
        plt.style.use("default")
        # tick_color = "k"  # color of the tick lines
        mtick_color = "k"  # color of the minor tick lines
        label_color = "k"  # color of the tick labels
        clabel_color = "k"  # color of the colorbar label

    # Set the fontstyle to Times New Roman
    font = {"family": "serif", "weight": "normal", "size": 10}
    plt.rc("font", **font)
    plt.rc("text", usetex=True)
    fig = plt.figure(num=None, figsize=fig_size, dpi=200, facecolor="w", edgecolor="k")
    fig.subplots_adjust(left=0.01, right=0.99, top=0.99, bottom=0.01, wspace=wspace, hspace=hspace)
    gs = gridspec.GridSpec(2, 2, width_ratios=[1, 1])

    # Set the font size for the axes
    label_size = 20  # fontsize for x and y labels
    t_label_size = 18  # fontsize for tick label
    c_label_size = 18  # fontsize for colorbar label
    ct_tick_size = 14  # fontsize for colorbar tick labels
    l_label_size = 14  # fontsize for legend label

    tick_len = 10  # length of the tick lines
    mtick_len = 7  # length of the minor tick lines
    tick_width = 1  # tick width in points
    mtick_width = 0.7  # minor tick width in points

    # box_style = dict(boxstyle="round", facecolor="wheat", alpha=0.5)
    if dark_mode:
        box_style = box_style
    else:
        box_style = dict(boxstyle="round", color="w", alpha=0.8, linewidth=1)
    y_vals = []
    x_intr_vals_list = []
    y_intr_vals_list = []
    dist_rc_list = []
    for i in range(len(image)):
        image_rotated = np.transpose(image[i])

        # Create the masked image from result for all the new processings
        # Find the number of rows in the original image
        n_rows, n_cols = image_rotated.shape

        # Make a grid of the data based on mumber of rows and columns
        X, Y = np.ogrid[:n_rows, :n_cols]

        # Find the central row and column
        c_row = int(n_rows / 2)
        c_col = int(n_cols / 2)
        # Find the distance of each pixel from the central pixel in terms of pixels
        dist_pxl = np.sqrt((X - c_row) ** 2 + (Y - c_col) ** 2)
        mask_image = dist_pxl > xrange[1] / dr

        if cmap_list is None:
            cmap_list = ["viridis", "viridis", "viridis", "viridis"]
        else:
            cmap_list = cmap_list
        if vmin is not None and vmax is not None:
            norm = plt.Normalize(vmin=vmin[i], vmax=vmax[i])
        else:
            norm = plt.Normalize()

        kwargs = {"sigmas": [sigma[i]], "black_ridges": False, "mode": mode, "alpha": 1}

        # Smoothen the image
        image_smooth = sp.ndimage.gaussian_filter(
            image_rotated, order=convolution_order[i], sigma=[5, 5], mode=mode
        )
        image_smooth_p = sp.ndimage.gaussian_filter(image_rotated, order=0, sigma=[5, 5], mode=mode)
        result = frangi(image_smooth, **kwargs)  # frangi, hessian, meijering, sato

        m_result = result.copy()
        m_result[mask_image] = np.nan
        new_image_rotated = image_rotated.copy()
        new_image_rotated[mask_image] = np.nan

        x_len = image_rotated.shape[0]
        y_len = image_rotated.shape[1]

        y_val = np.full(y_len, np.nan)
        y_vals.append(y_val)
        im_max_val = np.full(y_len, np.nan)
        for xx in range(y_len):
            try:
                y_val[xx] = np.nanargmax(m_result[:, xx]) * dr + yrange[0]
                im_max_val[xx] = np.nanargmax(new_image_rotated[:, xx]) * dr + yrange[0]
            except Exception:
                pass

        # TODO: Find a better way to do this
        j = i // 2
        k = i % 2

        axs1 = plt.subplot(gs[j, k])
        im1 = axs1.imshow(
            image_smooth_p,
            extent=[xrange[0], xrange[1], yrange[0], yrange[1]],
            origin="lower",
            cmap=cmap_list[i],
            norm=norm,
            interpolation=interpolation,
            alpha=1,
        )
        divider1 = make_axes_locatable(axs1)
        # Draw a circle of radius 15 (terminator) around the center of the image
        axs1.add_patch(plt.Circle((0, 0), radius=15, color="gray", fill=False, lw=0.5))

        # Take rolling average of the y_val array
        y_val_avg = np.full(len(y_val), np.nan)
        im_max_val_avg = np.full(len(y_val), np.nan)

        r_a_l = 5
        for xx in range(len(y_val)):
            y_val_avg[xx] = np.nanmean(y_val[max(0, xx - r_a_l) : min(len(y_val), xx + r_a_l)])
            im_max_val_avg[xx] = np.nanmean(
                im_max_val[max(0, xx - r_a_l) : min(len(y_val), xx + r_a_l)]
            )

        if draw_ridge:
            x_intr_vals = np.linspace(xrange[0], xrange[1], x_len)
            y_intr_vals = im_max_val_avg
            axs1.plot(x_intr_vals, y_intr_vals, color="aqua", ls="-", alpha=0.9)

        # Find the interpolation function corresponding to the x_vals and y_val_avg array
        line_intrp = line_fnc_der(x=np.linspace(xrange[0], xrange[1], x_len), y=im_max_val_avg)

        # Spacecraft position
        r0 = mms_sc_pos[:3]

        x_vals = np.linspace(xrange[0], xrange[1], x_len)

        curve = np.array([[x_vals[i], im_max_val_avg[i]] for i in range(len(x_vals))])
        points = np.array([r0[1:]])
        curves = np.array([curve], dtype=object)

        # compute unsigned distance
        dist_u = tt.basedists.distance(points, curves, argPnts=True)

        # Direction of the magnetosheath magnetic field at the position of the spacecraft
        # TODO: Check if this is same as direction/magnitude given by the Cooling model
        b_msh_dir = b_msh[:3] / np.linalg.norm(b_msh[:3])
        v_msh_dir = v_msh[:3] / np.linalg.norm(v_msh[:3])

        # Find the closest point on the reconnection line in the direction of the magnetosheath
        # magnetic field.
        xn = np.full(300, np.nan)
        yn = np.full(300, np.nan)
        for n in range(-150, 150):
            xn[50 + n] = r0[1] + n / 3 * b_msh_dir[1]
            yn[50 + n] = r0[2] + n / 3 * b_msh_dir[2]

        # Find the expected values of y-coordinate based on the x-coordinate, in the direction of
        # the magnetosheath magnetic field.
        yn_interp = line_intrp(xn)

        # Compute the distance between the spacecraft position and the coordinates computed in
        # previous step.
        dist_rn = np.abs(yn - yn_interp)
        # Find the index of the minimum distance
        min_dist_rn_idx = np.argmin(dist_rn)

        # Find the x- and y-coordinates corresponding to the minimum distance
        xn_rc = xn[min_dist_rn_idx]
        yn_rc = yn[min_dist_rn_idx]

        # Save the x- and y-coordinates of the reconnection line along with the spacecraft position
        x_intr_vals = [r0[1], xn_rc]
        y_intr_vals = [r0[2], yn_rc]

        # Append the x- and y-coordinates to the list in order to save to a file
        x_intr_vals_list.append(x_intr_vals)
        y_intr_vals_list.append(y_intr_vals)

        # Find the distance between the spacecraft position and the reconnection line
        dist_rc = dist_u[0]["UnsignedDistance"][0]
        x_y_point = dist_u[0]["ArgminPoints"][0]

        if dist_rc > xrange[1]:
            dist_rc = np.nan
        dist_rc_list.append(dist_rc)

        if i == 0:
            method_used = "shear"
        elif i == 1:
            method_used = "rx_en"
        elif i == 2:
            method_used = "va_cs"
        elif i == 3:
            method_used = "bisection"

        # Save the data to a text file
        # Create the rc-folder if it doesn"t exist
        if save_rc_file:
            if not os.path.exists(rc_folder):
                os.makedirs(rc_folder)
            var_list = (
                "mms_spc_num,date_from,date_to,spc_pos_x,spc_pos_y,spc_pos_z,"
                "b_msh_x,b_msh_y,b_msh_z,r_rc,method_used,b_imf_x,b_imf_y,"
                "b_imf_z,dipole,imf_clock_angle,p_dyn"
            )
            data_dict = {
                "mms_spc_num": mms_probe_num,
                "date_from": t_range[0],
                "date_to": t_range[1],
                "spc_pos_x": r0[0],
                "spc_pos_y": r0[1],
                "spc_pos_z": r0[2],
                "b_msh_x": b_msh[0],
                "b_msh_y": b_msh[1],
                "b_msh_z": b_msh[2],
                "r_rc": np.round(dist_rc, 3),
                "method_used": method_used,
                "b_imf_x": b_imf[0],
                "b_imf_y": b_imf[1],
                "b_imf_z": b_imf[2],
                "dipole": dipole_tilt_angle * 180 / np.pi,
                "imf_clock_angle": imf_clock_angle,
                "p_dyn": p_dyn,
            }
            # Add keys and data from df_jet_reversal to data_dict if those keys aren"t already
            # present in the dictionary
            try:
                for key in df_jet_reversal.keys():
                    if key not in data_dict.keys():
                        data_dict[key] = df_jet_reversal[key]
                        # Add the key to the variable list
                        var_list += "," + key
            except Exception:
                pass
            # Save data to the csv file using tab delimiter
            os.makedirs(rc_folder, exist_ok=True)
            if not os.path.exists(rc_folder + rc_file_name):
                with open(rc_folder + rc_file_name, "w") as f:
                    f.write(var_list + "\n")
                    f.close()
                    vprint(2, f"Created {rc_folder + rc_file_name} to store data")
            # Open file and append the relevant data
            with open(rc_folder + rc_file_name, "a") as f:
                for key in data_dict.keys():
                    try:
                        f.write(f"{np.round(data_dict[key], 3)},")
                    except Exception:
                        f.write(f"{data_dict[key]},")
                f.write("\n")
                f.close()
                vprint(2, f"Saved data to {rc_folder + rc_file_name}")

        # plot an arrow along the magnetosheath magnetic field direction
        axs1.arrow(
            r0[1],
            r0[2],
            5 * b_msh_dir[1],
            5 * b_msh_dir[2],
            head_width=0.4,
            head_length=0.7,
            fc="w",
            ec="r",
            linewidth=2,
            ls="-",
        )

        # Plot the arrow along the magnetosheath velocity direction
        axs1.arrow(
            r0[1],
            r0[2],
            5 * v_msh_dir[1],
            5 * v_msh_dir[2],
            head_width=0.4,
            head_length=0.7,
            fc="w",
            ec="b",
            linewidth=2,
            ls="-",
        )

        # Plot line connecting the spacecraft position and the reconnection line
        if ~np.isnan(dist_rc):
            # axs1.plot(x_intr_vals, y_intr_vals, "--", color="w", linewidth=2)
            axs1.plot([r0[1], x_y_point[0]], [r0[2], x_y_point[1]], "--", color="w", linewidth=2)
            distance = f"$R_{{\\rm rc}}$ = {dist_rc:.2f} $R_{{\\rm E}}$"
            axs1.text(
                x_intr_vals[0] - 6,
                y_intr_vals[0] + 2,
                distance,
                fontsize=l_label_size * 1.2,
                color="k",
                ha="left",
                va="bottom",
            )

        # Plot a horizontal line at x=0 and a vertical line at y=0
        axs1.axhline(0, color="k", linestyle="-", linewidth=0.5, alpha=0.5)
        axs1.axvline(0, color="k", linestyle="-", linewidth=0.5, alpha=0.5)

        if draw_patch:
            patch = patches.Circle(
                (0, 0), radius=xrange[1], transform=axs1.transData, fc="none", ec="k", lw=0.5
            )
            im1.set_clip_path(patch)
        axs1.add_patch(patch)
        if i == 0 or i == 2:
            axs1.set_ylabel(r"Z [GSM, $R_{\rm E}$]", fontsize=label_size, color=label_color)
        if i == 2 or i == 3:
            axs1.set_xlabel(r"Y [GSM, $R_{\rm E}$]", fontsize=label_size, color=label_color)
        if i == 1 or i == 3:
            axs1.set_ylabel(r"Z [GSM, $R_{\rm E}$]", fontsize=label_size, color=label_color)
            axs1.yaxis.set_label_position("right")

        if dark_mode:
            text_color = "white"
        else:
            text_color = "black"
        if i == 1:
            axs1.text(
                1.15,
                1.12,
                f"Model: {tsy_model}",
                horizontalalignment="right",
                verticalalignment="bottom",
                transform=axs1.transAxes,
                rotation=0,
                color=text_color,
                fontsize=l_label_size,
                bbox=box_style,
            )

        if i == 0:
            axs1.text(
                -0.18,
                1.12,
                f"MMS Position - [{mms_sc_pos[0]:.2f}, {mms_sc_pos[1]:.2f}, {mms_sc_pos[2]:.2f}] $R_E$ \n [GSM]",
                horizontalalignment="left",
                verticalalignment="bottom",
                transform=axs1.transAxes,
                rotation=0,
                color=text_color,
                fontsize=l_label_size,
                bbox=box_style,
            )

        # Define the location of the colorbar, it"s size relative to main figure and the padding
        # between the colorbar and the figure, the orientation the colorbar
        cax1 = divider1.append_axes("top", size="5%", pad=0.01)
        cbar1 = plt.colorbar(
            im1, cax=cax1, orientation="horizontal", ticks=None, fraction=0.05, pad=0.01
        )
        cbar1.ax.tick_params(
            axis="x",
            direction="in",
            top=True,
            labeltop=True,
            bottom=False,
            labelbottom=False,
            pad=0.01,
            labelsize=ct_tick_size,
            labelcolor=label_color,
        )
        cbar1.ax.xaxis.set_label_position("top")
        cbar1.ax.set_xlabel(f"{c_label[i]}", fontsize=c_label_size, color=clabel_color)
        # Draw the spacecraft position
        axs1.plot(mms_sc_pos[1], mms_sc_pos[2], "white", marker="$\\bigoplus$", ms=15, alpha=1)

        # Set tick label parameters
        if i == 0 or i == 2:
            axs1.tick_params(
                axis="both",
                direction="in",
                which="major",
                left=True,
                right=True,
                top=True,
                bottom=True,
                labelleft=True,
                labelright=False,
                labeltop=False,
                labelbottom=True,
                labelsize=t_label_size,
                length=tick_len,
                width=tick_width,
                labelcolor=label_color,
            )
        else:
            axs1.tick_params(
                axis="both",
                direction="in",
                which="major",
                left=True,
                right=True,
                top=True,
                bottom=True,
                labelleft=False,
                labelright=True,
                labeltop=False,
                labelbottom=True,
                labelsize=t_label_size,
                length=tick_len,
                width=tick_width,
                labelcolor=label_color,
            )
        # Write the timme range on the plot
        if i == 2:
            axs1.text(
                -0.17,
                -0.1,
                f"Clock Angle: {imf_clock_angle:.2f}$^\\circ$",
                horizontalalignment="left",
                verticalalignment="top",
                transform=axs1.transAxes,
                rotation=0,
                color=text_color,
                fontsize=l_label_size,
                bbox=box_style,
            )
        elif i == 3:
            axs1.text(
                1.17,
                -0.1,
                f"Dipole tilt: {dipole_tilt_angle * 180 / np.pi:.2f}$^\\circ$",
                horizontalalignment="right",
                verticalalignment="top",
                transform=axs1.transAxes,
                rotation=0,
                color=text_color,
                fontsize=l_label_size,
                bbox=box_style,
            )
        if i == 0:
            # Add a label "(a)" to the plot to indicate the panel number
            axs1.text(
                0.05,
                0.1,
                "(a)",
                horizontalalignment="left",
                verticalalignment="top",
                transform=axs1.transAxes,
                rotation=0,
                color=text_color,
                fontsize=1.2 * l_label_size,
            )
        elif i == 1:
            # Add a label "(b)" to the plot to indicate the panel number
            axs1.text(
                0.1,
                0.1,
                "(b)",
                horizontalalignment="right",
                verticalalignment="top",
                transform=axs1.transAxes,
                rotation=0,
                color=text_color,
                fontsize=1.2 * l_label_size,
            )
        elif i == 2:
            # Add a label "(c)" to the plot to indicate the panel number
            axs1.text(
                0.05,
                0.1,
                "(c)",
                horizontalalignment="left",
                verticalalignment="top",
                transform=axs1.transAxes,
                rotation=0,
                color=text_color,
                fontsize=1.2 * l_label_size,
            )
        elif i == 3:
            # Add a label "(d)" to the plot to indicate the panel number
            axs1.text(
                0.1,
                0.1,
                "(d)",
                horizontalalignment="right",
                verticalalignment="top",
                transform=axs1.transAxes,
                rotation=0,
                color=text_color,
                fontsize=1.2 * l_label_size,
            )

        # Show minor ticks
        axs1.minorticks_on()
        axs1.tick_params(
            axis="both",
            which="minor",
            direction="in",
            length=mtick_len,
            left=True,
            right=True,
            top=True,
            bottom=True,
            color=mtick_color,
            width=mtick_width,
        )
        # Set the number of ticks on the x-axis
        axs1.xaxis.set_major_locator(MaxNLocator(nbins=5, prune="lower"))
        # Set the number of ticks on the y-axis
        axs1.yaxis.set_major_locator(MaxNLocator(nbins=5, prune="lower"))

        # Setting the tickmarks labels in such a way that they don"t overlap
        plt.setp(axs1.get_xticklabels(), rotation=0, ha="right", va="top", visible=True)
        plt.setp(axs1.get_yticklabels(), rotation=0, va="center", visible=True)
        # Set the title of the plot
        fig.suptitle(
            f"Time range: {t_range[0]} - {t_range[1]} \n $B_{{\\rm {{imf}}}}$ = [{b_imf[0]:.2f}, {b_imf[1]:.2f}, {b_imf[2]:.2f}] nT",
            fontsize=label_size,
            color=text_color,
            y=title_y_pos,
            alpha=0.65,
        )

    if save_fig:
        try:
            # TODO: Add folder name as one of the path and make sure that the code creates the
            # folder. Gives out error if the folder can"t be created.
            temp1 = parser.parse(t_range[1]).strftime("%Y-%m-%d_%H-%M-%S")
            fig_time_range = f"{parser.parse(t_range[0]).strftime('%Y-%m-%d_%H-%M-%S')}_{temp1}"
            fig_folder = Path(
                f"figures/all_ridge_plots/{tsy_model}/{interpolation}"
                f"_interpolation_mms{mms_probe_num}/{fig_version}"
            )

            if not fig_folder.exists():
                fig_folder.mkdir(parents=True, exist_ok=True)
                vprint(2, "created folder : ", fig_folder)
            else:
                vprint(2, f"folder already exists: {fig_folder}\n")

            bbb = f"{b_imf[0]:.0f}_{b_imf[1]:.0f}_{b_imf[2]:.0f}"
            fig_name = fig_folder / f"ridge_plot_{fig_time_range}_{bbb}.{fig_format}"
            plt.savefig(str(fig_name), bbox_inches="tight", pad_inches=0.05, format=fig_format, dpi=200)
            abs_fig_name = os.path.abspath(str(fig_name))
            vprint(2, f"Figure saved as {abs_fig_name}")
        except Exception as e:
            vprint(1, e, color="red")
            vprint(2, "Figure not saved, folder does not exist. Create folder figures")
            vprint(2, "try again")
            # pass
        # plt.close()
    plt.close()

    # Build dictionary of R_rc values to return
    dist_rc_dict = {}
    for i, dist_rc in enumerate(dist_rc_list):
        method_used = c_label[i] if c_label[i] else f"model_{i}"
        key_name = f"r_rc_{method_used}"
        dist_rc_dict[key_name] = np.round(dist_rc, 3) if not np.isnan(dist_rc) else np.nan

    return y_vals, x_intr_vals_list, y_intr_vals_list, dist_rc_dict

ridge_finder_multiple_interactive(image=[None, None, None, None], convolution_order=[1, 1, 1, 1], t_range=['2016-12-24 15:08:00', '2016-12-24 15:12:00'], dt=5, b_imf=[-5, 0, 0], b_msh=[-5, 0, 0], v_msh=[-200, 50, 50], xrange=[-15.1, 15], yrange=[-15.1, 15], mms_probe_num='1', mms_sc_pos=[0, 0], dr=0.5, dipole_tilt_angle=None, p_dyn=None, imf_clock_angle=None, np_imf=None, v_imf=[None, None, None], sym_h=None, sigma=[2.2, 2.2, 2.2, 2.2], mode='nearest', alpha=1.0, vmin=[None, None, None, None], vmax=[None, None, None, None], cmap_list=['Viridis', 'Viridis', 'Viridis', 'Viridis'], draw_patch=[True, True, True, True], draw_ridge=[True, True, True, True], save_fig=True, fig_name='new', fig_format='html', c_label=[None, None, None, None], wspace=0.1, hspace=0.1, fig_size=(10, 10), box_style=None, title_y_pos=0.95, interpolation='nearest', tsy_model='t96', dark_mode=True, rc_file_name='rc_file.csv', rc_folder='data', save_rc_file=False, fig_version='v001', df_jet_reversal=None, **kwargs)

Finds and extracts the continuous ridge (X-line) locations across multiple reconnection models and prepares them for interactive Plotly 3D visualization.

This function evaluates the specified reconnection models over a 2D spatial grid (Y, Z) and extracts the coordinates of the maximal regions (ridges), representing the theoretical reconnection X-lines.

Parameters:

Name	Type	Description	Default
image	list	List of image arrays corresponding to each model.	[None, None, None, None]
convolution_order	list	Order of the convolution applied.	[1, 1, 1, 1]
t_range	list	Time range of the event.	['2016-12-24 15:08:00', '2016-12-24 15:12:00']
mms_probe_num	str	Probe number used for plotting title.	'1'
dr	float	Grid resolution.	0.5
mms_sc_pos	list	Spacecraft position vector [X, Y, Z].	[0, 0]

Returns:

Type	Description
list	A list of plotly trace objects representing the 3D X-line.

Source code in src/rxn_location/rx_model_funcs.py

def ridge_finder_multiple_interactive(
    image=[None, None, None, None],
    convolution_order=[1, 1, 1, 1],
    t_range=["2016-12-24 15:08:00", "2016-12-24 15:12:00"],
    dt=5,
    b_imf=[-5, 0, 0],
    b_msh=[-5, 0, 0],
    v_msh=[-200, 50, 50],
    xrange=[-15.1, 15],
    yrange=[-15.1, 15],
    mms_probe_num="1",
    mms_sc_pos=[0, 0],
    dr=0.5,
    dipole_tilt_angle=None,
    p_dyn=None,
    imf_clock_angle=None,
    np_imf=None,
    v_imf=[None, None, None],
    sym_h=None,
    sigma=[2.2, 2.2, 2.2, 2.2],
    mode="nearest",
    alpha=1.0,
    vmin=[None, None, None, None],
    vmax=[None, None, None, None],
    cmap_list=["Viridis", "Viridis", "Viridis", "Viridis"],
    draw_patch=[True, True, True, True],
    draw_ridge=[True, True, True, True],
    save_fig=True,
    fig_name="new",
    fig_format="html",
    c_label=[None, None, None, None],
    wspace=0.1,
    hspace=0.1,
    fig_size=(10, 10),
    box_style=None,
    title_y_pos=0.95,
    interpolation="nearest",
    tsy_model="t96",
    dark_mode=True,
    rc_file_name="rc_file.csv",
    rc_folder="data",
    save_rc_file=False,
    fig_version="v001",
    df_jet_reversal=None,
    **kwargs,
):
    """
    Finds and extracts the continuous ridge (X-line) locations across multiple 
    reconnection models and prepares them for interactive Plotly 3D visualization.

    This function evaluates the specified reconnection models over a 2D spatial grid
    (Y, Z) and extracts the coordinates of the maximal regions (ridges), representing
    the theoretical reconnection X-lines.

    Parameters
    ----------
    image : list
        List of image arrays corresponding to each model.
    convolution_order : list
        Order of the convolution applied.
    t_range : list
        Time range of the event.
    mms_probe_num : str
        Probe number used for plotting title.
    dr : float
        Grid resolution.
    mms_sc_pos : list
        Spacecraft position vector [X, Y, Z].

    Returns
    -------
    list
        A list of plotly trace objects representing the 3D X-line.
    """
    if image is None:
        raise ValueError("No image given")

    if len(t_range) == 1:
        if isinstance(t_range[0], datetime.datetime):
            t_range_date = t_range[0]
        else:
            t_range_date = datetime.datetime.strptime(t_range[0], "%Y-%m-%d %H:%M:%S")
        t_range_date_min = t_range_date - datetime.timedelta(minutes=dt)
        t_range_date_max = t_range_date + datetime.timedelta(minutes=dt)
        t_range = [
            t_range_date_min.strftime("%Y-%m-%d %H:%M:%S"),
            t_range_date_max.strftime("%Y-%m-%d %H:%M:%S"),
        ]

    # Plotly layout
    fig = make_subplots(rows=2, cols=2, horizontal_spacing=0.20, vertical_spacing=0.08)

    dist_rc_list = []

    y_vals = []
    x_intr_vals_list = []
    y_intr_vals_list = []

    for i in range(len(image)):
        image_rotated = np.transpose(image[i])
        n_rows, n_cols = image_rotated.shape
        X, Y = np.ogrid[:n_rows, :n_cols]
        c_row = int(n_rows / 2)
        c_col = int(n_cols / 2)
        dist_pxl = np.sqrt((X - c_row) ** 2 + (Y - c_col) ** 2)
        mask_image = dist_pxl > xrange[1] / dr

        kwargs = {"sigmas": [sigma[i]], "black_ridges": False, "mode": mode, "alpha": 1}
        image_smooth = sp.ndimage.gaussian_filter(
            image_rotated, order=convolution_order[i], sigma=[5, 5], mode=mode
        )
        image_smooth_p = sp.ndimage.gaussian_filter(image_rotated, order=0, sigma=[5, 5], mode=mode)
        result = frangi(image_smooth, **kwargs)

        m_result = result.copy()
        m_result[mask_image] = np.nan
        new_image_rotated = image_rotated.copy()
        new_image_rotated[mask_image] = np.nan
        image_smooth_p = image_smooth_p.astype(float)
        image_smooth_p[mask_image] = np.nan

        x_len = image_rotated.shape[0]
        y_len = image_rotated.shape[1]
        y_val = np.full(y_len, np.nan)
        y_vals.append(y_val)
        im_max_val = np.full(y_len, np.nan)

        for xx in range(y_len):
            try:
                y_val[xx] = np.nanargmax(m_result[:, xx]) * dr + yrange[0]
                im_max_val[xx] = np.nanargmax(new_image_rotated[:, xx]) * dr + yrange[0]
            except Exception:
                pass

        row = (i // 2) + 1
        col = (i % 2) + 1

        x_grid = np.linspace(xrange[0], xrange[1], x_len)
        y_grid = np.linspace(yrange[0], yrange[1], y_len)

        # Calculate dynamic colorbar positions to align with new tighter layout
        cb_x = 0.42 if col == 1 else 1.02
        cb_y = 0.77 if row == 1 else 0.23

        # Plotly Heatmap
        fig.add_trace(
            go.Heatmap(
                z=image_smooth_p,
                x=x_grid,
                y=y_grid,
                colorscale=cmap_list[i],
                showscale=True,
                colorbar=dict(
                    title=c_label[i] if c_label[i] else "",
                    orientation="v",
                    x=cb_x,
                    y=cb_y,
                    len=0.46,
                    thickness=15,
                ),
            ),
            row=row,
            col=col,
        )

        # Add circle outline (Terminator)
        if draw_patch[i]:
            fig.add_shape(
                type="circle",
                xref=f"x{i+1}",
                yref=f"y{i+1}",
                x0=-15,
                y0=-15,
                x1=15,
                y1=15,
                line=dict(color="gray", width=1),
            )

        y_val_avg = np.full(len(y_val), np.nan)
        im_max_val_avg = np.full(len(y_val), np.nan)
        r_a_l = 5
        for xx in range(len(y_val)):
            y_val_avg[xx] = np.nanmean(y_val[max(0, xx - r_a_l) : min(len(y_val), xx + r_a_l)])
            im_max_val_avg[xx] = np.nanmean(
                im_max_val[max(0, xx - r_a_l) : min(len(y_val), xx + r_a_l)]
            )

        if draw_ridge[i]:
            x_intr_vals = x_grid
            y_intr_vals = im_max_val_avg
            fig.add_trace(
                go.Scatter(
                    x=x_intr_vals,
                    y=y_intr_vals,
                    mode="lines",
                    line=dict(color="aqua", width=2),
                    showlegend=False,
                ),
                row=row,
                col=col,
            )

        r0 = mms_sc_pos[:3]

        fig.add_trace(
            go.Scatter(
                x=[r0[1]],
                y=[r0[2]],
                mode="markers",
                marker=dict(symbol="circle-cross", size=15, color="white"),
                showlegend=False,
            ),
            row=row,
            col=col,
        )

        fig.update_xaxes(title_text="Y [GSM, R_E]", range=xrange, row=row, col=col)
        fig.update_yaxes(
            title_text="Z [GSM, R_E]",
            range=yrange,
            row=row,
            col=col,
            scaleanchor=f"x{i+1}",
            scaleratio=1,
        )

        # Execute the heavy math to compute distance and directions
        line_intrp = line_fnc_der(x=x_grid, y=im_max_val_avg)

        curve = np.array([[x_grid[j], im_max_val_avg[j]] for j in range(len(x_grid))])
        points = np.array([r0[1:]])
        curves = np.array([curve], dtype=object)

        dist_u = tt.basedists.distance(points, curves, argPnts=True)

        b_msh_dir = b_msh[:3] / np.linalg.norm(b_msh[:3])
        v_msh_dir = v_msh[:3] / np.linalg.norm(v_msh[:3])

        xn = np.full(300, np.nan)
        yn = np.full(300, np.nan)
        for n in range(-150, 150):
            xn[50 + n] = r0[1] + n / 3 * b_msh_dir[1]
            yn[50 + n] = r0[2] + n / 3 * b_msh_dir[2]

        yn_interp = line_intrp(xn)
        dist_rn = np.abs(yn - yn_interp)
        min_dist_rn_idx = np.argmin(dist_rn)
        xn_rc = xn[min_dist_rn_idx]
        yn_rc = yn[min_dist_rn_idx]

        x_intr_vals = [r0[1], xn_rc]
        y_intr_vals = [r0[2], yn_rc]
        x_intr_vals_list.append(x_intr_vals)
        y_intr_vals_list.append(y_intr_vals)

        dist_rc = dist_u[0]["UnsignedDistance"][0]
        x_y_point = dist_u[0]["ArgminPoints"][0]

        if dist_rc > xrange[1]:
            dist_rc = np.nan
        dist_rc_list.append(dist_rc)

        # CSV saving moved outside the loop
        # Plotly Annotations for Arrows
        fig.add_annotation(
            x=r0[1] + 5 * b_msh_dir[1],
            y=r0[2] + 5 * b_msh_dir[2],
            ax=r0[1],
            ay=r0[2],
            xref=f"x{i+1}",
            yref=f"y{i+1}",
            axref=f"x{i+1}",
            ayref=f"y{i+1}",
            showarrow=True,
            arrowhead=2,
            arrowsize=1,
            arrowwidth=2,
            arrowcolor="red",
        )

        fig.add_annotation(
            x=r0[1] + 5 * v_msh_dir[1],
            y=r0[2] + 5 * v_msh_dir[2],
            ax=r0[1],
            ay=r0[2],
            xref=f"x{i+1}",
            yref=f"y{i+1}",
            axref=f"x{i+1}",
            ayref=f"y{i+1}",
            showarrow=True,
            arrowhead=2,
            arrowsize=1,
            arrowwidth=2,
            arrowcolor="blue",
        )

        if not np.isnan(dist_rc):
            fig.add_trace(
                go.Scatter(
                    x=[r0[1], x_y_point[0]],
                    y=[r0[2], x_y_point[1]],
                    mode="lines",
                    line=dict(color="white", width=2, dash="dash"),
                    showlegend=False,
                ),
                row=row,
                col=col,
            )
            # Add text annotation
            fig.add_annotation(
                x=r0[1] - 0.5,
                y=r0[2] + 0.5,
                xref=f"x{i+1}",
                yref=f"y{i+1}",
                text=f"R_rc = {dist_rc:.2f} R_E",
                showarrow=False,
                font=dict(color="white", size=13),
                bgcolor="rgba(0,0,0,0.6)",
                bordercolor="white",
                borderwidth=1,
                xanchor="right",
                yanchor="bottom",
            )

        # Plotly lines for axes
        fig.add_hline(y=0, line_dash="solid", line_color="gray", line_width=0.5, row=row, col=col)
        fig.add_vline(x=0, line_dash="solid", line_color="gray", line_width=0.5, row=row, col=col)

    if save_rc_file:
        os.makedirs(rc_folder, exist_ok=True)
        var_list = (
            "mms_spc_num,date_from,date_to,spc_pos_x,spc_pos_y,spc_pos_z,"
            "b_msh_x,b_msh_y,b_msh_z,b_imf_x,b_imf_y,"
            "b_imf_z,dipole,imf_clock_angle,p_dyn"
        )

        data_dict = {
            "mms_spc_num": mms_probe_num,
            "date_from": t_range[0],
            "date_to": t_range[1],
            "spc_pos_x": mms_sc_pos[0],
            "spc_pos_y": mms_sc_pos[1],
            "spc_pos_z": mms_sc_pos[2],
            "b_msh_x": b_msh[0],
            "b_msh_y": b_msh[1],
            "b_msh_z": b_msh[2],
            "b_imf_x": b_imf[0],
            "b_imf_y": b_imf[1],
            "b_imf_z": b_imf[2],
            "dipole": dipole_tilt_angle * 180 / np.pi if dipole_tilt_angle is not None else np.nan,
            "imf_clock_angle": imf_clock_angle if imf_clock_angle is not None else np.nan,
            "p_dyn": p_dyn if p_dyn is not None else np.nan,
        }

        try:
            if df_jet_reversal is not None:
                for key in df_jet_reversal.keys():
                    if key not in data_dict.keys():
                        data_dict[key] = df_jet_reversal[key]
                        var_list += "," + key
        except Exception:
            pass

        for i, dist_rc in enumerate(dist_rc_list):
            method_used = c_label[i] if c_label[i] else f"model_{i}"
            key_name = f"r_rc_{method_used}"
            data_dict[key_name] = np.round(dist_rc, 3)
            var_list += "," + key_name

        rc_path = os.path.join(rc_folder, rc_file_name)
        if not os.path.exists(rc_path):
            with open(rc_path, "w") as f:
                f.write(var_list + "\n")

        with open(rc_path, "a") as f:
            for key in data_dict.keys():
                try:
                    f.write(f"{np.round(float(data_dict[key]), 3)},")
                except Exception:
                    f.write(f"{data_dict[key]},")
            f.write("\n")

    # Define toggle button based on initial dark_mode
    if dark_mode:
        initial_label = "🌙"
        args1 = [
            {"template": pio.templates["plotly_white"], "updatemenus[0].buttons[0].label": "☀️"}
        ]
        args2 = [
            {"template": pio.templates["plotly_dark"], "updatemenus[0].buttons[0].label": "🌙"}
        ]
    else:
        initial_label = "☀️"
        args1 = [
            {"template": pio.templates["plotly_dark"], "updatemenus[0].buttons[0].label": "🌙"}
        ]
        args2 = [
            {"template": pio.templates["plotly_white"], "updatemenus[0].buttons[0].label": "☀️"}
        ]

    # Add center text annotation
    info_text = (
        f"Time range: {t_range[0]} - {t_range[1]}<br>"
        f"Model: {tsy_model.upper()}<br>"
        f"B_imf: [{b_imf[0]:.2f}, {b_imf[1]:.2f}, {b_imf[2]:.2f}] nT<br>"
        f"MMS Pos: [{mms_sc_pos[0]:.2f}, {mms_sc_pos[1]:.2f}, {mms_sc_pos[2]:.2f}] R_E<br>"
        f"Dipole Tilt: {dipole_tilt_angle * 180 / np.pi:.2f}°<br>"
        f"Clock Angle: {imf_clock_angle:.2f}°"
    )

    fig.add_annotation(
        text=info_text,
        xref="paper",
        yref="paper",
        x=0.5,
        y=1.15,
        xanchor="center",
        yanchor="top",
        showarrow=False,
        font=dict(size=14),
        bordercolor="gray",
        borderwidth=1,
        bgcolor="rgba(0,0,0,0.5)" if dark_mode else "rgba(255,255,255,0.8)",
        align="center",
    )

    # Update layout for interactive plot
    fig.update_layout(
        template="plotly_dark" if dark_mode else "plotly_white",
        width=1200,
        height=1300,
        margin=dict(t=180),
        autosize=True,
        updatemenus=[
            dict(
                type="buttons",
                direction="left",
                buttons=list(
                    [dict(args=args1, args2=args2, label=initial_label, method="relayout")]
                ),
                pad={"r": 10, "t": 10},
                showactive=False,
                x=1.0,
                xanchor="right",
                y=1.15,
                yanchor="top",
            )
        ],
    )

    if save_fig:
        try:
            temp1 = parser.parse(t_range[1]).strftime("%Y-%m-%d_%H-%M-%S")
            fig_time_range = f"{parser.parse(t_range[0]).strftime('%Y-%m-%d_%H-%M-%S')}_{temp1}"
            fig_folder = Path(
                f"interactive_figures/all_ridge_plots/{tsy_model}/{interpolation}"
                f"_interpolation_mms{mms_probe_num}/{fig_version}"
            )

            if not fig_folder.exists():
                fig_folder.mkdir(parents=True, exist_ok=True)

            bbb = f"{b_imf[0]:.0f}_{b_imf[1]:.0f}_{b_imf[2]:.0f}"
            fig_name = fig_folder / f"ridge_plot_{fig_time_range}_{bbb}.html"
            fig.write_html(str(fig_name), default_width="100%", default_height="100%")
            abs_fig_name = os.path.abspath(str(fig_name))
            vprint(2, f"Interactive figure saved as {abs_fig_name}")
        except Exception as e:
            vprint(1, e, color="red")

    # Build dictionary of R_rc values to return
    dist_rc_dict = {}
    for i, dist_rc in enumerate(dist_rc_list):
        method_used = c_label[i] if c_label[i] else f"model_{i}"
        key_name = f"r_rc_{method_used}"
        dist_rc_dict[key_name] = np.round(dist_rc, 3) if not np.isnan(dist_rc) else np.nan

    return fig, dist_rc_dict

rx_model(probe=None, trange=['2016-12-24 15:08:00', '2016-12-24 15:12:00'], dt=5, omni_level='hro', mms_probe_num='3', model_type='t96', m_p=0.5, dr=0.5, min_max_val=15, y_min=None, y_max=None, z_min=None, z_max=None, save_data=False, latest_version=False, nprocesses=None)

This function computes the magnetosheath and magnetospheric magnetic fields using the T96 model and solar wind paramters and returns the value of those fields as well as the value of other parameters associated with dayside reconnection, such as shear angle, reconnection energy, exhaust velocity, and the bisection field based on bisection model, for both the magnetospheric and magnetosheath bisection fields.

Parameters:

Name	Type	Description	Default
probe	str	The probe to be used for the computation. Options: "mms1", "mms2", "mms3", "mms4"	None
trange	list or an array of length 2	The time range to use. Should in the format [start, end], where start and end times should be a string in the format "YYYY-MM-DD HH:MM:SS".	['2016-12-24 15:08:00', '2016-12-24 15:12:00']
dt	float	The time interval to use in case the "trange" has only one element.	5
omni_level	str	The omni-data level to be used for the computation. Default is "hro". Options: "hro", "hro1"	'hro'
mms_probe_num	optional	The MMS probe number to be used for the computation. Default is "3". Options: "1", "2", "3", "4"	'3'
model_type	str	The model type to be used for the computation. Default is "t96". Options: "t96", "t01"	't96'
m_p	float	Thickness of the magnetosphere in earth radius units . Default is 0.5 R_E.	0.5
dr	float	Finess of the grid for model computation. Default is 0.5 R_E.	0.5
min_max_val	float	The minimum and maximum values of the y and z-axis to be used for the computation. Default is 15, meaning the model will be computed for -15 < y < 15 and -15 < z < 15.	15
y_min	float	The minimum value of the y-axis to be used for the computation. Default is None, in which case y_min is set to -min_max_val.	None
y_max	float	The maximum value of the y-axis to be used for the computation. Default is None, in which case y_max is set to min_max_val.	None
z_min	float	The minimum value of the z-axis to be used for the computation. Default is None, in which case z_min is set to -min_max_val.	None
z_max	float	The maximum value of the z-axis to be used for the computation. Default is None, in which case z_max is set to min_max_val.	None
save_data	bool	Whether to save the data or not. Default is False. If True, the data will be saved in a "HDF5" file.	False
latest_version	bool	Whether to use the latest version of the data. Default is False.	False
nprocesses	int	The number of processes to use for the computation. Default is None, in which case the number of processes will be set to the number of cores in the system.	None

Returns:

Type	Description
dict	A dictionary containing the computed solar wind parameters, input model grid coordinates, and the resulting arrays for magnetospheric and magnetosheath fields, shear angle, reconnection energy, exhaust velocity, and bisection fields.

Source code in src/rxn_location/rx_model_funcs.py

def rx_model(
    probe=None,
    trange=["2016-12-24 15:08:00", "2016-12-24 15:12:00"],
    dt=5,
    omni_level="hro",
    mms_probe_num="3",
    model_type="t96",
    m_p=0.5,
    dr=0.5,
    min_max_val=15,
    y_min=None,
    y_max=None,
    z_min=None,
    z_max=None,
    save_data=False,
    latest_version=False,
    nprocesses=None,
):
    """
    This function computes the magnetosheath and magnetospheric magnetic fields using the T96 model
    and  solar wind paramters and returns the value of those fields as well as the value of other
    parameters associated with dayside reconnection, such as shear angle, reconnection energy,
    exhaust velocity, and the bisection field based on bisection model, for both the magnetospheric
    and magnetosheath bisection fields.

    Parameters
    ----------
    probe : str
        The probe to be used for the computation.
        Options: "mms1", "mms2", "mms3", "mms4"
    trange : list or an array of length 2
        The time range to use. Should in the format [start, end], where start and end times should
        be a string in the format "YYYY-MM-DD HH:MM:SS".
    dt : float
        The time interval to use in case the "trange" has only one element.
    omni_level : str, optional
        The omni-data level to be used for the computation. Default is "hro".
        Options: "hro", "hro1"
    mms_probe_num : str. optional
        The MMS probe number to be used for the computation. Default is "3".
        Options: "1", "2", "3", "4"
    model_type : str
        The model type to be used for the computation. Default is "t96".
        Options: "t96", "t01"
    m_p : float
        Thickness of the magnetosphere in earth radius units . Default is 0.5 R_E.
    dr : float
        Finess of the grid for model computation. Default is 0.5 R_E.
    min_max_val : float
        The minimum and maximum values of the y and z-axis to be used for the computation. Default
        is 15, meaning the model will be computed for -15 < y < 15 and -15 < z < 15.
    y_min : float, optional
        The minimum value of the y-axis to be used for the computation. Default is None, in which
        case y_min is set to -min_max_val.
    y_max : float, optional
        The maximum value of the y-axis to be used for the computation. Default is None, in which
        case y_max is set to min_max_val.
    z_min : float, optional
        The minimum value of the z-axis to be used for the computation. Default is None, in which
        case z_min is set to -min_max_val.
    z_max : float, optional
        The maximum value of the z-axis to be used for the computation. Default is None, in which
        case z_max is set to min_max_val.
    save_data : bool, optional
        Whether to save the data or not. Default is False. If True, the data will be saved in a
        "HDF5" file.
    latest_version : bool, optional
        Whether to use the latest version of the data. Default is False.
    nprocesses : int, optional
        The number of processes to use for the computation. Default is None, in which case the
        number of processes will be set to the number of cores in the system.

    Returns
    -------
    dict
        A dictionary containing the computed solar wind parameters, input model grid coordinates,
        and the resulting arrays for magnetospheric and magnetosheath fields, shear angle,
        reconnection energy, exhaust velocity, and bisection fields.
    """

    if y_min is None:
        y_min = -min_max_val
    if y_max is None:
        y_max = min_max_val
    if z_min is None:
        z_min = -min_max_val
    if z_max is None:
        z_max = min_max_val

    # If trange has only one element, then make it a list of length 2 with "dt" minutes of padding
    if len(trange) == 1:
        # Check if trange is a datetime object
        if isinstance(trange[0], datetime.datetime):
            trange_date = trange[0]
        else:
            trange_date = datetime.datetime.strptime(trange[0], "%Y-%m-%d %H:%M:%S")
        trange_date_min = trange_date - datetime.timedelta(minutes=dt)
        trange_date_max = trange_date + datetime.timedelta(minutes=dt)
        trange = [
            trange_date_min.strftime("%Y-%m-%d %H:%M:%S"),
            trange_date_max.strftime("%Y-%m-%d %H:%M:%S"),
        ]
        # Add timezones to trange (UTC)
        trange = [trange[0] + "Z", trange[1] + "Z"]

    # Get the solar wind parameters for the model
    sw_params = get_sw_params(
        probe=probe,
        omni_level=omni_level,
        trange=trange,
        mms_probe_num=mms_probe_num,
        latest_version=latest_version,
        verbose=True,
    )

    n_arr_y = int((y_max - y_min) / dr) + 1
    n_arr_z = int((z_max - z_min) / dr) + 1
    bx = np.full((n_arr_y, n_arr_z), np.nan)
    by = np.full((n_arr_y, n_arr_z), np.nan)
    bz = np.full((n_arr_y, n_arr_z), np.nan)
    b_msx = np.full((n_arr_y, n_arr_z), np.nan)
    b_msy = np.full((n_arr_y, n_arr_z), np.nan)
    b_msz = np.full((n_arr_y, n_arr_z), np.nan)
    shear = np.full((n_arr_y, n_arr_z), np.nan)
    rx_en = np.full((n_arr_y, n_arr_z), np.nan)
    va_cs = np.full((n_arr_y, n_arr_z), np.nan)
    bisec_msp = np.full((n_arr_y, n_arr_z), np.nan)
    bisec_msh = np.full((n_arr_y, n_arr_z), np.nan)

    x_shu = np.full((n_arr_y, n_arr_z), np.nan)
    y_shu = np.full((n_arr_y, n_arr_z), np.nan)
    z_shu = np.full((n_arr_y, n_arr_z), np.nan)

    b_msx = np.full((n_arr_y, n_arr_z), np.nan)
    b_msy = np.full((n_arr_y, n_arr_z), np.nan)
    b_msz = np.full((n_arr_y, n_arr_z), np.nan)

    # Shue et al.,1998, equation 10
    ro = (10.22 + 1.29 * np.tanh(0.184 * (sw_params["b_imf"][2] + 8.14))) * (
        sw_params["p_dyn"]
    ) ** (-1.0 / 6.6)

    # Shue et al.,1998, equation 11
    alpha = (0.58 - 0.007 * sw_params["b_imf"][2]) * (1 + 0.024 * np.log(sw_params["p_dyn"]))
    rmp = ro * (2 / (1 + np.cos(0.0))) ** alpha  # Stand off position of the magnetopause

    len_y = int((y_max - y_min) / dr) + 1
    len_z = int((z_max - z_min) / dr) + 1

    if nprocesses is None:
        p = mp.Pool()
    else:
        p = mp.Pool(processes=nprocesses)

    input = (
        (j, k, y_max, z_max, dr, m_p, ro, alpha, rmp, sw_params, model_type)
        for j in range(len_y)
        for k in range(len_z)
    )

    vprint(2, "Running the model \n")
    res = p.map(model_run, input)
    vprint(2, "Model run complete \n")

    p.close()
    p.join()

    for r in res:
        j = r[0]
        k = r[1]
        bx[j, k] = r[2]
        by[j, k] = r[3]
        bz[j, k] = r[4]

        shear[j, k] = r[5]
        rx_en[j, k] = r[6]
        va_cs[j, k] = r[7]
        bisec_msp[j, k] = r[8]
        bisec_msh[j, k] = r[9]

        x_shu[j, k] = r[10]
        y_shu[j, k] = r[11]
        z_shu[j, k] = r[12]

        b_msx[j, k] = r[13]
        b_msy[j, k] = r[14]
        b_msz[j, k] = r[15]

    if save_data:
        try:
            today_date = datetime.datetime.today().strftime("%Y-%m-%d")
            fn = f"data/all_data_rx_model_{dr}re_{m_p}mp_{model_type}_{today_date}.h5"
            data_file = hf.File(fn, "w")

            data_file.create_dataset("bx", data=bx)
            data_file.create_dataset("by", data=by)
            data_file.create_dataset("bz", data=bz)

            data_file.create_dataset("b_msx", data=b_msx)
            data_file.create_dataset("b_msy", data=b_msy)
            data_file.create_dataset("b_msz", data=b_msz)

            data_file.create_dataset("shear", data=shear)
            data_file.create_dataset("rx_en", data=rx_en)
            data_file.create_dataset("va_cs", data=va_cs)
            data_file.create_dataset("bisec_msp", data=bisec_msp)
            data_file.create_dataset("bisec_msh", data=bisec_msh)

            data_file.create_dataset("x_shu", data=x_shu)
            data_file.create_dataset("y_shu", data=y_shu)
            data_file.create_dataset("z_shu", data=z_shu)

            data_file.close()
            vprint(2, f"Date saved to file {fn} \n")
        except Exception as e:
            vprint(1, e, color="red")
            vprint(2, 
                f"Data not saved to file {fn}. Please make sure that file name is correctly"
                + " assigned and that the directory exists and you have write permissions"
            )

    return (
        bx,
        by,
        bz,
        shear,
        rx_en,
        va_cs,
        bisec_msp,
        bisec_msh,
        sw_params,
        x_shu,
        y_shu,
        z_shu,
        b_msx,
        b_msy,
        b_msz,
    )

target_fnc(r, r0, b_msh, line_fnc, line_intrp)

Objective function used to minimize the distance between the theoretical X-line curve and the interpolated magnetic surface.

Parameters:

Name	Type	Description	Default
r	float	Current radial distance guess.	required
r0	float	Initial reference radial distance.	required
b_msh	array - like	Magnetosheath magnetic field vector.	required
line_fnc	function	Function describing the theoretical X-line geometry.	required
line_intrp	interp1d	Interpolation function of the evaluated magnetic surface.	required

Returns:

Type	Description
float	The squared difference to be minimized.

Source code in src/rxn_location/rx_model_funcs.py

def target_fnc(r, r0, b_msh, line_fnc, line_intrp):
    """
    Objective function used to minimize the distance between the theoretical 
    X-line curve and the interpolated magnetic surface.

    Parameters
    ----------
    r : float
        Current radial distance guess.
    r0 : float
        Initial reference radial distance.
    b_msh : array-like
        Magnetosheath magnetic field vector.
    line_fnc : function
        Function describing the theoretical X-line geometry.
    line_intrp : scipy.interpolate.interp1d
        Interpolation function of the evaluated magnetic surface.

    Returns
    -------
    float
        The squared difference to be minimized.
    """
    p_line = line_fnc(r0=r0, b_msh=b_msh, r=r)
    z_surface = line_intrp(p_line[1])
    return np.sum((p_line[2] - z_surface) ** 2)

rxn_location.jet_reversal_check_function

check_jet_location(df_mms=None, jet_len=3, time_cadence_median=0.15, v_thresh=70, ind_msh=None, verbose=True, ind_crossing=None, date_obs=None, save_delta_v_int=False, save_delta_v_stat=False, dark_mode=False)

Detects the presence and characteristics of a reconnection jet within the provided MMS data.

Parameters:

Name	Type	Description	Default
df_mms	DataFrame	Data frame containing the MMS data	None
jet_len	float	Length of the time for which jet must have threshold velocity to be considered as a jet, both at the positive and negative side of the crossing point	3
time_cadence_median	float	Median time cadence of the data in seconds	0.15
v_thresh	float	Threshold velocity in km/s. Default is 70 km/s	70
ind_msh	ndarray	Indices of the MSH data	None
verbose	bool	If True, prints the jet location and the jet length	True
ind_crossing	ndarray	Indices of the crossing points	None
date_obs	str	Date of observation in the format 'YYYY-MM-DD'	None

Returns:

Name	Type	Description
jet_present	bool	True if jet is present, False otherwise
delta_v_min	float	Difference between the velocity with respect to a reference magnetosheath velocity and the time series from 1 minute before/after the jet center to the point just before/after the jet starts/ends
delta_v_max	float	Difference between the velocity with respect to a reference magnetosheath velocity and the time series from 1 minute before/after the jet center to the point just before/after the jet starts/ends
t_jet_center	datetime	Time of the jet center in datetime format
ind_jet_center	int	Index of the jet center in the dataframe
ind_jet_center_minus_1_min	int	Index of the point 1 minute before the jet center in the dataframe
ind_jet_center_plus_1_min	int	Index of the point 1 minute after the jet center in the dataframe
vp_lmn_diff_l	ndarray	Difference between the velocity with respect to a reference magnetosheath velocity and the 'l' component of the velocity

Source code in src/rxn_location/jet_reversal_check_function.py

def check_jet_location(
    df_mms=None,
    jet_len=3,
    time_cadence_median=0.15,
    v_thresh=70,
    ind_msh=None,
    verbose=True,
    ind_crossing=None,
    date_obs=None,
    save_delta_v_int=False,
    save_delta_v_stat=False,
    dark_mode=False,
):
    """
    Detects the presence and characteristics of a reconnection jet within the provided MMS data.

    Parameters
    ----------
    df_mms : pandas.DataFrame
        Data frame containing the MMS data
    jet_len : float
        Length of the time for which jet must have threshold velocity to be considered as a jet,
        both at the positive and negative side of the crossing point
    time_cadence_median : float
        Median time cadence of the data in seconds
    v_thresh : float
        Threshold velocity in km/s. Default is 70 km/s
    ind_msh : numpy.ndarray
        Indices of the MSH data
    verbose : bool
        If True, prints the jet location and the jet length
    ind_crossing : numpy.ndarray
        Indices of the crossing points
    date_obs : str
        Date of observation in the format 'YYYY-MM-DD'

    Returns
    -------
    jet_present : bool
        True if jet is present, False otherwise
    delta_v_min : float
        Difference between the velocity with respect to a reference magnetosheath velocity and the
        time series from 1 minute before/after the jet center to the point just before/after the jet
        starts/ends
    delta_v_max : float
        Difference between the velocity with respect to a reference magnetosheath velocity and the
        time series from 1 minute before/after the jet center to the point just before/after the jet
        starts/ends
    t_jet_center : datetime.datetime
        Time of the jet center in datetime format
    ind_jet_center : int
        Index of the jet center in the dataframe
    ind_jet_center_minus_1_min: int
        Index of the point 1 minute before the jet center in the dataframe
    ind_jet_center_plus_1_min: int
        Index of the point 1 minute after the jet center in the dataframe
    vp_lmn_diff_l: np.ndarray
        Difference between the velocity with respect to a reference magnetosheath velocity and the
        'l' component of the velocity
    """
    # Compute the number of points corresponding to jet_len
    n_points_jet = max(1, int(jet_len / time_cadence_median))

    # Define the distance within which maximum and minimum values of jet velocity must lie
    # to be considered as a jet (within 30 seconds)
    delta_jet_min_max_ind = int(60 / time_cadence_median)
    # Get the median value of the velocity corresponding to the magnetosheath in the lmn coordinate
    # system
    vp_lmn_vec_msh_median = np.array(
        [
            np.nanmedian(df_mms["vp_lmn_l"].iloc[ind_msh]),
            np.nanmedian(df_mms["vp_lmn_m"].iloc[ind_msh]),
            np.nanmedian(df_mms["vp_lmn_n"].iloc[ind_msh]),
        ]
    )

    # Subtract the median value of the velocity corresponding to the magnetosheath from the velocity
    # in the lmn coordinate system
    vp_lmn_diff_l = df_mms["vp_lmn_l"] - vp_lmn_vec_msh_median[0]

    # Add the difference in the velocity in the lmn coordinate system to the dataframe
    df_mms["vp_lmn_diff_l"] = vp_lmn_diff_l

    print("\nvp_lmn_diff_l added to the dataframe\n")

    # Find the index where vp_lmn_diff_l has the maximum value
    ind_jet_max = np.argmax(vp_lmn_diff_l)
    t_jet_max = df_mms.index[ind_jet_max]

    # Within 30 seconds of ind_jet_max find the index where vp_lmn_diff_l has the minimum value
    # NOTE: The shifted indices are used to avoid the cases where the minimum (maximum) value of
    # vp_lmn_diff_l is at the beginning (end) of the dataframe
    shifted_ind_jet_min = np.max([ind_jet_max - delta_jet_min_max_ind, 0])
    shifted_ind_jet_max = np.min(
        [ind_jet_max + delta_jet_min_max_ind, len(vp_lmn_diff_l)]
    )
    ind_jet_min = shifted_ind_jet_min + np.argmin(
        vp_lmn_diff_l[shifted_ind_jet_min:shifted_ind_jet_max]
    )
    t_jet_min = df_mms.index[ind_jet_min]

    if t_jet_max < t_jet_min:
        # Find the time centered between t_jet_max and t_jet_min, which are datetime objects
        t_jet_center = t_jet_max + (t_jet_min - t_jet_max) / 2
        ind_jet_center = np.argmin(np.abs(df_mms.index - t_jet_center))

        # Find the time difference between t_jet_max and t_jet_center
        delta_t_jet_max_center = t_jet_center - t_jet_max
        # Find the number of points corresponding to delta_t_jet_max_center
        delta_n_jet_max_center = int(
            delta_t_jet_max_center.total_seconds() / time_cadence_median
        )

        # Find the time difference between t_jet_min and t_jet_center
        delta_t_jet_min_center = t_jet_min - t_jet_center
        # Find the number of points corresponding to delta_t_jet_min_center
        delta_n_jet_min_center = int(
            delta_t_jet_min_center.total_seconds() / time_cadence_median
        )

        # Find the median value of vp_lmn_diff_l from 1 minute before t_jet_center until just before
        # the jet starts
        t_jet_center_minus_1_min = t_jet_center - datetime.timedelta(minutes=1)
        ind_jet_center_minus_1_min = np.argmin(
            np.abs(df_mms.index - t_jet_center_minus_1_min)
        )
        v_max_median = np.nanmedian(
            vp_lmn_diff_l[
                ind_jet_center_minus_1_min : (
                    ind_jet_center - 2 * delta_n_jet_max_center
                )
            ]
        )

        # Find the median value of vp_lmn_diff_l between t_jet_center and 1 minute after
        # t_jet_center
        t_jet_center_plus_1_min = t_jet_center + datetime.timedelta(minutes=1)
        ind_jet_center_plus_1_min = np.argmin(
            np.abs(df_mms.index - t_jet_center_plus_1_min)
        )
        v_min_median = np.nanmedian(
            vp_lmn_diff_l[
                (
                    ind_jet_center + 2 * delta_n_jet_min_center
                ) : ind_jet_center_plus_1_min
            ]
        )

        # Subtract v_max_median from vp_lmn_diff_l from t_jet_cener_minus_1_min to
        # t_jet_center
        delta_v_max = (
            vp_lmn_diff_l[ind_jet_center_minus_1_min:ind_jet_center] - v_max_median
        )

        # Subtract v_min_median from vp_lmn_diff_l from t_jet_center to t_jet_center_plus_1_min
        delta_v_min = (
            vp_lmn_diff_l[ind_jet_center:ind_jet_center_plus_1_min] - v_min_median
        )

        # Check if delta_v_max has a sustained value greater than v_thresh for n_points_jet
        # points
        ind_v_gt_vthresh = np.flatnonzero(
            np.convolve(
                delta_v_max >= v_thresh, np.ones(n_points_jet, dtype=int), "valid"
            )
            >= n_points_jet
        )
        # Check if delta_v_min has a sustained value less than -v_thresh for n_points_jet
        # points
        ind_v_lt_vthresh = np.flatnonzero(
            np.convolve(
                delta_v_min <= -v_thresh, np.ones(n_points_jet, dtype=int), "valid"
            )
            >= n_points_jet
        )
        # If both conditions are satisfied then a jet is found
        if len(ind_v_gt_vthresh) > 0 and len(ind_v_lt_vthresh) > 0:
            jet_detection = True
        else:
            jet_detection = False
            # ind_jet_center = []
            # delta_v_min = []
            # delta_v_max = []
    else:
        t_jet_center = t_jet_min + (t_jet_max - t_jet_min) / 2
        ind_jet_center = np.argmin(np.abs(df_mms.index - t_jet_center))

        # Find the time difference between t_jet_max and t_jet_center
        delta_t_jet_max_center = t_jet_max - t_jet_center
        # Find the number of points corresponding to delta_t_jet_max_center
        delta_n_jet_max_center = int(
            delta_t_jet_max_center.total_seconds() / time_cadence_median
        )

        # Find the time difference between t_jet_min and t_jet_center
        delta_t_jet_min_center = t_jet_center - t_jet_min
        # Find the number of points corresponding to delta_t_jet_min_center
        delta_n_jet_min_center = int(
            delta_t_jet_min_center.total_seconds() / time_cadence_median
        )

        # Find the median value of vp_lmn_diff_l from 1 minute before t_jet_center to just before
        # the jet starts
        t_jet_center_minus_1_min = t_jet_center - datetime.timedelta(minutes=1)
        ind_jet_center_minus_1_min = np.argmin(
            np.abs(df_mms.index - t_jet_center_minus_1_min)
        )
        v_min_median = np.nanmedian(
            vp_lmn_diff_l[
                ind_jet_center_minus_1_min : (
                    ind_jet_center - 2 * delta_n_jet_min_center
                )
            ]
        )

        # Find the median value of vp_lmn_diff_l between just following the jet and 1 minute after
        # t_jet_center
        t_jet_center_plus_1_min = t_jet_center + datetime.timedelta(minutes=1)
        ind_jet_center_plus_1_min = np.argmin(
            np.abs(df_mms.index - t_jet_center_plus_1_min)
        )
        v_max_median = np.nanmedian(
            vp_lmn_diff_l[
                (
                    ind_jet_center + 2 * delta_n_jet_max_center
                ) : ind_jet_center_plus_1_min
            ]
        )

        # Subtract v_min_median from vp_lmn_diff_l from t_jet_cener_minus_1_min to
        # t_jet_center
        delta_v_min = (
            vp_lmn_diff_l[ind_jet_center_minus_1_min:ind_jet_center] - v_min_median
        )

        # Subtract v_max_median from vp_lmn_diff_l from t_jet_center to t_jet_center_plus_1_min
        delta_v_max = (
            vp_lmn_diff_l[ind_jet_center:ind_jet_center_plus_1_min] - v_max_median
        )

        # Check if delta_v_min has a sustained value less than -v_thresh for n_points_jet
        # points
        ind_v_lt_vthresh = np.flatnonzero(
            np.convolve(
                delta_v_min <= -v_thresh, np.ones(n_points_jet, dtype=int), "valid"
            )
            >= n_points_jet
        )
        # Check if delta_v_max has a sustained value greater than v_thresh for n_points_jet
        # points
        ind_v_gt_vthresh = np.flatnonzero(
            np.convolve(
                delta_v_max >= v_thresh, np.ones(n_points_jet, dtype=int), "valid"
            )
            >= n_points_jet
        )

        # If both conditions are satisfied then a jet is found
        if len(ind_v_lt_vthresh) > 0 and len(ind_v_gt_vthresh) > 0:
            jet_detection = True
        else:
            jet_detection = False
            # ind_jet_center = []
            # delta_v_min = []
            # delta_v_max = []
    if verbose:
        if jet_detection:
            print(f"\n\033[1;32m Jet found at {t_jet_center}\033[0m \n")
        else:
            print(f"\n\033[1;31m No jet found at {t_jet_center}\033[0m \n")

    if save_delta_v_stat or save_delta_v_int:
        if jet_detection:
            save_folder = f"figures/jet_reversal_checks/check_{date_obs}/delta_v/jet/"
        else:
            save_folder = f"figures/jet_reversal_checks/check_{date_obs}/delta_v/no_jet/"

        if not os.path.exists(save_folder):
            os.makedirs(save_folder)

        # Matplotlib Static Plot
        if save_delta_v_stat:
            plt.figure(figsize=(10, 5))
            plt.plot(df_mms.index, vp_lmn_diff_l, "r-", alpha=0.3)
            plt.plot(df_mms.index[ind_jet_center:ind_jet_center_plus_1_min], delta_v_max, "b-", alpha=1, label="$\\Delta v_{\\rm max}$")
            plt.plot(df_mms.index[ind_jet_center_minus_1_min:ind_jet_center], delta_v_min, "g-", alpha=1, label="$\\Delta v_{\\rm min}$")
            plt.axvline(t_jet_center, color="k", linestyle="--", alpha=0.5)
            plt.axvline(t_jet_center_minus_1_min, color="m", linestyle="--", alpha=0.5)
            plt.axvline(t_jet_center_plus_1_min, color="g", linestyle="--", alpha=0.5)
            plt.axhline(v_thresh, color="k", linestyle="--", alpha=0.5)
            plt.axhline(-v_thresh, color="k", linestyle="--", alpha=0.5)
            plt.title(f"Different delta as a function of time at {t_jet_center.strftime('%Y-%m-%d %H:%M:%S')}")
            plt.xlabel("Time [UTC]")
            plt.ylabel("$\\Delta V$ [km/s]")
            plt.legend(loc=1)

            color = "g" if jet_detection else "r"
            plt.text(0.05, 0.95, ind_crossing, horizontalalignment="left", verticalalignment="top", transform=plt.gca().transAxes, color=color)

            fig_name = f"{save_folder}/delta_v_{t_jet_center.strftime('%Y-%m-%d_%H-%M-%S')}.png"
            plt.savefig(
                fig_name, dpi=300, bbox_inches="tight", pad_inches=0.1,
                transparent=True, facecolor="w", edgecolor="w", orientation="landscape",
            )
            plt.close()

        # Plotly Interactive Plot
        if save_delta_v_int:
            import plotly.graph_objects as go
            fig = go.Figure()
            fig.add_trace(go.Scatter(x=df_mms.index, y=vp_lmn_diff_l, mode='lines', line=dict(color='red'), opacity=0.3, name='vp_lmn_diff_l'))
            fig.add_trace(go.Scatter(x=df_mms.index[ind_jet_center:ind_jet_center_plus_1_min], y=delta_v_max, mode='lines', line=dict(color='blue'), name='Δv_max'))
            fig.add_trace(go.Scatter(x=df_mms.index[ind_jet_center_minus_1_min:ind_jet_center], y=delta_v_min, mode='lines', line=dict(color='green'), name='Δv_min'))

            fig.add_vline(x=t_jet_center, line_dash="dash", line_color="black", opacity=0.5)
            fig.add_vline(x=t_jet_center_minus_1_min, line_dash="dash", line_color="magenta", opacity=0.5)
            fig.add_vline(x=t_jet_center_plus_1_min, line_dash="dash", line_color="green", opacity=0.5)

            fig.add_hline(y=v_thresh, line_dash="dash", line_color="black", opacity=0.5)
            fig.add_hline(y=-v_thresh, line_dash="dash", line_color="black", opacity=0.5)

            # Add annotation for crossing index
            fig.add_annotation(
                x=0.05, y=0.95, xref="paper", yref="paper", text=str(ind_crossing), 
                showarrow=False, font=dict(color="green" if jet_detection else "red")
            )

            fig.update_layout(
                title=f"Different delta as a function of time at {t_jet_center.strftime('%Y-%m-%d %H:%M:%S')}",
                xaxis_title="Time [UTC]",
                yaxis_title="ΔV [km/s]",
                template="plotly_dark" if dark_mode else "plotly_white"
            )

            fig_name_html = f"{save_folder}/delta_v_{t_jet_center.strftime('%Y-%m-%d_%H-%M-%S')}.html"
            fig.write_html(fig_name_html)

    return (
        jet_detection,
        delta_v_min,
        delta_v_max,
        t_jet_center,
        ind_jet_center,
        ind_jet_center_minus_1_min,
        ind_jet_center_plus_1_min,
        vp_lmn_diff_l,
    )

check_msp_msh_location(df_mms=None, time_cadence_median=0.15, verbose=True)

Determine the data indices corresponding to the magnetosheath (MSH) and magnetosphere (MSP) based on proton density thresholds.

Parameters:

Name	Type	Description	Default
df_mms	DataFrame	Dataframe containing the mms data	None
time_cadence_median	float	Median time cadence of the data	0.15
verbose	bool	If True, print information	True

Returns:

Name	Type	Description
ind_range_msh	ndarray or None	Indices representing the location of the magnetosheath, or None if not found.
ind_range_msp	ndarray or None	Indices representing the location of the magnetosphere, or None if not found.

Source code in src/rxn_location/jet_reversal_check_function.py

def check_msp_msh_location(df_mms=None, time_cadence_median=0.15, verbose=True):
    """
    Determine the data indices corresponding to the magnetosheath (MSH) and magnetosphere (MSP) based on proton density thresholds.

    Parameters
    ----------
    df_mms : pandas.DataFrame
        Dataframe containing the mms data
    time_cadence_median : float
        Median time cadence of the data
    verbose : bool
        If True, print information

    Returns
    -------
    ind_range_msh : numpy.ndarray or None
        Indices representing the location of the magnetosheath, or None if not found.
    ind_range_msp : numpy.ndarray or None
        Indices representing the location of the magnetosphere, or None if not found.
    """
    # TODO: Check if threshold value of 5 and 10 are fine or if we need to decrease/increase it
    n_thresh_msp = 5
    n_thresh_msh = 10

    ind_range_msp = None
    ind_range_msh = None

    # Compute the median time difference between the points
    # NOTE: The reason for factor 10**9 is to convert the time difference to seconds from
    # nanoseconds which is the format in which df_mms.index.astype(np.int64) outputs time
    time_cadence_median = np.median(np.diff(df_mms.index.astype(np.int64) / 10**9))

    # TODO: Instead of one difference between indices, use the median value of differences between
    # all the indices
    n_points_msp = int(5 / time_cadence_median)
    n_points_msh = int(5 / time_cadence_median)

    # Find an interval of length at least 'n_points_msp_msh' where 'np' is greater than 'n_thresh'
    # on both sides of the minimum value
    np_msp_bool_array = df_mms["np"] < n_thresh_msp
    np_msh_bool_array = df_mms["np"] > n_thresh_msh

    # ind_np_msp_vals = np.flatnonzero(np.convolve(np_msp_bool_array > 0,
    #                                              np.ones(
    #                                                  n_points_msp, dtype=int),
    #                                              'valid') >= n_points_msp)

    result_msp = list(mit.run_length.encode(np_msp_bool_array))
    # Find the length of longest subsequence of True, and the location if that index in result
    max_true_count_msp = -1
    max_true_idx_msp = -1
    for idx, (val, count) in enumerate(result_msp):
        if val and max_true_count_msp < count:
            max_true_count_msp = count
            max_true_idx_msp = idx
    # Find total elements before and after the longest subsequence tuple
    elems_before_idx_msp = sum((idx[1] for idx in result_msp[:max_true_idx_msp]))
    # elems_after_idx_msp = sum((idx[1] for idx in result_msp[max_true_idx_msp + 1:]))

    # Check if the longest subsequence is greater than the threshold
    if max_true_count_msp >= n_points_msp:
        ind_min_msp = int(
            elems_before_idx_msp + (max_true_count_msp - n_points_msp) / 2
        )
        ind_max_msp = int(
            elems_before_idx_msp + (max_true_count_msp + n_points_msp) / 2
        )
        ind_range_msp = np.arange(ind_min_msp, ind_max_msp)

    # ind_np_msh_vals = np.flatnonzero(np.convolve(np_msh_bool_array > 0,
    #                                  np.ones(n_points_msh, dtype=int),
    #                                  'valid') >= n_points_msh)

    result_msh = list(mit.run_length.encode(np_msh_bool_array))
    # Find the length of longest subsequence of True, and the location if that index in result
    max_true_count_msh = -1
    max_true_idx_msh = -1
    for idx, (val, count) in enumerate(result_msh):
        if val and max_true_count_msh < count:
            max_true_count_msh = count
            max_true_idx_msh = idx
    # Find total elements before and after the longest subsequence tuple
    elems_before_idx_msh = sum((idx[1] for idx in result_msh[:max_true_idx_msh]))
    # Check if the longest subsequence is greater than the threshold
    if max_true_count_msh >= n_points_msh:
        ind_min_msh = int(
            elems_before_idx_msh + (max_true_count_msh - n_points_msh) / 2
        )
        ind_max_msh = int(
            elems_before_idx_msh + (max_true_count_msh + n_points_msh) / 2
        )
        ind_range_msh = np.arange(ind_min_msh, ind_max_msh)

    if verbose:
        try:
            print(f"ind_min_msp: {ind_min_msp}")
            print(f"ind_max_msp: {ind_max_msp}")
            print(f"ind_min_msh: {ind_min_msh}")
            print(f"ind_max_msh: {ind_max_msh}")
        except Exception:
            pass

    return ind_range_msp, ind_range_msh

jet_reversal_check(crossing_time=None, dt=90, probe=3, data_rate='fast', level='l2', coord_type='lmn', data_type='dis-moms', time_clip=True, latest_version=False, jet_len=5, figname=None, date_obs=None, fname='data/mms_jet_reversal_times.csv', error_file_log_name='data/mms_jet_reversal_check_error_log.csv', verbose=True, return_plotly_fig=False, dark_mode=False, save_jet_int=True, save_jet_stat=True, save_delta_v_int=False, save_delta_v_stat=False)

For a given crossing time and a given probe, the function finds out if MMS observed a jet during magnetopause crossing. If there was indeed a jet reversal, then the function saves that time to a csv file, along with the probe number, position of the spacecraft, and the time of the crossing.

Parameters:

Name	Type	Description	Default
crossing_time	datetime	The time of the crossing of the magnetopause.	None
dt	int	The time interval to look for a jet reversal.	90
probe	int	The probe number. Can be 1, 2, 3, or 4. Default is 3. (since the magnetopause crossing times are for MMS3)	3
data_rate	str	The data rate. Can be 'fast' or 'srvy' or 'brst. Default is 'fast'.	'fast'
level	str	The data level. Can be either 'l1' or 'l2'. Default is 'l2'.	'l2'
coord_type	str	The coordinate type. Can be either 'lmn' or 'gse'. Default is 'lmn'.	'lmn'
data_type	str	The data type. Can be either 'dis-moms' or 'des-moms'. Default is 'dis-moms'.	'dis-moms'
time_clip	bool	Whether or not to clip the data to the time range of the crossing time. Default is True.	True
latest_version	bool	Whether or not to use the latest version of the data. Default is True.	False
jet_len	int	The time length for which the jet should be observed. Default is 5.	5
figname	str	The name of the figure to be saved. Default is 'mms_jet_reversal_check'.	None
date_obs	str	The date of the observation. Default is None.	None
fname	str	The name of the csv file to save the data to. Default is 'mms_jet_reversal_times.csv'.	'data/mms_jet_reversal_times.csv'
error_file_log_name	str	The name of the csv file to save the error log to. Default is 'data/mms_jet_reversal_check_error_log.csv'.	'data/mms_jet_reversal_check_error_log.csv'
verbose	bool	Whether or not to print the status of the function. Default is True.	True

Returns:

Name	Type	Description
fig	Figure or Figure	The generated figure showing the MMS data and jet reversal.
jet_detection	bool	True if a jet was detected during the magnetopause crossing, False otherwise.
data_dict	dict	A dictionary containing extracted parameters at the jet center, including time, position, shear angle, and field parameters.

Source code in src/rxn_location/jet_reversal_check_function.py

def jet_reversal_check(
    crossing_time=None,
    dt=90,
    probe=3,
    data_rate="fast",
    level="l2",
    coord_type="lmn",
    data_type="dis-moms",
    time_clip=True,
    latest_version=False,
    jet_len=5,
    figname=None,
    date_obs=None,
    fname="data/mms_jet_reversal_times.csv",
    error_file_log_name="data/mms_jet_reversal_check_error_log.csv",
    verbose=True,
    return_plotly_fig=False,
    dark_mode=False,
    save_jet_int=True,
    save_jet_stat=True,
    save_delta_v_int=False,
    save_delta_v_stat=False,
):
    """
    For a given crossing time and a given probe, the function finds out if MMS observed a jet
    during magnetopause crossing. If there was indeed a jet reversal, then the function saves that
    time to a csv file, along with the probe number, position of the spacecraft, and the time of
    the crossing.

    Parameters
    ----------
    crossing_time : datetime.datetime
        The time of the crossing of the magnetopause.
    dt : int
        The time interval to look for a jet reversal.
    probe : int
        The probe number. Can be 1, 2, 3, or 4. Default is 3. (since the magnetopause crossing
        times are for MMS3)
    data_rate : str
        The data rate. Can be 'fast' or 'srvy' or 'brst. Default is 'fast'.
    level : str
        The data level. Can be either 'l1' or 'l2'. Default is 'l2'.
    coord_type : str
        The coordinate type. Can be either 'lmn' or 'gse'. Default is 'lmn'.
    data_type : str
        The data type. Can be either 'dis-moms' or 'des-moms'. Default is 'dis-moms'.
    time_clip : bool
        Whether or not to clip the data to the time range of the crossing time. Default is True.
    latest_version : bool
        Whether or not to use the latest version of the data. Default is True.
    jet_len : int
        The time length for which the jet should be observed. Default is 5.
    figname : str
        The name of the figure to be saved. Default is 'mms_jet_reversal_check'.
    date_obs : str
        The date of the observation. Default is None.
    fname : str
        The name of the csv file to save the data to. Default is 'mms_jet_reversal_times.csv'.
    error_file_log_name : str
        The name of the csv file to save the error log to. Default is
        'data/mms_jet_reversal_check_error_log.csv'.
    verbose : bool
        Whether or not to print the status of the function. Default is True.

    Returns
    -------
    fig : matplotlib.figure.Figure or plotly.graph_objects.Figure
        The generated figure showing the MMS data and jet reversal.
    jet_detection : bool
        True if a jet was detected during the magnetopause crossing, False otherwise.
    data_dict : dict
        A dictionary containing extracted parameters at the jet center, including time, position, shear angle, and field parameters.
    """

    # Define the absolute permeability of free space in m^2 kg^-1 s^-1
    mu_0 = 4 * np.pi * 1e-7

    # Define the Boltzmann constant in J K^-1
    k_B = 1.38064852e-23

    if date_obs is None and crossing_time is not None:
        date_obs = crossing_time.strftime("%Y%m%d")

    # Define the conversion factor for converting temperature from ev to K
    ev_to_K = 11604.525

    crossing_time_min = crossing_time - datetime.timedelta(seconds=dt)
    crossing_time_max = crossing_time + datetime.timedelta(seconds=dt)
    trange = [crossing_time_min, crossing_time_max]

    # Get the index corresponding to the crossing time in the data
    try:
        df_crossing_temp = pd.read_csv("data/brst_intervals.csv", index_col=False)
        df_crossing_temp["start_time"] = pd.to_datetime(df_crossing_temp["start_time"])
        ind_crossing = np.where(df_crossing_temp["start_time"] == crossing_time)[0][0]
    except Exception:
        ind_crossing = 0
    print("Crossing index:", ind_crossing)

    # Get the data from the FPI
    mms_fpi_varnames = [
        f"mms{probe}_dis_numberdensity_{data_rate}",
        f"mms{probe}_dis_bulkv_gse_{data_rate}",
        f"mms{probe}_dis_temppara_{data_rate}",
        f"mms{probe}_dis_tempperp_{data_rate}",
        f"mms{probe}_dis_energyspectr_omni_{data_rate}",
        f"mms{probe}_des_energyspectr_omni_{data_rate}",
    ]

    _ = mms.fpi(
        trange=trange,
        probe=probe,
        data_rate=data_rate,
        level=level,
        datatype=[data_type, "des-moms"] if isinstance(data_type, str) else data_type,
        varnames=mms_fpi_varnames,
        time_clip=time_clip,
        latest_version=latest_version,
    )

    mms_fpi_time_unix = spd.get_data(mms_fpi_varnames[0])[0]
    # Convert the time to a datetime object
    mms_fpi_time_local = pd.to_datetime(mms_fpi_time_unix, unit="s")
    mms_fpi_time = mms_fpi_time_local.tz_localize(None).tz_localize(pytz.utc)

    mms_fpi_numberdensity = spd.get_data(mms_fpi_varnames[0])[1]
    _ = spd.get_data(mms_fpi_varnames[1])[1:4][0]
    mms_fpi_temppara = spd.get_data(mms_fpi_varnames[2])[1]
    mms_fpi_tempperp = spd.get_data(mms_fpi_varnames[3])[1]

    # Store both the temperatures in the ptt
    spd.store_data(
        "Tp",
        data=[
            f"mms{probe}_dis_temppara_{data_rate}",
            f"mms{probe}_dis_tempperp_{data_rate}",
        ],
    )

    # Covert gse to gsm
    _ = spd.cotrans(
        name_in=f"mms{probe}_dis_bulkv_gse_{data_rate}",
        name_out=f"mms{probe}_dis_bulkv_gsm_{data_rate}",
        coord_in="gse",
        coord_out="gsm",
    )

    mms_fpi_bulkv_gsm = spd.get_data(f"mms{probe}_dis_bulkv_gsm_{data_rate}")[1:4][0]
    mms_fpi_bulkv_gse = spd.get_data(f"mms{probe}_dis_bulkv_gse_{data_rate}")[1:4][0]

    if coord_type == "lmn":
        # Convert gsm to lmn
        _ = mms_cotrans_lmn.mms_cotrans_lmn(
            name_in=f"mms{probe}_dis_bulkv_gsm_{data_rate}",
            name_out=f"mms{probe}_dis_bulkv_lmn_{data_rate}",
            gse=False,
            gsm=True,
            probe=str(probe),
            data_rate="srvy",
        )

        mms_fpi_bulkv_lmn = spd.get_data(f"mms{probe}_dis_bulkv_lmn_{data_rate}")[1:4][
            0
        ]

    # Create a dataframe with the FPI data
    if coord_type == "lmn":
        df_mms_fpi = pd.DataFrame(
            index=mms_fpi_time,
            data={
                "np": mms_fpi_numberdensity,
                "vp_lmn_l": mms_fpi_bulkv_lmn[:, 0],
                "vp_lmn_m": mms_fpi_bulkv_lmn[:, 1],
                "vp_lmn_n": mms_fpi_bulkv_lmn[:, 2],
                "tp_para": mms_fpi_temppara,
                "tp_perp": mms_fpi_tempperp,
            },
        )
    else:
        df_mms_fpi = pd.DataFrame(
            index=mms_fpi_time,
            data={
                "np": mms_fpi_numberdensity,
                "vp_gsm_x": mms_fpi_bulkv_gsm[:, 0],
                "vp_gsm_y": mms_fpi_bulkv_gsm[:, 1],
                "vp_gsm_z": mms_fpi_bulkv_gsm[:, 2],
                "vp_gse_x": mms_fpi_bulkv_gse[:, 0],
                "vp_gse_y": mms_fpi_bulkv_gse[:, 1],
                "vp_gse_z": mms_fpi_bulkv_gse[:, 2],
                "tp_para": mms_fpi_temppara,
                "tp_perp": mms_fpi_tempperp,
            },
        )

    # Make sure that the time indices are in increasing order
    df_mms_fpi = df_mms_fpi.sort_index()

    if verbose:
        print("\n\033[1;32m FPI dataframe created \033[0m \n")
        print(f"The fpi Datafram:\n {df_mms_fpi.head()}")

    # Add rolling median to the dataframe
    df_mms_fpi["np_rolling_median"] = (
        df_mms_fpi["np"].rolling("60s", center=True).median()
    )

    if coord_type == "lmn":
        df_mms_fpi["vp_lmn_l_rolling_median"] = (
            df_mms_fpi["vp_lmn_l"].rolling("60s", center=True).median()
        )
        df_mms_fpi["vp_lmn_m_rolling_median"] = (
            df_mms_fpi["vp_lmn_m"].rolling("60s", center=True).median()
        )
        df_mms_fpi["vp_lmn_n_rolling_median"] = (
            df_mms_fpi["vp_lmn_n"].rolling("60s", center=True).median()
        )

    else:
        df_mms_fpi["vp_gsm_x_rolling_median"] = (
            df_mms_fpi["vp_gsm_x"].rolling("60s", center=True).median()
        )
        df_mms_fpi["vp_gsm_y_rolling_median"] = (
            df_mms_fpi["vp_gsm_y"].rolling("60s", center=True).median()
        )
        df_mms_fpi["vp_gsm_z_rolling_median"] = (
            df_mms_fpi["vp_gsm_z"].rolling("60s", center=True).median()
        )

    df_mms_fpi["tp_para_rolling_median"] = (
        df_mms_fpi["tp_para"].rolling("60s", center=True).median()
    )
    df_mms_fpi["tp_perp_rolling_median"] = (
        df_mms_fpi["tp_perp"].rolling("60s", center=True).median()
    )

    if coord_type == "lmn":
        df_mms_fpi["vp_diff_x"] = (
            df_mms_fpi["vp_lmn_l"] - df_mms_fpi["vp_lmn_l_rolling_median"]
        )
        df_mms_fpi["vp_diff_y"] = (
            df_mms_fpi["vp_lmn_m"] - df_mms_fpi["vp_lmn_m_rolling_median"]
        )
        df_mms_fpi["vp_diff_z"] = (
            df_mms_fpi["vp_lmn_n"] - df_mms_fpi["vp_lmn_n_rolling_median"]
        )
        # df_mms_fpi['vp_diff_z'] = df_mms_fpi['vp_lmn_l'] - np.nanmean(df_mms_fpi['vp_lmn_l'])
    else:
        df_mms_fpi["vp_diff_x"] = (
            df_mms_fpi["vp_gsm_x"] - df_mms_fpi["vp_gsm_x_rolling_median"]
        )
        df_mms_fpi["vp_diff_y"] = (
            df_mms_fpi["vp_gsm_y"] - df_mms_fpi["vp_gsm_y_rolling_median"]
        )
        df_mms_fpi["vp_diff_z"] = (
            df_mms_fpi["vp_gsm_z"] - df_mms_fpi["vp_gsm_z_rolling_median"]
        )

    # Get the data from the FGM
    _ = mms.fgm(
        trange=trange,
        probe=probe,
        time_clip=time_clip,
        latest_version=False,
        get_fgm_ephemeris=True,
    )
    # Get the time corresponding to the FGM data
    mms_fgm_time = spd.get_data(f"mms{probe}_fgm_b_gsm_srvy_{level}")[0]
    # Convert the time to a datetime object
    mms_fgm_time = pd.to_datetime(mms_fgm_time, unit="s")
    mms_fgm_time = mms_fgm_time.tz_localize(None).tz_localize(pytz.utc)

    mms_fgm_b_gsm = spd.get_data(f"mms{probe}_fgm_b_gsm_srvy_{level}")[1:4][0]
    mms_fgm_b_gse = spd.get_data(f"mms{probe}_fgm_b_gse_srvy_{level}")[1:4][0]
    mms_fgm_r_gsm = spd.get_data(f"mms{probe}_fgm_r_gsm_srvy_{level}")[1:4][0]

    if coord_type == "lmn":
        _ = mms_cotrans_lmn.mms_cotrans_lmn(
            name_in=f"mms{probe}_fgm_b_gsm_srvy_{level}_bvec",
            name_out=f"mms{probe}_fgm_b_lmn_srvy_{level}",
            gsm=True,
            probe=str(probe),
            data_rate="srvy",
        )

        mms_fgm_b_lmn = spd.get_data(f"mms{probe}_fgm_b_lmn_srvy_{level}")[1:4][0]

    # Create a dataframe with the FGM data
    if coord_type == "lmn":
        df_mms_fgm = pd.DataFrame(
            index=mms_fgm_time,
            data={
                "b_lmn_l": mms_fgm_b_lmn[:, 0],
                "b_lmn_m": mms_fgm_b_lmn[:, 1],
                "b_lmn_n": mms_fgm_b_lmn[:, 2],
            },
        )
    else:
        df_mms_fgm = pd.DataFrame(
            index=mms_fgm_time,
            data={
                "b_gsm_x": mms_fgm_b_gsm[:, 0],
                "b_gsm_y": mms_fgm_b_gsm[:, 1],
                "b_gsm_z": mms_fgm_b_gsm[:, 2],
                "b_gse_x": mms_fgm_b_gse[:, 0],
                "b_gse_y": mms_fgm_b_gse[:, 1],
                "b_gse_z": mms_fgm_b_gse[:, 2],
            },
        )

    # Make sure all time indices are in increasing order
    df_mms_fgm = df_mms_fgm.sort_index()

    if verbose:
        print("\n\033[1;32m FGM dataframe created \033[0m \n")
        print(f"FGM Dataframe: \n{df_mms_fgm.head()} \n")
    # Merge the two dataframes
    df_mms = pd.merge_asof(df_mms_fpi, df_mms_fgm, left_index=True, right_index=True)

    # Compute the median time difference between the points
    # NOTE: The reason for factor 10**9 is to convert the time difference to seconds from
    # nanoseconds which is the format in which df_mms.index.astype(np.int64) outputs time
    time_cadence_median = np.median(np.diff(df_mms.index.astype(np.int64) / 10**9))

    # Make the time index timezone aware
    try:
        df_mms.index = df_mms.index.tz_localize(None).tz_localize(pytz.utc)
    except Exception:
        if verbose:
            print("\033[1;31m Timezone conversion failed \033[0m \n")

    # Check if magnetosheath and magnetopshere are present and if they are then get the indices
    # corresponding to the magnetosheath and magnetosphere
    ind_range_msp, ind_range_msh = check_msp_msh_location(
        df_mms=df_mms, time_cadence_median=time_cadence_median, verbose=verbose
    )

    if ind_range_msp is None or ind_range_msh is None:
        if verbose:
            print(
                f"\n\033[1;31m Magnetosphere or magnetosheath region not found for date {crossing_time}\033[0m \n"
            )
        return None

    (
        jet_detection,
        delta_v_min,
        delta_v_max,
        t_jet_center,
        ind_jet_center,
        ind_jet_center_minus_1_min,
        ind_jet_center_plus_1_min,
        vp_lmn_diff_l,
    ) = check_jet_location(
        df_mms=df_mms,
        jet_len=jet_len,
        v_thresh=70,
        ind_msh=ind_range_msh,
        time_cadence_median=time_cadence_median,
        verbose=verbose,
        ind_crossing=ind_crossing,
        date_obs=date_obs,
        save_delta_v_int=save_delta_v_int,
        save_delta_v_stat=save_delta_v_stat,
        dark_mode=dark_mode,
    )

    # Add vp_lmn_diff_l to ptt FIRST so the pseudo-variable can find it!
    spd.store_data("vp_lmn_diff_l", data={"x": mms_fpi_time_unix, "y": vp_lmn_diff_l})

    if jet_detection:
        t_delta_v_min = mms_fpi_time_unix[ind_jet_center_minus_1_min:ind_jet_center]
        delta_v_min_data = vp_lmn_diff_l[ind_jet_center_minus_1_min:ind_jet_center]

        t_delta_v_max = mms_fpi_time_unix[ind_jet_center:ind_jet_center_plus_1_min]
        delta_v_max_data = vp_lmn_diff_l[ind_jet_center:ind_jet_center_plus_1_min]

        spd.store_data("delta_v_min", data={"x": t_delta_v_min, "y": delta_v_min_data})
        spd.store_data("delta_v_max", data={"x": t_delta_v_max, "y": delta_v_max_data})

        spd.store_data(
            "delta_v_vp_lmn_diff_l",
            data=["vp_lmn_diff_l", "delta_v_min", "delta_v_max"],
        )

    # Get different parameters for magnetosphere and magnetosheath
    np_msp = df_mms["np"].iloc[ind_range_msp] * 1e6  # Convert to m^-3 from cm^-3
    np_msh = df_mms["np"].iloc[ind_range_msh] * 1e6  # Convert to m^-3 from cm^-3

    if coord_type == "lmn":

        vp_lmn_vec_msp = (
            np.array(
                [
                    df_mms["vp_lmn_n"].iloc[ind_range_msp],
                    df_mms["vp_lmn_m"].iloc[ind_range_msp],
                    df_mms["vp_lmn_l"].iloc[ind_range_msp],
                ]
            )
            * 1e3
        )  # Convert to m/s
        vp_lmn_vec_msp = vp_lmn_vec_msp.T

        vp_lmn_vec_msh = (
            np.array(
                [
                    df_mms["vp_lmn_n"].iloc[ind_range_msh],
                    df_mms["vp_lmn_m"].iloc[ind_range_msh],
                    df_mms["vp_lmn_l"].iloc[ind_range_msh],
                ]
            )
            * 1e3
        )  # Convert to m/s
        vp_lmn_vec_msh = vp_lmn_vec_msh.T

        b_lmn_vec_msp = (
            np.array(
                [
                    df_mms["b_lmn_n"].iloc[ind_range_msp],
                    df_mms["b_lmn_m"].iloc[ind_range_msp],
                    df_mms["b_lmn_l"].iloc[ind_range_msp],
                ]
            )
            * 1e-9
        )  # Convert to T from nT
        b_lmn_vec_msp = b_lmn_vec_msp.T
        b_lmn_vec_msh = (
            np.array(
                [
                    df_mms["b_lmn_n"].iloc[ind_range_msh],
                    df_mms["b_lmn_m"].iloc[ind_range_msh],
                    df_mms["b_lmn_l"].iloc[ind_range_msh],
                ]
            )
            * 1e-9
        )  # Convert to T from nT
        b_lmn_vec_msh = b_lmn_vec_msh.T

        # Get the mean and median values of np, vp, and b for the magnetosphere and magnetosheath
        np_msp_median = np.nanmedian(np_msp) / 1e6  # Convert to cm^-3 from m^-3
        np_msh_median = np.nanmedian(np_msh) / 1e6  # Convert to cm^-3 from m^-3
        np_msp_mean = np.nanmean(np_msp) / 1e6  # Convert to cm^-3 from m^-3
        np_msh_mean = np.nanmean(np_msh) / 1e6  # Convert to cm^-3 from m^-3

        vp_lmn_vec_msp_median = np.nanmedian(vp_lmn_vec_msp, axis=0)  # km/sec
        vp_lmn_vec_msh_median = np.nanmedian(vp_lmn_vec_msh, axis=0)  # km/sec
        vp_lmn_vec_msp_mean = np.nanmean(vp_lmn_vec_msp, axis=0)  # km/sec
        vp_lmn_vec_msh_mean = np.nanmean(vp_lmn_vec_msh, axis=0)  # km/sec

        b_lmn_vec_msp_median = np.nanmedian(b_lmn_vec_msp, axis=0)  # T
        b_lmn_vec_msh_median = np.nanmedian(b_lmn_vec_msh, axis=0)  # T
        b_lmn_vec_msp_mean = np.nanmean(b_lmn_vec_msp, axis=0)  # T
        b_lmn_vec_msh_mean = np.nanmean(b_lmn_vec_msh, axis=0)  # T

        # Get the angle between the two vectors
        angle_b_lmn_vec_msp_msh_median = (
            np.arccos(
                np.dot(b_lmn_vec_msp_median, b_lmn_vec_msh_median)
                / (
                    np.linalg.norm(b_lmn_vec_msp_median)
                    * np.linalg.norm(b_lmn_vec_msh_median)
                )
            )
            * 180
            / np.pi
        )
    else:
        vp_gse_vec_msp = (
            np.array(
                [
                    df_mms["vp_gse_x"].iloc[ind_range_msp],
                    df_mms["vp_gse_y"].iloc[ind_range_msp],
                    df_mms["vp_gse_z"].iloc[ind_range_msp],
                ]
            )
            * 1e3
        )  # Convert km/s==> m/s
        vp_gse_vec_msp = vp_gse_vec_msp.T
        vp_gse_vec_msh = (
            np.array(
                [
                    df_mms["vp_gse_x"].iloc[ind_range_msh],
                    df_mms["vp_gse_y"].iloc[ind_range_msh],
                    df_mms["vp_gse_z"].iloc[ind_range_msh],
                ]
            )
            * 1e3
        )  # Convert km/s==> m/s
        vp_gse_vec_msh = vp_gse_vec_msh.T
        b_gse_vec_msp = (
            np.array(
                [
                    df_mms["b_gse_x"].iloc[ind_range_msp],
                    df_mms["b_gse_y"].iloc[ind_range_msp],
                    df_mms["b_gse_z"].iloc[ind_range_msp],
                ]
            )
            * 1e-9
        )  # Convert to T from nT
        b_gse_vec_msp = b_gse_vec_msp.T
        b_gse_vec_msh = (
            np.array(
                [
                    df_mms["b_gse_x"].iloc[ind_range_msh],
                    df_mms["b_gse_y"].iloc[ind_range_msh],
                    df_mms["b_gse_z"].iloc[ind_range_msh],
                ]
            )
            * 1e-9
        )  # Convert to T from nT
        b_gse_vec_msh = b_gse_vec_msh.T

    tp_para_msp = (
        df_mms["tp_para"].iloc[ind_range_msp] * ev_to_K
    )  # Convert to K from ev
    tp_para_msh = (
        df_mms["tp_para"].iloc[ind_range_msh] * ev_to_K
    )  # Convert to K from ev
    tp_perp_msp = (
        df_mms["tp_perp"].iloc[ind_range_msp] * ev_to_K
    )  # Convert to K from ev
    tp_perp_msh = (
        df_mms["tp_perp"].iloc[ind_range_msh] * ev_to_K
    )  # Convert to K from ev

    # Get the mean and median values of temperature for the magnetosphere and magnetosheath
    tp_para_msp_median = np.nanmedian(tp_para_msp)
    tp_para_msh_median = np.nanmedian(tp_para_msh)
    tp_para_msp_mean = np.nanmean(tp_para_msp)
    tp_para_msh_mean = np.nanmean(tp_para_msh)
    tp_perp_msp_median = np.nanmedian(tp_perp_msp)
    tp_perp_msh_median = np.nanmedian(tp_perp_msh)
    tp_perp_msp_mean = np.nanmean(tp_perp_msp)
    tp_perp_msh_mean = np.nanmean(tp_perp_msh)

    # Compute the magnetosheath beta value
    beta_msh_mean = (
        2
        * mu_0
        * np_msh_mean
        * 1e6
        * k_B
        * (2 * tp_para_msh_mean + tp_perp_msh_mean)
        / (3 * np.linalg.norm(b_lmn_vec_msh_mean) ** 2)
    )
    beta_msp_mean = (
        2
        * mu_0
        * np_msp_mean
        * 1e6
        * k_B
        * (2 * tp_para_msp_mean + tp_perp_msp_mean)
        / (3 * np.linalg.norm(b_lmn_vec_msp_mean) ** 2)
    )

    # Check if within 2 minutes of crossing time the values went above and below the threshold
    # If ind_vals is not empty, then append the crossing time to the csv file
    # if len(ind_vals) > 0:
    data_dict = None
    if jet_detection:
        # Position of MMS in GSM coordinates in earth radii (r_e) units
        r_e = 6378.137  # Earth radius in km
        # _ = spd.mms.mec
        # mms_sc_pos = spd.get_data(mms_mec_varnames[0])[1:3][0][0] / r_e
        x = np.nanmean(mms_fgm_r_gsm[:, 0]) / r_e
        y = np.nanmean(mms_fgm_r_gsm[:, 1]) / r_e
        z = np.nanmean(mms_fgm_r_gsm[:, 2]) / r_e
        r_yz = np.sqrt(y**2 + z**2)  # Projection distance in yz plane.

        # TODO: Add magnetic beta for sheath
        # List of variables to be saved in the csv file
        var_list = (
            "Date,Probe,jet_detection,x_gsm,y_gsm,z_gsm,r_spc,"
            "jet_time,ind_min_msp,ind_max_msp,ind_min_msh,ind_max_msh,ind_jet_center,"
            "angle_b_lmn_vec_msp_msh_median,b_lmn_vec_msp_mean_n,b_lmn_vec_msp_mean_m,"
            "b_lmn_vec_msp_mean_l,b_lmn_vec_msp_median_n,b_lmn_vec_msp_median_m,"
            "b_lmn_vec_msp_median_l,b_lmn_vec_msh_mean_n,b_lmn_vec_msh_mean_m,"
            "b_lmn_vec_msh_mean_l,b_lmn_vec_msh_median_n,b_lmn_vec_msh_median_m,"
            "b_lmn_vec_msh_median_l,np_msp_median,np_msp_mean,np_msh_median,np_msh_mean,"
            "vp_lmn_vec_msp_mean_n,vp_lmn_vec_msp_mean_m,vp_lmn_vec_msp_mean_l,"
            "vp_lmn_vec_msp_median_n,vp_lmn_vec_msp_median_m,vp_lmn_vec_msp_median_l,"
            "vp_lmn_vec_msh_mean_n,vp_lmn_vec_msh_mean_m,vp_lmn_vec_msh_mean_l,"
            "vp_lmn_vec_msh_median_n,vp_lmn_vec_msh_median_m,vp_lmn_vec_msh_median_l,"
            "tp_para_msp_median,tp_para_msh_median,tp_para_msp_mean,tp_para_msh_mean,"
            "tp_perp_msp_median,tp_perp_msh_median,tp_perp_msp_mean,tp_perp_msh_mean,"
            "beta_msh_mean,beta_msp_mean"
        )

        data_dict = {
            "Date": crossing_time,
            "Probe": probe,
            "jet_detection": jet_detection,
            "x_gsm": np.round(x, 3),
            "y_gsm": np.round(y, 3),
            "z_gsm": np.round(z, 3),
            "r_spc": np.round(r_yz, 3),
            "jet_time": t_jet_center,
            "ind_min_msp": ind_range_msp[0],
            "ind_max_msp": ind_range_msp[-1],
            "ind_min_msh": ind_range_msh[0],
            "ind_max_msh": ind_range_msh[-1],
            "ind_jet_center": ind_jet_center,
            "angle_b_lmn_vec_msp_msh_median": np.round(
                angle_b_lmn_vec_msp_msh_median, 3
            ),
            "b_lmn_vec_msp_mean_n": np.round(b_lmn_vec_msp_mean[0] * 1e9, 3),
            "b_lmn_vec_msp_mean_m": np.round(b_lmn_vec_msp_mean[1] * 1e9, 3),
            "b_lmn_vec_msp_mean_l": np.round(b_lmn_vec_msp_mean[2] * 1e9, 3),
            "b_lmn_vec_msp_median_n": np.round(b_lmn_vec_msp_median[0] * 1e9, 3),
            "b_lmn_vec_msp_median_m": np.round(b_lmn_vec_msp_median[1] * 1e9, 3),
            "b_lmn_vec_msp_median_l": np.round(b_lmn_vec_msp_median[2] * 1e9, 3),
            "b_lmn_vec_msh_mean_n": np.round(b_lmn_vec_msh_median[0] * 1e9, 3),
            "b_lmn_vec_msh_mean_m": np.round(b_lmn_vec_msh_median[1] * 1e9, 3),
            "b_lmn_vec_msh_mean_l": np.round(b_lmn_vec_msh_median[2] * 1e9, 3),
            "b_lmn_vec_msh_median_n": np.round(b_lmn_vec_msh_median[0] * 1e9, 3),
            "b_lmn_vec_msh_median_m": np.round(b_lmn_vec_msh_median[1] * 1e9, 3),
            "b_lmn_vec_msh_median_l": np.round(b_lmn_vec_msh_median[2] * 1e9, 3),
            "np_msp_median": np.round(np_msp_median, 3),
            "np_msp_mean": np.round(np_msp_mean, 3),
            "np_msh_median": np.round(np_msh_median, 3),
            "np_msh_mean": np.round(np_msh_mean, 3),
            "vp_lmn_vec_msp_mean_n": np.round(vp_lmn_vec_msp_mean[0] / 1e3, 3),
            "vp_lmn_vec_msp_mean_m": np.round(vp_lmn_vec_msp_mean[1] / 1e3, 3),
            "vp_lmn_vec_msp_mean_l": np.round(vp_lmn_vec_msp_mean[2] / 1e3, 3),
            "vp_lmn_vec_msp_median_n": np.round(vp_lmn_vec_msp_median[0] / 1e3, 3),
            "vp_lmn_vec_msp_median_m": np.round(vp_lmn_vec_msp_median[1] / 1e3, 3),
            "vp_lmn_vec_msp_median_l": np.round(vp_lmn_vec_msp_median[2] / 1e3, 3),
            "vp_lmn_vec_msh_mean_n": np.round(vp_lmn_vec_msh_mean[0] / 1e3, 3),
            "vp_lmn_vec_msh_mean_m": np.round(vp_lmn_vec_msh_mean[1] / 1e3, 3),
            "vp_lmn_vec_msh_mean_l": np.round(vp_lmn_vec_msh_mean[2] / 1e3, 3),
            "vp_lmn_vec_msh_median_n": np.round(vp_lmn_vec_msh_median[0] / 1e3, 3),
            "vp_lmn_vec_msh_median_m": np.round(vp_lmn_vec_msh_median[1] / 1e3, 3),
            "vp_lmn_vec_msh_median_l": np.round(vp_lmn_vec_msh_median[2] / 1e3, 3),
            "tp_para_msp_median": np.round(tp_para_msp_median, 3),
            "tp_para_msh_median": np.round(tp_para_msh_median, 3),
            "tp_para_msp_mean": np.round(tp_para_msp_mean, 3),
            "tp_para_msh_mean": np.round(tp_para_msh_mean, 3),
            "tp_perp_msp_median": np.round(tp_perp_msp_median, 3),
            "tp_perp_msh_median": np.round(tp_perp_msh_median, 3),
            "tp_perp_msp_mean": np.round(tp_perp_msp_mean, 3),
            "tp_perp_msh_mean": np.round(tp_perp_msh_mean, 3),
            "beta_msh_mean": np.round(beta_msh_mean, 3),
            "beta_msp_mean": np.round(beta_msp_mean, 3),
        }

        if x > -5 and x < 12 and r_yz < 12:
            # Check if the file exists, if not then create it
            os.makedirs(os.path.dirname(fname), exist_ok=True)
            if not os.path.isfile(fname):
                with open(fname, "w") as f:
                    f.write(var_list + "\n")
                    if verbose:
                        print(f"File {fname} created")
            # Check if the file already contains the data corresponding to the crossing time
            if os.path.isfile(fname):
                df_added_list = pd.read_csv(fname, sep=",", index_col=False)
                if not np.any(df_added_list["Date"].values == str(crossing_time)):
                    with open(fname, "a") as f:
                        for key in data_dict.keys():
                            f.write(f"{data_dict[key]},")
                        f.write("\n")
                    if verbose:
                        print(f"Data for all keys written to file {fname}")
                    f.close()

    fig = tplot_fnc(
        probe=probe,
        data_rate=data_rate,
        level=level,
        df_mms=df_mms,
        ind_range_msp=ind_range_msp,
        ind_range_msh=ind_range_msh,
        t_jet_center=t_jet_center,
        jet_detection=jet_detection,
        ind_crossing=ind_crossing,
        shear_val=int(angle_b_lmn_vec_msp_msh_median),
        date_obs=date_obs,
        return_plotly_fig=return_plotly_fig,
        dark_mode=dark_mode,
        figname=figname,
        save_pytplot=save_jet_stat,
    )

    return fig, jet_detection, data_dict

tplot_fnc(probe=3, data_rate='brst', level='l2', df_mms=None, ind_range_msp=None, ind_range_msh=None, t_jet_center=None, jet_detection=False, ind_crossing=None, shear_val=None, date_obs=None, return_plotly_fig=False, dark_mode=False, figname=None, save_pytplot=True)

Plot the data from the MMS spacecraft along with jet detection results

Parameters:

Name	Type	Description	Default
probe	int	The probe number. Default is 3	3
data_rate	str	The data rate. Default is 'brst'	'brst'
df_mms	DataFrame	The dataframe containing the MMS data	None
ind_range_msp	ndarray	The indices of the MSP	None
ind_range_msh	ndarray	The indices of the MSH	None
t_jet_center	datetime	The time of the jet center	None
jet_detection	bool	The status of the jet detection	False
ind_crossing	int	The index of the crossing date	None
shear_val	float	The shear value between the magnetosheath and the magnetosphere	None
date_obs	datetime	The observation date	None
return_plotly_fig	bool	Whether to return a Plotly figure instead of saving a tplot. Default is False.	False
dark_mode	bool	Whether to use dark mode for the Plotly figure. Default is False.	False
figname	str	The full filename to save the tplot figure to. If None, a default name is generated.	None

Returns:

Type	Description
None

Source code in src/rxn_location/jet_reversal_check_function.py

def tplot_fnc(
    probe=3,
    data_rate="brst",
    level="l2",
    df_mms=None,
    ind_range_msp=None,
    ind_range_msh=None,
    t_jet_center=None,
    jet_detection=False,
    ind_crossing=None,
    shear_val=None,
    date_obs=None,
    return_plotly_fig=False,
    dark_mode=False,
    figname=None,
    save_pytplot=True,
):
    """
    Plot the data from the MMS spacecraft along with jet detection results

    Parameters
    ----------
    probe : int
        The probe number. Default is 3
    data_rate : str
        The data rate. Default is 'brst'
    df_mms : pandas.DataFrame
        The dataframe containing the MMS data
    ind_range_msp : numpy.ndarray
        The indices of the MSP
    ind_range_msh : numpy.ndarray
        The indices of the MSH
    t_jet_center : datetime.datetime
        The time of the jet center
    jet_detection : bool
        The status of the jet detection
    ind_crossing : int
        The index of the crossing date
    shear_val : float
        The shear value between the magnetosheath and the magnetosphere
    date_obs : datetime.datetime
        The observation date
    return_plotly_fig : bool
        Whether to return a Plotly figure instead of saving a tplot. Default is False.
    dark_mode : bool
        Whether to use dark mode for the Plotly figure. Default is False.
    figname : str
        The full filename to save the tplot figure to. If None, a default name is generated.

    Returns
    -------
    None
    """
    # Set the fontstyle to Times New Roman
    font = {"family": "serif", "weight": "normal", "size": 12}
    plt.rc("font", **font)
    plt.rc("text", usetex=False)
    # Set tick parameters such that ticks are visible on all sides
    plt.rcParams["xtick.direction"] = "in"
    plt.rcParams["ytick.direction"] = "in"
    plt.rcParams["xtick.top"] = True
    plt.rcParams["ytick.right"] = True

    # Set the tick length and widths
    plt.rcParams["xtick.major.size"] = 5
    plt.rcParams["xtick.major.width"] = 1
    plt.rcParams["ytick.major.size"] = 5
    plt.rcParams["ytick.major.width"] = 1
    plt.rcParams["xtick.minor.size"] = 3
    plt.rcParams["xtick.minor.width"] = 1
    plt.rcParams["ytick.minor.size"] = 3
    plt.rcParams["ytick.minor.width"] = 1

    # Define the size of the figure
    plt.rcParams["figure.figsize"] = [10, 12]

    tplot_global_options = {
        # "show_all_axes": True,
        # "black_background": True,
        # "crosshair": True,
        # "vertical_spacing": 0,
        # "wsize": [1080, 1080],
    }
    for key in tplot_global_options:
        spd.tplot_options(key, tplot_global_options[key])

    keys_to_plot = [
        f"mms{probe}_dis_energyspectr_omni_{data_rate}",
        # f'mms{probe}_des_energyspectr_omni_{data_rate}',
        f"mms{probe}_fgm_b_lmn_srvy_{level}",
        f"mms{probe}_dis_numberdensity_{data_rate}",
        "Tp",
        f"mms{probe}_dis_bulkv_lmn_{data_rate}",
        "delta_v_vp_lmn_diff_l" if jet_detection else "vp_lmn_diff_l",
    ]

    ion_energy_spectr_dict_option = {
        "Colormap": "viridis_r",
        "ylog": True,
        "zlog": True,
        "ytitle": "$i^+$",
        "ysubtitle": "[eV]",
        "ztitle": "$\\rm{[keV/(cm^2\\,s\\,sr\\ keV)]}$",
        # 'title': f"Jet Detection for MMS{probe} {data_rate} data",
    }

    electron_energy_spectr_dict_option = {
        "Colormap": "Spectral_r",
        "ylog": True,
        "zlog": True,
        "ytitle": "$e^-$",
        "ysubtitle": "[eV]",
        "ztitle": "Counts",
    }

    number_density_dict_option = {
        "ytitle": "$n_{\\rm i}$",
        "ysubtitle": "[cm$^{-3}$]",
        "ylog": True,
    }

    tp_dict_option = {
        "ylog": True,
        "color": ["red", "blue"],
        "linestyle": "-",
        "legend_names": ["$\\parallel$", "$\\perp$"],
        "ytitle": "$T_{\\rm i}$",
        "ysubtitle": "[eV]",
    }

    b_dict_option = {
        "ytitle": "$B$",
        "ysubtitle": "[nT]",
        "color": ["blue", "green", "red"],
        "legend_names": ["L", "M", "N"],
        "linestyle": "-",
    }

    bulkv_dict_option = {
        "ytitle": "$v_{\\rm i}$",
        "ysubtitle": "[km/s]",
        "xtitle": "Time",
        "xsubtitle": "(HH:MM) [UTC]",
        "color": ["blue", "green", "red"],
        "linestyle": "-",
        "legend_names": ["L", "M", "N"],
    }

    print(f"\n Jet detection: {jet_detection} \n")

    # Define the limits of the delta v plot
    dv_vpl_min = np.nanmin(spd.get_data("vp_lmn_diff_l")[1:])
    dv_vpl_max = np.nanmax(spd.get_data("vp_lmn_diff_l")[1:])

    if jet_detection:
        dv_min_min = np.nanmin(spd.get_data("delta_v_min")[1:])
        dv_min_max = np.nanmax(spd.get_data("delta_v_min")[1:])
        dv_max_min = np.nanmin(spd.get_data("delta_v_max")[1:])
        dv_max_max = np.nanmax(spd.get_data("delta_v_max")[1:])

        dv_min = 1.1 * min(dv_min_min, dv_max_min, dv_vpl_min)
        dv_max = 1.1 * max(dv_min_max, dv_max_max, dv_vpl_max)

        delta_v_dict_option = {
            "color": ["blue", "green", "red"],
            "linestyle": ["-", "--", "--"],
            "yrange": [dv_min, dv_max],
            "ytitle": "$\\Delta v$",
            "ysubtitle": "[km/s]",
            "xtitle": "Time",
            "xsubtitle": "(HH:MM) [UTC]",
            "legend_names": [
                "$\\Delta \\rm v_{\\rm i,L}$",
                "$\\Delta \\rm v_{\\rm min}$",
                "$\\Delta \\rm v_{\\rm max}$",
            ],
        }
    else:
        dv_min = 1.1 * dv_vpl_min
        dv_max = 1.1 * dv_vpl_max

        delta_v_dict_option = {
            "color": ["blue"],
            "linestyle": ["-"],
            "yrange": [dv_min, dv_max],
            "ytitle": "$\\Delta v$",
            "ysubtitle": "[km/s]",
            "xtitle": "Time",
            "xsubtitle": "(HH:MM) [UTC]",
            "legend_names": ["$\\Delta \\rm v_{\\rm i,L}$"],
        }

    msp_time = df_mms.index[ind_range_msp[int(len(ind_range_msp) / 2)]]
    msh_time = df_mms.index[ind_range_msh[int(len(ind_range_msh) / 2)]]

    # Convert msp_time to UNIX time
    msp_time_unix = msp_time.timestamp()
    msh_time_unix = msh_time.timestamp()
    t_jet_center_unix = t_jet_center.timestamp()
    spd.timebar(msp_time_unix, databar=False, color="red", dash=True, thick=2)
    spd.timebar(msh_time_unix, databar=False, color="blue", dash=True, thick=2)
    spd.timebar(t_jet_center_unix, databar=False, color="green", dash=True, thick=1)

    spd.options(
        f"mms{probe}_dis_energyspectr_omni_{data_rate}",
        opt_dict=ion_energy_spectr_dict_option,
    )
    spd.options(
        f"mms{probe}_des_energyspectr_omni_{data_rate}",
        opt_dict=electron_energy_spectr_dict_option,
    )
    spd.options(
        f"mms{probe}_dis_numberdensity_{data_rate}", opt_dict=number_density_dict_option
    )

    spd.options("Tp", opt_dict=tp_dict_option)
    spd.options(f"mms{probe}_fgm_b_lmn_srvy_{level}", opt_dict=b_dict_option)
    spd.options(f"mms{probe}_dis_bulkv_lmn_{data_rate}", opt_dict=bulkv_dict_option)
    spd.options(
        "delta_v_vp_lmn_diff_l" if jet_detection else "vp_lmn_diff_l",
        opt_dict=delta_v_dict_option,
    )

    if figname is None:
        if jet_detection:
            folder_name = (
                f"figures/jet_reversal_checks/check_{date_obs}/{data_rate}/jet"
            )
            if not os.path.exists(folder_name):
                os.makedirs(folder_name)
        else:
            folder_name = (
                f"figures/jet_reversal_checks/check_{date_obs}/{data_rate}/no_jet"
            )
            if not os.path.exists(folder_name):
                os.makedirs(folder_name)

        figname = (
            f"{folder_name}/mms{probe}_{t_jet_center.strftime('%Y%m%d_%H%M')}_"
            + f"{str(ind_crossing).zfill(5)}_{shear_val}s_"
        )

    if save_pytplot:
        spd.tplot(keys_to_plot, save_png=figname, display=False)

    if return_plotly_fig:
        from rxn_location.plotly_utils import convert_tplot_to_plotly

        fig = convert_tplot_to_plotly(keys_to_plot, dark_mode=dark_mode)

        # Add timebars using go.Scatter vertical lines
        for unix_t, color, label in [
            (msp_time_unix, "red", "Magnetopause"),
            (t_jet_center_unix, "green", "Jet Center"),
            (msh_time_unix, "blue", "Magnetosheath"),
        ]:
            t_dt = datetime.datetime.fromtimestamp(unix_t, datetime.timezone.utc)
            fig.add_vline(
                x=t_dt, 
                line_dash="dash", 
                line_color=color, 
                line_width=2,
                annotation_text=label,
                annotation_position="top right",
                annotation_font_size=12,
                annotation_font_color=color
            )

        fig.update_layout(title=f"Jet Detection for MMS{probe} {data_rate} data")
        return fig
    else:
        return None

Data Management

rxn_location.master_jet_list

Master Jet List Management Module

Provides persistent storage and management of observed reconnection jet events. The master list is stored as a JSON file at ~/.rxn_location_master_jets.json, enabling deduplication (2-minute window), export, and quick re-run functionality.

add_jet(entries, data_dict, crossing_time, params, window_minutes=2)

Add a jet to the master list if no duplicate exists within the time window.

Parameters:

Name	Type	Description	Default
entries	list of dict	The current master list (modified in place if added).	required
data_dict	dict	The data dictionary returned by jet_reversal_check().	required
crossing_time	datetime	The user-input crossing time.	required
params	dict	The sidebar parameters used for this run (probe, dt, data_rate, etc.).	required
window_minutes	int or float	Deduplication window in minutes (default 2).	2

Returns:

Type	Description
tuple of (bool, dict or None)	(True, None) if the jet was added as new. (False, existing_entry) if a duplicate was found.

Source code in src/rxn_location/master_jet_list.py

def add_jet(entries, data_dict, crossing_time, params, window_minutes=2):
    """
    Add a jet to the master list if no duplicate exists within the time window.

    Parameters
    ----------
    entries : list of dict
        The current master list (modified in place if added).
    data_dict : dict
        The data dictionary returned by jet_reversal_check().
    crossing_time : datetime.datetime
        The user-input crossing time.
    params : dict
        The sidebar parameters used for this run (probe, dt, data_rate, etc.).
    window_minutes : int or float
        Deduplication window in minutes (default 2).

    Returns
    -------
    tuple of (bool, dict or None)
        (True, None) if the jet was added as new.
        (False, existing_entry) if a duplicate was found.
    """
    jet_time_str = str(data_dict.get("jet_time", crossing_time))

    existing = find_nearby_jet(entries, jet_time_str, window_minutes=window_minutes)
    if existing is not None:
        return False, existing

    # Build a serializable entry with rounded jet_time (nearest second)
    jet_time_rounded = _round_time_to_seconds(jet_time_str)
    crossing_rounded = _round_time_to_seconds(str(crossing_time))

    entry = {
        "jet_time": jet_time_rounded,
        "crossing_time": crossing_rounded,
        "added_at": datetime.datetime.now(pytz.utc).strftime("%Y-%m-%d %H:%M:%S"),
    }

    # Add sidebar parameters
    for key in [
        "mms_probe",
        "dt",
        "jet_len",
        "data_rate",
        "level",
        "coord_type",
        "time_clip",
        "t_delta",
        "max_attempts",
        "tsy_model",
        "recon_models",
        "omni_level",
        "m_p",
        "dr",
        "limits",
    ]:
        if key in params:
            entry[key] = params[key]

    # Add all data_dict fields with intelligent rounding
    for key, value in data_dict.items():
        rounded = _round_for_storage(key, value)
        entry[f"data_{key}"] = _to_serializable(rounded)

    entries.append(entry)
    return True, None

delete_jets(entries, indices)

Remove jets from the master list by index.

Parameters:

Name	Type	Description	Default
entries	list of dict	The current master list (modified in place).	required
indices	list of int	Sorted list of indices to remove.	required

Returns:

Type	Description
list of dict	The modified master list with entries removed.

Source code in src/rxn_location/master_jet_list.py

def delete_jets(entries, indices):
    """
    Remove jets from the master list by index.

    Parameters
    ----------
    entries : list of dict
        The current master list (modified in place).
    indices : list of int
        Sorted list of indices to remove.

    Returns
    -------
    list of dict
        The modified master list with entries removed.
    """
    # Remove in reverse order to preserve indices
    for idx in sorted(indices, reverse=True):
        if 0 <= idx < len(entries):
            entries.pop(idx)
    return entries

export_to_csv(entries, path='master_jets.csv')

Export the master jet list to a CSV file.

Parameters:

Name	Type	Description	Default
entries	list of dict	The master list entries.	required
path	str	Output file path.	'master_jets.csv'

Returns:

Type	Description
bytes	The CSV content as bytes for download.

Source code in src/rxn_location/master_jet_list.py

def export_to_csv(entries, path="master_jets.csv"):
    """
    Export the master jet list to a CSV file.

    Parameters
    ----------
    entries : list of dict
        The master list entries.
    path : str
        Output file path.

    Returns
    -------
    bytes
        The CSV content as bytes for download.
    """
    if not entries:
        return b""
    df = pd.DataFrame(entries)
    return df.to_csv(index=False).encode("utf-8")

export_to_json(entries)

Export the master jet list to JSON bytes.

Parameters:

Name	Type	Description	Default
entries	list of dict	The master list entries.	required

Returns:

Type	Description
bytes	The JSON content as bytes for download.

Source code in src/rxn_location/master_jet_list.py

def export_to_json(entries):
    """
    Export the master jet list to JSON bytes.

    Parameters
    ----------
    entries : list of dict
        The master list entries.

    Returns
    -------
    bytes
        The JSON content as bytes for download.
    """
    return json.dumps(entries, indent=2, default=str).encode("utf-8")

export_to_pickle(entries)

Export the master jet list to pickle bytes.

Parameters:

Name	Type	Description	Default
entries	list of dict	The master list entries.	required

Returns:

Type	Description
bytes	The pickled content as bytes for download.

Source code in src/rxn_location/master_jet_list.py

def export_to_pickle(entries):
    """
    Export the master jet list to pickle bytes.

    Parameters
    ----------
    entries : list of dict
        The master list entries.

    Returns
    -------
    bytes
        The pickled content as bytes for download.
    """
    return pickle.dumps(entries)

find_nearby_jet(entries, time, window_minutes=2)

Find a jet in the master list within a time window of the given time.

Parameters:

Name	Type	Description	Default
entries	list of dict	The master list entries.	required
time	datetime or str	The time to search around.	required
window_minutes	int or float	The half-window size in minutes (default 2).	2

Returns:

Type	Description
dict or None	The matching entry if found within the window, otherwise None.

Source code in src/rxn_location/master_jet_list.py

def find_nearby_jet(entries, time, window_minutes=2):
    """
    Find a jet in the master list within a time window of the given time.

    Parameters
    ----------
    entries : list of dict
        The master list entries.
    time : datetime.datetime or str
        The time to search around.
    window_minutes : int or float
        The half-window size in minutes (default 2).

    Returns
    -------
    dict or None
        The matching entry if found within the window, otherwise None.
    """
    target = _parse_time(time)
    window = datetime.timedelta(minutes=window_minutes)

    closest_entry = None
    closest_delta = None

    for entry in entries:
        try:
            entry_time = _parse_time(entry["jet_time"])
            delta = abs(entry_time - target)
            if delta <= window:
                if closest_delta is None or delta < closest_delta:
                    closest_entry = entry
                    closest_delta = delta
        except (KeyError, ValueError):
            continue

    return closest_entry

load_master_list(path=None)

Load the master jet list from disk.

Parameters:

Name	Type	Description	Default
path	str or Path	Path to the JSON file. Defaults to ~/.rxn_location_master_jets.json.	None

Returns:

Type	Description
list of dict	Each dict represents one observed jet event with its parameters.

Source code in src/rxn_location/master_jet_list.py

def load_master_list(path=None):
    """
    Load the master jet list from disk.

    Parameters
    ----------
    path : str or Path, optional
        Path to the JSON file. Defaults to ~/.rxn_location_master_jets.json.

    Returns
    -------
    list of dict
        Each dict represents one observed jet event with its parameters.
    """
    if path is None:
        path = DEFAULT_MASTER_LIST_PATH
    path = Path(path)
    if path.exists():
        try:
            with open(path, "r") as f:
                return json.load(f)
        except (json.JSONDecodeError, IOError):
            return []
    return []

master_list_to_stats_csv(entries, time_start=None, time_end=None)

Convert master list entries to a CSV string compatible with the statistics plots.

Filters entries to the given time range and reconstructs the wide-format CSV that generate_statistics_plots() expects.

Parameters:

Name	Type	Description	Default
entries	list of dict	The master list entries.	required
time_start	datetime	Start of time range filter.	None
time_end	datetime	End of time range filter.	None

Returns:

Type	Description
DataFrame or None	A DataFrame compatible with the statistics plotting functions, or None if no matching entries.

Source code in src/rxn_location/master_jet_list.py

def master_list_to_stats_csv(entries, time_start=None, time_end=None):
    """
    Convert master list entries to a CSV string compatible with the statistics plots.

    Filters entries to the given time range and reconstructs the wide-format CSV
    that generate_statistics_plots() expects.

    Parameters
    ----------
    entries : list of dict
        The master list entries.
    time_start : datetime.datetime, optional
        Start of time range filter.
    time_end : datetime.datetime, optional
        End of time range filter.

    Returns
    -------
    pandas.DataFrame or None
        A DataFrame compatible with the statistics plotting functions, or None if no matching entries.
    """
    if not entries:
        return None

    filtered = entries
    if time_start is not None:
        time_start = _parse_time(time_start)
        filtered = [
            e
            for e in filtered
            if _parse_time(e.get("jet_time", e.get("crossing_time"))) >= time_start
        ]
    if time_end is not None:
        time_end = _parse_time(time_end)
        filtered = [
            e
            for e in filtered
            if _parse_time(e.get("jet_time", e.get("crossing_time"))) <= time_end
        ]

    if not filtered:
        return None

    df = pd.DataFrame(filtered)

    # Rename columns to remove 'data_' prefix for compatibility with stats plots
    rename_map = {c: c[5:] for c in df.columns if c.startswith("data_")}
    df.rename(columns=rename_map, inplace=True)

    # Duplicate the dataframe 4 times to populate the standard 2x2 plot grid
    models = ["shear", "rx_en", "va_cs", "bisection"]

    # We map the standard model names to the exact keys output by the ridge finder.
    # We stripped "data_" earlier, so the column will be "r_rc_Shear", "r_rc_Reconnection Energy", etc.
    model_mapping = {
        "shear": "r_rc_Shear",
        "bisection": "r_rc_Bisection Field",
        "rx_en": "r_rc_Reconnection Energy",
        "va_cs": "r_rc_Exhaust Velocity",
    }

    df_list = []
    for model in models:
        df_copy = df.copy()
        df_copy["method_used"] = model

        # Check if the theoretical R_rc was computed and saved for this model
        target_col = model_mapping[model]

        if target_col in df.columns:
            df_copy["r_rc"] = df[target_col]
        else:
            df_copy["r_rc"] = np.nan

        df_list.append(df_copy)

    return pd.concat(df_list, ignore_index=True)

save_master_list(entries, path=None)

Save the master jet list to disk.

Parameters:

Name	Type	Description	Default
entries	list of dict	The master list entries to persist.	required
path	str or Path	Path to the JSON file. Defaults to ~/.rxn_location_master_jets.json.	None

Source code in src/rxn_location/master_jet_list.py

def save_master_list(entries, path=None):
    """
    Save the master jet list to disk.

    Parameters
    ----------
    entries : list of dict
        The master list entries to persist.
    path : str or Path, optional
        Path to the JSON file. Defaults to ~/.rxn_location_master_jets.json.
    """
    if path is None:
        path = DEFAULT_MASTER_LIST_PATH
    path = Path(path)
    try:
        with open(path, "w") as f:
            json.dump(entries, f, indent=2, default=str)
    except IOError:
        pass

rxn_location.presets

Presets Management Module

Allows users to save and load sidebar parameter configurations as named presets. Presets are stored in a JSON file at ~/.rxn_location_presets.json.

delete_preset(name, path=None)

Delete a named preset.

Parameters:

Name	Type	Description	Default
name	str	The preset name to remove.	required
path	str or Path	Path to the presets JSON file.	None

Returns:

Type	Description
bool	True if the preset was found and deleted, False otherwise.

Source code in src/rxn_location/presets.py

def delete_preset(name, path=None):
    """
    Delete a named preset.

    Parameters
    ----------
    name : str
        The preset name to remove.
    path : str or Path, optional
        Path to the presets JSON file.

    Returns
    -------
    bool
        True if the preset was found and deleted, False otherwise.
    """
    presets = load_presets(path)
    if name in presets:
        del presets[name]
        save_presets(presets, path)
        return True
    return False

get_current_params(session_state)

Extract the current sidebar parameters from Streamlit session state.

Parameters:

Name	Type	Description	Default
session_state	session_state	The current Streamlit session state.	required

Returns:

Type	Description
dict	A dictionary of all sidebar parameters.

Source code in src/rxn_location/presets.py

def get_current_params(session_state):
    """
    Extract the current sidebar parameters from Streamlit session state.

    Parameters
    ----------
    session_state : streamlit.session_state
        The current Streamlit session state.

    Returns
    -------
    dict
        A dictionary of all sidebar parameters.
    """
    return {
        "crossing_time_str": session_state.get(
            "crossing_time_str", "2015-09-02 16:45:00"
        ),
        "mms_probe": session_state.get("preset_mms_probe", 3),
        "dt": session_state.get("preset_dt", 300),
        "jet_len": session_state.get("preset_jet_len", 3),
        "data_rate": session_state.get("preset_data_rate", "brst"),
        "level": session_state.get("preset_level", "l2"),
        "coord_type": session_state.get("preset_coord_type", "lmn"),
        "time_clip": session_state.get("preset_time_clip", True),
        "t_delta": session_state.get("preset_t_delta", 10),
        "max_attempts": session_state.get("preset_max_attempts", 5),
        "tsy_model": session_state.get("preset_tsy_model", "t96"),
        "recon_models": session_state.get(
            "preset_recon_models",
            ["shear", "bisection", "reconnection energy", "exhaust velocity"],
        ),
        "omni_level": session_state.get("preset_omni_level", "hro_1min"),
        "m_p": session_state.get("preset_m_p", 1.0),
        "dr": session_state.get("preset_dr", 0.5),
        "limits": session_state.get("preset_limits", 20.0),
    }

list_preset_names(path=None)

List all available preset names.

Parameters:

Name	Type	Description	Default
path	str or Path	Path to the presets JSON file.	None

Returns:

Type	Description
list of str	Sorted list of preset names.

Source code in src/rxn_location/presets.py

def list_preset_names(path=None):
    """
    List all available preset names.

    Parameters
    ----------
    path : str or Path, optional
        Path to the presets JSON file.

    Returns
    -------
    list of str
        Sorted list of preset names.
    """
    presets = load_presets(path)
    return sorted(presets.keys())

load_presets(path=None)

Load all saved presets from disk.

Parameters:

Name	Type	Description	Default
path	str or Path	Path to the presets JSON file. Defaults to ~/.rxn_location_presets.json.	None

Returns:

Type	Description
dict	A dictionary mapping preset names to their parameter dictionaries.

Source code in src/rxn_location/presets.py

def load_presets(path=None):
    """
    Load all saved presets from disk.

    Parameters
    ----------
    path : str or Path, optional
        Path to the presets JSON file. Defaults to ~/.rxn_location_presets.json.

    Returns
    -------
    dict
        A dictionary mapping preset names to their parameter dictionaries.
    """
    if path is None:
        path = DEFAULT_PRESETS_PATH
    path = Path(path)
    if path.exists():
        try:
            with open(path, "r") as f:
                return json.load(f)
        except (json.JSONDecodeError, IOError):
            return {}
    return {}

save_preset(name, params, path=None)

Save or update a single named preset.

Parameters:

Name	Type	Description	Default
name	str	The preset name.	required
params	dict	The parameter dictionary to save.	required
path	str or Path	Path to the presets JSON file.	None

Source code in src/rxn_location/presets.py

def save_preset(name, params, path=None):
    """
    Save or update a single named preset.

    Parameters
    ----------
    name : str
        The preset name.
    params : dict
        The parameter dictionary to save.
    path : str or Path, optional
        Path to the presets JSON file.
    """
    presets = load_presets(path)
    presets[name] = params
    save_presets(presets, path)

save_presets(presets, path=None)

Save all presets to disk.

Parameters:

Name	Type	Description	Default
presets	dict	The full presets dictionary to persist.	required
path	str or Path	Path to the presets JSON file.	None

Source code in src/rxn_location/presets.py

def save_presets(presets, path=None):
    """
    Save all presets to disk.

    Parameters
    ----------
    presets : dict
        The full presets dictionary to persist.
    path : str or Path, optional
        Path to the presets JSON file.
    """
    if path is None:
        path = DEFAULT_PRESETS_PATH
    path = Path(path)
    try:
        with open(path, "w") as f:
            json.dump(presets, f, indent=2)
    except IOError:
        pass

rxn_location.logger

set_verbosity(level)

Configure global verbosity level: 3: Everything prints (including PySPEDAS). 2: Silence PySPEDAS/PyTplot, keep standard app logs. 1: Silence detailed app logs (only print major milestones/errors). 0: Silence completely.

Source code in src/rxn_location/logger.py

def set_verbosity(level):
    """
    Configure global verbosity level:
    3: Everything prints (including PySPEDAS).
    2: Silence PySPEDAS/PyTplot, keep standard app logs.
    1: Silence detailed app logs (only print major milestones/errors).
    0: Silence completely.
    """
    global _VERBOSITY
    _VERBOSITY = level

    if level <= 2:
        logging.getLogger("pyspedas").setLevel(logging.ERROR)
        logging.getLogger("pytplot").setLevel(logging.ERROR)
        logging.getLogger("urllib3").setLevel(logging.ERROR)
    else:
        logging.getLogger("pyspedas").setLevel(logging.INFO)
        logging.getLogger("pytplot").setLevel(logging.INFO)

vprint(level, message, color=None)

Print a message if the current verbosity is >= level. Available colors: red, green, yellow, blue, magenta, cyan, bold

Source code in src/rxn_location/logger.py

def vprint(level, message, color=None):
    """
    Print a message if the current verbosity is >= `level`.
    Available colors: red, green, yellow, blue, magenta, cyan, bold
    """
    if _VERBOSITY >= level:
        if color:
            colors = {
                "red": "\033[91m",
                "green": "\033[92m",
                "yellow": "\033[93m",
                "blue": "\033[94m",
                "magenta": "\033[95m",
                "cyan": "\033[96m",
                "bold": "\033[1m",
            }
            reset = "\033[0m"
            c = colors.get(color.lower(), "")
            print(f"{c}{message}{reset}")
            sys.stdout.flush()
        else:
            print(message)
            sys.stdout.flush()

Statistical Plotting

rxn_location.rc_stats_fncs

convert_wide_to_long(df)

Converts the new wide-format CSV into the legacy long-format for plotting compatibility.

This ensures that existing plotting code works seamlessly by melting the DataFrame so that each crossing is represented across multiple rows based on the theoretical models (shear, rx_en, etc.).

Parameters:

Name	Type	Description	Default
df	DataFrame	The input DataFrame, which can be in either wide or long format.	required

Returns:

Name	Type	Description
df_long	DataFrame	The output DataFrame converted to long format. If the input is already in long format, it is returned unmodified.

Source code in src/rxn_location/rc_stats_fncs.py

def convert_wide_to_long(df):
    """
    Converts the new wide-format CSV into the legacy long-format for plotting compatibility.

    This ensures that existing plotting code works seamlessly by melting the DataFrame so that
    each crossing is represented across multiple rows based on the theoretical models (shear, rx_en, etc.).

    Parameters
    ----------
    df : pandas.DataFrame
        The input DataFrame, which can be in either wide or long format.

    Returns
    -------
    df_long : pandas.DataFrame
        The output DataFrame converted to long format. If the input is already in long format,
        it is returned unmodified.
    """
    rc_cols = [col for col in df.columns if col.startswith("r_rc_")]

    if len(rc_cols) > 0:
        # It's a wide format dataframe
        id_vars = [
            col
            for col in df.columns
            if not col.startswith("r_rc_") and col != "r_rc" and col != "method_used"
        ]

        # Drop 'r_rc' if it already exists, as melt will create it and crash otherwise
        if "r_rc" in df.columns:
            df = df.drop(columns=["r_rc"])

        # Melt it
        df_long = pd.melt(
            df,
            id_vars=id_vars,
            value_vars=rc_cols,
            var_name="method_used",
            value_name="r_rc",
        )

        # Clean up the method_used column
        df_long["method_used"] = df_long["method_used"].str.replace("r_rc_", "")

        # Map readable names back to internal abbreviations for backward compatibility
        name_map = {
            "Shear": "shear",
            "Bisection Field": "bisection",
            "Reconnection Energy": "rx_en",
            "Exhaust Velocity": "va_cs",
        }
        df_long["method_used"] = df_long["method_used"].replace(name_map)

        return df_long
    else:
        # It's already in the legacy long format or doesn't have r_rc columns
        return df

plot_hist(file_name, fig_size=(6, 6), dark_mode=True, bins=8, fig_folder='figures', fig_name='new', fig_format='pdf', histtype='step', linewidth=1, cut_type='jet', r_lim=[0, 15], density=False, return_fig=False)

Generate histograms of reconnection location parameters.

Parameters:

Name	Type	Description	Default
file_name	str	Path to the CSV file containing the data.	required
fig_size	tuple	Figure size as (width, height). Default is (6, 6).	(6, 6)
dark_mode	bool	If True, applies a dark mode theme to the plot. Default is True.	True
bins	int	Number of bins for the histogram. Default is 8.	8
fig_folder	str	Directory where the figure will be saved. Default is "figures".	'figures'
fig_name	str	Base name of the saved figure. Default is "new".	'new'
fig_format	str	Format of the saved figure (e.g., 'pdf', 'png'). Default is "pdf".	'pdf'
histtype	str	Type of histogram to draw (e.g., 'step', 'bar'). Default is "step".	'step'
linewidth	int	Width of the histogram lines. Default is 1.	1
cut_type	str	Data filtering criteria (e.g., 'jet', 'bz_neg', 'bz_pos'). Default is "jet".	'jet'
r_lim	list of float	Limits for the x-axis representing the r_rc parameter. Default is [0, 15].	[0, 15]
density	bool	If True, plot a probability density rather than counts. Default is False.	False
return_fig	bool	If True, returns the matplotlib Figure object along with the dataframes. Default is False.	False

Returns:

Name	Type	Description
fig	Figure(optional)	The generated figure object, only returned if return_fig is True.
df_shear	DataFrame	Data frame containing data filtered by the shear model.
df_rx_en	DataFrame	Data frame containing data filtered by the reconnection energy model.
df_va_cs	DataFrame	Data frame containing data filtered by the exhaust velocity model.
df_bisec	DataFrame	Data frame containing data filtered by the bisection model.

Source code in src/rxn_location/rc_stats_fncs.py

def plot_hist(
    file_name,
    fig_size=(6, 6),
    dark_mode=True,
    bins=8,
    fig_folder="figures",
    fig_name="new",
    fig_format="pdf",
    histtype="step",
    linewidth=1,
    cut_type="jet",
    r_lim=[0, 15],
    density=False,
    return_fig=False,
):
    """
    Generate histograms of reconnection location parameters.

    Parameters
    ----------
    file_name : str
        Path to the CSV file containing the data.
    fig_size : tuple, optional
        Figure size as (width, height). Default is (6, 6).
    dark_mode : bool, optional
        If True, applies a dark mode theme to the plot. Default is True.
    bins : int, optional
        Number of bins for the histogram. Default is 8.
    fig_folder : str, optional
        Directory where the figure will be saved. Default is "figures".
    fig_name : str, optional
        Base name of the saved figure. Default is "new".
    fig_format : str, optional
        Format of the saved figure (e.g., 'pdf', 'png'). Default is "pdf".
    histtype : str, optional
        Type of histogram to draw (e.g., 'step', 'bar'). Default is "step".
    linewidth : int, optional
        Width of the histogram lines. Default is 1.
    cut_type : str, optional
        Data filtering criteria (e.g., 'jet', 'bz_neg', 'bz_pos'). Default is "jet".
    r_lim : list of float, optional
        Limits for the x-axis representing the r_rc parameter. Default is [0, 15].
    density : bool, optional
        If True, plot a probability density rather than counts. Default is False.
    return_fig : bool, optional
        If True, returns the matplotlib Figure object along with the dataframes. Default is False.

    Returns
    -------
    fig : matplotlib.figure.Figure (optional)
        The generated figure object, only returned if `return_fig` is True.
    df_shear : pandas.DataFrame
        Data frame containing data filtered by the shear model.
    df_rx_en : pandas.DataFrame
        Data frame containing data filtered by the reconnection energy model.
    df_va_cs : pandas.DataFrame
        Data frame containing data filtered by the exhaust velocity model.
    df_bisec : pandas.DataFrame
        Data frame containing data filtered by the bisection model.
    """
    df = pd.read_csv(file_name, index_col=False)
    df = convert_wide_to_long(df)

    # Set date_from as index
    df = df.set_index("date_from")

    df_shear = df[df.method_used == "shear"].copy()
    df_rx_en = df[df.method_used == "rx_en"].copy()
    df_va_cs = df[df.method_used == "va_cs"].copy()
    df_bisec = df[df.method_used == "bisection"].copy()

    def safe_cone_angle(df_n):
        """
        Safely calculates the solar wind cone angle from the IMF vector.

        The cone angle is defined as the angle between the Sun-Earth line (X-axis) 
        and the IMF vector. It safely handles invalid or zero magnetic field magnitudes.

        Parameters
        ----------
        df_n : pandas.DataFrame
            DataFrame containing the IMF components `b_imf_x`, `b_imf_y`, and `b_imf_z`.

        Returns
        -------
        float
            The cone angle in degrees. Returns NaN if the inputs are invalid.
        """
        if (
            all(c in df_n.columns for c in ["b_imf_x", "b_imf_y", "b_imf_z"])
            and len(df_n) > 0
        ):
            b_mag = np.sqrt(df_n.b_imf_x**2 + df_n.b_imf_y**2 + df_n.b_imf_z**2)
            cosang = (df_n.b_imf_x / b_mag).where(b_mag > 0).clip(-1, 1)
            return np.degrees(np.arccos(cosang))
        return pd.Series(dtype=float)

    cone_angle_shear = safe_cone_angle(df_shear)
    cone_angle_rx_en = safe_cone_angle(df_rx_en)
    cone_angle_va_cs = safe_cone_angle(df_va_cs)
    cone_angle_bisec = safe_cone_angle(df_bisec)

    df_shear["cone_angle"] = cone_angle_shear
    df_rx_en["cone_angle"] = cone_angle_rx_en
    df_va_cs["cone_angle"] = cone_angle_va_cs
    df_bisec["cone_angle"] = cone_angle_bisec

    if cut_type == "bz_neg":
        df_shear = df_shear[df_shear["b_imf_z"] < 0]
        df_rx_en = df_rx_en[df_rx_en["b_imf_z"] < 0]
        df_va_cs = df_va_cs[df_va_cs["b_imf_z"] < 0]
        df_bisec = df_bisec[df_bisec["b_imf_z"] < 0]

    if cut_type == "bz_pos":
        df_shear = df_shear[df_shear["b_imf_z"] > 0]
        df_rx_en = df_rx_en[df_rx_en["b_imf_z"] > 0]
        df_va_cs = df_va_cs[df_va_cs["b_imf_z"] > 0]
        df_bisec = df_bisec[df_bisec["b_imf_z"] > 0]

    if cut_type == "cone_angle":
        # Select all data where cone angle is greater than 36.87 degrees
        df_shear = df_shear[(df_shear.cone_angle <= 36.87)]
        df_rx_en = df_rx_en[(df_rx_en.cone_angle <= 36.87)]
        df_va_cs = df_va_cs[(df_va_cs.cone_angle <= 36.87)]
        df_bisec = df_bisec[(df_bisec.cone_angle <= 36.87)]

    if cut_type == "cone_and_bz_neg":
        df_shear = df_shear[
            (df_shear.cone_angle >= 36.87)
            & (df_shear["b_imf_z"] < 0)
            & (df_shear.cone_angle <= 150)
        ]
        df_rx_en = df_rx_en[
            (df_rx_en.cone_angle >= 36.87)
            & (df_rx_en["b_imf_z"] < 0)
            & (df_rx_en.cone_angle <= 150)
        ]
        df_va_cs = df_va_cs[
            (df_va_cs.cone_angle >= 36.87)
            & (df_va_cs["b_imf_z"] < 0)
            & (df_va_cs.cone_angle <= 150)
        ]
        df_bisec = df_bisec[
            (df_bisec.cone_angle >= 36.87)
            & (df_bisec["b_imf_z"] < 0)
            & (df_bisec.cone_angle <= 150)
        ]

    if dark_mode:
        plt.style.use("dark_background")
        # tick_color = 'w'  # color of the tick lines
        mtick_color = "w"  # color of the minor tick lines
        label_color = "w"  # color of the tick labels
        # clabel_color = 'w'  # color of the colorbar label
    else:
        plt.style.use("default")
        # tick_color = 'k'  # color of the tick lines
        mtick_color = "k"  # color of the minor tick lines
        label_color = "k"  # color of the tick labels
        # clabel_color = 'k'  # color of the colorbar label

    # Set the fontstyle to Times New Roman
    font = {"family": "serif", "weight": "normal", "size": 10}
    plt.rc("font", **font)
    # plt.rc('text', usetex=True)

    if not return_fig:
        plt.close("all")

    r_lim_val = 12
    shear_r_rc_mean = df_shear[df_shear.r_rc < r_lim_val]["r_rc"].mean()
    rx_en_r_rc_mean = df_rx_en[df_rx_en.r_rc < r_lim_val]["r_rc"].mean()
    va_cs_r_rc_mean = df_va_cs[df_va_cs.r_rc < r_lim_val]["r_rc"].mean()
    bisec_r_rc_mean = df_bisec[df_bisec.r_rc < r_lim_val]["r_rc"].mean()
    y_label = "Density" if density else "Counts"
    facecolor = "k" if dark_mode else "w"
    edgecolor = "w" if dark_mode else "k"
    fig = plt.figure(
        num=None,
        figsize=fig_size,
        dpi=200,
        facecolor=facecolor,
        edgecolor=edgecolor,
    )
    fig.subplots_adjust(
        left=0.01, right=0.99, top=0.99, bottom=0.01, wspace=0.0, hspace=0.0
    )
    gs = gridspec.GridSpec(2, 2, width_ratios=[1, 1])

    # Plot the histogram of the shear data
    axs1 = plt.subplot(gs[0, 0])
    axs1.hist(
        df_shear.r_rc,
        bins=bins,
        range=(0, 15),
        color="#1f77b4",
        alpha=0.5,
        density=density,
        histtype=histtype,
        linewidth=linewidth,
    )
    # Plot the median of the shear data and add atext to the line
    axs1.axvline(shear_r_rc_mean, color="#1f77b4", linestyle="--", linewidth=2)
    axs1.text(
        shear_r_rc_mean + 0.2,
        axs1.get_ylim()[1] * 0.2,
        "$R_{{\\rm{{rc}}}}$ = {:.2f}$R_{{\\rm E}}$".format(shear_r_rc_mean),
        fontsize=1.1 * t_label_size,
        color=label_color,
    )
    axs1.text(
        0.05,
        0.9,
        "(a)",
        fontsize=1.1 * t_label_size,
        color=label_color,
        transform=axs1.transAxes,
    )
    axs1.set_xlim(r_lim[0], r_lim[1])
    axs1.set_xscale("linear")
    # axs1.set_xlabel(r'$r_{rc}$', fontsize=label_size, color=label_color, labelpad=label_pad)
    axs1.set_ylabel(y_label, fontsize=label_size, color=label_color, labelpad=label_pad)

    # Plot the histogram of the rx_en data
    axs2 = plt.subplot(gs[0, 1])
    axs2.hist(
        df_rx_en.r_rc,
        bins=bins,
        range=(0, 15),
        color="#ff7f0e",
        alpha=0.5,
        density=density,
        histtype=histtype,
        linewidth=linewidth,
    )
    # Plot the median of the rx_en data and add atext to the line
    axs2.axvline(rx_en_r_rc_mean, color="#ff7f0e", linestyle="--", linewidth=2)
    axs2.text(
        rx_en_r_rc_mean + 0.2,
        axs2.get_ylim()[1] * 0.2,
        "$R_{{\\rm{{rc}}}}$ = {:.2f}$R_{{\\rm E}}$".format(rx_en_r_rc_mean),
        fontsize=1.1 * t_label_size,
        color=label_color,
    )
    axs2.text(
        0.05,
        0.9,
        "(b)",
        fontsize=1.1 * t_label_size,
        color=label_color,
        transform=axs2.transAxes,
    )
    axs2.set_xlim(r_lim[0], r_lim[1])
    axs2.set_xscale("linear")
    # axs2.set_xlabel(r'$r_{rc}$', fontsize=label_size, color=label_color, labelpad=label_pad)
    axs2.set_ylabel(y_label, fontsize=label_size, color=label_color, labelpad=label_pad)
    axs2.yaxis.set_label_position("right")

    # Plot the histogram of the va_cs data
    axs3 = plt.subplot(gs[1, 0])
    axs3.hist(
        df_va_cs.r_rc,
        bins=bins,
        range=(0, 15),
        color="#2ca02c",
        alpha=0.5,
        density=density,
        histtype=histtype,
        linewidth=linewidth,
    )
    # Plot the median of the va_cs data and add atext to the line
    axs3.axvline(va_cs_r_rc_mean, color="#2ca02c", linestyle="--", linewidth=2)
    axs3.text(
        va_cs_r_rc_mean + 0.2,
        axs3.get_ylim()[1] * 0.2,
        "$R_{{\\rm{{rc}}}}$ = {:.2f}$R_{{\\rm E}}$".format(va_cs_r_rc_mean),
        fontsize=1.1 * t_label_size,
        color=label_color,
    )
    axs3.text(
        0.05,
        0.9,
        "(c)",
        fontsize=1.1 * t_label_size,
        color=label_color,
        transform=axs3.transAxes,
    )
    axs3.set_xlim(r_lim[0], r_lim[1])
    axs3.set_xscale("linear")
    axs3.set_xlabel(
        r"$R_{\rm {rc}} [R_{{\rm E}}]$",
        fontsize=label_size,
        color=label_color,
        labelpad=label_pad,
    )
    axs3.set_ylabel(y_label, fontsize=label_size, color=label_color, labelpad=label_pad)

    # Plot the histogram of the bisection data
    axs4 = plt.subplot(gs[1, 1])
    axs4.hist(
        df_bisec.r_rc,
        bins=bins,
        range=(0, 15),
        color="#d62728",
        alpha=0.5,
        density=density,
        histtype=histtype,
        linewidth=linewidth,
    )
    # Plot the median of the bisection data and add atext to the line
    axs4.axvline(bisec_r_rc_mean, color="#d62728", linestyle="--", linewidth=2)
    axs4.text(
        bisec_r_rc_mean + 0.2,
        axs4.get_ylim()[1] * 0.2,
        "$R_{{\\rm{{rc}}}}$ = {:.2f}$R_{{\\rm E}}$".format(bisec_r_rc_mean),
        fontsize=1.1 * t_label_size,
        color=label_color,
    )
    axs4.text(
        0.05,
        0.9,
        "(d)",
        fontsize=1.1 * t_label_size,
        color=label_color,
        transform=axs4.transAxes,
    )
    axs4.set_xlim(r_lim[0], r_lim[1])
    axs4.set_xscale("linear")
    axs4.set_xlabel(
        r"$R_{\rm {rc}} [R_{{\rm E}}]$",
        fontsize=label_size,
        color=label_color,
        labelpad=label_pad,
    )
    axs4.set_ylabel(y_label, fontsize=label_size, color=label_color, labelpad=label_pad)
    axs4.yaxis.set_label_position("right")

    # Set the tick parameters
    axs1.tick_params(
        axis="both",
        direction="in",
        which="major",
        left=True,
        right=True,
        top=True,
        bottom=True,
        labelleft=True,
        labelright=False,
        labeltop=False,
        labelbottom=False,
        labelsize=t_label_size,
        length=tick_len,
        width=tick_width,
        labelcolor=label_color,
    )

    axs2.tick_params(
        axis="both",
        direction="in",
        which="major",
        left=True,
        right=True,
        top=True,
        bottom=True,
        labelleft=False,
        labelright=True,
        labeltop=False,
        labelbottom=False,
        labelsize=t_label_size,
        length=tick_len,
        width=tick_width,
        labelcolor=label_color,
    )

    axs3.tick_params(
        axis="both",
        direction="in",
        which="major",
        left=True,
        right=True,
        top=True,
        bottom=True,
        labelleft=True,
        labelright=False,
        labeltop=False,
        labelbottom=True,
        labelsize=t_label_size,
        length=tick_len,
        width=tick_width,
        labelcolor=label_color,
    )

    axs4.tick_params(
        axis="both",
        direction="in",
        which="major",
        left=True,
        right=True,
        top=True,
        bottom=True,
        labelleft=False,
        labelright=True,
        labeltop=False,
        labelbottom=True,
        labelsize=t_label_size,
        length=tick_len,
        width=tick_width,
        labelcolor=label_color,
    )

    axs1.text(
        0.95,
        0.95,
        "Shear",
        ha="right",
        va="top",
        transform=axs1.transAxes,
        fontsize=c_label_size,
        color=label_color,
    )
    axs2.text(
        0.95,
        0.95,
        "Reconnection\n Energy",
        ha="right",
        va="top",
        transform=axs2.transAxes,
        fontsize=c_label_size,
        color=label_color,
    )
    axs3.text(
        0.95,
        0.95,
        "Exhaust\n Velocity",
        ha="right",
        va="top",
        transform=axs3.transAxes,
        fontsize=c_label_size,
        color=label_color,
    )
    axs4.text(
        0.95,
        0.95,
        "Bisection",
        ha="right",
        va="top",
        transform=axs4.transAxes,
        fontsize=c_label_size,
        color=label_color,
    )

    # Show minor ticks
    axs1.minorticks_on()
    axs1.tick_params(
        axis="both",
        which="minor",
        direction="in",
        length=mtick_len,
        left=True,
        right=True,
        top=True,
        bottom=True,
        color=mtick_color,
        width=mtick_width,
    )
    axs2.minorticks_on()
    axs2.tick_params(
        axis="both",
        which="minor",
        direction="in",
        length=mtick_len,
        left=True,
        right=True,
        top=True,
        bottom=True,
        color=mtick_color,
        width=mtick_width,
    )
    axs3.minorticks_on()
    axs3.tick_params(
        axis="both",
        which="minor",
        direction="in",
        length=mtick_len,
        left=True,
        right=True,
        top=True,
        bottom=True,
        color=mtick_color,
        width=mtick_width,
    )
    axs4.minorticks_on()
    axs4.tick_params(
        axis="both",
        which="minor",
        direction="in",
        length=mtick_len,
        left=True,
        right=True,
        top=True,
        bottom=True,
        color=mtick_color,
        width=mtick_width,
    )

    # Set the number of ticks on the x-and y-axis
    axs1.xaxis.set_major_locator(MaxNLocator(nbins=5, prune="lower"))
    axs1.yaxis.set_major_locator(MaxNLocator(nbins=5, prune="lower"))
    axs2.xaxis.set_major_locator(MaxNLocator(nbins=5, prune="lower"))
    axs2.yaxis.set_major_locator(MaxNLocator(nbins=5, prune="lower"))
    axs3.xaxis.set_major_locator(MaxNLocator(nbins=5, prune="lower"))
    axs3.yaxis.set_major_locator(MaxNLocator(nbins=5, prune="lower"))
    axs4.xaxis.set_major_locator(MaxNLocator(nbins=5, prune="lower"))
    axs4.yaxis.set_major_locator(MaxNLocator(nbins=5, prune="lower"))

    # Setting the tickmarks labels in such a way that they don't overlap
    plt.setp(axs1.get_xticklabels(), rotation=0, ha="right", va="top", visible=True)
    plt.setp(axs1.get_yticklabels(), rotation=0, va="center", visible=True)

    if dark_mode:
        transparent = False
    else:
        transparent = True

    # Save the figure
    if not os.path.exists(fig_folder):
        os.makedirs(fig_folder)
    fig_name = f"{fig_folder}/{fig_name}_{cut_type}.{fig_format}"
    plt.savefig(
        fig_name,
        bbox_inches="tight",
        pad_inches=0.05,
        dpi=200,
        transparent=transparent,
        format=fig_format,
    )
    if not return_fig:
        plt.close()
    print(f"Figure saved as {fig_name} in {fig_format} format in {fig_folder}")

    if return_fig:
        return fig, df_shear, df_rx_en, df_va_cs, df_bisec
    return df_shear, df_rx_en, df_va_cs, df_bisec

rxn_location.app_seaborn_plots

Seaborn-based statistical plots for the Streamlit GUI.

Adapted from seaborn_plots_fncs.py and SeabornFig2Grid.py for in-app rendering. These functions are designed to return matplotlib Figure objects for use with st.pyplot(), without saving to disk.

generate_seaborn_jointplots(df_full, x_key='r_rc', y_key='b_imf_z', x_label='Reconnection Distance [R_E]', y_label='IMF Bz [nT]', x_lim=None, y_lim=None, dark_mode=False, marker_size_var='None', nbins=(40, 40), figsize=(16, 16))

Generate a 2×2 grid of seaborn joint-plots (one per reconnection model) from a statistics DataFrame.

Parameters:

Name	Type	Description	Default
df_full	DataFrame	Full reconnection statistics dataframe containing a method_used column used to split data into the four models.	required
x_key	str	Column names for the x- and y-axis variables.	'r_rc'
y_key	str	Column names for the x- and y-axis variables.	'r_rc'
x_label	str	Human-readable axis labels.	'Reconnection Distance [R_E]'
y_label	str	Human-readable axis labels.	'Reconnection Distance [R_E]'
x_lim	tuple of float	Axis limits. Auto-calculated when None.	None
y_lim	tuple of float	Axis limits. Auto-calculated when None.	None
dark_mode	bool	Apply dark background styling.	False
marker_size_var	str	Scale marker size by the specified column name.	'None'
nbins	tuple of int	Number of bins for (x, y) marginal histograms.	(40, 40)
figsize	tuple of int	Overall figure size.	(16, 16)

Returns:

Name	Type	Description
fig	Figure	The composited 2×2 figure ready for st.pyplot(fig).

Source code in src/rxn_location/app_seaborn_plots.py

def generate_seaborn_jointplots(
    df_full,
    x_key="r_rc",
    y_key="b_imf_z",
    x_label="Reconnection Distance [R_E]",
    y_label="IMF Bz [nT]",
    x_lim=None,
    y_lim=None,
    dark_mode=False,
    marker_size_var="None",
    nbins=(40, 40),
    figsize=(16, 16),
):
    """
    Generate a 2×2 grid of seaborn joint-plots (one per reconnection model)
    from a statistics DataFrame.

    Parameters
    ----------
    df_full : pandas.DataFrame
        Full reconnection statistics dataframe containing a ``method_used``
        column used to split data into the four models.
    x_key, y_key : str
        Column names for the x- and y-axis variables.
    x_label, y_label : str
        Human-readable axis labels.
    x_lim, y_lim : tuple of float, optional
        Axis limits.  Auto-calculated when *None*.
    dark_mode : bool
        Apply dark background styling.
    marker_size_var : str
        Scale marker size by the specified column name.
    nbins : tuple of int
        Number of bins for (x, y) marginal histograms.
    figsize : tuple of int
        Overall figure size.

    Returns
    -------
    fig : matplotlib.figure.Figure
        The composited 2×2 figure ready for ``st.pyplot(fig)``.
    """

    # --- Set style ---
    if dark_mode:
        plt.style.use("dark_background")
        plt.rcParams.update({"axes.facecolor": "black", "figure.facecolor": "black"})
    else:
        plt.style.use("default")

    font = {"family": "serif", "weight": "normal", "size": 10}
    plt.rc("font", **font)
    # Do NOT enable usetex — Streamlit servers rarely have LaTeX installed
    plt.rc("text", usetex=False)

    # --- Split by model ---
    model_map = {
        "shear": "Shear",
        "rx_en": "Reconnection-Energy",
        "va_cs": "Exhaust-Velocity",
        "bisection": "Bisection",
    }
    # Normalise the method_used values coming from the CSV
    method_col = df_full["method_used"].str.strip().str.lower()

    # Build sub-dataframes in canonical order
    canonical_order = ["shear", "rx_en", "va_cs", "bisection"]
    # Create a mapping from possible CSV values → canonical key
    alias_map = {
        "shear": "shear",
        "reconnection energy": "rx_en",
        "reconnection-energy": "rx_en",
        "rx_en": "rx_en",
        "exhaust velocity": "va_cs",
        "exhaust-velocity": "va_cs",
        "va_cs": "va_cs",
        "bisection": "bisection",
        "bisection field": "bisection",
    }

    df_dict = {}
    for canon in canonical_order:
        matching_aliases = [k for k, v in alias_map.items() if v == canon]
        mask = method_col.isin(matching_aliases)
        df_dict[canon] = df_full[mask].copy()

    df_list = [df_dict[c] for c in canonical_order]
    data_type = [model_map[c] for c in canonical_order]
    color_list = ["#1f77b4", "#ff7f0e", "#2ca02c", "#d62728"]

    # --- Compute derived columns if needed ---
    # Cone angle
    if "cone_angle" not in df_full.columns:
        # Use shear df as reference for b_imf components
        ref = df_list[0] if len(df_list[0]) > 0 else df_full
        if all(c in ref.columns for c in ("b_imf_x", "b_imf_y", "b_imf_z")):
            bmag = np.sqrt(
                ref["b_imf_x"] ** 2 + ref["b_imf_y"] ** 2 + ref["b_imf_z"] ** 2
            )
            cone = np.arccos(ref["b_imf_x"] / bmag) * 180 / np.pi
            for dfn in df_list:
                if len(dfn) > 0:
                    dfn["cone_angle"] = cone.reindex(dfn.index).values

    # By/|B| ratio
    if "bb" not in df_full.columns:
        ref = df_list[0] if len(df_list[0]) > 0 else df_full
        if all(c in ref.columns for c in ("b_imf_x", "b_imf_y", "b_imf_z")):
            bmag = np.sqrt(
                ref["b_imf_x"] ** 2 + ref["b_imf_y"] ** 2 + ref["b_imf_z"] ** 2
            )
            bb_ratio = ref["b_imf_y"] / bmag
            for dfn in df_list:
                if len(dfn) > 0:
                    dfn["bb"] = bb_ratio.reindex(dfn.index).values

    # Helper for model-specific Rc keys
    model_map_rc = {
        "shear": "r_rc_Shear",
        "bisection": "r_rc_Bisection Field",
        "rx_en": "r_rc_Reconnection Energy",
        "va_cs": "r_rc_Exhaust Velocity",
    }

    def get_actual_key(base_key, canon_model, test_df):
        """
        Helper to map user-friendly label names back to their raw master jet list keys.
        """
        if base_key == "r_rc":
            target = model_map_rc.get(canon_model)
            if target:
                if f"data_{target}" in test_df.columns:
                    return f"data_{target}"
                if target in test_df.columns:
                    return target
        return base_key

    # --- Auto-calculate limits ---
    if x_lim is None:
        all_x = []
        for c, df in df_dict.items():
            k = get_actual_key(x_key, c, df_full)
            if len(df) > 0 and k in df.columns:
                all_x.append(df[k].dropna().values)
        if all_x:
            combined = np.concatenate(all_x)
            x_lim = (float(np.nanmin(combined)), float(np.nanmax(combined)))
        else:
            x_lim = (0, 1)
        if x_lim[0] == x_lim[1]:
            x_lim = (x_lim[0] - 1, x_lim[1] + 1)

    if y_lim is None:
        all_y = []
        for c, df in df_dict.items():
            k = get_actual_key(y_key, c, df_full)
            if len(df) > 0 and k in df.columns:
                all_y.append(df[k].dropna().values)
        if all_y:
            combined = np.concatenate(all_y)
            y_lim = (float(np.nanmin(combined)), float(np.nanmax(combined)))
        else:
            y_lim = (0, 1)
        if y_lim[0] == y_lim[1]:
            y_lim = (y_lim[0] - 1, y_lim[1] + 1)

    # --- Build individual JointGrid panels ---
    panels = []
    bins = [
        np.linspace(x_lim[0], x_lim[1], nbins[0]),
        np.linspace(y_lim[0], y_lim[1], nbins[1]),
    ]

    for i, (canon_model, df) in enumerate(df_dict.items()):
        act_x = get_actual_key(x_key, canon_model, df_full)
        act_y = get_actual_key(y_key, canon_model, df_full)
        act_marker = get_actual_key(marker_size_var, canon_model, df_full)

        if len(df) == 0 or act_x not in df.columns or act_y not in df.columns:
            panels.append(None)
            continue

        spearman = (
            df[act_y].corr(df[act_x], method="spearman")
            if act_x in df.columns and marker_size_var == "None"
            else None
        )

        msz = 100
        if marker_size_var != "None" and act_marker in df.columns:
            if marker_size_var == "r_rc":
                msz = 3 * df[act_marker].values ** 2
            else:
                vals = df[act_marker].values
                v_min = np.nanmin(vals)
                v_max = np.nanmax(vals)
                if v_max > v_min:
                    norm_vals = (vals - v_min) / (v_max - v_min)
                    msz = 20 + norm_vals * 280
                else:
                    msz = 100

        jg = _kde_plot_panel(
            df=df,
            x=act_x,
            y=act_y,
            x_label=x_label,
            y_label=y_label,
            data_type=data_type[i],
            xlim=x_lim,
            ylim=y_lim,
            color=color_list[i],
            marker_size=msz,
            spearman=spearman,
            bins=bins,
            dark_mode=dark_mode,
            marker_size_var=marker_size_var,
        )
        panels.append(jg)

    # --- Composite into a single 2×2 figure ---
    fig = plt.figure(figsize=figsize)
    gs = gridspec.GridSpec(2, 2)

    for idx, panel in enumerate(panels):
        if panel is not None:
            _SeabornFig2Grid(panel, fig, gs[idx])

    gs.tight_layout(fig, rect=[0.02, 0.02, 0.98, 0.98])
    gs.update(top=0.95, bottom=0.08, left=0.1, right=0.95, hspace=0.1, wspace=0.3)

    # Close any stray panel figures to prevent memory leaks
    plt.close("all")
    # Re-register our composited figure so it stays alive
    # (plt.close("all") would have closed it too)
    # Instead, we just return it — the caller will render with st.pyplot(fig)

    return fig

rxn_location.app_stats_plots_interactive

Interactive Statistics Plots Module

Generates interactive Plotly figures for the statistics mode of the GUI. This is the interactive counterpart to app_stats_plots.py (which generates static matplotlib figures). All plots are rendered as Plotly figures that support zoom, pan, hover tooltips, and data selection.

generate_interactive_plots(csv_file_path, dark_mode=False, selected_plots=None, x_var='b_imf_z', y_var='r_rc', marker_size_var=None)

Generate interactive Plotly figures from a statistics CSV file.

Parameters:

Name	Type	Description	Default
csv_file_path	str	Path to the CSV file containing statistics data.	required
dark_mode	bool	If True, use a dark background template.	False
selected_plots	list of str	Which plot types to generate. Options: "Histograms", "KDE Plots", "2D Histograms", "Scatter Plots", "MMS Location Scatter Plot".	None
x_var	str	Column name for the x-axis variable.	'b_imf_z'
y_var	str	Column name for the y-axis variable.	'r_rc'

Returns:

Name	Type	Description
figures	dict	Mapping of plot title to Plotly Figure object.
err	str or None	Error message if something went wrong, else None.

Source code in src/rxn_location/app_stats_plots_interactive.py

def generate_interactive_plots(
    csv_file_path,
    dark_mode=False,
    selected_plots=None,
    x_var="b_imf_z",
    y_var="r_rc",
    marker_size_var=None,
):
    """
    Generate interactive Plotly figures from a statistics CSV file.

    Parameters
    ----------
    csv_file_path : str
        Path to the CSV file containing statistics data.
    dark_mode : bool
        If True, use a dark background template.
    selected_plots : list of str
        Which plot types to generate. Options: "Histograms", "KDE Plots",
        "2D Histograms", "Scatter Plots", "MMS Location Scatter Plot".
    x_var : str
        Column name for the x-axis variable.
    y_var : str
        Column name for the y-axis variable.

    Returns
    -------
    figures : dict
        Mapping of plot title to Plotly Figure object.
    err : str or None
        Error message if something went wrong, else None.
    """
    if selected_plots is None:
        selected_plots = ["Scatter Plots"]

    figures = {}
    template = "plotly_dark" if dark_mode else "plotly_white"

    # Label mappings
    var_labels = {
        "b_imf_z": "IMF Bz [nT]",
        "b_imf_x": "IMF Bx [nT]",
        "b_imf_y": "IMF By [nT]",
        "imf_clock_angle": "IMF Clock Angle [°]",
        "cone_angle": "Cone Angle [°]",
        "p_dyn": "Dynamic Pressure [nPa]",
        "msh_msp_shear": "Shear Angle [°]",
        "r_rc": "Reconnection Distance [Rₑ]",
        "delta_beta": "Δβ",
    }

    x_label = var_labels.get(x_var, x_var.replace("_", " ").title())
    y_label = var_labels.get(y_var, y_var.replace("_", " ").title())

    # Read and prepare data
    try:
        df = pd.read_csv(csv_file_path, index_col=False)
        df = rcsf.convert_wide_to_long(df)
    except Exception as e:
        return figures, f"Error reading CSV: {e}"

    # Split by model
    model_names = ["Shear", "Reconnection Energy", "Exhaust Velocity", "Bisection"]
    model_filters = [
        ["shear", "Shear"],
        ["rx_en", "Reconnection Energy"],
        ["va_cs", "Exhaust Velocity"],
        ["bisection", "Bisection Field"],
    ]
    color_list = ["#1f77b4", "#ff7f0e", "#2ca02c", "#d62728"]

    df_list = []
    for filters in model_filters:
        df_model = df[df.method_used.isin(filters)].copy()
        # Compute derived columns if possible
        if len(df_model) > 0:
            if "imf_clock_angle" in df_model.columns:
                df_model["imf_clock_angle"] = df_model["imf_clock_angle"].apply(
                    lambda x: 360 - x if x > 180 else x
                )
            if all(c in df_model.columns for c in ["b_imf_x", "b_imf_y", "b_imf_z"]):
                mag = np.sqrt(
                    df_model.b_imf_x**2 + df_model.b_imf_y**2 + df_model.b_imf_z**2
                )
                df_model["cone_angle"] = np.arccos(df_model.b_imf_x / mag) * 180 / np.pi
                df_model["bb"] = df_model.b_imf_y / mag
        df_list.append(df_model)

    # ── Histograms ──────────────────────────────────────────────────────
    if "Histograms" in selected_plots:
        try:
            fig = make_subplots(
                rows=2,
                cols=2,
                subplot_titles=model_names,
                vertical_spacing=0.12,
                horizontal_spacing=0.08,
            )
            for i, (df_n, name) in enumerate(zip(df_list, model_names)):
                row, col = i // 2 + 1, i % 2 + 1
                if len(df_n) > 0 and "r_rc" in df_n.columns:
                    vals = df_n["r_rc"].dropna()
                    if len(vals) > 0:
                        fig.add_trace(
                            go.Histogram(
                                x=vals,
                                name=name,
                                marker_color=color_list[i],
                                opacity=0.8,
                                showlegend=False,
                            ),
                            row=row,
                            col=col,
                        )
                        fig.update_xaxes(
                            title_text="Reconnection Distance [Rₑ]" if row == 2 else "",
                            row=row,
                            col=col,
                        )
                        fig.update_yaxes(
                            title_text="Count" if col == 1 else "",
                            row=row,
                            col=col,
                        )

            fig.update_layout(
                template=template, paper_bgcolor="black" if dark_mode else "white", plot_bgcolor="black" if dark_mode else "white", font=dict(color="white" if dark_mode else "black"),
                height=700,
                title_text="Reconnection Distance Histograms",
                title_x=0.5,
            )
            figures["Histograms"] = fig
        except Exception as e:
            print(f"Interactive Histogram error: {e}")

    # ── KDE Plots ───────────────────────────────────────────────────────
    if "KDE Plots" in selected_plots:
        import plotly.figure_factory as ff
        try:
            fig = make_subplots(
                rows=2,
                cols=2,
                subplot_titles=model_names,
                vertical_spacing=0.12,
                horizontal_spacing=0.08,
            )
            for i, (df_n, name) in enumerate(zip(df_list, model_names)):
                row, col = i // 2 + 1, i % 2 + 1
                if len(df_n) > 0 and x_var in df_n.columns:
                    vals = df_n[x_var].dropna()
                    if len(vals) > 1:
                        # Create distplot just for the KDE curve
                        try:
                            kde_fig = ff.create_distplot([vals], [name], show_hist=False, show_rug=False, colors=[color_list[i]])
                            for trace in kde_fig['data']:
                                fig.add_trace(trace, row=row, col=col)
                        except Exception as e:
                            print(f"Skipping KDE for {name} due to error: {e}")

                fig.update_xaxes(
                    title_text=x_label if row == 2 else "",
                    row=row,
                    col=col,
                )
                fig.update_yaxes(
                    title_text="Density" if col == 1 else "",
                    row=row,
                    col=col,
                )

            fig.update_layout(
                template=template, paper_bgcolor="black" if dark_mode else "white", plot_bgcolor="black" if dark_mode else "white", font=dict(color="white" if dark_mode else "black"),
                height=700,
                title_text=f"KDE Plots: {x_label}",
                title_x=0.5,
                showlegend=False,
            )
            figures["KDE Plots"] = fig
        except Exception as e:
            print(f"Interactive KDE Plot error: {e}")

    # ── 2D Histograms (Heatmaps) ────────────────────────────────────────
    if "2D Histograms" in selected_plots:
        try:
            colorscales = ["Viridis", "Cividis", "Plasma", "Magma"]
            fig = make_subplots(
                rows=2,
                cols=2,
                subplot_titles=model_names,
                vertical_spacing=0.12,
                horizontal_spacing=0.10,
            )
            for i, (df_n, name) in enumerate(zip(df_list, model_names)):
                row, col = i // 2 + 1, i % 2 + 1
                if len(df_n) > 0 and x_var in df_n.columns and y_var in df_n.columns:
                    fig.add_trace(
                        go.Histogram2d(
                            x=df_n[x_var],
                            y=df_n[y_var],
                            colorscale=colorscales[i],
                            nbinsx=25,
                            nbinsy=25,
                            showscale=True,
                            name=name,
                        ),
                        row=row,
                        col=col,
                    )
                fig.update_xaxes(
                    title_text=x_label if row == 2 else "",
                    row=row,
                    col=col,
                )
                fig.update_yaxes(
                    title_text=y_label if col == 1 else "",
                    row=row,
                    col=col,
                )

            fig.update_layout(
                template=template, paper_bgcolor="black" if dark_mode else "white", plot_bgcolor="black" if dark_mode else "white", font=dict(color="white" if dark_mode else "black"),
                height=700,
                title_text=f"2D Histograms: {x_label} vs {y_label}",
                title_x=0.5,
            )
            figures["2D Histograms"] = fig
        except Exception as e:
            print(f"Interactive 2D Histogram error: {e}")

    # ── Scatter Plots ───────────────────────────────────────────────────
    if "Scatter Plots" in selected_plots:
        try:
            fig = make_subplots(
                rows=2,
                cols=2,
                subplot_titles=model_names,
                vertical_spacing=0.12,
                horizontal_spacing=0.08,
            )

            for i, (df_n, name) in enumerate(zip(df_list, model_names)):
                row, col = i // 2 + 1, i % 2 + 1
                if len(df_n) > 0 and x_var in df_n.columns and y_var in df_n.columns:
                    # Build hover text with key metadata
                    hover_parts = [
                        f"<b>{name}</b>",
                        f"{x_label}: %{{x:.3f}}",
                        f"{y_label}: %{{y:.3f}}",
                    ]
                    # Add jet_time if available
                    custom_data = []
                    extra_hover = []
                    if "jet_time" in df_n.columns:
                        custom_data.append(df_n["jet_time"].astype(str).values)
                        extra_hover.append("Jet Time: %{customdata[0]}")
                    if "Date" in df_n.columns:
                        custom_data.append(df_n["Date"].astype(str).values)
                        extra_hover.append(
                            "Date: %{customdata[" + str(len(custom_data) - 1) + "]}"
                        )

                    hovertemplate = (
                        "<br>".join(hover_parts + extra_hover) + "<extra></extra>"
                    )

                    scatter_kwargs = dict(
                        x=df_n[x_var],
                        y=df_n[y_var],
                        mode="markers",
                        name=name,
                        marker=dict(
                            color=color_list[i],
                            size=7,
                            opacity=0.7,
                            line=dict(
                                width=0.5, color="white" if dark_mode else "black"
                            ),
                        ),
                        hovertemplate=hovertemplate,
                        showlegend=False,
                    )
                    if custom_data:
                        scatter_kwargs["customdata"] = np.column_stack(custom_data)

                    fig.add_trace(go.Scatter(**scatter_kwargs), row=row, col=col)

                fig.update_xaxes(
                    title_text=x_label if row == 2 else "",
                    row=row,
                    col=col,
                )
                fig.update_yaxes(
                    title_text=y_label if col == 1 else "",
                    row=row,
                    col=col,
                )

            fig.update_layout(
                template=template, paper_bgcolor="black" if dark_mode else "white", plot_bgcolor="black" if dark_mode else "white", font=dict(color="white" if dark_mode else "black"),
                height=700,
                title_text=f"Scatter: {x_label} vs {y_label}",
                title_x=0.5,
            )
            figures["Scatter Plots"] = fig
        except Exception as e:
            print(f"Interactive Scatter Plot error: {e}")

    # ── MMS Location Scatter Plot ───────────────────────────────────────
    if "MMS Location Scatter Plot" in selected_plots:
        try:
            pos_cols_xy = ["spc_pos_x", "spc_pos_y"]
            pos_cols_yz = ["spc_pos_y", "spc_pos_z"]
            # Also check for master-list column names
            alt_pos_cols_xy = ["x_gsm", "y_gsm"]
            alt_pos_cols_yz = ["y_gsm", "z_gsm"]

            fig = make_subplots(
                rows=1,
                cols=2,
                subplot_titles=["Spacecraft Y vs X", "Spacecraft Z vs Y"],
                horizontal_spacing=0.10,
            )
            for i, (df_n, name) in enumerate(zip(df_list, model_names)):
                if len(df_n) == 0:
                    continue

                # Determine which column names exist
                if all(c in df_n.columns for c in pos_cols_xy):
                    x_col, y_col = pos_cols_xy
                    y2_col, z_col = pos_cols_yz
                elif all(c in df_n.columns for c in alt_pos_cols_xy):
                    x_col, y_col = alt_pos_cols_xy
                    y2_col, z_col = alt_pos_cols_yz
                else:
                    continue

                fig.add_trace(
                    go.Scatter(
                        x=df_n[x_col],
                        y=df_n[y_col],
                        mode="markers",
                        name=name,
                        marker=dict(color=color_list[i], size=6, opacity=0.7),
                        legendgroup=name,
                        showlegend=True,
                    ),
                    row=1,
                    col=1,
                )
                fig.add_trace(
                    go.Scatter(
                        x=df_n[y2_col],
                        y=df_n[z_col],
                        mode="markers",
                        name=name,
                        marker=dict(color=color_list[i], size=6, opacity=0.7),
                        legendgroup=name,
                        showlegend=False,
                    ),
                    row=1,
                    col=2,
                )

            fig.update_xaxes(title_text="MMS Position X [Rₑ]", row=1, col=1)
            fig.update_yaxes(title_text="MMS Position Y [Rₑ]", row=1, col=1)
            fig.update_xaxes(title_text="MMS Position Y [Rₑ]", row=1, col=2)
            fig.update_yaxes(title_text="MMS Position Z [Rₑ]", row=1, col=2)

            fig.update_layout(
                template=template, paper_bgcolor="black" if dark_mode else "white", plot_bgcolor="black" if dark_mode else "white", font=dict(color="white" if dark_mode else "black"),
                height=500,
                title_text="MMS Spacecraft Positions",
                title_x=0.5,
            )
            figures["MMS Location Scatter Plot"] = fig
        except Exception as e:
            print(f"Interactive MMS Location Plot error: {e}")

    # ── Interactive Joint-Plots ─────────────────────────────────────────
    if "Seaborn Joint-Plots" in selected_plots:
        try:
            # We can't do a 2x2 grid of joint plots easily in Plotly Subplots, 
            # so we create 4 separate figures and return them as a list under this key.
            # The frontend (app.py) will render them in a 2x2 st.columns grid.
            joint_figures = []
            import plotly.express as px

            for i, (df_n, name) in enumerate(zip(df_list, model_names)):
                if len(df_n) > 0 and x_var in df_n.columns and y_var in df_n.columns:
                    # Filter out NaNs
                    df_clean = df_n.dropna(subset=[x_var, y_var]).copy()

                    scatter_kwargs = {
                        "x": x_var,
                        "y": y_var,
                        "marginal_x": "histogram",
                        "marginal_y": "histogram",
                        "color_discrete_sequence": [color_list[i]],
                        "title": name
                    }
                    if marker_size_var and marker_size_var != "None" and marker_size_var in df_clean.columns:
                        scatter_kwargs["size"] = marker_size_var
                        # Ensure sizes are positive and non-zero for Plotly
                        df_clean[marker_size_var] = df_clean[marker_size_var].abs() + 0.1

                    if len(df_clean) > 0:
                        fig = px.scatter(df_clean, **scatter_kwargs)
                        fig.update_layout(
                            template=template, paper_bgcolor="black" if dark_mode else "white", plot_bgcolor="black" if dark_mode else "white", font=dict(color="white" if dark_mode else "black"),
                            xaxis_title=x_label,
                            yaxis_title=y_label,
                            height=500
                        )
                        joint_figures.append(fig)
                    else:
                        joint_figures.append(None)
                else:
                    joint_figures.append(None)

            figures["Seaborn Joint-Plots"] = joint_figures
        except Exception as e:
            print(f"Interactive Joint-Plots error: {e}")

    return figures, None

rxn_location.plotly_utils

convert_tplot_to_plotly(keys_to_plot, dark_mode=False)

Converts a list of tplot variables into a Plotly Figure.

Source code in src/rxn_location/plotly_utils.py

def convert_tplot_to_plotly(keys_to_plot, dark_mode=False):
    """
    Converts a list of tplot variables into a Plotly Figure.
    """
    num_plots = len(keys_to_plot)
    fig = make_subplots(
        rows=num_plots, cols=1, shared_xaxes=True, vertical_spacing=0.02
    )

    layout_updates = {}

    for i, key in enumerate(keys_to_plot):
        meta = spd.get_data(key, metadata=True)
        plot_opts = meta.get("plot_options", {}) if meta else {}

        yaxis_opt = plot_opts.get("yaxis_opt", {})
        zaxis_opt = plot_opts.get("zaxis_opt", {})

        y_title = format_latex(yaxis_opt.get("axis_label", key))
        y_subtitle = format_latex(yaxis_opt.get("axis_subtitle", ""))
        if y_subtitle:
            y_title += f"<br>{y_subtitle}"

        y_type = "log" if yaxis_opt.get("y_axis_type") == "log" else "linear"
        z_type = "log" if zaxis_opt.get("z_axis_type") == "log" else "linear"

        # Calculate vertical center for colorbar/legend
        y_center = 1 - (i + 0.5) / num_plots

        # Determine all variables to plot in this panel (handles pseudo-variables)
        vars_to_plot = plot_opts.get("overplots_mpl", [])
        if not vars_to_plot:
            vars_to_plot = [key]

        legend_names = plot_opts.get("line_opt", {}).get("legend_names", [])
        legend_idx = 0

        for var_name in vars_to_plot:
            data = spd.get_data(var_name)
            if data is None:
                continue

            times = [
                datetime.datetime.fromtimestamp(t, datetime.timezone.utc)
                for t in data.times
            ]

            # Check if it's a spectrogram (2D data with v)
            if hasattr(data, "v") or hasattr(data, "v1") or hasattr(data, "v2"):
                v_data = (
                    data.v
                    if hasattr(data, "v")
                    else (data.v1 if hasattr(data, "v1") else data.v2)
                )
                y_data = data.y

                z_title = format_latex(zaxis_opt.get("axis_label", ""))

                # For Plotly heatmap, z is the 2D array, x is time, y is v_data
                fig.add_trace(
                    go.Heatmap(
                        x=times,
                        y=v_data if v_data.ndim == 1 else v_data[0],
                        z=y_data.T,
                        colorscale="Viridis",
                        showscale=True,
                        colorbar=dict(
                            title=z_title,
                            len=1.0 / num_plots,
                            y=y_center,
                            yanchor="middle",
                            x=1.02,
                        ),
                    ),
                    row=i + 1,
                    col=1,
                )

            else:
                # Line plot
                y_data = data.y

                # Configure a separate legend for this subplot
                legend_id = f"legend{i+1}" if i > 0 else "legend"
                layout_updates[legend_id] = dict(
                    y=y_center,
                    yanchor="middle",
                    x=1.02,
                    xanchor="left",
                    tracegroupgap=0,
                )

                if y_data.ndim == 1:
                    name = format_latex(
                        legend_names[legend_idx]
                        if legend_idx < len(legend_names)
                        else var_name
                    )
                    legend_idx += 1
                    fig.add_trace(
                        go.Scatter(
                            x=times, y=y_data, mode="lines", name=name, legend=legend_id
                        ),
                        row=i + 1,
                        col=1,
                    )
                else:
                    for col in range(y_data.shape[1]):
                        name = format_latex(
                            legend_names[legend_idx]
                            if legend_idx < len(legend_names)
                            else f"{var_name}_{col}"
                        )
                        legend_idx += 1
                        fig.add_trace(
                            go.Scatter(
                                x=times,
                                y=y_data[:, col],
                                mode="lines",
                                name=name,
                                legend=legend_id,
                            ),
                            row=i + 1,
                            col=1,
                        )

        if y_type == "log":
            fig.update_yaxes(
                type=y_type,
                title_text=y_title,
                exponentformat="power",
                dtick=1,
                row=i + 1,
                col=1,
            )
        else:
            fig.update_yaxes(type=y_type, title_text=y_title, row=i + 1, col=1)

    template = "plotly_dark" if dark_mode else "plotly_white"
    layout_updates.update(
        dict(height=250 * num_plots, template=template, showlegend=True)
    )
    fig.update_layout(**layout_updates)

    return fig

format_latex(text)

Formats text strings into LaTeX math strings for Plotly annotations.

Parameters:

Name	Type	Description	Default
text	str	The raw text string, which may contain underscores or greek letter names.	required

Returns:

Type	Description
str	The formatted string with Plotly-friendly HTML (e.g., /) for rendering.

Source code in src/rxn_location/plotly_utils.py

def format_latex(text):
    """
    Formats text strings into LaTeX math strings for Plotly annotations.

    Parameters
    ----------
    text : str
        The raw text string, which may contain underscores or greek letter names.

    Returns
    -------
    str
        The formatted string with Plotly-friendly HTML (e.g., <sub>/<sup>) for rendering.
    """
    if not text:
        return text
    text = str(text)

    is_latex = "$" in text

    # Remove math mode $
    text = text.replace("$", "")
    # Remove \rm
    text = text.replace("\\rm ", "")
    text = text.replace("\\rm", "")

    # Replace common symbols
    text = text.replace("\\Delta ", "Δ")
    text = text.replace("\\Delta", "Δ")
    text = text.replace("\\parallel ", "∥")
    text = text.replace("\\parallel", "∥")
    text = text.replace("\\perp ", "⟂")
    text = text.replace("\\perp", "⟂")

    if is_latex:
        # Handle subscripts: _{...} or _a
        text = re.sub(r"_\{([^}]+)\}", r"<sub>\1</sub>", text)
        text = re.sub(r"_([a-zA-Z0-9])", r"<sub>\1</sub>", text)

        # Handle superscripts: ^{...} or ^a
        text = re.sub(r"\^\{([^}]+)\}", r"<sup>\1</sup>", text)
        text = re.sub(r"\^([a-zA-Z0-9])", r"<sup>\1</sup>", text)

    return text

Command Line Interface

rxn_location.batch_statistics_cli

main()

The main entry point for the batch statistics processing script.

This function loads the requested configurations, reads the input times list, loops through each target timestamp, runs the jet reversal checks and reconnection models, and logs valid detections to the master jet list. Finally, it can optionally trigger Seaborn statistical plots.

Source code in src/rxn_location/batch_statistics_cli.py

def main():
    """
    The main entry point for the batch statistics processing script.

    This function loads the requested configurations, reads the input times list,
    loops through each target timestamp, runs the jet reversal checks and 
    reconnection models, and logs valid detections to the master jet list.
    Finally, it can optionally trigger Seaborn statistical plots.
    """
    args = parse_args()

    # Print default/passed parameters using Tabulate
    param_table = [
        ["Input File", args.input],
        ["MMS Probe", args.probe],
        ["Tsyganenko Model", args.tsy_model],
        ["Data Rate", args.data_rate],
        ["Plot Format", args.format],
        ["Output Directory", args.outdir],
        ["Stats CSV Name", args.csv_name],
        ["Jet Window (dt)", args.dt],
        ["Jet Length Threshold", args.jet_len],
        ["MP Standoff (m_p)", args.m_p],
        ["Grid Resolution (dr)", args.dr],
        ["Grid Limits", args.limits],
        ["Verbosity", args.verbosity],
        ["Retry Δt (min)", args.t_delta],
        ["Max Retries", args.max_retries],
    ]
    vprint(1, "\n" + "=" * 60, color="cyan")
    vprint(1, "RXN Location: Batch Statistics Processing", color="bold")
    vprint(1, "=" * 60, color="cyan")
    vprint(
        1,
        tabulate(param_table, headers=["Parameter", "Value"], tablefmt="fancy_grid"),
        color="cyan",
    )
    vprint(1, "=" * 60 + "\n", color="cyan")

    # Configure global verbosity
    set_verbosity(args.verbosity)

    if not os.path.exists(args.input):
        vprint(1, f"Error: Input file '{args.input}' not found.", color="red")
        sys.exit(1)

    os.makedirs(args.outdir, exist_ok=True)

    # Read times
    file_times = []
    with open(args.input, "r") as f:
        for line in f:
            line = line.strip()
            if not line or line.startswith(","):
                continue
            first_col = line.split(",")[0]
            try:
                t = parser.parse(first_col)
                if t.tzinfo is None:
                    t = t.replace(tzinfo=pytz.utc)
                file_times.append(t)
            except Exception:
                pass

    if not file_times:
        vprint(1, "Error: No valid timestamps found in the input file.", color="red")
        sys.exit(1)

    vprint(1, f"Found {len(file_times)} valid timestamps to process.\n", color="green")

    master_jets = mjl.load_master_list()

    # Define models to run
    models_to_run = ["shear", "reconnection energy", "exhaust velocity", "bisection"]

    for idx, c_time_original in enumerate(file_times):
        time_str = c_time_original.strftime("%Y-%m-%d %H:%M:%S")
        vprint(1, f"[{idx+1}/{len(file_times)}] Processing: {time_str}", color="blue")

        # Check if time is already in master list
        already_processed = False
        for entry in master_jets:
            try:
                ct = entry.get("crossing_time")
                if ct:
                    ct_dt = mjl._parse_time(ct)
                    if abs(ct_dt - c_time_original) <= datetime.timedelta(seconds=2):
                        already_processed = True
                        break
            except Exception:
                pass

        if not already_processed:
            existing = mjl.find_nearby_jet(
                master_jets, c_time_original, window_minutes=2
            )
            if existing is not None:
                already_processed = True

        if already_processed:
            vprint(
                1,
                f"  -> Skipping: Time {time_str} is already present in the Master List.",
                color="yellow",
            )
            continue

        # 1. Jet Reversal Check (with retry logic)
        # Try the original time first, then shift by ±t_delta minutes
        fig, s1, det = None, False, None
        c_time = c_time_original  # the time that actually found the jet

        # Build list of times to try: original, then ±1*delta, ±2*delta, ...
        times_to_try = [c_time_original]
        for i in range(1, args.max_retries + 1):
            times_to_try.append(
                c_time_original + datetime.timedelta(minutes=i * args.t_delta)
            )
            times_to_try.append(
                c_time_original - datetime.timedelta(minutes=i * args.t_delta)
            )

        for attempt_idx, attempt_time in enumerate(times_to_try):
            time_str_safe = attempt_time.strftime("%Y-%m-%d_%H%M%S")
            jet_dir = os.path.join(args.outdir, "jet_reversal_checks", f"mms{args.probe}")
            os.makedirs(jet_dir, exist_ok=True)
            jet_plot_filename = os.path.join(
                jet_dir, f"jet_plot_{time_str_safe}.{args.format}"
            )

            try:
                res = jet_reversal_check(
                    attempt_time,
                    probe=args.probe,
                    data_rate=args.data_rate,
                    dt=args.dt,
                    jet_len=args.jet_len,
                    return_plotly_fig=(args.plot_engine == "plotly"),
                    figname=(
                        jet_plot_filename if args.plot_engine == "matplotlib" else None
                    ),
                )
            except Exception as e:
                if attempt_idx == 0:
                    vprint(1, f"  -> Error running Jet Check: {e}", color="red")
                continue

            if res is None:
                continue

            fig_try, s1_try, det_try = res

            if s1_try and det_try is not None:
                fig, s1, det = fig_try, s1_try, det_try
                c_time = attempt_time
                if attempt_idx > 0:
                    vprint(
                        1,
                        f"  -> Jet found at shifted time {attempt_time.strftime('%Y-%m-%d %H:%M:%S')} "
                        f"(attempt {attempt_idx + 1}/{len(times_to_try)})",
                        color="yellow",
                    )
                break
            elif attempt_idx == 0:
                vprint(
                    2,
                    f"  -> No jet at original time, searching nearby times (±{args.t_delta}min steps)...",
                    color="magenta",
                )

        if not s1 or det is None:
            vprint(
                1,
                f"  -> No jet found at {time_str} or within ±{args.max_retries * args.t_delta}min window.",
                color="red",
            )
            continue

        vprint(1, f"  -> Jet detected!", color="green")

        # Save the jet plot
        if args.plot_engine == "plotly" and fig is not None:
            try:
                if args.format == "html":
                    fig.write_html(jet_plot_filename)
                else:
                    fig.write_image(jet_plot_filename)

                abs_jet_plot_filename = os.path.abspath(jet_plot_filename)
                vprint(
                    2,
                    f"  -> Saved Jet Plot to {abs_jet_plot_filename}",
                    color="magenta",
                )
            except Exception as e:
                vprint(1, f"  -> Failed to save Jet Plot: {e}", color="red")
        elif args.plot_engine == "matplotlib":
            abs_jet_plot_filename = os.path.abspath(jet_plot_filename)
            vprint(
                2, f"  -> Saved Jet Plot to {abs_jet_plot_filename}", color="magenta"
            )

        # 2. Reconnection Models (use exact jet time from detection)
        exact_jet_time = det.get("jet_time")
        if exact_jet_time is None:
            exact_jet_time = c_time
        jet_time_str = exact_jet_time.strftime("%Y-%m-%d %H:%M:%S")
        try:
            model_inputs = {
                "trange": [jet_time_str],
                "probe": None,
                "omni_level": "hro",
                "mms_probe_num": args.probe,
                "model_type": args.tsy_model.lower(),
                "m_p": args.m_p,
                "dr": args.dr,
                "min_max_val": args.limits,
            }

            res = rx_model(**model_inputs)
            if not res:
                vprint(1, f"  -> Model generation failed for {time_str}", color="red")
                continue

            sw_params = res[8]

            images, c_labels = [], []
            model_mapping = {
                "shear": {"var_idx": 3, "label": "Shear"},
                "bisection": {"var_idx": 6, "label": "Bisection Field"},
                "reconnection energy": {"var_idx": 4, "label": "Reconnection Energy"},
                "exhaust velocity": {"var_idx": 5, "label": "Exhaust Velocity"},
            }

            for m in models_to_run:
                var_idx = model_mapping[m]["var_idx"]
                raw_data = res[var_idx]
                norm_data = (raw_data - np.nanmin(raw_data)) / (
                    np.nanmax(raw_data) - np.nanmin(raw_data)
                )
                images.append(norm_data)
                c_labels.append(model_mapping[m]["label"])

            figure_inputs = {
                "image": images,
                "convolution_order": [1] * len(images),
                "t_range": [jet_time_str],
                "b_imf": np.round(sw_params["b_imf"], 2),
                "b_msh": np.round(sw_params["mms_b_gsm"], 2),
                "xrange": [-args.limits, args.limits],
                "yrange": [-args.limits, args.limits],
                "mms_probe_num": str(args.probe),
                "mms_sc_pos": np.round(sw_params["mms_sc_pos"], 2),
                "dr": args.dr,
                "dipole_tilt_angle": sw_params["ps"],
                "p_dyn": np.round(sw_params["p_dyn"], 2),
                "imf_clock_angle": sw_params["imf_clock_angle"],
                "np_imf": np.round(sw_params["np"], 2),
                "v_imf": np.round(sw_params["v_imf"], 2),
                "sym_h": np.round(sw_params["sym_h"], 2),
                "sigma": [2] * len(images),
                "mode": "nearest",
                "alpha": 1,
                "vmin": [0] * len(images),
                "vmax": [1] * len(images),
                "cmap_list": ["viridis", "cividis", "plasma", "magma"][: len(images)],
                "draw_patch": [True] * len(images),
                "draw_ridge": [True] * len(images),
                "save_fig": True,
                "fig_name": f"recon_models_{time_str.replace(' ', '_').replace(':', '')}",
                "fig_format": args.format,
                "c_label": c_labels,
                "wspace": 0.0,
                "hspace": 0.20,
                "fig_size": (8.775, 10) if args.format == "html" else (10, 8),
                "box_style": dict(boxstyle="round", color="k", alpha=0.8),
                "title_y_pos": 1.09,
                "interpolation": "None",
                "tsy_model": args.tsy_model.lower(),
                "dark_mode": False,
                "save_rc_file": True,
                "rc_file_name": args.csv_name,
                "rc_folder": args.outdir,
                "df_jet_reversal": det,
                "b_grids": (
                    (res[0], res[1], res[2], res[12], res[13], res[14])
                    if "res" in locals() and res
                    else None
                ),
            }

            vprint(2, "  -> Generating Reconnection Models...", color="magenta")
            if args.format == "html":
                _, dist_rc_dict = ridge_finder_multiple_interactive(**figure_inputs)
            else:
                _, _, _, dist_rc_dict = ridge_finder_multiple(**figure_inputs)

            if dist_rc_dict:
                det.update(dist_rc_dict)

            vprint(2, "  -> Models generated and saved!", color="magenta")

        except Exception as e:
            vprint(1, f"  -> Error running Models: {e}", color="red")
            continue

        # 3. Save to Master List
        run_params = {
            "mms_probe": int(args.probe),
            "dt": args.dt,
            "jet_len": args.jet_len,
            "data_rate": args.data_rate,
            "level": "l2",
            "coord_type": "lmn",
            "time_clip": True,
            "t_delta": 10,
            "max_attempts": 5,
            "tsy_model": args.tsy_model,
            "recon_models": models_to_run,
            "omni_level": "hro",
            "m_p": args.m_p,
            "dr": args.dr,
            "limits": args.limits,
        }

        # Pull latest OMNI SW params for contextual Master List logging
        try:
            trange_date_min = exact_jet_time - datetime.timedelta(minutes=30)
            trange_date_max = exact_jet_time + datetime.timedelta(minutes=30)
            trange_min_str = trange_date_min.strftime("%Y-%m-%d %H:%M:%S") + "Z"
            trange_max_str = trange_date_max.strftime("%Y-%m-%d %H:%M:%S") + "Z"

            sw_params = get_sw_params(
                trange=[trange_min_str, trange_max_str],
                omni_level="hro",
                mms_probe_num=str(args.probe),
            )
            if sw_params is not None:
                det["sw_b_imf_gsm_x"] = sw_params["b_imf"][0]
                det["sw_b_imf_gsm_y"] = sw_params["b_imf"][1]
                det["sw_b_imf_gsm_z"] = sw_params["b_imf"][2]
                det["sw_v_imf_gse_x"] = sw_params["v_imf"][0]
                det["sw_v_imf_gse_y"] = sw_params["v_imf"][1]
                det["sw_v_imf_gse_z"] = sw_params["v_imf"][2]
                det["sw_np"] = sw_params["np"]
                det["sw_tp"] = sw_params["t_p"]
                det["sw_sym_h"] = sw_params["sym_h"]
                det["sw_clock_angle"] = sw_params["imf_clock_angle"]
                det["sw_p_dyn"] = sw_params["p_dyn"]

                # Compute cone angle
                import math

                bx, by, bz = (
                    sw_params["b_imf"][0],
                    sw_params["b_imf"][1],
                    sw_params["b_imf"][2],
                )
                b_mag = math.sqrt(bx**2 + by**2 + bz**2)
                if b_mag > 0:
                    det["sw_cone_angle"] = math.acos(bx / b_mag) * 180 / math.pi
        except Exception:
            pass

        was_added, existing = mjl.add_jet(
            master_jets, det, c_time, run_params, window_minutes=2
        )

        if was_added:
            mjl.save_master_list(master_jets)
            vprint(1, "  -> Logged to Master Jet List.", color="green")
        else:
            vprint(
                2,
                f"  -> Jet already exists in Master List (near {existing.get('jet_time')}).",
                color="yellow",
            )

    if args.plot_seaborn:
        vprint(1, "\n" + "=" * 60, color="cyan")
        vprint(
            1,
            "Generating Seaborn joint-plots for master list parameters...",
            color="bold",
        )
        vprint(1, "=" * 60, color="cyan")
        df_stats = mjl.master_list_to_stats_csv(master_jets)
        if df_stats is not None and not df_stats.empty:
            from rxn_location.app_seaborn_plots import generate_seaborn_jointplots

            try:
                _ = generate_seaborn_jointplots(df_full=df_stats, dark_mode=False)
                vprint(1, "  -> Done generating Seaborn plots.", color="green")
            except Exception as e:
                vprint(1, f"  -> Error generating Seaborn plot: {e}", color="red")
        else:
            vprint(1, "  -> Not enough data to generate Seaborn plots.", color="yellow")

    vprint(
        1, f"\nBatch processing complete! Output saved to {args.outdir}/", color="green"
    )

parse_args()

Parses command line arguments for the batch statistics script.

Returns:

Type	Description
Namespace	The parsed arguments namespace containing all configuration options.

Source code in src/rxn_location/batch_statistics_cli.py

def parse_args():
    """
    Parses command line arguments for the batch statistics script.

    Returns
    -------
    argparse.Namespace
        The parsed arguments namespace containing all configuration options.
    """
    class CommentedArgumentParser(argparse.ArgumentParser):
        """
        Custom ArgumentParser that strips out inline comments from configuration files.
        """
        def convert_arg_line_to_args(self, arg_line):
            """
            Reads a single line from a configuration file, stripping comments and yielding valid arguments.
            """
            arg_line = arg_line.split("#")[0].strip()
            if arg_line:
                yield arg_line

    parser = CommentedArgumentParser(
        description="Batch process jet reversal checks and reconnection models from a list of times.",
        fromfile_prefix_chars="@",
    )
    parser.add_argument(
        "-i",
        "--input",
        required=True,
        help="Path to input .txt or .csv file with a list of times.",
    )
    parser.add_argument(
        "--probe", type=str, default="3", help="MMS probe number (1-4). Default: 3"
    )
    parser.add_argument(
        "--tsy_model",
        type=str,
        default="T96",
        help="Tsyganenko model to use. Default: T96",
    )
    parser.add_argument(
        "--data_rate",
        type=str,
        default="fast",
        help="MMS data rate (fast/brst). Default: fast",
    )
    parser.add_argument(
        "--format",
        choices=["html", "pdf", "png"],
        default="html",
        help="Format for saved plots (html for interactive Plotly, pdf/png for static Matplotlib). Default: html",
    )
    parser.add_argument(
        "--outdir",
        type=str,
        default="./figures",
        help="Output directory for plots. Default: ./figures",
    )
    parser.add_argument(
        "--csv_name",
        type=str,
        default="batch_reconnection_stats.csv",
        help="Name of output stats CSV.",
    )
    parser.add_argument(
        "--dt",
        type=int,
        default=300,
        help="Time window around crossing time to search for jets (s). Default: 300",
    )
    parser.add_argument(
        "--jet_len",
        type=int,
        default=3,
        help="Minimum number of data points for a valid jet. Default: 3",
    )
    parser.add_argument(
        "--m_p",
        type=float,
        default=0.5,
        help="Magnetopause standoff distance scaling. Default: 0.5",
    )
    parser.add_argument(
        "--dr",
        type=float,
        default=0.25,
        help="Spatial resolution for 3D model grid (Re). Default: 0.25",
    )
    parser.add_argument(
        "--limits",
        type=int,
        default=20,
        help="X, Y, Z boundaries for 3D grid (Re). Default: 20",
    )
    parser.add_argument(
        "--verbosity",
        type=int,
        default=2,
        choices=[0, 1, 2, 3],
        help="Verbosity level (0-3). 0: silent, 1: important only, 2: standard (no PySPEDAS), 3: all. Default: 2",
    )
    parser.add_argument(
        "--t_delta",
        type=int,
        default=2,
        help="Minutes to shift the crossing time on each retry when searching for jets. Default: 2",
    )
    parser.add_argument(
        "--max_retries",
        type=int,
        default=5,
        help="Maximum number of time-shifted retries when initial jet check fails. Default: 5",
    )
    parser.add_argument(
        "--plot-seaborn",
        action="store_true",
        help="Generate final statistical Seaborn joint-plots from the Master List.",
    )
    parser.add_argument(
        "--plot_engine",
        type=str,
        default="plotly",
        choices=["plotly", "matplotlib"],
        help="Engine used to plot the jet reversal check (plotly or matplotlib). Default: plotly",
    )
    return parser.parse_args()