diff --git a/CHANGELOG.md b/CHANGELOG.md
index f179173..b41e4e5 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -10,6 +10,7 @@
 
 * [Issue 255](https://github.com/MassimoCimmino/pygfunction/issues/255) - Default to an `orientation` of `0.` when `tilt` is `0.` in `boreholes.Borehole` class.
 * [Issue 274](https://github.com/MassimoCimmino/pygfunction/issues/274) - Fixed scalar assignment from ndim-1 array. It is deprecated as of `numpy` version `1.25`. Only ndim-0 arrays can be treated as scalars.
+* [Issue 285](https://github.com/MassimoCimmino/pygfunction/issues/285) - Use `numpy.complex128` instead of `numpy.cfloat`. This is to comply with backward-incompatible changes introduced in `numpy` version `2.0`.
 
 ## Version 2.2.2 (2023-01-09)
 
diff --git a/pygfunction/pipes.py b/pygfunction/pipes.py
index 5a93d08..c3650ef 100644
--- a/pygfunction/pipes.py
+++ b/pygfunction/pipes.py
@@ -3056,7 +3056,7 @@ def multipole(pos, r_out, r_b, k_s, k_g, R_fp, T_b, q_p, J,
     # Multipoles (EQ. 38)
     # -------------------
     if J > 0:
-        P = np.zeros((n_p, J), dtype=np.cfloat)
+        P = np.zeros((n_p, J), dtype=np.complex128)
         coeff = -np.array([[(1 - (k+1)*beta_m)/(1 + (k+1)*beta_m)
                            for k in range(J)] for beta_m in beta_p])
         while eps_max > eps and it < it_max:
@@ -3161,7 +3161,7 @@ def _F_mk(q_p, P, n_p, J, r_b, r_out, z, pikg, sigma):
         Matrix F_mk from Claesson and Hellstrom (2011), EQ. 34.
 
     """
-    F = np.zeros((n_p, J), dtype=np.cfloat)
+    F = np.zeros((n_p, J), dtype=np.complex128)
     dz = np.add.outer(z, -z) + np.eye(n_p)
     zz = np.multiply.outer(z, np.conj(z))
     for k in range(J):