diff --git a/src/FOC_reduction.py b/src/FOC_reduction.py
index 864e7ab..d3e6ac6 100755
--- a/src/FOC_reduction.py
+++ b/src/FOC_reduction.py
@@ -192,19 +192,20 @@ def main():
     P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P = proj_red.compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers)
 
     ## Step 4:
-    # Save image to FITS.
-    Stokes_test = proj_fits.save_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P, headers, figname+figtype, data_folder=data_folder, return_hdul=True)
+    # crop to desired region of interest (roi)
+#    stokescrop = proj_plots.crop_map(copy.deepcopy(stokes_test), copy.deepcopy(data_mask), snrp_cut=snrp_cut, snri_cut=snri_cut)
+#    stokescrop.run()
+#    stokes_crop, data_mask = stokescrop.crop()
 
     ## Step 5:
-    # Crop to desired Region of Interest (ROI)
-#    StokesCrop = proj_plots.crop_map(copy.deepcopy(Stokes_test), copy.deepcopy(data_mask), SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut)
-#    StokesCrop.run()
-#    Stokes_crop, data_mask = StokesCrop.crop()
+    # Save image to FITS.
+    Stokes_test = proj_fits.save_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P, headers, figname+figtype, data_folder=data_folder, return_hdul=True)
 
     # Plot polarization map (Background is either total Flux, Polarization degree or Polarization degree error).
     proj_plots.polarization_map(copy.deepcopy(Stokes_test), data_mask, rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype, plots_folder=plots_folder, display=None)
     proj_plots.polarization_map(copy.deepcopy(Stokes_test), data_mask, rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_P_flux", plots_folder=plots_folder, display='Pol_Flux')
     proj_plots.polarization_map(copy.deepcopy(Stokes_test), data_mask, rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_P", plots_folder=plots_folder, display='Pol_deg')
+    proj_plots.polarization_map(copy.deepcopy(Stokes_test), data_mask, rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_I_err", plots_folder=plots_folder, display='I_err')
     proj_plots.polarization_map(copy.deepcopy(Stokes_test), data_mask, rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_P_err", plots_folder=plots_folder, display='Pol_deg_err')
     proj_plots.polarization_map(copy.deepcopy(Stokes_test), data_mask, rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_SNRi", plots_folder=plots_folder, display='SNRi')
     proj_plots.polarization_map(copy.deepcopy(Stokes_test), data_mask, rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_SNRp", plots_folder=plots_folder, display='SNRp')
diff --git a/src/lib/plots.py b/src/lib/plots.py
index 51c472a..ccb077d 100755
--- a/src/lib/plots.py
+++ b/src/lib/plots.py
@@ -372,8 +372,15 @@ def polarization_map(Stokes, data_mask, rectangle=None, SNRp_cut=3., SNRi_cut=30
     elif display.lower() in ['s_p','pol_err','pol_deg_err']:
         # Display polarization degree error map
         vmin, vmax = 0., 10.
-        im = ax.imshow(pol_err.data*100.,extent=[-pol_err.data.shape[1]/2.,pol_err.data.shape[1]/2.,-pol_err.data.shape[0]/2.,pol_err.data.shape[0]/2.], vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.)
+        p_err = pol_err.data.copy()
+        p_err[p_err*100. > vmax] = np.nan
+        im = ax.imshow(p_err*100.,extent=[-pol_err.data.shape[1]/2.,pol_err.data.shape[1]/2.,-pol_err.data.shape[0]/2.,pol_err.data.shape[0]/2.], vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.)
         cbar = plt.colorbar(im, cax=cbar_ax, label=r"$\sigma_P$ [%]")
+    elif display.lower() in ['s_i','i_err']:
+        # Display intensity error map
+        vmin, vmax = 0., np.max(np.sqrt(stk_cov.data[0,0][stk_cov.data[0,0] > 0.])*convert_flux)
+        im = ax.imshow(np.sqrt(stk_cov.data[0,0])*convert_flux,extent=[-stkI.data.shape[1]/2.,stkI.data.shape[1]/2.,-stkI.data.shape[0]/2.,stkI.data.shape[0]/2.], vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.)
+        cbar = plt.colorbar(im, cax=cbar_ax, label=r"$\sigma_I$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]")
     elif display.lower() in ['snr','snri']:
         # Display I_stokes signal-to-noise map
         vmin, vmax = 0., np.max(SNRi[SNRi > 0.])
diff --git a/src/lib/reduction.py b/src/lib/reduction.py
index 988b043..a349bb0 100755
--- a/src/lib/reduction.py
+++ b/src/lib/reduction.py
@@ -863,13 +863,13 @@ def smooth_data(data_array, error_array, data_mask, headers, FWHM=1.,
                 for c in range(image.shape[1]):
                     dist_rc = np.where(data_mask, fmax, np.sqrt((r-xx)**2+(c-yy)**2))
                     g_rc = np.exp(-0.5*(dist_rc/stdev)**2)/(2.*np.pi*stdev**2)
-                    smoothed[i][r,c] = (1.-data_mask[r,c])*np.sum(data_array*weight*g_rc)/np.sum(weight*g_rc)
-                    error[i][r,c] = np.sqrt(np.sum(weight*g_rc**2))/np.sum(weight*g_rc)
+                    smoothed[i][r,c] = (1.-data_mask[r,c])*np.sum(image*g_rc)
+                    error[i][r,c] = np.sqrt(np.sum(error_array[i]*g_rc**2))
 
             # Nan handling
-            error[i][np.isnan(smoothed)] = 0.
-            smoothed[i][np.isnan(smoothed)] = 0.
-            error[i][np.isnan(error)] = 0.
+            error[i][np.isnan(smoothed[i])] = 0.
+            smoothed[i][np.isnan(smoothed[i])] = 0.
+            error[i][np.isnan(error[i])] = 0.
 
     else:
         raise ValueError("{} is not a valid smoothing option".format(smoothing))