diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 2195e6f1..3182f9e2 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -17,6 +17,7 @@ on:
 env:
   CI: true
   UV_SYSTEM_PYTHON: 1
+  COVERAGE_CORE: sysmon
 
 jobs:
   build:
diff --git a/doc/Changelog.md b/doc/Changelog.md
index 3a99d64b..362c4a81 100644
--- a/doc/Changelog.md
+++ b/doc/Changelog.md
@@ -19,6 +19,12 @@ og:description: See what's new in the latest release of Roseau Load Flow !
 
 ## Unreleased
 
+- {gh-pr}`304` Add top-level modules `rlf.constants` and `rlf.types`. The old modules in the `utils`
+  package are deprecated and will be removed in a future release. The `utils` package is for internal
+  use only and should not be considered stable.
+- {gh-pr}`304` Add top-level module `rlf.sym` for symmetrical components utilities. The `sym_to_phasor`,
+  `phasor_to_sym` and `series_phasor_to_sym` functions are moved from the `rlf.converters` module to
+  this module. The old functions are deprecated and will be removed in a future release.
 - {gh-pr}`303` Fix missing `voltage_level` in `en.res_buses_voltages` when the buses define nominal
   voltage but not voltage limits.
 - {gh-pr}`303` Add `rlf.SQRT3` constant for the square root of 3. It can be useful for the conversion
diff --git a/doc/conf.py b/doc/conf.py
index d6a0fc73..9d769891 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -132,7 +132,16 @@
 
 # -- Options for AutoAPI -------------------------------------------------
 autoapi_dirs = ["../roseau"]
-autoapi_ignore = ["**/tests/**", "**/conftest.py", "__about__.py"]
+autoapi_ignore = [
+    # Tests
+    "**/tests/**",
+    "**/conftest.py",
+    # Cruft
+    "**/roseau/load_flow/__about__.py",
+    # Internal utilities
+    "**/roseau/load_flow/utils/**",
+    "**/roseau/load_flow/io/**",
+]
 autoapi_options = ["members", "show-inheritance", "show-module-summary", "imported-members"]
 autoapi_python_class_content = "both"  # without this, the __init__ docstring is not shown
 autoapi_python_use_implicit_namespaces = True
diff --git a/doc/usage/Extras.md b/doc/usage/Extras.md
index adf822df..71e67b63 100644
--- a/doc/usage/Extras.md
+++ b/doc/usage/Extras.md
@@ -48,18 +48,18 @@ If we take the example network from the [Getting Started page](gs-creating-netwo
 
 As there are no transformers between the two buses, they all belong to the same cluster.
 
-## Conversion to symmetrical components
+## Symmetrical components
 
-{mod}`roseau.load_flow.converters` contains helpers to convert between phasor and symmetrical
-components. For example, to convert a phasor voltage to symmetrical components:
+{mod}`roseau.load_flow.sym` contains helpers to work with symmetrical components. For example, to
+convert a phasor voltage to symmetrical components:
 
 ```pycon
 >>> import numpy as np
->>> from roseau.load_flow.converters import phasor_to_sym, sym_to_phasor
+>>> import roseau.load_flow as rlf
 >>> v = 230 * np.exp([0, -2j * np.pi / 3, 2j * np.pi / 3])
 >>> v
 array([ 230.  +0.j        , -115.-199.18584287j, -115.+199.18584287j])
->>> v_sym = phasor_to_sym(v)
+>>> v_sym = rlf.sym.phasor_to_sym(v)
 >>> v_sym
 array([[ 8.52651283e-14-1.42108547e-14j],
        [ 2.30000000e+02+4.19109192e-14j],
@@ -80,8 +80,7 @@ You can also convert pandas Series to symmetrical components. If we take the exa
 [Getting Started](Getting_Started.md) page:
 
 ```pycon
->>> from roseau.load_flow.converters import series_phasor_to_sym
->>> series_phasor_to_sym(en.res_buses_voltages["voltage"])
+>>> rlf.sym.series_phasor_to_sym(en.res_buses_voltages["voltage"])
 bus_id  sequence
 lb      zero        8.526513e-14-1.421085e-14j
         pos         2.219282e+02+4.167975e-14j
@@ -92,14 +91,14 @@ sb      zero        9.947598e-14-1.421085e-14j
 Name: voltage, dtype: complex128
 ```
 
-_Roseau Load Flow_ also provides useful helpers to create three-phase balanced quantities by only
+The `rlf.sym` module also provides useful helpers to create three-phase balanced quantities by only
 providing the magnitude of the quantities. For example, to create a three-phase balanced positive
 sequence voltage:
 
 ```pycon
 >>> import numpy as np
 >>> import roseau.load_flow as rlf
->>> V = 230 * rlf.PositiveSequence
+>>> V = 230 * rlf.sym.PositiveSequence
 >>> V
 array([ 230.  +0.j        , -115.-199.18584287j, -115.+199.18584287j])
 >>> np.abs(V)
@@ -108,23 +107,24 @@ array([230., 230., 230.])
 array([   0., -120.,  120.])
 ```
 
-Similarly, you can use `rlf.NegativeSequence` and `rlf.ZeroSequence` to create negative-sequence
-and zero-sequence quantities respectively.
+Similarly, you can use `rlf.sym.NegativeSequence` and `rlf.sym.ZeroSequence` to create negative-sequence
+and zero-sequence quantities respectively. Because these are so common, you can also access them
+directly from the top-level module as `rlf.PositiveSequence`, etc.
 
 ## Potentials to voltages conversion
 
-{mod}`roseau.load_flow.converters` also contains helpers to convert a vector of potentials to a
-vector of voltages. Example:
+{mod}`roseau.load_flow.converters` contains helpers to convert a vector of potentials to a vector of
+voltages. Example:
 
 ```pycon
 >>> import numpy as np
->>> from roseau.load_flow.converters import calculate_voltages, calculate_voltage_phases
+>>> import roseau.load_flow as rlf
 >>> potentials = 230 * np.array([1, np.exp(-2j * np.pi / 3), np.exp(2j * np.pi / 3), 0])
 >>> potentials
 array([ 230.  +0.j        , -115.-199.18584287j, -115.+199.18584287j,
           0.  +0.j        ])
 >>> phases = "abcn"
->>> calculate_voltages(potentials, phases)
+>>> rlf.converters.calculate_voltages(potentials, phases)
 array([ 230.  +0.j        , -115.-199.18584287j, -115.+199.18584287j]) <Unit('volt')>
 ```
 
@@ -132,40 +132,46 @@ Because the phases include the neutral, the voltages calculated are phase-to-neu
 You can also calculate phase-to-phase voltages by omitting the neutral:
 
 ```pycon
->>> calculate_voltages(potentials[:-1], phases[:-1])
+>>> rlf.converters.calculate_voltages(potentials[:-1], phases[:-1])
 array([ 345.+199.18584287j,    0.-398.37168574j, -345.+199.18584287j]) <Unit('volt')>
 ```
 
 To get the phases of the voltage, you can use `calculate_voltage_phases`:
 
 ```pycon
->>> calculate_voltage_phases(phases)
+>>> rlf.converters.calculate_voltage_phases(phases)
 ['an', 'bn', 'cn']
 ```
 
 Of course these functions work with arbitrary phases:
 
 ```pycon
->>> calculate_voltages(potentials[:2], phases[:2])
+>>> rlf.converters.calculate_voltages(potentials[:2], phases[:2])
 array([345.+199.18584287j]) <Unit('volt')>
->>> calculate_voltage_phases(phases[:2])
+>>> rlf.converters.calculate_voltage_phases(phases[:2])
 ['ab']
->>> calculate_voltage_phases("abc")
+>>> rlf.converters.calculate_voltage_phases("abc")
 ['ab', 'bc', 'ca']
->>> calculate_voltage_phases("bc")
+>>> rlf.converters.calculate_voltage_phases("bc")
 ['bc']
->>> calculate_voltage_phases("bcn")
+>>> rlf.converters.calculate_voltage_phases("bcn")
 ['bn', 'cn']
 ```
 
 ## Constants
 
-{mod}`roseau.load_flow.utils.constants` contains some common constants like the resistivity
-and permeability of common materials in addition to other useful constants. Please refer to
-the module documentation for more details.
+{mod}`roseau.load_flow.constants` contains some common mathematical and physical constants like the
+resistivity and permeability of common materials in addition to other useful constants. Please refer
+to the module documentation for more details. An enumeration of available materials can be found in
+the {mod}`roseau.load_flow.types` module.
+
+Some commonly used constants can be accessed directly from the top-level module for convenience.
+Notable top-level constants:
 
-An enumeration of available materials can be found in the {mod}`roseau.load_flow.utils.types`
-module.
+- `rlf.SQRT3`: the square root of 3. Useful for converting between phase-to-phase and phase-to-neutral
+  voltages.
+- `rlf.ALPHA`: the alpha constant. Rotates a complex number by 120°.
+- `rlf.ALPHA2`: the alpha constant squared. Rotates a complex number by 240° (or -120°).
 
 ## Voltage unbalance
 
diff --git a/pyproject.toml b/pyproject.toml
index 1b65adfb..449a3bc0 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -110,7 +110,7 @@ docstring-code-format = true
 "*.ipynb" = ["E402", "F403", "F405"]
 
 [tool.coverage.run]
-branch = true
+branch = false
 omit = ["roseau/load_flow/__about__.py"]
 plugins = ["coverage_conditional_plugin"]
 
diff --git a/roseau/load_flow/__init__.py b/roseau/load_flow/__init__.py
index 92fd0467..fe080fd2 100644
--- a/roseau/load_flow/__init__.py
+++ b/roseau/load_flow/__init__.py
@@ -9,7 +9,7 @@
 import importlib.metadata
 from typing import Any
 
-from roseau.load_flow import converters
+from roseau.load_flow import constants, converters, plotting, sym
 from roseau.load_flow.__about__ import (
     __authors__,
     __copyright__,
@@ -20,6 +20,7 @@
     __status__,
     __url__,
 )
+from roseau.load_flow.constants import SQRT3
 from roseau.load_flow.exceptions import RoseauLoadFlowException, RoseauLoadFlowExceptionCode
 from roseau.load_flow.license import License, activate_license, deactivate_license, get_license
 from roseau.load_flow.models import (
@@ -43,10 +44,10 @@
     VoltageSource,
 )
 from roseau.load_flow.network import ElectricalNetwork
+from roseau.load_flow.sym import ALPHA, ALPHA2, NegativeSequence, PositiveSequence, ZeroSequence
+from roseau.load_flow.types import Insulator, LineType, Material
 from roseau.load_flow.units import Q_, ureg
-from roseau.load_flow.utils import Insulator, LineType, Material, constants
-from roseau.load_flow.utils._versions import show_versions
-from roseau.load_flow.utils.constants import ALPHA, ALPHA2, SQRT3, NegativeSequence, PositiveSequence, ZeroSequence
+from roseau.load_flow.utils import show_versions
 
 __version__ = importlib.metadata.version("roseau-load-flow")
 
@@ -63,6 +64,8 @@
     "show_versions",
     "converters",
     "constants",
+    "plotting",
+    "sym",
     # Electrical Network
     "ElectricalNetwork",
     # Buses
@@ -119,12 +122,12 @@ def __getattr__(name: str) -> Any:
     if name in deprecated_classes and name not in globals():
         import warnings
 
-        from roseau.load_flow.utils._exceptions import find_stack_level
+        from roseau.load_flow.utils.exceptions import find_stack_level
 
         new_class = deprecated_classes[name]
         warnings.warn(
             f"The `{name}` class is deprecated. Use `{new_class.__name__}` instead.",
-            category=DeprecationWarning,
+            category=FutureWarning,
             stacklevel=find_stack_level(),
         )
         globals()[name] = new_class
diff --git a/roseau/load_flow/_compat.py b/roseau/load_flow/_compat.py
index 4f571cfe..2f46cd8d 100644
--- a/roseau/load_flow/_compat.py
+++ b/roseau/load_flow/_compat.py
@@ -1,3 +1,4 @@
+# Cannot move to utils because of circular imports
 import sys
 from enum import Enum
 
diff --git a/roseau/load_flow/_solvers.py b/roseau/load_flow/_solvers.py
index 755c046a..73f77ed2 100644
--- a/roseau/load_flow/_solvers.py
+++ b/roseau/load_flow/_solvers.py
@@ -9,7 +9,7 @@
 from roseau.load_flow.exceptions import RoseauLoadFlowException, RoseauLoadFlowExceptionCode
 from roseau.load_flow.license import activate_license, get_license
 from roseau.load_flow.typing import JsonDict, Solver
-from roseau.load_flow.utils._exceptions import find_stack_level
+from roseau.load_flow.utils import find_stack_level
 from roseau.load_flow_engine.cy_engine import CyAbstractSolver, CyBackwardForward, CyNewton, CyNewtonGoldstein
 
 logger = logging.getLogger(__name__)
diff --git a/roseau/load_flow/_wrapper.py b/roseau/load_flow/_wrapper.py
deleted file mode 100644
index 59158479..00000000
--- a/roseau/load_flow/_wrapper.py
+++ /dev/null
@@ -1,148 +0,0 @@
-import functools
-from collections.abc import Callable, Iterable, MutableSequence
-from inspect import Parameter, Signature, signature
-from itertools import zip_longest
-from typing import Any, TypeVar
-
-from pint import Quantity, Unit
-from pint.registry import UnitRegistry
-from pint.util import to_units_container
-
-T = TypeVar("T")
-FuncT = TypeVar("FuncT", bound=Callable)
-
-
-def _parse_wrap_args(args: Iterable[str | Unit | None]) -> Callable:
-    """Create a converter function for the wrapper"""
-    # _to_units_container
-    args_as_uc = [to_units_container(arg) for arg in args]
-
-    # Check for references in args, remove None values
-    unit_args_ndx = {ndx for ndx, arg in enumerate(args_as_uc) if arg is not None}
-
-    def _converter(ureg: "UnitRegistry", sig: "Signature", values: "list[Any]", kw: "dict[Any]"):
-        len_initial_values = len(values)
-
-        # pack kwargs
-        for i, param_name in enumerate(sig.parameters):
-            if i >= len_initial_values:
-                values.append(kw[param_name])
-
-        # convert arguments
-        for ndx in unit_args_ndx:
-            value = values[ndx]
-            if isinstance(value, ureg.Quantity):
-                values[ndx] = ureg.convert(value.magnitude, value.units, args_as_uc[ndx])
-            elif isinstance(value, MutableSequence):
-                for i, val in enumerate(value):
-                    if isinstance(val, ureg.Quantity):
-                        value[i] = ureg.convert(val.magnitude, val.units, args_as_uc[ndx])
-
-        # unpack kwargs
-        for i, param_name in enumerate(sig.parameters):
-            if i >= len_initial_values:
-                kw[param_name] = values[i]
-
-        return values[:len_initial_values], kw
-
-    return _converter
-
-
-def _apply_defaults(sig: Signature, args: tuple[Any], kwargs: dict[str, Any]) -> tuple[list[Any], dict[str, Any]]:
-    """Apply default keyword arguments.
-
-    Named keywords may have been left blank. This function applies the default
-    values so that every argument is defined.
-    """
-    n = len(args)
-    for i, param in enumerate(sig.parameters.values()):
-        if i >= n and param.default != Parameter.empty and param.name not in kwargs:
-            kwargs[param.name] = param.default
-    return list(args), kwargs
-
-
-def wraps(
-    ureg: UnitRegistry,
-    ret: str | Unit | Iterable[str | Unit | None] | None,
-    args: str | Unit | Iterable[str | Unit | None] | None,
-) -> Callable[[FuncT], FuncT]:
-    """Wraps a function to become pint-aware.
-
-    Use it when a function requires a numerical value but in some specific
-    units. The wrapper function will take a pint quantity, convert to the units
-    specified in `args` and then call the wrapped function with the resulting
-    magnitude.
-
-    The value returned by the wrapped function will be converted to the units
-    specified in `ret`.
-
-    Args:
-        ureg:
-            A UnitRegistry instance.
-
-        ret:
-            Units of each of the return values. Use `None` to skip argument conversion.
-
-        args:
-             Units of each of the input arguments. Use `None` to skip argument conversion.
-
-    Returns:
-        The wrapper function.
-
-    Raises:
-        TypeError
-            if the number of given arguments does not match the number of function parameters.
-            if any of the provided arguments is not a unit a string or Quantity
-    """
-    if not isinstance(args, (list, tuple)):
-        args = (args,)
-
-    for arg in args:
-        if arg is not None and not isinstance(arg, (ureg.Unit, str)):
-            raise TypeError(f"wraps arguments must by of type str or Unit, not {type(arg)} ({arg})")
-
-    converter = _parse_wrap_args(args)
-
-    is_ret_container = isinstance(ret, (list, tuple))
-    if is_ret_container:
-        for arg in ret:
-            if arg is not None and not isinstance(arg, (ureg.Unit, str)):
-                raise TypeError(f"wraps 'ret' argument must by of type str or Unit, not {type(arg)} ({arg})")
-        ret = ret.__class__([to_units_container(arg, ureg) for arg in ret])
-    else:
-        if ret is not None and not isinstance(ret, (ureg.Unit, str)):
-            raise TypeError(f"wraps 'ret' argument must by of type str or Unit, not {type(ret)} ({ret})")
-        ret = to_units_container(ret, ureg)
-
-    def decorator(func: "Callable[..., Any]") -> "Callable[..., Quantity]":
-        sig = signature(func)
-        count_params = len(sig.parameters)
-        if len(args) != count_params:
-            raise TypeError(f"{func.__name__} takes {count_params} parameters, but {len(args)} units were passed")
-
-        assigned = tuple(attr for attr in functools.WRAPPER_ASSIGNMENTS if hasattr(func, attr))
-        updated = tuple(attr for attr in functools.WRAPPER_UPDATES if hasattr(func, attr))
-
-        @functools.wraps(func, assigned=assigned, updated=updated)
-        def wrapper(*values, **kw) -> "Quantity":
-            values, kw = _apply_defaults(sig, values, kw)
-
-            # In principle, the values are used as is
-            # When then extract the magnitudes when needed.
-            new_values, new_kw = converter(ureg, sig, values, kw)
-
-            result = func(*new_values, **new_kw)
-
-            if is_ret_container:
-                return ret.__class__(
-                    res if unit is None else ureg.Quantity(res, unit) for unit, res in zip_longest(ret, result)
-                )
-
-            if ret is None:
-                return result
-
-            return ureg.Quantity(result, ret)
-
-        return wrapper
-
-    return decorator
diff --git a/roseau/load_flow/conftest.py b/roseau/load_flow/conftest.py
index 66d4dcf6..068b9e9a 100644
--- a/roseau/load_flow/conftest.py
+++ b/roseau/load_flow/conftest.py
@@ -6,7 +6,7 @@
 import pytest
 from _pytest.monkeypatch import MonkeyPatch
 
-import roseau
+import roseau.load_flow
 from roseau.load_flow.utils.log import set_logging_config
 
 HERE = Path(__file__).parent.expanduser().absolute()
@@ -63,6 +63,8 @@ def bar(*args, **kwargs):  # pragma: no-cover
                 for f in filenames:
                     if not f.endswith(".py"):
                         continue
+                    if f in ("constants.py", "types.py") and base_module == "roseau.load_flow.utils":
+                        continue  # TODO: Remove when deprecated modules are removed
                     module = importlib.import_module(f"{base_module}.{f.removesuffix('.py')}")
                     for _, klass in inspect.getmembers(
                         module,
diff --git a/roseau/load_flow/constants.py b/roseau/load_flow/constants.py
new file mode 100644
index 00000000..cff001cb
--- /dev/null
+++ b/roseau/load_flow/constants.py
@@ -0,0 +1,130 @@
+"""Physical and mathematical constants used in the load flow calculations."""
+
+import cmath
+import math
+from typing import Final
+
+from roseau.load_flow.types import Insulator, Material
+from roseau.load_flow.units import Q_
+
+__all__ = [
+    # Mathematical constants
+    "SQRT3",
+    "ALPHA",
+    "ALPHA2",
+    "PI",
+    # Physical constants
+    "MU_0",
+    "MU_R",
+    "EPSILON_0",
+    "EPSILON_R",
+    "F",
+    "OMEGA",
+    "RHO",
+    "TAN_D",
+    "DELTA_P",
+]
+
+SQRT3: Final = math.sqrt(3)
+"""The square root of 3."""
+
+ALPHA: Final = cmath.exp(2 / 3 * cmath.pi * 1j)
+"""complex: Phasor rotation operator :math:`\\alpha`, which rotates a phasor vector counterclockwise
+by 120 degrees when multiplied by it."""
+
+ALPHA2: Final = ALPHA**2
+"""complex: Phasor rotation operator :math:`\\alpha^2`, which rotates a phasor vector clockwise by
+120 degrees when multiplied by it."""
+
+PI: Final = math.pi
+"""The famous mathematical constant :math:`\\pi = 3.141592\\ldots`."""
+
+MU_0: Final = Q_(1.25663706212e-6, "H/m")
+"""Magnetic permeability of the vacuum :math:`\\mu_0 = 4 \\pi \\times 10^{-7}` (H/m)."""
+
+EPSILON_0: Final = Q_(8.8541878128e-12, "F/m")
+"""Vacuum permittivity :math:`\\varepsilon_0 =  8.8541878128 \\times 10^{-12}` (F/m)."""
+
+F: Final = Q_(50.0, "Hz")
+"""Network frequency :math:`f = 50` (Hz)."""
+
+OMEGA: Final = Q_(2 * PI * F, "rad/s")
+"""Angular frequency :math:`\\omega = 2 \\pi f` (rad/s)."""
+
+RHO: Final[dict[Material, Q_[float]]] = {
+    Material.CU: Q_(1.7241e-8, "ohm*m"),  # IEC 60287-1-1 Table 1
+    Material.AL: Q_(2.8264e-8, "ohm*m"),  # IEC 60287-1-1 Table 1
+    Material.AM: Q_(3.26e-8, "ohm*m"),  # verified
+    Material.AA: Q_(4.0587e-8, "ohm*m"),  # verified (approx. AS 3607 ACSR/GZ)
+    Material.LA: Q_(3.26e-8, "ohm*m"),
+}
+"""Resistivity of common conductor materials (Ohm.m)."""
+
+MU_R: Final[dict[Material, Q_[float]]] = {
+    Material.CU: Q_(0.9999935849131266),
+    Material.AL: Q_(1.0000222328028834),
+    Material.AM: Q_(0.9999705074463784),
+    Material.AA: Q_(1.0000222328028834),  # ==AL
+    Material.LA: Q_(0.9999705074463784),  # ==AM
+}
+"""Relative magnetic permeability of common conductor materials."""
+
+DELTA_P: Final[dict[Material, Q_[float]]] = {
+    Material.CU: Q_(9.33, "mm"),
+    Material.AL: Q_(11.95, "mm"),
+    Material.AM: Q_(12.85, "mm"),
+    Material.AA: Q_(14.34, "mm"),
+    Material.LA: Q_(12.85, "mm"),
+}
+"""Skin depth of common conductor materials :math:`\\sqrt{\\dfrac{\\rho}{\\pi f \\mu_r \\mu_0}}` (mm)."""
+# Skin depth is the depth at which the current density is reduced to 1/e (~37%) of the surface value.
+# Generated with:
+# ---------------
+# def delta_p(rho, mu_r):
+#     return np.sqrt(rho / (PI * F * mu_r * MU_0))
+# for material in Material:
+#     print(material, delta_p(RHO[material], MU_R[material]).m_as("mm"))
+
+TAN_D: Final[dict[Insulator, Q_[float]]] = {
+    Insulator.PVC: Q_(1000e-4),
+    Insulator.HDPE: Q_(10e-4),
+    Insulator.MDPE: Q_(10e-4),
+    Insulator.LDPE: Q_(10e-4),
+    Insulator.XLPE: Q_(40e-4),
+    Insulator.EPR: Q_(200e-4),
+    Insulator.IP: Q_(100e-4),
+    Insulator.NONE: Q_(0),
+}
+"""Loss angles of common insulator materials according to the IEC 60287 standard."""
+# IEC 60287-1-1 Table 3. We only include the MV values.
+
+EPSILON_R: Final[dict[Insulator, Q_[float]]] = {
+    Insulator.PVC: Q_(8.0),
+    Insulator.HDPE: Q_(2.3),
+    Insulator.MDPE: Q_(2.3),
+    Insulator.LDPE: Q_(2.3),
+    Insulator.XLPE: Q_(2.5),
+    Insulator.EPR: Q_(3.0),
+    Insulator.IP: Q_(4.0),
+    Insulator.NONE: Q_(1.0),
+}
+"""Relative permittivity of common insulator materials according to the IEC 60287 standard."""
+# IEC 60287-1-1 Table 3. We only include the MV values.
+
+
+def __getattr__(name: str):
+    import warnings
+
+    from roseau.load_flow.utils import find_stack_level
+
+    if name in ("PositiveSequence", "NegativeSequence", "ZeroSequence"):
+        warnings.warn(
+            f"'rlf.constants.{name}' is deprecated. Use 'rlf.sym.{name}' instead.",
+            category=FutureWarning,
+            stacklevel=find_stack_level(),
+        )
+        from roseau.load_flow import sym
+
+        return getattr(sym, name)
+
+    raise AttributeError(f"module '{__name__}' has no attribute '{name}'")
diff --git a/roseau/load_flow/converters.py b/roseau/load_flow/converters.py
index 54f2fec2..3c39f060 100644
--- a/roseau/load_flow/converters.py
+++ b/roseau/load_flow/converters.py
@@ -3,106 +3,19 @@
 
 Available functions:
 
-* convert between phasor and symmetrical components
 * convert potentials to voltages
 """
 
 import logging
 
 import numpy as np
-import pandas as pd
 
 from roseau.load_flow.exceptions import RoseauLoadFlowException, RoseauLoadFlowExceptionCode
 from roseau.load_flow.typing import ComplexArray, ComplexArrayLike1D
 from roseau.load_flow.units import Q_, ureg_wraps
-from roseau.load_flow.utils.constants import NegativeSequence, PositiveSequence, ZeroSequence
-from roseau.load_flow.utils.types import SequenceDtype
 
 logger = logging.getLogger(__name__)
 
-A = np.array([ZeroSequence, PositiveSequence, NegativeSequence], dtype=np.complex128)
-"""numpy.ndarray[complex]: "A" matrix: transformation matrix from phasor to symmetrical components."""
-
-_A_INV = np.linalg.inv(A)
-_SEQ_INDEX = pd.CategoricalIndex(["zero", "pos", "neg"], name="sequence", dtype=SequenceDtype)
-
-
-def phasor_to_sym(v_abc: ComplexArrayLike1D) -> ComplexArray:
-    """Compute the symmetrical components `(0, +, -)` from the phasor components `(a, b, c)`."""
-    v_abc_array = np.asarray(v_abc)
-    orig_shape = v_abc_array.shape
-    v_012 = _A_INV @ v_abc_array.reshape((3, 1))
-    return v_012.reshape(orig_shape)
-
-
-def sym_to_phasor(v_012: ComplexArrayLike1D) -> ComplexArray:
-    """Compute the phasor components `(a, b, c)` from the symmetrical components `(0, +, -)`."""
-    v_012_array = np.asarray(v_012)
-    orig_shape = v_012_array.shape
-    v_abc = A @ v_012_array.reshape((3, 1))
-    return v_abc.reshape(orig_shape)
-
-
-def series_phasor_to_sym(s_abc: pd.Series) -> pd.Series:
-    """Compute the symmetrical components `(0, +, -)` from the phasor components `(a, b, c)` of a series.
-
-    Args:
-        s_abc:
-            Series of phasor components (voltage, current, ...). The series must have a
-            multi-index with a `'phase'` level containing the phases in order (a -> b -> c).
-
-    Returns:
-        Series of the symmetrical components representing the input phasor series. The series has
-        a multi-index with the phase level replaced by a `'sequence'` level of values
-        `('zero', 'pos', 'neg')`.
-
-    Example:
-        Say we have a pandas series of three-phase voltages of every bus in the network:
-
-        >>> voltage
-        bus_id  phase
-        vs      an       200000000000.0+0.00000000j
-                bn      -10000.000000-17320.508076j
-                cn      -10000.000000+17320.508076j
-        bus     an       19999.00000095+0.00000000j
-                bn       -9999.975000-17320.464775j
-                cn       -9999.975000+17320.464775j
-        Name: voltage, dtype: complex128
-
-        We can get the `zero`, `positive`, and `negative` sequences of the voltage using:
-
-        >>> voltage_sym_components = series_phasor_to_sym(voltage)
-        >>> voltage_sym_components
-        bus_id  sequence
-        bus     zero        3.183231e-12-9.094947e-13j
-                pos         1.999995e+04+3.283594e-12j
-                neg        -1.796870e-07-2.728484e-12j
-        vs      zero        5.002221e-12-9.094947e-13j
-                pos         2.000000e+04+3.283596e-12j
-                neg        -1.796880e-07-1.818989e-12j
-        Name: voltage, dtype: complex128
-
-        We can now access each sequence of the symmetrical components individually:
-
-        >>> voltage_sym_components.loc[:, "zero"]  # get zero sequence values
-        bus_id
-        bus     3.183231e-12-9.094947e-13j
-        vs      5.002221e-12-9.094947e-13j
-        Name: voltage, dtype: complex128
-
-    """
-    if not isinstance(s_abc, pd.Series):
-        raise TypeError("Input must be a pandas Series.")
-    if not isinstance(s_abc.index, pd.MultiIndex):
-        raise ValueError("Input series must have a MultiIndex.")
-    if "phase" not in s_abc.index.names:
-        raise ValueError("Input series must have a 'phase' level in the MultiIndex.")
-    level_names = [name for name in s_abc.index.names if name != "phase"]
-    s_012 = s_abc.groupby(level=level_names, sort=False).apply(
-        lambda s: pd.Series(_A_INV @ s, index=_SEQ_INDEX, dtype=np.complex128)
-    )
-    return s_012
-
 
 def _calculate_voltages(potentials: ComplexArray, phases: str) -> ComplexArray:
     if len(potentials) != len(phases):
@@ -199,3 +112,21 @@ def calculate_voltage_phases(phases: str) -> list[str]:
         msg = f"Invalid phases '{phases}'. Must be one of {', '.join(_VOLTAGE_PHASES_CACHE)}."
         logger.error(msg)
         raise RoseauLoadFlowException(msg, code=RoseauLoadFlowExceptionCode.BAD_PHASE) from None
+
+
+def __getattr__(name: str):
+    import warnings
+
+    from roseau.load_flow.utils import find_stack_level
+
+    if name in ("phasor_to_sym", "sym_to_phasor", "series_phasor_to_sym"):
+        warnings.warn(
+            f"'rlf.converters.{name}' is deprecated. Use 'rlf.sym.{name}' instead.",
+            category=FutureWarning,
+            stacklevel=find_stack_level(),
+        )
+        from roseau.load_flow import sym
+
+        return getattr(sym, name)
+
+    raise AttributeError(f"module '{__name__}' has no attribute '{name}'")
diff --git a/roseau/load_flow/io/dgs/__init__.py b/roseau/load_flow/io/dgs/__init__.py
index 2a789a10..8d83e38f 100644
--- a/roseau/load_flow/io/dgs/__init__.py
+++ b/roseau/load_flow/io/dgs/__init__.py
@@ -33,7 +33,7 @@
     VoltageSource,
 )
 from roseau.load_flow.typing import Id, StrPath
-from roseau.load_flow.utils._exceptions import find_stack_level
+from roseau.load_flow.utils import find_stack_level
 
 logger = logging.getLogger(__name__)
 
diff --git a/roseau/load_flow/io/dgs/constants.py b/roseau/load_flow/io/dgs/constants.py
index 1193a33d..fe6a33c6 100644
--- a/roseau/load_flow/io/dgs/constants.py
+++ b/roseau/load_flow/io/dgs/constants.py
@@ -1,6 +1,6 @@
 from typing import Final, Literal, TypeAlias
 
-from roseau.load_flow.utils import Insulator, LineType, Material
+from roseau.load_flow.types import Insulator, LineType, Material
 
 # Lines
 LINE_TYPES: Final[dict[int, LineType]] = {
diff --git a/roseau/load_flow/io/dgs/lines.py b/roseau/load_flow/io/dgs/lines.py
index fb6b0d46..8305e30c 100644
--- a/roseau/load_flow/io/dgs/lines.py
+++ b/roseau/load_flow/io/dgs/lines.py
@@ -8,7 +8,7 @@
 from roseau.load_flow.io.dgs.constants import INSULATORS, LINE_TYPES, MATERIALS
 from roseau.load_flow.models import Bus, Ground, Line, LineParameters
 from roseau.load_flow.typing import Id
-from roseau.load_flow.utils._exceptions import find_stack_level
+from roseau.load_flow.utils import find_stack_level
 
 logger = logging.getLogger(__name__)
 
diff --git a/roseau/load_flow/io/dgs/switches.py b/roseau/load_flow/io/dgs/switches.py
index 0ae9ffeb..fa9f289f 100644
--- a/roseau/load_flow/io/dgs/switches.py
+++ b/roseau/load_flow/io/dgs/switches.py
@@ -7,7 +7,7 @@
 from roseau.load_flow.exceptions import RoseauLoadFlowException, RoseauLoadFlowExceptionCode
 from roseau.load_flow.models import Bus, Switch
 from roseau.load_flow.typing import Id
-from roseau.load_flow.utils._exceptions import find_stack_level
+from roseau.load_flow.utils import find_stack_level
 
 logger = logging.getLogger(__name__)
 
diff --git a/roseau/load_flow/io/dict.py b/roseau/load_flow/io/dict.py
index ff221dcb..e43ba51f 100644
--- a/roseau/load_flow/io/dict.py
+++ b/roseau/load_flow/io/dict.py
@@ -28,7 +28,7 @@
     VoltageSource,
 )
 from roseau.load_flow.typing import Id, JsonDict
-from roseau.load_flow.utils._exceptions import find_stack_level
+from roseau.load_flow.utils import find_stack_level
 
 if TYPE_CHECKING:
     from roseau.load_flow.network import ElectricalNetwork
diff --git a/roseau/load_flow/io/tests/test_dict.py b/roseau/load_flow/io/tests/test_dict.py
index 8ccf8125..e5f9534d 100644
--- a/roseau/load_flow/io/tests/test_dict.py
+++ b/roseau/load_flow/io/tests/test_dict.py
@@ -22,7 +22,7 @@
 )
 from roseau.load_flow.network import ElectricalNetwork
 from roseau.load_flow.testing import assert_json_close
-from roseau.load_flow.utils import Insulator, LineType, Material
+from roseau.load_flow.types import Insulator, LineType, Material
 
 # Store the expected hashes of the files that should not be modified
 EXPECTED_HASHES = {
diff --git a/roseau/load_flow/models/buses.py b/roseau/load_flow/models/buses.py
index e9c79763..bf596936 100644
--- a/roseau/load_flow/models/buses.py
+++ b/roseau/load_flow/models/buses.py
@@ -9,12 +9,13 @@
 from shapely.geometry.base import BaseGeometry
 from typing_extensions import Self
 
-from roseau.load_flow.converters import _calculate_voltages, calculate_voltage_phases, phasor_to_sym
+from roseau.load_flow.converters import _calculate_voltages, calculate_voltage_phases
 from roseau.load_flow.exceptions import RoseauLoadFlowException, RoseauLoadFlowExceptionCode
 from roseau.load_flow.models.core import Element
+from roseau.load_flow.sym import phasor_to_sym
 from roseau.load_flow.typing import ComplexArray, ComplexArrayLike1D, FloatArray, Id, JsonDict
 from roseau.load_flow.units import Q_, ureg_wraps
-from roseau.load_flow.utils._exceptions import find_stack_level
+from roseau.load_flow.utils import find_stack_level
 from roseau.load_flow_engine.cy_engine import CyBus
 
 logger = logging.getLogger(__name__)
diff --git a/roseau/load_flow/models/core.py b/roseau/load_flow/models/core.py
index 0679b5b2..c75afeee 100644
--- a/roseau/load_flow/models/core.py
+++ b/roseau/load_flow/models/core.py
@@ -9,8 +9,7 @@
 
 from roseau.load_flow.exceptions import RoseauLoadFlowException, RoseauLoadFlowExceptionCode
 from roseau.load_flow.typing import Id
-from roseau.load_flow.utils import Identifiable, JsonMixin
-from roseau.load_flow.utils._exceptions import find_stack_level
+from roseau.load_flow.utils import Identifiable, JsonMixin, find_stack_level
 from roseau.load_flow_engine.cy_engine import CyElement
 
 if TYPE_CHECKING:
diff --git a/roseau/load_flow/models/lines/parameters.py b/roseau/load_flow/models/lines/parameters.py
index 10ec501b..42066637 100644
--- a/roseau/load_flow/models/lines/parameters.py
+++ b/roseau/load_flow/models/lines/parameters.py
@@ -13,7 +13,9 @@
 from typing_extensions import Self
 
 from roseau.load_flow._compat import StrEnum
+from roseau.load_flow.constants import EPSILON_0, EPSILON_R, MU_0, OMEGA, PI, RHO, TAN_D, F
 from roseau.load_flow.exceptions import RoseauLoadFlowException, RoseauLoadFlowExceptionCode
+from roseau.load_flow.types import Insulator, LineType, Material
 from roseau.load_flow.typing import (
     ComplexArray,
     ComplexArrayLike2D,
@@ -24,23 +26,7 @@
     JsonDict,
 )
 from roseau.load_flow.units import Q_, ureg_wraps
-from roseau.load_flow.utils import (
-    EPSILON_0,
-    EPSILON_R,
-    MU_0,
-    OMEGA,
-    PI,
-    RHO,
-    TAN_D,
-    CatalogueMixin,
-    F,
-    Identifiable,
-    Insulator,
-    JsonMixin,
-    LineType,
-    Material,
-)
-from roseau.load_flow.utils._exceptions import find_stack_level
+from roseau.load_flow.utils import CatalogueMixin, Identifiable, JsonMixin, find_stack_level
 
 logger = logging.getLogger(__name__)
 
@@ -56,9 +42,9 @@
     LineType.UNDERGROUND: Insulator.PVC,
 }
 
-MaterialArray: TypeAlias = NDArray[Material]
-InsulatorArray: TypeAlias = NDArray[Insulator]
-_StrEnumType: TypeAlias = TypeVar("_StrEnumType", bound=StrEnum)
+MaterialArray: TypeAlias = NDArray[np.object_]
+InsulatorArray: TypeAlias = NDArray[np.object_]
+_StrEnumType = TypeVar("_StrEnumType", bound=StrEnum)
 
 
 class LineParameters(Identifiable, JsonMixin, CatalogueMixin[pd.DataFrame]):
diff --git a/roseau/load_flow/models/loads/flexible_parameters.py b/roseau/load_flow/models/loads/flexible_parameters.py
index bf94195b..4fd87322 100644
--- a/roseau/load_flow/models/loads/flexible_parameters.py
+++ b/roseau/load_flow/models/loads/flexible_parameters.py
@@ -17,8 +17,7 @@
     ProjectionType,
 )
 from roseau.load_flow.units import Q_, ureg_wraps
-from roseau.load_flow.utils import JsonMixin, _optional_deps
-from roseau.load_flow.utils._exceptions import find_stack_level
+from roseau.load_flow.utils import JsonMixin, find_stack_level, optional_deps
 from roseau.load_flow_engine.cy_engine import CyControl, CyFlexibleParameter, CyProjection
 
 logger = logging.getLogger(__name__)
@@ -1168,7 +1167,7 @@ def plot_pq(
             The axis on which the plot has been drawn and the resulting flexible powers (the input if not `None` else
             the computed values).
         """
-        plt = _optional_deps.pyplot  # this line first for better error handling
+        plt = optional_deps.pyplot  # this line first for better error handling
         from matplotlib import colormaps, patheffects
 
         # Get the axes
@@ -1262,7 +1261,7 @@ def plot_control_p(
             The axis on which the plot has been drawn and the resulting flexible powers (the input if not `None` else
             the computed values).
         """
-        plt = _optional_deps.pyplot
+        plt = optional_deps.pyplot
 
         # Get the axes
         if ax is None:
@@ -1312,7 +1311,7 @@ def plot_control_q(
             The axis on which the plot has been drawn and the resulting flexible powers (the input if not `None` else
             the computed values).
         """
-        plt = _optional_deps.pyplot
+        plt = optional_deps.pyplot
 
         # Get the axes
         if ax is None:
diff --git a/roseau/load_flow/models/loads/loads.py b/roseau/load_flow/models/loads/loads.py
index 435f098c..67b1547e 100644
--- a/roseau/load_flow/models/loads/loads.py
+++ b/roseau/load_flow/models/loads/loads.py
@@ -13,7 +13,7 @@
 from roseau.load_flow.models.loads.flexible_parameters import FlexibleParameter
 from roseau.load_flow.typing import ComplexArray, ComplexScalarOrArrayLike1D, Id, JsonDict
 from roseau.load_flow.units import Q_, ureg_wraps
-from roseau.load_flow.utils._exceptions import find_stack_level
+from roseau.load_flow.utils import find_stack_level
 from roseau.load_flow_engine.cy_engine import (
     CyAdmittanceLoad,
     CyCurrentLoad,
diff --git a/roseau/load_flow/models/potential_refs.py b/roseau/load_flow/models/potential_refs.py
index 4e9d3a23..5e9e4f02 100644
--- a/roseau/load_flow/models/potential_refs.py
+++ b/roseau/load_flow/models/potential_refs.py
@@ -10,7 +10,7 @@
 from roseau.load_flow.models.grounds import Ground
 from roseau.load_flow.typing import Id, JsonDict
 from roseau.load_flow.units import Q_, ureg_wraps
-from roseau.load_flow.utils._exceptions import find_stack_level
+from roseau.load_flow.utils import find_stack_level
 from roseau.load_flow_engine.cy_engine import CyDeltaPotentialRef, CyPotentialRef
 
 logger = logging.getLogger(__name__)
diff --git a/roseau/load_flow/models/sources.py b/roseau/load_flow/models/sources.py
index f706819d..55665060 100644
--- a/roseau/load_flow/models/sources.py
+++ b/roseau/load_flow/models/sources.py
@@ -10,10 +10,10 @@
 from roseau.load_flow.exceptions import RoseauLoadFlowException, RoseauLoadFlowExceptionCode
 from roseau.load_flow.models.buses import Bus
 from roseau.load_flow.models.core import Element
+from roseau.load_flow.sym import PositiveSequence
 from roseau.load_flow.typing import ComplexArray, ComplexScalarOrArrayLike1D, Id, JsonDict
 from roseau.load_flow.units import Q_, ureg_wraps
-from roseau.load_flow.utils._exceptions import find_stack_level
-from roseau.load_flow.utils.constants import PositiveSequence
+from roseau.load_flow.utils import find_stack_level
 from roseau.load_flow_engine.cy_engine import CyDeltaVoltageSource, CyVoltageSource
 
 logger = logging.getLogger(__name__)
diff --git a/roseau/load_flow/models/tests/test_line_parameters.py b/roseau/load_flow/models/tests/test_line_parameters.py
index 5791109b..a5f9a5b4 100644
--- a/roseau/load_flow/models/tests/test_line_parameters.py
+++ b/roseau/load_flow/models/tests/test_line_parameters.py
@@ -8,8 +8,8 @@
 
 from roseau.load_flow.exceptions import RoseauLoadFlowException, RoseauLoadFlowExceptionCode
 from roseau.load_flow.models import Bus, Ground, Line, LineParameters
+from roseau.load_flow.types import Insulator, LineType, Material
 from roseau.load_flow.units import Q_
-from roseau.load_flow.utils import Insulator, LineType, Material
 
 
 def test_line_parameters():
@@ -214,12 +214,12 @@ def test_from_geometry():
         external_diameter=0.04,
     )
     np.testing.assert_allclose(
-        lp.y_shunt.m.imag * 4,  # because InsulatorType.IP has 4x epsilon_r
+        lp.y_shunt.m.imag * 4,  # because Insulator.IP has 4x epsilon_r
         LineParameters.from_geometry(
             id="test",
             line_type=lp.line_type,
             material=lp.materials[0],
-            insulator=Insulator.IP,  # 4x epsilon_r of InsulatorType.NONE
+            insulator=Insulator.IP,  # 4x epsilon_r of Insulator.NONE
             section=lp.sections[0],
             height=-0.5,
             external_diameter=0.04,
diff --git a/roseau/load_flow/models/tests/test_lines.py b/roseau/load_flow/models/tests/test_lines.py
index 4dfd4a6d..cfb1775d 100644
--- a/roseau/load_flow/models/tests/test_lines.py
+++ b/roseau/load_flow/models/tests/test_lines.py
@@ -4,8 +4,8 @@
 
 from roseau.load_flow.exceptions import RoseauLoadFlowException, RoseauLoadFlowExceptionCode
 from roseau.load_flow.models import Bus, Ground, Line, LineParameters
+from roseau.load_flow.sym import PositiveSequence as PosSeq
 from roseau.load_flow.units import Q_
-from roseau.load_flow.utils import PositiveSequence as PosSeq
 
 
 def test_lines_length():
diff --git a/roseau/load_flow/network.py b/roseau/load_flow/network.py
index 17294f98..eef94584 100644
--- a/roseau/load_flow/network.py
+++ b/roseau/load_flow/network.py
@@ -22,6 +22,7 @@
 from typing_extensions import Self
 
 from roseau.load_flow._solvers import AbstractSolver
+from roseau.load_flow.constants import SQRT3
 from roseau.load_flow.exceptions import RoseauLoadFlowException, RoseauLoadFlowExceptionCode
 from roseau.load_flow.io import network_from_dgs, network_from_dict, network_to_dict
 from roseau.load_flow.models import (
@@ -40,9 +41,15 @@
     VoltageSource,
 )
 from roseau.load_flow.typing import Id, JsonDict, MapOrSeq, Solver, StrPath
-from roseau.load_flow.utils import SQRT3, CatalogueMixin, JsonMixin, _optional_deps
-from roseau.load_flow.utils._exceptions import find_stack_level
-from roseau.load_flow.utils.types import _DTYPES, LoadTypeDtype, VoltagePhaseDtype
+from roseau.load_flow.utils import (
+    DTYPES,
+    CatalogueMixin,
+    JsonMixin,
+    LoadTypeDtype,
+    VoltagePhaseDtype,
+    find_stack_level,
+    optional_deps,
+)
 from roseau.load_flow_engine.cy_engine import CyElectricalNetwork
 
 if TYPE_CHECKING:
@@ -462,7 +469,7 @@ def to_graph(self) -> "Graph":
             This method requires *networkx* to be installed. You can install it with the ``"graph"``
             extra if you are using pip: ``pip install "roseau-load-flow[graph]"``.
         """
-        nx = _optional_deps.networkx
+        nx = optional_deps.networkx
         graph = nx.Graph()
         for bus in self.buses.values():
             graph.add_node(bus.id, geom=bus.geometry)
@@ -667,7 +674,7 @@ def res_buses(self) -> pd.DataFrame:
         """
         self._check_valid_results()
         res_dict = {"bus_id": [], "phase": [], "potential": []}
-        dtypes = {c: _DTYPES[c] for c in res_dict}
+        dtypes = {c: DTYPES[c] for c in res_dict}
         for bus_id, bus in self.buses.items():
             for potential, phase in zip(bus._res_potentials_getter(warning=False), bus.phases, strict=True):
                 res_dict["bus_id"].append(bus_id)
@@ -707,7 +714,7 @@ def res_buses_voltages(self) -> pd.DataFrame:
             "max_voltage_level": [],
             "nominal_voltage": [],
         }
-        dtypes = {c: _DTYPES[c] for c in voltages_dict} | {"phase": VoltagePhaseDtype}
+        dtypes = {c: DTYPES[c] for c in voltages_dict} | {"phase": VoltagePhaseDtype}
         for bus_id, bus in self.buses.items():
             nominal_voltage = bus._nominal_voltage
             min_voltage_level = bus._min_voltage_level
@@ -808,7 +815,7 @@ def res_lines(self) -> pd.DataFrame:
             "max_loading": [],
             "ampacity": [],
         }
-        dtypes = {c: _DTYPES[c] for c in res_dict}
+        dtypes = {c: DTYPES[c] for c in res_dict}
         for line in self.lines.values():
             currents1, currents2 = line._res_currents_getter(warning=False)
             potentials1, potentials2 = line._res_potentials_getter(warning=False)
@@ -905,7 +912,7 @@ def res_transformers(self) -> pd.DataFrame:
             "max_loading": [],
             "sn": [],
         }
-        dtypes = {c: _DTYPES[c] for c in res_dict}
+        dtypes = {c: DTYPES[c] for c in res_dict}
         for transformer in self.transformers.values():
             currents1, currents2 = transformer._res_currents_getter(warning=False)
             potentials1, potentials2 = transformer._res_potentials_getter(warning=False)
@@ -972,7 +979,7 @@ def res_switches(self) -> pd.DataFrame:
             "potential1": [],
             "potential2": [],
         }
-        dtypes = {c: _DTYPES[c] for c in res_dict}
+        dtypes = {c: DTYPES[c] for c in res_dict}
         for switch in self.switches.values():
             if not isinstance(switch, Switch):
                 continue
@@ -1009,7 +1016,7 @@ def res_loads(self) -> pd.DataFrame:
         """
         self._check_valid_results()
         res_dict = {"load_id": [], "phase": [], "type": [], "current": [], "power": [], "potential": []}
-        dtypes = {c: _DTYPES[c] for c in res_dict} | {"type": LoadTypeDtype}
+        dtypes = {c: DTYPES[c] for c in res_dict} | {"type": LoadTypeDtype}
         for load_id, load in self.loads.items():
             currents = load._res_currents_getter(warning=False)
             potentials = load._res_potentials_getter(warning=False)
@@ -1038,7 +1045,7 @@ def res_loads_voltages(self) -> pd.DataFrame:
         """
         self._check_valid_results()
         voltages_dict = {"load_id": [], "phase": [], "type": [], "voltage": []}
-        dtypes = {c: _DTYPES[c] for c in voltages_dict} | {"phase": VoltagePhaseDtype, "type": LoadTypeDtype}
+        dtypes = {c: DTYPES[c] for c in voltages_dict} | {"phase": VoltagePhaseDtype, "type": LoadTypeDtype}
         for load_id, load in self.loads.items():
             for voltage, phase in zip(load._res_voltages_getter(warning=False), load.voltage_phases, strict=True):
                 voltages_dict["load_id"].append(load_id)
@@ -1064,7 +1071,7 @@ def res_loads_flexible_powers(self) -> pd.DataFrame:
         """
         self._check_valid_results()
         loads_dict = {"load_id": [], "phase": [], "flexible_power": []}
-        dtypes = {c: _DTYPES[c] for c in loads_dict} | {"phase": VoltagePhaseDtype}
+        dtypes = {c: DTYPES[c] for c in loads_dict} | {"phase": VoltagePhaseDtype}
         for load_id, load in self.loads.items():
             if not (isinstance(load, PowerLoad) and load.is_flexible):
                 continue
@@ -1091,7 +1098,7 @@ def res_sources(self) -> pd.DataFrame:
         """
         self._check_valid_results()
         res_dict = {"source_id": [], "phase": [], "current": [], "power": [], "potential": []}
-        dtypes = {c: _DTYPES[c] for c in res_dict}
+        dtypes = {c: DTYPES[c] for c in res_dict}
         for source_id, source in self.sources.items():
             currents = source._res_currents_getter(warning=False)
             potentials = source._res_potentials_getter(warning=False)
@@ -1115,7 +1122,7 @@ def res_grounds(self) -> pd.DataFrame:
         """
         self._check_valid_results()
         res_dict = {"ground_id": [], "potential": []}
-        dtypes = {c: _DTYPES[c] for c in res_dict}
+        dtypes = {c: DTYPES[c] for c in res_dict}
         for ground in self.grounds.values():
             potential = ground._res_potential_getter(warning=False)
             res_dict["ground_id"].append(ground.id)
@@ -1134,7 +1141,7 @@ def res_potential_refs(self) -> pd.DataFrame:
         """
         self._check_valid_results()
         res_dict = {"potential_ref_id": [], "current": []}
-        dtypes = {c: _DTYPES[c] for c in res_dict}
+        dtypes = {c: DTYPES[c] for c in res_dict}
         for p_ref in self.potential_refs.values():
             current = p_ref._res_current_getter(warning=False)
             res_dict["potential_ref_id"].append(p_ref.id)
diff --git a/roseau/load_flow/plotting.py b/roseau/load_flow/plotting.py
index 8df052d2..ddc8919f 100644
--- a/roseau/load_flow/plotting.py
+++ b/roseau/load_flow/plotting.py
@@ -5,7 +5,9 @@
 
 import numpy as np
 
-import roseau.load_flow as rlf
+from roseau.load_flow.models import AbstractLoad, Bus, VoltageSource
+from roseau.load_flow.network import ElectricalNetwork
+from roseau.load_flow.sym import phasor_to_sym
 from roseau.load_flow.typing import Complex, ComplexArray, Float
 
 if TYPE_CHECKING:
@@ -76,9 +78,7 @@ def _draw_voltage_phasor(ax: "Axes", potential1: Complex, potential2: Complex, c
 #
 # Phasor plotting functions
 #
-def plot_voltage_phasors(
-    element: rlf.Bus | rlf.AbstractLoad | rlf.VoltageSource, *, ax: "Axes | None" = None
-) -> "Axes":
+def plot_voltage_phasors(element: Bus | AbstractLoad | VoltageSource, *, ax: "Axes | None" = None) -> "Axes":
     """Plot the voltage phasors of a bus, load, or voltage source.
 
     Args:
@@ -91,7 +91,7 @@ def plot_voltage_phasors(
     Returns:
         The axes with the plot.
     """
-    from roseau.load_flow.utils._optional_deps import pyplot as plt
+    from roseau.load_flow.utils.optional_deps import pyplot as plt
 
     if ax is None:
         ax = plt.gca()
@@ -121,9 +121,7 @@ def plot_voltage_phasors(
     return ax
 
 
-def plot_symmetrical_voltages(
-    element: rlf.Bus | rlf.AbstractLoad | rlf.VoltageSource, *, ax: "Axes | None" = None
-) -> "Axes":
+def plot_symmetrical_voltages(element: Bus | AbstractLoad | VoltageSource, *, ax: "Axes | None" = None) -> "Axes":
     """Plot the symmetrical voltages of a bus, load, or voltage source.
 
     Args:
@@ -137,13 +135,13 @@ def plot_symmetrical_voltages(
     Returns:
         The axes with the plot.
     """
-    from roseau.load_flow.utils._optional_deps import pyplot as plt
+    from roseau.load_flow.utils.optional_deps import pyplot as plt
 
     if element.phases not in {"abc", "abcn"}:
         raise ValueError("The element must have 'abc' or 'abcn' phases.")
     if ax is None:
         ax = plt.gca()
-    voltages_sym = rlf.converters.phasor_to_sym(element.res_voltages.m)
+    voltages_sym = phasor_to_sym(element.res_voltages.m)
     _configure_axes(ax, voltages_sym)
     ax.set_title(f"{element.id} (symmetrical)")
     for sequence, voltage in zip(("zero", "pos", "neg"), voltages_sym, strict=True):
@@ -158,7 +156,7 @@ def plot_symmetrical_voltages(
 # Map plotting functions
 #
 def plot_interactive_map(
-    network: rlf.ElectricalNetwork,
+    network: ElectricalNetwork,
     *,
     style_color: str = "#234e83",
     highlight_color: str = "#cad40e",
diff --git a/roseau/load_flow/sym.py b/roseau/load_flow/sym.py
new file mode 100644
index 00000000..26eeb0e8
--- /dev/null
+++ b/roseau/load_flow/sym.py
@@ -0,0 +1,110 @@
+"""Symmetrical components utilities."""
+
+from typing import Final
+
+import numpy as np
+import pandas as pd
+
+from roseau.load_flow.constants import ALPHA, ALPHA2
+from roseau.load_flow.typing import ComplexArray, ComplexArrayLike1D
+from roseau.load_flow.utils.dtypes import SequenceDtype
+
+__all__ = [
+    "A",
+    "PositiveSequence",
+    "NegativeSequence",
+    "ZeroSequence",
+    "phasor_to_sym",
+    "sym_to_phasor",
+    "series_phasor_to_sym",
+]
+
+PositiveSequence: Final = np.array([1, ALPHA2, ALPHA], dtype=np.complex128)
+"""numpy.ndarray[complex]: Unit vector of positive-sequence components of a three-phase system."""
+NegativeSequence: Final = np.array([1, ALPHA, ALPHA2], dtype=np.complex128)
+"""numpy.ndarray[complex]: Unit vector of negative-sequence components of a three-phase system."""
+ZeroSequence: Final = np.array([1, 1, 1], dtype=np.complex128)
+"""numpy.ndarray[complex]: Unit vector of zero-sequence components of a three-phase system."""
+
+A: Final = np.array([ZeroSequence, PositiveSequence, NegativeSequence], dtype=np.complex128)
+"""numpy.ndarray[complex]: `A` matrix: transformation matrix from phasor to symmetrical components."""
+
+_A_INV = np.linalg.inv(A)
+_SEQ_INDEX = pd.CategoricalIndex(["zero", "pos", "neg"], name="sequence", dtype=SequenceDtype)
+
+
+def phasor_to_sym(v_abc: ComplexArrayLike1D) -> ComplexArray:
+    """Compute the symmetrical components `(0, +, -)` from the phasor components `(a, b, c)`."""
+    v_abc_array = np.asarray(v_abc)
+    orig_shape = v_abc_array.shape
+    v_012 = _A_INV @ v_abc_array.reshape((3, 1))
+    return v_012.reshape(orig_shape)
+
+
+def sym_to_phasor(v_012: ComplexArrayLike1D) -> ComplexArray:
+    """Compute the phasor components `(a, b, c)` from the symmetrical components `(0, +, -)`."""
+    v_012_array = np.asarray(v_012)
+    orig_shape = v_012_array.shape
+    v_abc = A @ v_012_array.reshape((3, 1))
+    return v_abc.reshape(orig_shape)
+
+
+def series_phasor_to_sym(s_abc: pd.Series) -> pd.Series:
+    """Compute the symmetrical components `(0, +, -)` from the phasor components `(a, b, c)` of a series.
+
+    Args:
+        s_abc:
+            Series of phasor components (voltage, current, ...). The series must have a
+            multi-index with a `'phase'` level containing the phases in order (a -> b -> c).
+
+    Returns:
+        Series of the symmetrical components representing the input phasor series. The series has
+        a multi-index with the phase level replaced by a `'sequence'` level of values
+        `('zero', 'pos', 'neg')`.
+
+    Example:
+        Say we have a pandas series of three-phase voltages of every bus in the network:
+
+        >>> voltage
+        bus_id  phase
+        vs      an       200000000000.0+0.00000000j
+                bn      -10000.000000-17320.508076j
+                cn      -10000.000000+17320.508076j
+        bus     an       19999.00000095+0.00000000j
+                bn       -9999.975000-17320.464775j
+                cn       -9999.975000+17320.464775j
+        Name: voltage, dtype: complex128
+
+        We can get the `zero`, `positive`, and `negative` sequences of the voltage using:
+
+        >>> voltage_sym_components = series_phasor_to_sym(voltage)
+        >>> voltage_sym_components
+        bus_id  sequence
+        bus     zero        3.183231e-12-9.094947e-13j
+                pos         1.999995e+04+3.283594e-12j
+                neg        -1.796870e-07-2.728484e-12j
+        vs      zero        5.002221e-12-9.094947e-13j
+                pos         2.000000e+04+3.283596e-12j
+                neg        -1.796880e-07-1.818989e-12j
+        Name: voltage, dtype: complex128
+
+        We can now access each sequence of the symmetrical components individually:
+
+        >>> voltage_sym_components.loc[:, "zero"]  # get zero sequence values
+        bus_id
+        bus     3.183231e-12-9.094947e-13j
+        vs      5.002221e-12-9.094947e-13j
+        Name: voltage, dtype: complex128
+
+    """
+    if not isinstance(s_abc, pd.Series):
+        raise TypeError("Input must be a pandas Series.")
+    if not isinstance(s_abc.index, pd.MultiIndex):
+        raise ValueError("Input series must have a MultiIndex.")
+    if "phase" not in s_abc.index.names:
+        raise ValueError("Input series must have a 'phase' level in the MultiIndex.")
+    level_names = [name for name in s_abc.index.names if name != "phase"]
+    s_012 = s_abc.groupby(level=level_names, sort=False).apply(
+        lambda s: pd.Series(_A_INV @ s, index=_SEQ_INDEX, dtype=np.complex128)
+    )
+    return s_012
diff --git a/roseau/load_flow/tests/test_converters.py b/roseau/load_flow/tests/test_converters.py
index 80e42ba2..aee34610 100644
--- a/roseau/load_flow/tests/test_converters.py
+++ b/roseau/load_flow/tests/test_converters.py
@@ -1,118 +1,7 @@
 import numpy as np
-import pandas as pd
-from pandas.testing import assert_series_equal
 
-from roseau.load_flow.converters import calculate_voltages, phasor_to_sym, series_phasor_to_sym, sym_to_phasor
+from roseau.load_flow.converters import calculate_voltages
 from roseau.load_flow.units import Q_, ureg
-from roseau.load_flow.utils import SequenceDtype
-
-
-def test_phasor_to_sym():
-    # Tests verified with https://phillipopambuh.info/portfolio/calculator--symmetrical_components.html
-    va = 230 + 0j
-    vb = 230 * np.e ** (1j * 4 * np.pi / 3)
-    vc = 230 * np.e ** (1j * 2 * np.pi / 3)
-
-    # Test balanced direct system: positive sequence
-    expected = np.array([0, 230, 0], dtype=complex)
-    assert np.allclose(phasor_to_sym([va, vb, vc]), expected)
-    # Also test numpy array input with different shapes
-    assert np.allclose(phasor_to_sym(np.array([va, vb, vc])), expected)
-    assert np.allclose(phasor_to_sym(np.array([[va], [vb], [vc]])), expected.reshape((3, 1)))
-
-    # Test balanced indirect system: negative sequence
-    expected = np.array([0, 0, 230], dtype=complex)
-    assert np.allclose(phasor_to_sym([va, vc, vb]), expected)
-
-    # Test unbalanced system: zero sequence
-    expected = np.array([230, 0, 0], dtype=complex)
-    assert np.allclose(phasor_to_sym([va, va, va]), expected)
-
-    # Test unbalanced system: general case
-    va = 200 + 0j
-    expected = np.array([10 * np.exp(1j * np.pi), 220, 10 * np.exp(1j * np.pi)], dtype=complex)
-    assert np.allclose(phasor_to_sym([va, vb, vc]), expected)
-
-
-def test_sym_to_phasor():
-    # Tests verified with https://phillipopambuh.info/portfolio/calculator--symmetrical_components.html
-    va = 230 + 0j
-    vb = 230 * np.e ** (1j * 4 * np.pi / 3)
-    vc = 230 * np.e ** (1j * 2 * np.pi / 3)
-
-    # Test balanced direct system: positive sequence
-    expected = np.array([va, vb, vc], dtype=complex)
-    assert np.allclose(sym_to_phasor([0, va, 0]), expected)
-    # Also test numpy array input with different shapes
-    assert np.allclose(sym_to_phasor(np.array([0, va, 0])), expected)
-    assert np.allclose(sym_to_phasor(np.array([[0], [va], [0]])), expected.reshape((3, 1)))
-
-    # Test balanced indirect system: negative sequence
-    expected = np.array([va, vc, vb], dtype=complex)
-    assert np.allclose(sym_to_phasor([0, 0, va]), expected)
-
-    # Test unbalanced system: zero sequence
-    expected = np.array([va, va, va], dtype=complex)
-    assert np.allclose(sym_to_phasor([va, 0, 0]), expected)
-
-    # Test unbalanced system: general case
-    va = 200 + 0j
-    expected = np.array([va, vb, vc], dtype=complex)
-    assert np.allclose(sym_to_phasor([10 * np.e ** (1j * np.pi), 220, 10 * np.e ** (1j * np.pi)]), expected)
-
-
-def test_phasor_sym_roundtrip():
-    va = 230 + 0j
-    vb = 230 * np.e ** (1j * 4 * np.pi / 3)
-    vc = 230 * np.e ** (1j * 2 * np.pi / 3)
-
-    # Test balanced direct system: positive sequence
-    assert np.allclose(sym_to_phasor(phasor_to_sym([va, vb, vc])), np.array([va, vb, vc]))
-
-    # Test balanced indirect system: negative sequence
-    assert np.allclose(sym_to_phasor(phasor_to_sym([va, vc, vb])), np.array([va, vc, vb]))
-
-    # Test unbalanced system: zero sequence
-    assert np.allclose(sym_to_phasor(phasor_to_sym([va, va, va])), np.array([va, va, va]))
-
-    # Test unbalanced system: general case
-    va = 200 + 0j
-    assert np.allclose(sym_to_phasor(phasor_to_sym([va, vb, vc])), np.array([va, vb, vc]))
-
-
-def test_series_phasor_to_sym():
-    va = 230 + 0j
-    vb = 230 * np.exp(1j * 4 * np.pi / 3)
-    vc = 230 * np.exp(1j * 2 * np.pi / 3)
-
-    # Test with different phases per bus, different systems, different magnitudes
-    # fmt: off
-    voltage_data = {
-        # Direct system (positive sequence)
-        ("bus1", "a"): va, ("bus1", "b"): vb, ("bus1", "c"): vc,
-        # Indirect system (negative sequence)
-        ("bus2", "an"): va / 2, ("bus2", "bn"): vc / 2, ("bus2", "cn"): vb / 2,
-        # Unbalanced system (zero sequence)
-        ("bus3", "ab"): va, ("bus3", "bc"): va, ("bus3", "ca"): va,
-    }
-    expected_sym_data = [
-        0, va, 0,  # Direct system (positive sequence)
-        0, 0, va / 2,  # Indirect system (negative sequence)
-        va, 0, 0,  # Unbalanced system (zero sequence)
-    ]
-    # fmt: on
-    expected_sym_index = pd.MultiIndex.from_arrays(
-        [
-            pd.Index(["bus1", "bus1", "bus1", "bus2", "bus2", "bus2", "bus3", "bus3", "bus3"]),
-            pd.CategoricalIndex(
-                ["zero", "pos", "neg", "zero", "pos", "neg", "zero", "pos", "neg"], dtype=SequenceDtype
-            ),
-        ],
-        names=["bus_id", "sequence"],
-    )
-    voltage = pd.Series(voltage_data, name="voltage").rename_axis(index=["bus_id", "phase"])
-    expected = pd.Series(data=expected_sym_data, index=expected_sym_index, name="voltage")
-    assert_series_equal(series_phasor_to_sym(voltage), expected, check_exact=False)
 
 
 def test_calculate_voltages():
diff --git a/roseau/load_flow/tests/test_deprecations.py b/roseau/load_flow/tests/test_deprecations.py
new file mode 100644
index 00000000..dc73cf31
--- /dev/null
+++ b/roseau/load_flow/tests/test_deprecations.py
@@ -0,0 +1,86 @@
+import importlib
+
+import pytest
+
+from roseau.load_flow import Insulator, Material, constants, converters, utils
+
+
+def test_moved_names():
+    # utils -> constants
+    for name in (
+        "ALPHA",
+        "ALPHA2",
+        "PI",
+        "MU_0",
+        "EPSILON_0",
+        "F",
+        "OMEGA",
+        "RHO",
+        "MU_R",
+        "DELTA_P",
+        "TAN_D",
+        "EPSILON_R",
+    ):
+        with pytest.warns(
+            FutureWarning,
+            match=f"Importing {name} from 'roseau.load_flow.utils' is deprecated. Use 'rlf.constants.{name}' instead.",
+        ):
+            getattr(utils, name)
+    # utils -> sym
+    for name in ("PositiveSequence", "NegativeSequence", "ZeroSequence"):
+        with pytest.warns(
+            FutureWarning,
+            match=f"Importing {name} from 'roseau.load_flow.utils' is deprecated. Use 'rlf.sym.{name}' instead.",
+        ):
+            getattr(utils, name)
+    # utils -> types
+    for name in ("LineType", "Material", "Insulator"):
+        with pytest.warns(
+            FutureWarning,
+            match=f"Importing {name} from 'roseau.load_flow.utils' is deprecated. Use 'rlf.types.{name}' instead.",
+        ):
+            getattr(utils, name)
+    # constants -> sym
+    for name in ("PositiveSequence", "NegativeSequence", "ZeroSequence"):
+        with pytest.warns(
+            FutureWarning,
+            match=f"'rlf.constants.{name}' is deprecated. Use 'rlf.sym.{name}' instead.",
+        ):
+            getattr(constants, name)
+    # converters -> sym
+    for name in ("phasor_to_sym", "sym_to_phasor", "series_phasor_to_sym"):
+        with pytest.warns(
+            FutureWarning,
+            match=f"'rlf.converters.{name}' is deprecated. Use 'rlf.sym.{name}' instead.",
+        ):
+            getattr(converters, name)
+
+
+def test_renamed_classes():
+    with pytest.warns(FutureWarning, match="The `ConductorType` class is deprecated. Use `Material` instead."):
+        from roseau.load_flow import ConductorType
+    assert ConductorType is Material
+    with pytest.warns(FutureWarning, match="The `InsulatorType` class is deprecated. Use `Insulator` instead."):
+        from roseau.load_flow import InsulatorType
+    assert InsulatorType is Insulator
+
+
+def test_non_existent_name():
+    from roseau import load_flow
+
+    for mod in (utils, constants, converters, load_flow):
+        with pytest.raises(AttributeError, match=f"module '{mod.__name__}' has no attribute 'non_existent_name'"):
+            _ = mod.non_existent_name
+
+
+def test_deprecated_utils_modules():
+    with pytest.warns(
+        FutureWarning,
+        match="Module 'roseau.load_flow.utils.constants' is deprecated. Use 'rlf.constants' directly instead.",
+    ):
+        importlib.import_module("roseau.load_flow.utils.constants")
+    with pytest.warns(
+        FutureWarning,
+        match="Module 'roseau.load_flow.utils.types' is deprecated. Use 'rlf.types' directly instead.",
+    ):
+        importlib.import_module("roseau.load_flow.utils.types")
diff --git a/roseau/load_flow/tests/test_sym.py b/roseau/load_flow/tests/test_sym.py
new file mode 100644
index 00000000..40382560
--- /dev/null
+++ b/roseau/load_flow/tests/test_sym.py
@@ -0,0 +1,114 @@
+import numpy as np
+import pandas as pd
+from pandas.testing import assert_series_equal
+
+from roseau.load_flow.sym import phasor_to_sym, series_phasor_to_sym, sym_to_phasor
+from roseau.load_flow.utils import SequenceDtype
+
+
+def test_phasor_to_sym():
+    # Tests verified with https://phillipopambuh.info/portfolio/calculator--symmetrical_components.html
+    va = 230 + 0j
+    vb = 230 * np.e ** (1j * 4 * np.pi / 3)
+    vc = 230 * np.e ** (1j * 2 * np.pi / 3)
+
+    # Test balanced direct system: positive sequence
+    expected = np.array([0, 230, 0], dtype=complex)
+    assert np.allclose(phasor_to_sym([va, vb, vc]), expected)
+    # Also test numpy array input with different shapes
+    assert np.allclose(phasor_to_sym(np.array([va, vb, vc])), expected)
+    assert np.allclose(phasor_to_sym(np.array([[va], [vb], [vc]])), expected.reshape((3, 1)))
+
+    # Test balanced indirect system: negative sequence
+    expected = np.array([0, 0, 230], dtype=complex)
+    assert np.allclose(phasor_to_sym([va, vc, vb]), expected)
+
+    # Test unbalanced system: zero sequence
+    expected = np.array([230, 0, 0], dtype=complex)
+    assert np.allclose(phasor_to_sym([va, va, va]), expected)
+
+    # Test unbalanced system: general case
+    va = 200 + 0j
+    expected = np.array([10 * np.exp(1j * np.pi), 220, 10 * np.exp(1j * np.pi)], dtype=complex)
+    assert np.allclose(phasor_to_sym([va, vb, vc]), expected)
+
+
+def test_sym_to_phasor():
+    # Tests verified with https://phillipopambuh.info/portfolio/calculator--symmetrical_components.html
+    va = 230 + 0j
+    vb = 230 * np.e ** (1j * 4 * np.pi / 3)
+    vc = 230 * np.e ** (1j * 2 * np.pi / 3)
+
+    # Test balanced direct system: positive sequence
+    expected = np.array([va, vb, vc], dtype=complex)
+    assert np.allclose(sym_to_phasor([0, va, 0]), expected)
+    # Also test numpy array input with different shapes
+    assert np.allclose(sym_to_phasor(np.array([0, va, 0])), expected)
+    assert np.allclose(sym_to_phasor(np.array([[0], [va], [0]])), expected.reshape((3, 1)))
+
+    # Test balanced indirect system: negative sequence
+    expected = np.array([va, vc, vb], dtype=complex)
+    assert np.allclose(sym_to_phasor([0, 0, va]), expected)
+
+    # Test unbalanced system: zero sequence
+    expected = np.array([va, va, va], dtype=complex)
+    assert np.allclose(sym_to_phasor([va, 0, 0]), expected)
+
+    # Test unbalanced system: general case
+    va = 200 + 0j
+    expected = np.array([va, vb, vc], dtype=complex)
+    assert np.allclose(sym_to_phasor([10 * np.e ** (1j * np.pi), 220, 10 * np.e ** (1j * np.pi)]), expected)
+
+
+def test_phasor_sym_roundtrip():
+    va = 230 + 0j
+    vb = 230 * np.e ** (1j * 4 * np.pi / 3)
+    vc = 230 * np.e ** (1j * 2 * np.pi / 3)
+
+    # Test balanced direct system: positive sequence
+    assert np.allclose(sym_to_phasor(phasor_to_sym([va, vb, vc])), np.array([va, vb, vc]))
+
+    # Test balanced indirect system: negative sequence
+    assert np.allclose(sym_to_phasor(phasor_to_sym([va, vc, vb])), np.array([va, vc, vb]))
+
+    # Test unbalanced system: zero sequence
+    assert np.allclose(sym_to_phasor(phasor_to_sym([va, va, va])), np.array([va, va, va]))
+
+    # Test unbalanced system: general case
+    va = 200 + 0j
+    assert np.allclose(sym_to_phasor(phasor_to_sym([va, vb, vc])), np.array([va, vb, vc]))
+
+
+def test_series_phasor_to_sym():
+    va = 230 + 0j
+    vb = 230 * np.exp(1j * 4 * np.pi / 3)
+    vc = 230 * np.exp(1j * 2 * np.pi / 3)
+
+    # Test with different phases per bus, different systems, different magnitudes
+    # fmt: off
+    voltage_data = {
+        # Direct system (positive sequence)
+        ("bus1", "a"): va, ("bus1", "b"): vb, ("bus1", "c"): vc,
+        # Indirect system (negative sequence)
+        ("bus2", "an"): va / 2, ("bus2", "bn"): vc / 2, ("bus2", "cn"): vb / 2,
+        # Unbalanced system (zero sequence)
+        ("bus3", "ab"): va, ("bus3", "bc"): va, ("bus3", "ca"): va,
+    }
+    expected_sym_data = [
+        0, va, 0,  # Direct system (positive sequence)
+        0, 0, va / 2,  # Indirect system (negative sequence)
+        va, 0, 0,  # Unbalanced system (zero sequence)
+    ]
+    # fmt: on
+    expected_sym_index = pd.MultiIndex.from_arrays(
+        [
+            pd.Index(["bus1", "bus1", "bus1", "bus2", "bus2", "bus2", "bus3", "bus3", "bus3"]),
+            pd.CategoricalIndex(
+                ["zero", "pos", "neg", "zero", "pos", "neg", "zero", "pos", "neg"], dtype=SequenceDtype
+            ),
+        ],
+        names=["bus_id", "sequence"],
+    )
+    voltage = pd.Series(voltage_data, name="voltage").rename_axis(index=["bus_id", "phase"])
+    expected = pd.Series(data=expected_sym_data, index=expected_sym_index, name="voltage")
+    assert_series_equal(series_phasor_to_sym(voltage), expected, check_exact=False)
diff --git a/roseau/load_flow/utils/tests/test_types.py b/roseau/load_flow/tests/test_types.py
similarity index 82%
rename from roseau/load_flow/utils/tests/test_types.py
rename to roseau/load_flow/tests/test_types.py
index 4a3c6e11..0e7203f3 100644
--- a/roseau/load_flow/utils/tests/test_types.py
+++ b/roseau/load_flow/tests/test_types.py
@@ -1,12 +1,24 @@
 import pytest
 
+from roseau.load_flow.constants import DELTA_P, EPSILON_R, MU_R, RHO, TAN_D
 from roseau.load_flow.exceptions import RoseauLoadFlowException, RoseauLoadFlowExceptionCode
-from roseau.load_flow.utils.types import Insulator, LineType, Material
+from roseau.load_flow.types import Insulator, LineType, Material
 
 TYPES = [Material, Insulator, LineType]
 TYPES_IDS = [x.__name__ for x in TYPES]
 
 
+def test_types_of_constants():
+    for x in Material:
+        assert x in MU_R
+        assert x in RHO
+        assert x in DELTA_P
+
+    for x in Insulator:
+        assert x in TAN_D
+        assert x in EPSILON_R
+
+
 @pytest.mark.parametrize(scope="module", argnames="t", argvalues=TYPES, ids=TYPES_IDS)
 def test_types_basic(t):
     for x in t:
diff --git a/roseau/load_flow/types.py b/roseau/load_flow/types.py
new file mode 100644
index 00000000..a648ecd7
--- /dev/null
+++ b/roseau/load_flow/types.py
@@ -0,0 +1,133 @@
+import logging
+from enum import auto
+
+from roseau.load_flow._compat import StrEnum
+from roseau.load_flow.exceptions import RoseauLoadFlowException, RoseauLoadFlowExceptionCode
+
+# The local logger
+logger = logging.getLogger(__name__)
+
+
+class LineType(StrEnum):
+    """The type of a line."""
+
+    OVERHEAD = auto()
+    """An overhead line that can be vertically or horizontally configured -- Fr = Aérien."""
+    UNDERGROUND = auto()
+    """An underground or a submarine cable -- Fr = Souterrain/Sous-Marin."""
+    TWISTED = auto()
+    """A twisted line commonly known as Aerial Cable or Aerial Bundled Conductor (ABC) -- Fr = Torsadé."""
+
+    # aliases
+    O = OVERHEAD  # noqa: E741
+    U = UNDERGROUND
+    T = TWISTED
+
+    @classmethod
+    def _missing_(cls, value: object) -> "LineType | None":
+        if isinstance(value, str):
+            try:
+                return cls[value.upper()]
+            except KeyError:
+                pass
+        msg = f"{value!r} cannot be converted into a LineType."
+        logger.error(msg)
+        raise RoseauLoadFlowException(msg, RoseauLoadFlowExceptionCode.BAD_LINE_TYPE)
+
+    def code(self) -> str:
+        """A code that can be used in line type names."""
+        return self.name[0]
+
+
+class Material(StrEnum):
+    """The type of the material of the conductor."""
+
+    CU = auto()
+    """Copper -- Fr = Cuivre."""
+    AL = auto()
+    """All Aluminum Conductor (AAC) -- Fr = Aluminium."""
+    AM = auto()
+    """All Aluminum Alloy Conductor (AAAC) -- Fr = Almélec."""
+    AA = auto()
+    """Aluminum Conductor Steel Reinforced (ACSR) -- Fr = Alu-Acier."""
+    LA = auto()
+    """Aluminum Alloy Conductor Steel Reinforced (AACSR) -- Fr = Almélec-Acier."""
+
+    # Aliases
+    AAC = AL  # 1350-H19 (Standard Round of Compact Round)
+    """All Aluminum Conductor (AAC) -- Fr = Aluminium."""
+    # AAC/TW  # 1380-H19 (Trapezoidal Wire)
+
+    AAAC = AM
+    """All Aluminum Alloy Conductor (AAAC) -- Fr = Almélec."""
+    # Aluminum alloy 6201-T81.
+    # Concentric-lay-stranded
+    # conforms to ASTM Specification B-399
+    # Applications: Overhead
+
+    ACSR = AA
+    """Aluminum Conductor Steel Reinforced (ACSR) -- Fr = Alu-Acier."""
+    # Aluminum alloy 1350-H-19
+    # Applications: Bare overhead transmission cable and primary and secondary distribution cable
+
+    AACSR = LA
+    """Aluminum Alloy Conductor Steel Reinforced (AACSR) -- Fr = Almélec-Acier."""
+
+    @classmethod
+    def _missing_(cls, value: object) -> "Material":
+        if isinstance(value, str):
+            try:
+                return cls[value.upper()]
+            except KeyError:
+                pass
+        msg = f"{value!r} cannot be converted into a Material."
+        logger.error(msg)
+        raise RoseauLoadFlowException(msg, RoseauLoadFlowExceptionCode.BAD_MATERIAL)
+
+    def code(self) -> str:
+        """A code that can be used in conductor type names."""
+        return self.name
+
+
+class Insulator(StrEnum):
+    """The type of the insulator for a wire."""
+
+    NONE = auto()
+    """No insulation."""
+
+    # General insulators (IEC 60287)
+    HDPE = auto()
+    """High-Density PolyEthylene (HDPE) insulation."""
+    MDPE = auto()
+    """Medium-Density PolyEthylene (MDPE) insulation."""
+    LDPE = auto()
+    """Low-Density PolyEthylene (LDPE) insulation."""
+    XLPE = auto()
+    """Cross-linked polyethylene (XLPE) insulation."""
+    EPR = auto()
+    """Ethylene-Propylene Rubber (EPR) insulation."""
+    PVC = auto()
+    """PolyVinyl Chloride (PVC) insulation."""
+    IP = auto()
+    """Impregnated Paper (IP) insulation."""
+
+    # Aliases
+    PEX = XLPE
+    """Alias -- Cross-linked polyethylene (XLPE) insulation."""
+    PE = MDPE
+    """Alias -- Medium-Density PolyEthylene (MDPE) insulation."""
+
+    @classmethod
+    def _missing_(cls, value: object) -> "Insulator":
+        if isinstance(value, str):
+            try:
+                return cls[value.upper()]
+            except KeyError:
+                pass
+        msg = f"{value!r} cannot be converted into a Insulator."
+        logger.error(msg)
+        raise RoseauLoadFlowException(msg, RoseauLoadFlowExceptionCode.BAD_INSULATOR)
+
+    def code(self) -> str:
+        """A code that can be used in insulator type names."""
+        return self.name
diff --git a/roseau/load_flow/units.py b/roseau/load_flow/units.py
index fd36893b..e38f9524 100644
--- a/roseau/load_flow/units.py
+++ b/roseau/load_flow/units.py
@@ -21,21 +21,22 @@
 .. _pint: https://pint.readthedocs.io/en/stable/getting/overview.html
 """
 
-from collections.abc import Callable, Iterable, Sequence
+import functools
+from collections.abc import Callable, Iterable, MutableSequence, Sequence
 from decimal import Decimal
 from fractions import Fraction
+from inspect import Parameter, Signature, signature
+from itertools import zip_longest
 from types import GenericAlias
 from typing import TYPE_CHECKING, Any, TypeAlias, TypeVar, overload
 
 import numpy as np
 from numpy.typing import NDArray
-from pint import Unit, UnitRegistry
 from pint.facets.numpy.quantity import NumpyQuantity
-from pint.util import UnitsContainer
+from pint.registry import Unit, UnitRegistry
+from pint.util import UnitsContainer, to_units_container
 
-from roseau.load_flow._wrapper import wraps
-
-FuncT = TypeVar("FuncT", bound=Callable)
+__all__ = ["ureg", "Q_", "ureg_wraps"]
 
 ureg: UnitRegistry = UnitRegistry(
     preprocessors=[
@@ -101,12 +102,11 @@ def __getattr__(self, name: str) -> Any: ...  # attributes of the magnitude are
     ureg.Quantity.__class_getitem__ = classmethod(GenericAlias)
     globals()["Q_"] = ureg.Quantity  # Use globals() to trick PyCharm
 
+FuncT = TypeVar("FuncT", bound=Callable)
+OptionalUnits: TypeAlias = str | Unit | None | tuple[str | Unit | None, ...] | list[str | Unit | None]
+
 
-def ureg_wraps(
-    ret: str | Unit | None | Iterable[str | Unit | None],
-    args: str | Unit | None | Iterable[str | Unit | None],
-    strict: bool = True,
-) -> Callable[[FuncT], FuncT]:
+def ureg_wraps(ret: OptionalUnits, args: OptionalUnits, strict: bool = True) -> Callable[[FuncT], FuncT]:
     """Wraps a function to become pint-aware.
 
     Args:
@@ -120,3 +120,135 @@ def ureg_wraps(
             Indicates that only quantities are accepted. (Default value = True)
     """
     return wraps(ureg, ret, args)
+
+
+def _parse_wrap_args(args: Iterable[str | Unit | None]) -> Callable:
+    """Create a converter function for the wrapper"""
+    # _to_units_container
+    args_as_uc = [to_units_container(arg) for arg in args]
+
+    # Check for references in args, remove None values
+    unit_args_ndx = {ndx for ndx, arg in enumerate(args_as_uc) if arg is not None}
+
+    def _converter(ureg: "UnitRegistry", sig: "Signature", values: "list[Any]", kw: "dict[Any]"):
+        len_initial_values = len(values)
+
+        # pack kwargs
+        for i, param_name in enumerate(sig.parameters):
+            if i >= len_initial_values:
+                values.append(kw[param_name])
+
+        # convert arguments
+        for ndx in unit_args_ndx:
+            value = values[ndx]
+            if isinstance(value, ureg.Quantity):
+                values[ndx] = ureg.convert(value.magnitude, value.units, args_as_uc[ndx])
+            elif isinstance(value, MutableSequence):
+                for i, val in enumerate(value):
+                    if isinstance(val, ureg.Quantity):
+                        value[i] = ureg.convert(val.magnitude, val.units, args_as_uc[ndx])
+
+        # unpack kwargs
+        for i, param_name in enumerate(sig.parameters):
+            if i >= len_initial_values:
+                kw[param_name] = values[i]
+
+        return values[:len_initial_values], kw
+
+    return _converter
+
+
+def _apply_defaults(sig: Signature, args: tuple[Any], kwargs: dict[str, Any]) -> tuple[list[Any], dict[str, Any]]:
+    """Apply default keyword arguments.
+
+    Named keywords may have been left blank. This function applies the default
+    values so that every argument is defined.
+    """
+    n = len(args)
+    for i, param in enumerate(sig.parameters.values()):
+        if i >= n and param.default != Parameter.empty and param.name not in kwargs:
+            kwargs[param.name] = param.default
+    return list(args), kwargs
+
+
+def wraps(ureg: UnitRegistry, ret: OptionalUnits, args: OptionalUnits) -> Callable[[FuncT], FuncT]:
+    """Wraps a function to become pint-aware.
+
+    Use it when a function requires a numerical value but in some specific
+    units. The wrapper function will take a pint quantity, convert to the units
+    specified in `args` and then call the wrapped function with the resulting
+    magnitude.
+
+    The value returned by the wrapped function will be converted to the units
+    specified in `ret`.
+
+    Args:
+        ureg:
+            A UnitRegistry instance.
+
+        ret:
+            Units of each of the return values. Use `None` to skip argument conversion.
+
+        args:
+             Units of each of the input arguments. Use `None` to skip argument conversion.
+
+    Returns:
+        The wrapper function.
+
+    Raises:
+        TypeError
+            if the number of given arguments does not match the number of function parameters.
+            if any of the provided arguments is not a unit a string or Quantity
+    """
+    if not isinstance(args, (list, tuple)):
+        args = (args,)
+
+    for arg in args:
+        if arg is not None and not isinstance(arg, (ureg.Unit, str)):
+            raise TypeError(f"wraps arguments must by of type str or Unit, not {type(arg)} ({arg})")
+
+    converter = _parse_wrap_args(args)
+
+    is_ret_container = isinstance(ret, (list, tuple))
+    if is_ret_container:
+        for arg in ret:
+            if arg is not None and not isinstance(arg, (ureg.Unit, str)):
+                raise TypeError(f"wraps 'ret' argument must by of type str or Unit, not {type(arg)} ({arg})")
+        ret = ret.__class__([to_units_container(arg, ureg) for arg in ret])
+    else:
+        if ret is not None and not isinstance(ret, (ureg.Unit, str)):
+            raise TypeError(f"wraps 'ret' argument must by of type str or Unit, not {type(ret)} ({ret})")
+        ret = to_units_container(ret, ureg)
+
+    def decorator(func):
+        sig = signature(func)
+        count_params = len(sig.parameters)
+        if len(args) != count_params:
+            raise TypeError(f"{func.__name__} takes {count_params} parameters, but {len(args)} units were passed")
+
+        assigned = tuple(attr for attr in functools.WRAPPER_ASSIGNMENTS if hasattr(func, attr))
+        updated = tuple(attr for attr in functools.WRAPPER_UPDATES if hasattr(func, attr))
+
+        @functools.wraps(func, assigned=assigned, updated=updated)
+        def wrapper(*values, **kw):
+            values, kw = _apply_defaults(sig, values, kw)
+
+            # In principle, the values are used as is
+            # When then extract the magnitudes when needed.
+            new_values, new_kw = converter(ureg, sig, values, kw)
+
+            result = func(*new_values, **new_kw)
+
+            if is_ret_container:
+                return ret.__class__(
+                    res if unit is None else ureg.Quantity(res, unit) for unit, res in zip_longest(ret, result)
+                )
+
+            if ret is None:
+                return result
+
+            return ureg.Quantity(result, ret)
+
+        return wrapper
+
+    return decorator
diff --git a/roseau/load_flow/utils/__init__.py b/roseau/load_flow/utils/__init__.py
index d3a09ab6..efb42a2f 100644
--- a/roseau/load_flow/utils/__init__.py
+++ b/roseau/load_flow/utils/__init__.py
@@ -1,67 +1,86 @@
 """
-This module contains utility classes and functions for Roseau Load Flow.
+Internal utility classes and functions for Roseau Load Flow.
 """
 
-from roseau.load_flow.utils.constants import (
-    ALPHA,
-    ALPHA2,
-    DELTA_P,
-    EPSILON_0,
-    EPSILON_R,
-    MU_0,
-    MU_R,
-    OMEGA,
-    PI,
-    RHO,
-    SQRT3,
-    TAN_D,
-    F,
-    NegativeSequence,
-    PositiveSequence,
-    ZeroSequence,
-)
-from roseau.load_flow.utils.mixins import CatalogueMixin, Identifiable, JsonMixin
-from roseau.load_flow.utils.types import (
+from roseau.load_flow.utils.dtypes import (
+    DTYPES,
     BranchTypeDtype,
-    Insulator,
-    LineType,
     LoadTypeDtype,
-    Material,
     PhaseDtype,
     SequenceDtype,
     VoltagePhaseDtype,
 )
+from roseau.load_flow.utils.exceptions import find_stack_level
+from roseau.load_flow.utils.log import set_logging_config
+from roseau.load_flow.utils.mixins import CatalogueMixin, Identifiable, JsonMixin
+from roseau.load_flow.utils.versions import show_versions
 
 __all__ = [
-    # Constants
-    "DELTA_P",
-    "EPSILON_0",
-    "EPSILON_R",
-    "F",
-    "MU_0",
-    "MU_R",
-    "OMEGA",
-    "PI",
-    "RHO",
-    "SQRT3",
-    "TAN_D",
-    "ALPHA",
-    "ALPHA2",
-    "PositiveSequence",
-    "NegativeSequence",
-    "ZeroSequence",
     # Mixins
     "Identifiable",
     "JsonMixin",
     "CatalogueMixin",
-    # Types
-    "LineType",
-    "Material",
-    "Insulator",
-    # Dtypes
-    "PhaseDtype",
-    "VoltagePhaseDtype",
+    # Exceptions and warnings
+    "find_stack_level",
+    # DTypes
+    "DTYPES",
     "BranchTypeDtype",
     "LoadTypeDtype",
+    "PhaseDtype",
     "SequenceDtype",
+    "VoltagePhaseDtype",
+    # Versions
+    "show_versions",
+    # Logging
+    "set_logging_config",
 ]
+
+
+def __getattr__(name: str):
+    import warnings
+
+    deprecation_template = (
+        "Importing {name} from 'roseau.load_flow.utils' is deprecated. Use 'rlf.{module}.{name}' instead."
+    )
+    if name in (
+        "ALPHA",
+        "ALPHA2",
+        "PI",
+        "SQRT3",
+        "DELTA_P",
+        "EPSILON_0",
+        "EPSILON_R",
+        "F",
+        "MU_0",
+        "MU_R",
+        "OMEGA",
+        "RHO",
+        "TAN_D",
+    ):
+        warnings.warn(
+            deprecation_template.format(name=name, module="constants"),
+            category=FutureWarning,
+            stacklevel=find_stack_level(),
+        )
+        from roseau.load_flow import constants
+
+        return getattr(constants, name)
+    elif name in ("PositiveSequence", "NegativeSequence", "ZeroSequence"):
+        warnings.warn(
+            deprecation_template.format(name=name, module="sym"),
+            category=FutureWarning,
+            stacklevel=find_stack_level(),
+        )
+        from roseau.load_flow import sym
+
+        return getattr(sym, name)
+    elif name in ("LineType", "Material", "Insulator"):
+        warnings.warn(
+            deprecation_template.format(name=name, module="types"),
+            category=FutureWarning,
+            stacklevel=find_stack_level(),
+        )
+        from roseau.load_flow import types
+
+        return getattr(types, name)
+    raise AttributeError(f"module '{__name__}' has no attribute '{name}'")
diff --git a/roseau/load_flow/utils/constants.py b/roseau/load_flow/utils/constants.py
index 1db870ab..44bd90a5 100644
--- a/roseau/load_flow/utils/constants.py
+++ b/roseau/load_flow/utils/constants.py
@@ -1,101 +1,10 @@
-import cmath
-import math
-from typing import Final
-
-import numpy as np
-
-from roseau.load_flow.units import Q_
-from roseau.load_flow.utils.types import Insulator, Material
-
-SQRT3: Final = math.sqrt(3)
-"""The square root of 3."""
-
-PI: Final = math.pi
-"""The famous mathematical constant :math:`\\pi = 3.141592\\ldots`."""
-
-ALPHA: Final = cmath.exp(2 / 3 * PI * 1j)
-"""complex: Phasor rotation operator :math:`\\alpha`, which rotates a phasor vector counterclockwise
-by 120 degrees when multiplied by it."""
-
-ALPHA2: Final = ALPHA**2
-"""complex: Phasor rotation operator :math:`\\alpha^2`, which rotates a phasor vector clockwise by
-120 degrees when multiplied by it."""
-
-PositiveSequence: Final = np.array([1, ALPHA2, ALPHA], dtype=np.complex128)
-"""numpy.ndarray[complex]: Unit positive sequence components of a three-phase system."""
-NegativeSequence: Final = np.array([1, ALPHA, ALPHA2], dtype=np.complex128)
-"""numpy.ndarray[complex]: Unit negative sequence components of a three-phase system."""
-ZeroSequence: Final = np.array([1, 1, 1], dtype=np.complex128)
-"""numpy.ndarray[complex]: Unit zero sequence components of a three-phase system."""
-
-MU_0: Final = Q_(1.25663706212e-6, "H/m")
-"""Magnetic permeability of the vacuum :math:`\\mu_0 = 4 \\pi \\times 10^{-7}` (H/m)."""
-
-EPSILON_0: Final = Q_(8.8541878128e-12, "F/m")
-"""Vacuum permittivity :math:`\\varepsilon_0 =  8.8541878128 \\times 10^{-12}` (F/m)."""
-
-F: Final = Q_(50.0, "Hz")
-"""Network frequency :math:`f = 50` (Hz)."""
-
-OMEGA: Final = Q_(2 * PI * F, "rad/s")
-"""Angular frequency :math:`\\omega = 2 \\pi f` (rad/s)."""
-
-RHO: Final[dict[Material, Q_[float]]] = {
-    Material.CU: Q_(1.7241e-8, "ohm*m"),  # IEC 60287-1-1 Table 1
-    Material.AL: Q_(2.8264e-8, "ohm*m"),  # IEC 60287-1-1 Table 1
-    Material.AM: Q_(3.26e-8, "ohm*m"),  # verified
-    Material.AA: Q_(4.0587e-8, "ohm*m"),  # verified (approx. AS 3607 ACSR/GZ)
-    Material.LA: Q_(3.26e-8, "ohm*m"),
-}
-"""Resistivity of common conductor materials (Ohm.m)."""
-
-MU_R: Final[dict[Material, Q_[float]]] = {
-    Material.CU: Q_(0.9999935849131266),
-    Material.AL: Q_(1.0000222328028834),
-    Material.AM: Q_(0.9999705074463784),
-    Material.AA: Q_(1.0000222328028834),  # ==AL
-    Material.LA: Q_(0.9999705074463784),  # ==AM
-}
-"""Relative magnetic permeability of common conductor materials."""
-
-DELTA_P: Final[dict[Material, Q_[float]]] = {
-    Material.CU: Q_(9.33, "mm"),
-    Material.AL: Q_(11.95, "mm"),
-    Material.AM: Q_(12.85, "mm"),
-    Material.AA: Q_(14.34, "mm"),
-    Material.LA: Q_(12.85, "mm"),
-}
-"""Skin depth of common conductor materials :math:`\\sqrt{\\dfrac{\\rho}{\\pi f \\mu_r \\mu_0}}` (mm)."""
-# Skin depth is the depth at which the current density is reduced to 1/e (~37%) of the surface value.
-# Generated with:
-# ---------------
-# def delta_p(rho, mu_r):
-#     return np.sqrt(rho / (PI * F * mu_r * MU_0))
-# for material in Material:
-#     print(material, delta_p(RHO[material], MU_R[material]).m_as("mm"))
-
-TAN_D: Final[dict[Insulator, Q_[float]]] = {
-    Insulator.PVC: Q_(1000e-4),
-    Insulator.HDPE: Q_(10e-4),
-    Insulator.MDPE: Q_(10e-4),
-    Insulator.LDPE: Q_(10e-4),
-    Insulator.XLPE: Q_(40e-4),
-    Insulator.EPR: Q_(200e-4),
-    Insulator.IP: Q_(100e-4),
-    Insulator.NONE: Q_(0),
-}
-"""Loss angles of common insulator materials according to the IEC 60287 standard."""
-# IEC 60287-1-1 Table 3. We only include the MV values.
-
-EPSILON_R: Final[dict[Insulator, Q_[float]]] = {
-    Insulator.PVC: Q_(8.0),
-    Insulator.HDPE: Q_(2.3),
-    Insulator.MDPE: Q_(2.3),
-    Insulator.LDPE: Q_(2.3),
-    Insulator.XLPE: Q_(2.5),
-    Insulator.EPR: Q_(3.0),
-    Insulator.IP: Q_(4.0),
-    Insulator.NONE: Q_(1.0),
-}
-"""Relative permittivity of common insulator materials according to the IEC 60287 standard."""
-# IEC 60287-1-1 Table 3. We only include the MV values.
+import warnings
+
+warnings.warn(
+    "Module 'roseau.load_flow.utils.constants' is deprecated. Use 'rlf.constants' directly instead.",
+    category=FutureWarning,
+    stacklevel=2,
+)
+# ruff: noqa: E402, F401
+from roseau.load_flow.constants import DELTA_P, EPSILON_0, EPSILON_R, MU_0, MU_R, OMEGA, PI, RHO, TAN_D, F
+from roseau.load_flow.sym import ALPHA, ALPHA2, NegativeSequence, PositiveSequence, ZeroSequence
diff --git a/roseau/load_flow/utils/_doc_utils.py b/roseau/load_flow/utils/doc_utils.py
similarity index 100%
rename from roseau/load_flow/utils/_doc_utils.py
rename to roseau/load_flow/utils/doc_utils.py
diff --git a/roseau/load_flow/utils/dtypes.py b/roseau/load_flow/utils/dtypes.py
new file mode 100644
index 00000000..d5690375
--- /dev/null
+++ b/roseau/load_flow/utils/dtypes.py
@@ -0,0 +1,55 @@
+from typing import Final
+
+import pandas as pd
+
+# pandas dtypes used in the data frames
+PhaseDtype: Final = pd.CategoricalDtype(categories=["a", "b", "c", "n"], ordered=True)
+"""Categorical data type used for the phase of potentials, currents, powers, etc."""
+VoltagePhaseDtype: Final = pd.CategoricalDtype(categories=["an", "bn", "cn", "ab", "bc", "ca"], ordered=True)
+"""Categorical data type used for the phase of voltages and flexible powers only."""
+BranchTypeDtype: Final = pd.CategoricalDtype(categories=["line", "transformer", "switch"], ordered=True)
+"""Categorical data type used for branch types."""
+LoadTypeDtype: Final = pd.CategoricalDtype(categories=["power", "current", "impedance"], ordered=True)
+"""Categorical data type used for load types."""
+SequenceDtype: Final = pd.CategoricalDtype(categories=["zero", "pos", "neg"], ordered=True)
+"""Categorical data type used for symmetrical components."""
+DTYPES: Final = {
+    "bus_id": object,
+    "branch_id": object,
+    "transformer_id": object,
+    "line_id": object,
+    "switch_id": object,
+    "load_id": object,
+    "source_id": object,
+    "ground_id": object,
+    "potential_ref_id": object,
+    "type": object,
+    "phase": PhaseDtype,
+    "current": complex,
+    "current1": complex,
+    "current2": complex,
+    "power": complex,
+    "flexible_power": complex,
+    "power1": complex,
+    "power2": complex,
+    "potential": complex,
+    "potential1": complex,
+    "potential2": complex,
+    "voltage": complex,
+    "voltage1": complex,
+    "voltage2": complex,
+    "series_losses": complex,
+    "shunt_losses": complex,
+    "series_current": complex,
+    "max_current": float,
+    "loading": float,
+    "max_loading": float,
+    "sn": float,
+    "ampacity": float,
+    "voltage_level": float,
+    "nominal_voltage": float,
+    "min_voltage_level": float,
+    "max_voltage_level": float,
+    "violated": pd.BooleanDtype(),
+    "flexible": pd.BooleanDtype(),
+}
diff --git a/roseau/load_flow/utils/_exceptions.py b/roseau/load_flow/utils/exceptions.py
similarity index 100%
rename from roseau/load_flow/utils/_exceptions.py
rename to roseau/load_flow/utils/exceptions.py
diff --git a/roseau/load_flow/utils/_optional_deps.py b/roseau/load_flow/utils/optional_deps.py
similarity index 100%
rename from roseau/load_flow/utils/_optional_deps.py
rename to roseau/load_flow/utils/optional_deps.py
diff --git a/roseau/load_flow/utils/tests/test_constants.py b/roseau/load_flow/utils/tests/test_constants.py
deleted file mode 100644
index 09b2cdc8..00000000
--- a/roseau/load_flow/utils/tests/test_constants.py
+++ /dev/null
@@ -1,13 +0,0 @@
-from roseau.load_flow.utils.constants import DELTA_P, EPSILON_R, MU_R, RHO, TAN_D
-from roseau.load_flow.utils.types import Insulator, Material
-
-
-def test_constants():
-    for x in Material:
-        assert x in MU_R
-        assert x in RHO
-        assert x in DELTA_P
-
-    for x in Insulator:
-        assert x in TAN_D
-        assert x in EPSILON_R
diff --git a/roseau/load_flow/utils/tests/test_optional_deps.py b/roseau/load_flow/utils/tests/test_optional_deps.py
index 63228654..1df02ddc 100644
--- a/roseau/load_flow/utils/tests/test_optional_deps.py
+++ b/roseau/load_flow/utils/tests/test_optional_deps.py
@@ -1,20 +1,20 @@
 import pytest
 
 from roseau.load_flow.exceptions import RoseauLoadFlowException, RoseauLoadFlowExceptionCode
-from roseau.load_flow.utils import _optional_deps
+from roseau.load_flow.utils import optional_deps
 
 
 def test_optional_deps(monkeypatch):
     # Working deps
-    plt = _optional_deps.pyplot
+    plt = optional_deps.pyplot
     assert plt.__class__.__name__ == "module"
-    networkx = _optional_deps.networkx
+    networkx = optional_deps.networkx
     assert networkx.__class__.__name__ == "module"
 
     # Fail
     with pytest.raises(AttributeError) as e:
-        _optional_deps.toto  # noqa: B018
-    assert e.value.args[0] == "module roseau.load_flow.utils._optional_deps has no attribute 'toto'"
+        optional_deps.toto  # noqa: B018
+    assert e.value.args[0] == "module roseau.load_flow.utils.optional_deps has no attribute 'toto'"
 
     # Module not installed
     def fake_import(name, *args, **kwargs):
@@ -24,7 +24,7 @@ def fake_import(name, *args, **kwargs):
         m.setattr("builtins.__import__", fake_import)
 
         with pytest.raises(RoseauLoadFlowException) as e:
-            _optional_deps.pyplot  # noqa: B018
+            optional_deps.pyplot  # noqa: B018
         assert (
             e.value.msg == 'matplotlib is required for plotting. Install it with the "plot" extra using '
             '`pip install -U "roseau-load-flow[plot]"`'
@@ -32,7 +32,7 @@ def fake_import(name, *args, **kwargs):
         assert e.value.code == RoseauLoadFlowExceptionCode.IMPORT_ERROR
 
         with pytest.raises(RoseauLoadFlowException) as e:
-            _optional_deps.networkx  # noqa: B018
+            optional_deps.networkx  # noqa: B018
         assert (
             e.value.msg == 'networkx is not installed. Install it with the "graph" extra using `pip install -U '
             '"roseau-load-flow[graph]"`'
diff --git a/roseau/load_flow/utils/tests/test_versions.py b/roseau/load_flow/utils/tests/test_versions.py
index 646d5e24..03ec57ce 100644
--- a/roseau/load_flow/utils/tests/test_versions.py
+++ b/roseau/load_flow/utils/tests/test_versions.py
@@ -1,7 +1,7 @@
 import re
 
 from roseau.load_flow import show_versions
-from roseau.load_flow.utils._versions import _get_dependency_info, _get_sys_info
+from roseau.load_flow.utils.versions import _get_dependency_info, _get_sys_info
 
 
 def test_versions(capsys):
diff --git a/roseau/load_flow/utils/types.py b/roseau/load_flow/utils/types.py
index 109f95e0..6353ccba 100644
--- a/roseau/load_flow/utils/types.py
+++ b/roseau/load_flow/utils/types.py
@@ -1,189 +1,11 @@
-import logging
-from enum import auto
-from typing import Final
-
-import pandas as pd
-
-from roseau.load_flow._compat import StrEnum
-from roseau.load_flow.exceptions import RoseauLoadFlowException, RoseauLoadFlowExceptionCode
-
-# The local logger
-logger = logging.getLogger(__name__)
-
-
-# pandas dtypes used in the data frames
-PhaseDtype = pd.CategoricalDtype(categories=["a", "b", "c", "n"], ordered=True)
-"""Categorical data type used for the phase of potentials, currents, powers, etc."""
-VoltagePhaseDtype = pd.CategoricalDtype(categories=["an", "bn", "cn", "ab", "bc", "ca"], ordered=True)
-"""Categorical data type used for the phase of voltages and flexible powers only."""
-BranchTypeDtype = pd.CategoricalDtype(categories=["line", "transformer", "switch"], ordered=True)
-"""Categorical data type used for branch types."""
-LoadTypeDtype = pd.CategoricalDtype(categories=["power", "current", "impedance"], ordered=True)
-"""Categorical data type used for load types."""
-SequenceDtype = pd.CategoricalDtype(categories=["zero", "pos", "neg"], ordered=True)
-"""Categorical data type used for symmetrical components."""
-_DTYPES: Final = {
-    "bus_id": object,
-    "branch_id": object,
-    "transformer_id": object,
-    "line_id": object,
-    "switch_id": object,
-    "load_id": object,
-    "source_id": object,
-    "ground_id": object,
-    "potential_ref_id": object,
-    "type": object,
-    "phase": PhaseDtype,
-    "current": complex,
-    "current1": complex,
-    "current2": complex,
-    "power": complex,
-    "flexible_power": complex,
-    "power1": complex,
-    "power2": complex,
-    "potential": complex,
-    "potential1": complex,
-    "potential2": complex,
-    "voltage": complex,
-    "voltage1": complex,
-    "voltage2": complex,
-    "series_losses": complex,
-    "shunt_losses": complex,
-    "series_current": complex,
-    "max_current": float,
-    "loading": float,
-    "max_loading": float,
-    "sn": float,
-    "ampacity": float,
-    "voltage_level": float,
-    "nominal_voltage": float,
-    "min_voltage_level": float,
-    "max_voltage_level": float,
-    "violated": pd.BooleanDtype(),
-    "flexible": pd.BooleanDtype(),
-}
-
-
-class LineType(StrEnum):
-    """The type of a line."""
-
-    OVERHEAD = auto()
-    """An overhead line that can be vertically or horizontally configured -- Fr = Aérien."""
-    UNDERGROUND = auto()
-    """An underground or a submarine cable -- Fr = Souterrain/Sous-Marin."""
-    TWISTED = auto()
-    """A twisted line commonly known as Aerial Cable or Aerial Bundled Conductor (ABC) -- Fr = Torsadé."""
-
-    # aliases
-    O = OVERHEAD  # noqa: E741
-    U = UNDERGROUND
-    T = TWISTED
-
-    @classmethod
-    def _missing_(cls, value: object) -> "LineType | None":
-        if isinstance(value, str):
-            try:
-                return cls[value.upper()]
-            except KeyError:
-                pass
-        msg = f"{value!r} cannot be converted into a LineType."
-        logger.error(msg)
-        raise RoseauLoadFlowException(msg, RoseauLoadFlowExceptionCode.BAD_LINE_TYPE)
-
-    def code(self) -> str:
-        """A code that can be used in line type names."""
-        return self.name[0]
-
-
-class Material(StrEnum):
-    """The type of the material of the conductor."""
-
-    CU = auto()
-    """Copper -- Fr = Cuivre."""
-    AL = auto()
-    """All Aluminum Conductor (AAC) -- Fr = Aluminium."""
-    AM = auto()
-    """All Aluminum Alloy Conductor (AAAC) -- Fr = Almélec."""
-    AA = auto()
-    """Aluminum Conductor Steel Reinforced (ACSR) -- Fr = Alu-Acier."""
-    LA = auto()
-    """Aluminum Alloy Conductor Steel Reinforced (AACSR) -- Fr = Almélec-Acier."""
-
-    # Aliases
-    AAC = AL  # 1350-H19 (Standard Round of Compact Round)
-    """All Aluminum Conductor (AAC) -- Fr = Aluminium."""
-    # AAC/TW  # 1380-H19 (Trapezoidal Wire)
-
-    AAAC = AM
-    """All Aluminum Alloy Conductor (AAAC) -- Fr = Almélec."""
-    # Aluminum alloy 6201-T81.
-    # Concentric-lay-stranded
-    # conforms to ASTM Specification B-399
-    # Applications: Overhead
-
-    ACSR = AA
-    """Aluminum Conductor Steel Reinforced (ACSR) -- Fr = Alu-Acier."""
-    # Aluminum alloy 1350-H-19
-    # Applications: Bare overhead transmission cable and primary and secondary distribution cable
-
-    AACSR = LA
-    """Aluminum Alloy Conductor Steel Reinforced (AACSR) -- Fr = Almélec-Acier."""
-
-    @classmethod
-    def _missing_(cls, value: object) -> "Material":
-        if isinstance(value, str):
-            try:
-                return cls[value.upper()]
-            except KeyError:
-                pass
-        msg = f"{value!r} cannot be converted into a Material."
-        logger.error(msg)
-        raise RoseauLoadFlowException(msg, RoseauLoadFlowExceptionCode.BAD_MATERIAL)
-
-    def code(self) -> str:
-        """A code that can be used in conductor type names."""
-        return self.name
-
-
-class Insulator(StrEnum):
-    """The type of the insulator for a wire."""
-
-    NONE = auto()
-    """No insulation."""
-
-    # General insulators (IEC 60287)
-    HDPE = auto()
-    """High-Density PolyEthylene (HDPE) insulation."""
-    MDPE = auto()
-    """Medium-Density PolyEthylene (MDPE) insulation."""
-    LDPE = auto()
-    """Low-Density PolyEthylene (LDPE) insulation."""
-    XLPE = auto()
-    """Cross-linked polyethylene (XLPE) insulation."""
-    EPR = auto()
-    """Ethylene-Propylene Rubber (EPR) insulation."""
-    PVC = auto()
-    """PolyVinyl Chloride (PVC) insulation."""
-    IP = auto()
-    """Impregnated Paper (IP) insulation."""
-
-    # Aliases
-    PEX = XLPE
-    """Alias -- Cross-linked polyethylene (XLPE) insulation."""
-    PE = MDPE
-    """Alias -- Medium-Density PolyEthylene (MDPE) insulation."""
-
-    @classmethod
-    def _missing_(cls, value: object) -> "Insulator":
-        if isinstance(value, str):
-            try:
-                return cls[value.upper()]
-            except KeyError:
-                pass
-        msg = f"{value!r} cannot be converted into a Insulator."
-        logger.error(msg)
-        raise RoseauLoadFlowException(msg, RoseauLoadFlowExceptionCode.BAD_INSULATOR)
-
-    def code(self) -> str:
-        """A code that can be used in insulator type names."""
-        return self.name
+import warnings
+
+warnings.warn(
+    "Module 'roseau.load_flow.utils.types' is deprecated. Use 'rlf.types' directly instead.",
+    category=FutureWarning,
+    stacklevel=2,
+)
+# ruff: noqa: E402, F401
+from roseau.load_flow.types import Insulator, LineType, Material
+from roseau.load_flow.utils.dtypes import DTYPES as _DTYPES
+from roseau.load_flow.utils.dtypes import BranchTypeDtype, LoadTypeDtype, PhaseDtype, SequenceDtype, VoltagePhaseDtype
diff --git a/roseau/load_flow/utils/_versions.py b/roseau/load_flow/utils/versions.py
similarity index 100%
rename from roseau/load_flow/utils/_versions.py
rename to roseau/load_flow/utils/versions.py