diff --git a/gprMax/cython/pml_build.pyx b/gprMax/cython/pml_build.pyx
new file mode 100644
index 00000000..421ca1a0
--- /dev/null
+++ b/gprMax/cython/pml_build.pyx
@@ -0,0 +1,70 @@
+# Copyright (C) 2015-2023: The University of Edinburgh, United Kingdom
+#                 Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
+#
+# This file is part of gprMax.
+#
+# gprMax is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# gprMax is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
+
+
+cpdef pml_average_er_mr(
+    str dir,
+    int s,
+    G
+):
+    """Calculates average permittivity and permeability for building PML 
+        (based on underlying material er and mr from solid array).
+
+    Args:
+        dir: string identifier for direction of PML.
+        s: int for starting cell of PML.
+        G: FDTDGrid class describing a grid in a model.
+    """
+
+    sumer = 0  # Sum of relative permittivities in PML slab
+    summr = 0  # Sum of relative permeabilities in PML slab
+
+    if dir == 'x':
+        for j in range(G.ny):
+            for k in range(G.nz):
+                numID = G.solid[s, j, k]
+                material = [x for x in G.materials if x.numID == numID]
+                material = material[0]
+                sumer += material.er
+                summr += material.mr
+        averageer = sumer / (G.ny * G.nz)
+        averagemr = summr / (G.ny * G.nz)
+
+    elif dir == 'y':
+        for i in range(G.nx):
+            for k in range(G.nz):
+                numID = G.solid[i, s, k]
+                material = [x for x in G.materials if x.numID == numID]
+                material = material[0]
+                sumer += material.er
+                summr += material.mr
+        averageer = sumer / (G.nx * G.nz)
+        averagemr = summr / (G.nx * G.nz)
+
+    elif dir == 'z':
+        for i in range(G.nx):
+            for j in range(G.ny):
+                numID = G.solid[i, j, s]
+                material = [x for x in G.materials if x.numID == numID]
+                material = material[0]
+                sumer += material.er
+                summr += material.mr
+        averageer = sumer / (G.nx * G.ny)
+        averagemr = summr / (G.nx * G.ny)
+
+    return averageer, averagemr
diff --git a/gprMax/pml.py b/gprMax/pml.py
index d92cb69e..a320c271 100644
--- a/gprMax/pml.py
+++ b/gprMax/pml.py
@@ -22,6 +22,8 @@ import numpy as np
 
 import gprMax.config as config
 
+from .cython.pml_build import pml_average_er_mr
+
 
 class CFSParameter:
     """Individual CFS parameter (e.g. alpha, kappa, or sigma)."""
@@ -609,9 +611,6 @@ def build_pml(G, key, value):
     elif config.sim_config.general['solver'] == 'opencl':
         pml_type = OpenCLPML
 
-    sumer = 0  # Sum of relative permittivities in PML slab
-    summr = 0  # Sum of relative permeabilities in PML slab
-
     if key[0] == 'x':
         if key == 'x0':
             pml = pml_type(G, ID=key, direction='xminus', 
@@ -621,14 +620,7 @@ def build_pml(G, key, value):
                            xs=G.nx - value, xf=G.nx, yf=G.ny, zf=G.nz)
         pml.CFS = G.pmls['cfs']
         G.pmls['slabs'].append(pml)
-        for j in range(G.ny):
-            for k in range(G.nz):
-                numID = G.solid[pml.xs, j, k]
-                material = next(x for x in G.materials if x.numID == numID)
-                sumer += material.er
-                summr += material.mr
-        averageer = sumer / (G.ny * G.nz)
-        averagemr = summr / (G.ny * G.nz)
+        averageer, averagemr = pml_average_er_mr('x', pml.xs, G)
 
     elif key[0] == 'y':
         if key == 'y0':
@@ -639,14 +631,7 @@ def build_pml(G, key, value):
                            ys=G.ny - value, xf=G.nx, yf=G.ny, zf=G.nz)
         pml.CFS = G.pmls['cfs']
         G.pmls['slabs'].append(pml)
-        for i in range(G.nx):
-            for k in range(G.nz):
-                numID = G.solid[i, pml.ys, k]
-                material = next(x for x in G.materials if x.numID == numID)
-                sumer += material.er
-                summr += material.mr
-        averageer = sumer / (G.nx * G.nz)
-        averagemr = summr / (G.nx * G.nz)
+        averageer, averagemr = pml_average_er_mr('y', pml.ys, G)
 
     elif key[0] == 'z':
         if key == 'z0':
@@ -657,13 +642,6 @@ def build_pml(G, key, value):
                            zs=G.nz - value, xf=G.nx, yf=G.ny, zf=G.nz)
         pml.CFS = G.pmls['cfs']
         G.pmls['slabs'].append(pml)
-        for i in range(G.nx):
-            for j in range(G.ny):
-                numID = G.solid[i, j, pml.zs]
-                material = next(x for x in G.materials if x.numID == numID)
-                sumer += material.er
-                summr += material.mr
-        averageer = sumer / (G.nx * G.ny)
-        averagemr = summr / (G.nx * G.ny)
+        averageer, averagemr = pml_average_er_mr('z', pml.zs, G)
 
     pml.calculate_update_coeffs(averageer, averagemr)