diff --git a/meanas/fdtd/energy.py b/meanas/fdtd/energy.py
index 26fb036..ffbbf80 100644
--- a/meanas/fdtd/energy.py
+++ b/meanas/fdtd/energy.py
@@ -4,23 +4,31 @@ import numpy
 from .. import dx_lists_t, field_t, field_updater
 
 
-def poynting(e, h):
-    s = (numpy.roll(e[1], -1, axis=0) * h[2] - numpy.roll(e[2], -1, axis=0) * h[1],
-         numpy.roll(e[2], -1, axis=1) * h[0] - numpy.roll(e[0], -1, axis=1) * h[2],
-         numpy.roll(e[0], -1, axis=2) * h[1] - numpy.roll(e[1], -1, axis=2) * h[0])
-    return numpy.array(s)
-
-
-def poynting_divergence(s=None, *, e=None, h=None, dxes=None): # TODO dxes
+def poynting(e, h, dxes=None):
     if dxes is None:
         dxes = tuple(tuple(numpy.ones(1) for _ in range(3)) for _ in range(2))
 
-    if s is None:
-        s = poynting(e, h)
+    ex = e[0] * dxes[0][0][:, None, None]
+    ey = e[1] * dxes[0][1][None, :, None]
+    ez = e[2] * dxes[0][2][None, None, :]
+    hx = h[0] * dxes[1][0][:, None, None]
+    hy = h[1] * dxes[1][1][None, :, None]
+    hz = h[2] * dxes[1][2][None, None, :]
 
-    ds = ((s[0] - numpy.roll(s[0], 1, axis=0)) / numpy.sqrt(dxes[0][0] * dxes[1][0])[:, None, None] +
-          (s[1] - numpy.roll(s[1], 1, axis=1)) / numpy.sqrt(dxes[0][1] * dxes[1][1])[None, :, None] +
-          (s[2] - numpy.roll(s[2], 1, axis=2)) / numpy.sqrt(dxes[0][2] * dxes[1][2])[None, None, :] )
+    s = numpy.empty_like(e)
+    s[0] = numpy.roll(ey, -1, axis=0) * hz - numpy.roll(ez, -1, axis=0) * hy
+    s[1] = numpy.roll(ez, -1, axis=1) * hx - numpy.roll(ex, -1, axis=1) * hz
+    s[2] = numpy.roll(ex, -1, axis=2) * hy - numpy.roll(ey, -1, axis=2) * hx
+    return s
+
+
+def poynting_divergence(s=None, *, e=None, h=None, dxes=None): # TODO dxes
+    if s is None:
+        s = poynting(e, h, dxes=dxes)
+
+    ds = ((s[0] - numpy.roll(s[0], 1, axis=0)) +
+          (s[1] - numpy.roll(s[1], 1, axis=1)) +
+          (s[2] - numpy.roll(s[2], 1, axis=2)) )
     return ds
 
 
diff --git a/meanas/test/test_fdtd.py b/meanas/test/test_fdtd.py
index 08b678e..1b185f5 100644
--- a/meanas/test/test_fdtd.py
+++ b/meanas/test/test_fdtd.py
@@ -77,16 +77,14 @@ def test_poynting_divergence(sim):
     args = {'dxes': sim.dxes,
             'epsilon': sim.epsilon}
 
-    dV = numpy.prod(numpy.meshgrid(*sim.dxes[0], indexing='ij'), axis=0)
-
     u_eprev = None
     for ii in range(1, 8):
         u_hstep = fdtd.energy_hstep(e0=sim.es[ii-1], h1=sim.hs[ii], e2=sim.es[ii],     **args)
         u_estep = fdtd.energy_estep(h0=sim.hs[ii],   e1=sim.es[ii], h2=sim.hs[ii + 1], **args)
-        delta_j_B = fdtd.delta_energy_j(j0=sim.js[ii], e1=sim.es[ii],   dxes=sim.dxes)
+        delta_j_B = fdtd.delta_energy_j(j0=sim.js[ii], e1=sim.es[ii], dxes=sim.dxes)
 
         du_half_h2e = u_estep - u_hstep - delta_j_B
-        div_s_h2e = sim.dt * fdtd.poynting_divergence(e=sim.es[ii], h=sim.hs[ii], dxes=sim.dxes) * dV
+        div_s_h2e = sim.dt * fdtd.poynting_divergence(e=sim.es[ii], h=sim.hs[ii], dxes=sim.dxes)
         assert_fields_close(du_half_h2e, -div_s_h2e)
 
         if u_eprev is None:
@@ -97,7 +95,7 @@ def test_poynting_divergence(sim):
         delta_j_A = fdtd.delta_energy_j(j0=sim.js[ii], e1=sim.es[ii-1], dxes=sim.dxes)
 
         du_half_e2h = u_hstep - u_eprev - delta_j_A
-        div_s_e2h = sim.dt * fdtd.poynting_divergence(e=sim.es[ii-1], h=sim.hs[ii], dxes=sim.dxes) * dV
+        div_s_e2h = sim.dt * fdtd.poynting_divergence(e=sim.es[ii-1], h=sim.hs[ii], dxes=sim.dxes)
         assert_fields_close(du_half_e2h, -div_s_e2h)
         u_eprev = u_estep
 
@@ -123,10 +121,7 @@ def test_poynting_planes(sim):
         u_hstep = fdtd.energy_hstep(e0=sim.es[ii-1], h1=sim.hs[ii], e2=sim.es[ii],     **args)
         u_estep = fdtd.energy_estep(h0=sim.hs[ii],   e1=sim.es[ii], h2=sim.hs[ii + 1], **args)
 
-        s_h2e = -fdtd.poynting(e=sim.es[ii], h=sim.hs[ii]) * sim.dt
-        s_h2e[0] *= sim.dxes[0][1][None, :, None] * sim.dxes[0][2][None, None, :]
-        s_h2e[1] *= sim.dxes[0][0][:, None, None] * sim.dxes[0][2][None, None, :]
-        s_h2e[2] *= sim.dxes[0][0][:, None, None] * sim.dxes[0][1][None, :, None]
+        s_h2e = -fdtd.poynting(e=sim.es[ii], h=sim.hs[ii], dxes=sim.dxes) * sim.dt
         planes = [s_h2e[px].sum(), -s_h2e[mx].sum(),
                   s_h2e[py].sum(), -s_h2e[my].sum(),
                   s_h2e[pz].sum(), -s_h2e[mz].sum()]
@@ -135,10 +130,7 @@ def test_poynting_planes(sim):
             u_eprev = u_estep
             continue
 
-        s_e2h = -fdtd.poynting(e=sim.es[ii - 1], h=sim.hs[ii]) * sim.dt
-        s_e2h[0] *= sim.dxes[0][1][None, :, None] * sim.dxes[0][2][None, None, :]
-        s_e2h[1] *= sim.dxes[0][0][:, None, None] * sim.dxes[0][2][None, None, :]
-        s_e2h[2] *= sim.dxes[0][0][:, None, None] * sim.dxes[0][1][None, :, None]
+        s_e2h = -fdtd.poynting(e=sim.es[ii - 1], h=sim.hs[ii], dxes=sim.dxes) * sim.dt
         planes = [s_e2h[px].sum(), -s_e2h[mx].sum(),
                   s_e2h[py].sum(), -s_e2h[my].sum(),
                   s_e2h[pz].sum(), -s_e2h[mz].sum()]