From 6502c0c87931fe61a094ea8f398cc6160c1f0a40 Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Thu, 13 Nov 2025 23:41:52 +0800 Subject: [PATCH] FEAT: update lamb1 interfaces --- pygrt/C_extension/include/grt/lamb/lamb1.h | 5 ++- .../C_extension/include/grt/lamb/lamb_util.h | 12 +------ pygrt/C_extension/src/lamb/grt_lamb1.c | 36 ++++++++----------- pygrt/C_extension/src/lamb/lamb1.c | 17 ++++----- pygrt/C_extension/src/lamb/lamb_util.c | 5 --- pygrt/c_interfaces.py | 2 +- pygrt/utils.py | 11 +++--- 7 files changed, 32 insertions(+), 56 deletions(-) diff --git a/pygrt/C_extension/include/grt/lamb/lamb1.h b/pygrt/C_extension/include/grt/lamb/lamb1.h index 7e65ba6..cbd2f57 100644 --- a/pygrt/C_extension/include/grt/lamb/lamb1.h +++ b/pygrt/C_extension/include/grt/lamb/lamb1.h @@ -17,8 +17,7 @@ /** * 使用广义闭合解求解第一类 Lamb 问题 * - * @param[in] alpha P 波速度 - * @param[in] beta S 波速度 + * @param[in] nu 泊松比， (0, 0.5) * @param[in] ts 归一化时间序列 * @param[in] nt 时间序列点数 * @param[in] azimuth 方位角，单位度 @@ -26,4 +25,4 @@ * */ void grt_solve_lamb1( - const MYREAL alpha, const MYREAL beta, const MYREAL *ts, const int nt, const MYREAL azimuth, MYREAL (*u)[3][3]); \ No newline at end of file + const MYREAL nu, const MYREAL *ts, const int nt, const MYREAL azimuth, MYREAL (*u)[3][3]); \ No newline at end of file diff --git a/pygrt/C_extension/include/grt/lamb/lamb_util.h b/pygrt/C_extension/include/grt/lamb/lamb_util.h index e1b2d57..eca4762 100644 --- a/pygrt/C_extension/include/grt/lamb/lamb_util.h +++ b/pygrt/C_extension/include/grt/lamb/lamb_util.h @@ -31,14 +31,4 @@ void grt_roots3(const MYREAL a, const MYREAL b, const MYREAL c, MYCOMPLEX y3[3]) * @return 多项式结果 * */ -MYCOMPLEX grt_evalpoly2(const MYCOMPLEX *C, const int n, const MYCOMPLEX y, const int offset); - - -/** - * 根据 Vs/Vp 比值，计算泊松比， \f$ \nu = \dfrac{1}{2} \dfrac{1 - 2k^2}{1 - k^2}, k=\dfrac{V_S}{V_P} \f$ - * - * @param[in] k Vs/Vp - * @return 泊松比 - */ - -MYREAL grt_possion_ratio(const MYREAL k); \ No newline at end of file +MYCOMPLEX grt_evalpoly2(const MYCOMPLEX *C, const int n, const MYCOMPLEX y, const int offset); \ No newline at end of file diff --git a/pygrt/C_extension/src/lamb/grt_lamb1.c b/pygrt/C_extension/src/lamb/grt_lamb1.c index efb5fc6..be3615c 100644 --- a/pygrt/C_extension/src/lamb/grt_lamb1.c +++ b/pygrt/C_extension/src/lamb/grt_lamb1.c @@ -21,9 +21,8 @@ typedef struct { /** 模型参数 */ struct { bool active; - MYREAL alpha; ///< P 波速度 - MYREAL beta; ///< S 波速度 - } M; + MYREAL nu; ///< 泊松比 + } P; /** 归一化时间序列 */ struct { @@ -62,11 +61,11 @@ printf("\n" "\n\n" "Usage:\n" "----------------------------------------------------------------\n" -" grt lamb1 -M/ -T//
-A\n" +" grt lamb1 -P -T//
-A\n" "\n\n" "Options:\n" "----------------------------------------------------------------\n" -" -M/ Vp and Vs of the halfspace.\n" +" -P Possion ratio of the halfspace, (0, 0.5).\n" "\n" " -T//
\n" " Dimensionless time.\n" @@ -80,7 +79,7 @@ printf("\n" "\n\n" "Examples:\n" "----------------------------------------------------------------\n" -" grt lamb1 -M8.0/4.62 -T0/2/1e-3 -A30\n" +" grt lamb1 -P0.25 -T0/2/1e-3 -A30\n" "\n\n\n" ); } @@ -91,21 +90,16 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ char* command = Ctrl->name; int opt; - while ((opt = getopt(argc, argv, ":M:T:A:h")) != -1) { + while ((opt = getopt(argc, argv, ":P:T:A:h")) != -1) { switch (opt) { - // 模型参数， -Malpha/beta/rho - case 'M': - Ctrl->M.active = true; - if(2 != sscanf(optarg, "%lf/%lf", &Ctrl->M.alpha, &Ctrl->M.beta)){ - GRTBadOptionError(command, M, ""); + // 模型参数， -P + case 'P': + Ctrl->P.active = true; + if(1 != sscanf(optarg, "%lf", &Ctrl->P.nu)){ + GRTBadOptionError(command, P, ""); } - if(Ctrl->M.alpha <= 0.0 || Ctrl->M.beta <= 0.0){ - GRTBadOptionError(command, M, "only support positive values."); - } - // 检查泊松比范围 - MYREAL nu = grt_possion_ratio(Ctrl->M.beta/Ctrl->M.alpha); - if(nu >= 0.5 || nu <= 0.0){ - GRTBadOptionError(command, M, "possion ratio (%lf) is out of bound.", nu); + if(Ctrl->P.nu <= 0.0 || Ctrl->P.nu >= 0.5){ + GRTBadOptionError(command, P, "possion ratio (%lf) is out of bound.", Ctrl->P.nu); } break; @@ -152,7 +146,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ // 检查必须设置的参数是否有设置 GRTCheckOptionSet(command, argc > 1); - GRTCheckOptionActive(command, Ctrl, M); + GRTCheckOptionActive(command, Ctrl, P); GRTCheckOptionActive(command, Ctrl, T); GRTCheckOptionActive(command, Ctrl, A); } @@ -170,7 +164,7 @@ int lamb1_main(int argc, char **argv){ getopt_from_command(Ctrl, argc, argv); // 求解，输出到标准输出 - grt_solve_lamb1(Ctrl->M.alpha, Ctrl->M.beta, Ctrl->T.ts, Ctrl->T.nt, Ctrl->A.azimuth, NULL); + grt_solve_lamb1(Ctrl->P.nu, Ctrl->T.ts, Ctrl->T.nt, Ctrl->A.azimuth, NULL); free_Ctrl(Ctrl); return EXIT_SUCCESS; diff --git a/pygrt/C_extension/src/lamb/lamb1.c b/pygrt/C_extension/src/lamb/lamb1.c index bfcc79d..2dff2f0 100644 --- a/pygrt/C_extension/src/lamb/lamb1.c +++ b/pygrt/C_extension/src/lamb/lamb1.c @@ -19,8 +19,6 @@ typedef struct { // 常量，不随时间变化 - MYREAL alpha; - MYREAL beta; MYREAL nu; MYREAL k; ///< k = beta/alpha; MYREAL kk; ///< kk = k*k; @@ -761,8 +759,13 @@ static void build_R(MYREAL tbar, VARS *V, MYREAL u[3][3]) void grt_solve_lamb1( - const MYREAL alpha, const MYREAL beta, const MYREAL *ts, const int nt, const MYREAL azimuth, MYREAL (*u)[3][3]) + const MYREAL nu, const MYREAL *ts, const int nt, const MYREAL azimuth, MYREAL (*u)[3][3]) { + // 检查泊松比范围 + if(nu <= 0.0 || nu >= 0.5){ + GRTRaiseError("possion ratio (%lf) is out of bound.", nu); + } + // 根据情况判断是打印在屏幕还是记录到内存中 bool isprint = (u == NULL); @@ -787,11 +790,9 @@ void grt_solve_lamb1( // 初始化相关变量 VARS V0 = {0}; VARS *V = &V0; - V0.alpha = alpha; - V0.beta = beta; - V0.k = beta/alpha; - V0.kk = V0.k*V0.k; - V0.nu = grt_possion_ratio(V0.k); + V0.kk = 0.5 * (1.0 - 2.0*nu)/(1.0 - nu); + V0.k = sqrt(V0.kk); + V0.nu = nu; V0.kpkp = 1.0 - V0.kk; V0.kp = sqrt(V0.kpkp); V0.h1 = 2.0 * V0.kk - 1.0; diff --git a/pygrt/C_extension/src/lamb/lamb_util.c b/pygrt/C_extension/src/lamb/lamb_util.c index 132ad68..0d32bd7 100644 --- a/pygrt/C_extension/src/lamb/lamb_util.c +++ b/pygrt/C_extension/src/lamb/lamb_util.c @@ -46,9 +46,4 @@ MYCOMPLEX grt_evalpoly2(const MYCOMPLEX *C, const int n, const MYCOMPLEX y, cons p *= y; } return res; -} - -MYREAL grt_possion_ratio(const MYREAL k) -{ - return 0.5 * (1.0 - 2.0*k*k)/(1.0 - k*k); } \ No newline at end of file diff --git a/pygrt/c_interfaces.py b/pygrt/c_interfaces.py index 24ba7ef..b06e141 100755 --- a/pygrt/c_interfaces.py +++ b/pygrt/c_interfaces.py @@ -159,5 +159,5 @@ def set_num_threads(n): """使用广义闭合解求解第一类 Lamb 问题""" C_grt_solve_lamb1.restype = None C_grt_solve_lamb1.argtypes = [ - REAL, REAL, PREAL, c_int, REAL, PREAL + REAL, PREAL, c_int, REAL, PREAL ] \ No newline at end of file diff --git a/pygrt/utils.py b/pygrt/utils.py index f7e120c..37c0cbb 100755 --- a/pygrt/utils.py +++ b/pygrt/utils.py @@ -1566,23 +1566,20 @@ def plot_statsdata_ptam(statsdata1:np.ndarray, statsdata2:np.ndarray, statsdata_ -def solve_lamb1(alpha:float, beta:float, ts:np.ndarray, azimuth:float): +def solve_lamb1(nu:float, ts:np.ndarray, azimuth:float): r""" 使用广义闭合解求解第一类 Lamb 问题，参考： 张海明, 冯禧 著. 2024. 地震学中的 Lamb 问题（下）. 科学出版社 - :param alpha: P 波速度 - :param beta: S 波速度 + :param nu: 泊松比， (0, 0.5) :param ts: 归一化时间序列 :math:`\bar{t}` ，:math:`\bar{t}=\dfrac{t}{r/\beta}` :param azimuth: 方位角，单位度 - :return: 形状为 (nt, 3, 3) 的归一化解，距离物理解还需乘 :math:`\pi^2 \mu r` + :return: 形状为 (nt, 3, 3) 的归一化解，距离物理解还需除 :math:`\pi^2 \mu r` """ # 检查数据 - if alpha <= 0.0 or beta <= 0.0: - raise ValueError("alpha and beta should be positive.") if np.any(ts < 0.0): raise ValueError("ts should be nonnegative.") if azimuth < 0.0 or azimuth > 360.0: @@ -1594,7 +1591,7 @@ def solve_lamb1(alpha:float, beta:float, ts:np.ndarray, azimuth:float): nt = len(ts) u = np.zeros((nt, 3, 3), dtype=NPCT_REAL_TYPE) - C_grt_solve_lamb1(alpha, beta, npct.as_ctypes(ts), nt, azimuth, npct.as_ctypes(u.ravel())) + C_grt_solve_lamb1(nu, npct.as_ctypes(ts), nt, azimuth, npct.as_ctypes(u.ravel())) return u