Magnetic Field Utilities

el_paso.processing.magnetic_field_utils

Classes

el_paso.processing.magnetic_field_utils.Coords

A class to handle coordinate transformations using the IRBEM library.

This class provides a Pythonic interface for converting coordinates between different systems (e.g., GEO, GSE, GSM) by calling the underlying Fortran routines of the IRBEM-LIB.

Source code in el_paso/processing/magnetic_field_utils/irbem.py

class Coords:
    """A class to handle coordinate transformations using the IRBEM library.

    This class provides a Pythonic interface for converting coordinates between
    different systems (e.g., GEO, GSE, GSM) by calling the underlying
    Fortran routines of the IRBEM-LIB.
    """

    def __init__(self, *, lib_path: Optional[str | Path] = DEFAULT_LIBIRBEM_PATH) -> None:
        """Initializes the Coords object and loads the IRBEM shared library.

        Args:
            lib_path (Optional[str | Path]): The path to the IRBEM shared library file.
                                                    Defaults to DEFAULT_LIBIRBEM_PATH.
        """
        if isinstance(lib_path, str):
            lib_path = Path(lib_path)

        self.irbem_obj_path = lib_path

        self.irbem_obj_path, self._irbem_obj = _load_shared_object(self.irbem_obj_path)

    def transform(
        self,
        time: list[datetime] | list[str] | datetime | str,
        pos: NDArray[np.floating],
        sysaxes_in: int | str,
        sysaxes_out: int | str,
    ) -> NDArray[np.float64]:
        """Transforms coordinates from one system to another.

        This method is a vectorized wrapper around the `coord_trans_vec1_` Fortran subroutine,
        which handles multiple input positions and times efficiently.

        Args:
            time (list[datetime] | list[str] | datetime | str): A single datetime object or a list of
                                                                datetime objects or ISO-formatted strings.
            pos (NDArray[np.floating]): An array of input positions. The shape should be (N, 3), where
                                        N is the number of time steps.
            sysaxes_in (int | str): The integer code or string name of the input coordinate system
                                    (e.g., 0 or 'GDZ' for GEI/GEO).
            sysaxes_out (int | str): The integer code or string name of the output coordinate system.

        Returns:
            NDArray[np.float64]: A NumPy array of the transformed positions with shape (N, 3).
        """
        if isinstance(time, (datetime)):
            time = [time]
        if isinstance(time, (str)):
            time = [time]

        pos = np.atleast_2d(pos)

        c_pos_in = ((ctypes.c_double * 3) * len(time))()
        c_pos_out = ((ctypes.c_double * 3) * len(time))()
        c_iyear, c_idoy, c_ut = self._convert_to_c_times(time)
        c_sys_in = self._get_c_sysaxes(sysaxes_in)
        c_sys_out = self._get_c_sysaxes(sysaxes_out)
        c_ntime = ctypes.c_int(len(time))

        for it, ix in itertools.product(range(pos.shape[0]), range(3)):
            c_pos_in[it][ix] = ctypes.c_double(pos[it, ix])

        self._irbem_obj.coord_trans_vec1_(
            ctypes.byref(c_ntime),
            ctypes.byref(c_sys_in),
            ctypes.byref(c_sys_out),
            ctypes.byref(c_iyear),
            ctypes.byref(c_idoy),
            ctypes.byref(c_ut),
            ctypes.byref(c_pos_in),
            ctypes.byref(c_pos_out),
        )
        return np.array(c_pos_out)

    def _convert_to_c_times(
        self, time: list[datetime] | list[str] | datetime | str
    ) -> tuple[ctypes.Array[ctypes.c_int], ctypes.Array[ctypes.c_int], ctypes.Array[ctypes.c_double]]:
        if isinstance(time, (datetime)):
            time = [time]
        if isinstance(time, (str)):
            time = [time]

        iyear = (ctypes.c_int * len(time))()
        idoy = (ctypes.c_int * len(time))()
        ut = (ctypes.c_double * len(time))()

        if isinstance(time[0], str):
            time = typing.cast("list[str]", time)
            time = [dateutil.parser.parse(t) for t in time]

        assert isinstance(time[0], datetime)
        time = typing.cast("list[datetime]", time)

        for it, t in enumerate(time):
            iyear[it] = ctypes.c_int(t.year)
            idoy[it] = ctypes.c_int(t.timetuple().tm_yday)
            ut[it] = ctypes.c_double(3600 * t.hour + 60 * t.minute + t.second)

        return iyear, idoy, ut

    def _get_c_sysaxes(self, sysaxes: int | str) -> ctypes.c_int:
        if isinstance(sysaxes, str):
            assert sysaxes.upper() in SYSAXES_STR_TO_INT, (
                "ERROR: Unknown coordinate system! Choose from GDZ, GEO, GSM, GSE, SM, GEI, MAG, SPH, RLL."
            )
            return ctypes.c_int(SYSAXES_STR_TO_INT[sysaxes])
        if isinstance(sysaxes, int):
            return ctypes.c_int(sysaxes)
        msg = "Error, coordinate axis can only be a string or int!"
        raise ValueError(msg)

Methods:

__init__

__init__

Initializes the Coords object and loads the IRBEM shared library.

Parameters:

Name	Type	Description	Default
lib_path	Optional[str | Path]	The path to the IRBEM shared library file. Defaults to DEFAULT_LIBIRBEM_PATH.	DEFAULT_LIBIRBEM_PATH

Source code in el_paso/processing/magnetic_field_utils/irbem.py

def __init__(self, *, lib_path: Optional[str | Path] = DEFAULT_LIBIRBEM_PATH) -> None:
    """Initializes the Coords object and loads the IRBEM shared library.

    Args:
        lib_path (Optional[str | Path]): The path to the IRBEM shared library file.
                                                Defaults to DEFAULT_LIBIRBEM_PATH.
    """
    if isinstance(lib_path, str):
        lib_path = Path(lib_path)

    self.irbem_obj_path = lib_path

    self.irbem_obj_path, self._irbem_obj = _load_shared_object(self.irbem_obj_path)

transform

transform

Transforms coordinates from one system to another.

This method is a vectorized wrapper around the coord_trans_vec1_ Fortran subroutine, which handles multiple input positions and times efficiently.

Parameters:

Name	Type	Description	Default
time	list[datetime] | list[str] | datetime | str	A single datetime object or a list of datetime objects or ISO-formatted strings.	required
pos	NDArray[floating]	An array of input positions. The shape should be (N, 3), where N is the number of time steps.	required
sysaxes_in	int | str	The integer code or string name of the input coordinate system (e.g., 0 or 'GDZ' for GEI/GEO).	required
sysaxes_out	int | str	The integer code or string name of the output coordinate system.	required

Returns:

Type	Description
NDArray[float64]	NDArray[np.float64]: A NumPy array of the transformed positions with shape (N, 3).

Source code in el_paso/processing/magnetic_field_utils/irbem.py

def transform(
    self,
    time: list[datetime] | list[str] | datetime | str,
    pos: NDArray[np.floating],
    sysaxes_in: int | str,
    sysaxes_out: int | str,
) -> NDArray[np.float64]:
    """Transforms coordinates from one system to another.

    This method is a vectorized wrapper around the `coord_trans_vec1_` Fortran subroutine,
    which handles multiple input positions and times efficiently.

    Args:
        time (list[datetime] | list[str] | datetime | str): A single datetime object or a list of
                                                            datetime objects or ISO-formatted strings.
        pos (NDArray[np.floating]): An array of input positions. The shape should be (N, 3), where
                                    N is the number of time steps.
        sysaxes_in (int | str): The integer code or string name of the input coordinate system
                                (e.g., 0 or 'GDZ' for GEI/GEO).
        sysaxes_out (int | str): The integer code or string name of the output coordinate system.

    Returns:
        NDArray[np.float64]: A NumPy array of the transformed positions with shape (N, 3).
    """
    if isinstance(time, (datetime)):
        time = [time]
    if isinstance(time, (str)):
        time = [time]

    pos = np.atleast_2d(pos)

    c_pos_in = ((ctypes.c_double * 3) * len(time))()
    c_pos_out = ((ctypes.c_double * 3) * len(time))()
    c_iyear, c_idoy, c_ut = self._convert_to_c_times(time)
    c_sys_in = self._get_c_sysaxes(sysaxes_in)
    c_sys_out = self._get_c_sysaxes(sysaxes_out)
    c_ntime = ctypes.c_int(len(time))

    for it, ix in itertools.product(range(pos.shape[0]), range(3)):
        c_pos_in[it][ix] = ctypes.c_double(pos[it, ix])

    self._irbem_obj.coord_trans_vec1_(
        ctypes.byref(c_ntime),
        ctypes.byref(c_sys_in),
        ctypes.byref(c_sys_out),
        ctypes.byref(c_iyear),
        ctypes.byref(c_idoy),
        ctypes.byref(c_ut),
        ctypes.byref(c_pos_in),
        ctypes.byref(c_pos_out),
    )
    return np.array(c_pos_out)

el_paso.processing.magnetic_field_utils.IrbemInput dataclass

A data class to hold all necessary input parameters for IRBEM calculations.

Attributes:

Name	Type	Description
magnetic_field	MagneticField	The magnetic field model to be used.
maginput	dict[MagInputKeys, NDArray[float64]]	A dictionary of magnetic field input parameters required by IRBEM (e.g., Kp, Dst).
irbem_options	list[int]	A list of integer options to configure the IRBEM library's behavior.
num_cores	int	The number of CPU cores to use for parallel processing. Defaults to 4.
irbem_lib_path	str | Path	The file path to the compiled IRBEM library. Defaults to the 'libirbem.so' located in the same directory as the el_paso package.

Source code in el_paso/processing/magnetic_field_utils/magnetic_field_functions.py

@dataclass
class IrbemInput:
    """A data class to hold all necessary input parameters for IRBEM calculations.

    Attributes:
        magnetic_field (MagneticField): The magnetic field model to be used.
        maginput (dict[MagInputKeys, NDArray[np.float64]]): A dictionary of
            magnetic field input parameters required by IRBEM (e.g., Kp, Dst).
        irbem_options (list[int]): A list of integer options to configure the
            IRBEM library's behavior.
        num_cores (int): The number of CPU cores to use for parallel processing.
            Defaults to 4.
        irbem_lib_path (str|Path): The file path to the compiled IRBEM library.
            Defaults to the 'libirbem.so' located in the same directory as the el_paso package.
    """

    magnetic_field: MagneticField
    maginput: dict[MagInputKeys, NDArray[np.float64]]
    irbem_options: list[int]
    num_cores: int = 4
    irbem_lib_path: str | Path = str(Path(ep.__file__).parent / "libirbem.so")

el_paso.processing.magnetic_field_utils.IrbemOutput

Bases: NamedTuple

A named tuple to represent the output of a single IRBEM calculation.

This is used for intermediate results in parallel processing.

Attributes:

Name	Type	Description
arr	NDArray[float64]	The calculated data as a NumPy array.
unit	UnitBase	The unit of the calculated data.

Source code in el_paso/processing/magnetic_field_utils/magnetic_field_functions.py

class IrbemOutput(NamedTuple):
    """A named tuple to represent the output of a single IRBEM calculation.

    This is used for intermediate results in parallel processing.

    Attributes:
        arr (NDArray[np.float64]): The calculated data as a NumPy array.
        unit (u.UnitBase): The unit of the calculated data.
    """

    arr: NDArray[np.float64]
    unit: u.UnitBase

el_paso.processing.magnetic_field_utils.MagneticField

Bases: Enum

Enum for magnetic field models.

Source code in el_paso/processing/magnetic_field_utils/mag_field_enum.py

class MagneticField(Enum):
    """Enum for magnetic field models."""

    T89 = "T89"
    T01 = "T01"
    T01s = "T01s"
    TS04 = "TS04"
    TS05 = "TS05"
    T04s = "T04s"
    T96 = "T96"
    OP77Q = "OP77Q"
    OP77 = "OP77"

    def kext(self) -> kext:
        """Returns the kext value for the magnetic field model."""
        return _magnetic_field_str_to_kext(self.value)

    @classmethod
    def _missing_(cls, value: object) -> None:
        msg = "{!r} is not a valid {}.  Valid types: {}".format(
            value,
            cls.__name__,
            ", ".join([repr(m.value) for m in cls]),
        )
        raise ValueError(msg)

Methods:

kext

kext

Returns the kext value for the magnetic field model.

Source code in el_paso/processing/magnetic_field_utils/mag_field_enum.py

def kext(self) -> kext:
    """Returns the kext value for the magnetic field model."""
    return _magnetic_field_str_to_kext(self.value)

Functions:

el_paso.processing.magnetic_field_utils.create_var_name

create_var_name

Creates a standardized variable name combining the variable type and magnetic field model.

This function is used internally to generate consistent variable names for all output quantities based on the magnetic field model used.

Parameters:

Name	Type	Description	Default
var_type	MagFieldVarTypes	The type of the magnetic field-related variable (e.g., "B_eq", "MLT").	required
mag_field	MagneticField	The specific magnetic field model used for the calculation.	required

Returns:

Name	Type	Description
str	str	The concatenated and standardized variable name.

Source code in el_paso/processing/magnetic_field_utils/magnetic_field_functions.py

def create_var_name(var_type: MagFieldVarTypes, mag_field: MagneticField) -> str:
    """Creates a standardized variable name combining the variable type and magnetic field model.

    This function is used internally to generate consistent variable names
    for all output quantities based on the magnetic field model used.

    Args:
        var_type (MagFieldVarTypes): The type of the magnetic field-related variable
                                     (e.g., "B_eq", "MLT").
        mag_field (MagneticField): The specific magnetic field model used for the calculation.

    Returns:
        str: The concatenated and standardized variable name.
    """
    return var_type + "_" + mag_field.value

el_paso.processing.magnetic_field_utils.get_Lstar

get_Lstar

Calculates Lm, Lstar, and the third adiabatic invariant (J).

This function computes Lm and Lstar for each satellite position and local pitch angle. These are crucial parameters for characterizing particle drift shells in the magnetosphere. The calculation is parallelized for performance.

Parameters:

Name	Type	Description	Default
xgeo_var	Variable	The variable containing satellite position data in GEO coordinates.	required
time_var	Variable	The variable containing the timestamps.	required
pa_local_var	Variable	The variable containing the local pitch angle data.	required
irbem_input	IrbemInput	A data class with all required IRBEM input parameters.	required

Returns:

Type	Description
dict[str, Variable]	dict[str, ep.Variable]: A dictionary containing the calculated Lm, Lstar, and XJ variables.

Source code in el_paso/processing/magnetic_field_utils/magnetic_field_functions.py

@timed_function()
def get_Lstar(
    xgeo_var: ep.Variable, time_var: ep.Variable, pa_local_var: ep.Variable, irbem_input: IrbemInput
) -> dict[str, ep.Variable]:
    """Calculates Lm, Lstar, and the third adiabatic invariant (J).

    This function computes Lm and Lstar for each satellite position and local
    pitch angle. These are crucial parameters for characterizing particle drift shells
    in the magnetosphere. The calculation is parallelized for performance.

    Args:
        xgeo_var (ep.Variable): The variable containing satellite position data in GEO coordinates.
        time_var (ep.Variable): The variable containing the timestamps.
        pa_local_var (ep.Variable): The variable containing the local pitch angle data.
        irbem_input (IrbemInput): A data class with all required IRBEM input parameters.

    Returns:
        dict[str, ep.Variable]: A dictionary containing the calculated `Lm`, `Lstar`, and `XJ` variables.
    """
    logger.info("\tCalculating Lstar and J ...")

    timestamps = time_var.get_data(ep.units.posixtime)
    x_geo = xgeo_var.get_data(ep.units.RE)
    pa_local = pa_local_var.get_data(u.deg)

    datetimes = [datetime.fromtimestamp(t, tz=timezone.utc) for t in timestamps]
    sysaxes = ep.IRBEM_SYSAXIS_GEO

    x_geo = x_geo.astype(np.float64)
    pa_local = pa_local.astype(np.float64)
    irbem_input.maginput = {key: arr.astype(np.float64) for key, arr in irbem_input.maginput.items()}

    if len(datetimes) != len(irbem_input.maginput["Kp"]):
        msg = (
            f"Encountered size mismatch for Kp: len of Kp data: {len(irbem_input.maginput['Kp'])}, "
            f"requested len: {len(datetimes)}"
        )
        raise ValueError(msg)
    if len(datetimes) != len(x_geo):
        msg = f"Encountered size mismatch for x_geo: len of x_geo data: {len(x_geo)}, requested len: {len(datetimes)}"
        raise ValueError(msg)
    if len(datetimes) != len(pa_local):
        msg = (
            f"Encountered size mismatch for pa_local: len of pa_local data: {len(pa_local)}, "
            f"requested len: {len(datetimes)}"
        )
        raise ValueError(msg)

    kext = irbem_input.magnetic_field.kext()

    irbem_args = (irbem_input.irbem_lib_path, irbem_input.irbem_options, kext, sysaxes)

    parallel_func = partial(
        _make_lstar_shell_splitting_parallel,
        irbem_args,
        x_geo,
        datetimes,
        irbem_input.maginput,
        pa_local,
    )

    with Pool(processes=irbem_input.num_cores) as pool:
        chunksize = max(1, len(datetimes) // irbem_input.num_cores // 4)  # same as default
        rs = pool.map_async(parallel_func, range(len(datetimes)), chunksize=chunksize)
        show_process_bar_for_map_async(rs, chunksize)

    # write async results into one array
    Lm = np.empty_like(pa_local)
    Lstar = np.empty_like(pa_local)
    xj = np.empty_like(pa_local)

    results = rs.get()

    for i in range(len(datetimes)):
        Lm[i, :] = results[i][0]
        Lstar[i, :] = results[i][1]
        xj[i, :] = results[i][2]

    # replace bad values with nan
    for arr in [Lm, Lstar, xj]:
        arr[arr < 0] = np.nan
        if not np.any(np.isfinite(arr)) and irbem_input.irbem_options[0] != 0:
            msg = (
                "Lstar calculation failed! All points are NaNs! Hints for debugging:\n"
                "1) The calculation can fail for very low pitch-angles, where particles"
                "are not actually trapped\n"
                "2) Make sure your equatorial pitch angles and xGEO positions are correct\n"
                "3) Check other magnetic field outputs like equatorial magnetic fields."
                "If they are also NaN, the maginput to IRBEM might be wrong and needs debugging."
            )
            raise ValueError(msg)

    Lm_var = ep.Variable(data=Lm.astype(np.float64), original_unit=u.dimensionless_unscaled)
    Lm_var.metadata.add_processing_note(
        f"Calculated Lm using IRBEM model {irbem_input.magnetic_field} with options {irbem_input.irbem_options}."
    )

    Lstar_var = ep.Variable(data=Lstar.astype(np.float64), original_unit=u.dimensionless_unscaled)
    Lstar_var.metadata.add_processing_note(
        f"Calculated Lstar using IRBEM model {irbem_input.magnetic_field} with options {irbem_input.irbem_options}."
    )

    XJ_var = ep.Variable(data=xj.astype(np.float64), original_unit=ep.units.RE)
    XJ_var.metadata.add_processing_note(
        f"Calculated XJ using IRBEM model {irbem_input.magnetic_field} with options {irbem_input.irbem_options}."
    )

    return {
        create_var_name("L_m", irbem_input.magnetic_field): Lm_var,
        create_var_name("L_star", irbem_input.magnetic_field): Lstar_var,
        create_var_name("I", irbem_input.magnetic_field): XJ_var,
    }

el_paso.processing.magnetic_field_utils.get_MLT

get_MLT

Calculates the magnetic local time (MLT).

Parameters:

Name	Type	Description	Default
xgeo_var	Variable	The variable containing satellite position data in GEO coordinates.	required
time_var	Variable	The variable containing the timestamps.	required
irbem_input	IrbemInput	A data class with all required IRBEM input parameters.	required

Returns:

Type	Description
dict[str, Variable]	dict[str, ep.Variable]: A dictionary containing the calculated MLT variable.

Source code in el_paso/processing/magnetic_field_utils/magnetic_field_functions.py

@timed_function()
def get_MLT(xgeo_var: ep.Variable, time_var: ep.Variable, irbem_input: IrbemInput) -> dict[str, ep.Variable]:
    """Calculates the magnetic local time (MLT).

    Args:
        xgeo_var (ep.Variable): The variable containing satellite position data in GEO coordinates.
        time_var (ep.Variable): The variable containing the timestamps.
        irbem_input (IrbemInput): A data class with all required IRBEM input parameters.

    Returns:
        dict[str, ep.Variable]: A dictionary containing the calculated `MLT` variable.
    """
    logger.info("\tCalculating magnetic local time ...")

    timestamps = time_var.get_data(ep.units.posixtime)
    x_geo = xgeo_var.get_data(ep.units.RE)

    datetimes = [datetime.fromtimestamp(t, tz=timezone.utc) for t in timestamps]
    sysaxes = ep.IRBEM_SYSAXIS_GEO

    # Ensure xGEO and maginput are floating-point arrays
    x_geo = x_geo.astype(np.float64)

    if len(datetimes) != len(x_geo):
        msg = f"Encountered size mismatch for x_geo: len of x_geo data: {len(x_geo)}, requested len: {len(datetimes)}"
        raise ValueError(msg)

    kext = irbem_input.magnetic_field.kext()

    model = MagFields(
        lib_path=irbem_input.irbem_lib_path,
        options=irbem_input.irbem_options,
        kext=kext,
        sysaxes=sysaxes,
    )

    mlt_output = np.empty_like(datetimes)

    for i in range(len(datetimes)):
        x_dict: dict[Literal["x1", "x2", "x3"], np.floating] = {
            "x1": x_geo[i, 0],
            "x2": x_geo[i, 1],
            "x3": x_geo[i, 2],
        }
        mlt_output[i] = model.get_mlt(datetimes[i], x_dict)

    mlt_output = mlt_output.astype(np.float64)

    var = ep.Variable(data=mlt_output, original_unit=u.hour)
    var.metadata.add_processing_note(
        f"Calculated MLT using IRBEM model {irbem_input.magnetic_field} with options {irbem_input.irbem_options}."
    )

    return {create_var_name("MLT", irbem_input.magnetic_field): var}

el_paso.processing.magnetic_field_utils.get_footpoint_atmosphere

get_footpoint_atmosphere

Calculates the magnetic field strength at the atmospheric foot point.

This function uses parallel processing to calculate the magnetic field strength at the atmospheric foot point (100 km altitude) for each satellite position. It returns the result as a dictionary of el_paso.Variable.

Parameters:

Name	Type	Description	Default
xgeo_var	Variable	The variable containing satellite position data in GEO coordinates.	required
time_var	Variable	The variable containing the timestamps.	required
irbem_input	IrbemInput	A data class with all required IRBEM input parameters.	required

Returns:

Type	Description
dict[str, Variable]	dict[str, ep.Variable]: A dictionary containing the calculated B_fofl variable.

Source code in el_paso/processing/magnetic_field_utils/magnetic_field_functions.py

@timed_function()
def get_footpoint_atmosphere(
    xgeo_var: ep.Variable, time_var: ep.Variable, irbem_input: IrbemInput
) -> dict[str, ep.Variable]:
    """Calculates the magnetic field strength at the atmospheric foot point.

    This function uses parallel processing to calculate the magnetic field strength
    at the atmospheric foot point (100 km altitude) for each satellite position.
    It returns the result as a dictionary of `el_paso.Variable`.

    Args:
        xgeo_var (ep.Variable): The variable containing satellite position data in GEO coordinates.
        time_var (ep.Variable): The variable containing the timestamps.
        irbem_input (IrbemInput): A data class with all required IRBEM input parameters.

    Returns:
        dict[str, ep.Variable]: A dictionary containing the calculated `B_fofl` variable.
    """
    logger.info("\tCalculating magnetic foot point at the atmosphere ...")

    timestamps = time_var.get_data(ep.units.posixtime)
    x_geo = xgeo_var.get_data(ep.units.RE)

    datetimes = [datetime.fromtimestamp(t, tz=timezone.utc) for t in timestamps]
    sysaxes = ep.IRBEM_SYSAXIS_GEO

    x_geo = x_geo.astype(np.float64)

    if len(datetimes) != len(x_geo):
        msg = f"Encountered size mismatch for x_geo: len of x_geo data: {len(x_geo)}, requested len: {len(datetimes)}"
        raise ValueError(msg)

    kext = irbem_input.magnetic_field.kext()

    irbem_args = (irbem_input.irbem_lib_path, irbem_input.irbem_options, kext, sysaxes)

    parallel_func = partial(_get_footpoint_atmosphere_parallel, irbem_args, x_geo, datetimes, irbem_input.maginput)

    with Pool(processes=irbem_input.num_cores) as pool:
        chunksize = max(1, len(datetimes) // irbem_input.num_cores // 4)  # same as default
        rs = pool.map_async(parallel_func, range(len(datetimes)), chunksize=chunksize)
        show_process_bar_for_map_async(rs, chunksize)

    # write async results into one array
    B_foot = np.empty_like(datetimes)

    results = rs.get()

    for i in range(len(datetimes)):
        B_foot[i] = results[i]

    B_foot[B_foot == FORTRAN_BAD_VALUE] = np.nan

    var = ep.Variable(data=B_foot.astype(np.float64), original_unit=u.nT)
    var.metadata.add_processing_note(
        f"Calculated foot point at the atmosphere using IRBEM model {irbem_input.magnetic_field} "
        f"with options {irbem_input.irbem_options}."
    )

    return {create_var_name("B_fofl", irbem_input.magnetic_field): var}

el_paso.processing.magnetic_field_utils.get_local_B_field

get_local_B_field

Calculates the local magnetic field magnitude.

Parameters:

Name	Type	Description	Default
xgeo_var	Variable	The variable containing satellite position data in GEO coordinates.	required
time_var	Variable	The variable containing the timestamps.	required
irbem_input	IrbemInput	A data class with all required IRBEM input parameters.	required

Returns:

Type	Description
dict[str, Variable]	dict[str, ep.Variable]: A dictionary containing the calculated B_local variable.

Source code in el_paso/processing/magnetic_field_utils/magnetic_field_functions.py

@timed_function()
def get_local_B_field(xgeo_var: ep.Variable, time_var: ep.Variable, irbem_input: IrbemInput) -> dict[str, ep.Variable]:
    """Calculates the local magnetic field magnitude.

    Args:
        xgeo_var (ep.Variable): The variable containing satellite position data in GEO coordinates.
        time_var (ep.Variable): The variable containing the timestamps.
        irbem_input (IrbemInput): A data class with all required IRBEM input parameters.

    Returns:
        dict[str, ep.Variable]: A dictionary containing the calculated `B_local` variable.
    """
    logger.info("\tCalculating local magnetic field values ...")

    timestamps = time_var.get_data(ep.units.posixtime)
    x_geo = xgeo_var.get_data(ep.units.RE)

    datetimes = [datetime.fromtimestamp(t, tz=timezone.utc) for t in timestamps]
    sysaxes = ep.IRBEM_SYSAXIS_GEO

    # Define Fortran bad value as a float
    fortran_bad_value = np.float64(-1.0e31)
    # Ensure x_geo and maginput are floating-point arrays
    x_geo = x_geo.astype(np.float64)
    for key in irbem_input.maginput:
        irbem_input.maginput[key] = np.array(irbem_input.maginput[key], dtype=np.float64)

    if len(datetimes) != len(irbem_input.maginput["Kp"]):
        msg = (
            f"Encountered size mismatch for Kp: len of Kp data: {len(irbem_input.maginput['Kp'])}, "
            f"requested len: {len(datetimes)}"
        )
        raise ValueError(msg)
    if len(datetimes) != len(x_geo):
        msg = f"Encountered size mismatch for x_geo: len of x_geo data: {len(x_geo)}, requested len: {len(datetimes)}"
        raise ValueError(msg)

    x_dict: dict[Literal["x1", "x2", "x3"], NDArray[np.floating]] = {
        "x1": x_geo[:, 0],
        "x2": x_geo[:, 1],
        "x3": x_geo[:, 2],
    }
    kext = irbem_input.magnetic_field.kext()

    model = MagFields(
        lib_path=irbem_input.irbem_lib_path,
        options=irbem_input.irbem_options,
        kext=kext,
        sysaxes=sysaxes,
    )

    field_multi_output = model.get_field_multi(datetimes, x_dict, irbem_input.maginput)

    # replace bad values with nan
    field_multi_output.bgeo[field_multi_output.bgeo == fortran_bad_value] = np.nan
    field_multi_output.blocal[field_multi_output.blocal == fortran_bad_value] = np.nan

    b_local_var = ep.Variable(data=field_multi_output.blocal, original_unit=u.nT)
    return {create_var_name("B_Calc", irbem_input.magnetic_field): b_local_var}

el_paso.processing.magnetic_field_utils.get_magequator

get_magequator

Calculates the magnetic field strength and radial distance at the magnetic equator.

This function uses parallel processing to efficiently compute magnetic field and position properties at the magnetic equator for a given set of satellite positions over time. It returns the results as a dictionary of el_paso.Variable objects.

Parameters:

Name	Type	Description	Default
xgeo_var	Variable	The variable containing satellite position data in GEO coordinates.	required
time_var	Variable	The variable containing the timestamps for the data.	required
irbem_input	IrbemInput	A data class with all required IRBEM input parameters.	required

Returns:

Type	Description
dict[str, Variable]	dict[str, ep.Variable]: A dictionary containing the calculated B_eq, R_eq, and xGEO_eq variables.

Source code in el_paso/processing/magnetic_field_utils/magnetic_field_functions.py

@timed_function()
def get_magequator(xgeo_var: ep.Variable, time_var: ep.Variable, irbem_input: IrbemInput) -> dict[str, ep.Variable]:
    """Calculates the magnetic field strength and radial distance at the magnetic equator.

    This function uses parallel processing to efficiently compute magnetic field
    and position properties at the magnetic equator for a given set of satellite
    positions over time. It returns the results as a dictionary of `el_paso.Variable` objects.

    Args:
        xgeo_var (ep.Variable): The variable containing satellite position data in GEO coordinates.
        time_var (ep.Variable): The variable containing the timestamps for the data.
        irbem_input (IrbemInput): A data class with all required IRBEM input parameters.

    Returns:
        dict[str, ep.Variable]: A dictionary containing the calculated `B_eq`, `R_eq`, and
                                `xGEO_eq` variables.
    """
    logger.info("\tCalculating magnetic field and radial distance at the equator ...")

    timestamps = time_var.get_data(ep.units.posixtime)
    x_geo = xgeo_var.get_data(ep.units.RE)

    datetimes = [datetime.fromtimestamp(t, tz=timezone.utc) for t in timestamps]
    sysaxes = ep.IRBEM_SYSAXIS_GEO

    x_geo = x_geo.astype(np.float64)

    if len(datetimes) != len(x_geo):
        msg = f"Encountered size mismatch for x_geo: len of x_geo data: {len(x_geo)}, requested len: {len(datetimes)}"
        raise ValueError(msg)

    kext = irbem_input.magnetic_field.kext()

    irbem_args = (irbem_input.irbem_lib_path, irbem_input.irbem_options, kext, sysaxes)

    parallel_func = partial(_get_magequator_parallel, irbem_args, x_geo, datetimes, irbem_input.maginput)

    with Pool(processes=irbem_input.num_cores) as pool:
        chunksize = max(1, len(datetimes) // irbem_input.num_cores // 4)  # same as default
        rs = pool.map_async(parallel_func, range(len(datetimes)), chunksize=chunksize)
        show_process_bar_for_map_async(rs, chunksize)

    # write async results into one array
    B_eq = np.empty_like(datetimes)
    x_geo_min = np.empty_like(x_geo)

    results = rs.get()

    for i in range(len(datetimes)):
        B_eq[i] = results[i][0]
        x_geo_min[i] = results[i][1]

    B_eq[B_eq == FORTRAN_BAD_VALUE] = np.nan
    x_geo_min[x_geo_min == FORTRAN_BAD_VALUE] = np.nan

    B_eq_var = ep.Variable(data=B_eq.astype(np.float64), original_unit=u.nT)
    B_eq_var.metadata.add_processing_note(
        f"Calculated magnetic field at the equator using IRBEM model {irbem_input.magnetic_field} "
        f"with options {irbem_input.irbem_options}."
    )

    x_geo_var = ep.Variable(data=x_geo_min.astype(np.float64), original_unit=ep.units.RE)
    x_geo_var.metadata.add_processing_note(
        f"Calculated radial distance at the equator using IRBEM model {irbem_input.magnetic_field} "
        f"with options {irbem_input.irbem_options}."
    )

    # add radial distance field in SM coordinates
    x_gsm = Coords(lib_path=irbem_input.irbem_lib_path).transform(
        datetimes,
        x_geo_min,
        ep.IRBEM_SYSAXIS_GEO,
        ep.IRBEM_SYSAXIS_GSM,
    )

    R_eq_var = ep.Variable(
        data=np.linalg.norm(x_gsm, ord=2, axis=1).astype(np.float64),
        original_unit=ep.units.RE,
    )
    R_eq_var.metadata.add_processing_note(
        f"Calculated radial distance at the equator in GSM coordinates using IRBEM model {irbem_input.magnetic_field} "
        f"with options {irbem_input.irbem_options}."
    )

    p_gsm = np.arctan2(x_gsm[:, 1], x_gsm[:, 0])
    mlt_gsm = ((p_gsm * 12 / np.pi) + 12) % 24

    mlt_eq_var = ep.Variable(
        data=mlt_gsm.astype(np.float64),
        original_unit=u.hour,
    )
    mlt_eq_var.metadata.add_processing_note(
        "Calculated magnetic local time at the equator in GSM coordinates using "
        f"IRBEM model {irbem_input.magnetic_field} with options {irbem_input.irbem_options}."
    )

    return {
        create_var_name("B_Eq", irbem_input.magnetic_field): B_eq_var,
        create_var_name("R_Eq", irbem_input.magnetic_field): R_eq_var,
        create_var_name("MLT_Eq", irbem_input.magnetic_field): mlt_eq_var,
        create_var_name("xGEO_Eq", irbem_input.magnetic_field): x_geo_var,
    }

el_paso.processing.magnetic_field_utils.get_mirror_point

get_mirror_point

Calculates the magnetic field strength at the mirror point.

This function computes the magnetic field strength at the mirror point for a given set of local pitch angles.

Parameters:

Name	Type	Description	Default
xgeo_var	Variable	The variable containing satellite position data in GEO coordinates.	required
time_var	Variable	The variable containing the timestamps.	required
pa_local_var	Variable	The variable containing the local pitch angle data.	required
irbem_input	IrbemInput	A data class with all required IRBEM input parameters.	required

Returns:

Type	Description
dict[str, Variable]	dict[str, ep.Variable]: A dictionary containing the calculated B_mirr variable.

Source code in el_paso/processing/magnetic_field_utils/magnetic_field_functions.py

@timed_function()
def get_mirror_point(
    xgeo_var: ep.Variable, time_var: ep.Variable, pa_local_var: ep.Variable, irbem_input: IrbemInput
) -> dict[str, ep.Variable]:
    """Calculates the magnetic field strength at the mirror point.

    This function computes the magnetic field strength at the mirror point for a
    given set of local pitch angles.

    Args:
        xgeo_var (ep.Variable): The variable containing satellite position data in GEO coordinates.
        time_var (ep.Variable): The variable containing the timestamps.
        pa_local_var (ep.Variable): The variable containing the local pitch angle data.
        irbem_input (IrbemInput): A data class with all required IRBEM input parameters.

    Returns:
        dict[str, ep.Variable]: A dictionary containing the calculated `B_mirr` variable.
    """
    logger.info("\tCalculating mirror points ...")

    timestamps = time_var.get_data(ep.units.posixtime)
    x_geo = xgeo_var.get_data(ep.units.RE)
    pa_local = pa_local_var.get_data(u.deg)

    datetimes = [datetime.fromtimestamp(t, tz=timezone.utc) for t in timestamps]
    sysaxes = ep.IRBEM_SYSAXIS_GEO

    x_geo = x_geo.astype(np.float64)
    pa_local = pa_local.astype(np.float64)
    irbem_input.maginput = {key: arr.astype(np.float64) for key, arr in irbem_input.maginput.items()}

    if len(datetimes) != len(irbem_input.maginput["Kp"]):
        msg = (
            f"Encountered size mismatch for Kp: len of Kp data: {len(irbem_input.maginput['Kp'])}, "
            f"requested len: {len(datetimes)}"
        )
        raise ValueError(msg)
    if len(datetimes) != len(x_geo):
        msg = f"Encountered size mismatch for x_geo: len of x_geo data: {len(x_geo)}, requested len: {len(datetimes)}"
        raise ValueError(msg)
    if len(datetimes) != len(pa_local):
        msg = (
            f"Encountered size mismatch for pa_local: len of pa_local data: {len(pa_local)}, "
            f"requested len: {len(datetimes)}"
        )
        raise ValueError(msg)

    kext = irbem_input.magnetic_field.kext()

    irbem_args = (irbem_input.irbem_lib_path, irbem_input.irbem_options, kext, sysaxes)

    parallel_func = partial(_get_mirror_point_parallel, irbem_args, x_geo, datetimes, irbem_input.maginput, pa_local)

    with Pool(processes=irbem_input.num_cores) as pool:
        chunksize = max(1, len(datetimes) // irbem_input.num_cores // 4)  # same as default
        rs = pool.map_async(parallel_func, range(len(datetimes)), chunksize=chunksize)
        show_process_bar_for_map_async(rs, chunksize)

    # write async results into one array
    mirror_point_output = np.empty_like(pa_local)

    results = rs.get()

    for i in range(len(datetimes)):
        mirror_point_output[i, :] = results[i]

    # replace bad values with nan
    mirror_point_output[mirror_point_output < 0] = np.nan

    var = ep.Variable(data=mirror_point_output.astype(np.float64), original_unit=u.nT)
    var.metadata.add_processing_note(
        f"Calculated mirror points using IRBEM model {irbem_input.magnetic_field} "
        f"with options {irbem_input.irbem_options}."
    )

    return {create_var_name("B_mirr", irbem_input.magnetic_field): var}