From 2cedd772a0de9f9e7f3922481400ecbed1663af2 Mon Sep 17 00:00:00 2001
From: Dengda98 <zhudengda@mail.iggcas.ac.cn>
Date: Wed, 5 Aug 2026 12:26:23 +0800
Subject: [PATCH] FIX: properly handle zero epicentral distance without
 clamping r

---
 pygrt/C_extension/include/grt/common/const.h  |  4 +-
 pygrt/C_extension/include/grt/common/coord.h  |  4 +-
 .../include/grt/dynamic/dyn_postprocess.h     |  3 +
 .../include/grt/static/static_postprocess.h   |  3 +
 pygrt/C_extension/src/common/coord.c          | 19 ++++-
 pygrt/C_extension/src/dynamic/grt_greenfn.c   | 14 +++-
 pygrt/C_extension/src/dynamic/grt_rotation.c  |  4 +-
 pygrt/C_extension/src/dynamic/grt_strain.c    |  6 +-
 pygrt/C_extension/src/dynamic/grt_stress.c    |  7 +-
 pygrt/C_extension/src/dynamic/grt_syn.c       | 36 ++++++--
 pygrt/C_extension/src/integral/dcm.c          |  3 +-
 pygrt/C_extension/src/integral/dwm.c          | 11 +--
 pygrt/C_extension/src/integral/fim.c          | 13 ++-
 pygrt/C_extension/src/integral/k_integ.c      | 50 ++++++-----
 pygrt/C_extension/src/integral/ptam.c         | 10 +--
 pygrt/C_extension/src/integral/safim.c        | 12 +--
 .../src/static/grt_static_greenfn.c           | 14 +++-
 .../src/static/grt_static_rotation.c          |  6 +-
 .../src/static/grt_static_strain.c            |  8 +-
 .../src/static/grt_static_stress.c            |  8 +-
 pygrt/C_extension/src/static/grt_static_syn.c | 25 ++++--
 pygrt/c_interfaces.py                         | 82 ++++++-------------
 pygrt/c_structures.py                         |  2 +
 pygrt/pymod.py                                | 21 +++--
 pygrt/utils.py                                | 75 +++++++----------
 25 files changed, 235 insertions(+), 205 deletions(-)

diff --git a/pygrt/C_extension/include/grt/common/const.h b/pygrt/C_extension/include/grt/common/const.h
index cb6742e0..905c6b7f 100755
--- a/pygrt/C_extension/include/grt/common/const.h
+++ b/pygrt/C_extension/include/grt/common/const.h
@@ -59,8 +59,8 @@ typedef double complex cplx_t;
 #define GCC_ALWAYS_INLINE __attribute__((always_inline))  ///< gcc编译器不改动内联函数
 
 #define GRT_SWAP(type, a, b) { type temp = a; a = b; b = temp; } ///< 交换两个变量的值
-#define GRT_MIN_DISTANCE    1e-5   ///< 最小震中距，用于限制
-#define GRT_IS_SMALLE_DISTANCE(r)  ((r) <= GRT_MIN_DISTANCE)   ///< 判断是否是过小的震中距 
+#define GRT_ZERO_DISTANCE   1e-8   ///< 判定为零震中距的阈值 (km)
+#define GRT_IS_ZERO(r)  ((r) <= GRT_ZERO_DISTANCE)   ///< 判断震中距是否为零（或数值上视为零）
 
 #define GRT_STRING_FMT "%18s"  ///< 字符串输出格式
 #define GRT_REAL_FMT "%18.8e"  ///< 浮点数输出格式
diff --git a/pygrt/C_extension/include/grt/common/coord.h b/pygrt/C_extension/include/grt/common/coord.h
index 96756bcb..0724c6f5 100644
--- a/pygrt/C_extension/include/grt/common/coord.h
+++ b/pygrt/C_extension/include/grt/common/coord.h
@@ -52,7 +52,7 @@ void grt_rot_zxy2zrt_symtensor2odr(real_t theta, real_t A[6]);
  * 
  * @param[in]       theta      r轴相对x轴的旋转弧度
  * @param[in,out]   u          柱坐标下的位移矢量
- * @param[in,out]   upar       柱坐标下的位移空间偏导
- * @param[in]       r          r轴坐标
+ * @param[in,out]   upar       柱坐标下的位移空间偏导（第三行已是 (1/r)∂_θ 有限部分）
+ * @param[in]       r          r 坐标 (cm)；r=0 时联络项 u/r 改用 ∂_r u
  */
 void grt_rot_zrt2zxy_upar(const real_t theta, real_t u[3], real_t upar[3][3], const real_t r);
\ No newline at end of file
diff --git a/pygrt/C_extension/include/grt/dynamic/dyn_postprocess.h b/pygrt/C_extension/include/grt/dynamic/dyn_postprocess.h
index db06c83b..459d66e5 100644
--- a/pygrt/C_extension/include/grt/dynamic/dyn_postprocess.h
+++ b/pygrt/C_extension/include/grt/dynamic/dyn_postprocess.h
@@ -12,6 +12,9 @@
  * 由动态位移偏导合成应变张量。
  * 数组布局：u[分量][采样点]、upar[偏导方向][分量][采样点]、
  * res[第二分量][第一分量][采样点]。
+ *
+ * ZRT 联络项：r≠0 用 u/r；r=0 改用 ∂_r u（upar[1][*]），
+ * 与 syn 中轴点处 (1/r)∂_θ 有限部分配套。
  */
 void grt_compute_strain(
     size_t npts, float dist, float *const u[GRT_CHANNEL_NUM],
diff --git a/pygrt/C_extension/include/grt/static/static_postprocess.h b/pygrt/C_extension/include/grt/static/static_postprocess.h
index 36b83838..c54b857b 100644
--- a/pygrt/C_extension/include/grt/static/static_postprocess.h
+++ b/pygrt/C_extension/include/grt/static/static_postprocess.h
@@ -13,6 +13,9 @@
  *
  * 数组布局：u[分量][点]、upar[偏导方向][分量][点]、
  * res[第二分量][第一分量][点]。仅写入 res 的上三角分量。
+ *
+ * ZRT 联络项：r≠0 用 u/r；r=0 改用 ∂_r u（upar[1][*]），
+ * 与 syn 中轴点处 (1/r)∂_θ 有限部分配套。
  */
 void grt_static_compute_stress(
     size_t nx, size_t ny, const real_t *xs, const real_t *ys,
diff --git a/pygrt/C_extension/src/common/coord.c b/pygrt/C_extension/src/common/coord.c
index 012b2edc..02f625b7 100644
--- a/pygrt/C_extension/src/common/coord.c
+++ b/pygrt/C_extension/src/common/coord.c
@@ -66,6 +66,17 @@ void grt_rot_zrt2zxy_upar(const real_t theta, real_t u[3], real_t upar[3][3], co
     real_t cct = ct*ct;
     real_t sct = st*ct;
 
+    // 变换含联络项 u_r/r、u_θ/r（r 单位 cm）。
+    // r=0: u/r 联络项改用 ∂_r u_r、∂_r u_θ（s11,s12），与 syn 中 (1/r)∂_θ 有限部分配套。
+    real_t u1_over_r, u2_over_r;
+    if(GRT_IS_ZERO(r * 1e-5)){  // cm → km 后再判零
+        u1_over_r = s11;
+        u2_over_r = s12;
+    } else {
+        u1_over_r = u1/r;
+        u2_over_r = u2/r;
+    }
+
     //           uz       ux       uy
     //  ∂z
     //  ∂x
@@ -82,17 +93,17 @@ void grt_rot_zrt2zxy_upar(const real_t theta, real_t u[3], real_t upar[3][3], co
     // ∂ uz / ∂ x
     upar[1][0] = s10*ct - s20*st;
     // ∂ ux / ∂ x
-    upar[1][1] = s11*cct + s22*sst - (s12+s21)*sct + u1*sst/r + u2*sct/r;
+    upar[1][1] = s11*cct + s22*sst - (s12+s21)*sct + u1_over_r*sst + u2_over_r*sct;
     // ∂ uy / ∂ x
-    upar[1][2] = s12*cct - s21*sst + (s11-s22)*sct - u1*sct/r + u2*sst/r;
+    upar[1][2] = s12*cct - s21*sst + (s11-s22)*sct - u1_over_r*sct + u2_over_r*sst;
 
 
     // ∂ uz / ∂ y
     upar[2][0] = s10*st + s20*ct;
     // ∂ ux / ∂ y
-    upar[2][1] = s21*cct - s12*sst + (s11-s22)*sct - u1*sct/r - u2*cct/r;
+    upar[2][1] = s21*cct - s12*sst + (s11-s22)*sct - u1_over_r*sct - u2_over_r*cct;
     // ∂ uy / ∂ y
-    upar[2][2] = s22*cct + s11*sst + (s12+s21)*sct + u1*cct/r - u2*sct/r;
+    upar[2][2] = s22*cct + s11*sst + (s12+s21)*sct + u1_over_r*cct - u2_over_r*sct;
 
 
     // 转矢量
diff --git a/pygrt/C_extension/src/dynamic/grt_greenfn.c b/pygrt/C_extension/src/dynamic/grt_greenfn.c
index 645b2c83..ec46d947 100644
--- a/pygrt/C_extension/src/dynamic/grt_greenfn.c
+++ b/pygrt/C_extension/src/dynamic/grt_greenfn.c
@@ -295,7 +295,8 @@ printf("\n"
 "                   <length> will be determined automatically\n"
 "                   in program with the criterion (Bouchon, 1980).\n"
 "                 + manually set one POSITIVE <length>, e.g. -L20\n"
-"                 For FIM or SAFIM:\n"
+"                 For FIM or SAFIM (large epicentral distance only;\n"
+"                 zero epicentral distance is not allowed):\n"
 "                 + +l<Flength> defines the dk of the FIM.\n"
 "                 + +a<Ftol> defines the tolerance of the SAFIM.\n"
 "                   you can't set both.\n"
@@ -845,6 +846,15 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){
     fclose(fp);
     GRT_SAFE_FREE_PTR(dummy);
 
+    // FIM/SAFIM 面向远震中距，公式含 1/r、1/√r，不能用于 r=0（含 kcut 分段）
+    if(Ctrl->L.FIM.active || Ctrl->L.SAFIM.active){
+        for(size_t ir=0; ir<Ctrl->R.nr; ++ir){
+            if(GRT_IS_ZERO(Ctrl->R.rs[ir])){
+                GRTBadOptionError(L, "FIM/SAFIM cannot be used with zero epicentral distance.");
+            }
+        }
+    }
+
 }
 
 
@@ -884,7 +894,7 @@ int greenfn_main(int argc, char **argv) {
     Ctrl->N.winT = Ctrl->N.nt*Ctrl->N.dt;
 
     // 最大震中距
-    real_t rmax = Ctrl->R.rs[grt_findMax_real_t(Ctrl->R.rs, Ctrl->R.nr)];   
+    real_t rmax = Ctrl->R.rs[grt_findMax_real_t(Ctrl->R.rs, Ctrl->R.nr)];
 
     // 时窗最大截止时刻
     real_t tmax = 0.0;
diff --git a/pygrt/C_extension/src/dynamic/grt_rotation.c b/pygrt/C_extension/src/dynamic/grt_rotation.c
index c218df9d..d983e579 100644
--- a/pygrt/C_extension/src/dynamic/grt_rotation.c
+++ b/pygrt/C_extension/src/dynamic/grt_rotation.c
@@ -59,8 +59,8 @@ void grt_compute_rotation(
     const char *chs = rot2ZNE ? GRT_ZNE_CODES : GRT_ZRT_CODES;
 
     for(size_t i=0; i<npts; ++i){
-        // ZRT 联络项 u_θ/r，1e-5: km→cm
-        float ut_over_r = u[2][i] / dist * 1e-5f;
+        // 联络项 u_θ/r（1e-5: km→cm）：r≠0 用 u_θ/r；r=0 改用 ∂_r u_θ
+        float ut_over_r = GRT_IS_ZERO(dist) ? upar[1][2][i] : (u[2][i] / dist * 1e-5f);
 
         for(int c=0; c<GRT_CHANNEL_NUM; ++c){
             for(int c2=c+1; c2<GRT_CHANNEL_NUM; ++c2){
diff --git a/pygrt/C_extension/src/dynamic/grt_strain.c b/pygrt/C_extension/src/dynamic/grt_strain.c
index e51f0fa4..2d098c85 100644
--- a/pygrt/C_extension/src/dynamic/grt_strain.c
+++ b/pygrt/C_extension/src/dynamic/grt_strain.c
@@ -57,9 +57,9 @@ void grt_compute_strain(
     const char *chs = rot2ZNE ? GRT_ZNE_CODES : GRT_ZRT_CODES;
 
     for(size_t i=0; i<npts; ++i){
-        // ZRT 联络项 u/r，1e-5: km→cm
-        float ur_over_r = u[1][i] / dist * 1e-5f;
-        float ut_over_r = u[2][i] / dist * 1e-5f;
+        // 联络项（1e-5: km→cm）：r≠0 用 u/r；r=0 改用 ∂_r u，与 syn 中 (1/r)∂_θ 有限部分配套
+        float ur_over_r = GRT_IS_ZERO(dist) ? upar[1][1][i] : (u[1][i] / dist * 1e-5f);
+        float ut_over_r = GRT_IS_ZERO(dist) ? upar[1][2][i] : (u[2][i] / dist * 1e-5f);
 
         for(int c=0; c<GRT_CHANNEL_NUM; ++c){
             for(int c2=c; c2<GRT_CHANNEL_NUM; ++c2){
diff --git a/pygrt/C_extension/src/dynamic/grt_stress.c b/pygrt/C_extension/src/dynamic/grt_stress.c
index ca3d0431..e1a40388 100644
--- a/pygrt/C_extension/src/dynamic/grt_stress.c
+++ b/pygrt/C_extension/src/dynamic/grt_stress.c
@@ -82,12 +82,13 @@ void grt_compute_stress(
         for(size_t i=0; i<nf; ++i)  lam_ukk[i] += fwd->W_f[i];
     }
 
-    // ZRT 联络项 u/r（1e-5: km→cm）；时域先算好再 FFT，避免频域再分支缩放
+    // 联络项（1e-5: km→cm）：r≠0 用 u/r；r=0 改用 ∂_r u，与 syn 中 (1/r)∂_θ 有限部分配套
+    // 时域先算好 ur_over_r / ut_over_r，再 FFT，避免频域再分支缩放
     float *ur_over_r = (float *)malloc(sizeof(float)*npts);
     float *ut_over_r = (float *)malloc(sizeof(float)*npts);
     for(size_t i=0; i<npts; ++i){
-        ur_over_r[i] = u[1][i] / dist * 1e-5f;
-        ut_over_r[i] = u[2][i] / dist * 1e-5f;
+        ur_over_r[i] = GRT_IS_ZERO(dist) ? upar[1][1][i] : (u[1][i] / dist * 1e-5f);
+        ut_over_r[i] = GRT_IS_ZERO(dist) ? upar[1][2][i] : (u[2][i] / dist * 1e-5f);
     }
 
     if(!rot2ZNE){
diff --git a/pygrt/C_extension/src/dynamic/grt_syn.c b/pygrt/C_extension/src/dynamic/grt_syn.c
index 933c34be..a705e8b6 100644
--- a/pygrt/C_extension/src/dynamic/grt_syn.c
+++ b/pygrt/C_extension/src/dynamic/grt_syn.c
@@ -140,6 +140,8 @@ printf("\n"
 "    -G<grn_path>  Green's Functions output directory of module `greenfn`.\n"
 "\n"
 "    -A<azimuth>   Azimuth in degree, from source to station.\n"
+"                  Ignored (forced to 0°) when Green's Functions\n"
+"                  have zero epicentral distance.\n"
 "\n"
 "    -S[u]<scale>  Scale factor to all kinds of source. \n"
 "                  + For Explosion, Shear and Moment Tensor,\n"
@@ -535,17 +537,27 @@ static void syn_accum_from_gf(
  * syn_upar[偏导方向][分量][采样点]。gf_uiz/gf_uir 在 calc_upar=false
  * 时可传 NULL；单个分量指针为 NULL 时跳过该道。
  *
- * @param[in]       azrad     方位角（弧度）
+ * r=0 时强制 *azrad=0（e_r→N、e_θ→E）并告警；
+ * 并用 ∂_r 格林函数合成 (1/r)∂_θ 的有限部分（见函数内注释）。
+ *
+ * @param[in,out]  azrad     方位角（弧度）；r=0 时写回 0
  */
 void grt_syn_from_gf(
     size_t npts, float dist,
     const pfloatChnlGrid gf, const pfloatChnlGrid gf_uiz, const pfloatChnlGrid gf_uir,
-    GRT_SYN_TYPE computeType, real_t M0, real_t VpVs_ratio, real_t azrad,
+    GRT_SYN_TYPE computeType, real_t M0, real_t VpVs_ratio, real_t *azrad,
     const real_t mchn[GRT_MECHANISM_NUM],
     bool rot2ZNE, bool calc_upar,
     float *const syn[GRT_CHANNEL_NUM], float *const syn_upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM])
 {
-    const real_t az = azrad;
+    // r=0：方位角无定义，约定 e_r→N、e_θ→E，故强制 az=0
+    if(GRT_IS_ZERO(dist)){
+        GRTRaiseWarning(
+            "Zero epicentral distance: azimuth is ignored "
+            "(forced to 0°; e_r→N, e_θ→E).");
+        *azrad = 0.0;
+    }
+    const real_t az = *azrad;
 
     int calcUTypes = calc_upar ? 4 : 1;
     realChnlGrid srcRadi = {0};
@@ -563,14 +575,16 @@ void grt_syn_from_gf(
     for(int ityp = 0; ityp < calcUTypes; ++ityp){
         real_t upar_scale = 1.0;
         // 求位移空间导数时，需调整比例系数（1e-5: km→cm）
-        // ZRT 协变导数拆两步：此处算 (1/r)∂_θ u，后处理再补 ±u/r。
+        // ZRT 协变导数拆两步：此处合成 (1/r)∂_θ u，后处理再补 ±u/r。
+        // r=0 时协变组合仍有限，两项分别换成有限极限（见下方 gf 与后处理）。
         if(ityp > 0){
             switch (GRT_ZRT_CODES[ityp-1]){
                 case 'Z': case 'R':
                     upar_scale = 1e-5;
                     break;
+                // (1/r)∂_θ：r≠0 时 scale∝1/r；r=0 时改用 ∂_r GF，scale 仅留 km→cm
                 case 'T':
-                    upar_scale = 1e-5 / dist;
+                    upar_scale = GRT_IS_ZERO(dist) ? 1e-5 : (1e-5 / dist);
                     break;
                 default:
                     break;
@@ -582,6 +596,11 @@ void grt_syn_from_gf(
             up = gf_uiz;
         } else if(ityp == 2){
             up = gf_uir;
+        } else if(ityp == 3 && GRT_IS_ZERO(dist)){
+            // r=0: 用 ∂_r GF（par_θ 辐射因子）合成 (1/r)∂_θ 的有限部分；
+            // 后处理中 u/r 联络项改用 ∂_r u，二者合并得有限直角/柱坐标导数。
+            // 对 u_z：m≥1 ⇒ u_z(0)=0，lim u_z/r=∂_r u_z；m=0 无 ∂_θ；无 ±u_z/r 联络项。
+            up = gf_uir;
         }
 
         memset(srcRadi, 0, sizeof(srcRadi));
@@ -752,9 +771,14 @@ int syn_main(int argc, char **argv)
     grt_syn_from_gf(
         (size_t)npts, Ctrl->dist,
         gf, calc_upar ? gf_uiz : NULL, calc_upar ? gf_uir : NULL,
-        Ctrl->computeType, Ctrl->S.M0, Ctrl->VpVs_ratio, Ctrl->A.azrad, Ctrl->mchn,
+        Ctrl->computeType, Ctrl->S.M0, Ctrl->VpVs_ratio, &Ctrl->A.azrad, Ctrl->mchn,
         false, calc_upar, syn, syn_upar);
 
+    // C 可能因 r=0 强制 azrad=0，同步方位角头段
+    Ctrl->A.azimuth = Ctrl->A.azrad / DEG1;
+    Ctrl->A.backazimuth = Ctrl->A.azimuth + 180.0;
+    if(Ctrl->A.backazimuth >= 360.0) Ctrl->A.backazimuth -= 360.0;
+
     // 时间函数 / 积分 / 微分（在旋转前，与历史行为一致）
     SACTRACE *tfsac = NULL;
     if(Ctrl->D.active){
diff --git a/pygrt/C_extension/src/integral/dcm.c b/pygrt/C_extension/src/integral/dcm.c
index 6ee8bb7b..7849c824 100644
--- a/pygrt/C_extension/src/integral/dcm.c
+++ b/pygrt/C_extension/src/integral/dcm.c
@@ -34,7 +34,8 @@ void grt_dcm_correction(size_t nr, real_t *rs, real_t dk, real_t kcut, K_INTEG *
 {
     for(size_t ir = 0; ir < nr; ++ir){
         real_t r = rs[ir];
-        if(GRT_IS_SMALLE_DISTANCE(r)) continue;
+        // 修正系数含 1/r；r=0 跳过（近场极限已在 k_integ 的 Bessel 极限中处理）
+        if(GRT_IS_ZERO(r)) continue;
 
         real_t c = 1.0 / r;
 
diff --git a/pygrt/C_extension/src/integral/dwm.c b/pygrt/C_extension/src/integral/dwm.c
index 8714ec8b..a5e2a218 100644
--- a/pygrt/C_extension/src/integral/dwm.c
+++ b/pygrt/C_extension/src/integral/dwm.c
@@ -70,13 +70,14 @@ real_t grt_discrete_integ(
         for(size_t ir = 0; ir < nr; ++ir){
             if(iendkrs[ir]) continue; // 该震中距下的波数k积分已收敛
 
-            // 跳过奇异点
-            if(depsrc == deprcv && GRT_IS_SMALLE_DISTANCE(rs[ir]))  continue;
+            // 震源-接收点重合且 r=0：格林函数奇异
+            if(depsrc == deprcv && GRT_IS_ZERO(rs[ir]))  continue;
 
             memset(K->SUM, 0, sizeof(cplxIntegGrid));
             
             // 计算被积函数一项 F(k,w)Jm(kr)k
-            grt_int_Pk(k, rs[ir], (K->applyDCM && GRT_IS_SMALLE_DISTANCE(rs[ir]))? K->QWV_raw : K->QWV, false, K->SUM);
+            // r=0 时 DCM 解析修正被跳过，故改用未扣除近场的 QWV_raw；近场由 k_integ 极限处理
+            grt_int_Pk(k, rs[ir], (K->applyDCM && GRT_IS_ZERO(rs[ir]))? K->QWV_raw : K->QWV, false, K->SUM);
             
             iendk0 = true;
 
@@ -102,7 +103,7 @@ real_t grt_discrete_integ(
             if(K->calc_upar){
                 // ------------------------------- ui_z -----------------------------------
                 // 计算被积函数一项 F(k,w)Jm(kr)k
-                grt_int_Pk(k, rs[ir], (K->applyDCM && GRT_IS_SMALLE_DISTANCE(rs[ir]))? K->QWVz_raw : K->QWVz, false, K->SUM);
+                grt_int_Pk(k, rs[ir], (K->applyDCM && GRT_IS_ZERO(rs[ir]))? K->QWVz_raw : K->QWVz, false, K->SUM);
                 
                 // keps不参与计算位移空间导数的积分，背后逻辑认为u收敛，则uiz也收敛
                 GRT_LOOP_IntegGrid(im, v){
@@ -111,7 +112,7 @@ real_t grt_discrete_integ(
 
                 // ------------------------------- ui_r -----------------------------------
                 // 计算被积函数一项 F(k,w)Jm(kr)k
-                grt_int_Pk(k, rs[ir], (K->applyDCM && GRT_IS_SMALLE_DISTANCE(rs[ir]))? K->QWV_raw : K->QWV, true, K->SUM);
+                grt_int_Pk(k, rs[ir], (K->applyDCM && GRT_IS_ZERO(rs[ir]))? K->QWV_raw : K->QWV, true, K->SUM);
                 
                 // keps不参与计算位移空间导数的积分，背后逻辑认为u收敛，则uiz也收敛
                 GRT_LOOP_IntegGrid(im, v){
diff --git a/pygrt/C_extension/src/integral/fim.c b/pygrt/C_extension/src/integral/fim.c
index 7194e13b..4407fdd9 100755
--- a/pygrt/C_extension/src/integral/fim.c
+++ b/pygrt/C_extension/src/integral/fim.c
@@ -29,9 +29,6 @@ real_t grt_linear_filon_integ(
 {   
     if(k0 + dk0 >= kmax)  return k0;
 
-    real_t depsrc = mstat->mod1d->depsrc;
-    real_t deprcv = mstat->mod1d->deprcv;
-    
     // 从0开始，存储第二部分Filon积分的结果
     K_INTEG *K2 = grt_init_K_INTEG(K->calc_upar, nr);
 
@@ -72,8 +69,8 @@ real_t grt_linear_filon_integ(
         for(size_t ir = 0; ir < nr; ++ir){
             if(iendkrs[ir]) continue; // 该震中距下的波数k积分已收敛
 
-            // 跳过奇异点
-            if(depsrc == deprcv && GRT_IS_SMALLE_DISTANCE(rs[ir]))  continue;
+            // 防御：r=0 应在参数入口已拒绝（FIM 不适于零震中距）
+            if(GRT_IS_ZERO(rs[ir]))  continue;
 
             memset(K2->SUM, 0, sizeof(cplxIntegGrid));
             
@@ -180,8 +177,8 @@ real_t grt_linear_filon_integ(
         for(size_t ir = 0; ir < nr; ++ir){
             real_t r = rs[ir];
 
-            // 跳过奇异点
-            if(depsrc == deprcv && GRT_IS_SMALLE_DISTANCE(rs[ir]))  continue;
+            // 防御：r=0 应在参数入口已拒绝（FIM 不适于零震中距）
+            if(GRT_IS_ZERO(rs[ir]))  continue;
 
             // Gc
             grt_int_Pk_filon(k0N, r, true, K->QWV, false, SUM_Gc[iik]);
@@ -232,6 +229,8 @@ real_t grt_linear_filon_integ(
     // 乘上总系数 sqrt(2.0/(PI*r)) / dk0,  除dks0是在该函数外还会再乘dk0, 并将结果加到原数组中
     for(size_t ir = 0; ir < nr; ++ir){
         real_t r = rs[ir];
+        // 防御：r=0 应在参数入口已拒绝（总系数含 1/√r）
+        if(GRT_IS_ZERO(r))  continue;
         real_t tmp = sqrt(2.0/(PI*r)) / dk0;
 
         GRT_LOOP_IntegGrid(im, v){
diff --git a/pygrt/C_extension/src/integral/k_integ.c b/pygrt/C_extension/src/integral/k_integ.c
index 33ddf781..5e840b08 100644
--- a/pygrt/C_extension/src/integral/k_integ.c
+++ b/pygrt/C_extension/src/integral/k_integ.c
@@ -86,34 +86,44 @@ void grt_int_Pk(real_t k, real_t r, const cplxChnlGrid QWV, bool calc_uir, cplxI
 {
     real_t bjmk[GRT_MORDER_MAX+1] = {0};
     real_t kr = k*r;
-    real_t kr_inv = 1.0/kr;
     real_t kcoef = k;
 
     real_t Jmcoef[GRT_MORDER_MAX+1] = {0};
 
-    grt_bessel012(kr, &bjmk[0], &bjmk[1], &bjmk[2]); 
-    if(calc_uir){
-        real_t bjmk0[GRT_MORDER_MAX+1] = {0};
-        for(int i=0; i<=GRT_MORDER_MAX; ++i)  bjmk0[i] = bjmk[i];
-
-        if(GRT_IS_SMALLE_DISTANCE(r)){
-            bjmk[0] = - bjmk0[1];
-            bjmk[1] = bjmk0[0] - 0.5;
-            bjmk[2] = bjmk0[1] - 0.25*kr;
-            Jmcoef[1] = - 0.125 * kr;
+    // r=0 ⇒ kr≡0（任意 k）。近场项含 J_m(kr)/r、∂_r(J_m/r)，直接取 x=kr→0 极限，
+    // 避免 1/kr。仅当 GRT_IS_ZERO(r) 时成立，不可对一般小 r、大 k 套用。
+    if(GRT_IS_ZERO(r)){
+        if(calc_uir){
+            // bjmk ← J_m'(0): J0'=0, J1'=1/2, J2'=0
+            bjmk[0] = 0.0;
+            bjmk[1] = 0.5;
+            bjmk[2] = 0.0;
+            // Jmcoef ← lim (J_m' - J_m/x)/x ；×k² 后即 lim ∂_r(J_m/r)
+            // m=1 → 0, m=2 → 1/8
+            Jmcoef[1] = 0.0;
             Jmcoef[2] = 0.125;
+            kcoef = k*k;
         } else {
-            grt_besselp012(kr, &bjmk[0], &bjmk[1], &bjmk[2]); 
-            for(int i=1; i<=GRT_MORDER_MAX; ++i)  Jmcoef[i] = kr_inv * (-kr_inv * bjmk0[i] + bjmk[i]);
-        }
-
-        kcoef = k*k;
-    } 
-    else {
-        if(GRT_IS_SMALLE_DISTANCE(r)){
+            // bjmk ← J_m(0): J0=1, J1=J2=0
+            bjmk[0] = 1.0;
+            bjmk[1] = 0.0;
+            bjmk[2] = 0.0;
+            // Jmcoef ← lim J_m(x)/x ；×k 后即 lim J_m/r
+            // m=1 → 1/2, m=2 → 0
             Jmcoef[1] = 0.5;
-            Jmcoef[2] = 0.125*kr;
+            Jmcoef[2] = 0.0;
+        }
+    } else {
+        grt_bessel012(kr, &bjmk[0], &bjmk[1], &bjmk[2]);
+        if(calc_uir){
+            real_t bjmk0[GRT_MORDER_MAX+1] = {0};
+            for(int i=0; i<=GRT_MORDER_MAX; ++i)  bjmk0[i] = bjmk[i];
+            real_t kr_inv = 1.0/kr;
+            grt_besselp012(kr, &bjmk[0], &bjmk[1], &bjmk[2]);
+            for(int i=1; i<=GRT_MORDER_MAX; ++i)  Jmcoef[i] = kr_inv * (-kr_inv * bjmk0[i] + bjmk[i]);
+            kcoef = k*k;
         } else {
+            real_t kr_inv = 1.0/kr;
             for(int i=1; i<=GRT_MORDER_MAX; ++i)  Jmcoef[i] = bjmk[i]*kr_inv;
         }
     }
diff --git a/pygrt/C_extension/src/integral/ptam.c b/pygrt/C_extension/src/integral/ptam.c
index 385c660f..f869200e 100755
--- a/pygrt/C_extension/src/integral/ptam.c
+++ b/pygrt/C_extension/src/integral/ptam.c
@@ -259,14 +259,11 @@ void grt_PTA_method(
     #undef __ARR
     #undef __CALLOC_ARRAY
 
-    real_t depsrc = mstat->mod1d->depsrc;
-    real_t deprcv = mstat->mod1d->deprcv;
-
     // 对于PTAM，不同震中距使用不同dk
     for(size_t ir = 0; ir < nr; ++ir){
 
-        // 跳过奇异点
-        if(depsrc == deprcv && GRT_IS_SMALLE_DISTANCE(rs[ir]))  continue;
+        // 零震中距：PTAM 步长含 1/r；震源-接收点重合时格林函数奇异
+        if(GRT_IS_ZERO(rs[ir]))  continue;
 
 
         real_t dk = PI/((GRT_PTAM_WAITS_MAX-1)*rs[ir]); 
@@ -316,8 +313,7 @@ void grt_PTA_method(
 
     // 做缩减序列，赋值最终解
     for(size_t ir = 0; ir < nr; ++ir){
-        // 跳过奇异点
-        if(depsrc == deprcv && GRT_IS_SMALLE_DISTANCE(rs[ir]))  continue;
+        if(GRT_IS_ZERO(rs[ir]))  continue;
 
         // 记录到文件
         if(ptam_fstatsnr != NULL)  grt_write_stats_ptam(ptam_fstatsnr[ir][1], Kpt[ir], (K->calc_upar)? Fpt_uiz[ir] : Fpt[ir]);
diff --git a/pygrt/C_extension/src/integral/safim.c b/pygrt/C_extension/src/integral/safim.c
index d64caa06..b817540b 100644
--- a/pygrt/C_extension/src/integral/safim.c
+++ b/pygrt/C_extension/src/integral/safim.c
@@ -232,13 +232,13 @@ static bool check_fit(
  * 以下实际拟合的二次函数是 sqrt(k)*F(k,w), 这样积分时可以避免计算超越函数
  * 
  */
-static void interv_integ(const KInterval *ptKitv, size_t nr, real_t *rs, K_INTEG *K, real_t depsrc, real_t deprcv)
+static void interv_integ(const KInterval *ptKitv, size_t nr, real_t *rs, K_INTEG *K)
 {
     // 震中距rs循环
     for(size_t ir = 0; ir < nr; ++ir){
 
-        // 跳过奇异点
-        if(depsrc == deprcv && GRT_IS_SMALLE_DISTANCE(rs[ir]))  continue;
+        // 防御：r=0 应在参数入口已拒绝（SAFIM 不适于零震中距）
+        if(GRT_IS_ZERO(rs[ir]))  continue;
 
 
         memset(K->SUM, 0, sizeof(cplxIntegGrid));
@@ -285,8 +285,6 @@ real_t grt_sa_filon_integ(
     MODEL1D_STATE *mstat, real_t k0, real_t dk0, real_t tol, real_t kmax, real_t kref,
     size_t nr, real_t *rs, K_INTEG *K, FILE *fstats, GRT_KernelFunc kerfunc)
 {   
-    real_t depsrc = mstat->mod1d->depsrc;
-    real_t deprcv = mstat->mod1d->deprcv;
     
     real_t kmin = k0 + dk0;
     if(kmin >= kmax)  return k0;
@@ -433,12 +431,14 @@ real_t grt_sa_filon_integ(
                 }
             }
             // 计算积分
-            interv_integ(&Kitv, nr, rs, K2, depsrc, deprcv);
+            interv_integ(&Kitv, nr, rs, K2);
         }
     } // END sampling
 
     // 乘上总系数 sqrt(2.0/(PI*r)) / dk0,  除dks0是在该函数外还会再乘dk0, 并将结果加到原数组中
     for(size_t ir = 0; ir < nr; ++ir){
+        // 防御：r=0 应在参数入口已拒绝（总系数含 1/√r）
+        if(GRT_IS_ZERO(rs[ir]))  continue;
         real_t tmp = sqrt(2.0/(PI*rs[ir])) / dk0;
 
         GRT_LOOP_IntegGrid(im, v){
diff --git a/pygrt/C_extension/src/static/grt_static_greenfn.c b/pygrt/C_extension/src/static/grt_static_greenfn.c
index 6dc7c73b..f48cedf1 100644
--- a/pygrt/C_extension/src/static/grt_static_greenfn.c
+++ b/pygrt/C_extension/src/static/grt_static_greenfn.c
@@ -195,7 +195,8 @@ printf("\n"
 "                 + (default) not set or set 0.0. \n"
 "                   <length> will be %.1f.\n", GRT_GREENFN_L_LENGTH); printf(
 "                 + manually set one POSITIVE <length>, e.g. -L20\n"
-"                 For FIM or SAFIM:\n"
+"                 For FIM or SAFIM (large epicentral distance only;\n"
+"                 zero epicentral distance is not allowed):\n"
 "                 + +l<Flength> defines the dk of the FIM.\n"
 "                 + +a<Ftol> defines the tolerance of the SAFIM.\n"
 "                   you can't set both.\n"
@@ -555,7 +556,16 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){
     Ctrl->rs = (real_t*)calloc(Ctrl->nr, sizeof(real_t));
     for(size_t ix=0; ix<Ctrl->X.nx; ++ix){
         for(size_t iy=0; iy<Ctrl->Y.ny; ++iy){
-            Ctrl->rs[iy + ix*Ctrl->Y.ny] = GRT_MAX(hypot(Ctrl->X.xs[ix], Ctrl->Y.ys[iy]), GRT_MIN_DISTANCE);  // 避免0震中距
+            Ctrl->rs[iy + ix*Ctrl->Y.ny] = hypot(Ctrl->X.xs[ix], Ctrl->Y.ys[iy]);
+        }
+    }
+
+    // FIM/SAFIM 面向远震中距，公式含 1/r、1/√r，不能用于 r=0（含 kcut 分段）
+    if(Ctrl->L.FIM.active || Ctrl->L.SAFIM.active){
+        for(size_t ir=0; ir<Ctrl->nr; ++ir){
+            if(GRT_IS_ZERO(Ctrl->rs[ir])){
+                GRTBadOptionError(L, "FIM/SAFIM cannot be used with zero epicentral distance.");
+            }
         }
     }
 
diff --git a/pygrt/C_extension/src/static/grt_static_rotation.c b/pygrt/C_extension/src/static/grt_static_rotation.c
index 78099ba6..8d01a498 100644
--- a/pygrt/C_extension/src/static/grt_static_rotation.c
+++ b/pygrt/C_extension/src/static/grt_static_rotation.c
@@ -60,9 +60,9 @@ void grt_static_compute_rotation(
     for(size_t ix=0; ix<nx; ++ix){
         for(size_t iy=0; iy<ny; ++iy){
             size_t ir = iy + ix*ny;
-            real_t dist = GRT_MAX(hypot(xs[ix], ys[iy]), GRT_MIN_DISTANCE);
-            // ZRT 联络项 u_θ/r，1e-5: km→cm
-            real_t ut_over_r = u[2][ir] / dist * 1e-5;
+            real_t dist = hypot(xs[ix], ys[iy]);
+            // 联络项 u_θ/r（1e-5: km→cm）：r≠0 用 u_θ/r；r=0 改用 ∂_r u_θ
+            real_t ut_over_r = GRT_IS_ZERO(dist) ? upar[1][2][ir] : (u[2][ir] / dist * 1e-5);
 
             for(int c=0; c<GRT_CHANNEL_NUM; ++c){
                 for(int c2=c+1; c2<GRT_CHANNEL_NUM; ++c2){
diff --git a/pygrt/C_extension/src/static/grt_static_strain.c b/pygrt/C_extension/src/static/grt_static_strain.c
index 94f5d631..f76dca19 100644
--- a/pygrt/C_extension/src/static/grt_static_strain.c
+++ b/pygrt/C_extension/src/static/grt_static_strain.c
@@ -58,10 +58,10 @@ void grt_static_compute_strain(
     for(size_t ix=0; ix<nx; ++ix){
         for(size_t iy=0; iy<ny; ++iy){
             size_t ir = iy + ix*ny;
-            real_t dist = GRT_MAX(hypot(xs[ix], ys[iy]), GRT_MIN_DISTANCE);
-            // ZRT 联络项 u/r，1e-5: km→cm
-            real_t ur_over_r = u[1][ir] / dist * 1e-5;
-            real_t ut_over_r = u[2][ir] / dist * 1e-5;
+            real_t dist = hypot(xs[ix], ys[iy]);
+            // 联络项（1e-5: km→cm）：r≠0 用 u/r；r=0 改用 ∂_r u，与 syn 中 (1/r)∂_θ 有限部分配套
+            real_t ur_over_r = GRT_IS_ZERO(dist) ? upar[1][1][ir] : (u[1][ir] / dist * 1e-5);
+            real_t ut_over_r = GRT_IS_ZERO(dist) ? upar[1][2][ir] : (u[2][ir] / dist * 1e-5);
 
             for(int c=0; c<GRT_CHANNEL_NUM; ++c){
                 for(int c2=c; c2<GRT_CHANNEL_NUM; ++c2){
diff --git a/pygrt/C_extension/src/static/grt_static_stress.c b/pygrt/C_extension/src/static/grt_static_stress.c
index 85780dc3..d4c1e176 100644
--- a/pygrt/C_extension/src/static/grt_static_stress.c
+++ b/pygrt/C_extension/src/static/grt_static_stress.c
@@ -60,10 +60,10 @@ void grt_static_compute_stress(
     for(size_t ix=0; ix<nx; ++ix){
         for(size_t iy=0; iy<ny; ++iy){
             size_t ir = iy + ix*ny;
-            real_t dist = GRT_MAX(hypot(xs[ix], ys[iy]), GRT_MIN_DISTANCE);
-            // ZRT 联络项 u/r，1e-5: km→cm
-            real_t ur_over_r = u[1][ir] / dist * 1e-5;
-            real_t ut_over_r = u[2][ir] / dist * 1e-5;
+            real_t dist = hypot(xs[ix], ys[iy]);
+            // 联络项（1e-5: km→cm）：r≠0 用 u/r；r=0 改用 ∂_r u，与 syn 中 (1/r)∂_θ 有限部分配套
+            real_t ur_over_r = GRT_IS_ZERO(dist) ? upar[1][1][ir] : (u[1][ir] / dist * 1e-5);
+            real_t ut_over_r = GRT_IS_ZERO(dist) ? upar[1][2][ir] : (u[2][ir] / dist * 1e-5);
 
             real_t lam_ukk = upar[0][0][ir] + upar[1][1][ir] + upar[2][2][ir];
             if(!rot2ZNE)  lam_ukk += ur_over_r;
diff --git a/pygrt/C_extension/src/static/grt_static_syn.c b/pygrt/C_extension/src/static/grt_static_syn.c
index 43a22051..0dbcf6f3 100644
--- a/pygrt/C_extension/src/static/grt_static_syn.c
+++ b/pygrt/C_extension/src/static/grt_static_syn.c
@@ -408,13 +408,16 @@ static void static_syn_from_gf_one(
     for(int ityp = 0; ityp < calcUTypes; ++ityp){
         real_t upar_scale = 1.0;
         // 求位移空间导数时，需调整比例系数（1e-5: km→cm）
+        // ZRT 协变导数拆两步：此处合成 (1/r)∂_θ u，后处理再补 ±u/r。
+        // r=0 时协变组合仍有限，两项分别换成有限极限（见下方 up 与后处理）。
         if(ityp > 0){
             switch (GRT_ZRT_CODES[ityp-1]){
                 case 'Z': case 'R':
                     upar_scale = 1e-5;
                     break;
+                // (1/r)∂_θ：r≠0 时 scale∝1/r；r=0 时改用 ∂_r GF，scale 仅留 km→cm
                 case 'T':
-                    upar_scale = 1e-5 / dist0;
+                    upar_scale = GRT_IS_ZERO(dist0) ? 1e-5 : (1e-5 / dist0);
                     break;
                 default:
                     break;
@@ -426,6 +429,11 @@ static void static_syn_from_gf_one(
             up = uiz;
         } else if(ityp == 2){
             up = uir;
+        } else if(ityp == 3 && GRT_IS_ZERO(dist0)){
+            // r=0: 用 ∂_r GF（par_θ 辐射因子）合成 (1/r)∂_θ 的有限部分；
+            // 后处理中 u/r 联络项改用 ∂_r u，二者合并得有限直角/柱坐标导数。
+            // 对 u_z：m≥1 ⇒ u_z(0)=0，lim u_z/r=∂_r u_z；m=0 无 ∂_θ；无 ±u_z/r 联络项。
+            up = uir;
         }
 
         memset(tmpsyn, 0, sizeof(tmpsyn));
@@ -448,6 +456,7 @@ static void static_syn_from_gf_one(
 
     if(rot2ZNE){
         if(calc_upar){
+            // dist0: km→cm；r=0 时 coord 内 u/r 联络项改用 ∂_r u
             grt_rot_zrt2zxy_upar(azrad, syn, syn_upar, dist0 * 1e5);
         } else {
             grt_rot_zxy2zrt_vec(-azrad, syn);
@@ -460,7 +469,7 @@ static void static_syn_from_gf_one(
  * 由静态格林函数合成三分量位移场（及可选空间偏导）。
  *
  * 对应动态解的 grt_syn_from_gf。输入为原 XY 网格上的格林函数，
- * 可插值到新 XY 网格。
+ * 可插值到新 XY 网格；r=0 时强制方位角为 0（e_r→N、e_θ→E）。
  *
  * 数组布局：u[震中距点][震源][分量]、syn[新网格点][分量]、
  * syn_upar[新网格点][偏导方向][分量]。uiz/uir 在 calc_upar=false 时可传 NULL。
@@ -482,7 +491,7 @@ void grt_static_syn_from_gf(
     for(size_t ix = 0; ix < nx0; ++ix){
         for(size_t iy = 0; iy < ny0; ++iy){
             size_t idx = iy + ix*ny0;
-            rs0[idx] = GRT_MAX(hypot(xs0[ix], ys0[iy]), GRT_MIN_DISTANCE);
+            rs0[idx] = hypot(xs0[ix], ys0[iy]);
             sort_rs0_idx[idx] = idx;
         }
     }
@@ -515,18 +524,18 @@ void grt_static_syn_from_gf(
             real_t y = ys[iy];
 
             // 震中距
-            real_t dist = GRT_MAX(hypot(x, y), GRT_MIN_DISTANCE);
+            real_t dist = hypot(x, y);
 
-            // 方位角
-            real_t azrad = atan2(y, x);
+            // 方位角；r=0 时 atan2(0,0) 无定义，约定 e_r→N、e_θ→E ⇒ az=0
+            real_t azrad = GRT_IS_ZERO(dist) ? 0.0 : atan2(y, x);
 
             size_t ir = iy + ix * ny;
 
             memset(syn[ir], 0, sizeof(syn[ir]));
             memset(syn_upar[ir], 0, sizeof(syn_upar[ir]));
 
-            // 检查是否越界
-            bool r_OutofBound = (dist < sort_rs0[0] || dist > sort_rs0[nr0-1] + 1e-8);
+            // 检查是否越界（允许查询点为精确的零震中距）
+            bool r_OutofBound = (dist < sort_rs0[0] - 1e-8 || dist > sort_rs0[nr0-1] + 1e-8);
             if(r_OutofBound){
                 GRTRaiseWarning("(x, y)=(%.3e, %.3e) is out of distance bounds, skip.", x, y);
                 continue;
diff --git a/pygrt/c_interfaces.py b/pygrt/c_interfaces.py
index 2b3fcfe9..87e4d06c 100755
--- a/pygrt/c_interfaces.py
+++ b/pygrt/c_interfaces.py
@@ -17,7 +17,6 @@
 
 FPOINTER = POINTER(c_float)
 IPOINTER = POINTER(c_int)
-DPOINTER = POINTER(c_double)
 
 
 libgrt = cdll.LoadLibrary(
@@ -44,6 +43,8 @@
     c_char_p
 ]
 
+
+
 C_grt_static_syn_from_gf = libgrt.grt_static_syn_from_gf
 """由静态格林函数合成三分量位移场（及可选空间偏导），可插值到新 XY 网格"""
 C_grt_static_syn_from_gf.restype = None
@@ -58,14 +59,6 @@
     POINTER(REAL * CHANNEL_NUM), POINTER((REAL * CHANNEL_NUM) * CHANNEL_NUM),
 ]
 
-# 通道维指针容器，对应 C 侧常用布局：
-#   T *arr[GRT_CHANNEL_NUM]                      → *_CHNL_PTRS
-#   T *arr[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]      → *_CHNL_PTR_MAT
-FLOAT_CHNL_PTRS = FPOINTER * CHANNEL_NUM
-FLOAT_CHNL_PTR_MAT = FLOAT_CHNL_PTRS * CHANNEL_NUM
-REAL_CHNL_PTRS = PREAL * CHANNEL_NUM
-REAL_CHNL_PTR_MAT = REAL_CHNL_PTRS * CHANNEL_NUM
-
 C_grt_syn_from_gf = libgrt.grt_syn_from_gf
 """由动态格林函数合成三分量地震图（及可选空间偏导）"""
 C_grt_syn_from_gf.restype = None
@@ -74,7 +67,7 @@
     POINTER(FPOINTER * CHANNEL_NUM),
     POINTER(FPOINTER * CHANNEL_NUM),
     POINTER(FPOINTER * CHANNEL_NUM),
-    c_int, REAL, REAL, REAL, REAL * MECHANISM_NUM,
+    c_int, REAL, REAL, PREAL, REAL * MECHANISM_NUM,
     c_bool, c_bool,
     FPOINTER * CHANNEL_NUM, POINTER(FPOINTER * CHANNEL_NUM),
 ]
@@ -85,23 +78,29 @@
 C_grt_compute_stress.restype = None
 C_grt_compute_stress.argtypes = [
     c_size_t, c_float, c_float, c_float, c_float, c_float, c_float, c_float,
-    FLOAT_CHNL_PTRS, POINTER(FLOAT_CHNL_PTRS), POINTER(FLOAT_CHNL_PTRS), c_bool,
+    FPOINTER * CHANNEL_NUM,
+    POINTER(FPOINTER * CHANNEL_NUM), POINTER(FPOINTER * CHANNEL_NUM),
+    c_bool,
 ]
 
 C_grt_compute_strain = libgrt.grt_compute_strain
 """由动态位移偏导合成应变张量"""
 C_grt_compute_strain.restype = None
 C_grt_compute_strain.argtypes = [
-    c_size_t, c_float, FLOAT_CHNL_PTRS, POINTER(FLOAT_CHNL_PTRS),
-    POINTER(FLOAT_CHNL_PTRS), c_bool,
+    c_size_t, c_float,
+    FPOINTER * CHANNEL_NUM,
+    POINTER(FPOINTER * CHANNEL_NUM), POINTER(FPOINTER * CHANNEL_NUM),
+    c_bool,
 ]
 
 C_grt_compute_rotation = libgrt.grt_compute_rotation
 """由动态位移偏导合成旋转张量"""
 C_grt_compute_rotation.restype = None
 C_grt_compute_rotation.argtypes = [
-    c_size_t, c_float, FLOAT_CHNL_PTRS, POINTER(FLOAT_CHNL_PTRS),
-    POINTER(FLOAT_CHNL_PTRS), c_bool,
+    c_size_t, c_float,
+    FPOINTER * CHANNEL_NUM,
+    POINTER(FPOINTER * CHANNEL_NUM), POINTER(FPOINTER * CHANNEL_NUM),
+    c_bool,
 ]
 
 C_grt_static_compute_stress = libgrt.grt_static_compute_stress
@@ -109,7 +108,8 @@
 C_grt_static_compute_stress.restype = None
 C_grt_static_compute_stress.argtypes = [
     c_size_t, c_size_t, PREAL, PREAL,
-    REAL_CHNL_PTRS, POINTER(REAL_CHNL_PTRS), POINTER(REAL_CHNL_PTRS),
+    PREAL * CHANNEL_NUM,
+    POINTER(PREAL * CHANNEL_NUM), POINTER(PREAL * CHANNEL_NUM),
     c_bool, REAL, REAL,
 ]
 
@@ -118,7 +118,9 @@
 C_grt_static_compute_strain.restype = None
 C_grt_static_compute_strain.argtypes = [
     c_size_t, c_size_t, PREAL, PREAL,
-    REAL_CHNL_PTRS, POINTER(REAL_CHNL_PTRS), POINTER(REAL_CHNL_PTRS), c_bool,
+    PREAL * CHANNEL_NUM,
+    POINTER(PREAL * CHANNEL_NUM), POINTER(PREAL * CHANNEL_NUM),
+    c_bool,
 ]
 
 C_grt_static_compute_rotation = libgrt.grt_static_compute_rotation
@@ -126,13 +128,15 @@
 C_grt_static_compute_rotation.restype = None
 C_grt_static_compute_rotation.argtypes = [
     c_size_t, c_size_t, PREAL, PREAL,
-    REAL_CHNL_PTRS, POINTER(REAL_CHNL_PTRS), POINTER(REAL_CHNL_PTRS), c_bool,
+    PREAL * CHANNEL_NUM,
+    POINTER(PREAL * CHANNEL_NUM), POINTER(PREAL * CHANNEL_NUM),
+    c_bool,
 ]
 
 
 C_grt_set_num_threads = libgrt.grt_set_num_threads
 """设置多线程数"""
-C_grt_set_num_threads.restype = None 
+C_grt_set_num_threads.restype = None
 C_grt_set_num_threads.argtypes = [c_int]
 
 
@@ -193,46 +197,6 @@ def set_num_threads(n):
 C_grt_get_ricker_wave.argtypes = [c_float, c_float, IPOINTER]
 
 
-# -------------------------------------------------------------------
-#                      C函数定义的旋转函数
-# -------------------------------------------------------------------
-C_grt_rot_zxy2zrt_vec = libgrt.grt_rot_zxy2zrt_vec
-"""直角坐标zxy到柱坐标zrt的矢量旋转"""
-C_grt_rot_zxy2zrt_vec.restype = None
-C_grt_rot_zxy2zrt_vec.argtypes = [c_double, DPOINTER]  # double, double[3]
-
-C_grt_rot_zxy2zrt_symtensor2odr = libgrt.grt_rot_zxy2zrt_symtensor2odr
-"""直角坐标zxy到柱坐标zrt的二阶对称张量旋转"""
-C_grt_rot_zxy2zrt_symtensor2odr.restype = None
-C_grt_rot_zxy2zrt_symtensor2odr.argtypes = [c_double, DPOINTER]  # double, double[6]
-
-C_grt_rot_zrt2zxy_upar = libgrt.grt_rot_zrt2zxy_upar
-"""柱坐标下的位移偏导 ∂u(z,r,t)/∂(z,r,t) 转到 直角坐标 ∂u(z,x,y)/∂(z,x,y)"""
-C_grt_rot_zrt2zxy_upar.restype = None
-C_grt_rot_zrt2zxy_upar.argtypes = [c_double, DPOINTER, DPOINTER, c_double]  # double, double[3], double[3][3], double
-
-
-# -------------------------------------------------------------------
-#                      C函数定义的衰减函数
-# -------------------------------------------------------------------
-C_grt_attenuation_law = libgrt.grt_attenuation_law
-"""品质因子Q 对 波速的影响"""
-C_grt_attenuation_law.restype = CPLX
-C_grt_attenuation_law.argtypes = [REAL, CPLX, CPLX]  # double, cplx, cplx
-
-
-# -------------------------------------------------------------------
-#                      C 函数定义的方向因子
-# -------------------------------------------------------------------
-C_grt_set_source_radiation = libgrt.grt_set_source_radiation
-C_grt_set_source_radiation.restype = None
-C_grt_set_source_radiation.argtypes = [
-    (REAL*CHANNEL_NUM)*SRC_M_NUM, c_int, c_bool,
-    REAL, REAL, REAL, REAL, REAL*MECHANISM_NUM
-]
-
-
-
 # -------------------------------------------------------------------
 #                      使用 C 函数求解 Lamb 问题
 # -------------------------------------------------------------------
diff --git a/pygrt/c_structures.py b/pygrt/c_structures.py
index 158eda07..e1129ce6 100755
--- a/pygrt/c_structures.py
+++ b/pygrt/c_structures.py
@@ -24,6 +24,7 @@
     "qwvchs",
     "MECHANISM_NUM",
     "MIN_DEPTH_GAP_SRC_RCV",
+    "ZERO_DISTANCE",
 
     "NPCT_REAL_TYPE",
     "NPCT_CMPLX_TYPE",
@@ -52,6 +53,7 @@
 qwvchs = ['q', 'w', 'v']
 MECHANISM_NUM = 6
 MIN_DEPTH_GAP_SRC_RCV = 0.1
+ZERO_DISTANCE = 1e-8  # 与 C 侧 GRT_ZERO_DISTANCE 一致
 
 NPCT_REAL_TYPE = 'f8'
 NPCT_CMPLX_TYPE = 'c16'
diff --git a/pygrt/pymod.py b/pygrt/pymod.py
index accb5f29..cc3394fa 100755
--- a/pygrt/pymod.py
+++ b/pygrt/pymod.py
@@ -219,6 +219,10 @@ def _get_grn_spectra(
         # 只能设置一种filon积分方法
         if safilonTol > 0.0 and filonLength > 0.0:
             raise ValueError(f"You should only set one of filonLength and safilonTol.")
+
+        # FIM/SAFIM 面向远震中距，不能用于 r=0（含 filonCut 分段）
+        if (filonLength > 0.0 or safilonTol > 0.0) and np.any(np.asarray(distarr) <= ZERO_DISTANCE):
+            raise ValueError("FIM/SAFIM cannot be used with zero epicentral distance.")
         
         nf = nt//2+1 
         df = 1/(nt*dt)
@@ -246,13 +250,10 @@ def _get_grn_spectra(
         # 虚频率 
         wI = zeta * np.pi/(nt*dt)
 
-        # 避免绝对0震中距 
         nrs = len(distarr)
         for ir in range(nrs):
             if(distarr[ir] < 0.0):
                 raise ValueError(f"r({distarr[ir]}) < 0")
-            elif(distarr[ir] == 0.0):
-                distarr[ir] = 1e-5 
 
         # 最大震中距
         rmax = np.max(distarr)
@@ -505,8 +506,8 @@ def compute_grn(
                                      or source and receiver are at the same depth, DCM is applied in Auto mode.
             :param    use_kmax_ref:   directly use kmax_ref as kmax, without amplitude search
             :param    Length:        integration step `dk=2\pi / (L*rmax)`, see Bouchon (1981) and 张海明 (2021) for the criterion, default set automatically.
-            :param    filonLength:   integration step of Fixed-Interval Filon's Integration Method
-            :param    safilonTol:    precision of Self-Adaptive Filon's Integration Method
+            :param    filonLength:   integration step of Fixed-Interval Filon's Integration Method (large distance only; not for r=0)
+            :param    safilonTol:    precision of Self-Adaptive Filon's Integration Method (large distance only; not for r=0)
             :param    filonCut:      The splitting point of DWM and (SA)FIM, k*=<filonCut>/rmax, default is 0
             :param    converg_method:   The method of explicit convergence, you can set "AUTO", "NONE", "DCM" or "PTAM". Default use "AUTO".
             :param    skipImagComps:    skip the amplitude compensation from imaginary frequency.
@@ -578,8 +579,8 @@ def compute_static_grn(
                                         or source and receiver are at the same depth, DCM is applied in Auto mode.
             :param       use_kmax_ref:   directly use kmax_ref as kmax, without amplitude search
             :param       Length:        integration step `dk=2\pi / (L*rmax)`, default L=15
-            :param       filonLength:   integration step of Fixed-Interval Filon's Integration Method
-            :param       safilonTol:    precision of Self-Adaptive Filon's Integration Method
+            :param       filonLength:   integration step of Fixed-Interval Filon's Integration Method (large distance only; not for r=0)
+            :param       safilonTol:    precision of Self-Adaptive Filon's Integration Method (large distance only; not for r=0)
             :param       filonCut:      The splitting point of DWM and (SA)FIM, k*=<filonCut>/rmax, default is 0
             :param    converg_method:   The method of explicit convergence, you can set "AUTO", "NONE", "DCM" or "PTAM". Default use "AUTO".
             :param       calc_upar:     whether calculate the spatial derivatives of displacements.
@@ -638,8 +639,12 @@ def compute_static_grn(
         rs = np.zeros((nr,), dtype=NPCT_REAL_TYPE)
         for iy in range(ny):
             for ix in range(nx):
-                rs[ix + iy*nx] = max(np.hypot(xarr[ix], yarr[iy]), 1e-5)
+                rs[ix + iy*nx] = np.hypot(xarr[ix], yarr[iy])
         c_rs = npct.as_ctypes(rs)
+
+        # FIM/SAFIM 面向远震中距，不能用于 r=0（含 filonCut 分段）
+        if (filonLength > 0.0 or safilonTol > 0.0) and np.any(rs <= ZERO_DISTANCE):
+            raise ValueError("FIM/SAFIM cannot be used with zero epicentral distance.")
         
         # 设置波数积分间隔
         if Length == 0.0:
diff --git a/pygrt/utils.py b/pygrt/utils.py
index 74fd4377..1679a415 100755
--- a/pygrt/utils.py
+++ b/pygrt/utils.py
@@ -6,9 +6,8 @@
     该文件包含一些数据处理操作上的补充:   
         1、剪切源、张裂源、单力源、爆炸源、矩张量源 通过格林函数合成理论地震图的函数\n
         2、Stream类型的时域卷积、微分、积分 (基于numpy和scipy)    \n
-        3、Stream类型写到本地sac文件，自定义名称    \n
-        4、读取波数积分和峰谷平均法过程文件  \n
-        5、其它辅助函数  \n
+        3、读取波数积分和峰谷平均法过程文件  \n
+        4、其它辅助函数  \n
 
 """
 
@@ -97,14 +96,10 @@ def _gen_syn_from_gf(st:Stream, calc_upar:bool, compute_type:GRT_SYN_TYPE, M0:fl
     src_vb = st[0].stats.sac['user7']
     VpVs_ratio = float(src_va / src_vb)
 
-    azrad = float(np.deg2rad(az))
+    azrad = REAL(np.deg2rad(az))
     dist = float(st[0].stats.sac['dist'])
     npts = int(st[0].stats.npts)
 
-    baz = 180 + az
-    if baz > 360:
-        baz -= 360
-
     FPtrs = FPOINTER * CHANNEL_NUM
     FGrid = FPtrs * SRC_M_NUM
     FMat = FPtrs * CHANNEL_NUM
@@ -165,6 +160,12 @@ def _src_chnl_ptrs(hold):
     )
     # ========================================================================
 
+    # C 可能因 r=0 强制 azrad=0，同步方位角头段
+    az = float(np.rad2deg(azrad.value))
+    baz = az + 180.0
+    if baz >= 360.0:
+        baz -= 360.0
+
     out_chs = ZNEchs if ZNE else ZRTchs
     stall = Stream()
     for c, ch in enumerate(out_chs):
@@ -279,35 +280,6 @@ def _set_source_mechanism(
     return mchn
 
 
-def _set_source_radi(
-    par_theta:bool, coef:float, compute_type:GRT_SYN_TYPE, M0:float, azrad:float, **kwargs):
-    r"""
-        使用 C 函数计算不同震源的方向因子矩阵
-
-        :param    par_theta:       是否求对theta的偏导
-        :param    coef:            比例系数
-        :param    compute_type:    计算震源类型
-        :param    M0:              地震矩
-        :param    azrad:           方位角(弧度)
-
-        - 其他参数根据计算类型可选:
-            - 单力源需要: fZ, fN, fE, 
-            - 剪切源需要: strike, dip, rake
-            - 矩张量源需要: MT=(Mxx, Mxy, Mxz, Myy, Myz, Mzz)
-    """
-
-    VpVs_ratio = kwargs['VpVs_ratio'] if 'VpVs_ratio' in kwargs else 0.0
-
-    src_coef = np.zeros((SRC_M_NUM, CHANNEL_NUM), dtype=NPCT_REAL_TYPE)
-    mchn = _set_source_mechanism(compute_type, **kwargs)
-
-    C_grt_set_source_radiation(
-        npct.as_ctypes(src_coef), compute_type.value, par_theta, 
-        M0, coef, VpVs_ratio, azrad, npct.as_ctypes(mchn))
-    
-    return src_coef
-
-
 def gen_syn_from_gf_DC(st:Union[Stream,dict], M0:float, strike:float, dip:float, rake:float, az:float=-999, ZNE=False, calc_upar:bool=False, **kwargs):
     '''
         Shear source, the unit of angles is all degrees(°)
@@ -498,6 +470,9 @@ def _compute_strain_rotation(st_syn:Stream, Type:str):
 
 def _prepare_dynamic_postprocess_arrays(st_syn:Stream, chs:List[str], npts:int):
     """收集动态位移/偏导数组，并构造 ctypes 通道指针表。"""
+    FPtrs = FPOINTER * CHANNEL_NUM
+    FMat = FPtrs * CHANNEL_NUM
+
     def data(channel:str):
         st = st_syn.select(channel=channel)
         if len(st) == 0:
@@ -508,23 +483,29 @@ def data(channel:str):
 
     u = [data(c) for c in chs]
     upar = [[data(f"{d.lower()}{c}") for c in chs] for d in chs]
-    u_ptrs = FLOAT_CHNL_PTRS(*(arr.ctypes.data_as(FPOINTER) for arr in u))
-    upar_ptrs = FLOAT_CHNL_PTR_MAT(*(
-        FLOAT_CHNL_PTRS(*(arr.ctypes.data_as(FPOINTER) for arr in row)) for row in upar))
+    u_ptrs = FPtrs(*(arr.ctypes.data_as(FPOINTER) for arr in u))
+    upar_ptrs = FMat(*(
+        FPtrs(*(arr.ctypes.data_as(FPOINTER) for arr in row)) for row in upar))
     return u, upar, u_ptrs, upar_ptrs
 
 
 def _prepare_dynamic_postprocess_result(npts:int):
     """分配动态后处理结果数组及其 ctypes 通道指针表。"""
+    FPtrs = FPOINTER * CHANNEL_NUM
+    FMat = FPtrs * CHANNEL_NUM
+
     resarr = np.zeros((CHANNEL_NUM, CHANNEL_NUM, npts), dtype=np.float32)
-    res_ptrs = FLOAT_CHNL_PTR_MAT(*(
-        FLOAT_CHNL_PTRS(*(resarr[c2, c1].ctypes.data_as(FPOINTER) for c1 in range(CHANNEL_NUM)))
+    res_ptrs = FMat(*(
+        FPtrs(*(resarr[c2, c1].ctypes.data_as(FPOINTER) for c1 in range(CHANNEL_NUM)))
         for c2 in range(CHANNEL_NUM)))
     return resarr, res_ptrs
 
 
 def _prepare_static_postprocess_arrays(syn:dict, chs:List[str]):
     """整理静态后处理所需的连续内存数组与 ctypes 通道指针表。"""
+    RPtrs = PREAL * CHANNEL_NUM
+    RMat = RPtrs * CHANNEL_NUM
+
     xarr = np.ascontiguousarray(syn['_xarr'], dtype=np.float64)
     yarr = np.ascontiguousarray(syn['_yarr'], dtype=np.float64)
     if xarr.ndim != 1 or yarr.ndim != 1:
@@ -540,12 +521,12 @@ def _prepare_static_postprocess_arrays(syn:dict, chs:List[str]):
             arr.shape != expected_shape for row in upar for arr in row):
         raise ValueError("Static displacement and derivative arrays must match '_xarr'/'_yarr'.")
 
-    u_ptrs = REAL_CHNL_PTRS(*(arr.ctypes.data_as(PREAL) for arr in u))
-    upar_ptrs = REAL_CHNL_PTR_MAT(*(
-        REAL_CHNL_PTRS(*(arr.ctypes.data_as(PREAL) for arr in row)) for row in upar))
+    u_ptrs = RPtrs(*(arr.ctypes.data_as(PREAL) for arr in u))
+    upar_ptrs = RMat(*(
+        RPtrs(*(arr.ctypes.data_as(PREAL) for arr in row)) for row in upar))
     resarr = np.zeros((CHANNEL_NUM, CHANNEL_NUM, *expected_shape), dtype=np.float64)
-    res_ptrs = REAL_CHNL_PTR_MAT(*(
-        REAL_CHNL_PTRS(*(resarr[c2, c1].ctypes.data_as(PREAL) for c1 in range(CHANNEL_NUM)))
+    res_ptrs = RMat(*(
+        RPtrs(*(resarr[c2, c1].ctypes.data_as(PREAL) for c1 in range(CHANNEL_NUM)))
         for c2 in range(CHANNEL_NUM)))
     return xarr, yarr, u, upar, u_ptrs, upar_ptrs, resarr, res_ptrs