DOC: add lamb1 (the first-kind Lamb problem)

docs/source/Lamb_problem/index.rst

@@ -0,0 +1,28 @@
+Lamb 问题
+================
+
+:Author: Zhu Dengda
+:Email:  zhudengda@mail.iggcas.ac.cn
+
+-----------------------------------------------------------
+
+北京大学张海明教授已在其编著的 **《地震学中的 Lamb 问题》上下两册** 中对三类 Lamb 问题进行了详细论述，
+并分别在上下两册中给出了频域解和时域解，过程之详细令人叹为观止。这里不再对 Lamb 问题及其解法做过多介绍，
+详见张海明老师的书（强烈推荐！！）
+
+**PyGRT** 程序包本身是在频域求解，这里我对照下册书也编程实现了 Lamb 问题的时域解（目前包括第一类 Lamb 问题，后续会再扩展），
+方便进行对比验证。以下简单示范用法以及绘图。
+
+.. note::
+
+    程序输入的时间为无量纲时间，输出的位移为和阶跃函数卷积后的无量纲位移。
+
+    + 无量纲时间 :math:`\bar{t} = \dfrac{t}{T_S}` ，其中 :math:`T_S = \dfrac{r}{\beta}`
+    + 无量纲位移 :math:`\bar{\mathbf{G}}^H = \pi^2 \mu r \mathbf{G}^H` ，上标 :math:`H` 表示已和阶跃函数卷积。
+
+
+.. toctree::
+   :maxdepth: 1
+
+   lamb1
+

docs/source/Lamb_problem/lamb1.rst

@@ -0,0 +1,54 @@
+第一类 Lamb 问题
+===================
+
+:Author: Zhu Dengda
+:Email:  zhudengda@mail.iggcas.ac.cn
+
+-----------------------------------------------------------
+
+第一类 Lamb 问题是指，在半空间模型中，源点和场点均位于地表，求解场点记录到的位移。
+
+.. tabs::  
+
+   .. group-tab:: C
+
+      C 程序 :command:`grt` 提供了模块 :doc:`/Module/lamb1` 求解第一类 Lamb 问题。
+
+      .. literalinclude:: run/run.sh
+         :language: bash
+         :start-after: BEGIN LAMB1
+         :end-before: END LAMB1
+
+      使用重定向将结果保存到文件 *lamb1.txt* 中，其内容格式类似于
+
+      .. literalinclude:: run/head_lamb1
+         :language: text
+
+   .. group-tab:: Python    
+
+      Python 提供了函数 :func:`solve_lamb1() <pygrt.utils.solve_lamb1>` 求解第一类 Lamb 问题。
+
+      .. literalinclude:: run/lamb1_plot_time.py
+         :language: python
+         :start-after: BEGIN LAMB1
+         :end-before: END LAMB1
+
+最后绘制计算得到的格林函数。
+
+:download:`lamb1_plot_time.py <run/lamb1_plot_time.py>`
+
+.. image:: run/lamb1_time.png
+   :align: center
+
+
+频域解和时域解的对比
+-------------------------------
+
+对比观察可发现频域解在波形突变出有明显的 Gibbs 效应。
+
+:download:`lamb1_plot_freq_time.py <run/lamb1_plot_freq_time.py>`
+
+.. image:: run/lamb1_compare_freq_time.png
+   :align: center
+
+

docs/source/Lamb_problem/run/lamb1_plot_freq_time.py

@@ -0,0 +1,86 @@
+import pygrt
+import numpy as np
+import matplotlib.pyplot as plt
+
+plt.rcParams.update({
+    "font.sans-serif": "Times New Roman",
+    "mathtext.fontset": "cm"
+})
+
+# 定义半无限空间模型 
+Vp = 8     # km/s
+Vs = 4.62  # km/s
+Rho = 3.3  # g/cm^3
+
+# 泊松比
+nu = 0.5 * (1 - 2*(Vs/Vp)**2) / (1 - (Vs/Vp)**2)
+
+# 模型数组，半无限空间
+modarr = np.array([
+    [0.,  Vp, Vs, Rho, 9e10, 9e10],
+])
+
+# 剪切模量
+mu = Vs**2 * Rho
+
+depsrc = 0.0  # 震源深度 km
+deprcv = 0.0  # 台站深度 km
+
+rs = np.array([10]) # 震中距数组，km
+
+nt = 1000     # 总点数，不要求2的幂次
+dt = 0.005    # 采样时间间隔(s)
+
+idx = 0      # 震中距索引
+r = rs[idx]
+
+t = np.arange(0, nt)*dt * Vs/r 
+
+
+pymod = pygrt.PyModel1D(modarr, depsrc, deprcv)  # 整理好的模型对象
+# 计算格林函数频谱
+st = pymod.compute_grn(
+    distarr=rs, 
+    nt=nt, 
+    dt=dt, 
+)[0]
+
+# 卷积阶跃函数
+pygrt.utils.stream_integral(st)
+
+
+# 时域解
+u = pygrt.utils.solve_lamb1(nu, t, 0).reshape(-1, 9)
+u = u[:, [0,2,4,6,8]]
+u.shape
+
+# 频域解
+v = np.zeros_like(u)
+coef = np.pi * np.pi * mu * r 
+v[:,0] = st.select(channel='HFR')[0].data * coef
+v[:,1] = st.select(channel='VFR')[0].data * coef
+v[:,2] = st.select(channel='HFT')[0].data * coef
+v[:,3] = st.select(channel='HFZ')[0].data * coef * (-1)
+v[:,4] = st.select(channel='VFZ')[0].data * coef * (-1)
+
+fig, axs = plt.subplots(5, 2, figsize=(10, 10), sharex=True)
+labels = [r"$\bar{{G}}^H_{11}$", r"$\bar{{G}}^H_{13}$", r"$\bar{{G}}^H_{22}$", r"$\bar{{G}}^H_{31}$", r"$\bar{{G}}^H_{33}$"]
+for i in range(5):
+    ax = axs[i, 0]
+    ax.plot(t, u[:,i])
+    ax.set_ylim(-2, 2)
+    ax.set_xlim(0, 2)
+    ax.grid()
+    ax.text(0.05, 0.92, labels[i], transform=ax.transAxes, ha='left', va='top', fontsize=12)
+
+    ax = axs[i, 1]
+    ax.plot(t, v[:,i])
+    ax.set_ylim(-2, 2)
+    ax.set_xlim(0, 2)
+    ax.grid()
+    ax.text(0.05, 0.92, labels[i], transform=ax.transAxes, ha='left', va='top', fontsize=12)
+
+axs[0,0].set_title("From Time-Domain")
+axs[0,1].set_title("From Frequency-Domain")
+
+fig.savefig("lamb1_compare_freq_time.png", dpi=100, bbox_inches='tight')

docs/source/Lamb_problem/run/lamb1_plot_time.py

@@ -0,0 +1,21 @@
+# BEGIN LAMB1
+import pygrt
+import numpy as np
+
+ts = np.arange(0, 2, 1e-4)
+u = pygrt.utils.solve_lamb1(0.25, ts, 30)
+# END LAMB1
+
+import matplotlib.pyplot as plt
+
+fig, axs = plt.subplots(3, 3, figsize=(10, 5), sharex=True)
+for i in range(3):
+    for j in range(3):
+        ax = axs[i, j]
+        ax.plot(ts, u[:, i, j])
+        ax.set_xlim(0, 2)
+        ax.set_ylim(-2, 2)
+
+        ax.text(0.1, 0.9, rf"$\bar{{G}}^H_{{{i+1}{j+1}}}$", transform=ax.transAxes, ha='left', va='top', fontsize=12)
+
+fig.savefig("lamb1_time.png", dpi=100, bbox_inches='tight')

docs/source/Lamb_problem/run/run.sh

@@ -0,0 +1,9 @@
+#!/bin/bash
+
+# -----------------------------------------------------------------------------------
+# BEGIN LAMB1
+grt lamb1 -P0.25 -T0/2/1e-4 -A30 > lamb1.txt
+# END LAMB1
+# -----------------------------------------------------------------------------------
+
+head -n 10 lamb1.txt > head_lamb1

docs/source/Module/index.rst

@@ -95,6 +95,7 @@
     - :doc:`sac2asc`
     - :doc:`ker2asc`
     - :doc:`travt`
+    - :doc:`lamb1`
 
 
 .. toctree::
@@ -114,4 +115,5 @@
     sac2asc
     ker2asc
     travt
+    lamb1
 

docs/source/Module/lamb1.rst

@@ -0,0 +1,50 @@
+.. include:: common_OPTs.rst_
+
+
+lamb1
+==============
+
+:简介: 使用广义闭合解求解第一类 Lamb 问题
+
+语法
+-----------
+
+**grt lamb1**
+|-P|\ *nu*
+|-T|\ *t1/t2/dt*
+|-A|\ *azimuth*
+[ **-h** ]
+
+
+描述
+--------
+
+第一类 Lamb1 问题指在半空间中，当源点和场点均位于地表时，求场点记录到的位移。
+**lamb1** 模块实现的理论基础来源于《地震学中的 Lamb 问题（下）》。
+结果为无量纲位移（距离物理解还需除 :math:`\pi^2 \mu r` ）和阶跃函数卷积后的结果，
+输出到标准输出。
+
+必选选项
+----------
+
+.. _-P:
+
+**-P**\ *nu*
+    半空间的泊松比 *nu*, 要求范围在 (0, 0.5) 。
+
+.. _-T:
+
+**-T**\ *t1/t2/dt*
+    无量纲时间序列 :math:`\bar{t}` ，根据起点时刻 *t1*, 终点时刻 *t2* 和间隔 *dt* 定义， 
+    :math:`\bar{t} = \dfrac{t}{T_S}` ，其中 :math:`T_S = \dfrac{r}{\beta}` ，:math:`r` 为震中距，:math:`\beta` 为 S 波速度。
+
+.. _-A:
+
+**-A**\ *azimuth*
+    方位角，单位为度。
+
+
+示例
+-------
+
+详见 :doc:`/Lamb_problem/lamb1` 。

docs/source/index.rst

@@ -120,6 +120,7 @@
    :hidden:
    :maxdepth: 1
 
+   Lamb_problem/index
    API/api
    changelog
    copyright

docs/source/run_all.sh

@@ -60,4 +60,13 @@ if [[ $1 == '7' || $1 == '' ]]; then
     python run.py
     python plot.py
     cd -
+fi
+
+if [[ $1 == '8' || $1 == '' ]]; then
+    cd Lamb_problem/run
+    chmod +x *.sh
+    ./run.sh
+    python lamb1_plot_time.py
+    python lamb1_plot_freq_time.py
+    cd -
 fi