diff --git a/machine_design/designs/computation.py b/machine_design/designs/computation.py
index 1b76a3a..425101d 100644
--- a/machine_design/designs/computation.py
+++ b/machine_design/designs/computation.py
@@ -13,6 +13,7 @@
 class ComputationBase(ABC):
     def __init__(self, geometry: GeometryBase) -> None:
         self.geometry = geometry
+        self.phases = "ABCDEFGHIJK"[: self.geometry.n_phases]
         self.setup_name = "Setup1"
         self.rotor_mesh = None
         self.set_oper_params()
@@ -60,6 +61,24 @@ def compute_phase_groups(self):
         Q = int(g.geom_params["SlotNumber"])
         return phase_groups(Q, g.PolePairs, g.n_phases, belt_offset=g.belt_offset)
 
+    def compute_solution_expressions(self):
+        return [
+            "Moving1.Position",
+            "Moving1.Torque",
+            *[f"FluxLinkage(Phase{p})" for p in self.phases],
+            *[f"InducedVoltage(Phase{p})" for p in self.phases],
+            *[f"InputCurrent(Phase{p})" for p in self.phases],
+            *[f"L(Phase{x},Phase{y})" for x in self.phases for y in self.phases],
+        ]
+
+    def extract_phase_results(self, solutions):
+        # SI units
+        flux_phases = np.stack([self.extract_expression(solutions, f"FluxLinkage(Phase{p})") for p in self.phases], axis=-1)
+        vind_phases = np.stack([self.extract_expression(solutions, f"InducedVoltage(Phase{p})") for p in self.phases], axis=-1)
+        current_phases = np.stack([self.extract_expression(solutions, f"InputCurrent(Phase{p})") for p in self.phases], axis=-1)
+        L_raw = [np.stack([self.extract_expression(solutions, f"L(Phase{x},Phase{y})") for y in self.phases], axis=-1) for x in self.phases]
+        return flux_phases, vind_phases, current_phases, L_raw
+
     def push_variables(self, m2d: Maxwell2d) -> None:
         for k, v in self.oper_params.items():
             m2d[k] = v
diff --git a/motors/synrm_3f_36s.py b/motors/synrm_3f_36s.py
index c054e41..2d68dac 100644
--- a/motors/synrm_3f_36s.py
+++ b/motors/synrm_3f_36s.py
@@ -190,14 +190,7 @@ def set_oper_params(self):
         self.set_initpos()
 
     def set_solution_expressions(self):
-        self.solution_expressions = [
-            "Moving1.Position",
-            "Moving1.Torque",
-            *[f"FluxLinkage(Phase{p})" for p in "ABC"],
-            *[f"InducedVoltage(Phase{p})" for p in "ABC"],
-            *[f"InputCurrent(Phase{p})" for p in "ABC"],
-            *[f"L(Phase{x},Phase{y})" for x in "ABC" for y in "ABC"],
-        ]
+        self.solution_expressions = self.compute_solution_expressions()
 
     def set_output_vars(self):
         self.output_vars = {}
@@ -259,20 +252,17 @@ def extract_results(self, solutions):
         torque = self.extract_expression(solutions, "Moving1.Torque")
 
         theta_el = (position - np.deg2rad(self.InitPos)) * self.geometry.PolePairs
-        flux_phases = np.stack([self.extract_expression(solutions, f"FluxLinkage(Phase{p})") for p in "ABC"], axis=-1)
-        vind_phases = np.stack([self.extract_expression(solutions, f"InducedVoltage(Phase{p})") for p in "ABC"], axis=-1)
-        current_phases = np.stack([self.extract_expression(solutions, f"InputCurrent(Phase{p})") for p in "ABC"], axis=-1)
+
+        flux_phases, vind_phases, current_phases, L_raw = self.extract_phase_results(solutions)
 
         Flux_d, Flux_q = to_dq(flux_phases, theta_el, harmonic=1)
         Ui_d, Ui_q = to_dq(vind_phases, theta_el, harmonic=1)
         I_d, I_q = to_dq(current_phases, theta_el, harmonic=1)
 
-        L_raw = [np.stack([self.extract_expression(solutions, f"L(Phase{x},Phase{y})") for y in "ABC"], axis=-1) for x in "ABC"]
-
-        L_d_row = np.zeros((len(position), 3))
-        L_q_row = np.zeros((len(position), 3))
+        L_d_row = np.zeros((len(position), self.geometry.n_phases))
+        L_q_row = np.zeros((len(position), self.geometry.n_phases))
         for i, L_row in enumerate(L_raw):
-            L_d_row[:, i], L_q_row[:, i] = (v * m for v, m in zip(to_dq(L_row, theta_el, harmonic=1), (3 / 2, -3 / 2)))
+            L_d_row[:, i], L_q_row[:, i] = (v * m for v, m in zip(to_dq(L_row, theta_el, harmonic=1), (self.geometry.n_phases / 2, -self.geometry.n_phases / 2)))
 
         L_d, _ = to_dq(L_d_row, theta_el, harmonic=1)
         _, L_q_raw = to_dq(L_q_row, theta_el, harmonic=1)
diff --git a/motors/synrm_5f_40s.py b/motors/synrm_5f_40s.py
index 99256a7..4d38e62 100644
--- a/motors/synrm_5f_40s.py
+++ b/motors/synrm_5f_40s.py
@@ -60,14 +60,7 @@ def set_oper_params(self):
         self.set_initpos()
 
     def set_solution_expressions(self):
-        self.solution_expressions = [
-            "Moving1.Position",
-            "Moving1.Torque",
-            *[f"FluxLinkage(Phase{p})" for p in "ABCDE"],
-            *[f"InducedVoltage(Phase{p})" for p in "ABCDE"],
-            *[f"InputCurrent(Phase{p})" for p in "ABCDE"],
-            *[f"L(Phase{x},Phase{y})" for x in "ABCDE" for y in "ABCDE"],
-        ]
+        self.solution_expressions = self.compute_solution_expressions()
 
     def set_output_vars(self):
         self.output_vars = {}
@@ -231,9 +224,8 @@ def extract_results(self, solutions):
         torque = self.extract_expression(solutions, "Moving1.Torque")
 
         theta_el = electrical_angle(position, np.deg2rad(self.InitPos), self.geometry.PolePairs, self.RotSign, degrees=False)
-        flux_phases = np.stack([self.extract_expression(solutions, f"FluxLinkage(Phase{p})") for p in "ABCDE"], axis=-1)
-        vind_phases = np.stack([self.extract_expression(solutions, f"InducedVoltage(Phase{p})") for p in "ABCDE"], axis=-1)
-        current_phases = np.stack([self.extract_expression(solutions, f"InputCurrent(Phase{p})") for p in "ABCDE"], axis=-1)
+
+        flux_phases, vind_phases, current_phases, L_raw = self.extract_phase_results(solutions)
 
         Flux_d1, Flux_q1 = to_dq(flux_phases, theta_el, harmonic=1)
         Flux_d3, Flux_q3 = to_dq(flux_phases, theta_el, harmonic=3)
@@ -248,9 +240,9 @@ def extract_results(self, solutions):
         epsI3 = np.atan2(self.Iq3, self.Id3)
 
         time = np.array(solutions.primary_sweep_values)
-        dI_dt_phases = np.zeros((len(time), 5))
-        for k in range(5):
-            phase_offset = -2 * np.pi * k / 5  # 5 angles, 0deg, -72deg, -144deg, -216deg, -288deg
+        dI_dt_phases = np.zeros((len(time), self.geometry.n_phases))
+        for k in range(self.geometry.n_phases):
+            phase_offset = -2 * np.pi * k / self.geometry.n_phases
             dI_dt_phases[:, k] = -Im1 * self.w * np.sin(self.w * time + phase_offset + epsI1 - np.pi) - Im3 * 3 * self.w * np.sin(3 * (self.w * time + phase_offset) + epsI3 - np.pi)
 
         V_phases = vind_phases + self.Rstat * current_phases + self.Lew * dI_dt_phases
@@ -258,15 +250,13 @@ def extract_results(self, solutions):
         V_d1, V_q1 = to_dq(V_phases, theta_el, harmonic=1)
         V_d3, V_q3 = to_dq(V_phases, theta_el, harmonic=3)
 
-        L_raw = [np.stack([self.extract_expression(solutions, f"L(Phase{x},Phase{y})") for y in "ABCDE"], axis=-1) for x in "ABCDE"]
-
-        L_d1_row = np.zeros((len(time), 5))
-        L_q1_row = np.zeros((len(time), 5))
-        L_d3_row = np.zeros((len(time), 5))
-        L_q3_row = np.zeros((len(time), 5))
+        L_d1_row = np.zeros((len(time), self.geometry.n_phases))
+        L_q1_row = np.zeros((len(time), self.geometry.n_phases))
+        L_d3_row = np.zeros((len(time), self.geometry.n_phases))
+        L_q3_row = np.zeros((len(time), self.geometry.n_phases))
         for i, L_row in enumerate(L_raw):
-            L_d1_row[:, i], L_q1_row[:, i] = (v * 5 / 2 for v in to_dq(L_row, theta_el, harmonic=1))
-            L_d3_row[:, i], L_q3_row[:, i] = (v * 5 / 2 for v in to_dq(L_row, theta_el, harmonic=3))
+            L_d1_row[:, i], L_q1_row[:, i] = (v * self.geometry.n_phases / 2 for v in to_dq(L_row, theta_el, harmonic=1))
+            L_d3_row[:, i], L_q3_row[:, i] = (v * self.geometry.n_phases / 2 for v in to_dq(L_row, theta_el, harmonic=3))
 
         Ld1, Ld1q1 = to_dq(L_d1_row, theta_el, harmonic=1)
         Lq1d1, Lq1 = to_dq(L_q1_row, theta_el, harmonic=1)
@@ -282,7 +272,7 @@ def extract_results(self, solutions):
         Flux_e_d3 = Flux_d3 - (Ld3d1 * I_d1 + Ld3q1 * I_q1 + Ld3 * I_d3 + Ld3q3 * I_q3)
         Flux_e_q3 = Flux_q3 - (Lq3d1 * I_d1 + Lq3q1 * I_q1 + Lq3d3 * I_d3 + Lq3 * I_q3)
 
-        # Torque_dq = 5 / 2 * self.geometry.PolePairs * ((Flux_d1 * I_q1 - Flux_q1 * I_d1) + 3 * (Flux_d3 * I_q3 - Flux_q3 * I_d3))
+        # Torque_dq = self.geometry.n_phases / 2 * self.geometry.PolePairs * ((Flux_d1 * I_q1 - Flux_q1 * I_d1) + 3 * (Flux_d3 * I_q3 - Flux_q3 * I_d3))
 
         out = {
             "V_d1": V_d1,
diff --git a/tests/test_computation.py b/tests/test_computation.py
index f675f33..37c34b1 100644
--- a/tests/test_computation.py
+++ b/tests/test_computation.py
@@ -1,3 +1,5 @@
+from types import SimpleNamespace
+
 import pytest
 
 from machine_design.designs.computation import ComputationBase
@@ -35,19 +37,19 @@ def extract_results(self, solutions):
 
 
 def test_no_reserved_names_does_not_raise():
-    _DummyComputation(object(), {"Im": "1A"}, {"pos": "Moving1.Position"})
+    _DummyComputation(SimpleNamespace(n_phases=3), {"Im": "1A"}, {"pos": "Moving1.Position"})
 
 
 def test_reserved_name_in_oper_params_raises():
     with pytest.raises(ValueError, match="F"):
-        _DummyComputation(object(), {"F": "50Hz"}, {})
+        _DummyComputation(SimpleNamespace(n_phases=3), {"F": "50Hz"}, {})
 
 
 def test_reserved_name_in_output_vars_raises():
     with pytest.raises(ValueError, match="Time"):
-        _DummyComputation(object(), {}, {"Time": "0s"})
+        _DummyComputation(SimpleNamespace(n_phases=3), {}, {"Time": "0s"})
 
 
 def test_reserved_name_check_is_case_insensitive():
     with pytest.raises(ValueError, match="f"):
-        _DummyComputation(object(), {"f": "50Hz"}, {})
+        _DummyComputation(SimpleNamespace(n_phases=3), {"f": "50Hz"}, {})