FIX: update default factor in kmax formula and update related doc

docs/source/Advanced/integ_converg/run_dcm/plot_depth_kernel.py

@@ -13,10 +13,9 @@ def gray_mapping(x):
     return (0.7 - x / 0.7) ** 0.8
 
 def _plot_one(ax:Axes, pattern:str, ktype:str):
-    norm = 0.0
-    yLst = []
-    evdpLst = []
+    plotLst = []
     maxnorm = 0
+    kmax = 0.0
     for path in glob.glob(pattern):
         data = pygrt.utils.read_statsfile(path)
 
@@ -26,16 +25,16 @@ def _plot_one(ax:Axes, pattern:str, ktype:str):
         Farr = data[ktype]
 
         norm = np.linalg.norm(np.real(Farr), np.inf)
-        yLst.append(np.real(Farr))
-        evdpLst.append(evdp)
+        plotLst.append((karr, np.real(Farr), evdp))
         maxnorm = max(maxnorm, norm)
+        kmax = max(kmax, karr[-1])
 
     ax.axhline(y=0, color='k', ls='--', linewidth=1)
 
-    for y, evdp in zip(yLst, evdpLst):
+    for karr, y, evdp in plotLst:
         ax.plot(karr, y/maxnorm * 1e8, c=str(gray_mapping(evdp)), zorder=int(evdp*1e6))
 
-    ax.set_xlim(0, karr[-1])
+    ax.set_xlim(0, kmax * 0.4)
     ax.grid(linewidth=0.4)
 
 

docs/source/Advanced/integ_converg/run_ptam/run.py

@@ -12,7 +12,7 @@
 distarr = [5,8,10]
 # 设置 converg_method='PTAM' 进行收敛
 stgrnLst = pymod.compute_grn(
-    distarr=distarr, nt=500, dt=0.02, converg_method='PTAM',
+    distarr=distarr, nt=500, dt=0.02, converg_method='PTAM', k0=2, ampk=1.2, use_kmax_ref=True,
     statsfile=f"pygrtstats_{depsrc}_{deprcv}", statsidxs=[50,100]
 )
 # END DEPSRC 0.0 DGRN
@@ -46,7 +46,7 @@
 
 xarr = np.array([2.0])
 yarr = np.array([2.0])
-static_grn = pymod.compute_static_grn(xarr, yarr, converg_method='PTAM', statsfile=f"static_pygrtstats_{depsrc}_{deprcv}")
+static_grn = pymod.compute_static_grn(xarr, yarr, converg_method='PTAM', statsfile=f"static_pygrtstats_{depsrc}_{deprcv}", k0=3, use_kmax_ref=True)
 
 ir = 0
 statsdata1, statsdata2, ptamdata, dist = pygrt.utils.read_statsfile_ptam(f"static_pygrtstats_{depsrc}_{deprcv}/PTAM_{ir:04d}_*/PTAM")

docs/source/Advanced/k_integ/kmax.rst

@@ -10,22 +10,38 @@
 
    P_m(\omega) = \Delta k \sum_{j=0}^{\infty} F_m(k_j,\omega)J_m(k_j r)k_j
 
-其中 :math:`\Delta k = 2\pi/L, k_j=j\Delta k`，:math:`L` 为特征长度，即 :math:`L` 控制波数积分的积分间隔，默认根据 :doc:`/Tutorial/dynamic/gfunc` 部分介绍的约束条件自动确定。
+其中 :math:`\Delta k = 2\pi/L, k_j=j\Delta k`，:math:`L` 为特征长度，即 :math:`L` 控制波数积分的积分间隔，
+默认根据 :doc:`/Tutorial/dynamic/gfunc` 部分介绍的约束条件自动确定。
 
-程序会取一个较大值 :math:`k_{\text{max}}` 作为波数积分的上限，经验公式为
+程序首先根据经验公式确定波数积分搜索区间的上界 :math:`k_{\text{max,ref}}` 。对动态全波解，
 
 .. math:: 
 
-    k_{\text{max}} = \sqrt{k_0^2 + \left(s \cdot \dfrac{\omega}{v_{\text{min}}}\right)^2}
+    k_{\text{max,ref}} = \sqrt{\left(k_0 \cdot \dfrac{\pi}{h_s}\right)^2 + \left(s \cdot \dfrac{\omega}{v_{\text{min}}}\right)^2}
 
 其中 
 
-+ :math:`k_0` 为零频的波数积分上限，默认为 :math:`\dfrac{5\pi}{h_s}`， :math:`h_s` 为震源和场点的深度差（绝对值），限制最小为1km，若实际深度差小于该值，则自动使用直接收敛法（DCM）；
++ :math:`k_0` 为零频项的系数，默认为 50，程序内部使用 :math:`k_0 \cdot \pi / h_s` ；
+  :math:`h_s=\max(|z_s-z_r|, 0.1)` km 为震源和场点的深度差；
 + :math:`\omega` 为角频率；
-+ :math:`v_{\text{min}}` 为参考最小速度，默认取自模型中的最小速度，且限制在0.1km/s以上。
-+ :math:`s` 为放大系数，默认为1.15；
++ :math:`v_{\text{min}}` 为参考最小速度，默认取自模型中的最小速度，
+  且限制在 0.1 km/s 以上；
++ :math:`s` 为放大系数（``ampk``），默认为 2.0。
 
-程序支持自动判断积分收敛 (|yao1983|) 。当所有积分满足如下表达式时，自动退出波数循环（若达到 :math:`k_{\text{max}}` 则强制退出 ），
+对静态解，:math:`k_{\text{max,ref}} = k_0 \cdot \pi / h_s` 。
+
+默认情况下，程序在 :math:`[\Delta k, k_{\text{max,ref}}]` 内基于核函数振幅搜索
+实际积分上限 :math:`k_{\text{max}}` 。
+同深度时判断核函数是否逼近常数，异深度时判断振幅是否衰减至 0 。
+若指定 ``use_kmax_ref=True`` （C 模块 **+f**），则直接使用 :math:`k_{\text{max,ref}}` 作为
+:math:`k_{\text{max}}` 。
+
+若振幅搜索达到 :math:`k_{\text{max,ref}}` 仍未满足收敛判据，
+程序在默认（Auto）模式下将自动启用直接收敛法（DCM）处理积分收敛。
+当震源与场点完全同深度时，也会自动使用 DCM 。
+
+程序还支持提前判断积分收敛 (|yao1983|) 。
+当所有积分满足如下表达式时，自动退出波数循环（若达到 :math:`k_{\text{max}}` 则强制退出 ）。
 
 .. math:: 
 
@@ -46,13 +62,18 @@
 
     .. group-tab:: Python
 
-        :func:`compute_grn() <pygrt.pymod.PyModel1D.compute_grn>` 函数支持以下可选参数来控制波数积分，具体说明详见API。
+        :func:`compute_grn() <pygrt.pymod.PyModel1D.compute_grn>` 函数支持以下可选参数来控制波数积分，
+        具体说明详见API。
 
-        + ``k0:float``, 对应公式中的 :math:`k_0` 的系数，默认为5 
-        + ``ampk:float``, 对应公式中的 :math:`s` ，默认为1.15 
-        + ``keps:float`` 对应公式中的 :math:`\epsilon`，默认为-1（不使用）
+        + ``k0:float``, 对应公式中零频项的系数 :math:`k_0` ，默认为 50 
+        + ``ampk:float``, 对应公式中的 :math:`s` ，默认为 2.0 
+        + ``keps:float`` 对应公式中的 :math:`\epsilon`，默认为 -1（不使用）
+        + ``use_kmax_ref:bool`` 为 True 时直接使用 :math:`k_{\text{max,ref}}` 作为
+          积分上限
 
-        :func:`compute_static_grn() <pygrt.pymod.PyModel1D.compute_static_grn>` 函数支持以下可选参数来控制波数积分，参数与上面对应，具体说明详见API。
+        :func:`compute_static_grn() <pygrt.pymod.PyModel1D.compute_static_grn>` 函数支持以下可选参数来控制波数积分，
+        参数与上面对应，具体说明详见API。
 
         + ``k0:float`` 
         + ``keps:float`` 
+        + ``use_kmax_ref:bool``

docs/source/Advanced/kernel_old/run/run.py

@@ -12,7 +12,7 @@
 # 不指定statsidx，默认输出全部频率点的积分过程文件
 # vmin_ref 显式给定参考速度（用于定义波数积分上限），避免使用PTAM
 # Length 给定波数积分间隔dk
-_ = pymod.compute_grn(distarr=[1], nt=500, dt=0.02, vmin_ref=0.1, Length=20, k0_is_fixed=True, converg_method='none', statsfile="pygrtstats")
+_ = pymod.compute_grn(distarr=[1], nt=500, dt=0.02, vmin_ref=0.1, Length=20, use_kmax_ref=True, converg_method='none', statsfile="pygrtstats")
 # END GRN
 # -----------------------------------------------------------------
 

docs/source/Gallery/ex17/run.sh

@@ -13,8 +13,8 @@ cat > halfspace_Q <<EOF
 EOF
 
 # 强衰减介质需要更高的积分上限
-grt greenfn -Mhalfspace   -N2000/0.001+a -D0/0  -R5  -OGRN     -Cd
-grt greenfn -Mhalfspace_Q -N2000/0.001+a -D0/0  -R5  -OGRN_Q   -Cd
+grt greenfn -Mhalfspace   -N2000/0.001+a -D0/0  -R5  -OGRN  
+grt greenfn -Mhalfspace_Q -N2000/0.001+a -D0/0  -R5  -OGRN_Q
 python plot.py
 
 cp compare.svg cover.svg

docs/source/Module/explain_-Cconverg.rst_

@@ -1,7 +1,8 @@
 .. _-C:
 
 **-C**\ **d|p|n**
-    设置波数积分收敛方法。默认当 :math:`|z_s - z_r| <= 1.0` km 时，自动使用直接收敛法。
+    设置波数积分收敛方法。默认当震源与场点完全同深度时，或振幅搜索达到 :math:`k_{\text{max,ref}}` 仍未收敛时，
+    自动使用直接收敛法。
     支持以下子选项：
 
     + **d** - :doc:`/Advanced/integ_converg/dcm` 

docs/source/Module/explain_-D.rst_

@@ -1,5 +1,4 @@
 .. _-D:
 
 **-D**\ *depsrc/deprcv*
-    震源深度 *depsrc* (km) 和台站深度 *deprcv* (km)。
-    如果是在使用波数积分法求解格林函数，则当二者深度差小于 1 km 时，自动使用快速收敛算法。
+    震源深度 *depsrc* (km) 和台站深度 *deprcv* (km)。

docs/source/Module/greenfn.rst

@@ -22,7 +22,7 @@ greenfn
 [ |-H|\ *f1/f2* ]
 [ |-L|\ *length*\ [**+l**\ *Flength*][**+a**\ *Ftol*][**+o**\ *offset*] ]
 [ |-C|\ **d|p|n** ]
-[ |-K|\ [**+k**\ *k0*][**+s**\ *ampk*][**+e**\ *keps*][**+v**\ *vmin*] ]
+[ |-K|\ [**+k**\ *k0*][**+f**][**+s**\ *ampk*][**+e**\ *keps*][**+v**\ *vmin*] ]
 [ |-E|\ [**p**]\ *t0*\ [/*v0*] ]
 [ |-P|\ *nthreads* ]
 [ |-G|\ **e|v|h|s** ]
@@ -98,18 +98,26 @@ greenfn
 
 .. _-K:
 
-**-K**\ [**+k**\ *k0*][**+s**\ *ampk*][**+e**\ *keps*][**+v**\ *vmin*]
-    控制波数积分上限
+**-K**\ [**+k**\ *k0*][**+f**][**+s**\ *ampk*][**+e**\ *keps*][**+v**\ *vmin*]
+    控制波数积分搜索区间的上界 :math:`k_{\text{max,ref}}`
 
     .. math::
 
-        k_{\text{max}} = \sqrt{ k_0 \cdot \dfrac{\pi}{\Delta h} + \textit{<ampk>} \cdot \left(\dfrac{\omega}{v_{\text{min}}}\right)^2}
-
-    + **+k**\ *k0* - 控制零频的积分上限 [5.0]，其中深度差 :math:`\Delta h = \max(|z_s - z_r|, 1.0)` 。
-    + **+s**\ *ampk* - 放大倍数 [1.15] 。
-    + **+e**\ *keps* - 用于判断提前结束波数积分的收敛精度[0.0, 默认不使用]，详见
-      Yao and Harkrider (1983) 和 :doc:`/Advanced/k_integ/kmax` 。
-    + **+v**\ *vmin* - 参考最小速度，默认 :math:`\max{(\min\limits_{i} (\alpha_i \cup \beta_i), 0.1)}` 。
+        k_{\text{max,ref}} = \sqrt{ \left(k_0 \cdot \dfrac{\pi}{\Delta h}\right)^2 + \textit{<ampk>}^2 \cdot \left(\dfrac{\omega}{v_{\text{min}}}\right)^2 }
+
+    程序在 :math:`[\Delta k, k_{\text{max,ref}}]` 内基于核函数振幅搜索实际积分上限 :math:`k_{\text{max}}` 。
+    若搜索达到 :math:`k_{\text{max,ref}}` 仍未收敛，或震源与场点完全同深度时，
+    默认模式下将自动启用 DCM 。
+
+    + **+k**\ *k0* - 零频项系数 [50.0]，
+      其中深度差 :math:`\Delta h = \max(|z_s - z_r|, 0.1)` 。
+    + **+f** - 直接使用 :math:`k_{\text{max,ref}}` 作为积分上限，
+      不进行振幅搜索。
+    + **+s**\ *ampk* - 放大倍数 [2.0] 。
+    + **+e**\ *keps* - 用于判断提前结束波数积分的收敛精度[0.0, 默认不使用]，
+      详见 Yao and Harkrider (1983) 和 :doc:`/Advanced/k_integ/kmax` 。
+    + **+v**\ *vmin* - 参考最小速度，
+      默认 :math:`\max{(\min\limits_{i} (\alpha_i \cup \beta_i), 0.1)}` 。
 
 .. include:: explain_-Cconverg.rst_
 

docs/source/Module/static_greenfn.rst

@@ -22,7 +22,7 @@ static_greenfn
 [ |-B|\ **f|F|r|R|h|H** ]
 [ |-L|\ *length*\ [**+l**\ *Flength*][**+a**\ *Ftol*][**+o**\ *offset*] ]
 [ |-C|\ **d|p|n** ]
-[ |-K|\ [**+k**\ *k0*][**+e**\ *keps*] ]
+[ |-K|\ [**+k**\ *k0*][**+f**][**+e**\ *keps*] ]
 [ |-S| ]
 [ **-e** ]
 [ **-h** ]
@@ -66,12 +66,19 @@ static_greenfn
 
 .. _-K:
 
-**-K**\ [**+k**\ *k0*][**+e**\ *keps*]
-    控制波数积分上限 :math:`k_0 \cdot \dfrac{\pi}{\Delta h}`
+**-K**\ [**+k**\ *k0*][**+f**][**+e**\ *keps*]
+    控制波数积分搜索区间的上界 :math:`k_{\text{max,ref}} = k_0 \cdot \dfrac{\pi}{\Delta h}`
 
-    + **+k**\ *k0* - 控制零频的积分上限 [5.0]，其中深度差 :math:`\Delta h = \max(|z_s - z_r|, 1.0)` 。
-    + **+e**\ *keps* - 用于判断提前结束波数积分的收敛精度[0.0, 默认不使用]，详见
-      Yao and Harkrider (1983) 和 :doc:`/Advanced/k_integ/kmax` 。
+    程序在 :math:`[\Delta k, k_{\text{max,ref}}]` 内基于核函数振幅搜索实际积分上限 :math:`k_{\text{max}}` 。
+    若搜索达到 :math:`k_{\text{max,ref}}` 仍未收敛，或震源与场点完全同深度时，
+    默认模式下将自动启用 DCM 。
+
+    + **+k**\ *k0* - 零频项系数 [50.0]，
+      其中深度差 :math:`\Delta h = \max(|z_s - z_r|, 0.1)` 。
+    + **+f** - 直接使用 :math:`k_{\text{max,ref}}` 作为积分上限，
+      不进行振幅搜索。
+    + **+e**\ *keps* - 用于判断提前结束波数积分的收敛精度[0.0, 默认不使用]，
+      详见 Yao and Harkrider (1983) 和 :doc:`/Advanced/k_integ/kmax` 。
 
 .. include:: explain_-Cconverg.rst_
 

pygrt/C_extension/include/grt/common/const.h

@@ -55,7 +55,7 @@ typedef double complex cplx_t;
 #define DEG1 0.017453292519943295  ///< \f$ \frac{\pi}{180} \f$ 
 #define GOLDEN_RATIO 0.6180339887498949  ///< \f$ \frac{\sqrt{5}-1}{2} \f$
 
-#define GRT_MIN_DEPTH_GAP_SRC_RCV  1.0  ///< 震源和台站的最小深度差（不做绝对限制，仅用于参考波数积分上限，以及判断是否需要其它收敛方法）
+#define GRT_MIN_DEPTH_GAP_SRC_RCV  0.1  ///< 震源和台站的最小深度差（不做绝对限制，仅用于 kmax_ref 中 hs 的下限）
 #define GCC_ALWAYS_INLINE __attribute__((always_inline))  ///< gcc编译器不改动内联函数
 
 #define GRT_SWAP(type, a, b) { type temp = a; a = b; b = temp; } ///< 交换两个变量的值

pygrt/C_extension/include/grt/integral/integ_process.h

@@ -34,15 +34,19 @@ typedef enum {
 
 // 描述不同波数积分方法的结构体
 typedef struct {
-    real_t k0;      ///< 波数积分的上限 \f$ \tilde{k_{max}}=\sqrt{(k_{0}*\pi/hs)^2 + (ampk*w/vmin_{ref})^2} \f$ ，k循环必须退出, hs=max(震源和台站深度差,1.0) 
-    bool k0_is_fixed;  ///< 固定 k0，默认在程序中自动调整 k0
-    real_t ampk;    ///< 影响波数k积分上限的系数
-    real_t keps;    ///< 波数积分的收敛条件，要求在某震中距下所有格林函数都收敛，为负数代表不提前判断收敛，按照波数积分上限进行积分 
-    real_t vmin;    ///< 参考最小速度，用于定义波数积分的上限
+    real_t k0;      ///< 用户参数 k0 经 \f$ \pi/hs \f$ 缩放后的零频项，
+                    ///< 与 \f$ ampk*\omega/vmin \f$ 共同确定搜索上界 kmax_ref；
+                    ///< hs=max(震源和台站深度差, 0.1)
+    bool use_kmax_ref;  ///< 为 true 时直接将 kmax_ref 作为积分上限，不进行基于振幅的搜索
+    real_t ampk;    ///< 影响 kmax_ref 中频率相关项的系数，默认 2.0
+    real_t keps;    ///< 波数积分的收敛条件，要求在某震中距下所有格林函数都收敛；
+                    ///< 为负数代表不提前判断收敛
+    real_t vmin;    ///< 参考最小速度，用于定义 kmax_ref
 
     real_t kcut;    ///< 波数积分和Filon积分的分割点
 
-    real_t kmax;    ///< 全局波数最大值，程序运行中会随频率变动
+    real_t kmax;    ///< 实际波数积分上限，默认在 [dk, kmax_ref] 内由振幅搜索确定，
+                    ///< 程序运行中会随频率变动
 
     real_t dk;      ///< DWM 的波数积分间隔
 

pygrt/C_extension/include/grt/integral/kmax.h

@@ -18,14 +18,14 @@
  *
  * 搜索在 log(k) 上等间距推进（自适应几何因子），目标约 40 步覆盖 [kmax_init, kmax_ref]
  * 在 kmax_low 之前不判断收敛
- * 同深度使用当前核函数与前一采样点的差值 dF 相对于搜索过程中的最大振幅 Fmax 判断逼近常数；
- * 异深度用振幅相对峰值衰减判断逼近 0
+ * 同深度使用当前核函数与前一采样点的差值 dF 相对于搜索过程中的最大振幅 Fmax
+ * 判断逼近常数；异深度用振幅相对峰值衰减判断逼近 0
  * 
  * @param[in,out]  mstat        已设置频率的模型状态
  * @param[in]      kerfunc      待检查的核函数
  * @param[in]      kmax_init    扫描的初始波数
  * @param[in]      kmax_low     允许判断收敛的最低波数
- * @param[in]      kmax_ref     最大上限
+ * @param[in]      kmax_ref     搜索区间的上界（经验公式）
  * @param[out]     Ncount       估计过程中计算核函数的次数，可为 NULL
  * @return     估计的 kmax
  */

pygrt/C_extension/src/dynamic/grn.c

@@ -150,9 +150,9 @@ void grt_integ_grn_spec(MODEL1D *mod1d, K_INTEG_PROCESS *Kproc, GRNSPEC *grn, co
 
         // ===================================================================================
         //                          Wavenumber Integration
-        // 每个频率均直接根据动态核函数估计积分上限
+        // 每个频率根据 kmax_ref 与核函数振幅搜索实际积分上限
         real_t kmax_ref = hypot(local_Kproc->k0, local_Kproc->ampk * w / local_Kproc->vmin);
-        if(local_Kproc->k0_is_fixed){
+        if(local_Kproc->use_kmax_ref){
             local_Kproc->kmax = kmax_ref;
             size_t nk = floor(local_Kproc->kmax / local_Kproc->dk) + 1;
             #pragma omp critical(grn_console)