Source code for custEM.misc.synthetic_definitions

# -*- 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)