FEAT: multi-depth static Green's functions with finite-fault and arbitrary-receiver synthesis

docs/source/Gallery/ex15/plot.py

@@ -50,8 +50,8 @@
 ax.plot(distarr, dynamic_Z, 'k', label='Dynamic Z')
 ax.plot(distarr, dynamic_R, 'k', label='Dynamic R')
 
-ax.plot(distarr, static_grn['variables']['EXZ']['data'][0] * coef, 'ro', ms=4, label='Static Z')
-ax.plot(distarr, static_grn['variables']['EXR']['data'][0] * coef, 'bo', ms=4, label='Static R')
+ax.plot(distarr, static_grn['variables']['EXZ']['data'][0, 0, 0] * coef, 'ro', ms=4, label='Static Z')
+ax.plot(distarr, static_grn['variables']['EXR']['data'][0, 0, 0] * coef, 'bo', ms=4, label='Static R')
 
 ax.set_xlim(0, 50)
 ax.set_xlabel('Distance (km)')

docs/source/Gallery/ex15/plot_all.py

@@ -62,12 +62,12 @@
     ms = 2
     lw = 0.8
     ax.plot(distarr, dynamic_Z, 'k', lw=lw, label='Dynamic Z')
-    ax.plot(distarr, static_grn['variables'][f'{src}Z']['data'][0] * coef, 'ro', ms=ms, label='Static Z')
+    ax.plot(distarr, static_grn['variables'][f'{src}Z']['data'][0, 0, 0] * coef, 'ro', ms=ms, label='Static Z')
     ax.plot(distarr, dynamic_R, 'k', lw=lw, label='Dynamic R')
-    ax.plot(distarr, static_grn['variables'][f'{src}R']['data'][0] * coef, 'bo', ms=ms, label='Static R')
+    ax.plot(distarr, static_grn['variables'][f'{src}R']['data'][0, 0, 0] * coef, 'bo', ms=ms, label='Static R')
     if src not in ['EX', 'VF', 'DD']:
         ax.plot(distarr, dynamic_T, 'k', lw=lw, label='Dynamic T')
-        ax.plot(distarr, static_grn['variables'][f'{src}T']['data'][0] * coef, 'go', ms=ms, label='Static T')
+        ax.plot(distarr, static_grn['variables'][f'{src}T']['data'][0, 0, 0] * coef, 'go', ms=ms, label='Static T')
 
     ax.set_xlim(0, 50)
     ax.set_xlabel('Distance (km)')

docs/source/Gallery/ex16/plot.py

@@ -115,8 +115,8 @@
     ax.set_ymargin(0.3)
 
     ax = axs2[i]
-    ax.plot(easts, static1['variables'][chLst[i]]['data'][0], **prop1)
-    ax.plot(easts, static2['variables'][chLst[i]]['data'][0] * sgn, **prop2)
+    ax.plot(easts, static1['variables'][chLst[i]]['data'][0, 0, 0], **prop1)
+    ax.plot(easts, static2['variables'][chLst[i]]['data'][0, 0, 0] * sgn, **prop2)
     ax.text(0.96, 0.9, chLst[i], transform=ax.transAxes, 
             ha='right', va='top', bbox=dict(fc='w'))
     ax.ticklabel_format(axis='y', style='sci', scilimits=(0,0))

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

@@ -24,7 +24,7 @@
 print(static_grn.keys())
 # dict_keys(['dimensions', 'variables', 'attributes'])
 print(list(static_grn["variables"].keys()))
-# ['north', 'east', 'EXZ', 'EXR', ...]
+# ['model', 'depsrc', 'deprcv', 'north', 'east', 'src_va', ..., 'EXZ', 'EXR', ...]
 # END GRN
 # ---------------------------------------------------------------------------------
 

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

@@ -118,7 +118,7 @@ grt static syn -S1e24 -M33/50 -N -Gstgrn.nc -Ostsyn_ts.nc
 gmt begin syn_ts pdf
     gmtplot_static stsyn_ts.nc -Si0.03c
     gmt meca -Sz0.5c <<EOF
-0 0 2 $(python tension2mt.py 33 50 stsyn_ts.nc) 24
+0 0 2 $(python tension2mt.py 33 50 stgrn.nc) 24
 EOF
     gmt colorbar -Bx+l"Z (cm)"
 gmt end
@@ -134,7 +134,7 @@ grt static syn -S1e24 -M33/90 -N -Gstgrn.nc -Ostsyn_ts2.nc
 gmt begin syn_ts2 pdf
     gmtplot_static stsyn_ts2.nc -Si0.03c
     gmt meca -Sz0.5c <<EOF
-0 0 2 $(python tension2mt.py 33 90 stsyn_ts2.nc) 24
+0 0 2 $(python tension2mt.py 33 90 stgrn.nc) 24
 EOF
     gmt colorbar -Bx+l"Z (cm)"
 gmt end

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

@@ -20,10 +20,10 @@
 
 vvec = nvec.copy()
 
-# 从 nc 文件中读取 src_va, src_vb
+# 单震源深度：直接取该层的 src_va / src_vb
 with netcdf_file(ncfile, mmap=False) as f:
-    src_va = f._attributes['src_va']
-    src_vb = f._attributes['src_vb']
+    src_va = float(f.variables["src_va"][0])
+    src_vb = float(f.variables["src_vb"][0])
 
 M = ((src_va/src_vb)**2 - 2) * np.eye(3) * np.sum(nvec * vvec) + (np.einsum('i,j', nvec, vvec) + np.einsum('i,j', vvec, nvec))
 

pygrt/C_extension/include/grt.h

@@ -26,6 +26,7 @@
 #include "grt/common/colorstr.h"
 #include "grt/common/const.h"
 #include "grt/common/coord.h"
+#include "grt/common/finite_fault.h"
 #include "grt/common/logo.h"
 #include "grt/common/matrix.h"
 #include "grt/common/model.h"
@@ -56,7 +57,9 @@
 #include "grt/modal/secular.h"
 
 
+#include "grt/static/recv_points.h"
 #include "grt/static/static_grn.h"
+#include "grt/static/stgrnlib.h"
 #include "grt/static/static_layer.h"
 #include "grt/static/static_source.h"
 

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

@@ -122,6 +122,7 @@ typedef double complex cplx_t;
 #define GRT_MORDER_MAX   2    ///< 2, 代码中阶数m的最大值
 #define GRT_SRC_M_NUM    6   ///< 6, 代码中不同震源、不同阶数的个数
 #define GRT_MECHANISM_NUM   6   ///<  6, 描述震源机制的最多参数
+#define GRT_MODARR_NCOL  6   ///< 6, 模型矩阵列数：Thk, Va, Vb, Rho, Qa, Qb
 
 #define GRT_PTAM_PT_MAX   36         ///< 36， 最后统计波峰波谷的目标数量
 #define GRT_PTAM_WINDOW_SIZE  3      ///< 3，  使用连续点数判断是否为波峰或波谷

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

@@ -0,0 +1,105 @@
+/**
+ * @file   finite_fault.h
+ * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn)
+ * @date   2026-08
+ *
+ * Coulomb 格式有限断层：读入、衍生量与几何剖分
+ *
+ */
+
+#pragma once
+
+#include "grt/common/const.h"
+
+/** Coulomb 程序格式的有限断层（及衍生量） */
+typedef struct {
+    real_t east_begin;
+    real_t north_begin;
+    real_t east_end;
+    real_t north_end;
+    real_t right_lateral;   ///< m
+    real_t reverse;         ///< m
+    real_t dip;             ///< degree
+    real_t top;             ///< km
+    real_t bot;             ///< km
+
+    // 衍生量，由 grt_finite_fault_set_derived 填充
+    real_t strike;          ///< degree
+    real_t rake;            ///< degree
+    real_t slip;            ///< m, hypot(right_lateral, reverse)
+} FINITE_FAULT;
+
+/** 剖分后的单个子断层（点源几何） */
+typedef struct {
+    real_t east;      ///< 中心 east (km)
+    real_t north;     ///< 中心 north (km)
+    real_t depsrc;    ///< 中心深度 (km)
+    real_t width;     ///< 沿倾向边长 (km)
+    real_t length;    ///< 沿走向边长 (km)
+    real_t potency;   ///< 矩势 (cm^3) = slip(m)*width*length*1e12
+} FINITE_SUBFAULT;
+
+/**
+ * 由 Coulomb 原始字段计算 strike / rake / slip
+ *
+ * @param[in,out]  f   有限断层结构体
+ */
+void grt_finite_fault_set_derived(FINITE_FAULT *f);
+
+/**
+ * 读取 Coulomb 格式有限断层文件
+ *
+ * 跳过前两行表头，逐行解析 11 列中的断层字段，并对每条调用
+ * grt_finite_fault_set_derived。要求 dip ∈ (0, 90] 且 bot > top
+ * 调用方负责 grt_finite_fault_free
+ *
+ * @param[in]   path     文件路径
+ * @param[out]  nfault   读入的断层段数
+ * @return      新分配的 FINITE_FAULT 数组，失败则报错退出
+ */
+FINITE_FAULT *grt_finite_fault_load_coulomb(const char *path, size_t *nfault);
+
+/**
+ * 释放 grt_finite_fault_load_coulomb 返回的数组
+ *
+ * @param[in,out]  faults   断层数组，可为 NULL
+ */
+void grt_finite_fault_free(FINITE_FAULT *faults);
+
+/**
+ * 由断层几何与 dL/dW 得到倾向/走向长度及剖分数
+ *
+ * 要求 dip ∈ (0, 90]、bot > top，且 dL/dW > 0
+ *
+ * @param[in]   fault  有限断层（需已有 dip / 端点坐标）
+ * @param[in]   dL     沿走向剖分间隔 (km)
+ * @param[in]   dW     沿倾向剖分间隔 (km)
+ * @param[out]  W      沿倾向长度 (km)
+ * @param[out]  L      沿走向长度 (km)
+ * @param[out]  nW     倾向方向子断层数
+ * @param[out]  nL     走向方向子断层数
+ */
+void grt_finite_fault_subdiv(
+    const FINITE_FAULT *fault, real_t dL, real_t dW,
+    real_t *W, real_t *L, size_t *nW, size_t *nL);
+
+/**
+ * 计算 (iW, iL) 子断层中心几何与矩势
+ *
+ * 要求 fault 已 set_derived，且 W/L/nW/nL 由 grt_finite_fault_subdiv 给出
+ * 末块可短于 dL/dW，中心取该块中点：i*d + size/2
+ *
+ * @param[in]   fault  有限断层
+ * @param[in]   dL     沿走向剖分间隔 (km)
+ * @param[in]   dW     沿倾向剖分间隔 (km)
+ * @param[in]   W      沿倾向总长 (km)
+ * @param[in]   L      沿走向总长 (km)
+ * @param[in]   iW     倾向方向子断层索引 [0, nW)
+ * @param[in]   iL     走向方向子断层索引 [0, nL)
+ * @param[out]  sub    子断层几何
+ */
+void grt_finite_fault_subfault(
+    const FINITE_FAULT *fault,
+    real_t dL, real_t dW, real_t W, real_t L,
+    size_t iW, size_t iL,
+    FINITE_SUBFAULT *sub);

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

@@ -128,6 +128,32 @@ void grt_free_mod1d(MODEL1D *mod1d);
  */
 MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, real_t deprcv, bool allowLiquid);
 
+/**
+ * 从模型文件读取 ``nlayer × 6`` 矩阵（Thk, Va, Vb, Rho, Qa, Qb）
+ * 若首列以层顶深度给出则先转为厚度；末层厚度保持文件原值
+ * 
+ * @param[in]     modelpath      模型文件路径
+ * @param[out]    nlayer         层数
+ * @param[in]     allowLiquid    是否允许液体层（Vs==0）
+ * @return        新分配的矩阵，调用方负责 free
+ */
+real_t (* grt_read_modarr_from_file(
+    const char *modelpath, size_t *nlayer, bool allowLiquid))[GRT_MODARR_NCOL];
+
+/**
+ * 由 ``nlayer × 6`` 模型矩阵按深度查层介质
+ * 末层视为半空间；恰落在层界面时取上层
+ * 
+ * @param[in]     nlayer    层数
+ * @param[in]     modarr    模型矩阵，每行 Thk/Va/Vb/Rho/Qa/Qb
+ * @param[in]     depth     深度 (km)，须 >= 0
+ * @param[out]    va        P 波速 (km/s)，可为 NULL
+ * @param[out]    vb        S 波速 (km/s)，可为 NULL
+ * @param[out]    rho       密度 (g/cm^3)，可为 NULL
+ */
+void grt_modarr_medium_at_depth(
+    size_t nlayer, const real_t (*modarr)[GRT_MODARR_NCOL],
+    real_t depth, real_t *va, real_t *vb, real_t *rho);
 
 /**
  * 设置模型的边界条件，并对底界面做检查

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

@@ -157,3 +157,21 @@ typedef int (*grt_compare_fn)(const void *a, const void *b);
 int grt_argsort(
     const void *base, size_t n, size_t element_size,
     grt_compare_fn compare, size_t *indices);
+
+/**
+ * 在升序数组上定位 1D 线性插值的左右邻点与权重
+ *
+ * n==1 时要求 |q-x[0]| 在容差内，此时 i0=i1=0、w=0；
+ * q 越界（超出端点容差）返回 false
+ *
+ * @param[in]   x    升序坐标数组
+ * @param[in]   n    长度
+ * @param[in]   q    查询值
+ * @param[out]  i0   左邻索引
+ * @param[out]  i1   右邻索引
+ * @param[out]  w    权重，使值 = (1-w)*x[i0] + w*x[i1]
+ * @return      是否在范围内
+ */
+bool grt_locateLinearInterp(
+    const real_t *x, size_t n, real_t q,
+    size_t *i0, size_t *i1, real_t *w);

pygrt/C_extension/include/grt/static/recv_points.h

@@ -0,0 +1,79 @@
+/**
+ * @file   recv_points.h
+ * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn)
+ * @date   2026-08
+ *
+ * 静态 syn / 后处理用的接收点列表
+ * 网格 (-X/-Y 或延用库坐标) 与任意点文件 (-Q) 均展开为点列
+ *
+ */
+
+#pragma once
+
+#include <stdbool.h>
+#include <stddef.h>
+
+#include "grt/common/const.h"
+
+/** layout 字符串：写入 nc 全局属性，读端据此分支 */
+#define GRT_RECV_LAYOUT_GRID   "grid"
+#define GRT_RECV_LAYOUT_POINTS "points"
+
+/**
+ * 接收点坐标列表（三个分量各自连续存放）
+ *
+ * norths/easts/depths 长度均为 npts，单位 km
+ * is_grid 为真时可由 nnorth/neast 还原二维网格，ipt = ieast + inorth*neast
+ */
+typedef struct {
+    size_t npts;
+    real_t *norths;
+    real_t *easts;
+    real_t *depths;
+
+    bool is_grid;
+    size_t nnorth;
+    size_t neast;
+} GRT_RECV_POINTS;
+
+/**
+ * 由 north/east 轴与单一深度展开为点列（is_grid=true）
+ *
+ * @param[in]   nnorth   north 方向点数
+ * @param[in]   norths   north 坐标 (km)
+ * @param[in]   neast    east 方向点数
+ * @param[in]   easts    east 坐标 (km)
+ * @param[in]   depth    接收深度 (km)
+ * @return      新分配的 GRT_RECV_POINTS*，调用方负责 grt_recv_points_free
+ */
+GRT_RECV_POINTS *grt_recv_points_from_grid(
+    size_t nnorth, const real_t *norths,
+    size_t neast,  const real_t *easts,
+    real_t depth);
+
+/**
+ * 从 ASCII 文件读任意接收点（is_grid=false）
+ *
+ * 每行 north east depth (km)，# 开头为注释
+ *
+ * @param[in]   path   文件路径
+ * @return      新分配的 GRT_RECV_POINTS*
+ */
+GRT_RECV_POINTS *grt_recv_points_from_file(const char *path);
+
+/**
+ * 释放 GRT_RECV_POINTS（含 norths/easts/depths）
+ *
+ * @param[in,out]  pts   可为 NULL
+ */
+void grt_recv_points_free(GRT_RECV_POINTS *pts);
+
+/**
+ * 从已打开的 nc 判断是否为 points 布局
+ *
+ * 优先读全局属性 layout；若无属性则看是否存在 point 维
+ *
+ * @param[in]   ncid   已打开的 nc id
+ * @return      true 表示 points，false 表示 grid
+ */
+bool grt_recv_nc_is_points(int ncid);

pygrt/C_extension/include/grt/static/static_grn.h

@@ -41,4 +41,4 @@ void grt_integ_static_grn(
     realChnlGrid grn[nr],
     realChnlGrid grn_uiz[nr],
     realChnlGrid grn_uir[nr],
-    const char *statsstr);
+    const char *statsstr);