From afce860f4286ae0c0f3816356fdfb57259e61480 Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Tue, 4 Aug 2026 22:57:12 +0800 Subject: [PATCH] REFAC: replace sqrt(x*x+y*y) with hypot for 2D distance --- pygrt/C_extension/src/common/travt.c | 2 +- pygrt/C_extension/src/dynamic/grt_greenfn.c | 2 +- pygrt/C_extension/src/modal/grt_modsum.c | 2 +- pygrt/C_extension/src/modal/secular.c | 4 ++-- pygrt/C_extension/src/static/grt_static_greenfn.c | 2 +- pygrt/pymod.py | 4 ++-- 6 files changed, 8 insertions(+), 8 deletions(-) diff --git a/pygrt/C_extension/src/common/travt.c b/pygrt/C_extension/src/common/travt.c index c722411f..1595faa2 100644 --- a/pygrt/C_extension/src/common/travt.c +++ b/pygrt/C_extension/src/common/travt.c @@ -83,7 +83,7 @@ real_t grt_compute_travt1d( // ========================================================= // ------------------- 同层直达波 ---------------------- if(imax - imin == 1){ // 位于同一物理层 - travt = sqrt(dist*dist + depdif*depdif) / vsrc; + travt = hypot(dist, depdif) / vsrc; // printf("direct wave in same layer, travt=%f\n", travt); } else { diff --git a/pygrt/C_extension/src/dynamic/grt_greenfn.c b/pygrt/C_extension/src/dynamic/grt_greenfn.c index 7cd95166..645b2c83 100644 --- a/pygrt/C_extension/src/dynamic/grt_greenfn.c +++ b/pygrt/C_extension/src/dynamic/grt_greenfn.c @@ -1046,7 +1046,7 @@ int greenfn_main(int argc, char **argv) { real_t delayT = 0.0; if (! Ctrl->E.refFirstP){ delayT = Ctrl->E.delayT0; - if(Ctrl->E.delayV0 > 0.0) delayT += sqrt( GRT_SQUARE(dist) + GRT_SQUARE(Ctrl->D.deprcv - Ctrl->D.depsrc) ) / Ctrl->E.delayV0; + if(Ctrl->E.delayV0 > 0.0) delayT += hypot( dist, Ctrl->D.deprcv - Ctrl->D.depsrc ) / Ctrl->E.delayV0; } else { delayT = Ctrl->E.delayT0 + travtPS[ir][0]; } diff --git a/pygrt/C_extension/src/modal/grt_modsum.c b/pygrt/C_extension/src/modal/grt_modsum.c index 510c0b93..280342bd 100644 --- a/pygrt/C_extension/src/modal/grt_modsum.c +++ b/pygrt/C_extension/src/modal/grt_modsum.c @@ -623,7 +623,7 @@ int modsum_main(int argc, char **argv){ real_t delayT = 0.0; if (! Ctrl->E.refFirstP){ delayT = Ctrl->E.delayT0; - if(Ctrl->E.delayV0 > 0.0) delayT += sqrt( GRT_SQUARE(dist) + GRT_SQUARE(Ctrl->D.deprcv - Ctrl->D.depsrc) ) / Ctrl->E.delayV0; + if(Ctrl->E.delayV0 > 0.0) delayT += hypot( dist, Ctrl->D.deprcv - Ctrl->D.depsrc ) / Ctrl->E.delayV0; } else { delayT = Ctrl->E.delayT0 + sac->hd.t0; } diff --git a/pygrt/C_extension/src/modal/secular.c b/pygrt/C_extension/src/modal/secular.c index b5fb32f1..e941bca6 100644 --- a/pygrt/C_extension/src/modal/secular.c +++ b/pygrt/C_extension/src/modal/secular.c @@ -379,8 +379,8 @@ void grt_secular_function_potential_Rayl( // 返回对应的垂直波函数 if(ppot != NULL){ // 假设一个比例 - ppot[2] = - Det[0][1] / sqrt( GRT_SQUARE(fabs(Det[0][0])) + GRT_SQUARE(fabs(Det[0][1])) ); - ppot[3] = + Det[0][0] / sqrt( GRT_SQUARE(fabs(Det[0][0])) + GRT_SQUARE(fabs(Det[0][1])) ); + ppot[2] = - Det[0][1] / hypot( fabs(Det[0][0]), fabs(Det[0][1]) ); + ppot[3] = + Det[0][0] / hypot( fabs(Det[0][0]), fabs(Det[0][1]) ); grt_cmat2x1_mul(mstat->M_BL.RD, ppot+2, ppot); } } diff --git a/pygrt/C_extension/src/static/grt_static_greenfn.c b/pygrt/C_extension/src/static/grt_static_greenfn.c index 22c7a08c..6dc7c73b 100644 --- a/pygrt/C_extension/src/static/grt_static_greenfn.c +++ b/pygrt/C_extension/src/static/grt_static_greenfn.c @@ -555,7 +555,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ Ctrl->rs = (real_t*)calloc(Ctrl->nr, sizeof(real_t)); for(size_t ix=0; ixX.nx; ++ix){ for(size_t iy=0; iyY.ny; ++iy){ - Ctrl->rs[iy + ix*Ctrl->Y.ny] = GRT_MAX(sqrt(GRT_SQUARE(Ctrl->X.xs[ix]) + GRT_SQUARE(Ctrl->Y.ys[iy])), GRT_MIN_DISTANCE); // 避免0震中距 + Ctrl->rs[iy + ix*Ctrl->Y.ny] = GRT_MAX(hypot(Ctrl->X.xs[ix], Ctrl->Y.ys[iy]), GRT_MIN_DISTANCE); // 避免0震中距 } } diff --git a/pygrt/pymod.py b/pygrt/pymod.py index b1046beb..accb5f29 100755 --- a/pygrt/pymod.py +++ b/pygrt/pymod.py @@ -406,7 +406,7 @@ def _get_stream_from_grn_spectra( # 计算延迟 delayT = delayT0 if delayV0 > 0.0: - delayT += np.sqrt(dist**2 + (deprcv-depsrc)**2)/delayV0 + delayT += np.hypot(dist, deprcv-depsrc)/delayV0 # 计算走时 travtP, travtS = self.compute_travt1d(dist) @@ -638,7 +638,7 @@ def compute_static_grn( rs = np.zeros((nr,), dtype=NPCT_REAL_TYPE) for iy in range(ny): for ix in range(nx): - rs[ix + iy*nx] = max(np.sqrt(xarr[ix]**2 + yarr[iy]**2), 1e-5) + rs[ix + iy*nx] = max(np.hypot(xarr[ix], yarr[iy]), 1e-5) c_rs = npct.as_ctypes(rs) # 设置波数积分间隔