diff --git a/ccdproc/core.py b/ccdproc/core.py
index 0c321092..d5888663 100644
--- a/ccdproc/core.py
+++ b/ccdproc/core.py
@@ -438,7 +438,9 @@ def gain_correct(ccd, gain, gain_unit=None):
 
 @log_to_metadata
 def flat_correct(ccd, flat, min_value=None):
-    """Correct the image for flatfielding.
+    """Correct the image for flat fielding.
+
+    The flat field image is normalized by its mean before flat correcting.
 
     Parameters
     ----------
@@ -448,31 +450,31 @@ def flat_correct(ccd, flat, min_value=None):
     flat : `~ccdproc.ccddata.CCDData`
         Flatfield to apply to the data
 
-       min_value : None or float
-          Minimum value for flat field.  The value can either be None and no
-          minimum value is applied to the flat or specified by a float which
-          will replace all values in the flat by the min_value.
+    min_value : None or float
+        Minimum value for flat field.  The value can either be None and no
+        minimum value is applied to the flat or specified by a float which
+        will replace all values in the flat by the min_value.
 
     {log}
 
     Returns
     -------
     ccd :  `~ccdproc.ccddata.CCDData`
-      CCDData object with flat corrected
+        CCDData object with flat corrected
     """
     #Use the min_value to replace any values in the flat
+    use_flat = flat
     if min_value is not None:
-       flat.data[flat.data < min_value] = min_value
-
-    # normalize the flat
-    flat.data = flat.data / flat.data.mean()
-    if flat.uncertainty is not None:
-        flat.uncertainty.array = flat.uncertainty.array / flat.data.mean()
+        flat_min = flat.copy()
+        flat_min.data[flat_min.data < min_value] = min_value
+        use_flat = flat_min
 
     # divide through the flat
-    ccd.divide(flat)
-
-    return ccd
+    flat_corrected = ccd.divide(use_flat)
+    # multiply by the mean of the flat
+    flat_corrected = flat_corrected.multiply(use_flat.data.mean() *
+                                             use_flat.unit)
+    return flat_corrected
 
 
 def sigma_func(arr):
diff --git a/ccdproc/tests/test_ccdproc.py b/ccdproc/tests/test_ccdproc.py
index 1d7d986c..8886923f 100644
--- a/ccdproc/tests/test_ccdproc.py
+++ b/ccdproc/tests/test_ccdproc.py
@@ -321,28 +321,39 @@ def test_subtract_dark_fails(ccd_data):
 def test_flat_correct(ccd_data):
     size = ccd_data.shape[0]
     orig_mean = ccd_data.data.mean()
-    # create the flat
-    data = 2 * np.ones((size, size))
+    # create the flat, with some scatter
+    data = 2 * np.random.normal(loc=1.0, scale=0.05, size=(size, size))
     flat = CCDData(data, meta=fits.header.Header(), unit=ccd_data.unit)
     flat_data = flat_correct(ccd_data, flat)
 
     #check that the flat was normalized
-    assert flat_data.data.mean() == ccd_data.data.mean()
-    assert (ccd_data.data / flat_data.data == flat.data).all()
+    # Should be the case that flat * flat_data = ccd_data * flat.data.mean
+    # if the normalization was done correctly.
+    np.testing.assert_almost_equal((flat_data.data * flat.data).mean(),
+                                   ccd_data.data.mean() * flat.data.mean())
+    np.testing.assert_allclose(ccd_data.data / flat_data.data,
+                               flat.data / flat.data.mean())
+
 
 # test for flat correction with min_value
 @pytest.mark.data_scale(10)
 def test_flat_correct_min_value(ccd_data):
     size = ccd_data.shape[0]
-    orig_mean = ccd_data.data.mean()
     # create the flat
-    data = 2 * np.ones((size, size))
+    data = 2 * np.random.normal(loc=1.0, scale=0.05, size=(size, size))
     flat = CCDData(data, meta=fits.header.Header(), unit=ccd_data.unit)
-    flat_data = flat_correct(ccd_data, flat, min_value = 4)
-
+    flat_orig_data = flat.data.copy()
+    min_value = 2.1  # should replace some, but not all, values
+    flat_data = flat_correct(ccd_data, flat, min_value=min_value)
+    flat_with_min = flat.copy()
+    flat_with_min.data[flat_with_min.data < min_value] = min_value
     #check that the flat was normalized
-    assert flat_data.data.mean() == ccd_data.data.mean()
-    assert (ccd_data.data / flat_data.data == flat.data).all()
+    np.testing.assert_almost_equal((flat_data.data * flat_with_min.data).mean(),
+                                   ccd_data.data.mean() * flat_with_min.data.mean())
+    np.testing.assert_allclose(ccd_data.data / flat_data.data,
+                               flat_with_min.data / flat_with_min.data.mean())
+    # Test that flat is not modified.
+    assert (flat_orig_data == flat.data).all()
 
 
 # test for variance and for flat correction