Skip to content

Commit cab9d59

Browse files
Dengda98Dengda98
authored
FEAT: update lamb1 interfaces (#118)
1 parent 7217d3d commit cab9d59

File tree

7 files changed

+32
-56
lines changed

7 files changed

+32
-56
lines changed

‎pygrt/C_extension/include/grt/lamb/lamb1.h‎

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -17,13 +17,12 @@
1717
/**
1818
* 使用广义闭合解求解第一类 Lamb 问题
1919
*
20-
* @param[in] alpha P 波速度
21-
* @param[in] beta S 波速度
20+
* @param[in] nu 泊松比， (0, 0.5)
2221
* @param[in] ts 归一化时间序列
2322
* @param[in] nt 时间序列点数
2423
* @param[in] azimuth 方位角，单位度
2524
* @param[out] u 记录结果的指针，如果为NULL则输出到标准输出
2625
*
2726
*/
2827
void grt_solve_lamb1(
29-
const MYREAL alpha, const MYREAL beta, const MYREAL *ts, const int nt, const MYREAL azimuth, MYREAL (*u)[3][3]);
28+
const MYREAL nu, const MYREAL *ts, const int nt, const MYREAL azimuth, MYREAL (*u)[3][3]);

‎pygrt/C_extension/include/grt/lamb/lamb_util.h‎

Lines changed: 1 addition & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -31,14 +31,4 @@ void grt_roots3(const MYREAL a, const MYREAL b, const MYREAL c, MYCOMPLEX y3[3])
3131
* @return 多项式结果
3232
*
3333
*/
34-
MYCOMPLEX grt_evalpoly2(const MYCOMPLEX *C, const int n, const MYCOMPLEX y, const int offset);
35-
36-
37-
/**
38-
* 根据 Vs/Vp 比值，计算泊松比， \f$ \nu = \dfrac{1}{2} \dfrac{1 - 2k^2}{1 - k^2}, k=\dfrac{V_S}{V_P} \f$
39-
*
40-
* @param[in] k Vs/Vp
41-
* @return 泊松比
42-
*/
43-
44-
MYREAL grt_possion_ratio(const MYREAL k);
34+
MYCOMPLEX grt_evalpoly2(const MYCOMPLEX *C, const int n, const MYCOMPLEX y, const int offset);

‎pygrt/C_extension/src/lamb/grt_lamb1.c‎

Lines changed: 15 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -21,9 +21,8 @@ typedef struct {
2121
/** 模型参数 */
2222
struct {
2323
bool active;
24-
MYREAL alpha; ///< P 波速度
25-
MYREAL beta; ///< S 波速度
26-
} M;
24+
MYREAL nu; ///< 泊松比
25+
} P;
2726

2827
/** 归一化时间序列 */
2928
struct {
@@ -62,11 +61,11 @@ printf("\n"
6261
"\n\n"
6362
"Usage:\n"
6463
"----------------------------------------------------------------\n"
65-
" grt lamb1 -M<Vp>/<Vs> -T<t1>/<t2>/<dt> -A<azimuth>\n"
64+
" grt lamb1 -P<nu> -T<t1>/<t2>/<dt> -A<azimuth>\n"
6665
"\n\n"
6766
"Options:\n"
6867
"----------------------------------------------------------------\n"
69-
" -M<Vp>/<Vs> Vp and Vs of the halfspace.\n"
68+
" -P<nu> Possion ratio of the halfspace, (0, 0.5).\n"
7069
"\n"
7170
" -T<t1>/<t2>/<dt>\n"
7271
" Dimensionless time.\n"
@@ -80,7 +79,7 @@ printf("\n"
8079
"\n\n"
8180
"Examples:\n"
8281
"----------------------------------------------------------------\n"
83-
" grt lamb1 -M8.0/4.62 -T0/2/1e-3 -A30\n"
82+
" grt lamb1 -P0.25 -T0/2/1e-3 -A30\n"
8483
"\n\n\n"
8584
);
8685
}
@@ -91,21 +90,16 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){
9190
char* command = Ctrl->name;
9291
int opt;
9392

94-
while ((opt = getopt(argc, argv, ":M:T:A:h")) != -1) {
93+
while ((opt = getopt(argc, argv, ":P:T:A:h")) != -1) {
9594
switch (opt) {
96-
// 模型参数， -Malpha/beta/rho
97-
case 'M':
98-
Ctrl->M.active = true;
99-
if(2 != sscanf(optarg, "%lf/%lf", &Ctrl->M.alpha, &Ctrl->M.beta)){
100-
GRTBadOptionError(command, M, "");
95+
// 模型参数， -P<nu>
96+
case 'P':
97+
Ctrl->P.active = true;
98+
if(1 != sscanf(optarg, "%lf", &Ctrl->P.nu)){
99+
GRTBadOptionError(command, P, "");
101100
}
102-
if(Ctrl->M.alpha <= 0.0 || Ctrl->M.beta <= 0.0){
103-
GRTBadOptionError(command, M, "only support positive values.");
104-
}
105-
// 检查泊松比范围
106-
MYREAL nu = grt_possion_ratio(Ctrl->M.beta/Ctrl->M.alpha);
107-
if(nu >= 0.5 || nu <= 0.0){
108-
GRTBadOptionError(command, M, "possion ratio (%lf) is out of bound.", nu);
101+
if(Ctrl->P.nu <= 0.0 || Ctrl->P.nu >= 0.5){
102+
GRTBadOptionError(command, P, "possion ratio (%lf) is out of bound.", Ctrl->P.nu);
109103
}
110104
break;
111105

@@ -152,7 +146,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){
152146

153147
// 检查必须设置的参数是否有设置
154148
GRTCheckOptionSet(command, argc > 1);
155-
GRTCheckOptionActive(command, Ctrl, M);
149+
GRTCheckOptionActive(command, Ctrl, P);
156150
GRTCheckOptionActive(command, Ctrl, T);
157151
GRTCheckOptionActive(command, Ctrl, A);
158152
}
@@ -170,7 +164,7 @@ int lamb1_main(int argc, char **argv){
170164
getopt_from_command(Ctrl, argc, argv);
171165

172166
// 求解，输出到标准输出
173-
grt_solve_lamb1(Ctrl->M.alpha, Ctrl->M.beta, Ctrl->T.ts, Ctrl->T.nt, Ctrl->A.azimuth, NULL);
167+
grt_solve_lamb1(Ctrl->P.nu, Ctrl->T.ts, Ctrl->T.nt, Ctrl->A.azimuth, NULL);
174168

175169
free_Ctrl(Ctrl);
176170
return EXIT_SUCCESS;

‎pygrt/C_extension/src/lamb/lamb1.c‎

Lines changed: 9 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -19,8 +19,6 @@
1919

2020
typedef struct {
2121
// 常量，不随时间变化
22-
MYREAL alpha;
23-
MYREAL beta;
2422
MYREAL nu;
2523
MYREAL k; ///< k = beta/alpha;
2624
MYREAL kk; ///< kk = k*k;
@@ -761,8 +759,13 @@ static void build_R(MYREAL tbar, VARS *V, MYREAL u[3][3])
761759

762760

763761
void grt_solve_lamb1(
764-
const MYREAL alpha, const MYREAL beta, const MYREAL *ts, const int nt, const MYREAL azimuth, MYREAL (*u)[3][3])
762+
const MYREAL nu, const MYREAL *ts, const int nt, const MYREAL azimuth, MYREAL (*u)[3][3])
765763
{
764+
// 检查泊松比范围
765+
if(nu <= 0.0 || nu >= 0.5){
766+
GRTRaiseError("possion ratio (%lf) is out of bound.", nu);
767+
}
768+
766769
// 根据情况判断是打印在屏幕还是记录到内存中
767770
bool isprint = (u == NULL);
768771

@@ -787,11 +790,9 @@ void grt_solve_lamb1(
787790
// 初始化相关变量
788791
VARS V0 = {0};
789792
VARS *V = &V0;
790-
V0.alpha = alpha;
791-
V0.beta = beta;
792-
V0.k = beta/alpha;
793-
V0.kk = V0.k*V0.k;
794-
V0.nu = grt_possion_ratio(V0.k);
793+
V0.kk = 0.5 * (1.0 - 2.0*nu)/(1.0 - nu);
794+
V0.k = sqrt(V0.kk);
795+
V0.nu = nu;
795796
V0.kpkp = 1.0 - V0.kk;
796797
V0.kp = sqrt(V0.kpkp);
797798
V0.h1 = 2.0 * V0.kk - 1.0;

‎pygrt/C_extension/src/lamb/lamb_util.c‎

Lines changed: 0 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -46,9 +46,4 @@ MYCOMPLEX grt_evalpoly2(const MYCOMPLEX *C, const int n, const MYCOMPLEX y, cons
4646
p *= y;
4747
}
4848
return res;
49-
}
50-
51-
MYREAL grt_possion_ratio(const MYREAL k)
52-
{
53-
return 0.5 * (1.0 - 2.0*k*k)/(1.0 - k*k);
5449
}

‎pygrt/c_interfaces.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -159,5 +159,5 @@ def set_num_threads(n):
159159
"""使用广义闭合解求解第一类 Lamb 问题"""
160160
C_grt_solve_lamb1.restype = None
161161
C_grt_solve_lamb1.argtypes = [
162-
REAL, REAL, PREAL, c_int, REAL, PREAL
162+
REAL, PREAL, c_int, REAL, PREAL
163163
]

‎pygrt/utils.py‎

Lines changed: 4 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -1566,23 +1566,20 @@ def plot_statsdata_ptam(statsdata1:np.ndarray, statsdata2:np.ndarray, statsdata_
15661566

15671567

15681568

1569-
def solve_lamb1(alpha:float, beta:float, ts:np.ndarray, azimuth:float):
1569+
def solve_lamb1(nu:float, ts:np.ndarray, azimuth:float):
15701570
r"""
15711571
使用广义闭合解求解第一类 Lamb 问题，参考：
15721572
15731573
张海明, 冯禧 著. 2024. 地震学中的 Lamb 问题（下）. 科学出版社
15741574
1575-
:param alpha: P 波速度
1576-
:param beta: S 波速度
1575+
:param nu: 泊松比， (0, 0.5)
15771576
:param ts: 归一化时间序列 :math:`\bar{t}` ，:math:`\bar{t}=\dfrac{t}{r/\beta}`
15781577
:param azimuth: 方位角，单位度
15791578
1580-
:return: 形状为 (nt, 3, 3) 的归一化解，距离物理解还需乘 :math:`\pi^2 \mu r`
1579+
:return: 形状为 (nt, 3, 3) 的归一化解，距离物理解还需除 :math:`\pi^2 \mu r`
15811580
"""
15821581

15831582
# 检查数据
1584-
if alpha <= 0.0 or beta <= 0.0:
1585-
raise ValueError("alpha and beta should be positive.")
15861583
if np.any(ts < 0.0):
15871584
raise ValueError("ts should be nonnegative.")
15881585
if azimuth < 0.0 or azimuth > 360.0:
@@ -1594,7 +1591,7 @@ def solve_lamb1(alpha:float, beta:float, ts:np.ndarray, azimuth:float):
15941591
nt = len(ts)
15951592
u = np.zeros((nt, 3, 3), dtype=NPCT_REAL_TYPE)
15961593

1597-
C_grt_solve_lamb1(alpha, beta, npct.as_ctypes(ts), nt, azimuth, npct.as_ctypes(u.ravel()))
1594+
C_grt_solve_lamb1(nu, npct.as_ctypes(ts), nt, azimuth, npct.as_ctypes(u.ravel()))
15981595

15991596
return u
16001597

0 commit comments

Comments
 (0)