IP and Q clarifications (#239)

Author: prisae

CHANGELOG.rst

@@ -10,11 +10,16 @@ v2.5.x
 """"""
 
 
-*latest*
---------
+v2.5.1 IP/Q clarifications
+--------------------------
+
+**2025-03-12**
+
+
+- Modelling routines ``ip_and_q``:
 
-- Modelling routines: Add warning to ``ip_and_q`` if ``src_z!=rec_z``; not
-  thoroughly tested.
+  - Use the corresponding co-planar loop configuration for the primary field.
+  - Add warning to it that it is experimental.
 
 
 v2.5.0 In-phase and quadrature

empymod/model.py

@@ -1775,6 +1775,8 @@ def ip_and_q(**kwargs):
     - 1st: with `xdirect=None` for the secondary field Hs.
     - 2nd: with `xdirect=True`, and only the first value of all model
       parameters and an empty `depth`-parameter for the primary field Hp.
+      For the primary field, the co-planar field of the source configuration is
+      used.
 
     It then returns the real (in-phase) and imaginary (quadrature) components
     of the ratio Hs/Hp, scaled by `scale`.
@@ -1786,9 +1788,17 @@ def ip_and_q(**kwargs):
     This function is only implemented for frequency-domain data of magnetic
     sources and receivers.
 
-    Warning: Only the case where source and receiver are at the same height has
-    been tested; if you encounter any issues with other configurations, please
-    let us know!
+    .. note::
+
+        This function is experimental. Please let us know if you encounter any
+        issues when using this function!
+
+        Also note that for the primary field computation, the co-planar field
+        of the source direction is used. Hence:
+
+        - Hp for ab=44 is used for ab in [44, 54, 64],
+        - Hp for ab=55 is used for ab in [45, 55, 65],
+        - Hp for ab=66 is used for ab in [46, 56, 66].
 
 
     Parameters
@@ -1826,8 +1836,8 @@ def ip_and_q(**kwargs):
     # Warning that it has not been thoroughly tested.
     verb = kwargs.get('verb', 2)
     if verb > 0 and not np.allclose(kwargs['src'][2], kwargs['rec'][2]):
-        print("* WARNING :: `src_z != rec_z`: has not been tested; "
-              "please report back any issues you might encounter.")
+        print("* WARNING :: This function is experimental. Please let us know "
+              "if you encounter any issues when using this function!")
 
     # Get or set scale
     scale = kwargs.pop('scale', 1e3)
@@ -1838,9 +1848,13 @@ def ip_and_q(**kwargs):
     # Primary magnetic field
     new = copy.deepcopy(kwargs)
 
-    # For PERP, ab=[46;64], Hp is zero; instead use Hp of the HCP config.
-    # Frischknecht et al., 1991, p. 111; doi: 10.1190/1.9781560802686.ch3.
-    new['ab'] = new['ab'] if new['ab'] not in [46, 64] else 66
+    # For the primary field, use the co-planar field of source configuration.
+    if new['ab'] in [44, 54, 64]:
+        new['ab'] = 44
+    elif new['ab'] in [45, 55, 65]:
+        new['ab'] = 55
+    elif new['ab'] in [46, 56, 66]:
+        new['ab'] = 66
 
     # Take only the first value of each parameter, and set depth to empty.
     new['depth'] = []

tests/test_model.py

@@ -1079,21 +1079,21 @@ def test_ip_and_q(capsys):
 
     # Status Quo
     system = {
-        'src': [0, 0, -1],
-        'rec': [2, 0, -1],
-        'freqtime': [1.0, 100.0, 10000.0],
-        'ab': 66,
-        'verb': 1,
+        "src": [0, 0, -1],
+        "rec": [2, 0, -1],
+        "freqtime": [1.0, 100.0, 10000.0],
+        "ab": 66,
+        "verb": 1,
     }
 
     model1 = {
-        'depth': [0, 2, 5],
-        'res': [2e14, 50, 0.1, 50],
-        'aniso': [1, 1.2, 1, 1],
-        'epermH': [0, 1, 1.1, 1],
-        'epermV': [0, 1, 1, 1.1],
-        'mpermH': [1, 1, 1.5, 1],
-        'mpermV': [1, 1.1, 1, 1],
+        "depth": [0, 2, 5],
+        "res": [2e14, 50, 0.1, 50],
+        "aniso": [1, 1.2, 1, 1],
+        "epermH": [0, 1, 1.1, 1],
+        "epermV": [0, 1, 1, 1.1],
+        "mpermH": [1, 1, 1.5, 1],
+        "mpermV": [1, 1.1, 1, 1],
     }
 
     IP1, Q1 = model.ip_and_q(**model1, **system)
@@ -1104,26 +1104,34 @@ def test_ip_and_q(capsys):
     assert_allclose(IP1*1e3, IP2)
     assert_allclose(Q1*1e3, Q2)
 
+    IP3, Q3 = model.ip_and_q(**model1, **{**system, 'ab': 44})
+    assert_allclose(IP3, [-1.3563986, -1.38122726, -7.4560389])
+    assert_allclose(Q3, [-0.00299279869, -0.293418700, -3.40236993])
+
+    IP4, Q4 = model.ip_and_q(**model1, **{**system, 'ab': 55})
+    assert_allclose(IP4, [-10.74321713, -10.69222953, 3.62780845])
+    assert_allclose(Q4, [0.00663009305, 0.651225190, 9.52444857])
+
     # Test errors
     with pytest.raises(ValueError, match="Only implemented for magnetic"):
-        model.ip_and_q(**model1, **{**system, 'ab': 13})
+        model.ip_and_q(**model1, **{**system, "ab": 13})
 
     with pytest.raises(ValueError, match="Only implemented for frequency"):
-        model.ip_and_q(**model1, **{**system, 'signal': 0})
+        model.ip_and_q(**model1, **{**system, "signal": 0})
 
     with pytest.raises(ValueError, match="Only implemented for frequency"):
-        model.ip_and_q(**model1, **{**system, 'signal': 1})
+        model.ip_and_q(**model1, **{**system, "signal": 1})
 
     # Fullspace - no secondary field -> zeros
-    model2 = {'depth': [0], 'res': [50, 50]}
-    IP3, Q3 = model.ip_and_q(**model2, **system)
-    assert_allclose(IP3, 0.0)
-    assert_allclose(Q3, 0.0)
+    model2 = {"depth": [0], "res": [50, 50]}
+    IP5, Q5 = model.ip_and_q(**model2, **system)
+    assert_allclose(IP5, 0.0)
+    assert_allclose(Q5, 0.0)
 
     _, _ = capsys.readouterr()
-    _, _ = model.ip_and_q(**model1, **{**system, 'src': [0, 0, 0]})
+    _, _ = model.ip_and_q(**model1, **{**system, "src": [0, 0, 0]})
     out, _ = capsys.readouterr()
-    assert 'has not been tested' in out
+    assert "This function is experimental" in out
 
 
 def test_fem():