FEAT: support set new XY grid for static results

docs/source/Tutorial/static/run/run.py

@@ -9,7 +9,10 @@
 
 xarr = np.linspace(-3, 3, 41)
 yarr = np.linspace(-2.5, 2.5, 33)
-static_grn = pymod.compute_static_grn(xarr, yarr)
+# 可以设置 distarr 来指定震中距序列
+# static_grn = pymod.compute_static_grn(distarr=np.arange(0,10+1e-8,0.1))
+# 也可以设置 xarr 和 yarr 来指定 XY 网格
+static_grn = pymod.compute_static_grn(xarr=xarr, yarr=yarr)
 print(static_grn.keys())
 # dict_keys(['_xarr', '_yarr', '_src_va', '_src_vb', '_src_rho', '_rcv_va', '_rcv_vb', '_rcv_rho', 'EXZ', 'VFZ', 'DDZ', 'HFZ', 'DSZ', 'SSZ', 'EXR', 'VFR', 'DDR', 'HFR', 'DSR', 'SSR', 'HFT', 'DST', 'SST'])
 # END GRN
@@ -114,3 +117,15 @@ def plot_static(static_syn:dict, out:Union[str,None]=None):
 plot_static(static_syn, "syn_mt2.svg")
 # END SYN MT2
 # ---------------------------------------------------------------------------------
+
+
+# ---------------------------------------------------------------------------------
+# BEGIN NEW XY
+xarr2 = np.arange(-3, 3+1e-8, 0.2)
+yarr2 = np.arange(-2.5, 2.5+1e-8, 0.25)
+static_syn = pygrt.utils.gen_syn_from_gf_DC(static_grn, M0=1e24, strike=33, dip=90, rake=0, ZNE=True, xarr=xarr2, yarr=yarr2)
+print(static_syn.keys())
+# dict_keys(['_xarr', '_yarr', '_src_va', '_src_vb', '_src_rho', '_rcv_va', '_rcv_vb', '_rcv_rho', 'Z', 'N', 'E'])
+plot_static(static_syn, "synXY_dc2.svg")
+# END NEW XY
+# ---------------------------------------------------------------------------------

docs/source/Tutorial/static/run/run.sh

@@ -41,6 +41,9 @@ y1=-2.5
 y2=2.5
 dy=0.15
 # 输出到NetCDF网格文件
+# 可以使用 -R 来指定震中距序列
+# grt static greenfn -Mmilrow -D${depsrc}/${deprcv} -R0/10/0.1 -Ostgrn.nc
+# 也可以使用 -X 和 -Y 来指定 XY 网格
 grt static greenfn -Mmilrow -D${depsrc}/${deprcv} -X$x1/$x2/$dx -Y$y1/$y2/$dy -Ostgrn.nc
 # END GRN
 # ---------------------------------------------------------------------------------
@@ -160,7 +163,6 @@ grt static syn -S1e24 -T0/-0.2/0/0/0/0 -N -Gstgrn.nc -Ostsyn_mt2.nc
 # END SYN MT2
 # ---------------------------------------------------------------------------------
 
-
 # X Y depth mrr   mtt   mff   mrt   mrf   mtf   exp
 #           mzz   mxx   myy   mxz  -myz  -mxy
 gmt begin syn_mt2 pdf
@@ -171,6 +173,23 @@ EOF
     gmt colorbar -Bx+l"Z (cm)"
 gmt end
 
+
+# ---------------------------------------------------------------------------------
+# BEGIN NEW XY
+# 从网格文件中读取格林函数，再将合成结果写入新网格
+grt static syn -S1e24 -M33/90/0 -N -Gstgrn.nc -OstsynXY_dc2.nc -X-3/3/0.2 -Y-2.5/2.5/0.25
+# END NEW XY
+# ---------------------------------------------------------------------------------
+
+gmt begin synXY_dc2 pdf
+    gmtplot_static stsynXY_dc2.nc -Si0.03c
+    gmt meca -Sa0.5c <<EOF
+0 0 2 33 90 0 5
+EOF
+    gmt colorbar -Bx+l"Z (cm)"
+gmt end
+
+
 # 统一转为 svg
 for pdfname in $(ls *.pdf); do
     name=$(basename $pdfname .pdf)

docs/source/Tutorial/static/static_gfunc.rst

@@ -8,11 +8,20 @@ Python中计算静态格林函数的主函数为 :func:`compute_static_grn() <py
 
 建议先阅读完 :doc:`/Tutorial/dynamic/gfunc` 部分。静态情况与动态情况采取的计算方法一致，只是推导细节会有不同，详见 |yao2026| 。
 
+静态解模块对于传入“震中距”的方式及后续处理与动态解模块相比稍有不同。
+不论是动态解还是静态解，在模型、源点深度和场点深度确定的情况下，格林函数仅与 **和震中距** 相关。
+但考虑到静态解的普遍应用场景，程序在计算静态格林函数阶段有两种传入震中距的方式：
+
++ 指定二维的 XY 网格。坐标 (0,0) 为源点的水平投影，每个节点的震中距为 :math:`r_{ij} = \sqrt{x_i^2 + y_j^2}` 。实际计算中会对震中距自动去重以减少计算量。
++ 指定一维的震中距序列。在计算和结果保存上等同于指定从原点 (X=0.0) 沿东向 (Y/km) 指定对应的采样点。
+
+为了兼容性，结果的输出仍然保持以 XY 网格的形式。在 :doc:`static_syn` 阶段，可以指定新的 XY 网格，
+此时每个节点的格林函数会近似为最近震中距的格林函数。
 
 示例程序
 -----------
 
-假设在 :file:`milrow` 模型中，震源深度2km，接收点位于地表。北向(X/km)在[-3,3]范围内等距间隔 0.15 km 采样，东向(Y/km)在[-2.5,2.5]范围内等距间隔 0.15 km 采样，计算这些点上的静态格林函数。
+假设在 :file:`milrow` 模型中，震源深度2km，接收点位于地表。
 
 .. tabs::  
 

docs/source/Tutorial/static/static_syn.rst

@@ -220,4 +220,48 @@ Python中合成静态位移的主函数为 :func:`gen_syn_from_gf_*() <pygrt.uti
 
 .. figure:: run/syn_mt2.svg
     :width: 500px
-    :align: center
+    :align: center
+
+
+指定新的 XY 网格
+---------------------
+
+以上示例结果中， XY 网格位置都沿用了静态格林函数计算时传入的 XY 网格。
+
+程序也支持在合成阶段指定新的 XY 网格，此时每个节点将使用最近震中距的格林函数来进行合成，
+因此这是一种近似方法，程序会打印震中距之差的统计信息，以方便评估近似效果。
+但近似好处是一旦静态格林函数计算好，合成阶段使用新的网格也可以复用格林函数，无需再计算。
+
+以下以一个与上面相同的走滑断层作为示例进行计算，选取了间隔稍大的网格。
+
+.. tabs::  
+
+    .. group-tab:: C 
+
+        .. literalinclude:: run/run.sh
+            :language: bash
+            :start-after: BEGIN NEW XY
+            :end-before: END NEW XY
+
+    .. group-tab:: Python 
+
+        .. literalinclude:: run/run.py
+            :language: python
+            :start-after: BEGIN NEW XY
+            :end-before: END NEW XY
+
+.. grid:: 2
+
+    .. grid-item::
+
+        .. figure:: run/syn_dc2.svg
+            :width: 400px
+            :align: center
+
+    .. grid-item::
+
+        .. figure:: run/synXY_dc2.svg
+            :width: 400px
+            :align: center
+
+

pygrt/C_extension/src/static/grt_static_greenfn.c

@@ -149,9 +149,14 @@ printf("\n"
 "\n\n"
 "Usage:\n"
 "----------------------------------------------------------------\n"
-"    grt static greenfn -M<model> -D<depsrc>/<deprcv>  -X<x1>/<x2>/<dx> \n"
-"          -Y<y1>/<y2>/<dy>  -O<outgrid>  [-L<length>] [-C[d|p|n]] \n" 
-"           [-Bf|F|r|R|h|H] [-K[+k<k0>][+e<keps>]] [-S]  [-e]\n"
+"    grt static greenfn -M<model> -D<depsrc>/<deprcv>  -O<outgrid>  \n"
+"           [-X<x1>/<x2>/<dx>] [-Y<y1>/<y2>/<dy>] [-R<r1>/<r2>/<dr>]\n" 
+"           [-L<length>]   [-C[d|p|n]]  [-Bf|F|r|R|h|H] \n"
+"           [-K[+k<k0>][+e<keps>]] [-S]  [-e]\n"
+"\n"
+"    There're two ways to define the \"epicentral distances\":\n"
+"    1. set both -X and -Y to define a XY grid in advance.\n"
+"    2. simply set -R, which equal to \"-X0/0/1 -Y<r1>/<r2>/<nr>\".\n"
 "\n\n"
 "Options:\n"
 "----------------------------------------------------------------\n"
@@ -184,6 +189,9 @@ printf("\n"
 "                 <y2>: end coordinate (km).\n"
 "                 <dy>: sampling interval (km).\n"
 "\n"
+"    -R<r1>/<r2>/<dr>\n"
+"                 equal to \"-X0/0/1 -Y<r1>/<r2>/<nr>\". \n"
+"\n"
 "    -O<outgrid>  Filepath to output nc grid.\n"
 "\n"
 "    -L[<length>][+l<Flength>][+a<Ftol>][+o<offset>]\n"
@@ -241,7 +249,8 @@ printf("\n"
 "\n\n"
 "Examples:\n"
 "----------------------------------------------------------------\n"
-"    grt static greenfn -Mmilrow -D2/0 -X-10/10/20 -Y-10/10/20 -Ostgrn.nc\n"
+"    grt static greenfn -Mmilrow -D2/0 -X-10/10/1 -Y-10/10/1 -Ostgrn.nc\n"
+"    grt static greenfn -Mmilrow -D2/0 -R0/20/1 -Ostgrn.nc\n"
 "\n\n\n"
 );
 }
@@ -260,7 +269,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){
 
     int opt;
 
-    while ((opt = getopt(argc, argv, ":M:D:B:L:C:K:X:Y:O:Seh")) != -1) {
+    while ((opt = getopt(argc, argv, ":M:D:B:L:C:K:X:Y:R:O:Seh")) != -1) {
         switch (opt) {
             // 模型路径，其中每行分别为 
             //      厚度(km)  Vp(km/s)  Vs(km/s)  Rho(g/cm^3)  Qp   Qs
@@ -466,6 +475,37 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){
                 }
                 break;
 
+            // -R 算是别名，相当于 -X0/0/1, -Yy1/y2/dy, 以此方式指定一维震中距序列
+            case 'R':
+                Ctrl->X.active = Ctrl->Y.active = true;
+                {
+                    real_t a1, a2, delta;
+                    if(3 != sscanf(optarg, "%lf/%lf/%lf", &a1, &a2, &delta)){
+                        GRTBadOptionError(R, "");
+                    };
+                    if(delta <= 0){
+                        GRTBadOptionError(R, "Can't set nonpositive dr(%f)", delta);
+                    }
+                    // 特别检查 ys 不能是负数
+                    if(a1 < 0.0){
+                        GRTBadOptionError(R, "Can't set nonpositive r1(%f)", a1);
+                    }
+                    if(a1 > a2){
+                        GRTBadOptionError(R, "r1(%f) > r2(%f).", a1, a2);
+                    }
+
+                    Ctrl->Y.ny = floor((a2-a1)/delta) + 1;
+                    Ctrl->Y.ys = (real_t*)calloc(Ctrl->Y.ny, sizeof(real_t));
+                    for(size_t i=0; i<Ctrl->Y.ny; ++i){
+                        Ctrl->Y.ys[i] = a1 + delta*i;
+                    }
+
+                    Ctrl->X.nx = 1;
+                    Ctrl->X.xs = (real_t*)calloc(Ctrl->X.nx, sizeof(real_t));
+                    Ctrl->X.xs[0] = 0.0;
+                }
+                break;
+
             // 输出 nc 文件名
             case 'O':
                 Ctrl->O.active = true;