diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml
index f09bc24..ddc6467 100644
--- a/.github/workflows/python-package.yml
+++ b/.github/workflows/python-package.yml
@@ -56,7 +56,7 @@ jobs:
         # test with as many cores as possible for speed
         pytest legwork -n=auto --cov=./ --cov-report=xml
     - name: Get Coverage
-      uses: orgoro/coverage@v3.1
+      uses: orgoro/coverage@v3.2
       with:
         coverageFile: ./coverage.xml
         token: ${{ secrets.GITHUB_TOKEN }}
diff --git a/changelog.md b/changelog.md
index fc9aa92..9b906a6 100644
--- a/changelog.md
+++ b/changelog.md
@@ -178,4 +178,8 @@ having two when we could just set `f_dom=2 f_orb`
 *TW 22/06/24*
 - Update to Python 3.11 and various dependencies
     - Fix deprecations based on these updates
-- Add __repr__ and __len__ functions for the Source class and its subclasses
\ No newline at end of file
+- Add __repr__ and __len__ functions for the Source class and its subclasses
+
+## 0.5.1
+*TW 11/01/25*
+- [Issue [#116](https://github.com/TeamLEGWORK/LEGWORK/issues/116)] LEGWORK will now raise errors if invalid harmonics are supplied
\ No newline at end of file
diff --git a/legwork/_version.py b/legwork/_version.py
index 3d18726..dd9b22c 100644
--- a/legwork/_version.py
+++ b/legwork/_version.py
@@ -1 +1 @@
-__version__ = "0.5.0"
+__version__ = "0.5.1"
diff --git a/legwork/strain.py b/legwork/strain.py
index f602626..ca6d843 100644
--- a/legwork/strain.py
+++ b/legwork/strain.py
@@ -102,6 +102,10 @@ def h_0_n(m_c, f_orb, ecc, n, dist, position=None, polarisation=None, inclinatio
     arrayed_args, _ = utils.ensure_array(m_c, f_orb, ecc, n, dist)
     m_c, f_orb, ecc, n, dist = arrayed_args
 
+    # raise a value error if any n is less than 1
+    if np.any(n < 1):
+        raise ValueError("All harmonics must be greater than or equal to 1")
+
     # if one timestep then extend dimensions
     if f_orb.ndim != 2:
         f_orb = f_orb[:, np.newaxis]
@@ -190,6 +194,10 @@ def h_c_n(m_c, f_orb, ecc, n, dist, position=None, polarisation=None, inclinatio
     arrayed_args, _ = utils.ensure_array(m_c, f_orb, ecc, n, dist)
     m_c, f_orb, ecc, n, dist = arrayed_args
 
+    # raise a value error if any n is less than 1
+    if np.any(n < 1):
+        raise ValueError("All harmonics must be greater than or equal to 1")
+
     # if one timestep then extend dimensions
     if f_orb.ndim != 2:
         f_orb = f_orb[:, np.newaxis]
diff --git a/legwork/tests/test_strains.py b/legwork/tests/test_strains.py
index 75efa5a..0f54b0e 100644
--- a/legwork/tests/test_strains.py
+++ b/legwork/tests/test_strains.py
@@ -28,3 +28,18 @@ def test_strain_conversion(self):
         fn_dot = utils.fn_dot(m_c, f_orb, e, n)
 
         self.assertTrue(np.allclose(should_be_fn_dot, fn_dot))
+
+    def test_bad_harmonics(self):
+        """Make sure fn_dot fails when given harmonics < 1"""
+
+        n_values = 100000
+
+        m_c = np.random.uniform(0, 10, n_values) * u.Msun
+        dist = np.random.uniform(0, 30, n_values) * u.kpc
+        f_orb = 10**(np.random.uniform(-5, -1, n_values)) * u.Hz
+        e = np.random.uniform(0, 0.9, n_values)
+
+        with self.assertRaises(ValueError):
+            strain.h_0_n(m_c, f_orb, e, [0, 1, 2], dist)
+        with self.assertRaises(ValueError):
+            strain.h_c_n(m_c, f_orb, e, [0, 1, 2], dist)
diff --git a/legwork/tests/test_utils.py b/legwork/tests/test_utils.py
index 4b7d707..bbb6907 100644
--- a/legwork/tests/test_utils.py
+++ b/legwork/tests/test_utils.py
@@ -68,3 +68,9 @@ def integrand2(theta, phi, psi, inc):
              [0, np.pi / 2]])  # inc
 
         self.assertAlmostEqual(result, 0.12)
+
+    def test_bad_harmonics(self):
+        """Make sure fn_dot fails when given harmonics < 1"""
+
+        with self.assertRaises(ValueError):
+            utils.fn_dot(m_c=10*u.Msun, f_orb=10*u.mHz, e=0.1, n=0)
\ No newline at end of file
diff --git a/legwork/utils.py b/legwork/utils.py
index 24eba02..18b3460 100644
--- a/legwork/utils.py
+++ b/legwork/utils.py
@@ -253,6 +253,9 @@ def fn_dot(m_c, f_orb, e, n):
     fn_dot : `float/array`
         Rate of change of nth frequency
     """
+    if np.any(n < 1):
+        raise ValueError("All harmonics must be greater than or equal to 1")
+
     fn_dot = (48 * n) / (5 * np.pi) * (c.G * m_c)**(5/3) / c.c**5 * (2 * np.pi * f_orb)**(11/3) * peters_f(e)
 
     # simplify units if present