diff --git a/pygrt/C_extension/include/grt/common/RT_matrix.h b/pygrt/C_extension/include/grt/common/RT_matrix.h
index ed6478e..e5e46c7 100644
--- a/pygrt/C_extension/include/grt/common/RT_matrix.h
+++ b/pygrt/C_extension/include/grt/common/RT_matrix.h
@@ -39,6 +39,8 @@ typedef struct {
  * @param[out]   M    RT_MATRIX 指针
  */
 void grt_reset_RT_matrix(RT_MATRIX *M);
+void grt_reset_RT_matrix_PSV(RT_MATRIX *M);
+void grt_reset_RT_matrix_SH(RT_MATRIX *M);
 
 /** 
  * 合并 recursion_RD_PSV(SH) ，仅计算RD/RDL
diff --git a/pygrt/C_extension/include/grt/dynamic/layer.h b/pygrt/C_extension/include/grt/dynamic/layer.h
index 75af04b..f00aa7f 100755
--- a/pygrt/C_extension/include/grt/dynamic/layer.h
+++ b/pygrt/C_extension/include/grt/dynamic/layer.h
@@ -24,7 +24,8 @@
  * @param[in,out]     mod1d          模型结构体指针，结果保存在 M_top
  * 
  */
-void grt_topbound_RU(GRT_MODEL1D *mod1d);
+void grt_topbound_RU_PSV(GRT_MODEL1D *mod1d);
+void grt_topbound_RU_SH(GRT_MODEL1D *mod1d);
 
 /**
  * 计算不同边界条件下底界面的反射系数RD，其中自由表面的公式见(5.3.10-14) 
@@ -33,7 +34,8 @@ void grt_topbound_RU(GRT_MODEL1D *mod1d);
  * @param[in,out]     mod1d          模型结构体指针，结果保存在 M_bot
  * 
  */
-void grt_botbound_RD(GRT_MODEL1D *mod1d);
+void grt_botbound_RD_PSV(GRT_MODEL1D *mod1d);
+void grt_botbound_RD_SH(GRT_MODEL1D *mod1d);
 
 
 /**
@@ -125,6 +127,8 @@ void grt_RT_matrix_ss_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M
  * 
  */
 void grt_delay_RT_matrix(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M);
+void grt_delay_RT_matrix_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M);
+void grt_delay_RT_matrix_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M);
 
 /**
  * 为虚拟层的广义 R/T 矩阵添加时间延迟因子
diff --git a/pygrt/C_extension/include/grt/static/static_layer.h b/pygrt/C_extension/include/grt/static/static_layer.h
index 7be2e83..b1479e5 100644
--- a/pygrt/C_extension/include/grt/static/static_layer.h
+++ b/pygrt/C_extension/include/grt/static/static_layer.h
@@ -25,7 +25,8 @@
  * @param[in,out]     mod1d             模型结构体指针，结果保存在 M_top
  * 
  */
-void grt_static_topbound_RU(GRT_MODEL1D *mod1d);
+void grt_static_topbound_RU_PSV(GRT_MODEL1D *mod1d);
+void grt_static_topbound_RU_SH(GRT_MODEL1D *mod1d);
 
 /**
  * 计算不同边界条件下底界面的静态反射系数RD，其中自由表面的公式见(6.3.12)
@@ -33,7 +34,8 @@ void grt_static_topbound_RU(GRT_MODEL1D *mod1d);
  * @param[in,out]     mod1d             模型结构体指针，结果保存在 M_bot
  * 
  */
-void grt_static_botbound_RD(GRT_MODEL1D *mod1d);
+void grt_static_botbound_RD_PSV(GRT_MODEL1D *mod1d);
+void grt_static_botbound_RD_SH(GRT_MODEL1D *mod1d);
 
 /**
  * 计算接收点位置的 P-SV 静态接收矩阵，将波场转为位移，公式(6.3.35)
diff --git a/pygrt/C_extension/src/common/RT_matrix.c b/pygrt/C_extension/src/common/RT_matrix.c
index d0f7d23..a32bb6f 100644
--- a/pygrt/C_extension/src/common/RT_matrix.c
+++ b/pygrt/C_extension/src/common/RT_matrix.c
@@ -22,20 +22,26 @@
 
 void grt_reset_RT_matrix(RT_MATRIX *M)
 {
-    RT_MATRIX *M0 = &(RT_MATRIX){
-        .RD = GRT_INIT_ZERO_2x2_MATRIX,
-        .RU = GRT_INIT_ZERO_2x2_MATRIX,
-        .TD = GRT_INIT_IDENTITY_2x2_MATRIX,
-        .TU = GRT_INIT_IDENTITY_2x2_MATRIX,
-        .RDL = 0.0,
-        .RUL = 0.0,
-        .TDL = 1.0,
-        .TUL = 1.0,
-        .invT = GRT_INIT_ZERO_2x2_MATRIX,
-        .invTL = 0.0,
-        .stats = GRT_INVERSE_SUCCESS 
-    };
-    memcpy(M, M0, sizeof(*M));
+    memset(M, 0, sizeof(*M));
+    M->TD[0][0] = M->TD[1][1] = 
+    M->TU[0][0] = M->TU[1][1] = 
+    M->TDL = M->TUL = 1.0;
+    M->stats = GRT_INVERSE_SUCCESS;
+}
+
+void grt_reset_RT_matrix_PSV(RT_MATRIX *M)
+{
+    memset(M, 0, sizeof(*M));
+    M->TD[0][0] = M->TD[1][1] = 
+    M->TU[0][0] = M->TU[1][1] = 1.0;
+    M->stats = GRT_INVERSE_SUCCESS;
+}
+
+void grt_reset_RT_matrix_SH(RT_MATRIX *M)
+{
+    M->RDL = M->RUL = M->invTL = 0.0;
+    M->TDL = M->TUL = 1.0;
+    M->stats = GRT_INVERSE_SUCCESS;
 }
 
 void grt_recursion_RD(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MATRIX *M)
diff --git a/pygrt/C_extension/src/common/model.c b/pygrt/C_extension/src/common/model.c
index 1acc3c9..e80983d 100755
--- a/pygrt/C_extension/src/common/model.c
+++ b/pygrt/C_extension/src/common/model.c
@@ -384,8 +384,8 @@ GRT_MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, rea
         nlay++;
     }
 
-    mod1d->isrc = isrc;
-    mod1d->ircv = ircv;
+    mod1d->isrc = (depsrc >= 0.0)? isrc : nlay; // 使用明显越界的索引，表明没有插入震源层或台站层
+    mod1d->ircv = (deprcv >= 0.0)? ircv : nlay;
     mod1d->ircvup = ircvup;
     mod1d->n = nlay;
     mod1d->depsrc = depsrc;
diff --git a/pygrt/C_extension/src/dynamic/layer.c b/pygrt/C_extension/src/dynamic/layer.c
index 9e2bfde..919a31c 100644
--- a/pygrt/C_extension/src/dynamic/layer.c
+++ b/pygrt/C_extension/src/dynamic/layer.c
@@ -29,7 +29,7 @@
     T N##2 = mod1d->N[iy];\
 
 
-static int __freeBound_R(bool isLiquid, cplx_t cbcb0, cplx_t xa0, cplx_t xb0, real_t adsgn, cplx_t M[2][2], cplx_t *ML)
+static int __freeBound_R_PSV(bool isLiquid, cplx_t cbcb0, cplx_t xa0, cplx_t xb0, real_t adsgn, cplx_t M[2][2])
 {
     if( ! isLiquid){
         // 固体表面
@@ -47,18 +47,29 @@ static int __freeBound_R(bool isLiquid, cplx_t cbcb0, cplx_t xa0, cplx_t xb0, re
         M[0][1] = 2.0 * xb0 * (1.0 - 0.5*cbcb0) * Delta * adsgn;
         M[1][0] = 2.0 * xa0 * (1.0 - 0.5*cbcb0) * Delta * adsgn;
         M[1][1] = M[0][0];
-        *ML = 1.0;
     }
     else {
         // 液体表面
         M[0][0] = -1.0;
         M[1][1] = M[0][1] = M[1][0] = 0.0;
+    }
+    return GRT_INVERSE_SUCCESS;
+}
+
+static int __freeBound_R_SH(bool isLiquid, cplx_t *ML)
+{
+    if( ! isLiquid){
+        // 固体表面
+        *ML = 1.0;
+    }
+    else {
+        // 液体表面
         *ML = 0.0;
     }
     return GRT_INVERSE_SUCCESS;
 }
 
-static int __rigidBound_R(bool isLiquid, cplx_t xa0, cplx_t xb0, real_t adsgn, cplx_t M[2][2], cplx_t *ML)
+static int __rigidBound_R_PSV(bool isLiquid, cplx_t xa0, cplx_t xb0, real_t adsgn, cplx_t M[2][2])
 {
     if( ! isLiquid){
         // 固体表面
@@ -73,31 +84,60 @@ static int __rigidBound_R(bool isLiquid, cplx_t xa0, cplx_t xb0, real_t adsgn, c
         M[0][1] = 2.0 * xb0 * Delta * adsgn;
         M[1][0] = 2.0 * xa0 * Delta * adsgn;
         M[1][1] = M[0][0];
-        *ML = - 1.0;
     }
     else {
         // 液体表面
         M[0][0] = 1.0;
         M[1][1] = M[0][1] = M[1][0] = 0.0;
+    }
+    return GRT_INVERSE_SUCCESS;
+}
+
+static int __rigidBound_R_SH(bool isLiquid, cplx_t *ML)
+{
+    if( ! isLiquid){
+        // 固体表面
+        *ML = - 1.0;
+    }
+    else {
+        // 液体表面
         *ML = 0.0;
     }
     return GRT_INVERSE_SUCCESS;
 }
 
-void grt_topbound_RU(GRT_MODEL1D *mod1d)
+void grt_topbound_RU_PSV(GRT_MODEL1D *mod1d)
 {
     cplx_t cbcb = mod1d->cbcb[0];
     cplx_t xa = mod1d->xa[0];
     cplx_t xb = mod1d->xb[0];
     bool   isLiquid = mod1d->isLiquid[0];
     if(mod1d->topbound == GRT_BOUND_FREE){
-        mod1d->M_top.stats = __freeBound_R(isLiquid, cbcb, xa, xb, 1.0, mod1d->M_top.RU, &mod1d->M_top.RUL);
+        mod1d->M_top.stats = __freeBound_R_PSV(isLiquid, cbcb, xa, xb, 1.0, mod1d->M_top.RU);
     }
     else if(mod1d->topbound == GRT_BOUND_RIGID){
-        mod1d->M_top.stats = __rigidBound_R(isLiquid, xa, xb, 1.0, mod1d->M_top.RU, &mod1d->M_top.RUL);
+        mod1d->M_top.stats = __rigidBound_R_PSV(isLiquid, xa, xb, 1.0, mod1d->M_top.RU);
     }
     else if(mod1d->topbound == GRT_BOUND_HALFSPACE){
         memset(mod1d->M_top.RU, 0, sizeof(cplx_t)*4);
+        mod1d->M_top.stats = GRT_INVERSE_SUCCESS;
+    }
+    else{
+        GRTRaiseError("Wrong execution.");
+    }
+    // RU 不需要时延
+}
+
+void grt_topbound_RU_SH(GRT_MODEL1D *mod1d)
+{
+    bool isLiquid = mod1d->isLiquid[0];
+    if(mod1d->topbound == GRT_BOUND_FREE){
+        mod1d->M_top.stats = __freeBound_R_SH(isLiquid, &mod1d->M_top.RUL);
+    }
+    else if(mod1d->topbound == GRT_BOUND_RIGID){
+        mod1d->M_top.stats = __rigidBound_R_SH(isLiquid, &mod1d->M_top.RUL);
+    }
+    else if(mod1d->topbound == GRT_BOUND_HALFSPACE){
         mod1d->M_top.RUL = 0.0;
         mod1d->M_top.stats = GRT_INVERSE_SUCCESS;
     }
@@ -107,22 +147,21 @@ void grt_topbound_RU(GRT_MODEL1D *mod1d)
     // RU 不需要时延
 }
 
-void grt_botbound_RD(GRT_MODEL1D *mod1d)
+void grt_botbound_RD_PSV(GRT_MODEL1D *mod1d)
 {
     size_t nlay = mod1d->n;
     cplx_t cbcb = mod1d->cbcb[nlay-2];
     cplx_t xa = mod1d->xa[nlay-2];
     cplx_t xb = mod1d->xb[nlay-2];
-    cplx_t isLiquid = mod1d->isLiquid[nlay-2];
+    bool isLiquid = mod1d->isLiquid[nlay-2];
     if(mod1d->botbound == GRT_BOUND_FREE){
-        mod1d->M_bot.stats = __freeBound_R(isLiquid, cbcb, xa, xb, -1.0, mod1d->M_bot.RD, &mod1d->M_bot.RDL);
+        mod1d->M_bot.stats = __freeBound_R_PSV(isLiquid, cbcb, xa, xb, -1.0, mod1d->M_bot.RD);
     }
     else if(mod1d->botbound == GRT_BOUND_RIGID){
-        mod1d->M_bot.stats = __rigidBound_R(cbcb, xa, xb, -1.0, mod1d->M_bot.RD, &mod1d->M_bot.RDL);
+        mod1d->M_bot.stats = __rigidBound_R_PSV(isLiquid, xa, xb, -1.0, mod1d->M_bot.RD);
     }
     else if(mod1d->botbound == GRT_BOUND_HALFSPACE){
         grt_RT_matrix_PSV(mod1d, nlay-1, &mod1d->M_bot);
-        grt_RT_matrix_SH(mod1d, nlay-1, &mod1d->M_bot);
     }
     else{
         GRTRaiseError("Wrong execution.");
@@ -141,6 +180,33 @@ void grt_botbound_RD(GRT_MODEL1D *mod1d)
 
     mod1d->M_bot.RD[0][0] *= ex2a;   mod1d->M_bot.RD[0][1] *= exab;
     mod1d->M_bot.RD[1][0] *= exab;   mod1d->M_bot.RD[1][1] *= ex2b;
+}
+
+void grt_botbound_RD_SH(GRT_MODEL1D *mod1d)
+{
+    size_t nlay = mod1d->n;
+    cplx_t xb = mod1d->xb[nlay-2];
+    bool isLiquid = mod1d->isLiquid[nlay-2];
+    if(mod1d->botbound == GRT_BOUND_FREE){
+        mod1d->M_bot.stats = __freeBound_R_SH(isLiquid, &mod1d->M_bot.RDL);
+    }
+    else if(mod1d->botbound == GRT_BOUND_RIGID){
+        mod1d->M_bot.stats = __rigidBound_R_SH(isLiquid, &mod1d->M_bot.RDL);
+    }
+    else if(mod1d->botbound == GRT_BOUND_HALFSPACE){
+        grt_RT_matrix_SH(mod1d, nlay-1, &mod1d->M_bot);
+    }
+    else{
+        GRTRaiseError("Wrong execution.");
+    }
+
+    // 时延 RD
+    real_t k = mod1d->k;
+    real_t thk = mod1d->Thk[nlay-2];
+
+    cplx_t exb, ex2b;
+    exb = exp(- k*thk*xb);
+    ex2b = exb * exb;
 
     mod1d->M_bot.RDL *= ex2b;
 }
@@ -558,6 +624,45 @@ void grt_delay_RT_matrix(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M
     M->TUL *= exb;
 }
 
+void grt_delay_RT_matrix_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M)
+{
+    real_t thk = mod1d->Thk[iy-1];
+    cplx_t xa1 = mod1d->xa[iy-1];
+    cplx_t xb1 = mod1d->xb[iy-1];
+    real_t k = mod1d->k;
+
+    cplx_t exa, exb, ex2a, ex2b, exab;
+    exa = exp(- k*thk*xa1);
+    exb = exp(- k*thk*xb1);
+    ex2a = exa * exa;
+    ex2b = exb * exb;
+    exab = exa * exb;
+
+    M->RD[0][0] *= ex2a;   M->RD[0][1] *= exab;
+    M->RD[1][0] *= exab;   M->RD[1][1] *= ex2b;
+
+    M->TD[0][0] *= exa;    M->TD[0][1] *= exb;
+    M->TD[1][0] *= exa;    M->TD[1][1] *= exb;
+
+    M->TU[0][0] *= exa;    M->TU[0][1] *= exa;
+    M->TU[1][0] *= exb;    M->TU[1][1] *= exb;
+}
+
+void grt_delay_RT_matrix_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M)
+{
+    real_t thk = mod1d->Thk[iy-1];
+    cplx_t xb1 = mod1d->xb[iy-1];
+    real_t k = mod1d->k;
+
+    cplx_t exb, ex2b;
+    exb = exp(- k*thk*xb1);
+    ex2b = exb * exb;
+
+    M->RDL *= ex2b;
+    M->TDL *= exb;
+    M->TUL *= exb;
+}
+
 
 
 void grt_delay_GRT_matrix(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M)
diff --git a/pygrt/C_extension/src/integral/kernel_template.c_ b/pygrt/C_extension/src/integral/kernel_template.c_
index 1448a52..06b64a9 100644
--- a/pygrt/C_extension/src/integral/kernel_template.c_
+++ b/pygrt/C_extension/src/integral/kernel_template.c_
@@ -47,8 +47,10 @@
 #define __grt_RT_matrix_SH        CONCAT_PREFIX(RT_matrix_SH)
 #define __grt_delay_RT_matrix     CONCAT_PREFIX(delay_RT_matrix)
 #define __grt_source_coef         CONCAT_PREFIX(source_coef)
-#define __grt_topbound_RU         CONCAT_PREFIX(topbound_RU)
-#define __grt_botbound_RD         CONCAT_PREFIX(botbound_RD)
+#define __grt_topbound_RU_PSV     CONCAT_PREFIX(topbound_RU_PSV)
+#define __grt_topbound_RU_SH      CONCAT_PREFIX(topbound_RU_SH)
+#define __grt_botbound_RD_PSV     CONCAT_PREFIX(botbound_RD_PSV)
+#define __grt_botbound_RD_SH      CONCAT_PREFIX(botbound_RD_SH)
 #define __grt_wave2qwv_REV_PSV    CONCAT_PREFIX(wave2qwv_REV_PSV)
 #define __grt_wave2qwv_REV_SH     CONCAT_PREFIX(wave2qwv_REV_SH)
 #define __grt_wave2qwv_z_REV_PSV  CONCAT_PREFIX(wave2qwv_z_REV_PSV)
@@ -153,14 +155,16 @@ void __GRT_MATRIX_FUNC__(GRT_MODEL1D *mod1d, const real_t k)
     //===================================================================================
 
     // 递推RU_FA
-    __grt_topbound_RU(mod1d);
+    __grt_topbound_RU_PSV(mod1d);
+    __grt_topbound_RU_SH(mod1d);
     if(M_top->stats==GRT_INVERSE_FAILURE)  goto BEFORE_RETURN;
     grt_recursion_RU(M_top, M_FA, M_FA);
     if(M_FA->stats==GRT_INVERSE_FAILURE)  goto BEFORE_RETURN;
 
     // 递推RD_BL, 在模型定义阶段已检查震源/台站的虚拟层成为底层时，只能是半无限空间
     if(imax+1 < mod1d->n){ // 此 if 仅减少不必要的计算
-        __grt_botbound_RD(mod1d);
+        __grt_botbound_RD_PSV(mod1d);
+        __grt_botbound_RD_SH(mod1d);
         if(M_bot->stats==GRT_INVERSE_FAILURE)  goto BEFORE_RETURN;
         grt_recursion_RD(M_BL, M_bot, M_BL);
         if(M_BL->stats==GRT_INVERSE_FAILURE)  goto BEFORE_RETURN;
diff --git a/pygrt/C_extension/src/static/static_layer.c b/pygrt/C_extension/src/static/static_layer.c
index a614a62..8392499 100644
--- a/pygrt/C_extension/src/static/static_layer.c
+++ b/pygrt/C_extension/src/static/static_layer.c
@@ -26,7 +26,7 @@
     T N##2 = mod1d->N[iy];\
 
 
-static int __freeBound_R(real_t d, real_t k, cplx_t delta1, cplx_t M[2][2], cplx_t *ML)
+static int __freeBound_R_PSV(real_t d, real_t k, cplx_t delta1, cplx_t M[2][2])
 {
     if(d <= 0.0){
         // 公式(6.3.12)
@@ -38,11 +38,16 @@ static int __freeBound_R(real_t d, real_t k, cplx_t delta1, cplx_t M[2][2], cplx
         M[0][1] = delta1 * (-4.0*d*d*k*k - 1.0);
         M[1][0] = -1.0/delta1;
     }
+    return GRT_INVERSE_SUCCESS;
+}
+
+static int __freeBound_R_SH(cplx_t *ML)
+{
     *ML = 1.0;
     return GRT_INVERSE_SUCCESS;
 }
 
-static int __rigidBound_R(real_t d, real_t k, cplx_t delta1, cplx_t M[2][2], cplx_t *ML)
+static int __rigidBound_R_PSV(real_t d, real_t k, cplx_t delta1, cplx_t M[2][2])
 {
     if(d <= 0.0){
         // 公式(6.3.12)
@@ -53,22 +58,44 @@ static int __rigidBound_R(real_t d, real_t k, cplx_t delta1, cplx_t M[2][2], cpl
         M[0][1] = (4.0*GRT_SQUARE(delta1*d*k) + 1.0);
         M[1][0] = 1.0;
     }
+    return GRT_INVERSE_SUCCESS;
+}
+
+static int __rigidBound_R_SH(cplx_t *ML)
+{
     *ML = -1.0;
     return GRT_INVERSE_SUCCESS;
 }
 
-void grt_static_topbound_RU(GRT_MODEL1D *mod1d)
+void grt_static_topbound_RU_PSV(GRT_MODEL1D *mod1d)
 {
     cplx_t delta = mod1d->delta[0];
     real_t k = mod1d->k;
     if(mod1d->topbound == GRT_BOUND_FREE){
-        mod1d->M_top.stats = __freeBound_R(0.0, k, delta, mod1d->M_top.RU, &mod1d->M_top.RUL);
+        mod1d->M_top.stats = __freeBound_R_PSV(0.0, k, delta, mod1d->M_top.RU);
     }
     else if(mod1d->topbound == GRT_BOUND_RIGID){
-        mod1d->M_top.stats = __rigidBound_R(0.0, k, delta, mod1d->M_top.RU, &mod1d->M_top.RUL);
+        mod1d->M_top.stats = __rigidBound_R_PSV(0.0, k, delta, mod1d->M_top.RU);
     }
     else if(mod1d->topbound == GRT_BOUND_HALFSPACE){
         memset(mod1d->M_top.RU, 0, sizeof(cplx_t)*4);
+        mod1d->M_top.stats = GRT_INVERSE_SUCCESS;
+    }
+    else{
+        GRTRaiseError("Wrong execution.");
+    }
+    // RU 不需要时延
+}
+
+void grt_static_topbound_RU_SH(GRT_MODEL1D *mod1d)
+{
+    if(mod1d->topbound == GRT_BOUND_FREE){
+        mod1d->M_top.stats = __freeBound_R_SH(&mod1d->M_top.RUL);
+    }
+    else if(mod1d->topbound == GRT_BOUND_RIGID){
+        mod1d->M_top.stats = __rigidBound_R_SH(&mod1d->M_top.RUL);
+    }
+    else if(mod1d->topbound == GRT_BOUND_HALFSPACE){
         mod1d->M_top.RUL = 0.0;
         mod1d->M_top.stats = GRT_INVERSE_SUCCESS;
     }
@@ -78,21 +105,20 @@ void grt_static_topbound_RU(GRT_MODEL1D *mod1d)
     // RU 不需要时延
 }
 
-void grt_static_botbound_RD(GRT_MODEL1D *mod1d)
+void grt_static_botbound_RD_PSV(GRT_MODEL1D *mod1d)
 {
     size_t nlay = mod1d->n;
     cplx_t delta = mod1d->delta[nlay-2];
     real_t thk = mod1d->Thk[nlay-2];
     real_t k = mod1d->k;
     if(mod1d->botbound == GRT_BOUND_FREE){
-        mod1d->M_bot.stats = __freeBound_R(thk, k, delta, mod1d->M_bot.RD, &mod1d->M_bot.RDL);
+        mod1d->M_bot.stats = __freeBound_R_PSV(thk, k, delta, mod1d->M_bot.RD);
     }
     else if(mod1d->botbound == GRT_BOUND_RIGID){
-        mod1d->M_bot.stats = __rigidBound_R(thk, k, delta, mod1d->M_bot.RD, &mod1d->M_bot.RDL);
+        mod1d->M_bot.stats = __rigidBound_R_PSV(thk, k, delta, mod1d->M_bot.RD);
     }
     else if(mod1d->botbound == GRT_BOUND_HALFSPACE){
         grt_static_RT_matrix_PSV(mod1d, nlay-1, &mod1d->M_bot);
-        grt_static_RT_matrix_SH(mod1d, nlay-1, &mod1d->M_bot);
     }
     else{
         GRTRaiseError("Wrong execution.");
@@ -105,6 +131,30 @@ void grt_static_botbound_RD(GRT_MODEL1D *mod1d)
 
     mod1d->M_bot.RD[0][0] *= ex2;   mod1d->M_bot.RD[0][1] *= ex2;
     mod1d->M_bot.RD[1][0] *= ex2;   mod1d->M_bot.RD[1][1] *= ex2;
+}
+
+void grt_static_botbound_RD_SH(GRT_MODEL1D *mod1d)
+{
+    size_t nlay = mod1d->n;
+    real_t thk = mod1d->Thk[nlay-2];
+    real_t k = mod1d->k;
+    if(mod1d->botbound == GRT_BOUND_FREE){
+        mod1d->M_bot.stats = __freeBound_R_SH(&mod1d->M_bot.RDL);
+    }
+    else if(mod1d->botbound == GRT_BOUND_RIGID){
+        mod1d->M_bot.stats = __rigidBound_R_SH(&mod1d->M_bot.RDL);
+    }
+    else if(mod1d->botbound == GRT_BOUND_HALFSPACE){
+        grt_static_RT_matrix_SH(mod1d, nlay-1, &mod1d->M_bot);
+    }
+    else{
+        GRTRaiseError("Wrong execution.");
+    }
+
+    // 时延 RD
+    cplx_t ex, ex2;
+    ex = exp(- k*thk);
+    ex2 = ex * ex;
 
     mod1d->M_bot.RDL *= ex2;
 }