diff --git a/pygrt/C_extension/include/common/const.h b/pygrt/C_extension/include/common/const.h index e08233e9..ced9c3c6 100755 --- a/pygrt/C_extension/include/common/const.h +++ b/pygrt/C_extension/include/common/const.h @@ -70,6 +70,8 @@ typedef int MYINT; ///< 整数 #define MIN_DEPTH_GAP_SRC_RCV 1.0 ///< 震源和台站的最小深度差（不做绝对限制，仅用于参考波数积分上限） #define GCC_ALWAYS_INLINE __attribute__((always_inline)) ///< gcc编译器不改动内联函数 +#define GRT_SWAP(type, a, b) { type temp = a; a = b; b = temp; } ///< 交换两个变量的值 + #define GRT_STRING_FMT "%18s" ///< 字符串输出格式 #define GRT_REAL_FMT "%18.8e" ///< 浮点数输出格式 #define GRT_CMPLX_FMT "%18.8e%-+14.8ej" ///< 复数输出格式 diff --git a/pygrt/C_extension/include/common/model.h b/pygrt/C_extension/include/common/model.h index daa2b2c5..fa813ea2 100755 --- a/pygrt/C_extension/include/common/model.h +++ b/pygrt/C_extension/include/common/model.h @@ -143,8 +143,20 @@ PYMODEL1D * init_pymod(MYINT n); * @param[in] modelpath 模型文件路径 * @param[in] depsrc 震源深度 * @param[in] deprcv 接收深度 + * @param[in] allowLiquid 是否允许液体层 * * @return `PYMODEL1D` 结构体指针 * */ -PYMODEL1D * read_pymod_from_file(const char *command, const char *modelpath, double depsrc, double deprcv); \ No newline at end of file +PYMODEL1D * read_pymod_from_file(const char *command, const char *modelpath, double depsrc, double deprcv, bool allowLiquid); + + +/** + * 计算PYMODEL1D结构体中的最大最小速度（非零值） + * + * @param pymod (in)`PYMODEL1D` 结构体指针 + * @param vmin (out)最小速度 + * @param vmax (out)最大速度 + * + */ +void get_pymod_vmin_vmax(const PYMODEL1D *pymod, double *vmin, double *vmax); \ No newline at end of file diff --git a/pygrt/C_extension/include/dynamic/layer.h b/pygrt/C_extension/include/dynamic/layer.h index 108d5081..23c4faa6 100755 --- a/pygrt/C_extension/include/dynamic/layer.h +++ b/pygrt/C_extension/include/dynamic/layer.h @@ -89,6 +89,7 @@ void calc_uiz_R_EV( * @param[in] kbkb2 下层的S波水平波数的平方 \f$ k_b^2=(\frac{\omega}{V_b})^2 \f$ * @param[in] mu2 下层的剪切模量 * @param[in] thk 上层层厚 + * @param[in] omega 角频率 * @param[in] k 波数 * @param[out] RD P-SV 下传反射系数矩阵 * @param[out] RDL SH 下传反射系数 @@ -104,11 +105,33 @@ void calc_uiz_R_EV( void calc_RT_2x2( MYREAL Rho1, MYCOMPLEX xa1, MYCOMPLEX xb1, MYCOMPLEX kbkb1, MYCOMPLEX mu1, MYREAL Rho2, MYCOMPLEX xa2, MYCOMPLEX xb2, MYCOMPLEX kbkb2, MYCOMPLEX mu2, - MYREAL thk, MYREAL k, + MYREAL thk, MYCOMPLEX omega, MYREAL k, MYCOMPLEX RD[2][2], MYCOMPLEX *RDL, MYCOMPLEX RU[2][2], MYCOMPLEX *RUL, MYCOMPLEX TD[2][2], MYCOMPLEX *TDL, MYCOMPLEX TU[2][2], MYCOMPLEX *TUL, MYINT *stats); +/** 固-固 界面，函数参数与 calc_RT_2x2 函数一致 */ +void calc_RT_ss_2x2( + MYREAL Rho1, MYCOMPLEX xa1, MYCOMPLEX xb1, MYCOMPLEX kbkb1, MYCOMPLEX mu1, + MYREAL Rho2, MYCOMPLEX xa2, MYCOMPLEX xb2, MYCOMPLEX kbkb2, MYCOMPLEX mu2, + MYREAL thk, MYCOMPLEX omega, MYREAL k, + MYCOMPLEX RD[2][2], MYCOMPLEX *RDL, MYCOMPLEX RU[2][2], MYCOMPLEX *RUL, + MYCOMPLEX TD[2][2], MYCOMPLEX *TDL, MYCOMPLEX TU[2][2], MYCOMPLEX *TUL, MYINT *stats); +/** 液-液 界面，函数参数与 calc_RT_2x2 函数一致 */ +void calc_RT_ll_2x2( + MYREAL Rho1, MYCOMPLEX xa1, MYCOMPLEX xb1, MYCOMPLEX kbkb1, MYCOMPLEX mu1, + MYREAL Rho2, MYCOMPLEX xa2, MYCOMPLEX xb2, MYCOMPLEX kbkb2, MYCOMPLEX mu2, + MYREAL thk, MYCOMPLEX omega, MYREAL k, + MYCOMPLEX RD[2][2], MYCOMPLEX *RDL, MYCOMPLEX RU[2][2], MYCOMPLEX *RUL, + MYCOMPLEX TD[2][2], MYCOMPLEX *TDL, MYCOMPLEX TU[2][2], MYCOMPLEX *TUL, MYINT *stats); + +/** 液-固 界面，函数参数与 calc_RT_2x2 函数一致 */ +void calc_RT_ls_2x2( + MYREAL Rho1, MYCOMPLEX xa1, MYCOMPLEX xb1, MYCOMPLEX kbkb1, MYCOMPLEX mu1, + MYREAL Rho2, MYCOMPLEX xa2, MYCOMPLEX xb2, MYCOMPLEX kbkb2, MYCOMPLEX mu2, + MYREAL thk, MYCOMPLEX omega, MYREAL k, + MYCOMPLEX RD[2][2], MYCOMPLEX *RDL, MYCOMPLEX RU[2][2], MYCOMPLEX *RUL, + MYCOMPLEX TD[2][2], MYCOMPLEX *TDL, MYCOMPLEX TU[2][2], MYCOMPLEX *TUL, MYINT *stats); /** * 【未使用，仅用于代码测试】 diff --git a/pygrt/C_extension/src/common/model.c b/pygrt/C_extension/src/common/model.c index a4076f20..1a0aec9b 100755 --- a/pygrt/C_extension/src/common/model.c +++ b/pygrt/C_extension/src/common/model.c @@ -214,7 +214,7 @@ void update_mod1d_omega(MODEL1D *mod1d, MYCOMPLEX omega){ atnb = attenuation_law(lay->Qbinv, omega); ka0 = omega/(Va0*atna); - kb0 = omega/(Vb0*atnb); + kb0 = (Vb0>RZERO)? omega/(Vb0*atnb) : CZERO; lay->kaka = ka0*ka0; lay->kbkb = kb0*kb0; @@ -255,7 +255,7 @@ void realloc_pymod(PYMODEL1D *pymod, MYINT n){ } -PYMODEL1D * read_pymod_from_file(const char *command, const char *modelpath, double depsrc, double deprcv){ +PYMODEL1D * read_pymod_from_file(const char *command, const char *modelpath, double depsrc, double deprcv, bool allowLiquid){ FILE *fp; if((fp = fopen(modelpath, "r")) == NULL){ fprintf(stderr, "[%s] " BOLD_RED "Model file open error.\n" DEFAULT_RESTORE, command); @@ -303,8 +303,13 @@ PYMODEL1D * read_pymod_from_file(const char *command, const char *modelpath, dou return NULL; }; - if(va <= 0.0 || vb <= 0.0 || rho <= 0.0 || qa <= 0.0 || qb <= 0.0){ - fprintf(stderr, "[%s] " BOLD_RED "In line %d, nonpositive value is not supported.\n" DEFAULT_RESTORE, command, iline); + if(va <= 0.0 || rho <= 0.0 || qa <= 0.0 || qb <= 0.0){ + fprintf(stderr, "[%s] " BOLD_RED "In model file, line %d, nonpositive value is not supported.\n" DEFAULT_RESTORE, command, iline); + return NULL; + } + + if(!allowLiquid && vb <=0){ + fprintf(stderr, "[%s] " BOLD_RED "In model file, line %d, nonpositive Vs is not supported.\n" DEFAULT_RESTORE, command, iline); return NULL; } @@ -391,4 +396,18 @@ PYMODEL1D * read_pymod_from_file(const char *command, const char *modelpath, dou free(modarr); return pymod; +} + + +void get_pymod_vmin_vmax(const PYMODEL1D *pymod, double *vmin, double *vmax){ + *vmin = __DBL_MAX__; + *vmax = RZERO; + const MYREAL *Va = pymod->Va; + const MYREAL *Vb = pymod->Vb; + for(MYINT i=0; in; ++i){ + if(Va[i] < *vmin) *vmin = Va[i]; + if(Va[i] > *vmax) *vmax = Va[i]; + if(Vb[i] < *vmin && Vb[i] > RZERO) *vmin = Vb[i]; + if(Vb[i] > *vmax && Vb[i] > RZERO) *vmax = Vb[i]; + } } \ No newline at end of file diff --git a/pygrt/C_extension/src/dynamic/grt_main.c b/pygrt/C_extension/src/dynamic/grt_main.c index 0eb07f82..3fca6ca9 100644 --- a/pygrt/C_extension/src/dynamic/grt_main.c +++ b/pygrt/C_extension/src/dynamic/grt_main.c @@ -991,17 +991,12 @@ int main(int argc, char **argv) { getopt_from_command(argc, argv); // 读入模型文件 - if((pymod = read_pymod_from_file(command, s_modelpath, depsrc, deprcv)) ==NULL){ + if((pymod = read_pymod_from_file(command, s_modelpath, depsrc, deprcv, true)) ==NULL){ exit(EXIT_FAILURE); } // 最大最小速度 - vmax = pymod->Va[findMinMax_MYREAL(pymod->Va, pymod->n, true)]; - vmin = pymod->Vb[findMinMax_MYREAL(pymod->Vb, pymod->n, false)]; - if(vmin > vmax) { - double tmp; - tmp = vmin; vmin = vmax; vmax = tmp; - } + get_pymod_vmin_vmax(pymod, &vmin, &vmax); // 参考最小速度 if(vmin_ref == 0.0){ diff --git a/pygrt/C_extension/src/dynamic/layer.c b/pygrt/C_extension/src/dynamic/layer.c index 9a12e90a..cf8df9c3 100755 --- a/pygrt/C_extension/src/dynamic/layer.c +++ b/pygrt/C_extension/src/dynamic/layer.c @@ -12,35 +12,44 @@ */ #include +#include #include #include "dynamic/layer.h" #include "common/model.h" #include "common/prtdbg.h" #include "common/matrix.h" +#include "common/colorstr.h" void calc_R_tilt( MYCOMPLEX xa0, MYCOMPLEX xb0, MYCOMPLEX kbkb0, MYREAL k, MYCOMPLEX R_tilt[2][2], MYINT *stats) { - - // // 公式(5.3.10-14) - MYCOMPLEX Delta = RZERO; - MYREAL kk = k*k; - MYCOMPLEX kbkb_k2inv = kbkb0/kk; - MYCOMPLEX kbkb_k4inv = RQUART*kbkb_k2inv*kbkb_k2inv; - - // 对公式(5.3.10-14)进行重新整理，对浮点数友好一些 - Delta = -RONE + xa0*xb0 + kbkb_k2inv - kbkb_k4inv; - if(Delta == CZERO){ - *stats = INVERSE_FAILURE; - return; + if(kbkb0 != CZERO){ + // 固体表面 + // 公式(5.3.10-14) + MYCOMPLEX Delta = RZERO; + MYREAL kk = k*k; + MYCOMPLEX kbkb_k2inv = kbkb0/kk; + MYCOMPLEX kbkb_k4inv = RQUART*kbkb_k2inv*kbkb_k2inv; + + // 对公式(5.3.10-14)进行重新整理，对浮点数友好一些 + Delta = -RONE + xa0*xb0 + kbkb_k2inv - kbkb_k4inv; + if(Delta == CZERO){ + *stats = INVERSE_FAILURE; + return; + } + R_tilt[0][0] = (RONE + xa0*xb0 - kbkb_k2inv + kbkb_k4inv) / Delta; + R_tilt[0][1] = RTWO * xb0 * (RONE - RHALF*kbkb_k2inv) / Delta; + R_tilt[1][0] = RTWO * xa0 * (RONE - RHALF*kbkb_k2inv) / Delta; + R_tilt[1][1] = R_tilt[0][0]; + } + else { + // 液体表面 + R_tilt[0][0] = -RONE; + R_tilt[1][1] = R_tilt[0][1] = R_tilt[1][0] = CZERO; } - R_tilt[0][0] = (RONE + xa0*xb0 - kbkb_k2inv + kbkb_k4inv) / Delta; - R_tilt[0][1] = RTWO * xb0 * (RONE - RHALF*kbkb_k2inv) / Delta; - R_tilt[1][0] = RTWO * xa0 * (RONE - RHALF*kbkb_k2inv) / Delta; - R_tilt[1][1] = R_tilt[0][0]; } @@ -51,9 +60,24 @@ void calc_R_EV( const MYCOMPLEX R[2][2], MYCOMPLEX RL, MYCOMPLEX R_EV[2][2], MYCOMPLEX *R_EVL) { - // 公式(5.2.19) - MYCOMPLEX D11[2][2] = {{k, k*xb_rcv}, {k*xa_rcv, k}}; - MYCOMPLEX D12[2][2] = {{k, -k*xb_rcv}, {-k*xa_rcv, k}}; + MYCOMPLEX D11[2][2], D12[2][2]; + if(xb_rcv != CONE){ + // 位于固体层 + // 公式(5.2.19) + D11[0][0] = k; D11[0][1] = k*xb_rcv; + D11[1][0] = k*xa_rcv; D11[1][1] = k; + D12[0][0] = k; D12[0][1] = -k*xb_rcv; + D12[1][0] = -k*xa_rcv; D12[1][1] = k; + *R_EVL = (RONE + (RL))*k; + } else { + // 位于液体层 + D11[0][0] = k; D11[0][1] = RZERO; + D11[1][0] = k*xa_rcv; D11[1][1] = RZERO; + D12[0][0] = k; D12[0][1] = RZERO; + D12[1][0] = -k*xa_rcv; D12[1][1] = RZERO; + *R_EVL = RZERO; + } + // 公式(5.7.7,25) if(ircvup){// 震源更深 @@ -63,7 +87,6 @@ void calc_R_EV( cmat2x2_mul(D11, R, R_EV); cmat2x2_add(D12, R_EV, R_EV); } - *R_EVL = (RONE + (RL))*k; } @@ -83,26 +106,192 @@ void calc_uiz_R_EV( MYCOMPLEX D11[2][2] = {{ak, bb}, {aa, bk}}; MYCOMPLEX D12[2][2] = {{-ak, bb}, {aa, -bk}}; + if(ircvup){// 震源更深 + *R_EVL = (RONE - (RL))*bk; + } else { // 接收点更深 + *R_EVL = (RL - RONE)*bk; + } + + // 位于液体层 + if(xb_rcv == CONE){ + D11[0][1] = D11[1][1] = D12[0][1] = D12[1][1] = *R_EVL = RZERO; + } + // 公式(5.7.7,25) if(ircvup){// 震源更深 cmat2x2_mul(D12, R, R_EV); cmat2x2_add(D11, R_EV, R_EV); - *R_EVL = (RONE - (RL))*bk; } else { // 接收点更深 cmat2x2_mul(D11, R, R_EV); cmat2x2_add(D12, R_EV, R_EV); - *R_EVL = (RL - RONE)*bk; } } -void calc_RT_2x2( +void calc_RT_ll_2x2( MYREAL Rho1, MYCOMPLEX xa1, MYCOMPLEX xb1, MYCOMPLEX kbkb1, MYCOMPLEX mu1, MYREAL Rho2, MYCOMPLEX xa2, MYCOMPLEX xb2, MYCOMPLEX kbkb2, MYCOMPLEX mu2, MYREAL thk, // 使用上层的厚度 - MYREAL k, + MYCOMPLEX omega, MYREAL k, + MYCOMPLEX RD[2][2], MYCOMPLEX *RDL, MYCOMPLEX RU[2][2], MYCOMPLEX *RUL, + MYCOMPLEX TD[2][2], MYCOMPLEX *TDL, MYCOMPLEX TU[2][2], MYCOMPLEX *TUL, MYINT *stats) +{ + MYCOMPLEX exa, ex2a; + + exa = exp(-k*thk*xa1); + ex2a = exa * exa; + + + bool computeRayl = true; + bool computeLove = true; + if(RD==NULL || RU==NULL || TD==NULL || TU==NULL) computeRayl=false; + if(RDL==NULL || RUL==NULL || TDL==NULL || TUL==NULL) computeLove=false; + + MYCOMPLEX A = xa1*Rho2 + xa2*Rho1; + + if(computeRayl){ + RD[0][0] = (xa1*Rho2 - xa2*Rho1)/A * ex2a; + RD[0][1] = RD[1][0] = RD[1][1] = CZERO; + + RU[0][0] = (xa2*Rho1 - xa1*Rho2)/A; + RU[0][1] = RU[1][0] = RU[1][1] = CZERO; + + TD[0][0] = RTWO*xa1*Rho1/A * exa; + TD[0][1] = TD[1][0] = TD[1][1] = CZERO; + + TU[0][0] = RTWO*xa2*Rho2/A * exa; + TU[0][1] = TU[1][0] = TU[1][1] = CZERO; + } + + if(computeLove){ + *RDL = RZERO; + *RUL = RZERO; + *TDL = RZERO; + *TUL = RZERO; + } +} + + +void calc_RT_ls_2x2( + MYREAL Rho1, MYCOMPLEX xa1, MYCOMPLEX xb1, MYCOMPLEX kbkb1, MYCOMPLEX mu1, + MYREAL Rho2, MYCOMPLEX xa2, MYCOMPLEX xb2, MYCOMPLEX kbkb2, MYCOMPLEX mu2, + MYREAL thk, // 使用上层的厚度 + MYCOMPLEX omega, MYREAL k, + MYCOMPLEX RD[2][2], MYCOMPLEX *RDL, MYCOMPLEX RU[2][2], MYCOMPLEX *RUL, + MYCOMPLEX TD[2][2], MYCOMPLEX *TDL, MYCOMPLEX TU[2][2], MYCOMPLEX *TUL, MYINT *stats) +{ + // 后缀1表示上层的液体的物理参数，后缀2表示下层的固体的物理参数 + // 若mu2==0, 则下层为液体，参数需相互交换 + + // 延迟因子始终作用于上层 + MYCOMPLEX exa, exb, exab, ex2a, ex2b; + MYCOMPLEX tmp; + + exa = exp(-k*thk*xa1); + exb = exp(-k*thk*xb1); + + exab = exa * exb; + ex2a = exa * exa; + ex2b = exb * exb; + + bool computeRayl = true; + bool computeLove = true; + if(RD==NULL || RU==NULL || TD==NULL || TU==NULL) computeRayl=false; + if(RDL==NULL || RUL==NULL || TDL==NULL || TUL==NULL) computeLove=false; + + // 讨论液-固 or 固-液 + bool isfluidUp = (mu1 == CZERO); // 上层是否为液体 + MYINT sgn = 1; + if(isfluidUp && mu2 == CZERO){ + fprintf(stderr, BOLD_RED "Error: fluid-fluid interface is not allowed in calc_RT_fs_2x2\n" DEFAULT_RESTORE); + exit(EXIT_FAILURE); + } + + + // 使用指针 + MYCOMPLEX (*pRD)[2], *pRDL, (*pRU)[2], *pRUL; + MYCOMPLEX (*pTD)[2], *pTDL, (*pTU)[2], *pTUL; + if(isfluidUp){ + pRD = RD; pRDL = RDL; pRU = RU; pRUL = RUL; + pTD = TD; pTDL = TDL; pTU = TU; pTUL = TUL; + } else { + pRD = RU; pRDL = RUL; pRU = RD; pRUL = RDL; + pTD = TU; pTDL = TUL; pTU = TD; pTUL = TDL; + GRT_SWAP(MYREAL, Rho1, Rho2); + GRT_SWAP(MYCOMPLEX, xa1, xa2); + GRT_SWAP(MYCOMPLEX, xb1, xb2); + GRT_SWAP(MYCOMPLEX, kbkb1, kbkb2); + GRT_SWAP(MYCOMPLEX, mu1, mu2); + sgn = -1; + } + + + // 定义一些中间变量来简化运算和书写 + MYREAL k2 = k*k; + MYCOMPLEX lamka1k = Rho1*omega*omega/k2; + MYCOMPLEX kb2k = kbkb2/k2; + MYCOMPLEX Og2k = RONE - RHALF*kb2k; + MYCOMPLEX Og2k2 = Og2k*Og2k; + MYCOMPLEX A = RTWO*Og2k2*xa1*mu2 + RHALF*lamka1k*kb2k*xa2 - RTWO*mu2*xa1*xa2*xb2; + MYCOMPLEX B = RTWO*Og2k2*xa1*mu2 - RHALF*lamka1k*kb2k*xa2 + RTWO*mu2*xa1*xa2*xb2; + MYCOMPLEX C = RTWO*Og2k2*xa1*mu2 + RHALF*lamka1k*kb2k*xa2 + RTWO*mu2*xa1*xa2*xb2; + MYCOMPLEX D = RTWO*Og2k2*xa1*mu2 - RHALF*lamka1k*kb2k*xa2 - RTWO*mu2*xa1*xa2*xb2; + + if(A == CZERO){ + *stats = INVERSE_FAILURE; + return; + } + + // 按液体层在上层处理 + if(computeRayl){ + pRD[0][0] = D/A; + pRD[0][1] = pRD[1][0] = pRD[1][1] = CZERO; + + pRU[0][0] = - B/A; + pRU[0][1] = - RFOUR*Og2k*xa1*xb2*mu2/A * sgn; + pRU[1][0] = pRU[0][1]/xb2 * xa2; + pRU[1][1] = - C/A; + + pTD[0][0] = - RTWO*Og2k*xa1*lamka1k/A; pTD[0][1] = CZERO; + pTD[1][0] = pTD[0][0]/Og2k*xa2 * sgn; pTD[1][1] = CZERO; + + pTU[0][0] = - RTWO*Og2k*xa2*mu2*kb2k/A; pTU[0][1] = pTU[0][0]/Og2k*xb2 * sgn; + pTU[1][0] = pTU[1][1] = CZERO; + + // 回归数组，增加时移 + RD[0][0] *= ex2a; RD[0][1] *= exab; + RD[1][0] *= exab; RD[1][1] *= ex2b; + + TD[0][0] *= exa; TD[0][1] *= exb; + TD[1][0] *= exa; TD[1][1] *= exb; + + TU[0][0] *= exa; TU[0][1] *= exa; + TU[1][0] *= exb; TU[1][1] *= exb; + } + + if(computeLove){ + *pRDL = RZERO; + *pRUL = RONE; + *pTDL = RZERO; + *pTUL = RZERO; + + // 回归数组，增加时移 + // 特殊数值，仅需为潜在的RDL添加 + *RDL *= ex2b; + } + + +} + + + +void calc_RT_ss_2x2( + MYREAL Rho1, MYCOMPLEX xa1, MYCOMPLEX xb1, MYCOMPLEX kbkb1, MYCOMPLEX mu1, + MYREAL Rho2, MYCOMPLEX xa2, MYCOMPLEX xb2, MYCOMPLEX kbkb2, MYCOMPLEX mu2, + MYREAL thk, // 使用上层的厚度 + MYCOMPLEX omega, MYREAL k, MYCOMPLEX RD[2][2], MYCOMPLEX *RDL, MYCOMPLEX RU[2][2], MYCOMPLEX *RUL, MYCOMPLEX TD[2][2], MYCOMPLEX *TDL, MYCOMPLEX TU[2][2], MYCOMPLEX *TUL, MYINT *stats) { @@ -205,6 +394,39 @@ void calc_RT_2x2( } +void calc_RT_2x2( + MYREAL Rho1, MYCOMPLEX xa1, MYCOMPLEX xb1, MYCOMPLEX kbkb1, MYCOMPLEX mu1, + MYREAL Rho2, MYCOMPLEX xa2, MYCOMPLEX xb2, MYCOMPLEX kbkb2, MYCOMPLEX mu2, + MYREAL thk, // 使用上层的厚度 + MYCOMPLEX omega, MYREAL k, + MYCOMPLEX RD[2][2], MYCOMPLEX *RDL, MYCOMPLEX RU[2][2], MYCOMPLEX *RUL, + MYCOMPLEX TD[2][2], MYCOMPLEX *TDL, MYCOMPLEX TU[2][2], MYCOMPLEX *TUL, MYINT *stats) +{ + // 根据界面两侧的具体情况选择函数 + if(mu1 != CZERO && mu2 != CZERO){ + calc_RT_ss_2x2( + Rho1, xa1, xb1, kbkb1, mu1, + Rho2, xa2, xb2, kbkb2, mu2, + thk, omega, k, + RD, RDL, RU, RUL, TD, TDL, TU, TUL, stats); + } + else if(mu1 == CZERO && mu2 == CZERO){ + calc_RT_ll_2x2( + Rho1, xa1, xb1, kbkb1, mu1, + Rho2, xa2, xb2, kbkb2, mu2, + thk, omega, k, + RD, RDL, RU, RUL, TD, TDL, TU, TUL, stats); + } + else{ + calc_RT_ls_2x2( + Rho1, xa1, xb1, kbkb1, mu1, + Rho2, xa2, xb2, kbkb2, mu2, + thk, omega, k, + RD, RDL, RU, RUL, TD, TDL, TU, TUL, stats); + } +} + + void get_layer_D( MYCOMPLEX xa, MYCOMPLEX xb, MYCOMPLEX kbkb, MYCOMPLEX mu, diff --git a/pygrt/C_extension/src/dynamic/propagate.c b/pygrt/C_extension/src/dynamic/propagate.c index 115fdfda..26e6f91c 100755 --- a/pygrt/C_extension/src/dynamic/propagate.c +++ b/pygrt/C_extension/src/dynamic/propagate.c @@ -194,7 +194,7 @@ void kernel( mod1d_Rho0, mod1d_xa0, mod1d_xb0, mod1d_kbkb0, mod1d_mu0, mod1d_Rho1, mod1d_xa1, mod1d_xb1, mod1d_kbkb1, mod1d_mu1, mod1d_thk0, // 使用iy-1层的厚度 - k, + omega, k, RD, pRDL, RU, pRUL, TD, pTDL, TU, pTUL, stats); if(*stats==INVERSE_FAILURE) goto BEFORE_RETURN; diff --git a/pygrt/C_extension/src/static/stgrt_main.c b/pygrt/C_extension/src/static/stgrt_main.c index bcb95f25..5c9169c8 100644 --- a/pygrt/C_extension/src/static/stgrt_main.c +++ b/pygrt/C_extension/src/static/stgrt_main.c @@ -476,17 +476,12 @@ int main(int argc, char **argv){ getopt_from_command(argc, argv); // 读入模型文件 - if((pymod = read_pymod_from_file(command, s_modelpath, depsrc, deprcv)) ==NULL){ + if((pymod = read_pymod_from_file(command, s_modelpath, depsrc, deprcv, false)) ==NULL){ exit(EXIT_FAILURE); } // 最大最小速度 - vmax = pymod->Va[findMinMax_MYREAL(pymod->Va, pymod->n, true)]; - vmin = pymod->Vb[findMinMax_MYREAL(pymod->Vb, pymod->n, false)]; - if(vmin > vmax) { - double tmp; - tmp = vmin; vmin = vmax; vmax = tmp; - } + get_pymod_vmin_vmax(pymod, &vmin, &vmax); // 参考最小速度 if(vmin_ref == 0.0){ diff --git a/pygrt/C_extension/src/travt/grt_travt.c b/pygrt/C_extension/src/travt/grt_travt.c index 6d29e1d7..593a2083 100644 --- a/pygrt/C_extension/src/travt/grt_travt.c +++ b/pygrt/C_extension/src/travt/grt_travt.c @@ -205,7 +205,7 @@ int main(int argc, char **argv){ getopt_from_command(argc, argv); // 读入模型文件 - if((pymod = read_pymod_from_file(command, s_modelpath, depsrc, deprcv)) ==NULL){ + if((pymod = read_pymod_from_file(command, s_modelpath, depsrc, deprcv, true)) ==NULL){ exit(EXIT_FAILURE); } // print_pymod(pymod); diff --git a/pygrt/C_extension/src/travt/travt.c b/pygrt/C_extension/src/travt/travt.c index c43eb105..01e852ba 100644 --- a/pygrt/C_extension/src/travt/travt.c +++ b/pygrt/C_extension/src/travt/travt.c @@ -17,9 +17,30 @@ MYREAL compute_travt1d( - const MYREAL *Thk, const MYREAL *Vel, const int nlay, + const MYREAL *Thk, const MYREAL *Vel0, const int nlay, const int isrc, const int ircv, const MYREAL dist) { + // 以防速度数组中存在零速度的情况，这里新建数组以去除0速度 + MYREAL *Vel = (MYREAL*)malloc(sizeof(MYREAL)*nlay); + for(int i=0; i