Skip to content

REFAC: Split source coefs into P-SV and SH #67

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Aug 20, 2025
Merged

REFAC: Split source coefs into P-SV and SH #67

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 4 additions & 2 deletions pygrt/C_extension/include/common/recursion.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]);
12 changes: 9 additions & 3 deletions pygrt/C_extension/include/dynamic/source.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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$
Expand All @@ -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]);
9 changes: 5 additions & 4 deletions pygrt/C_extension/include/static/static_source.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]);

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]);
11 changes: 6 additions & 5 deletions pygrt/C_extension/src/common/recursion.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);

Expand Down
14 changes: 8 additions & 6 deletions pygrt/C_extension/src/dynamic/propagate.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -357,7 +359,7 @@ void kernel(
tmpRL2 = R_EVL * tmpRL;

for(MYINT i=0; i<SRC_M_NUM; ++i){
get_qwv(ircvup, tmp2x2, tmpRL2, RD_BL, RDL_BL, src_coef[i], QWV[i]);
get_qwv(ircvup, tmp2x2, tmpRL2, RD_BL, RDL_BL, src_coef_PSV[i], src_coef_SH[i], QWV[i]);
}


Expand All @@ -368,7 +370,7 @@ void kernel(
tmpRL2 = uiz_R_EVL * tmpRL;

for(MYINT i=0; i<SRC_M_NUM; ++i){
get_qwv(ircvup, tmp2x2_uiz, tmpRL2, RD_BL, RDL_BL, src_coef[i], QWV_uiz[i]);
get_qwv(ircvup, tmp2x2_uiz, tmpRL2, RD_BL, RDL_BL, src_coef_PSV[i], src_coef_SH[i], QWV_uiz[i]);
}
}
}
Expand Down Expand Up @@ -402,7 +404,7 @@ void kernel(
tmpRL2 = R_EVL * tmpRL;

for(MYINT i=0; i<SRC_M_NUM; ++i){
get_qwv(ircvup, tmp2x2, tmpRL2, RU_FA, RUL_FA, src_coef[i], QWV[i]);
get_qwv(ircvup, tmp2x2, tmpRL2, RU_FA, RUL_FA, src_coef_PSV[i], src_coef_SH[i], QWV[i]);
}


Expand All @@ -413,7 +415,7 @@ void kernel(
tmpRL2 = uiz_R_EVL * tmpRL;

for(MYINT i=0; i<SRC_M_NUM; ++i){
get_qwv(ircvup, tmp2x2_uiz, tmpRL2, RU_FA, RUL_FA, src_coef[i], QWV_uiz[i]);
get_qwv(ircvup, tmp2x2_uiz, tmpRL2, RU_FA, RUL_FA, src_coef_PSV[i], src_coef_SH[i], QWV_uiz[i]);
}
}

Expand Down
42 changes: 33 additions & 9 deletions pygrt/C_extension/src/dynamic/source.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,23 +18,20 @@
#include "common/prtdbg.h"




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])
{
// 先全部赋0
for(MYINT i=0; i<SRC_M_NUM; ++i){
for(MYINT j=0; j<QWV_NUM; ++j){
for(MYINT j=0; j<QWV_NUM-1; ++j){
for(MYINT p=0; p<2; ++p){
coef[i][j][p] = CZERO;
}
}
}


MYCOMPLEX a = k*src_xa;
MYCOMPLEX b = k*src_xb;
MYREAL kk = k*k;
Expand All @@ -53,7 +50,6 @@ void source_coef(
// 方向性因子包含水平力方向与震源台站连线方向的夹角
coef[2][0][0] = tmp = -k / a; coef[2][0][1] = tmp;
coef[2][1][0] = tmp = -RONE; coef[2][1][1] = - tmp;
coef[2][2][0] = tmp = src_kbkb / k / b; coef[2][2][1] = tmp;

// 剪切位错 (4.8.34)
// m=0
Expand All @@ -62,13 +58,41 @@ void source_coef(
// m=1
coef[4][0][0] = tmp = RTWO*k; coef[4][0][1] = - tmp;
coef[4][1][0] = tmp = (RTWO*kk - src_kbkb) / b; coef[4][1][1] = tmp;
coef[4][2][0] = tmp = - src_kbkb / k; coef[4][2][1] = - tmp;

// m=2
coef[5][0][0] = tmp = - kk / a; coef[5][0][1] = tmp;
coef[5][1][0] = tmp = - k; coef[5][1][1] = - tmp;
coef[5][2][0] = tmp = src_kbkb / b; coef[5][2][1] = tmp;

}


void source_coef_SH(
MYCOMPLEX src_xb, MYCOMPLEX src_kbkb,
MYREAL k,
MYCOMPLEX coef[SRC_M_NUM][2])
{
// 先全部赋0
for(MYINT i=0; i<SRC_M_NUM; ++i){
for(MYINT p=0; p<2; ++p){
coef[i][p] = CZERO;
}
}


MYCOMPLEX b = k*src_xb;
MYCOMPLEX tmp;

// 水平力源 (4.6.21,26), 这里可以把x1,x2方向的力转到r,theta方向
// 推导可发现，r方向的力形成P,SV波, theta方向的力形成SH波
// 方向性因子包含水平力方向与震源台站连线方向的夹角
coef[2][0] = tmp = src_kbkb / k / b; coef[2][1] = tmp;

// 剪切位错 (4.8.34)
// m=1
coef[4][0] = tmp = - src_kbkb / k; coef[4][1] = - tmp;

// m=2
coef[5][0] = tmp = src_kbkb / b; coef[5][1] = tmp;

}

Expand Down
14 changes: 8 additions & 6 deletions pygrt/C_extension/src/static/static_propagate.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -193,8 +193,10 @@ void static_kernel(


// 计算震源系数
MYCOMPLEX src_coef[SRC_M_NUM][QWV_NUM][2] = {0};
static_source_coef(src_delta, k, src_coef);
MYCOMPLEX src_coef_PSV[SRC_M_NUM][QWV_NUM-1][2] = {0};
MYCOMPLEX src_coef_SH[SRC_M_NUM][2] = {0};
static_source_coef_PSV(src_delta, k, src_coef_PSV);
static_source_coef_SH(k, src_coef_SH);

// 临时中转矩阵 (temperary)
MYCOMPLEX tmpR2[2][2], tmp2x2[2][2], tmpRL, tmp2x2_uiz[2][2], tmpRL_uiz;
Expand Down Expand Up @@ -237,7 +239,7 @@ void static_kernel(
tmpRL = R_EVL * invT / (RONE - RDL_BL * RUL_FB);

for(MYINT i=0; i<SRC_M_NUM; ++i){
get_qwv(ircvup, tmp2x2, tmpRL, RD_BL, RDL_BL, src_coef[i], QWV[i]);
get_qwv(ircvup, tmp2x2, tmpRL, RD_BL, RDL_BL, src_coef_PSV[i], src_coef_SH[i], QWV[i]);
}


Expand All @@ -248,7 +250,7 @@ void static_kernel(
tmpRL_uiz = tmpRL / R_EVL * uiz_R_EVL;

for(MYINT i=0; i<SRC_M_NUM; ++i){
get_qwv(ircvup, tmp2x2_uiz, tmpRL_uiz, RD_BL, RDL_BL, src_coef[i], QWV_uiz[i]);
get_qwv(ircvup, tmp2x2_uiz, tmpRL_uiz, RD_BL, RDL_BL, src_coef_PSV[i], src_coef_SH[i], QWV_uiz[i]);
}
}
}
Expand Down Expand Up @@ -279,7 +281,7 @@ void static_kernel(
tmpRL = R_EVL * invT / (RONE - RUL_FA * RDL_AL);

for(MYINT i=0; i<SRC_M_NUM; ++i){
get_qwv(ircvup, tmp2x2, tmpRL, RU_FA, RUL_FA, src_coef[i], QWV[i]);
get_qwv(ircvup, tmp2x2, tmpRL, RU_FA, RUL_FA, src_coef_PSV[i], src_coef_SH[i], QWV[i]);
}

if(calc_uiz){
Expand All @@ -289,7 +291,7 @@ void static_kernel(
tmpRL_uiz = tmpRL / R_EVL * uiz_R_EVL;

for(MYINT i=0; i<SRC_M_NUM; ++i){
get_qwv(ircvup, tmp2x2_uiz, tmpRL_uiz, RU_FA, RUL_FA, src_coef[i], QWV_uiz[i]);
get_qwv(ircvup, tmp2x2_uiz, tmpRL_uiz, RU_FA, RUL_FA, src_coef_PSV[i], src_coef_SH[i], QWV_uiz[i]);
}
}
} // END if
Expand Down
31 changes: 24 additions & 7 deletions pygrt/C_extension/src/static/static_source.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,13 +16,11 @@
#include "static/static_source.h"
#include "common/const.h"


void static_source_coef(
MYCOMPLEX delta, MYREAL k, MYCOMPLEX coef[SRC_M_NUM][QWV_NUM][2])
void static_source_coef_PSV(MYCOMPLEX delta, MYREAL k, MYCOMPLEX coef[SRC_M_NUM][QWV_NUM-1][2])
{
// 先全部赋0
for(MYINT i=0; i<SRC_M_NUM; ++i){
for(MYINT j=0; j<QWV_NUM; ++j){
for(MYINT j=0; j<QWV_NUM-1; ++j){
for(MYINT p=0; p<2; ++p){
coef[i][j][p] = CZERO;
}
Expand All @@ -42,7 +40,6 @@ void static_source_coef(
// 水平力源
coef[2][0][0] = tmp = RONE/(RTWO*A*k); coef[2][0][1] = tmp;
coef[2][1][0] = - tmp; coef[2][1][1] = - tmp;
coef[2][2][0] = tmp = -RONE/k; coef[2][2][1] = tmp;

// 剪切位错
// m=0
Expand All @@ -51,11 +48,31 @@ void static_source_coef(
// m=1
coef[4][0][0] = tmp = -delta/A; coef[4][0][1] = -tmp;
coef[4][1][0] = tmp = RONE/A; coef[4][1][1] = -tmp;
coef[4][2][0] = tmp = RONE; coef[4][2][1] = -tmp;
// m=2
coef[5][0][0] = tmp = RONE/(RTWO*A); coef[5][0][1] = tmp;
coef[5][1][0] = tmp = -RONE/(RTWO*A); coef[5][1][1] = tmp;
coef[5][2][0] = tmp = -RONE; coef[5][2][1] = tmp;
}


void static_source_coef_SH(MYREAL k, MYCOMPLEX coef[SRC_M_NUM][2])
{
// 先全部赋0
for(MYINT i=0; i<SRC_M_NUM; ++i){
for(MYINT p=0; p<2; ++p){
coef[i][p] = CZERO;
}
}

MYCOMPLEX tmp;

// 水平力源
coef[2][0] = tmp = -RONE/k; coef[2][1] = tmp;

// 剪切位错
// m=1
coef[4][0] = tmp = RONE; coef[4][1] = -tmp;
// m=2
coef[5][0] = tmp = -RONE; coef[5][1] = tmp;
}