add read_kernels_freqs() function, add "kernel frequency response" in Tutorial

docs/source/Tutorial/index.rst

@@ -83,7 +83,7 @@
    dynamic/strain_stress
    static/static_strain_stress
    integ_converg/integ_converg
-
+   kernel/kernel
 
 
 

docs/source/Tutorial/integ_converg/integ_converg.rst

@@ -1,3 +1,5 @@
+.. _integ_converg_rst:
+
 积分收敛性
 ===================
 
@@ -199,13 +201,18 @@ C和Python导出的核函数文件是一致的，底层调用的是相同的函
             :start-after: BEGIN SGRN
             :end-before: END SGRN
 
+----------------------------------
+
 
-+ 只使用离散波数积分
++ **只使用离散波数积分**
 
 .. image:: run/DC_20_0.1_static.png
     :align: center
 
-+ 使用峰谷平均法
+----------------------------------
+
+
++ **使用峰谷平均法**
 
 .. image:: run/DC_20_0.1_ptam_static.png
     :align: center

docs/source/Tutorial/kernel/kernel.rst

@@ -0,0 +1,82 @@
+核函数频谱图 
+=============
+
+:Author: Zhu Dengda
+:Email:  zhudengda@mail.iggcas.ac.cn
+
+
+-----------------------------------------------------------
+
+在 :ref:`integ_converg_rst` 部分介绍了程序可以导出积分过程中不同频率的核函数值，这使得我们可以绘制核函数的频谱图 [#]_ ，以观察其中的频散特征 。
+
+以薄层模型为例，
+
+.. literalinclude:: run/mod1
+    :language: text 
+
+.. [#] 感谢席超强博士 `@xichaoqiang <https://github.com/xichaoqiang>`_ 的建议。
+
+
+生成核函数文件
+----------------------
+
+.. note:: 
+
+    与之前的示例相比，这里添加了控制波数积分的参数，以显式控制核函数的采样间隔和范围。但这并非是计算格林函数的最佳参数，故此时格林函数计算结果很可能是不收敛的。
+
+
+.. tabs:: 
+
+    .. group-tab:: C 
+
+        .. literalinclude:: run/run.sh
+            :language: bash 
+            :start-after: BEGIN GRN
+            :end-before: END GRN
+
+
+    .. group-tab:: Python
+
+        .. literalinclude:: run/run.py
+            :language: python
+            :start-after: BEGIN GRN
+            :end-before: END GRN
+
+
+读取，采样，获得核函数频谱
+----------------------------------
+Python端提供了 :py:func:`pygrt.utils.read_kernels_freqs` 函数来完成所有频率的核函数读取以及从波数 ``k`` 插值到特定速度点 ``v`` 的工作 ( :math:`k=\omega/v` )，最后获得保存有核函数频谱的字典。
+
+.. literalinclude:: run/run.py
+    :language: python
+    :start-after: BEGIN read
+    :end-before: END read
+
+绘制核函数频谱图
+------------------------
+以下将使用Python进行图件绘制。
+
+.. literalinclude:: run/run.py
+    :language: python
+    :start-after: BEGIN plot
+    :end-before: END plot
+
+----------------------------------
+
++ **虚部**
+
+.. image:: run/imag.png
+    :align: center
+
+----------------------------------
+
++ **实部**
+
+.. image:: run/real.png
+    :align: center
+
+
+----------------------------------
+
+
+从图上可看到， **核函数频谱的虚部峰值正是频散曲线的位置** ，其中 :math:`q_m，w_m` 呈现的频散特征一致，只是不同震源的核函数幅值不同，这对应的是 **Rayleigh波频散** ，而 :math:`v_m` 对应的是 **Love波频散**。

docs/source/Tutorial/kernel/run/mod1

@@ -0,0 +1,4 @@
+0.01    1.5   0.18   1.78   1e4  1e4
+0.01    1.7   0.35   1.85   1e4  1e4
+0.02    1.6   0.25   1.80   1e4  1e4
+0.00    2.0   0.6    1.94   1e4  1e4

docs/source/Tutorial/kernel/run/run.py

@@ -0,0 +1,82 @@
+# -----------------------------------------------------------------
+# BEGIN GRN
+import numpy as np
+import matplotlib.pyplot as plt
+from typing import Union
+import pygrt 
+
+modarr = np.loadtxt("mod1")
+
+pymod = pygrt.PyModel1D(modarr, depsrc=0.01, deprcv=0.0)
+
+# 不指定statsidx，默认输出全部频率点的积分过程文件
+# vmin_ref 显式给定参考速度（用于定义波数积分上限），避免使用PTAM
+# Length 给定波数积分间隔dk
+_ = pymod.compute_grn(distarr=[1], nt=500, dt=0.02, vmin_ref=0.1, Length=20, statsfile="pygrtstats")
+# END GRN
+# -----------------------------------------------------------------
+
+# -----------------------------------------------------------------
+# BEGIN read
+# 指定待采样的速度数组
+vels = np.arange(0.1, 0.6, 0.001)
+
+# 读取所有频率的核函数，并插值到vels
+# 不指定ktypes，默认返回全部核函数，均以2D数组的形式保存，shape=(nfreqs, nvels)
+kerDct = pygrt.utils.read_kernels_freqs("pygrtstats", vels)
+print(kerDct.keys())
+# dict_keys(['_vels', '_freqs', 'EXP_q0', 'EXP_w0', 'VF_q0', 'VF_w0', 'HF_q1', 'HF_w1', 'HF_v1', 'DC_q0', 'DC_w0', 'DC_q1', 'DC_w1', 'DC_v1', 'DC_q2', 'DC_w2', 'DC_v2'])
+# END read
+# -----------------------------------------------------------------
+
+
+# -----------------------------------------------------------------
+# BEGIN plot
+# 绘制图像
+def plot_kernel(kerDct:dict, RorI:bool, out:Union[str,None]=None):
+    funcRorI = np.real if RorI else np.imag
+
+    ktypes = []
+    for key in kerDct:
+        if key[0] == '_':
+            continue
+        ktypes.append(key)
+
+    srctypes = ['EXP0', 'VF0', 'HF1', 'DC0', 'DC1', 'DC2']
+
+    vels = kerDct['_vels']
+    freqs = kerDct['_freqs']
+
+    fig = plt.figure(figsize=(12, 3*len(srctypes)))
+    gs = fig.add_gridspec(len(srctypes), 3)
+    qwvLst = ['q', 'w', 'v']
+    for ik, key in enumerate(ktypes):
+        srctype = key.split("_")[0]+key.split("_")[1][1:]
+        qwv = key.split("_")[1][0]
+
+        ax = fig.add_subplot(gs[srctypes.index(srctype), qwvLst.index(qwv)])
+
+        # 对不同速度间取归一化
+        data = kerDct[key].copy()
+        data[...] = data/np.max(np.abs(data), axis=1)[:,None]
+
+        pcm = ax.pcolormesh(freqs, vels, np.abs(funcRorI(data)).T, vmin=0, vmax=1, shading='nearest')
+        ax.set_xlabel("Frequency (Hz)")
+        ax.set_ylabel("Velocity (km/s)")
+        ax.set_title(key)
+        fig.colorbar(pcm, ax=ax)
+
+    if RorI:
+        fig.suptitle("Real parts of Kernels", fontsize=20, x=0.5, y=0.99)
+    else:
+        fig.suptitle("Imag parts of Kernels", fontsize=20, x=0.5, y=0.99)
+
+    if out is not None:
+        fig.tight_layout()
+        fig.savefig(out, dpi=100)
+
+
+plot_kernel(kerDct, False, "imag.png")
+plot_kernel(kerDct, True, "real.png")
+# END plot
+# -----------------------------------------------------------------

docs/source/Tutorial/kernel/run/run.sh

@@ -0,0 +1,10 @@
+#!/bin/bash 
+
+# -----------------------------------------------------------------
+# BEGIN GRN
+# -S-1 表示输出所有频率点的核函数
+# -V0.1 显式给定参考速度（用于定义波数积分上限），避免使用PTAM
+# -L20 定义波数积分间隔dk
+grt -Mmod1 -D0/0 -N500/0.02/0.8 -OGRN -R1 -V0.1 -S-1 -L20
+# END GRN
+# -----------------------------------------------------------------

docs/source/Tutorial/run_all.sh

@@ -42,4 +42,12 @@ if [[ $1 == '5' || $1 == '' ]]; then
     ./run.sh
     python run.py
     cd -
+fi
+
+if [[ $1 == '6' || $1 == '' ]]; then
+    cd kernel/run
+    chmod +x *.sh
+    ./run.sh
+    python run.py
+    cd -
 fi

pygrt/utils.py

@@ -23,6 +23,7 @@
 from scipy.signal import oaconvolve
 from scipy.fft import rfft, irfft
 from scipy.special import jv
+from scipy.interpolate import interpn
 import math 
 import os
 import glob
@@ -48,6 +49,7 @@
     "stream_diff",
     "stream_write_sac",
 
+    "read_kernels_freqs",
     "read_statsfile",
     "read_statsfile_ptam",
     "plot_statsdata",
@@ -1250,6 +1252,75 @@ def read_statsfile(statsfile:str):
     return data
 
 
+def read_kernels_freqs(statsdir:str, vels:np.ndarray, ktypes:Union[List[str],None]=None):
+    r"""
+        读取statsdir目录下所有频率（除了零频）的积分过程文件（K_开头），再将核函数线性插值到指定的速度数组vels
+
+        :param        statsdir:     存储积分过程文件的目录
+        :param        vels:         待插值的速度数组(km/s)，必须正序
+        :param        ktype:        指定返回一系列的核函数名称，如EXP_q0，DC_w2等，默认返回全部
+
+        :return:
+            - **kerDct**  -   字典格式的核函数插值结果
+    """
+
+    if not np.all(np.diff(vels) > 0):
+        raise ValueError("vels must be in ascending order.")
+
+    KLst = np.array(glob.glob(os.path.join(statsdir, "K_*")))
+    freqs = np.array([float(s.split("_")[-1]) for s in KLst])
+    # 根据freqs排序
+    _idx = np.argsort(freqs)
+    freqs[:] = freqs[_idx]
+    KLst[:] = KLst[_idx]
+    del _idx 
+
+    # 去除零频
+    if freqs[0] == 0.0:
+        freqs = freqs[1:]
+        KLst = KLst[1:]
+
+    kerDct = {}
+    kerDct['_vels'] = vels.copy()
+    kerDct['_freqs'] = freqs.copy()
+
+    for i in range(len(freqs)):
+        Kpath = KLst[i]
+        freq = freqs[i]
+        w = 2*np.pi*freq
+
+        data = read_statsfile(Kpath)
+        v = w/data['k']
+
+        # 检查v范围
+        v1 = np.min(v)
+        v2 = np.max(v)
+        if v1 > vels[0] or v2 < vels[-1]:
+            raise ValueError(f"In freq={freq:.5e}, minV={v1:.5e}, maxV={v2:.5e}, insufficient wavenumber samples"
+                              " to interpolate on vels.")
+
+        for key in data.dtype.names:
+            if key == 'k':
+                continue 
+            if (ktypes is not None) and (key not in ktypes):
+                continue 
+
+            if key not in kerDct.keys():
+                kerDct[key] = []
+
+            # 如果越界会报错
+            F = interpn((v,), data[key], vels)
+            kerDct[key].append(F)
+
+    # 将每个核函数结果拼成2D数组
+    for key in kerDct.keys():
+        if key[0] == '_':
+            continue
+        kerDct[key] = np.vstack(kerDct[key])
+
+    return kerDct
+
+
 def read_statsfile_ptam(statsfile:str):
     '''
         读取单个频率下峰谷平均法的记录文件