diff --git a/kwant/plotter.py b/kwant/plotter.py
index e160bddee28da7c8fd433376eaa80a7e26eb9ee9..03ca6260eda9a458c1092df24dd8615c4e673537 100644
--- a/kwant/plotter.py
+++ b/kwant/plotter.py
@@ -1872,11 +1872,7 @@ def interpolate_current(syst, current, relwidth=None, abswidth=None, n=9):
     else:
         width = abswidth
 
-    # TODO: Generalize 'factor' prefactor to arbitrary dimensions and remove
-    #       this check. This check is done here to keep changes local
-    if dim != 2:
-        raise ValueError("'interpolate_current' only works for 2D systems.")
-    factor = (3 / np.pi) / (width / 2)
+    # Define length scale in terms of the bump width.
     scale = 2 / width
     lens *= scale
 
@@ -1899,9 +1895,13 @@ def interpolate_current(syst, current, relwidth=None, abswidth=None, n=9):
     slices[:, :, 0] = np.floor((min_hops - bbox_min) * grid_density)
     slices[:, :, 1] = np.ceil((max_hops + 1.5*width - bbox_min) * grid_density)
 
-    # F is the antiderivative of the smoothing function
-    # f(ρ, z) = (1 - ρ^2 - z^2)^2 Θ(1 - ρ^2 - z^2), with respect to
-    # z,with F(ρ, -∞) = 0 and where Θ is the heaviside function.
+    # Define the smoothing function.
+
+    # We generate the smoothing function by convolving the current
+    # defined on a line between the two sites with
+    # f(ρ, z) = (1 - ρ^2 - z^2)^2 Θ(1 - ρ^2 - z^2), where ρ and z are
+    # cylindrical coords defined with respect to the hopping.
+    # 'F' is the result of the convolution.
     def F(rho, z):
         r = 1 - rho * rho
         r[r < 0] = 0
@@ -1912,6 +1912,20 @@ def interpolate_current(syst, current, relwidth=None, abswidth=None, n=9):
         mm = m * m
         return m * (mm * (mm/5 - (2/3) * rr) + rrrr) + (8 / 15) * rrrr * r
 
+    # We need to normalize the smoothing function so that it has unit cross
+    # section in the plane perpendicular to the hopping. This is equivalent
+    # to normalizing the integral of 'f' over the unit hypersphere to 1.
+    # The smoothing function goes as F(ρ) = (16/15) (1 - ρ^2)^(5/2) in the
+    # plane perpendicular to the hopping, so the cross section is:
+    # A_n = (16 / 15) * σ_n * ∫_0^1 ρ^(n-1) (1 - ρ^2)^(5/2) dρ
+    # where σ_n is the surface element prefactor (2 in 2D, 2π in 3D). Rather
+    # that calculate A_n every time, we hard code its value for 1, 2 and 3D.
+    if dim > 3:
+        raise ValueError("'interpolate_current' only works for systems "
+                         "with dimension <= 3.")
+    factors = [16 / 15, np.pi / 3, 32 * np.pi / 105]
+    cross_section_factor = factors[dim - 1]
+
     # Interpolate the field for each hopping.
     for i in range(len(current)):
 
@@ -1935,7 +1949,7 @@ def interpolate_current(syst, current, relwidth=None, abswidth=None, n=9):
         rho *= scale
 
         magns = F(rho, z) - F(rho, z - lens[i])
-        magns *= current[i] * factor
+        magns *= current[i] * scale / cross_section_factor
 
         field[field_slice] += dirs[i] * magns[..., None]