diff --git a/openpmd_viewer/addons/pic/lpa_diagnostics.py b/openpmd_viewer/addons/pic/lpa_diagnostics.py
index 7f5d4d7..4873fdb 100644
--- a/openpmd_viewer/addons/pic/lpa_diagnostics.py
+++ b/openpmd_viewer/addons/pic/lpa_diagnostics.py
@@ -675,6 +675,56 @@ def get_laser_envelope( self, t=None, iteration=None, pol=None,
         # Return the result
         return( envelope, info )
 
+
+    def get_electromagnetic_energy( self, t=None, iteration=None ):
+        """
+        Compute the total electromagnetic energy inside the box, i.e.
+
+        .. math::
+
+            \\mathcal{E}_{tot} = \\int dV\,\\left[ \\frac{\\epsilon_0}{2}\\boldsymbol{E}^2 + \\frac{1}{2\mu_0}\\boldsymbol{B}^2 \\right]
+
+        For simulations of high-intensity lasers propagating in underdense plasmas,
+        this is approximately equal to the laser energy (since the plasma fields are usually low)
+
+        Parameters:
+        -----------
+        t : float (in seconds), optional
+            Time at which to obtain the data (if this does not correspond to
+            an existing iteration, the closest existing iteration will be used)
+            Either `t` or `iteration` should be given by the user.
+
+        iteration : int
+            The iteration at which to obtain the data
+            Either `t` or `iteration` should be given by the user.
+
+        Returns:
+        --------
+        A float with the total electromagnetic energy inside the box (in Joules)
+        """
+        # Check that the required data is available
+        if 'E' not in self.avail_fields or 'B' not in self.avail_fields:
+            raise ValueError('The fields E and B are required to compute the electromagnetic energy.')
+        # Check that the fields have either the thetaMode or 3dcartesian geometry
+        if (self.fields_metadata['E']['geometry'] not in ['thetaMode', '3dcartesian']) or \
+              (self.fields_metadata['B']['geometry'] not in ['thetaMode', '3dcartesian']):
+            raise ValueError('The electromagnetic energy can only be computed for 3D Cartesian and cylindrical simulations.')
+
+        # For RZ: use `theta=None` to get the full 3D field
+        Ex, info = self.get_field('E', 'x', t=t, iteration=iteration, theta=None )
+        Ey, info = self.get_field('E', 'y', t=t, iteration=iteration, theta=None )
+        Ez, info = self.get_field('E', 'z', t=t, iteration=iteration, theta=None )
+        Bx, info = self.get_field('B', 'x', t=t, iteration=iteration, theta=None )
+        By, info = self.get_field('B', 'y', t=t, iteration=iteration, theta=None )
+        Bz, info = self.get_field('B', 'z', t=t, iteration=iteration, theta=None )
+
+        # Compute the energy
+        energy_density = const.epsilon_0/2.*(Ex**2 + Ey**2 + Ez**2) + 1./(2*const.mu_0)*(Bx**2 + By**2 + Bz**2)
+        volume = info.dx*info.dy*info.dz # Cell volume
+        E = energy_density.sum() * volume
+        return E
+
+
     def get_main_frequency( self, t=None, iteration=None, pol=None, m='all',
                             method='max'):
         """
diff --git a/tests/test_laser_energy.py b/tests/test_laser_energy.py
new file mode 100644
index 0000000..79f8c50
--- /dev/null
+++ b/tests/test_laser_energy.py
@@ -0,0 +1,57 @@
+"""
+This test file is part of the openPMD-viewer.
+
+It checks that the function `get_electromagnetic_energy` works correctly
+
+Usage:
+This file is meant to be run from the root directory of openPMD-viewer,
+by any of the following commands
+$ python tests/test_laser_energy.py
+$ py.test
+$ python -m pytest tests
+
+Copyright 2024, openPMD-viewer contributors
+Authors: Remi Lehe
+License: 3-Clause-BSD-LBNL
+"""
+from openpmd_viewer.addons import LpaDiagnostics
+from scipy.constants import c
+import numpy as np
+
+# Download required datasets
+import os
+def download_if_absent( dataset_name ):
+    "Function that downloads and decompress a chosen dataset"
+    if os.path.exists( dataset_name ) is False:
+        import wget, tarfile
+        tar_name = "%s.tar.gz" %dataset_name
+        url = "https://github.com/openPMD/openPMD-example-datasets/raw/draft/%s" %tar_name
+        wget.download(url, tar_name)
+        with tarfile.open( tar_name ) as tar_file:
+            tar_file.extractall()
+        os.remove( tar_name )
+download_if_absent( 'example-thetaMode' )
+
+def test_laser_energy():
+    """
+    Check that the function `get_electromagnetic_energy` gives the same result
+    as the engineering formula for a Gaussian pulse
+
+    E = 2.7e-5*(a0*w0/lambd)**2*(tau[fs])
+    """
+    ts = LpaDiagnostics('./example-thetaMode/hdf5')
+    iteration = 300
+
+    # Evaluate the laser energy using the engineering formula
+    a0 = ts.get_a0(iteration=iteration, pol='y')
+    w0 = ts.get_laser_waist(iteration=iteration, pol='y')
+    tau_fs = ts.get_ctau(iteration=iteration, pol='y') / c * 1e15
+    omega = ts.get_main_frequency(iteration=iteration, pol='y')
+    lambd = 2*np.pi*c/omega
+    E_eng = 2.7e-5*(a0*w0/lambd)**2*tau_fs
+
+    # Evaluate the laser energy using the function
+    E_func = ts.get_electromagnetic_energy(iteration=iteration)
+
+    # Compare the two results
+    assert np.isclose(E_eng, E_func, rtol=0.04)