diff --git a/src/load_distribution/load_distribution.py b/src/load_distribution/load_distribution.py
index 4bfab61..13c4ce5 100644
--- a/src/load_distribution/load_distribution.py
+++ b/src/load_distribution/load_distribution.py
@@ -215,21 +215,26 @@ def singularities_to_polygon(los: list[Singularity], xy: bool = False) -> Polygo
     # required shape, even if that means the exact x value is omitted (because we are
     # keeping the value immediately to the left and immediately to the right instead).
     x_acc = sorted(list(set(x_acc)))
-
     x_ord_count = Counter([round(x, 6) for x in x_acc])
     to_filter = []
     for key, count in x_ord_count.items():
         if count == 3:
             to_filter.append(key)
 
+    filtered_indexes = []
     if to_filter:
         rounded_x_acc = [round(x, 6) for x in x_acc]
         for filter_val in to_filter:
             index = rounded_x_acc.index(filter_val)
-            x_acc.pop(index)
+            filtered_indexes.append(index + 1)
+
+    filtered_x = []
+    for idx, x_val in enumerate(x_acc):
+        if idx not in filtered_indexes:
+            filtered_x.append(x_val)
 
     # Now, for every x value, compute the corresponding y value
-    y_acc = [sum([sing(x) for sing in sorted_sings]) for x in x_acc[:-1]]
+    y_acc = [sum([sing(x) for sing in sorted_sings]) for x in filtered_x[:-1]]
     # Always ends on 0.0
     y_acc.append(0.0)
     if xy:
@@ -237,9 +242,9 @@ def singularities_to_polygon(los: list[Singularity], xy: bool = False) -> Polygo
     else:
         precision = n if n else 2
         xy_acc = zip(
-            [round(x, precision) for x in x_acc], [round(y, precision) for y in y_acc]
+            [round(x, precision) for x in filtered_x],
+            [round(y, precision) for y in y_acc],
         )
-
         return Polygon(xy_acc)
 
 
diff --git a/tests/test_data.py b/tests/test_data.py
new file mode 100644
index 0000000..92af0d3
--- /dev/null
+++ b/tests/test_data.py
@@ -0,0 +1,319 @@
+from load_distribution import Singularity
+import numpy as np
+
+# An output from papermodels that caused an error in the singularities_to_polygon function due to lots of list popping
+
+big_list_o_singularities = [
+    Singularity(
+        x0=6.628,
+        x1=7.859,
+        m=np.float64(0.0),
+        y0=np.float64(-111.0726290117895),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=0.001,
+        x1=0.413,
+        m=np.float64(0.0),
+        y0=np.float64(1684.3833962328304),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=3.172,
+        x1=6.876,
+        m=np.float64(0.0),
+        y0=np.float64(-110.42669320082763),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=3.172,
+        x1=6.876,
+        m=np.float64(0.0),
+        y0=np.float64(-110.42669320082763),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=0.001,
+        x1=3.504,
+        m=np.float64(0.0),
+        y0=np.float64(-110.44792139070424),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=0.001,
+        x1=3.504,
+        m=np.float64(0.0),
+        y0=np.float64(-110.44792139070424),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=3.172,
+        x1=6.876,
+        m=np.float64(0.0),
+        y0=np.float64(-110.42669320082763),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=3.172,
+        x1=6.876,
+        m=np.float64(0.0),
+        y0=np.float64(-110.42669320082763),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=0.001,
+        x1=3.504,
+        m=np.float64(0.0),
+        y0=np.float64(-110.44792139070424),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=0.001,
+        x1=3.504,
+        m=np.float64(0.0),
+        y0=np.float64(-110.44792139070424),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=0.001,
+        x1=3.504,
+        m=np.float64(0.0),
+        y0=np.float64(-110.44792139070424),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=6.628,
+        x1=7.859,
+        m=np.float64(0.0),
+        y0=np.float64(-111.0726290117895),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=7.57,
+        x1=7.859,
+        m=np.float64(0.0),
+        y0=np.float64(1274.5157420687394),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=6.628,
+        x1=7.859,
+        m=np.float64(0.0),
+        y0=np.float64(-111.0726290117895),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=0.001,
+        x1=0.413,
+        m=np.float64(0.0),
+        y0=np.float64(1684.3833962328304),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=3.172,
+        x1=6.876,
+        m=np.float64(0.0),
+        y0=np.float64(-110.42669320082763),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=3.172,
+        x1=6.876,
+        m=np.float64(0.0),
+        y0=np.float64(-110.42669320082763),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=0.001,
+        x1=3.504,
+        m=np.float64(0.0),
+        y0=np.float64(-110.44792139070424),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=0.001,
+        x1=3.504,
+        m=np.float64(0.0),
+        y0=np.float64(-110.44792139070424),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=3.172,
+        x1=6.876,
+        m=np.float64(0.0),
+        y0=np.float64(-110.42669320082763),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=3.172,
+        x1=6.876,
+        m=np.float64(0.0),
+        y0=np.float64(-110.42669320082763),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=0.001,
+        x1=3.504,
+        m=np.float64(0.0),
+        y0=np.float64(-110.44792139070424),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=0.001,
+        x1=3.504,
+        m=np.float64(0.0),
+        y0=np.float64(-110.44792139070424),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=0.001,
+        x1=3.504,
+        m=np.float64(0.0),
+        y0=np.float64(-110.44792139070424),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=6.628,
+        x1=7.859,
+        m=np.float64(0.0),
+        y0=np.float64(-111.0726290117895),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=7.57,
+        x1=7.859,
+        m=np.float64(0.0),
+        y0=np.float64(1274.5157420687394),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=6.628,
+        x1=7.859,
+        m=np.float64(0.0),
+        y0=np.float64(-111.0726290117895),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=0.001,
+        x1=0.413,
+        m=np.float64(0.0),
+        y0=np.float64(1684.3833962328304),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=3.172,
+        x1=6.876,
+        m=np.float64(0.0),
+        y0=np.float64(-110.42669320082763),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=3.172,
+        x1=6.876,
+        m=np.float64(0.0),
+        y0=np.float64(-110.42669320082763),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=0.001,
+        x1=3.504,
+        m=np.float64(0.0),
+        y0=np.float64(-110.44792139070424),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=0.001,
+        x1=3.504,
+        m=np.float64(0.0),
+        y0=np.float64(-110.44792139070424),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=3.172,
+        x1=6.876,
+        m=np.float64(0.0),
+        y0=np.float64(-110.42669320082763),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=3.172,
+        x1=6.876,
+        m=np.float64(0.0),
+        y0=np.float64(-110.42669320082763),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=0.001,
+        x1=3.504,
+        m=np.float64(0.0),
+        y0=np.float64(-110.44792139070424),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=0.001,
+        x1=3.504,
+        m=np.float64(0.0),
+        y0=np.float64(-110.44792139070424),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=0.001,
+        x1=3.504,
+        m=np.float64(0.0),
+        y0=np.float64(-110.44792139070424),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=6.628,
+        x1=7.859,
+        m=np.float64(0.0),
+        y0=np.float64(-111.0726290117895),
+        precision=6,
+        eps=1e-12,
+    ),
+    Singularity(
+        x0=7.57,
+        x1=7.859,
+        m=np.float64(0.0),
+        y0=np.float64(1274.5157420687394),
+        precision=6,
+        eps=1e-12,
+    ),
+]
diff --git a/tests/test_load_distribution.py b/tests/test_load_distribution.py
index 610ec39..95ee7cc 100644
--- a/tests/test_load_distribution.py
+++ b/tests/test_load_distribution.py
@@ -1,4 +1,5 @@
 import load_distribution as ld
+import test_data as data
 from shapely import Polygon, GeometryCollection
 
 import pytest
@@ -236,6 +237,14 @@ def test_singularities_to_polygon():
         == "POLYGON ((0 0, 0.034368 0, 0.051573 0.5, 0.051573 0.5, 11.705 0.5, 11.705 0, 0 0))"
     )
 
+    # This test fails due to popping multiple items from a list where the indexes have changed
+    # because of the popping
+
+    assert (
+        (ld.singularities_to_polygon(data.big_list_o_singularities).wkt)
+        == "POLYGON ((0 0, 0.001 0, 0.001 3396.431373, 0.413 3396.431373, 0.413 -1656.718815, 3.172 -1656.718815, 3.172 -2981.839131, 3.504 -2981.839131, 3.504 -1325.120316, 6.628 -1325.120316, 6.628 -1991.55609, 6.876 -1991.55609, 6.876 -666.435774, 7.57 -666.435774, 7.57 3157.111452, 7.859 3157.111452, 7.859 0, 0 0))"
+    )
+
 
 def test_overlap_region_to_singularity():
     assert ld.overlap_region_to_singularity(