diff --git a/tweakwcs/linearfit.py b/tweakwcs/linearfit.py
index 01e2ba5..f9ca986 100644
--- a/tweakwcs/linearfit.py
+++ b/tweakwcs/linearfit.py
@@ -10,16 +10,18 @@
 """
 import logging
 import numbers
-from packaging.version import Version
 
-import numpy as np
 import astropy
-from astropy.modeling.fitting import LevMarLSQFitter
+import numpy as np
+from packaging.version import Version
+
 if Version(astropy.__version__) >= Version('5.1'):
     from astropy.modeling.fitting import fitter_to_model_params
+    from astropy.modeling.fitting import LMLSQFitter
 else:
     from astropy.modeling.fitting import (_fitter_to_model_params as
                                           fitter_to_model_params)
+    from astropy.modeling.fitting import LevMarLSQFitter as LMLSQFitter
 
 from . linalg import inv
 from . import __version__  # noqa: F401
@@ -27,8 +29,10 @@
 __author__ = 'Mihai Cara, Warren Hack'
 
 __all__ = [
-    'iter_linear_fit', 'build_fit_matrix', 'SUPPORTED_FITGEOM_MODES',
-    '_LevMarLSQFitter2x2'
+    'iter_linear_fit',
+    'build_fit_matrix',
+    'SUPPORTED_FITGEOM_MODES',
+    '_LevMarLSQFitter2x2',
 ]
 
 # Supported fitgeom modes and corresponding minobj
@@ -846,7 +850,7 @@ def build_fit_matrix(rot, scale=1):
     return matrix
 
 
-class _LevMarLSQFitter2x2(LevMarLSQFitter):
+class _LevMarLSQFitter2x2(LMLSQFitter):
     """ Performs fits of 2D vector-models to 2D reference points. """
     def objective_function(self, fps, *args):
         model, weights, inputs, meas = args
diff --git a/tweakwcs/tests/test_linearfit.py b/tweakwcs/tests/test_linearfit.py
index fafdc89..be8924d 100644
--- a/tweakwcs/tests/test_linearfit.py
+++ b/tweakwcs/tests/test_linearfit.py
@@ -7,9 +7,10 @@
 from itertools import product
 import math
 import sys
-import pytest
-import numpy as np
+
 from astropy.modeling.models import Shift, Rotation2D
+import numpy as np
+import pytest
 from tweakwcs import linearfit, linalg