# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# ### Model the DC field of a dipole array above the sea floor

# Here's a sketch:
# ![sketch](GE_DC_Sketch_B.png)

import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics import ert

# A few definitions
waterDepth = 40
ltx = 3.6 # length of transmitter dipole
lrx = 3.36 # length of receiver dipole
htx = 0.3
htx2 = 10.0
dztxrx = 0.8 # fixed distance
resWater = 0.2838 # North Sea water resistivity
current = 10.0 # current in A

# make a world with two layers, the first being 40m thick
world = mt.createWorld(start=[-100, -100], end=[100, 0], layers=[-waterDepth])
pg.show(world, markers=True);

# create the ERT data (cental to origin)
shm = pg.DataContainerERT()
shm.createSensor([-ltx/2, -waterDepth+htx])
shm.createSensor([+ltx/2, -waterDepth+htx])
shm.createSensor([-lrx/2, -waterDepth+htx+dztxrx])
shm.createSensor([+lrx/2, -waterDepth+htx+dztxrx])
shm.createFourPointData(0, 0, 1, 2, 3)
shm["k"] = ert.geometricFactors(shm, dim=2)
print(shm, shm["k"])

for el in shm.sensors():
    world.createNode(el, marker=-99)
    world.createNode(el+[0, 0.03])

mesh = mt.createMesh(world, quality=34.3)
print(mesh)
pg.show(mesh, markers=True, showMesh=True);

for resSubsurface in [1, 2, 3, 4]:
    resMap = [[1, resWater], [2, resSubsurface]]
    data = ert.simulate(mesh=mesh, scheme=shm, res=resMap)
    voltage = data["rhoa"] / data["k"] * current
    print(resSubsurface, data["rhoa"][0], voltage[0])

# +
# Now we do the same for higher floor distance
shm = pg.DataContainerERT()
shm.createSensor([-ltx/2, -waterDepth+htx2])
shm.createSensor([+ltx/2, -waterDepth+htx2])
shm.createSensor([-lrx/2, -waterDepth+htx2+dztxrx])
shm.createSensor([+lrx/2, -waterDepth+htx2+dztxrx])
shm.createFourPointData(0, 0, 1, 2, 3)
shm["k"] = ert.geometricFactors(shm, dim=2)
world = mt.createWorld(start=[-100, -100], end=[100, 0], layers=[-waterDepth])
for el in shm.sensors():
    world.createNode(el, marker=-99)
    world.createNode(el+[0, 0.03])

mesh = mt.createMesh(world, quality=34.3)
pg.show(mesh, markers=True, showMesh=True); # refined differently
# -

resSubsurface = 1.0  # calibration in water
resMap = [[1, resWater], [2, resSubsurface]]
data2 = ert.simulate(mesh=mesh, scheme=shm, res=resMap)
voltage = data2["rhoa"] / data2["k"] * current
print(resSubsurface, data2["rhoa"][0], voltage[0])

StarcoderdataJupyterScriptsDedupFiltered