diff --git a/bluemira/equilibria/flux_surfaces.py b/bluemira/equilibria/flux_surfaces.py
index 14a5a90bee..5b618092ab 100644
--- a/bluemira/equilibria/flux_surfaces.py
+++ b/bluemira/equilibria/flux_surfaces.py
@@ -449,7 +449,7 @@ def __init__(self, coords: Coordinates):
 
         self.alpha = None
 
-    def clip(self, first_wall: Coordinates):
+    def clip(self, first_wall: Coordinates) -> None:
         """
         Clip the PartialOpenFluxSurface to a first wall.
 
diff --git a/bluemira/radiation_transport/advective_transport.py b/bluemira/radiation_transport/advective_transport.py
index d99233bee4..a4d276ed17 100644
--- a/bluemira/radiation_transport/advective_transport.py
+++ b/bluemira/radiation_transport/advective_transport.py
@@ -8,21 +8,25 @@
 A simplified 2-D solver for calculating charged particle heat loads.
 """
 
-from copy import deepcopy
 from dataclasses import dataclass, fields
 
 import matplotlib.pyplot as plt
 import numpy as np
 from matplotlib.axes import Axes
+from numpy import typing as npt
 
 import bluemira.radiation_transport.flux_surfaces_maker as fsm
 from bluemira.base.constants import EPS
 from bluemira.base.look_and_feel import bluemira_warn
 from bluemira.display.plotter import Zorder, plot_coordinates
+from bluemira.equilibria.flux_surfaces import PartialOpenFluxSurface
 from bluemira.geometry.coordinates import Coordinates
 from bluemira.geometry.plane import BluemiraPlane
 from bluemira.radiation_transport.error import AdvectionTransportError
-from bluemira.radiation_transport.flux_surfaces_maker import _clip_flux_surfaces
+from bluemira.radiation_transport.flux_surfaces_maker import (
+    _clip_flux_surfaces,
+    _process_first_wall,
+)
 
 __all__ = ["ChargedParticleSolver"]
 
@@ -72,15 +76,17 @@ def __init__(
         self._o_point = o_points[0]
         z = self._o_point.z
         self._yz_plane = BluemiraPlane.from_3_points([0, 0, z], [1, 0, z], [1, 1, z])
+        self._process_first_wall = staticmethod(_process_first_wall)
 
     @property
-    def flux_surfaces(self):
+    def flux_surfaces(self) -> list[PartialOpenFluxSurface]:
         """
         All flux surfaces in the ChargedParticleSolver.
 
         Returns
         -------
-        flux_surfaces: List[PartialOpenFluxSurface]
+        flux_surfaces:
+            The list of partially open flux surfaces.
         """
         flux_surfaces = []
         for group in [
@@ -127,27 +133,30 @@ def _check_params(self):
             )
 
     @staticmethod
-    def _process_first_wall(first_wall):
-        """
-        Force working first wall geometry to be closed and counter-clockwise.
-        """
-        first_wall = deepcopy(first_wall)
-
-        if not first_wall.check_ccw(axis=[0, 1, 0]):
-            bluemira_warn(
-                "First wall should be oriented counter-clockwise. Reversing it."
-            )
-            first_wall.reverse()
-
-        if not first_wall.closed:
-            bluemira_warn("First wall should be a closed geometry. Closing it.")
-            first_wall.close()
-        return first_wall
-
-    @staticmethod
-    def _get_arrays(flux_surfaces):
+    def _get_arrays(
+        flux_surfaces,
+    ) -> tuple[
+        npt.NDArray[float],
+        npt.NDArray[float],
+        npt.NDArray[float],
+        npt.NDArray[float],
+        npt.NDArray[float],
+    ]:
         """
         Get arrays of flux surface values.
+
+        Returns
+        -------
+        x_mp: npt.NDArray[float]
+            the array of mid-plane intersection point x-coordinate for each flux surface
+        z_mp: npt.NDArray[float]
+            the array of mid-plane intersection point z-coordinate for each flux surface
+        x_fw: npt.NDArray[float]
+            the array of first-wall intersection point x-coordinate for each flux surface
+        z_fw: npt.NDArray[float]
+            the array of first-wall intersection point z-coordinate for each flux surface
+        alpha: npt.NDArray[float]
+            the array of alpha angle for each flux surface
         """
         x_mp = np.array([fs.x_start for fs in flux_surfaces])
         z_mp = np.array([fs.z_start for fs in flux_surfaces])
@@ -156,7 +165,7 @@ def _get_arrays(flux_surfaces):
         alpha = np.array([fs.alpha for fs in flux_surfaces])
         return x_mp, z_mp, x_fw, z_fw, alpha
 
-    def _make_flux_surfaces_ob(self):
+    def _make_flux_surfaces_ob(self) -> None:
         """
         Make the flux surfaces on the outboard.
         """
@@ -249,6 +258,16 @@ def analyse(self, first_wall: Coordinates):
     def _analyse_SN(self):
         """
         Calculation for the case of single nulls.
+
+        Returns
+        -------
+        x: np.array
+            The x coordinates of the flux surface intersections
+        z: np.array
+            The z coordinates of the flux surface intersections
+        heat_flux: np.array
+            The perpendicular heat fluxes at the intersection points [MW/m^2]
+
         """
         self._make_flux_surfaces_ob()
 
@@ -296,6 +315,15 @@ def _analyse_SN(self):
     def _analyse_DN(self):
         """
         Calculation for the case of double nulls.
+
+        Returns
+        -------
+        x: np.array
+            The x coordinates of the flux surface intersections
+        z: np.array
+            The z coordinates of the flux surface intersections
+        heat_flux: np.array
+            The perpendicular heat fluxes at the intersection points [MW/m^2]
         """
         self._make_flux_surfaces_ob()
         self._make_flux_surfaces_ib()
@@ -411,6 +439,11 @@ def _analyse_DN(self):
     def _q_par(self, x, dx, B, Bp, *, outboard=True):
         """
         Calculate the parallel power at the midplane.
+
+        Returns
+        -------
+        :
+            parallel power at the midplane
         """
         p_sol_near = self.params.P_sep_particle * self.params.f_p_sol_near
         p_sol_far = self.params.P_sep_particle * (1 - self.params.f_p_sol_near)
@@ -432,6 +465,11 @@ def _q_par(self, x, dx, B, Bp, *, outboard=True):
     def plot(self, ax: Axes = None, *, show=False) -> Axes:
         """
         Plot the ChargedParticleSolver results.
+
+        Returns
+        -------
+        :
+            The axes object on which the ChargedParticleSolver is plotted.
         """
         if ax is None:
             _, ax = plt.subplots()
@@ -472,6 +510,11 @@ def _make_params(config):
             Unsupported config type
         ValueError
             Unknown configuration parameters
+
+        Returns
+        -------
+        :
+            a ChargedParticleSolverParams object
         """
         if isinstance(config, dict):
             try:
diff --git a/bluemira/radiation_transport/flux_surfaces_maker.py b/bluemira/radiation_transport/flux_surfaces_maker.py
index 299051ec37..215644e1c3 100644
--- a/bluemira/radiation_transport/flux_surfaces_maker.py
+++ b/bluemira/radiation_transport/flux_surfaces_maker.py
@@ -11,13 +11,14 @@
 from copy import deepcopy
 
 import numpy as np
+from numpy import typing as npt
 
 from bluemira.base.constants import EPS
 from bluemira.base.look_and_feel import bluemira_warn
 from bluemira.equilibria.equilibrium import Equilibrium
 from bluemira.equilibria.find import find_flux_surface_through_point
 from bluemira.equilibria.find_legs import LegFlux, NumNull, SortSplit
-from bluemira.equilibria.flux_surfaces import OpenFluxSurface
+from bluemira.equilibria.flux_surfaces import OpenFluxSurface, PartialOpenFluxSurface
 from bluemira.geometry.coordinates import Coordinates, coords_plane_intersect
 from bluemira.geometry.plane import BluemiraPlane
 from bluemira.geometry.wire import BluemiraWire
@@ -28,18 +29,38 @@
 
 def analyse_first_wall_flux_surfaces(
     equilibrium: Equilibrium, first_wall: BluemiraWire, dx_mp: float
-):
+) -> tuple[
+    npt.NDArray[float],
+    npt.NDArray[float] | None,
+    list[PartialOpenFluxSurface],
+    float,
+    float | None,
+]:
     """
     A simplified charged particle transport model along open field lines.
 
     Parameters
     ----------
-    equilibrium:
+    equilibrium: Equilibrium
         The equilibrium defining flux surfaces.
     first_wall:
-    dx_mp:
+    dx_mp: float
         The midplane spatial resolution between flux surfaces [m]
         (default: 0.001).
+
+
+    Returns
+    -------
+    dx_omp: np.ndarray
+        The midplane spatial resolution between flux surfaces at the outboard [m]
+    dx_imp: np.ndarray or None
+        The midplane spatial resolution between flux surfaces at the inboard [m]
+    flux_surfaces: list[PartialOpenFluxSurface]
+        list of flux surfaces, all of which terminating at the first walls.
+    x_sep_omp: float
+        intersection between the separatrix outboard and mid-plane.
+    x_sep_imp: float or None
+        intersection between the separatrix inboard and mid-plane.
     """
     o_point = equilibrium.get_OX_points()[0][0]  # 1st o_point
     z = o_point.z
@@ -63,6 +84,11 @@ def analyse_first_wall_flux_surfaces(
 def _process_first_wall(first_wall: Coordinates) -> Coordinates:
     """
     Force working first wall geometry to be closed and counter-clockwise.
+
+    Returns
+    -------
+    first_wall: Coordinates
+        A closed first wall geometry, running counter clockwise.
     """
     first_wall = deepcopy(first_wall)
 
@@ -76,9 +102,20 @@ def _process_first_wall(first_wall: Coordinates) -> Coordinates:
     return first_wall
 
 
-def _analyse_SN(first_wall, dx_mp, equilibrium, o_point, yz_plane):
+def _analyse_SN(
+    first_wall, dx_mp, equilibrium, o_point, yz_plane
+) -> tuple[npt.NDArray[float], list[PartialOpenFluxSurface], float]:
     """
     Calculation for the case of single nulls.
+
+    Returns
+    -------
+    :
+        horizontal distances between outboard flux surfaces and outboard separatrix.
+    :
+        list of flux surfaces, all of which terminating at the first walls.
+    x_sep_omp: float
+        intersection between the separatrix outboard and mid-plane.
     """
     x_sep_omp, x_out_omp = _get_sep_out_intersection(
         equilibrium, first_wall, yz_plane, outboard=True
@@ -97,9 +134,26 @@ def _analyse_SN(first_wall, dx_mp, equilibrium, o_point, yz_plane):
     )
 
 
-def _analyse_DN(first_wall, dx_mp, equilibrium: Equilibrium, o_point, yz_plane):
+def _analyse_DN(
+    first_wall: Coordinates, dx_mp, equilibrium: Equilibrium, o_point, yz_plane
+) -> tuple[
+    npt.NDArray[float], npt.NDArray[float], list[PartialOpenFluxSurface], float, float
+]:
     """
     Calculation for the case of double nulls.
+
+    Returns
+    -------
+    :
+        horizontal distances between outboard flux surfaces and outboard separatrix.
+    :
+        horizontal distances between inboard flux surfaces and inboard separatrix.
+    :
+        list of flux surfaces, all of which terminating at the first walls.
+    x_sep_omp: float
+        intersection between the separatrix outboard and mid-plane.
+    x_sep_imp: float
+        intersection between the separatrix inboard and mid-plane.
     """
     x_sep_omp, x_out_omp = _get_sep_out_intersection(
         equilibrium, first_wall, yz_plane, outboard=True
@@ -126,10 +180,18 @@ def _analyse_DN(first_wall, dx_mp, equilibrium: Equilibrium, o_point, yz_plane):
     )
 
 
-def _clip_flux_surfaces(first_wall, flux_surfaces):
+def _clip_flux_surfaces(
+    first_wall: Coordinates, flux_surfaces: list[PartialOpenFluxSurface]
+) -> list[PartialOpenFluxSurface]:
     """
     Clip the flux surfaces to a first wall. Catch the cases where no intersections
     are found.
+
+    Returns
+    -------
+    flux_surfaces: list[PartialOpenFluxSurface]
+        A list of flux surface groups. Each group only contains flux surfaces that
+        intersect the first_wall.
     """
     for group in flux_surfaces:
         if group:
@@ -139,52 +201,91 @@ def _clip_flux_surfaces(first_wall, flux_surfaces):
                     # No intersection detected between flux surface and first wall
                     # Drop the flux surface from the group
                     group.pop(i)  # noqa: B909
+                    # TODO: fix this ^ B909 error
     return flux_surfaces
 
 
-def get_array_x_mp(flux_surfaces):
+def get_array_x_mp(flux_surfaces) -> npt.NDArray[float]:
     """
-    Get x_mp array of flux surface values.
+    Get the x-coordinate of the mid-plane intersection point for each flux surface.
+
+    Returns
+    -------
+    :
+        np.ndarray of mid-plane intersection point x-coordinate.
+
     """
     return np.array([fs.x_start for fs in flux_surfaces])
 
 
-def get_array_z_mp(flux_surfaces):
+def get_array_z_mp(flux_surfaces) -> npt.NDArray[float]:
     """
-    Get z_mp array of flux surface values.
+    Get the z-coordinate of the mid-plane intersection point for each flux surface.
+
+    Returns
+    -------
+    :
+        np.ndarray of mid-plane intersection point z-coordinate.
     """
     return np.array([fs.z_start for fs in flux_surfaces])
 
 
-def get_array_x_fw(flux_surfaces):
+def get_array_x_fw(flux_surfaces) -> npt.NDArray[float]:
     """
-    Get x_fw array of flux surface values.
+    Get the x-coordinate of the first-wall intersection point for each flux surface.
+
+    Returns
+    -------
+    :
+        np.ndarray of first-wall intersection point x-coordinate.
     """
     return np.array([fs.x_end for fs in flux_surfaces])
 
 
-def get_array_z_fw(flux_surfaces):
+def get_array_z_fw(flux_surfaces) -> npt.NDArray[float]:
     """
-    Get z_fw array of flux surface values.
+    Get the z-coordinate of the first-wall intersection point for each flux surface.
+
+    Returns
+    -------
+    :
+        np.ndarray of first-wall intersection point z-coordinate.
     """
     return np.array([fs.z_end for fs in flux_surfaces])
 
 
-def get_array_alpha(flux_surfaces):
+def get_array_alpha(flux_surfaces) -> npt.NDArray[float]:
     """
-    Get alpha angle array of flux surface values.
+    Get the alpha angle for each flux surface.
+
+    Returns
+    -------
+    :
+        np.ndarray of alpha.
     """
     return np.array([fs.alpha for fs in flux_surfaces])
 
 
-def _get_sep_out_intersection(eq: Equilibrium, first_wall, yz_plane, *, outboard=True):
+def _get_sep_out_intersection(
+    eq: Equilibrium, first_wall, yz_plane, *, outboard=True
+) -> tuple[float, float]:
     """
-    Find the middle and maximum outboard mid-plane psi norm values
+    Find the x-coordinate of where the mid-plane intersect the separatrix and at
+    the inboard/outboard.
 
     Raises
     ------
     RadiationTransportError
         Separatrix doesnt cross midplane
+
+    Returns
+    -------
+    x_sep_mp: float
+        the x-coordinate of the intersection point between the inboard-side separatrix
+        (outboard=False)/outboard-side separatrix (inboard=True) and the mid-plane.
+    x_out_mp: float
+        the x-coordinate of the intersection point between the inboard first wall
+        (outboard=False)/outboard first wall (outboard=True), and the mid-plane.
     """
     sep = LegFlux(eq)
 
@@ -220,7 +321,12 @@ def _get_sep_out_intersection(eq: Equilibrium, first_wall, yz_plane, *, outboard
 
 def _make_flux_surfaces(x, z, equilibrium, o_point, yz_plane):
     """
-    Make individual PartialOpenFluxSurfaces through a point.
+    Make individual PartialOpenFluxSurface through a point.
+
+    Returns
+    -------
+    :
+        The PartialOpenFluxSurface that passes through the point.
     """
     coords = find_flux_surface_through_point(
         equilibrium.x, equilibrium.z, equilibrium.psi(), x, z, equilibrium.psi(x, z)
@@ -232,9 +338,16 @@ def _make_flux_surfaces(x, z, equilibrium, o_point, yz_plane):
 
 def _make_flux_surfaces_ibob(
     dx_mp, equilibrium, o_point, yz_plane, x_sep_mp, x_out_mp, *, outboard: bool
-):
+) -> tuple[list]:
     """
     Make the flux surfaces on the inboard or outboard.
+
+    Returns
+    -------
+    flux_surfaces_lfs:
+        inboard flux surfaces
+    flux_surfaces_hfs:
+        outboard flux surfaces
     """
     sign = 1 if outboard else -1
 
diff --git a/bluemira/radiation_transport/midplane_temperature_density.py b/bluemira/radiation_transport/midplane_temperature_density.py
index 4d8eb4b0ba..2facfc8685 100644
--- a/bluemira/radiation_transport/midplane_temperature_density.py
+++ b/bluemira/radiation_transport/midplane_temperature_density.py
@@ -69,12 +69,18 @@ class MidplaneProfiles:
     """electron temperature at the mid-plane. Unit [keV]"""
 
 
-def midplane_profiles(params: ParameterFrameLike):
+def midplane_profiles(params: ParameterFrameLike) -> MidplaneProfiles:
     """
     Calculate the core radiation source profiles.
 
     Temperature and density are assumed to be constant along a
     single flux tube.
+
+    Returns
+    -------
+    :
+        Electron density and electron temperature as function of rho, from the magnetic
+        axis to the separatrix along the midplane.
     """
     params = make_parameter_frame(params, MidplaneProfilesParams, allow_unknown=True)
     rho_ped = (params.rho_ped_n.value + params.rho_ped_t.value) / 2
@@ -112,6 +118,11 @@ def collect_rho_core_values(
     n_points_mantle:
         no of discretisation points to separatrix
 
+    Returns
+    -------
+    :
+        radial coordinate rho [dimensionless].
+
     Notes
     -----
     The plasma bulk is divided into plasma core and plasma mantle according to rho
@@ -145,6 +156,10 @@ def core_electron_density_temperature_profile(
     rho_ped:
         dimensionless pedestal radius. Values between 0 and 1
 
+    Returns
+    -------
+    :
+        electron density and electron temperature, stored as a MidplaneProfiles object.
 
     Notes
     -----
diff --git a/bluemira/radiation_transport/radiation_profile.py b/bluemira/radiation_transport/radiation_profile.py
index 5d4c3b0036..2e98c2ba57 100644
--- a/bluemira/radiation_transport/radiation_profile.py
+++ b/bluemira/radiation_transport/radiation_profile.py
@@ -325,6 +325,11 @@ def mp_profile_plot(
         imp_name:
             impurity names
 
+        Returns
+        -------
+        ax: plt.Axes
+            axes on which the mid-plane radiation power distribution profile is plotted.
+
         Notes
         -----
         if rad_power.ndim > 1 the plot shows as many line as the number
@@ -368,6 +373,10 @@ def __init__(
         impurity_content,
         impurity_data,
     ):
+        """
+        In addition to `Radiation`, this class also includes the impurity data of all
+        gases except Argon.
+        """
         super().__init__(eq, params)
 
         self.H_content = impurity_content["H"]
@@ -457,7 +466,7 @@ def calculate_core_distribution(self) -> list[list[np.ndarray]]:
 
         return self.rad
 
-    def calculate_core_radiation_map(self):
+    def calculate_core_radiation_map(self) -> None:
         """
         Build core radiation map.
         """
@@ -470,7 +479,7 @@ def calculate_core_radiation_map(self):
 
     def radiation_distribution_plot(
         self, flux_tubes: np.ndarray, power_density: np.ndarray, ax=None
-    ):
+    ) -> plt.Axes:
         """
         2D plot of the core radiation power distribution.
 
@@ -481,6 +490,10 @@ def radiation_distribution_plot(
         power_density:
             arrays containing the power radiation density of the
             points lying on each flux tube [MW/m^3]
+
+        Returns
+        -------
+        ax: plt.Axes
         """
         if ax is None:
             fig, ax = plt.subplots()
@@ -1173,6 +1186,11 @@ def radiation_distribution_plot(self, flux_tubes, power_density, firstwall, ax=N
             expected len(flux_tubes) = len(power_density)
         firstwall: Grid
             first wall geometry
+
+        Returns
+        -------
+        ax: plt.Axes
+            Axes on which the 2D radiation power distribution is plotted.
         """
         if ax is None:
             fig, ax = plt.subplots()
@@ -1228,6 +1246,11 @@ def poloidal_distribution_plot(
         flux_property
             arrays containing the property values associated
             to the points of each flux tube.
+
+        Returns
+        -------
+        ax: plt.Axes
+            Axes on which the flux tubes' properties are plotted.
         """
         if ax is None:
             _, ax = plt.subplots()
@@ -1262,6 +1285,13 @@ def plot_t_vs_n(flux_tube, t_distribution, n_distribution, ax1=None):
         n_distribution: np.array([np.array])
             arrays containing the density values associated
             to the points of the flux tube.
+
+        Returns
+        -------
+        ax1: plt.Axes
+            Axes object on which the electron temperature is plotted
+        ax2: plt.Axes
+            Axes object on which the electron density is plotted
         """
         if ax1 is None:
             _, ax1 = plt.subplots()
@@ -1774,6 +1804,7 @@ def calculate_sol_radiation_map(self, lfs: np.ndarray, hfs: np.ndarray):
         self.x_tot = np.concatenate([flux_tube.coords.x for flux_tube in flux_tubes])
         self.z_tot = np.concatenate([flux_tube.coords.z for flux_tube in flux_tubes])
         self.rad_tot = np.concatenate(
+            # concatenate these two lists into a single float np.ndarray
             functools.reduce(operator.iadd, [self.total_rad_lfs, self.total_rad_hfs], [])
         )
 
@@ -1865,6 +1896,12 @@ def analyse(
         firstwall_geom:
             The closed first wall geometry
 
+        Returns
+        -------
+        self.core_rad: CoreRadiation
+            Core radiation source
+        self.sol_rad: ScrapeOffLayerRadiation
+            Scrape-off-layer radiation source
         """
         self.core_rad = CoreRadiation(
             self.eq,
@@ -2031,6 +2068,15 @@ def rad_map(self, firstwall_geom: Grid):
         Mapping all the radiation values associated to all the points
         as three arrays containing x coordinates, z coordinates and
         local radiated power density [MW/m^3]
+
+        Returns
+        -------
+        self.x_tot:
+            x-coordinates of flux tubes [m]
+        self.z_tot:
+            z-coordinates of flux tubes [m]
+        self.rad_tot:
+            total radiated power density [MW/m^3]
         """
         self.core_rad.calculate_core_radiation_map()
 
@@ -2046,9 +2092,14 @@ def rad_map(self, firstwall_geom: Grid):
 
         return self.x_tot, self.z_tot, self.rad_tot
 
-    def plot(self, ax=None):
+    def plot(self, ax=None) -> plt.Axes:
         """
         Plot the RadiationSolver results.
+
+        Returns
+        -------
+        ax: plt.Axes
+            The axes object on which radiation solver results are plotted.
         """
         if ax is None:
             fig, ax = plt.subplots()
diff --git a/bluemira/radiation_transport/radiation_tools.py b/bluemira/radiation_transport/radiation_tools.py
index d1206c9020..b7c691b66c 100644
--- a/bluemira/radiation_transport/radiation_tools.py
+++ b/bluemira/radiation_transport/radiation_tools.py
@@ -595,7 +595,7 @@ def radiative_loss_function_values(
 
 def radiative_loss_function_plot(
     t_ref: np.ndarray, lz_val: Iterable[np.ndarray], species: Iterable[str]
-):
+) -> plt.Axes:
     """
     Radiative loss function plot for a set of given impurities.
 
@@ -606,7 +606,12 @@ def radiative_loss_function_plot(
     l_z:
         radiative power loss reference [Wm^3]
     species:
-        name species
+        species names
+
+    Returns
+    -------
+    ax: plt.Axes
+        The axes object on which radiative loss function is plotted.
     """
     _fig, ax = plt.subplots()
     plt.title("Radiative loss functions vs Electron Temperature")
@@ -808,9 +813,16 @@ def filtering_in_or_out(
 
 def get_impurity_data(
     impurities_list: Iterable[str] = ("H", "He"), confinement_time_ms: float = 0.1
-):
+) -> dict[str, dict[str, tuple[np.ndarray, np.ndarray, np.ndarray]]]:
     """
     Function getting the PROCESS impurity data
+
+    Returns
+    -------
+    impurity_data:
+        The dictionary of impurities at the defined time, sorted by species, then sorted
+        by "T_ref" v.s. "L_ref", where "T_ref" = reference ion temperature [eV], "L_ref"
+        = the loss function value $L_z(n_e, T_e)$ [W m^3].
     """
     # This is a function
     imp_data_getter = get_code_interface("PROCESS").Solver.get_species_data
@@ -846,9 +858,14 @@ def detect_radiation(
     world: World,
     *,
     verbose: bool = False,
-):
+) -> DetectedRadiation:
     """
     To sample the wall and detect radiation
+
+    Returns
+    -------
+    :
+        DetectedRadiation object describing the radiation data.
     """
     # Storage lists for results
     power_density = []
@@ -945,6 +962,11 @@ def make_wall_detectors(
 ) -> list[WallDetector]:
     """
     To make the detectors on the wall
+
+    Returns
+    -------
+    wall_detectors:
+        list of WallDetectors
     """
     # number of detectors
     num = np.shape(wall_r)[0] - 2
@@ -1044,7 +1066,7 @@ def plot_radiation_loads(
     radiation_function, wall_detectors, wall_loads, plot_title, fw_shape
 ):
     """
-    To plot the radiation on the wall as MW/m^2
+    To plot the radiation on the wall as [MW/m^2]
     """
     min_r = min(fw_shape.x)
     max_r = max(fw_shape.x)
@@ -1137,8 +1159,19 @@ def solve(
         *,
         plot: bool = True,
         verbose: bool = False,
-    ):
-        """Solve first wall radiation problem"""
+    ) -> DetectedRadiation:
+        """
+        Solve first wall radiation problem
+
+        Plots
+        -----
+        Plot and show the radiation on the wall [MW/m^2], if plot=True.
+
+        Returns
+        -------
+        wall_loads: DetectedRadiation
+            Detected radiation data.
+        """
         shift = translate(0, 0, np.min(self.fw_shape.z))
         height = np.max(self.fw_shape.z) - np.min(self.fw_shape.z)
         rad_3d = AxisymmetricMapper(self.rad_source)
diff --git a/pyproject.toml b/pyproject.toml
index bb5940905d..b0b7b2dce6 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -363,9 +363,7 @@ ignore-names = [
 # "examples/*" = ["DOC201"]  # 51
 "bluemira/utilities/*" = ["DOC201"]  # 38
 "bluemira/structural/*" = ["DOC201"]  # 29
-"bluemira/radiation_transport/*" = ["DOC201"]  # 103
-"bluemira/plasma_physics/*" = ["DOC201"]  # 4
-"bluemira/optimisation/*" = ["DOC201"]  # 24
+"bluemira/radiation_transport/neutronics/*" = ["DOC201"] # 64
 "bluemira/mesh/*" = ["DOC201"]  # 10
 "bluemira/materials/*" = ["DOC201"]  # 62
 "bluemira/magnetostatics/*" = ["DOC201"]  # 45