Merge 7ffddd0 into e47127c

CHANGELOG.rst

@@ -13,7 +13,14 @@ v2.0.x
 latest
 ------
 
-- Documentation: Expanded note on FFTLog.
+- One can define new ``+np.infty`` as interface. Only use-case is to enforce a
+  coordinate system in a two-layer case with an interface at ``z`` (see example
+  coordinate system in the educational section of the gallery).
+
+- Documentation:
+
+  - Expanded note on FFTLog.
+  - Expanded note on coordinate system.
 
 - Maintenance:
 

empymod/utils.py

@@ -728,8 +728,13 @@ def check_model(depth, res, aniso, epermH, epermV, mpermH, mpermV, xdirect,
     # => The top-layer (-infinity to first interface) is layer 0.
     if depth.size == 0:
         depth = np.array([-np.infty, ])
-    elif depth[0] != -np.infty:
-        depth = np.insert(depth, 0, -np.infty)
+    else:
+        if depth[0] != -np.infty:
+            depth = np.r_[-np.infty, depth]
+
+        # Remove +np.infty (can be used to define 2-layer coordinate system).
+        if depth[-1] == np.infty:
+            depth = depth[:-1]
 
     # Check if the user provided a model for etaH/etaV/zetaH/zetaV
     if isinstance(res, dict):

examples/educational/coordinate_system.py

@@ -44,25 +44,17 @@
 coordinate systems:
 
 - The interfaces (``depth``) have to be defined continuously increasing or
-  decreasing, either from lowest to highest or the other way around. E.g., a
-  simple five-layer model with the sea-surface at 0 m, a 100 m water column,
-  and a target of 50 m 900 m below the seafloor can be defined in four ways:
+  decreasing. E.g., think of a simple five-layer model with the sea-surface at
+  0 m, a 100 m water column, and a target of 50 m 900 m below the seafloor,
+  with the following resistivities:
 
-  - ``[0, 100, 1000, 1050]`` -> LHS low to high
-  - ``[0, -100, -1000, -1050]`` -> RHS high to low
-  - ``[1050, 1000, 100, 0]`` -> LHS high to low
-  - ``[-1050, -1000, -100, 0]`` -> RHS low to high
+  - ``res = [1e12, 0.3, 1, 50, 1]``: air, water, background, target,
+    background.
 
-- The above point affects also all model parameters (``res``, ``aniso``,
-  ``{e;m}perm{H;V}``). E.g., for the above five-layer example this would be
+  We can now define the depths in two different ways:
 
-  - ``res = [1e12, 0.3, 1, 50, 1]`` -> LHS low to high
-  - ``res = [1e12, 0.3, 1, 50, 1]`` -> RHS high to low
-  - ``res = [1, 50, 1, 0.3, 1e12]`` -> LHS high to low
-  - ``res = [1, 50, 1, 0.3, 1e12]`` -> RHS low to high
-
-  Note that in a two-layer scenario the values are always taken as low-to-high
-  (as it is not possible to detect the direction from only one interface).
+  - ``depth = [0, 100, 1000, 1050]``: **LHS** (+z down)
+  - ``depth = [0, -100, -1000, -1050]``: **RHS** (+z up)
 
 - A source or a receiver *exactly on* a boundary is taken as being in the lower
   layer. Hence, if :math:`z_\rm{rec} = z_0`, where :math:`z_0` is the surface,
@@ -71,10 +63,51 @@
   receiver is taken as in the sea in the LHS, but as in the subsurface in the
   RHS. This can be avoided by never placing it exactly on a boundary, but
   slightly (e.g., 1 mm) in the layer where you want to have it.
-- In ``bipole``, the ``dip`` switches sign. Correspondingly in ``dipole``, the
-  ``ab``'s containing vertical directions switch the sign for each vertical
-  component.
-- Sign switches also occur for magnetic sources or receivers.
+
+- Sign switches:
+
+  - In ``bipole``, the ``dip`` switches sign.
+  - In ``dipole``, the ``ab``'s containing vertical directions switch the sign
+    for each vertical component.
+  - Sign switches also occur for magnetic sources or receivers.
+
+  These sign switches multiply; e.g., ``ab=34`` has no sign switch, as the
+  switch due to the vertical source and the switch due to the magnetic receiver
+  cancel each other. Here an overview of the sign switches (--) for ``dipole``
+  when changing from the default LHS to a RHS:
+
+  +---------------+-------+------+------+------+------+------+------+
+  |                       | electric  source   | magnetic source    |
+  +===============+=======+======+======+======+======+======+======+
+  |                       | **x**| **y**| **z**| **x**| **y**| **z**|
+  +---------------+-------+------+------+------+------+------+------+
+  |               | **x** |      |      |  --  |  --  |  --  |      |
+  + **electric**  +-------+------+------+------+------+------+------+
+  |               | **y** |      |      |  --  |  --  |  --  |      |
+  + **receiver**  +-------+------+------+------+------+------+------+
+  |               | **z** |  --  |  --  |      |      |      |  --  |
+  +---------------+-------+------+------+------+------+------+------+
+  |               | **x** |  --  |  --  |      |      |      |  --  |
+  + **magnetic**  +-------+------+------+------+------+------+------+
+  |               | **y** |  --  |  --  |      |      |      |  --  |
+  + **receiver**  +-------+------+------+------+------+------+------+
+  |               | **z** |      |      |  --  |  --  |  --  |      |
+  +---------------+-------+------+------+------+------+------+------+
+
+- The depths can also be defined the other way around, starting at +/-1050
+  going to 0. You have to change the order of all model parameters (``res``,
+  ``aniso``, ``{e;m}perm{H;V}``) accordingly.
+
+
+.. note::
+
+  In a two-layer scenario with only one ``depth`` it always assumes a **LHS**,
+  as it is not possible to detect the direction from only one interface. To
+  force any of the other system you can define ``+/-np.infty`` for the
+  down-most interface at the appropriate place:
+
+  - ``0`` or ``[0, np.infty]``: **LHS** (+z down; default)
+  - ``[0, -np.infty]``: **RHS** (+z up)
 
 
 In this example we first create a sketch of the LHS and RHS for visualization,

tests/test_model.py

@@ -700,6 +700,77 @@ def test_all_depths():
     assert_allclose(lhs_l2h, rhs_h2l)
 
 
+def test_coordinate_systems():
+    srcLHS = (0, 0, -10)
+    srcRHS = (0, 0, +10)
+
+    x = np.arange(1, 11)*1000
+    recLHS = (x, x, +3)
+    recRHS = (x, x, -3)
+
+    air, hs, tg = 2e14, 100, 1000
+    z0, z1, z2 = 0, 10, 20
+
+    inp = {'freqtime': 1, 'verb': 1}
+    inpLHS = {'src': srcLHS, 'rec': recLHS}
+    inpRHS = {'src': srcRHS, 'rec': recRHS}
+
+    for ab in [11, 31, 23, 33, 25, 35, 16, 66, 51, 61, 43, 63, 44, 65, 56, 66]:
+        # Sign switches occur for each z-component; each m-component
+        sign = 1
+        if ab % 10 > 3:   # If True: magnetic src
+            sign *= -1
+        if ab // 10 > 3:  # If True: magnetic rec
+            sign *= -1
+        if str(ab)[0] in ['3', '6']:  # Vertical source component
+            sign *= -1
+        if str(ab)[1] in ['3', '6']:  # Vertical receiver component
+            sign *= -1
+        inp['ab'] = ab
+
+        # # 2-layer case
+
+        # Default/original: LHS low to high
+        orig = dipole(depth=z0, res=[air, hs], **inpLHS, **inp)
+
+        # Alternatives LHS: low to high and high to low
+        LHSl2h = dipole(depth=[z0, np.infty], res=[air, hs], **inpLHS, **inp)
+        LHSh2l = dipole(depth=[np.infty, z0], res=[hs, air], **inpLHS, **inp)
+
+        assert_allclose(orig, LHSl2h)
+        assert_allclose(orig, LHSh2l)
+
+        # Alternatives LHS: low to high and high to low
+        RHSlth = sign*dipole(
+            depth=[-np.infty, -z0], res=[hs, air], **inpRHS, **inp)
+        RHSh2l = sign*dipole(
+            depth=[-z0, -np.infty], res=[air, hs], **inpRHS, **inp)
+
+        assert_allclose(orig, RHSlth)
+        assert_allclose(orig, RHSh2l)
+
+        # # 4-layer case
+
+        # Default/original: LHS low to high
+        orig = dipole(
+            depth=[z0, z1, z2], res=[air, hs, tg, hs], **inpLHS, **inp)
+
+        # Alternatives LHS: low to high and high to low
+        LHSh2l = dipole(
+            depth=[z2, z1, z0], res=[hs, tg, hs, air], **inpLHS, **inp)
+
+        assert_allclose(orig, LHSh2l)
+
+        # Alternatives LHS: low to high and high to low
+        RHSlth = sign*dipole(
+            depth=[-z2, -z1, -z0], res=[hs, tg, hs, air], **inpRHS, **inp)
+        RHSh2l = sign*dipole(
+            depth=[-z0, -z1, -z2], res=[air, hs, tg, hs], **inpRHS, **inp)
+
+        assert_allclose(orig, RHSlth)
+        assert_allclose(orig, RHSh2l)
+
+
 class TestLoop:
     # Loop is a subset of bipole, with a frequency-dependent factor at the
     # frequency level.