From 28af3f40088ea7f2eb33121f9595e3b72659727c Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Wed, 20 Aug 2025 17:58:25 +0800 Subject: [PATCH] REFAC: Split source coefs into P-SV and SH --- pygrt/C_extension/include/common/recursion.h | 6 ++- pygrt/C_extension/include/dynamic/source.h | 12 ++++-- .../include/static/static_source.h | 9 ++-- pygrt/C_extension/src/common/recursion.c | 11 ++--- pygrt/C_extension/src/dynamic/propagate.c | 14 ++++--- pygrt/C_extension/src/dynamic/source.c | 42 +++++++++++++++---- .../C_extension/src/static/static_propagate.c | 14 ++++--- pygrt/C_extension/src/static/static_source.c | 31 ++++++++++---- 8 files changed, 97 insertions(+), 42 deletions(-) diff --git a/pygrt/C_extension/include/common/recursion.h b/pygrt/C_extension/include/common/recursion.h index b4e7ac5..8cf62da 100644 --- a/pygrt/C_extension/include/common/recursion.h +++ b/pygrt/C_extension/include/common/recursion.h @@ -354,11 +354,13 @@ void recursion_RT_SH_imaginary( * @param[in] RL1 SH波， \f$ R_1\f$ * @param[in] R2 P-SV波，\f$\mathbf{R_2}\f$矩阵 * @param[in] RL2 SH波， \f$ R_2\f$ - * @param[in] coef 震源系数，\f$ P_m, SV_m, SH_m \f$ ，维度2表示下行波(p=0)和上行波(p=1) + * @param[in] coef_PSV P-SV 波震源系数，\f$ P_m, SV_m\f$ ，维度2表示下行波(p=0)和上行波(p=1) + * @param[in] coef_SH SH 波震源系数，\f$ SH_m \f$ ，维度2表示下行波(p=0)和上行波(p=1) * @param[out] qwv 最终通过矩阵传播计算出的在台站位置的\f$ q_m,w_m,v_m\f$ */ void get_qwv( bool ircvup, const MYCOMPLEX R1[2][2], MYCOMPLEX RL1, const MYCOMPLEX R2[2][2], MYCOMPLEX RL2, - const MYCOMPLEX coef[QWV_NUM][2], MYCOMPLEX qwv[QWV_NUM]); + const MYCOMPLEX coef_PSV[QWV_NUM-1][2], const MYCOMPLEX coef_SH[2], + MYCOMPLEX qwv[QWV_NUM]); diff --git a/pygrt/C_extension/include/dynamic/source.h b/pygrt/C_extension/include/dynamic/source.h index e5d244b..fc9fbae 100755 --- a/pygrt/C_extension/include/dynamic/source.h +++ b/pygrt/C_extension/include/dynamic/source.h @@ -18,7 +18,7 @@ /** * 根据公式(4.6.6)，(4.6.15)，(4.6.21,26)，(4.8.34)计算不同震源不同阶数的震源系数， * 数组形状代表在[i][j][p]时表示i类震源的 - * P(j=0),SV(j=1),SH(j=2)的震源系数(分别对应q,w,v)，且分为下行波(p=0)和上行波(p=1). + * P(j=0),SV(j=1)的震源系数(分别对应q,w)，且分为下行波(p=0)和上行波(p=1). * * @param[in] src_xa 震源层的P波归一化垂直波数 \f$ \sqrt{1 - (k_a/k)^2} \f$ * @param[in] src_xb 震源层的S波归一化垂直波数 \f$ \sqrt{1 - (k_b/k)^2} \f$ @@ -28,7 +28,13 @@ * @param[out] coef 震源系数 \f$ P_m, SV_m, SH_m \f$ * */ -void source_coef( +void source_coef_PSV( MYCOMPLEX src_xa, MYCOMPLEX src_xb, MYCOMPLEX src_kaka, MYCOMPLEX src_kbkb, MYREAL k, - MYCOMPLEX coef[SRC_M_NUM][QWV_NUM][2]); + MYCOMPLEX coef[SRC_M_NUM][QWV_NUM-1][2]); + +/* SH 波的震源系数，参数见 source_coef_PSV */ +void source_coef_SH( + MYCOMPLEX src_xb, MYCOMPLEX src_kbkb, + MYREAL k, + MYCOMPLEX coef[SRC_M_NUM][2]); diff --git a/pygrt/C_extension/include/static/static_source.h b/pygrt/C_extension/include/static/static_source.h index db82412..f24cdcf 100644 --- a/pygrt/C_extension/include/static/static_source.h +++ b/pygrt/C_extension/include/static/static_source.h @@ -16,12 +16,13 @@ * 计算不同震源的静态震源系数，文献/书中仅提供剪切源的震源系数，其它震源系数重新推导 * * 数组形状代表在[i][j][p]时表示i类震源的 - * P(j=0),SV(j=1),SH(j=2)的震源系数(分别对应q,w,v)，且分为下行波(p=0)和上行波(p=1). + * P(j=0),SV(j=1)的震源系数(分别对应q,w)，且分为下行波(p=0)和上行波(p=1). * * @param[in] delta 震源层的\f$ \Delta \f$ * @param[in] k 波数 * @param[out] coef 震源系数 \f$ P_m, SV_m, SH_m \f$ */ -void static_source_coef( - MYCOMPLEX delta, MYREAL k, MYCOMPLEX coef[SRC_M_NUM][QWV_NUM][2]); - \ No newline at end of file +void static_source_coef_PSV(MYCOMPLEX delta, MYREAL k, MYCOMPLEX coef[SRC_M_NUM][QWV_NUM-1][2]); + +/* SH 波的静态震源系数，参数见 static_source_coef_PSV */ +void static_source_coef_SH(MYREAL k, MYCOMPLEX coef[SRC_M_NUM][2]); \ No newline at end of file diff --git a/pygrt/C_extension/src/common/recursion.c b/pygrt/C_extension/src/common/recursion.c index a8b85f9..676ca5b 100644 --- a/pygrt/C_extension/src/common/recursion.c +++ b/pygrt/C_extension/src/common/recursion.c @@ -380,19 +380,20 @@ void get_qwv( bool ircvup, const MYCOMPLEX R1[2][2], MYCOMPLEX RL1, const MYCOMPLEX R2[2][2], MYCOMPLEX RL2, - const MYCOMPLEX coef[QWV_NUM][2], MYCOMPLEX qwv[QWV_NUM]) + const MYCOMPLEX coef_PSV[QWV_NUM-1][2], const MYCOMPLEX coef_SH[2], + MYCOMPLEX qwv[QWV_NUM]) { MYCOMPLEX qw0[2], qw1[2], v0; - MYCOMPLEX coefD[2] = {coef[0][0], coef[1][0]}; - MYCOMPLEX coefU[2] = {coef[0][1], coef[1][1]}; + MYCOMPLEX coefD[2] = {coef_PSV[0][0], coef_PSV[1][0]}; + MYCOMPLEX coefU[2] = {coef_PSV[0][1], coef_PSV[1][1]}; if(ircvup){ cmat2x1_mul(R2, coefD, qw0); qw0[0] += coefU[0]; qw0[1] += coefU[1]; - v0 = RL1 * (RL2*coef[2][0] + coef[2][1]); + v0 = RL1 * (RL2*coef_SH[0] + coef_SH[1]); } else { cmat2x1_mul(R2, coefU, qw0); qw0[0] += coefD[0]; qw0[1] += coefD[1]; - v0 = RL1 * (coef[2][0] + RL2*coef[2][1]); + v0 = RL1 * (coef_SH[0] + RL2*coef_SH[1]); } cmat2x1_mul(R1, qw0, qw1); diff --git a/pygrt/C_extension/src/dynamic/propagate.c b/pygrt/C_extension/src/dynamic/propagate.c index 7844af3..12b065d 100755 --- a/pygrt/C_extension/src/dynamic/propagate.c +++ b/pygrt/C_extension/src/dynamic/propagate.c @@ -273,8 +273,10 @@ void kernel( // 计算震源系数 - MYCOMPLEX src_coef[SRC_M_NUM][QWV_NUM][2] = {0}; - source_coef(src_xa, src_xb, src_kaka, src_kbkb, k, src_coef); + MYCOMPLEX src_coef_PSV[SRC_M_NUM][QWV_NUM-1][2] = {0}; + MYCOMPLEX src_coef_SH[SRC_M_NUM][2] = {0}; + source_coef_PSV(src_xa, src_xb, src_kaka, src_kbkb, k, src_coef_PSV); + source_coef_SH(src_xb, src_kbkb, k, src_coef_SH); // 临时中转矩阵 (temperary) MYCOMPLEX tmpR2[2][2], tmp2x2[2][2], tmpRL, tmp2x2_uiz[2][2], tmpRL2; @@ -357,7 +359,7 @@ void kernel( tmpRL2 = R_EVL * tmpRL; for(MYINT i=0; i