# Checkout www.pygimli.org for more examples
Refraction in 3D#
This example shows refracted ray paths in a three-dimensional vertical gradient medium.
Note
This is a placeholder/proof-of-concept. The code should be refactored partly to `tt.showRayPaths()`
import numpy as np
import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics import traveltime as tt
from pygimli.viewer.pv import drawSensors
import pyvista
Build mesh.
depth = 15
width = 30
plc = mt.createCube(size=[width, width, depth], pos=[0, 0, -depth/2], area=5)
n_sensors = 8
sensors = np.zeros((n_sensors, 3))
sensors[0, 0] = 15
sensors[0, 1] = -10
sensors[1:, 0] = -15
sensors[1:, 1] = np.linspace(-15, 15, n_sensors - 1)
for pos in sensors:
plc.createNode(pos)
mesh = mt.createMesh(plc)
mesh.createSecondaryNodes(1)
Create vertical gradient model.
vel = 300 + -pg.z(mesh.cellCenters()) * 100
pg.show(mesh, vel, label=pg.utils.unit("vel"), showMesh=True)
Set-up data container.
data = tt.createRAData(sensors)
data.markInvalid(data("s") > 1)
data.set("t", np.zeros(data.size()))
data.removeInvalid()
Do raytracing.
# fop = pg.core.TravelTimeDijkstraModelling(mesh, data)
fop = tt.TravelTimeModelling()
fop.setData(data)
fop.setMesh(mesh)
print(fop.mesh())
# This is to show single raypaths.
graph = fop.createGraph(1 / vel)
dij = tt.Dijkstra(graph)
dij.setStartNode(mesh.findNearestNode([15, -10, 0]))
rays = []
for receiver in sensors[1:]:
ni = dij.shortestPathTo(mesh.findNearestNode(receiver))
pos = mesh.positions(withSecNodes=True)[ni]
segs = np.zeros((len(pos), 3))
segs[:,0] = pg.x(pos)
segs[:,1] = pg.y(pos)
segs[:,2] = pg.z(pos)
rays.append(segs)
Plot final ray paths.
pl, _ = pg.show(mesh, style='wireframe', line_width=0.1,
hold=True)
drawSensors(pl, sensors, diam=0.5, color='yellow')
for ray in rays:
for i in range(len(ray) - 1):
start = tuple(ray[i])
stop = tuple(ray[i + 1])
line = pyvista.Line(start, stop)
pl.add_mesh(line, color='green', line_width=2)
pl.show()