diff --git a/pygrt/C_extension/include/grt/common/kernel.h b/pygrt/C_extension/include/grt/common/kernel.h index ff2ccef9..e62b2304 100644 --- a/pygrt/C_extension/include/grt/common/kernel.h +++ b/pygrt/C_extension/include/grt/common/kernel.h @@ -44,8 +44,7 @@ typedef void (*GRT_KernelFunc) ( #define __grt_RT_matrix_PSV CONCAT_PREFIX(RT_matrix_PSV) #define __grt_RT_matrix_SH CONCAT_PREFIX(RT_matrix_SH) #define __grt_delay_RT_matrix CONCAT_PREFIX(delay_RT_matrix) - #define __grt_source_coef_PSV CONCAT_PREFIX(source_coef_PSV) - #define __grt_source_coef_SH CONCAT_PREFIX(source_coef_SH) + #define __grt_source_coef CONCAT_PREFIX(source_coef) #define __grt_topfree_RU CONCAT_PREFIX(topfree_RU) #define __grt_wave2qwv_REV_PSV CONCAT_PREFIX(wave2qwv_REV_PSV) #define __grt_wave2qwv_REV_SH CONCAT_PREFIX(wave2qwv_REV_SH) diff --git a/pygrt/C_extension/include/grt/common/kernel_template.h b/pygrt/C_extension/include/grt/common/kernel_template.h index 6651f08b..9c44e777 100644 --- a/pygrt/C_extension/include/grt/common/kernel_template.h +++ b/pygrt/C_extension/include/grt/common/kernel_template.h @@ -58,28 +58,26 @@ * @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_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[in] coef 震源系数，\f$ P_m, SV_m，SH_m\f$ ，维度2表示下行波(p=0)和上行波(p=1) * @param[out] qwv 最终通过矩阵传播计算出的在台站位置的\f$ q_m,w_m,v_m\f$ */ inline void grt_construct_qwv( bool ircvup, const cplx_t R1[2][2], cplx_t RL1, const cplx_t R2[2][2], cplx_t RL2, - const cplx_t coef_PSV[GRT_QWV_NUM-1][2], const cplx_t coef_SH[2], - cplx_t qwv[GRT_QWV_NUM]) + const cplx_t coef[GRT_QWV_NUM][2], cplx_t qwv[GRT_QWV_NUM]) { cplx_t qw0[2], qw1[2], v0; - cplx_t coefD[2] = {coef_PSV[0][0], coef_PSV[1][0]}; - cplx_t coefU[2] = {coef_PSV[0][1], coef_PSV[1][1]}; + cplx_t coefD[2] = {coef[0][0], coef[1][0]}; + cplx_t coefU[2] = {coef[0][1], coef[1][1]}; if(ircvup){ grt_cmat2x1_mul(R2, coefD, qw0); qw0[0] += coefU[0]; qw0[1] += coefU[1]; - v0 = RL1 * (RL2*coef_SH[0] + coef_SH[1]); + v0 = RL1 * (RL2*coef[2][0] + coef[2][1]); } else { grt_cmat2x1_mul(R2, coefU, qw0); qw0[0] += coefD[0]; qw0[1] += coefD[1]; - v0 = RL1 * (coef_SH[0] + RL2*coef_SH[1]); + v0 = RL1 * (coef[2][0] + RL2*coef[2][1]); } grt_cmat2x1_mul(R1, qw0, qw1); @@ -256,10 +254,8 @@ static void __KERNEL_FUNC__( // 计算震源系数 - cplx_t src_coef_PSV[GRT_SRC_M_NUM][GRT_QWV_NUM-1][2] = {0}; - cplx_t src_coef_SH[GRT_SRC_M_NUM][2] = {0}; - __grt_source_coef_PSV(mod1d, src_coef_PSV); - __grt_source_coef_SH(mod1d, src_coef_SH); + cplx_t src_coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2] = {0}; + __grt_source_coef(mod1d, src_coef); // 临时中转矩阵 (temperary) cplx_t tmpR2[2][2], tmp2x2[2][2], tmpRL, tmp2x2_uiz[2][2], tmpRL2; @@ -295,7 +291,7 @@ static void __KERNEL_FUNC__( tmpRL2 = R_EVL * tmpRL; for(int i=0; iRD, M_BL->RDL, src_coef_PSV[i], src_coef_SH[i], QWV[i]); + grt_construct_qwv(ircvup, tmp2x2, tmpRL2, M_BL->RD, M_BL->RDL, src_coef[i], QWV[i]); } @@ -306,7 +302,7 @@ static void __KERNEL_FUNC__( tmpRL2 = uiz_R_EVL * tmpRL; for(int i=0; iRD, M_BL->RDL, src_coef_PSV[i], src_coef_SH[i], QWV_uiz[i]); + grt_construct_qwv(ircvup, tmp2x2_uiz, tmpRL2, M_BL->RD, M_BL->RDL, src_coef[i], QWV_uiz[i]); } } } @@ -334,7 +330,7 @@ static void __KERNEL_FUNC__( tmpRL2 = R_EVL * tmpRL; for(int i=0; iRU, M_FA->RUL, src_coef_PSV[i], src_coef_SH[i], QWV[i]); + grt_construct_qwv(ircvup, tmp2x2, tmpRL2, M_FA->RU, M_FA->RUL, src_coef[i], QWV[i]); } @@ -345,7 +341,7 @@ static void __KERNEL_FUNC__( tmpRL2 = uiz_R_EVL * tmpRL; for(int i=0; iRU, M_FA->RUL, src_coef_PSV[i], src_coef_SH[i], QWV_uiz[i]); + grt_construct_qwv(ircvup, tmp2x2_uiz, tmpRL2, M_FA->RU, M_FA->RUL, src_coef[i], QWV_uiz[i]); } } diff --git a/pygrt/C_extension/include/grt/dynamic/source.h b/pygrt/C_extension/include/grt/dynamic/source.h index e04c8377..f0ff9520 100755 --- a/pygrt/C_extension/include/grt/dynamic/source.h +++ b/pygrt/C_extension/include/grt/dynamic/source.h @@ -24,7 +24,10 @@ * @param[out] coef 震源系数 \f$ P_m, SV_m, SH_m \f$ * */ -void grt_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM-1][2]); +void grt_source_coef(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]); + +/* P-SV 波的震源系数 */ +void grt_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]); /* SH 波的震源系数 */ -void grt_source_coef_SH(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][2]); +void grt_source_coef_SH(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]); diff --git a/pygrt/C_extension/include/grt/static/static_source.h b/pygrt/C_extension/include/grt/static/static_source.h index cc11e54c..ef591cd6 100644 --- a/pygrt/C_extension/include/grt/static/static_source.h +++ b/pygrt/C_extension/include/grt/static/static_source.h @@ -22,7 +22,11 @@ * @param[in] mod1d 模型结构体指针 * @param[out] coef 震源系数 \f$ P_m, SV_m, SH_m \f$ */ -void grt_static_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM-1][2]); +void grt_static_source_coef(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]); + + +/* P-SV 波的静态震源系数 */ +void grt_static_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]); /* SH 波的静态震源系数 */ -void grt_static_source_coef_SH(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][2]); \ No newline at end of file +void grt_static_source_coef_SH(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]); \ No newline at end of file diff --git a/pygrt/C_extension/src/dynamic/source.c b/pygrt/C_extension/src/dynamic/source.c index 594724e9..4d96f068 100755 --- a/pygrt/C_extension/src/dynamic/source.c +++ b/pygrt/C_extension/src/dynamic/source.c @@ -11,93 +11,89 @@ #include #include +#include #include "grt/dynamic/source.h" -#include "grt/common/model.h" -#include "grt/common/matrix.h" -#include "grt/common/prtdbg.h" -void grt_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM-1][2]) +inline void _source_PSV( + const cplx_t xa, const cplx_t caca, + const cplx_t xb, const cplx_t cbcb, const real_t k, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]) { - // 先全部赋0 - for(int i=0; iisrc; - cplx_t xa = mod1d->xa[isrc]; - cplx_t caca = mod1d->caca[isrc]; - cplx_t xb = mod1d->xb[isrc]; - cplx_t cbcb = mod1d->cbcb[isrc]; - real_t k = mod1d->k; - cplx_t tmp; - // 爆炸源， 通过(4.9.8)的矩张量源公式，提取各向同性的量(M11+M22+M33)，-a+k^2/a -> ka^2/a - coef[0][0][0] = tmp = (caca / xa) * k; coef[0][0][1] = tmp; + coef[0][0][0] = tmp = (caca / xa) * k; coef[0][0][1] = tmp; // 垂直力源 (4.6.15) - coef[1][0][0] = tmp = -1.0; coef[1][0][1] = - tmp; + coef[1][0][0] = tmp = -1.0; coef[1][0][1] = - tmp; coef[1][1][0] = tmp = -1.0 / xb; coef[1][1][1] = tmp; - // 水平力源 (4.6.21,26), 这里可以把x1,x2方向的力转到r,theta方向 - // 推导可发现，r方向的力形成P,SV波, theta方向的力形成SH波 - // 方向性因子包含水平力方向与震源台站连线方向的夹角 + // 水平力源 (4.6.21,26) coef[2][0][0] = tmp = -1.0 / xa; coef[2][0][1] = tmp; - coef[2][1][0] = tmp = -1.0; coef[2][1][1] = - tmp; + coef[2][1][0] = tmp = -1.0; coef[2][1][1] = - tmp; // 剪切位错 (4.8.34) // m=0 coef[3][0][0] = tmp = ((2.0*caca - 3.0) / xa) * k; coef[3][0][1] = tmp; - coef[3][1][0] = tmp = -3.0*k; coef[3][1][1] = - tmp; + coef[3][1][0] = tmp = -3.0*k; coef[3][1][1] = - tmp; // m=1 coef[4][0][0] = tmp = 2.0*k; coef[4][0][1] = - tmp; coef[4][1][0] = tmp = ((2.0 - cbcb) / xb) * k; coef[4][1][1] = tmp; // m=2 - coef[5][0][0] = tmp = - (1.0 / xa) * k; coef[5][0][1] = tmp; + coef[5][0][0] = tmp = - (1.0 / xa) * k; coef[5][0][1] = tmp; coef[5][1][0] = tmp = - k; coef[5][1][1] = - tmp; } +inline void _source_SH(const cplx_t xb, const cplx_t cbcb, const real_t k, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]) +{ + cplx_t tmp; + + // 水平力源 (4.6.21,26) + coef[2][2][0] = tmp = cbcb / xb; coef[2][2][1] = tmp; + + // 剪切位错 (4.8.34) + // m=1 + coef[4][2][0] = tmp = - cbcb * k; coef[4][2][1] = - tmp; + + // m=2 + coef[5][2][0] = tmp = (cbcb / xb) * k; coef[5][2][1] = tmp; +} + -void grt_source_coef_SH(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][2]) +void grt_source_coef(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]) { // 先全部赋0 - for(int i=0; iisrc; + cplx_t xa = mod1d->xa[isrc]; + cplx_t caca = mod1d->caca[isrc]; cplx_t xb = mod1d->xb[isrc]; cplx_t cbcb = mod1d->cbcb[isrc]; real_t k = mod1d->k; - // cplx_t b = k*src_xb; - cplx_t tmp; - - // 水平力源 (4.6.21,26), 这里可以把x1,x2方向的力转到r,theta方向 - // 推导可发现，r方向的力形成P,SV波, theta方向的力形成SH波 - // 方向性因子包含水平力方向与震源台站连线方向的夹角 - coef[2][0] = tmp = cbcb / xb; coef[2][1] = tmp; + _source_PSV(xa, caca, xb, cbcb, k, coef); +} - // 剪切位错 (4.8.34) - // m=1 - coef[4][0] = tmp = - cbcb * k; coef[4][1] = - tmp; - // m=2 - coef[5][0] = tmp = (cbcb / xb) * k; coef[5][1] = tmp; +void grt_source_coef_SH(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]) +{ + size_t isrc = mod1d->isrc; + cplx_t xb = mod1d->xb[isrc]; + cplx_t cbcb = mod1d->cbcb[isrc]; + real_t k = mod1d->k; + _source_SH(xb, cbcb, k, coef); } - - diff --git a/pygrt/C_extension/src/static/static_layer.c b/pygrt/C_extension/src/static/static_layer.c index f71defaf..a8779c63 100644 --- a/pygrt/C_extension/src/static/static_layer.c +++ b/pygrt/C_extension/src/static/static_layer.c @@ -52,6 +52,7 @@ void grt_static_wave2qwv_REV_PSV(const GRT_MODEL1D *mod1d, const cplx_t R[2][2], void grt_static_wave2qwv_REV_SH(const GRT_MODEL1D *mod1d, cplx_t RL, cplx_t *R_EVL) { + (void)mod1d; // 暂停该参数 -Wunused-parameter 警告 *R_EVL = (1.0 + (RL)); } diff --git a/pygrt/C_extension/src/static/static_source.c b/pygrt/C_extension/src/static/static_source.c index 07a64cec..a1136b72 100644 --- a/pygrt/C_extension/src/static/static_source.c +++ b/pygrt/C_extension/src/static/static_source.c @@ -12,25 +12,13 @@ #include #include +#include #include "grt/static/static_source.h" -#include "grt/common/const.h" -void grt_static_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM-1][2]) -{ - // 先全部赋0 - for(int i=0; iisrc; - cplx_t delta = mod1d->delta[isrc]; - real_t k = mod1d->k; +inline void _source_PSV(const cplx_t delta, const real_t k, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]) +{ cplx_t tmp; cplx_t A = 1.0+delta; @@ -39,10 +27,10 @@ void grt_static_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_ // 垂直力源 coef[1][0][0] = tmp = -1.0/(2.0*A*k); coef[1][0][1] = - tmp; - coef[1][1][0] = tmp; coef[1][1][1] = - tmp; + coef[1][1][0] = tmp; coef[1][1][1] = - tmp; // 水平力源 - coef[2][0][0] = tmp = 1.0/(2.0*A*k); coef[2][0][1] = tmp; + coef[2][0][0] = tmp = 1.0/(2.0*A*k); coef[2][0][1] = tmp; coef[2][1][0] = - tmp; coef[2][1][1] = - tmp; // 剪切位错 @@ -50,7 +38,7 @@ void grt_static_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_ coef[3][0][0] = tmp = (-1.0+4.0*delta)/(2.0*A); coef[3][0][1] = tmp; coef[3][1][0] = tmp = -3.0/(2.0*A); coef[3][1][1] = tmp; // m=1 - coef[4][0][0] = tmp = -delta/A; coef[4][0][1] = -tmp; + coef[4][0][0] = tmp = -delta/A; coef[4][0][1] = -tmp; coef[4][1][0] = tmp = 1.0/A; coef[4][1][1] = -tmp; // m=2 coef[5][0][0] = tmp = 1.0/(2.0*A); coef[5][0][1] = tmp; @@ -58,27 +46,46 @@ void grt_static_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_ } -void grt_static_source_coef_SH(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][2]) +inline void _source_SH(const real_t k, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]) { - real_t k = mod1d->k; - - // 先全部赋0 - for(int i=0; iisrc; + cplx_t delta = mod1d->delta[isrc]; + real_t k = mod1d->k; + + _source_PSV(delta, k, coef); +} + + +void grt_static_source_coef_SH(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]) +{ + real_t k = mod1d->k; + + _source_SH(k, coef); }