# -*- coding: utf-8 -*-
"""
@author: Rochlitz.R
"""
import numpy as np
"""
This file contains several definitions of functions for specifying synthetic
topographies for modeling studies and the benchmark examples.
"""
[docs]
def flat_topo(x, y, h=0.):
"""
Flat topo specified as function for implementation tests.
"""
return(h)
[docs]
def almost_flat_topo(x, y, h=0.1):
"""
Almost flat topo specified as function for mesh generation tests.
"""
return(h + h * np.sin(1e-3 * y))
[docs]
def topo_func1(x, y, h=100.):
"""
Example cosine-topography in y-direction.
"""
return(h * np.cos(3e-3 * (y - 100.)))
[docs]
def topo_func2(x, y, h=200.):
"""
Slope in x-direction.
"""
z_offset = 200.
return(5e-2 * x + z_offset)
[docs]
def topo_func3(x, y, h=200.):
"""
Example sine-topography in y-direction.
"""
z_offset = -2000.
return(h * np.sin(1e-3 * y) + z_offset)
[docs]
def subtopo_func1(x, y, h=2e2):
"""
Example sine-shaped subsurface topography in x-direction.
"""
z_offset = -1000.
return(h * np.sin(4e4 * x) + z_offset)
[docs]
def subtopo_func2(x, y, h=1e-1):
"""
Slope of subsurface layer in y-direction.
"""
z_offset = -3000.
return(h * y + z_offset)
[docs]
def example_4_topo_cos_sin(x, y, h=200.):
"""
Example topography for frequency-domain example 4, x- and y-directed
cosine/sine topography.
"""
z_offset = 0.
return(h * np.cos(1e-3 * (x - 2.2e3)) + (h/2.) * np.sin(3e-3 * y) +
z_offset)
[docs]
def sample_topo_func(x, y, h=100.):
"""
Example sine topography for test scripts.
"""
z_offset = 250.
return(h * np.sin(1e-3 * x) + y * 1e-2 + z_offset)
[docs]
def pyramid_topo(x, y, h=1e3, slope=1., z_off=0., x_off=0., y_off=0.):
"""
Pyramid-shaped topography.
"""
z = np.zeros(len(x))
for j in range(len(x)):
if (x[j] + x_off) < h / slope and (x[j] + x_off) > -h / slope and\
(y[j] + y_off) < h / slope and (y[j] + y_off) > -h / slope:
z[j] = np.min([h - np.abs((x[j] + x_off) / slope),
h - np.abs((y[j] + y_off) / slope)])
return(z + z_off)
[docs]
def roof_topo(x, y, h=1e3, slope=1., z_off=0., x_off=0., ylim=[-1e3, 1e3]):
"""
Pyramid-shaped topography with cutted top (roof).
"""
z = np.zeros(len(x))
for j in range(len(x)):
if (x[j] + x_off) < h / slope and (x[j] + x_off) > -h / slope and\
y[j] >= ylim[0] and y[j] <= ylim[1]:
z[j] = h - np.abs((x[j] + x_off) / slope)
return(z + z_off)
[docs]
def tx_slope(x, y, h=1e1, z_off=0., x_off=50., flank=10., ylim=[-6e1, 6e1]):
"""
Topography description building a slight slope just beneath a loop_r Tx.
"""
slope = h / (2 * x_off)
z = np.zeros(len(x))
for j in range(len(x)):
if x[j] >= -x_off and x[j] <= x_off and \
y[j] >= ylim[0] and y[j] <= ylim[1]:
z[j] = h - np.abs((x[j] + x_off) * slope)
elif x[j] <= -x_off and x[j] >= -x_off - flank and\
y[j] >= ylim[0] and y[j] <= ylim[1]:
z[j] = h - np.abs((x[j] + x_off))
for j in range(len(x)):
if y[j] <= -x_off and y[j] >= -x_off - flank and\
x[j] >= ylim[0] and x[j] <= ylim[1]:
z[j] = np.min([h - np.abs((y[j] + x_off)), z[j]])
if y[j] >= x_off and y[j] <= x_off + flank and\
x[j] >= ylim[0] and x[j] <= ylim[1]:
z[j] = np.min([h - np.abs((y[j] - x_off)), z[j]])
# z[j] = h - np.abs((y[j] - x_off))
return(z + z_off)
[docs]
def valley(x, y, height=1000., z=0.):
"""
Example topography with a valley for testing bathy-tools.
"""
# return(-height * np.cos(stretch*1e-3 * np.sqrt(x**2)) + z +
# -height * np.sin(stretch*2e-3 * np.sqrt(y**2)) + z)
return(-height / 2. * np.cos(0.3e-3 * x) + z -
(np.cos(0.3e-3*y))*100.)
[docs]
def cone_bathy(x, y, h=1000., z=0.):
"""
Example topography with a cone hill for testing bathy-tools.
"""
return(h * np.cos(0.3e-3 * np.sqrt(x**2 + y**2)) + z)
[docs]
def anti_cone_bathy(x, y, h=1000., z=0.):
"""
Example topography with a cone hole for testing bathy-tools.
"""
return(-(h * np.cos(0.3e-3 * np.sqrt(x**2 + y**2)) + z))
[docs]
def double_cone_bathy(x, y, h=1000., sigma=1., zl=0.):
"""
Example topography with two cones for testing bathy-tools.
"""
z = np.zeros(len(x))
for j in range(len(x)):
if np.sqrt((x[j]**2 + y[j]**2)) < 1e2:
z[j] = +(2*h * np.cos(0.3e-3 * np.sqrt(x[j]**2 + y[j]**2)) + zl)
elif x[j] < 1e3 and x[j] > -1e3 and y[j] < 1e3 and y[j] > -1e3:
z[j] = -(h * np.cos(0.3e-3 * np.sqrt(x[j]**2 + y[j]**2)) + zl)
else:
z[j] = +(h * np.cos(0.3e-3 * np.sqrt(x[j]**2 + y[j]**2)) + zl)
return(z)
[docs]
def gaussian_topo(x, y, h=3e9, sigma=1e3, z=0., x_off=0., y_off=0.):
"""
Example topography with hill descriped as gaussian function.
"""
z = h/(sigma**2 * np.pi * 2.) * np.exp(
-((x - x_off)**2 / (2. * sigma**2) +
(y - y_off)**2 / (2. * sigma**2)))
z[z < 1.] = 0.
return(z)
[docs]
def gaussian_topo_hole(x, y, h=3e9, sigma=1e3, x_off=0., y_off=0.):
"""
Example topography with hole descriped as gaussian function.
"""
z = h/(sigma**2 * np.pi * 2.) * np.exp(
-((x - x_off)**2 / (2. * sigma**2) +
(y - y_off)**2 / (2. * sigma**2)))
z[z < 1.] = 0.
for j in range(len(x)):
if x[j] < 2.001e3 and x[j] > -2.001e3 and \
y[j] < 2.001e3 and y[j] > -2.001e3:
z[j] = 0.
return(z)
[docs]
def surface_anomly_outcrop_1():
"""
Example outcrop polygon for surface anomaly.
"""
from custEM.meshgen import meshgen_utils as mu
poly = mu.loop_r(start=[0.5e3, -5e2], stop=[0.6e3, 3e2], n_segs=14)
poly[2:7, 0] += np.array([12., 30., 40., 35., 20.])
poly[9:14, 0] += np.array([12., 30., 40., 35., 20.])
poly = mu.rotate(poly, np.deg2rad(20.), 'z')
return(poly)
[docs]
def surface_anomly_outcrop_2():
"""
Example outcrop polygon for surface anomaly.
"""
from custEM.meshgen import meshgen_utils as mu
poly = mu.loop_r(start=[-4e2, 0.3e3], stop=[-1e2, 1.2e3], n_segs=14)
poly[2:7, 0] += np.array([12., 30., 40., 35., 20.])
poly[9:14, 0] += np.array([12., 30., 40., 35., 20.])
poly = mu.rotate(poly, np.deg2rad(20.), 'z')
return(poly)
[docs]
def surface_anomly_outcrop_3():
"""
Example outcrop polygon for surface anomaly.
"""
from custEM.meshgen import meshgen_utils as mu
poly = mu.loop_r(start=[-3e3, -7e2], stop=[1.5e3, -6e2], nx=20, ny=3)
shift = np.append(np.linspace(5., 100., 15) * np.logspace(0., 0.6, 15),
np.linspace(100., 10., 4) * np.logspace(0.6, 0., 4))
poly[1:20, 1] += shift
poly[1:21, 1] += np.linspace(1., 80., 20)
poly[21:23, 1] += np.array([60., 30.])
poly[24:43, 1] += shift[::-1]
poly = mu.rotate(poly, np.deg2rad(20.), 'z')
return(poly)
[docs]
def example_3_line_tx(x=None):
"""
Deprecated, not used anymore.
"""
# z_offset = 250.
# return(np.sin(1e-3 * x) + y * 1e-2 + z_offset)
if x is None:
x = np.linspace(-5e2, 5e2, 100)
X = np.zeros((len(x), 3))
X[:, 1] = x
for j in range(len(x)):
if x[j] > 100.:
X[j, 0] = 50.
elif x[j] < -200.:
X[j, 0] = -100.
else:
X[j, 0] = 0.5 * X[j, 1]
return(X)
[docs]
def example_2_loop_tx():
"""
Crooked loop Tx description for frequency-domain example 2
"""
from custEM.meshgen import meshgen_utils as mu
loop = mu.loop_r(start=[-4e2, -6e2], stop=[4e2, 7e2], nx=10, ny=10)
np.random.seed(23)
loop[:, 0] += np.random.rand(len(loop)) * 50. - 25.
np.random.seed(42)
loop[:, 1] += np.random.rand(len(loop)) * 100. - 50.
loop[30:, 0] *= np.linspace(1., 2., 10)
return(loop)
[docs]
def borehole_tx():
"""
Description of Tx grounded in two different boreholes.
"""
borehole_tx = np.array([[-5e2, 0., -5e2],
[-5e2, 0., 0.],
[5e2, 0., 0.],
[5e2, 0., -5e2]])
return(borehole_tx)