三维情况下，球坐标系中复杂介质的走时场 及 射线追踪 (修正版)

在原有的示例中，\(\theta\) 的范围设置成[-90, 90]，这是不太合适的，尽管程序能计算出结果。

球坐标系中三个坐标的范围为：

\(r \in [0, +\infty ]\)

\(\theta \in [0, \pi]\)

\(\phi \in [0, 2\pi]\)

import pyfmm
import numpy as np
import matplotlib.pyplot as plt

# 三维情况
rarr = np.arange(20, 50, 0.2)
tarr = np.deg2rad(np.arange(0, 180, 0.5))
parr = np.deg2rad(np.arange(0, 5, 0.5))

# 慢度场
slw  = np.ones((len(rarr), len(tarr), len(parr)), dtype='f')

# 随意添加一些异常
slw -= 0.3*np.cos(3*tarr)[None,:,None]*np.sin(5*(rarr - np.min(rarr))/(np.max(rarr) - np.min(rarr)))[:,None,None]

# 慢度分布
fig, ax = plt.subplots(1, 1, subplot_kw=dict(projection='polar'))
pcm = ax.pcolormesh(tarr, rarr, slw[:, :, 6])
# ax.set_thetalim([-np.pi/2, np.pi/2])
ax.set_thetalim([0, np.pi])
ax.set_rorigin(0.0)
fig.colorbar(pcm)

<matplotlib.colorbar.Colorbar at 0x7ff41ba37820>

srcloc = [45, np.deg2rad(10.0), np.deg2rad(2.0)]

# FMM解
TT = pyfmm.travel_time_source(
    srcloc,
    rarr, tarr, parr, slw, sphcoord=True)

# 射线追踪
rcvloc = [45, np.deg2rad(178.0), np.deg2rad(4.0)]

travt, rays = pyfmm.raytracing(
    TT, srcloc, rcvloc, rarr, tarr, parr, 0.1, sphcoord=True)
print(f"traveltime={travt:.3f}")

traveltime=92.476

# 设置phi方向的索引，观察对应切片
idx = 4
fig, ax = plt.subplots(1, 1, subplot_kw=dict(projection='polar'))
cs = ax.contour(tarr, rarr, TT[:, :, idx], levels=30, linewidths=0.5)
ax.set_thetalim([0, np.pi])
ax.set_rorigin(0.0)
ax.clabel(cs)

ax.plot(rays[:,1], rays[:,0], c='r', lw=0.8)

[<matplotlib.lines.Line2D at 0x7ff418580e50>]

# Fast Sweeping
FSMTT = pyfmm.travel_time_source(
    srcloc,
    rarr, tarr, parr, slw, sphcoord=True, useFSM=True, FSMparallel=True)

travt, rays = pyfmm.raytracing(
    FSMTT, srcloc, rcvloc, rarr, tarr, parr, 0.1, sphcoord=True)
print(f"traveltime={travt:.3f}")

fig, ax = plt.subplots(1, 1, subplot_kw=dict(projection='polar'))
cs = ax.contour(tarr, rarr, FSMTT[:, :, idx], levels=30, linewidths=0.5)
ax.set_thetalim([0, np.pi])
ax.set_rorigin(0.0)
ax.clabel(cs)

ax.plot(rays[:,1], rays[:,0], c='r', lw=0.8)

[WARNING] For parallel FSM, maxLoops must set at least 2 (already changed).

traveltime=92.139

[<matplotlib.lines.Line2D at 0x7ff418618790>]