Add possibility for el. current density in bipole (#220)

CHANGELOG.rst

@@ -10,6 +10,14 @@ v2.3.x
 """"""
 
 
+v2.4.0 El. current density
+--------------------------
+
+**2024-10-25**
+
+The ``mrec`` keyword in ``empymod.bipole`` can now be set to ``j``, in which case the electric current density (A/m2) is returned.
+
+
 v2.3.3 pyproject.toml
 ---------------------
 

empymod/model.py

@@ -151,8 +151,13 @@ def bipole(src, rec, depth, res, freqtime, signal=None, aniso=None,
         #mpermH = #mpermV = #res. If mpermH is provided but not mpermV,
         isotropic behaviour is assumed.
 
-    msrc, mrec : bool, default: False
+    msrc, mrec : bool or string, default: False
         If True, source/receiver (msrc/mrec) is magnetic, else electric.
+        The receiver can also be set to `mrec='j'`, in which case the electric
+        current density is returned (with the approximation that the current
+        density is proportional to the electric field, J=σE). Only implemented
+        for isotropic resistivities and electric permittivities at receiver
+        level.
 
     srcpts, recpts : int, default: 1
         Number of integration points for bipole source/receiver:
@@ -209,7 +214,7 @@ def bipole(src, rec, depth, res, freqtime, signal=None, aniso=None,
             - If > 0, interpolation is used.
 
           - `diff_quad`: criteria when to swap to QUAD (only relevant if
-            pts_per_dec=-1) (default: 100)
+            pts_per_dec>0) (default: 100)
           - `a`: lower limit for QUAD (default: first interval from QWE)
           - `b`: upper limit for QUAD (default: last interval from QWE)
           - `limit`: limit for quad (default: maxint)
@@ -310,6 +315,7 @@ def bipole(src, rec, depth, res, freqtime, signal=None, aniso=None,
 
         - If rec is electric, returns E [V/m].
         - If rec is magnetic, returns H [A/m].
+        - If rec is `j`, returns J [A/m2].
 
         EMArray is a subclassed ndarray with `.pha` and `.amp` attributes
         (only relevant for frequency-domain data).
@@ -422,6 +428,13 @@ def bipole(src, rec, depth, res, freqtime, signal=None, aniso=None,
     src, nsrc, nsrcz, srcdipole = check_bipole(src, 'src')
     rec, nrec, nrecz, recdipole = check_bipole(rec, 'rec')
 
+    # Check if receiver is a `j`.
+    if mrec == 'j':
+        rec_j = True
+        mrec = False
+    else:
+        rec_j = False
+
     # === 3. EM-FIELD CALCULATION ============
 
     # Pre-allocate output EM array
@@ -487,6 +500,11 @@ def bipole(src, rec, depth, res, freqtime, signal=None, aniso=None,
                     # Get layer number in which rec resides
                     lrec, zrec = get_layer_nr(tirec, depth)
 
+                    # Check eta at receiver level (only isotropic implemented).
+                    if rec_j and verb > 0 and etaH[0, lrec] != etaV[0, lrec]:
+                        print("* WARNING :: `etaH != etaV` at receiver level, "
+                              "only `etaH` considered for e-current density.")
+
                     # Gather variables
                     finp = (off, angle, zsrc, zrec, lsrc, lrec, depth, freq,
                             etaH, etaV, zetaH, zetaV, xdirect, isfullspace, ht,
@@ -528,6 +546,10 @@ def bipole(src, rec, depth, res, freqtime, signal=None, aniso=None,
                 src_rec_w *= np.tile(rec_w, isrc)
             sEM *= src_rec_w
 
+            # Multiply with eta of the rec-layer if ecurrent.
+            if rec_j:
+                sEM *= etaH[:, lrec, None]
+
             # Add this src-rec signal
             if nrec == nrecz:
                 if nsrc == nsrcz:  # Case 1: Looped over each src and each rec
@@ -898,7 +920,8 @@ def loop(src, rec, depth, res, freqtime, signal=None, aniso=None, epermH=None,
 
         - True: Magnetic dipole receiver;
         - False: Electric dipole receiver;
-        - 'loop': Magnetic receiver consisting of an electric-wire loop.
+        - 'loop': Magnetic receiver consisting of an electric-wire loop. Only
+          implemented for isotropic magnetic permeability at loop levels.
 
     recpts : int, default: 1
         Number of integration points for bipole receiver:
@@ -1119,7 +1142,7 @@ def loop(src, rec, depth, res, freqtime, signal=None, aniso=None, epermH=None,
                 # Get layer number in which rec resides
                 lrec, zrec = get_layer_nr(tirec, depth)
 
-                # Check mu at receiver level.
+                # Check mu at receiver level (only isotropic implemented).
                 if rec_loop and verb > 0 and mpermH[lrec] != mpermV[lrec]:
                     print("* WARNING :: `mpermH != mpermV` at receiver level, "
                           "only `mpermH` considered for loop factor.")

tests/test_model.py

@@ -588,6 +588,28 @@ def test_shape(self):
         assert_allclose(out[:, 0], b)
         assert_allclose(out[:, 1], d)
 
+    def test_j(self, capsys):
+        # Compare a electric * sigma with j
+        inp = {
+            'src': [-25, 25, -25, 25, 100, 170.7107],
+            'rec': [8000, 200, 300, 0, 0],
+            'depth': [0, 250],
+            'res': [1e20, 0.3, 5],
+            'freqtime': 1,
+        }
+        efield = bipole(**inp)
+        ecurr = bipole(mrec='j', **inp)
+        assert_allclose(efield/5, ecurr)
+
+        # Check warning
+        inp['aniso'] = [1, 1, 2]
+        efield = bipole(**inp)
+        _, _ = capsys.readouterr()  # Empty it
+        ecurr = bipole(mrec='j', **inp)
+        out, _ = capsys.readouterr()
+        assert_allclose(efield/5, ecurr)
+        assert "* WARNING :: `etaH != etaV` at receiver level, " in out
+
 
 def test_dipole():
     # As this is a subset of bipole, just run two tests to ensure