diff --git a/docs/api_reference.rst b/docs/api_reference.rst
index c0efba1..8bc38d8 100644
--- a/docs/api_reference.rst
+++ b/docs/api_reference.rst
@@ -13,8 +13,6 @@ API
 .. automodapi:: petrofit.modeling.models
     :no-inheritance-diagram:
     :skip: Nuker2D
-    :skip: Moffat2D
-    :skip: EllipMoffat2D
     :skip: sersic_enclosed_model
     :skip: petrosian_model
 
diff --git a/petrofit/modeling/models.py b/petrofit/modeling/models.py
index 8891976..05604d7 100644
--- a/petrofit/modeling/models.py
+++ b/petrofit/modeling/models.py
@@ -547,10 +547,6 @@ class GenSersic2D(models.Sersic2D):
     c_0 : float
         Boxiness of elliptical isophote.
 
-    See Also
-    --------
-    Gaussian2D, Moffat2D
-
     Notes
     -----
     Model formula:
@@ -660,28 +656,6 @@ def CoreSersic2D(x, y, amplitude=1, r_eff=1, r_break=1, n=1, x_0=0, y_0=0,
     )
 
 
-@custom_model
-def Moffat2D(x, y, amplitude=1.0, x_0=0.0, y_0=0.0, gamma=1.0, alpha=1.0):
-    """Two dimensional Moffat function."""
-    rr_gg = ((x - x_0) ** 2 + (y - y_0) ** 2) / gamma ** 2
-    return amplitude * (1 + rr_gg) ** (-alpha)
-
-
-@custom_model
-def EllipMoffat2D(x, y, amplitude=1.0, x_0=0.0, y_0=0.0, gamma=1.0, alpha=1.0, ellip=0, theta=0, r=1):
-    """Two dimensional Moffat function."""
-
-    a, b = 1 * r, (1 - ellip) * r
-    cos_theta, sin_theta = np.cos(theta), np.sin(theta)
-    x_maj = (x - x_0) * cos_theta + (y - y_0) * sin_theta
-    x_min = -(x - x_0) * sin_theta + (y - y_0) * cos_theta
-    z = np.sqrt((x_maj / a) ** 2 + (x_min / b) ** 2)
-
-    rr_gg = (z) / gamma ** 2
-
-    return amplitude * (1 + rr_gg) ** (-alpha)
-
-
 def sersic_enclosed(
         r,
         amplitude,
diff --git a/petrofit/modeling/tests/test_fitting.py b/petrofit/modeling/tests/test_fitting.py
index add06b3..5c00af3 100644
--- a/petrofit/modeling/tests/test_fitting.py
+++ b/petrofit/modeling/tests/test_fitting.py
@@ -3,7 +3,7 @@
 from astropy.convolution import convolve
 from astropy.modeling import FittableModel, Parameter, custom_model, models
 
-from petrofit.modeling.models import Moffat2D, PSFModel, make_grid, PSFConvolvedModel2D
+from petrofit.modeling.models import PSFModel, make_grid, PSFConvolvedModel2D
 from petrofit.modeling import model_to_image, fit_model
 from petrofit.segmentation import masked_segm_image
 
@@ -39,7 +39,7 @@ def test_psfmodel():
 
     # Make a PSF
     x_grid, y_grid = make_grid(51, factor=1)
-    PSF = Moffat2D(x_0=25.0, y_0=25.0)(x_grid, y_grid)
+    PSF = models.Moffat2D(x_0=25.0, y_0=25.0)(x_grid, y_grid)
     PSF /= PSF.sum()
 
     # Make a PSF image using model image and PSF
@@ -110,7 +110,7 @@ def test_psf_convolved_image_model():
 
         # Make a PSF
         x_grid, y_grid = make_grid(51, factor=1)
-        PSF = Moffat2D(x_0=25.0, y_0=25.0)(x_grid, y_grid)
+        PSF = models.Moffat2D(x_0=25.0, y_0=25.0)(x_grid, y_grid)
         PSF /= PSF.sum()
 
         # Make a PSF image using model image and PSF