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 f035a60..ca40f0b 100644
--- a/petrofit/modeling/models.py
+++ b/petrofit/modeling/models.py
@@ -13,8 +13,7 @@
 from astropy.modeling import FittableModel, Parameter, custom_model, models
 
 __all__ = [
-    'get_default_sersic_bounds', 'make_grid', 'PSFConvolvedModel2D',
-    'Nuker2D', 'Moffat2D', 'EllipMoffat2D',
+    'get_default_sersic_bounds', 'make_grid', 'PSFConvolvedModel2D', 'Nuker2D',
     'sersic_enclosed', 'sersic_enclosed_inv', 'sersic_enclosed_model', 'PetroApprox',
     'petrosian_profile', 'petrosian_model', 'get_default_gen_sersic_bounds', 'GenSersic2D'
 ]
@@ -607,28 +606,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