Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
from SimPEG import maps
from SimPEG.potential_fields import gravity
# sphinx_gallery_thumbnail_number = 2
#############################################
# Defining Topography
# -------------------
#
# Surface topography is defined as an (N, 3) numpy array. We create it here but
# the topography could also be loaded from a file.
#
[x_topo, y_topo] = np.meshgrid(np.linspace(-200, 200, 41), np.linspace(-200, 200, 41))
z_topo = -15*np.exp(-(x_topo**2 + y_topo**2) / 80**2)
x_topo, y_topo, z_topo = mkvc(x_topo), mkvc(y_topo), mkvc(z_topo)
xyz_topo = np.c_[x_topo, y_topo, z_topo]
#############################################
# Defining the Survey
# -------------------
#
# Here, we define survey that will be used for the forward simulation. Gravity
# surveys are simple to create. The user only needs an (N, 3) array to define
# the xyz locations of the observation locations, and a list of field components
# which are to be measured.
#
# Define the observation locations as an (N, 3) numpy array or load them.
x = np.linspace(-80., 80., 17)
y = np.linspace(-80., 80., 17)
h = [(dh, nbc)]
mesh = TreeMesh([h, h, h], x0='CCC')
# Mesh refinement based on topography
mesh = refine_tree_xyz(
mesh, topo_xyz, octree_levels=[0, 0, 0, 1], method='surface', finalize=False
)
# Mesh refinement near transmitters and receivers
mesh = refine_tree_xyz(
mesh, rx_locs, octree_levels=[2, 4], method='radial', finalize=False
)
# Refine core mesh region
xp, yp, zp = np.meshgrid([-300., 300.], [-300., 300.], [-300., 0.])
xyz = np.c_[mkvc(xp), mkvc(yp), mkvc(zp)]
mesh = refine_tree_xyz(
mesh, xyz, octree_levels=[0, 2, 4], method='box', finalize=False
)
mesh.finalize()
###############################################################
# Create Resistivity Model and Mapping for OcTree Mesh
# ----------------------------------------------------
#
# Here, we define the electrical properties of the Earth as a resistivity
# model. The model consists of a long vertical conductive pipe within a more
# resistive background.
#
# Log-Resistivity in log[Ohm m]
from SimPEG import maps
from SimPEG.potential_fields import gravity
# sphinx_gallery_thumbnail_number = 2
#############################################
# Defining Topography
# -------------------
#
# Surface topography is defined as an (N, 3) numpy array. We create it here but
# the topography could also be loaded from a file.
#
[x_topo, y_topo] = np.meshgrid(np.linspace(-200, 200, 41), np.linspace(-200, 200, 41))
z_topo = -15*np.exp(-(x_topo**2 + y_topo**2) / 80**2)
x_topo, y_topo, z_topo = mkvc(x_topo), mkvc(y_topo), mkvc(z_topo)
xyz_topo = np.c_[x_topo, y_topo, z_topo]
#############################################
# Defining the Survey
# -------------------
#
# Here, we define survey that will be used for the forward simulation. Gravity
# surveys are simple to create. The user only needs an (N, 3) array to define
# the xyz locations of the observation locations, and a list of field components
# which are to be measured.
#
# Define the observation locations as an (N, 3) numpy array or load them
x = np.linspace(-80., 80., 17)
y = np.linspace(-80., 80., 17)
# sphinx_gallery_thumbnail_number = 2
#############################################
# Topography
# ----------
#
# Here we define surface topography as an (N, 3) numpy array. Topography could
# also be loaded from a file.
#
[x_topo, y_topo] = np.meshgrid(
np.linspace(-200, 200, 41), np.linspace(-200, 200, 41)
)
z_topo = -15*np.exp(-(x_topo**2 + y_topo**2) / 80**2)
x_topo, y_topo, z_topo = mkvc(x_topo), mkvc(y_topo), mkvc(z_topo)
xyz_topo = np.c_[x_topo, y_topo, z_topo]
#############################################
# Defining the Survey
# -------------------
#
# Here, we define survey that will be used for the simulation. Magnetic
# surveys are simple to create. The user only needs an (N, 3) array to define
# the xyz locations of the observation locations, the list of field components
# which are to be modeled and the properties of the Earth's field.
#
# Define the observation locations as an (N, 3) numpy array or load them.
x = np.linspace(-80., 80., 17)
y = np.linspace(-80., 80., 17)
x, y = np.meshgrid(x, y)
# result in a decrease in numerical accuracy for the predicted fields.
#
###############################################################
# Defining Topography
# -------------------
#
# Here we define surface topography as an (N, 3) numpy array. Topography could
# also be loaded from a file. Here we define flat topography, however more
# complex topographies can be considered.
#
xx, yy = np.meshgrid(np.linspace(-3000, 3000, 101), np.linspace(-3000, 3000, 101))
zz = np.zeros(np.shape(xx))
topo_xyz = np.c_[mkvc(xx), mkvc(yy), mkvc(zz)]
#####################################################################
# Create Airborne Survey
# ----------------------
#
# For this example, the survey consists of a uniform grid of airborne
# measurements. To save time, we will compute the response for a minimal
# range of time channels
#
# time channels
n_times = 6
tc = np.logspace(-4, -3, n_times)
# Defining transmitter locations
nbcz = 2**int(np.round(np.log(z_length/dz)/np.log(2.)))
# Define the base mesh
hx = [(dx, nbcx)]
hy = [(dy, nbcy)]
hz = [(dz, nbcz)]
mesh = TreeMesh([hx, hy, hz], x0='CCN')
# Refine based on surface topography
mesh = refine_tree_xyz(
mesh, topo, octree_levels=[2, 2], method='surface', finalize=False
)
# Refine box base on region of interest
xp, yp, zp = np.meshgrid([-100., 100.], [-100., 100.], [-80., 0.])
xyz = np.c_[mkvc(xp), mkvc(yp), mkvc(zp)]
mesh = refine_tree_xyz(
mesh, xyz, octree_levels=[2, 2], method='box', finalize=False
)
mesh.finalize()
##########################################################
# Create Magnetic Vector Intensity Model (MVI)
# --------------------------------------------
#
# Magnetic vector models are defined by three-component effective
# susceptibilities. To create a magnetic vector
# model, we must
#
# 1) Define the magnetic susceptibility for each cell. Then multiply by the
# ----------------------
#
# For this example, the survey consists of a uniform grid of airborne
# measurements. To save time, we will only compute the response for a single
# frequency.
#
# Frequencies being predicted
freq = 100.
# Defining transmitter locations
N = 11
xtx, ytx, ztx = np.meshgrid(
np.linspace(-200, 200, N), np.linspace(-200,200, N), [50]
)
src_locs = np.c_[mkvc(xtx), mkvc(ytx), mkvc(ztx)]
ntx = np.size(xtx)
# Define receiver locations
xrx, yrx, zrx = np.meshgrid(
np.linspace(-200, 200, N), np.linspace(-200,200, N), [30]
)
rx_locs = np.c_[mkvc(xrx), mkvc(yrx), mkvc(zrx)]
source_list = [] # Create empty list to store sources
# Each unique location and frequency defines a new transmitter
for ii in range(ntx):
# Define receivers of different type at each location
bzr = receivers.Point_bSecondary(rx_locs[ii, :], 'z', 'real')
bzi = receivers.Point_bSecondary(rx_locs[ii, :], 'z', 'imag')
# sphinx_gallery_thumbnail_number = 2
#############################################
# Topography
# ----------
#
# Surface topography is defined as an (N, 3) numpy array. We create it here but
# topography could also be loaded from a file.
#
[x_topo, y_topo] = np.meshgrid(
np.linspace(-200, 200, 41), np.linspace(-200, 200, 41)
)
z_topo = -15*np.exp(-(x_topo**2 + y_topo**2) / 80**2)
x_topo, y_topo, z_topo = mkvc(x_topo), mkvc(y_topo), mkvc(z_topo)
xyz_topo = np.c_[x_topo, y_topo, z_topo]
#############################################
# Defining the Survey
# -------------------
#
# Here, we define survey that will be used for the simulation. Magnetic
# surveys are simple to create. The user only needs an (N, 3) array to define
# the xyz locations of the observation locations, the list of field components
# which are to be modeled and the properties of the Earth's field.
#
# Define the observation locations as an (N, 3) numpy array or load them.
x = np.linspace(-80., 80., 17)
y = np.linspace(-80., 80., 17)
x, y = np.meshgrid(x, y)