FEAT: support spatial derivatives of displacements

pygrt/C_extension/include/common/bessel.h

@@ -19,4 +19,16 @@
  * @param bj2  (out) \f$ J_2(x) \f$
  * 
  */
-void bessel012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2);
+void bessel012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2);
+
+
+/**
+ * 计算Bessel函数的一阶导数 \f$ J_m^'(x), m=0,1,2 \f$ 
+ * 
+ * @param x          自变量 
+ * @param bj0  (inout) 传入\f$ J_0(x) \f$, 返回\f$ J_0^'(x) \f$
+ * @param bj1  (inout) 传入\f$ J_1(x) \f$, 返回\f$ J_1^'(x) \f$
+ * @param bj2  (inout) 传入\f$ J_2(x) \f$, 返回\f$ J_2^'(x) \f$
+ * 
+ */
+void besselp012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2);

pygrt/C_extension/include/common/const.h

@@ -68,6 +68,9 @@ typedef int MYINT;  ///< 整数
     #define PI2 6.2831853f        ///< \f$ 2\pi \f$
     #define HALFPI 1.5707963f   ///< \f$ \frac{\pi}{2} \f$
     #define QUARTERPI 0.78539816f   ///< \f$ \frac{\pi}{4} \f$
+    #define THREEQUARTERPI 2.35619450f  ///< \f$ \frac{3\pi}{4} \f$
+    #define FIVEQUARTERPI  3.92699082f  ///< \f$ \frac{5\pi}{4} \f$
+    #define SEVENQUARTERPI  5.49778714f   ///< \f$ \frac{7\pi}{4} \f$
     #define INV_SQRT_TWO 0.70710678f   ///< \f$ \frac{1}{\sqrt{2}} \f$ 
 
 #else 
@@ -102,6 +105,9 @@ typedef int MYINT;  ///< 整数
     #define PI2 6.283185307179586        ///< \f$ 2\pi \f$
     #define HALFPI 1.5707963267948966   ///< \f$ \frac{\pi}{2} \f$
     #define QUARTERPI 0.7853981633974483   ///< \f$ \frac{\pi}{4} \f$
+    #define THREEQUARTERPI 2.356194490192345  ///< \f$ \frac{3\pi}{4} \f$
+    #define FIVEQUARTERPI  3.9269908169872414  ///< \f$ \frac{5\pi}{4} \f$
+    #define SEVENQUARTERPI  5.497787143782138   ///< \f$ \frac{7\pi}{4} \f$
     #define INV_SQRT_TWO 0.7071067811865475   ///< \f$ \frac{1}{\sqrt{2}} \f$ 
 
 #endif

pygrt/C_extension/include/dynamic/dwm.h

@@ -32,6 +32,17 @@
  * @param  sum_VF_J[nr][3][4]   (out)垂直力源
  * @param  sum_HF_J[nr][3][4]   (out)水平力源
  * @param  sum_DC_J[nr][3][4]   (out)双力偶源
+ * 
+ * @param  calc_upar       (in)是否计算位移u的空间导数
+ * @param  sum_EXP_uiz_J[nr][3][4]  (out)爆炸源产生的ui_z(位移u对坐标z的偏导)，下同
+ * @param  sum_VF_uiz_J[nr][3][4]   (out)垂直力源
+ * @param  sum_HF_uiz_J[nr][3][4]   (out)水平力源
+ * @param  sum_DC_uiz_J[nr][3][4]   (out)双力偶源
+ * @param  sum_EXP_uir_J[nr][3][4]  (out)爆炸源产生的ui_r(位移u对坐标r的偏导)，下同
+ * @param  sum_VF_uir_J[nr][3][4]   (out)垂直力源
+ * @param  sum_HF_uir_J[nr][3][4]   (out)水平力源
+ * @param  sum_DC_uir_J[nr][3][4]   (out)双力偶源
+ * 
  * @param  fstats[nr])          (out)不同震中距的格林函数积分过程文件
  * 
  * @return  k        积分截至时的波数
@@ -41,4 +52,9 @@ MYREAL discrete_integ(
     MYINT nr, MYREAL *rs,
     MYCOMPLEX sum_EXP_J[nr][3][4], MYCOMPLEX sum_VF_J[nr][3][4],  
     MYCOMPLEX sum_HF_J[nr][3][4],  MYCOMPLEX sum_DC_J[nr][3][4],  
+    bool calc_upar,
+    MYCOMPLEX sum_EXP_uiz_J[nr][3][4], MYCOMPLEX sum_VF_uiz_J[nr][3][4],  
+    MYCOMPLEX sum_HF_uiz_J[nr][3][4],  MYCOMPLEX sum_DC_uiz_J[nr][3][4],  
+    MYCOMPLEX sum_EXP_uir_J[nr][3][4], MYCOMPLEX sum_VF_uir_J[nr][3][4],  
+    MYCOMPLEX sum_HF_uir_J[nr][3][4],  MYCOMPLEX sum_DC_uir_J[nr][3][4],  
     FILE *(fstats[nr]));

pygrt/C_extension/include/dynamic/fim.h

@@ -38,6 +38,17 @@
  * @param  sum_VF_J[nr][3][4]   (out)垂直力源
  * @param  sum_HF_J[nr][3][4]   (out)水平力源
  * @param  sum_DC_J[nr][3][4]   (out)双力偶源
+ * 
+ * @param  calc_upar       (in)是否计算位移u的空间导数
+ * @param  sum_EXP_uiz_J[nr][3][4]  (out)爆炸源产生的ui_z(位移u对坐标z的偏导)，下同
+ * @param  sum_VF_uiz_J[nr][3][4]   (out)垂直力源
+ * @param  sum_HF_uiz_J[nr][3][4]   (out)水平力源
+ * @param  sum_DC_uiz_J[nr][3][4]   (out)双力偶源
+ * @param  sum_EXP_uir_J[nr][3][4]  (out)爆炸源产生的ui_r(位移u对坐标r的偏导)，下同
+ * @param  sum_VF_uir_J[nr][3][4]   (out)垂直力源
+ * @param  sum_HF_uir_J[nr][3][4]   (out)水平力源
+ * @param  sum_DC_uir_J[nr][3][4]   (out)双力偶源
+ * 
  * @param  fstats[nr]           (out)不同震中距的格林函数积分过程文件
  * 
  * @return  k        积分截至时的波数
@@ -47,6 +58,11 @@ MYREAL linear_filon_integ(
     MYINT nr, MYREAL *rs,
     MYCOMPLEX sum_EXP_J[nr][3][4], MYCOMPLEX sum_VF_J[nr][3][4],  
     MYCOMPLEX sum_HF_J[nr][3][4],  MYCOMPLEX sum_DC_J[nr][3][4],  
+    bool calc_upar,
+    MYCOMPLEX sum_EXP_uiz_J[nr][3][4], MYCOMPLEX sum_VF_uiz_J[nr][3][4],  
+    MYCOMPLEX sum_HF_uiz_J[nr][3][4],  MYCOMPLEX sum_DC_uiz_J[nr][3][4],  
+    MYCOMPLEX sum_EXP_uir_J[nr][3][4], MYCOMPLEX sum_VF_uir_J[nr][3][4],  
+    MYCOMPLEX sum_HF_uir_J[nr][3][4],  MYCOMPLEX sum_DC_uir_J[nr][3][4],  
     FILE *(fstats[nr]));
 
 
@@ -63,6 +79,7 @@ MYREAL linear_filon_integ(
  * @param    VF_qwv[3][3]     (in)垂直力源核函数
  * @param    HF_qwv[3][3]     (in)水平力源核函数
  * @param    DC_qwv[3][3]     (in)双力偶源核函数 
+ * @param    calc_uir         (in)是否计算ui_r（位移u对坐标r的偏导）
  * @param    EXP_J[3][4]      (out)爆炸源，该dk区间内的积分值，下同
  * @param    VF_J[3][4]       (out)垂直力源
  * @param    HF_J[3][4]       (out)水平力源
@@ -73,6 +90,7 @@ void int_Pk_filon(
     MYREAL k, MYREAL r, 
     const MYCOMPLEX EXP_qwv[3][3], const MYCOMPLEX VF_qwv[3][3], 
     const MYCOMPLEX HF_qwv[3][3],  const MYCOMPLEX DC_qwv[3][3], 
+    bool calc_uir, 
     MYCOMPLEX EXP_J[3][4], MYCOMPLEX VF_J[3][4], 
     MYCOMPLEX HF_J[3][4],  MYCOMPLEX DC_J[3][4] );
 

pygrt/C_extension/include/dynamic/grt.h

@@ -49,6 +49,21 @@ void set_num_threads(int num_threads);
  * @param      DSgrn[nr][3]       (out)复数数组，90度倾滑的Z、R、T分量频谱结果
  * @param      SSgrn[nr][3]       (out)复数数组，90度走滑的Z、R、T分量频谱结果
  * 
+ * @param      calc_upar              (in)是否计算位移u的空间导数
+ * @param      EXPgrn_uiz[nr][2]      (out)复数数组，爆炸源产生的ui_z(位移u对坐标z的偏导)的Z、R分量频谱结果，下同
+ * @param      VFgrn_uiz[nr][2]       (out)复数数组，垂直力源的Z、R分量频谱结果
+ * @param      HFgrn_uiz[nr][3]       (out)复数数组，水平力源的Z、R、T分量频谱结果
+ * @param      DDgrn_uiz[nr][2]       (out)复数数组，45度倾滑的Z、R分量频谱结果
+ * @param      DSgrn_uiz[nr][3]       (out)复数数组，90度倾滑的Z、R、T分量频谱结果
+ * @param      SSgrn_uiz[nr][3]       (out)复数数组，90度走滑的Z、R、T分量频谱结果
+ * @param      EXPgrn_uir[nr][2]      (out)复数数组，爆炸源产生的ui_r(位移u对坐标r的偏导)的Z、R分量频谱结果，下同
+ * @param      VFgrn_uir[nr][2]       (out)复数数组，垂直力源的Z、R分量频谱结果
+ * @param      HFgrn_uir[nr][3]       (out)复数数组，水平力源的Z、R、T分量频谱结果
+ * @param      DDgrn_uir[nr][2]       (out)复数数组，45度倾滑的Z、R分量频谱结果
+ * @param      DSgrn_uir[nr][3]       (out)复数数组，90度倾滑的Z、R、T分量频谱结果
+ * @param      SSgrn_uir[nr][3]       (out)复数数组，90度走滑的Z、R、T分量频谱结果
+ * 
+ * 
  * @param      statsstr          (in) 积分结果输出目录
  * @param      nstatsidxs        (in) 输出积分结果的特定频点的个数
  * @param      statsidxs         (in) 特定频点的索引值
@@ -68,6 +83,20 @@ void integ_grn_spec_in_C(
     MYCOMPLEX *DScplx[nr][3],  // DSZ, DSR, DST 的实部和虚部 [DS: 90-dip slip]
     MYCOMPLEX *SScplx[nr][3],  // SSZ, SSR, SST 的实部和虚部 [SS: strike slip]
 
+    bool calc_upar,
+    MYCOMPLEX *EXPcplx_uiz[nr][2], // EXZ, EXR 的实部和虚部
+    MYCOMPLEX *VFcplx_uiz[nr][2],  // VFZ, VFR 的实部和虚部
+    MYCOMPLEX *HFcplx_uiz[nr][3],  // HFZ, HFR, HFT 的实部和虚部
+    MYCOMPLEX *DDcplx_uiz[nr][2],  // DDZ, DDR 的实部和虚部      [DD: 45-dip slip]
+    MYCOMPLEX *DScplx_uiz[nr][3],  // DSZ, DSR, DST 的实部和虚部 [DS: 90-dip slip]
+    MYCOMPLEX *SScplx_uiz[nr][3],  // SSZ, SSR, SST 的实部和虚部 [SS: strike slip]
+    MYCOMPLEX *EXPcplx_uir[nr][2], // EXZ, EXR 的实部和虚部
+    MYCOMPLEX *VFcplx_uir[nr][2],  // VFZ, VFR 的实部和虚部
+    MYCOMPLEX *HFcplx_uir[nr][3],  // HFZ, HFR, HFT 的实部和虚部
+    MYCOMPLEX *DDcplx_uir[nr][2],  // DDZ, DDR 的实部和虚部      [DD: 45-dip slip]
+    MYCOMPLEX *DScplx_uir[nr][3],  // DSZ, DSR, DST 的实部和虚部 [DS: 90-dip slip]
+    MYCOMPLEX *SScplx_uir[nr][3],  // SSZ, SSR, SST 的实部和虚部 [SS: strike slip]
+
     const char *statsstr, // 积分结果输出
     MYINT  nstatsidxs, // 仅输出特定频点
     MYINT *statsidxs
@@ -99,6 +128,20 @@ void integ_grn_spec_in_C(
  * @param      DSgrn[nr][3]       (out)`GRN` 结构体指针，90度倾滑的Z、R、T分量频谱结果
  * @param      SSgrn[nr][3]       (out)`GRN` 结构体指针，90度走滑的Z、R、T分量频谱结果
  * 
+ * @param      calc_upar              (in)是否计算位移u的空间导数
+ * @param      EXPgrn_uiz[nr][2]      (out)`GRN` 结构体指针，爆炸源产生的ui_z(位移u对坐标z的偏导)的Z、R分量频谱结果，下同
+ * @param      VFgrn_uiz[nr][2]       (out)`GRN` 结构体指针，垂直力源的Z、R分量频谱结果
+ * @param      HFgrn_uiz[nr][3]       (out)`GRN` 结构体指针，水平力源的Z、R、T分量频谱结果
+ * @param      DDgrn_uiz[nr][2]       (out)`GRN` 结构体指针，45度倾滑的Z、R分量频谱结果
+ * @param      DSgrn_uiz[nr][3]       (out)`GRN` 结构体指针，90度倾滑的Z、R、T分量频谱结果
+ * @param      SSgrn_uiz[nr][3]       (out)`GRN` 结构体指针，90度走滑的Z、R、T分量频谱结果
+ * @param      EXPgrn_uir[nr][2]      (out)`GRN` 结构体指针，爆炸源产生的ui_r(位移u对坐标r的偏导)的Z、R分量频谱结果，下同
+ * @param      VFgrn_uir[nr][2]       (out)`GRN` 结构体指针，垂直力源的Z、R分量频谱结果
+ * @param      HFgrn_uir[nr][3]       (out)`GRN` 结构体指针，水平力源的Z、R、T分量频谱结果
+ * @param      DDgrn_uir[nr][2]       (out)`GRN` 结构体指针，45度倾滑的Z、R分量频谱结果
+ * @param      DSgrn_uir[nr][3]       (out)`GRN` 结构体指针，90度倾滑的Z、R、T分量频谱结果
+ * @param      SSgrn_uir[nr][3]       (out)`GRN` 结构体指针，90度走滑的Z、R、T分量频谱结果
+ * 
  * @param      statsstr          (in) 积分结果输出目录
  * @param      nstatsidxs        (in) 输出积分结果的特定频点的个数
  * @param      statsidxs         (in) 特定频点的索引值
@@ -118,6 +161,20 @@ void integ_grn_spec(
     GRN *DSgrn[nr][3],  // DSZ, DSR, DST 的实部和虚部 [DS: 90-dip slip]
     GRN *SSgrn[nr][3],  // SSZ, SSR, SST 的实部和虚部 [SS: strike slip]
 
+    bool calc_upar,
+    GRN *EXPgrn_uiz[nr][2], // EXZ, EXR 的实部和虚部
+    GRN *VFgrn_uiz[nr][2],  // VFZ, VFR 的实部和虚部
+    GRN *HFgrn_uiz[nr][3],  // HFZ, HFR, HFT 的实部和虚部
+    GRN *DDgrn_uiz[nr][2],  // DDZ, DDR 的实部和虚部      [DD: 45-dip slip]
+    GRN *DSgrn_uiz[nr][3],  // DSZ, DSR, DST 的实部和虚部 [DS: 90-dip slip]
+    GRN *SSgrn_uiz[nr][3],  // SSZ, SSR, SST 的实部和虚部 [SS: strike slip]
+    GRN *EXPgrn_uir[nr][2], // EXZ, EXR 的实部和虚部
+    GRN *VFgrn_uir[nr][2],  // VFZ, VFR 的实部和虚部
+    GRN *HFgrn_uir[nr][3],  // HFZ, HFR, HFT 的实部和虚部
+    GRN *DDgrn_uir[nr][2],  // DDZ, DDR 的实部和虚部      [DD: 45-dip slip]
+    GRN *DSgrn_uir[nr][3],  // DSZ, DSR, DST 的实部和虚部 [DS: 90-dip slip]
+    GRN *SSgrn_uir[nr][3],  // SSZ, SSR, SST 的实部和虚部 [SS: strike slip]
+
     const char *statsstr, // 积分结果输出
     MYINT  nstatsidxs, // 仅输出特定频点
     MYINT *statsidxs

pygrt/C_extension/include/dynamic/layer.h

@@ -50,6 +50,13 @@ void calc_R_EV(
     MYCOMPLEX R_EV[2][2], MYCOMPLEX *R_EVL);
 
 
+void calc_uiz_R_EV(
+    MYCOMPLEX xa_rcv, MYCOMPLEX xb_rcv, bool ircvup,
+    MYREAL k, 
+    const MYCOMPLEX R[2][2], MYCOMPLEX RL, 
+    MYCOMPLEX R_EV[2][2], MYCOMPLEX *R_EVL);
+
+
 /**
  * 计算界面的反射系数RD/RDL/RU/RUL, 透射系数TD/TDL/TU/TUL, 包括时间延迟因子，
  * 后缀L表示SH波的系数, 其余表示P-SV波的系数, 根据公式(5.4.14)和(5.4.31)计算系数   

pygrt/C_extension/include/dynamic/propagate.h

@@ -88,11 +88,18 @@
  *  @param    VF_qwv[3][3]     (out)垂直力源核函数
  *  @param    HF_qwv[3][3]     (out)水平力源核函数
  *  @param    DC_qwv[3][3]     (out)双力偶源核函数 
+ *  @param    calc_uiz       (in)是否计算ui_z（位移u对坐标z的偏导）
+ *  @param    EXP_uiz_qwv[3][3]    (out)爆炸源产生的ui_z的核函数，下同
+ *  @param    VF_uiz_qwv[3][3]     (out)垂直力源核函数
+ *  @param    HF_uiz_qwv[3][3]     (out)水平力源核函数
+ *  @param    DC_uiz_qwv[3][3]     (out)双力偶源核函数 
  * 
  */
 void kernel(
     const MODEL1D *mod1d, MYCOMPLEX omega, MYREAL k,
-    MYCOMPLEX EXP_qwv[3][3], MYCOMPLEX VF_qwv[3][3], MYCOMPLEX HF_qwv[3][3], MYCOMPLEX DC_qwv[3][3]);
+    MYCOMPLEX EXP_qwv[3][3], MYCOMPLEX VF_qwv[3][3], MYCOMPLEX HF_qwv[3][3], MYCOMPLEX DC_qwv[3][3],
+    bool calc_uiz,
+    MYCOMPLEX EXP_uiz_qwv[3][3], MYCOMPLEX VF_uiz_qwv[3][3], MYCOMPLEX HF_uiz_qwv[3][3], MYCOMPLEX DC_uiz_qwv[3][3]);
 
 
 
@@ -108,6 +115,7 @@ void kernel(
  * @param    VF_qwv[3][3]     (in)垂直力源核函数
  * @param    HF_qwv[3][3]     (in)水平力源核函数
  * @param    DC_qwv[3][3]     (in)双力偶源核函数
+ * @param    calc_uir         (in)是否计算ui_r（位移u对坐标r的偏导）
  * @param    EXP_J[3][4]      (out)爆炸源，该dk区间内的积分值，下同
  * @param    VF_J[3][4]       (out)垂直力源
  * @param    HF_J[3][4]       (out)水平力源
@@ -118,6 +126,7 @@ void int_Pk(
     MYREAL k, MYREAL r, 
     const MYCOMPLEX EXP_qwv[3][3], const MYCOMPLEX VF_qwv[3][3], 
     const MYCOMPLEX HF_qwv[3][3],  const MYCOMPLEX DC_qwv[3][3], 
+    bool calc_uir,
     MYCOMPLEX EXP_J[3][4], MYCOMPLEX VF_J[3][4], 
     MYCOMPLEX HF_J[3][4],  MYCOMPLEX DC_J[3][4] );
 

pygrt/C_extension/include/dynamic/ptam.h

@@ -40,6 +40,17 @@
  * @param    sum_VF_J0[nr][3][4]    (out)垂直力源
  * @param    sum_HF_J0[nr][3][4]    (out)水平力源
  * @param    sum_DC_J0[nr][3][4]    (out)双力偶源
+ * 
+ * @param  calc_upar       (in)是否计算位移u的空间导数
+ * @param  sum_EXP_uiz_J0[nr][3][4]  (out)爆炸源产生的ui_z(位移u对坐标z的偏导)，下同
+ * @param  sum_VF_uiz_J0[nr][3][4]   (out)垂直力源
+ * @param  sum_HF_uiz_J0[nr][3][4]   (out)水平力源
+ * @param  sum_DC_uiz_J0[nr][3][4]   (out)双力偶源
+ * @param  sum_EXP_uir_J0[nr][3][4]  (out)爆炸源产生的ui_r(位移u对坐标r的偏导)，下同
+ * @param  sum_VF_uir_J0[nr][3][4]   (out)垂直力源
+ * @param  sum_HF_uir_J0[nr][3][4]   (out)水平力源
+ * @param  sum_DC_uir_J0[nr][3][4]   (out)双力偶源
+ * 
  * @param    fstats[nr]          (out)波数积分过程文件指针
  * @param    ptam_fstats[nr]     (out)峰谷平均法过程文件指针
  * 
@@ -50,6 +61,11 @@ void PTA_method(
     MYINT nr, MYREAL *rs, 
     MYCOMPLEX sum_EXP_J0[nr][3][4], MYCOMPLEX sum_VF_J0[nr][3][4],  
     MYCOMPLEX sum_HF_J0[nr][3][4],  MYCOMPLEX sum_DC_J0[nr][3][4],  
+    bool calc_upar,
+    MYCOMPLEX sum_EXP_uiz_J0[nr][3][4], MYCOMPLEX sum_VF_uiz_J0[nr][3][4],  
+    MYCOMPLEX sum_HF_uiz_J0[nr][3][4],  MYCOMPLEX sum_DC_uiz_J0[nr][3][4],  
+    MYCOMPLEX sum_EXP_uir_J0[nr][3][4], MYCOMPLEX sum_VF_uir_J0[nr][3][4],  
+    MYCOMPLEX sum_HF_uir_J0[nr][3][4],  MYCOMPLEX sum_DC_uir_J0[nr][3][4], 
     FILE *(fstats[nr]), FILE *(ptam_fstats[nr]));
 
 

pygrt/C_extension/src/common/bessel.c

@@ -16,3 +16,11 @@ void bessel012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2){
     *bj2 = JN(2, x);
 }
 
+void besselp012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2){
+    MYREAL j0=*bj0;
+    MYREAL j1=*bj1;
+    MYREAL j2=*bj2;
+    *bj0 = -j1;
+    *bj1 = j0 - RONE/x * j1;
+    *bj2 = j1 - RTWO/x * j2;
+}

pygrt/C_extension/src/dynamic/dwm.c

@@ -27,6 +27,11 @@ MYREAL discrete_integ(
     MYINT nr, MYREAL *rs,
     MYCOMPLEX sum_EXP_J[nr][3][4], MYCOMPLEX sum_VF_J[nr][3][4],  
     MYCOMPLEX sum_HF_J[nr][3][4],  MYCOMPLEX sum_DC_J[nr][3][4],  
+    bool calc_upar,
+    MYCOMPLEX sum_EXP_uiz_J[nr][3][4], MYCOMPLEX sum_VF_uiz_J[nr][3][4],  
+    MYCOMPLEX sum_HF_uiz_J[nr][3][4],  MYCOMPLEX sum_DC_uiz_J[nr][3][4],  
+    MYCOMPLEX sum_EXP_uir_J[nr][3][4], MYCOMPLEX sum_VF_uir_J[nr][3][4],  
+    MYCOMPLEX sum_HF_uir_J[nr][3][4],  MYCOMPLEX sum_DC_uir_J[nr][3][4],  
     FILE *(fstats[nr]))
 {
     MYCOMPLEX EXP_J[3][4], VF_J[3][4], HF_J[3][4],  DC_J[3][4];
@@ -40,6 +45,12 @@ MYREAL discrete_integ(
     MYCOMPLEX (*pVF_qwv)[3]  = (sum_VF_J!=NULL)?  VF_qwv  : NULL;
     MYCOMPLEX (*pHF_qwv)[3]  = (sum_HF_J!=NULL)?  HF_qwv  : NULL;
     MYCOMPLEX (*pDC_qwv)[3]  = (sum_DC_J!=NULL)?  DC_qwv  : NULL;
+
+    MYCOMPLEX EXP_uiz_qwv[3][3], VF_uiz_qwv[3][3], HF_uiz_qwv[3][3], DC_uiz_qwv[3][3]; 
+    MYCOMPLEX (*pEXP_uiz_qwv)[3] = (sum_EXP_uiz_J!=NULL)? EXP_uiz_qwv : NULL;
+    MYCOMPLEX (*pVF_uiz_qwv)[3]  = (sum_VF_uiz_J!=NULL)?  VF_uiz_qwv  : NULL;
+    MYCOMPLEX (*pHF_uiz_qwv)[3]  = (sum_HF_uiz_J!=NULL)?  HF_uiz_qwv  : NULL;
+    MYCOMPLEX (*pDC_uiz_qwv)[3]  = (sum_DC_uiz_J!=NULL)?  DC_uiz_qwv  : NULL;
 
     MYREAL k = 0.0;
     MYINT ik = 0;
@@ -61,7 +72,8 @@ MYREAL discrete_integ(
 
         // printf("w=%15.5e, ik=%d\n", CREAL(omega), ik);
         // 计算核函数 F(k, w)
-        kernel(mod1d, omega, k, pEXP_qwv, pVF_qwv, pHF_qwv, pDC_qwv); 
+        kernel(mod1d, omega, k, pEXP_qwv, pVF_qwv, pHF_qwv, pDC_qwv, 
+            calc_upar, pEXP_uiz_qwv, pVF_uiz_qwv, pHF_uiz_qwv, pDC_uiz_qwv); 
 
 
         // 震中距rs循环
@@ -77,8 +89,8 @@ MYREAL discrete_integ(
 
             // 计算被积函数一项 F(k,w)Jm(kr)k
             int_Pk(k, rs[ir], 
-                    pEXP_qwv, pVF_qwv, pHF_qwv, pDC_qwv,
-                    EXP_J, VF_J, HF_J, DC_J);
+                   pEXP_qwv, pVF_qwv, pHF_qwv, pDC_qwv, false,
+                   EXP_J, VF_J, HF_J, DC_J);
 
             // 记录积分结果
             if(fstats[ir]!=NULL){
@@ -114,7 +126,42 @@ MYREAL discrete_integ(
                 iendk = iendkrs[ir] = false;
             }
 
-
+
+            // ---------------- 位移空间导数，EXP_J, VF_J, HF_J, DC_J数组重复利用 --------------------------
+            if(calc_upar){
+                // ------------------------------- ui_z -----------------------------------
+                // 计算被积函数一项 F(k,w)Jm(kr)k
+                int_Pk(k, rs[ir], 
+                       pEXP_uiz_qwv, pVF_uiz_qwv, pHF_uiz_qwv, pDC_uiz_qwv, false,
+                       EXP_J, VF_J, HF_J, DC_J);
+
+                // keps不参与计算位移空间导数的积分，背后逻辑认为u收敛，则uiz也收敛
+                for(MYINT m=0; m<3; ++m){
+                    for(MYINT v=0; v<4; ++v){
+                        if(sum_EXP_uiz_J!=NULL) sum_EXP_uiz_J[ir][m][v] += EXP_J[m][v];
+                        if(sum_VF_uiz_J!=NULL)  sum_VF_uiz_J[ir][m][v]  += VF_J[m][v];
+                        if(sum_HF_uiz_J!=NULL)  sum_HF_uiz_J[ir][m][v]  += HF_J[m][v];
+                        if(sum_DC_uiz_J!=NULL)  sum_DC_uiz_J[ir][m][v]  += DC_J[m][v];
+                    }
+                }
+
+
+                // ------------------------------- ui_r -----------------------------------
+                // 计算被积函数一项 F(k,w)Jm(kr)k
+                int_Pk(k, rs[ir], 
+                       pEXP_qwv, pVF_qwv, pHF_qwv, pDC_qwv, true,
+                       EXP_J, VF_J, HF_J, DC_J);
+
+                // keps不参与计算位移空间导数的积分，背后逻辑认为u收敛，则uir也收敛
+                for(MYINT m=0; m<3; ++m){
+                    for(MYINT v=0; v<4; ++v){
+                        if(sum_EXP_uir_J!=NULL) sum_EXP_uir_J[ir][m][v] += EXP_J[m][v];
+                        if(sum_VF_uir_J!=NULL)  sum_VF_uir_J[ir][m][v]  += VF_J[m][v];
+                        if(sum_HF_uir_J!=NULL)  sum_HF_uir_J[ir][m][v]  += HF_J[m][v];
+                        if(sum_DC_uir_J!=NULL)  sum_DC_uir_J[ir][m][v]  += DC_J[m][v];
+                    }
+                }
+            } // END if calc_upar
 
         } // END rs loop