Skip to content

FEAT: update lamb1 interfaces #118

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
Nov 13, 2025
Merged

FEAT: update lamb1 interfaces #118

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
5 changes: 2 additions & 3 deletions pygrt/C_extension/include/grt/lamb/lamb1.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,13 +17,12 @@
/**
* 使用广义闭合解求解第一类 Lamb 问题
*
* @param[in] alpha P 波速度
* @param[in] beta S 波速度
* @param[in] nu 泊松比， (0, 0.5)
* @param[in] ts 归一化时间序列
* @param[in] nt 时间序列点数
* @param[in] azimuth 方位角，单位度
* @param[out] u 记录结果的指针，如果为NULL则输出到标准输出
*
*/
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]);
12 changes: 1 addition & 11 deletions pygrt/C_extension/include/grt/lamb/lamb_util.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
MYCOMPLEX grt_evalpoly2(const MYCOMPLEX *C, const int n, const MYCOMPLEX y, const int offset);
36 changes: 15 additions & 21 deletions pygrt/C_extension/src/lamb/grt_lamb1.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,9 +21,8 @@ typedef struct {
/** 模型参数 */
struct {
bool active;
MYREAL alpha; ///< P 波速度
MYREAL beta; ///< S 波速度
} M;
MYREAL nu; ///< 泊松比
} P;

/** 归一化时间序列 */
struct {
Expand Down Expand Up @@ -62,11 +61,11 @@ printf("\n"
"\n\n"
"Usage:\n"
"----------------------------------------------------------------\n"
" grt lamb1 -M<Vp>/<Vs> -T<t1>/<t2>/<dt> -A<azimuth>\n"
" grt lamb1 -P<nu> -T<t1>/<t2>/<dt> -A<azimuth>\n"
"\n\n"
"Options:\n"
"----------------------------------------------------------------\n"
" -M<Vp>/<Vs> Vp and Vs of the halfspace.\n"
" -P<nu> Possion ratio of the halfspace, (0, 0.5).\n"
"\n"
" -T<t1>/<t2>/<dt>\n"
" Dimensionless time.\n"
Expand All @@ -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"
);
}
Expand All @@ -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<nu>
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;

Expand Down Expand Up @@ -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);
}
Expand All @@ -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;
Expand Down
17 changes: 9 additions & 8 deletions pygrt/C_extension/src/lamb/lamb1.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,8 +19,6 @@

typedef struct {
// 常量，不随时间变化
MYREAL alpha;
MYREAL beta;
MYREAL nu;
MYREAL k; ///< k = beta/alpha;
MYREAL kk; ///< kk = k*k;
Expand Down Expand Up @@ -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);

Expand All @@ -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;
Expand Down
5 changes: 0 additions & 5 deletions pygrt/C_extension/src/lamb/lamb_util.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}
2 changes: 1 addition & 1 deletion pygrt/c_interfaces.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
]
11 changes: 4 additions & 7 deletions pygrt/utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand All @@ -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