From f0bb97ac6796414984ba3bde86eba33f997157f2 Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Wed, 3 Sep 2025 12:14:29 +0800 Subject: [PATCH 1/5] BUILD: update dependency --- pygrt/C_extension/Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygrt/C_extension/Makefile b/pygrt/C_extension/Makefile index 1c92f3cb..99ceea79 100644 --- a/pygrt/C_extension/Makefile +++ b/pygrt/C_extension/Makefile @@ -89,7 +89,7 @@ endif $(BUILD_DIR)/%.d: $(SRC_DIR)/%.c @mkdir -p $(shell dirname $@) @$(CC) $(CFLAGS) -MM $< > $@.$$$$; \ - sed 's,\($*\)\.o[ :]*,$(BUILD_DIR)/\1.o $@ : ,g' < $@.$$$$ > $@; \ + sed 's,\($(notdir $*)\)\.o[ :]*,$(shell dirname $@)/\1.o $@ : ,g' < $@.$$$$ > $@; \ rm -f $@.$$$$ # ---------------------------------------------------------------------- From 84d349e010023c2e00f12b071778f14eb9f97ec1 Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Wed, 3 Sep 2025 13:13:16 +0800 Subject: [PATCH 2/5] REFAC: rename macros and global variables with prefix `GRT` --- pygrt/C_extension/include/grt/common/const.h | 79 ++++---- pygrt/C_extension/include/grt/common/dwm.h | 6 +- pygrt/C_extension/include/grt/common/fim.h | 6 +- .../C_extension/include/grt/common/integral.h | 14 +- .../C_extension/include/grt/common/iostats.h | 6 +- pygrt/C_extension/include/grt/common/kernel.h | 4 +- pygrt/C_extension/include/grt/common/matrix.h | 12 +- pygrt/C_extension/include/grt/common/ptam.h | 8 +- .../include/grt/common/radiation.h | 4 +- .../include/grt/common/recursion.h | 13 +- pygrt/C_extension/include/grt/common/safim.h | 6 +- pygrt/C_extension/include/grt/common/util.h | 6 +- pygrt/C_extension/include/grt/dynamic/grn.h | 6 +- .../include/grt/dynamic/propagate.h | 4 +- .../C_extension/include/grt/dynamic/source.h | 4 +- .../include/grt/static/static_grn.h | 6 +- .../include/grt/static/static_propagate.h | 4 +- .../include/grt/static/static_source.h | 4 +- pygrt/C_extension/src/common/attenuation.c | 4 +- pygrt/C_extension/src/common/bessel.c | 4 +- pygrt/C_extension/src/common/const.c | 13 +- pygrt/C_extension/src/common/dwm.c | 38 ++-- pygrt/C_extension/src/common/fim.c | 82 ++++---- pygrt/C_extension/src/common/integral.c | 68 +++---- pygrt/C_extension/src/common/iostats.c | 56 +++--- pygrt/C_extension/src/common/model.c | 12 +- pygrt/C_extension/src/common/ptam.c | 137 +++++++------ pygrt/C_extension/src/common/quadratic.c | 2 +- pygrt/C_extension/src/common/radiation.c | 4 +- pygrt/C_extension/src/common/recursion.c | 56 +++--- pygrt/C_extension/src/common/safim.c | 114 +++++------ pygrt/C_extension/src/common/util.c | 30 +-- pygrt/C_extension/src/dynamic/grn.c | 68 +++---- pygrt/C_extension/src/dynamic/grt_greenfn.c | 22 +-- pygrt/C_extension/src/dynamic/grt_rotation.c | 2 +- pygrt/C_extension/src/dynamic/grt_strain.c | 2 +- pygrt/C_extension/src/dynamic/grt_stress.c | 2 +- pygrt/C_extension/src/dynamic/grt_syn.c | 20 +- pygrt/C_extension/src/dynamic/layer.c | 182 +++++++++--------- pygrt/C_extension/src/dynamic/propagate.c | 142 +++++++------- pygrt/C_extension/src/dynamic/source.c | 26 +-- .../src/static/grt_static_greenfn.c | 34 ++-- .../src/static/grt_static_rotation.c | 2 +- .../src/static/grt_static_strain.c | 2 +- .../src/static/grt_static_stress.c | 2 +- pygrt/C_extension/src/static/grt_static_syn.c | 40 ++-- pygrt/C_extension/src/static/static_grn.c | 54 +++--- pygrt/C_extension/src/static/static_layer.c | 38 ++-- .../C_extension/src/static/static_propagate.c | 94 ++++----- pygrt/C_extension/src/static/static_source.c | 38 ++-- 50 files changed, 797 insertions(+), 785 deletions(-) diff --git a/pygrt/C_extension/include/grt/common/const.h b/pygrt/C_extension/include/grt/common/const.h index 8d429668..6233c769 100755 --- a/pygrt/C_extension/include/grt/common/const.h +++ b/pygrt/C_extension/include/grt/common/const.h @@ -10,6 +10,7 @@ #include #include +#include #include "grt/common/checkerror.h" @@ -46,13 +47,6 @@ typedef int MYINT; ///< 整数 #endif // 常数 -#define RZERO 0.0 ///< 0.0 -#define RQUART 0.25 ///< 0.25 -#define RHALF 0.5 ///< 0.5 -#define RONE 1.0 ///< 1.0 -#define RTWO 2.0 ///< 2.0 -#define RTHREE 3.0 ///< 3.0 -#define RFOUR 4.0 ///< 4.0 #define RTWOTHIRD 0.6666666666666667 ///< 2/3 #define PI 3.141592653589793 ///< \f$ \pi \f$ #define PI2 6.283185307179586 ///< \f$ 2\pi \f$ @@ -63,21 +57,25 @@ typedef int MYINT; ///< 整数 #define SEVENQUARTERPI 5.497787143782138 ///< \f$ \frac{7\pi}{4} \f$ #define INV_SQRT_TWO 0.7071067811865475 ///< \f$ \frac{1}{\sqrt{2}} \f$ #define DEG1 0.017453292519943295 ///< \f$ \frac{\pi}{180} \f$ +#define GOLDEN_RATIO 0.6180339887498949 ///< \f$ \frac{\sqrt{5}-1}{2} \f$ -#define CZERO CMPLX(RZERO, RZERO) ///< 0.0 + j0.0 -#define CONE CMPLX(RONE, RZERO) ///< 1.0 + j0.0 -#define INIT_C_ZERO_2x2_MATRIX {{CZERO, CZERO}, {CZERO, CZERO}} ///< 初始化复数0矩阵 -#define INIT_C_IDENTITY_2x2_MATRIX {{CONE, CZERO}, {CZERO, CONE}} ///< 初始化复数单位阵 +#define GRT_INIT_ZERO_2x2_MATRIX {{0, 0}, {0, 0}} ///< 初始化复数0矩阵 +#define GRT_INIT_IDENTITY_2x2_MATRIX {{1, 0}, {0, 1}} ///< 初始化复数单位阵 -#define MIN_DEPTH_GAP_SRC_RCV 1.0 ///< 震源和台站的最小深度差(不做绝对限制,仅用于参考波数积分上限) +#define GRT_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" ///< 复数输出格式 +#define GRT_CMPLX_FMT "%18.8e%16.8e" ///< 复数输出格式 #define GRT_STR_CMPLX_FMT "%34s" ///< 与复数格式同长度的字符串输出格式 +#define GRT_STRING_LONG_FMT "%25s" ///< 字符串输出格式(更长) +#define GRT_REAL_LONG_FMT "%25.16e" ///< 浮点数输出格式(更长) +#define GRT_CMPLX_LONG_FMT "%25.16e %25.16e" ///< 复数输出格式(更长) +#define GRT_STR_CMPLX_LONG_FMT "%51s" ///< 与复数格式同长度的字符串输出格式(更长) +#define GRT_CMPLX_SPLIT(x) creal(x), cimag(x) ///< 用于打印复数时将实部虚部分开 #define GRT_MAX(a, b) ((a) > (b) ? (a) : (b)) ///< 求两者较大值 #define GRT_MIN(a, b) ((a) < (b) ? (a) : (b)) ///< 求两者较小值 @@ -89,6 +87,7 @@ typedef int MYINT; ///< 整数 #define GRT_SAFE_FREE_PTR(ptr) ({\ if(ptr!=NULL) {\ free(ptr);\ + ptr=NULL;\ }\ }) @@ -96,9 +95,11 @@ typedef int MYINT; ///< 整数 #define GRT_SAFE_FREE_PTR_ARRAY(ptr, count) ({\ if(ptr!=NULL){\ for(MYINT i=0; i /** 分别对应爆炸源(0阶),垂直力源(0阶),水平力源(1阶),剪切源(0,1,2阶) */ -const MYINT SRC_M_ORDERS[SRC_M_NUM] = {0, 0, 1, 0, 1, 2}; +const MYINT GRT_SRC_M_ORDERS[GRT_SRC_M_NUM] = {0, 0, 1, 0, 1, 2}; /** 不同震源类型使用的格林函数类型,0为Gij,1为格林函数导数Gij,k */ -const MYINT SRC_M_GTYPES[SRC_M_NUM] = {1, 0, 0, 1, 1, 1}; +const MYINT GRT_SRC_M_GTYPES[GRT_SRC_M_NUM] = {1, 0, 0, 1, 1, 1}; /** 不同震源,不同阶数的名称简写,用于命名 */ -const char *SRC_M_NAME_ABBR[SRC_M_NUM] = {"EX", "VF", "HF", "DD", "DS", "SS"}; +const char *GRT_SRC_M_NAME_ABBR[GRT_SRC_M_NUM] = {"EX", "VF", "HF", "DD", "DS", "SS"}; /** q, w, v 名称代号 */ -const char qwvchs[] = {'q', 'w', 'v'}; +const char GRT_QWV_CODES[] = {'q', 'w', 'v'}; /** ZRT三分量代号 */ -const char ZRTchs[] = {'Z', 'R', 'T'}; +const char GRT_ZRT_CODES[] = {'Z', 'R', 'T'}; /** ZNE三分量代号 */ -const char ZNEchs[] = {'Z', 'N', 'E'}; +const char GRT_ZNE_CODES[] = {'Z', 'N', 'E'}; diff --git a/pygrt/C_extension/src/common/dwm.c b/pygrt/C_extension/src/common/dwm.c index d151d552..c0445409 100644 --- a/pygrt/C_extension/src/common/dwm.c +++ b/pygrt/C_extension/src/common/dwm.c @@ -26,17 +26,17 @@ MYREAL grt_discrete_integ( const GRT_MODEL1D *mod1d, MYREAL dk, MYREAL kmax, MYREAL keps, MYCOMPLEX omega, MYINT nr, MYREAL *rs, - MYCOMPLEX sum_J[nr][SRC_M_NUM][INTEG_NUM], + MYCOMPLEX sum_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], bool calc_upar, - MYCOMPLEX sum_uiz_J[nr][SRC_M_NUM][INTEG_NUM], - MYCOMPLEX sum_uir_J[nr][SRC_M_NUM][INTEG_NUM], + MYCOMPLEX sum_uiz_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], + MYCOMPLEX sum_uir_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], FILE *fstats, GRT_KernelFunc kerfunc, MYINT *stats) { - MYCOMPLEX SUM[SRC_M_NUM][INTEG_NUM]; + MYCOMPLEX SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM]; // 不同震源不同阶数的核函数 F(k, w) - MYCOMPLEX QWV[SRC_M_NUM][QWV_NUM]; - MYCOMPLEX QWV_uiz[SRC_M_NUM][QWV_NUM]; + MYCOMPLEX QWV[GRT_SRC_M_NUM][GRT_QWV_NUM]; + MYCOMPLEX QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]; MYREAL k = 0.0; MYINT ik = 0; @@ -59,7 +59,7 @@ MYREAL grt_discrete_integ( // printf("w=%15.5e, ik=%d\n", creal(omega), ik); // 计算核函数 F(k, w) kerfunc(mod1d, omega, k, QWV, calc_upar, QWV_uiz, stats); - if(*stats==INVERSE_FAILURE) goto BEFORE_RETURN; + if(*stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; // 记录积分核函数 if(fstats!=NULL) grt_write_stats(fstats, k, QWV); @@ -69,9 +69,9 @@ MYREAL grt_discrete_integ( for(MYINT ir=0; ir RZERO){ + if(keps > 0.0){ iendkrs[ir] = iendk0; iendk = iendk && iendkrs[ir]; } else { @@ -107,8 +107,8 @@ MYREAL grt_discrete_integ( grt_int_Pk(k, rs[ir], QWV_uiz, false, SUM); // keps不参与计算位移空间导数的积分,背后逻辑认为u收敛,则uiz也收敛 - for(MYINT i=0; iQbinv[i] > 0.0)? grt_attenuation_law(mod1d->Qbinv[i], omega) : 1.0; ka0 = omega/(Va0*atna); - kb0 = (Vb0>RZERO)? omega/(Vb0*atnb) : CZERO; + kb0 = (Vb0>0.0)? omega/(Vb0*atnb) : 0.0; mod1d->kaka[i] = ka0*ka0; mod1d->kbkb[i] = kb0*kb0; @@ -376,15 +376,15 @@ GRT_MODEL1D * grt_read_mod1d_from_file(const char *command, const char *modelpat } -void grt_get_mod1d_vmin_vmax(const GRT_MODEL1D *mod1d, double *vmin, double *vmax){ - *vmin = __DBL_MAX__; - *vmax = RZERO; +void grt_get_mod1d_vmin_vmax(const GRT_MODEL1D *mod1d, MYREAL *vmin, MYREAL *vmax){ + *vmin = 9.0e30; + *vmax = 0.0; const MYREAL *Va = mod1d->Va; const MYREAL *Vb = mod1d->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]; + if(Vb[i] < *vmin && Vb[i] > 0.0) *vmin = Vb[i]; + if(Vb[i] > *vmax && Vb[i] > 0.0) *vmax = Vb[i]; } } \ No newline at end of file diff --git a/pygrt/C_extension/src/common/ptam.c b/pygrt/C_extension/src/common/ptam.c index e7cb1d5e..069c4a8c 100755 --- a/pygrt/C_extension/src/common/ptam.c +++ b/pygrt/C_extension/src/common/ptam.c @@ -43,21 +43,21 @@ */ static void process_peak_or_trough( MYINT ir, MYINT im, MYINT v, MYREAL k, MYREAL dk, - MYCOMPLEX (*J3)[PTAM_WINDOW_SIZE][SRC_M_NUM][INTEG_NUM], MYREAL (*Kpt)[SRC_M_NUM][INTEG_NUM][PTAM_MAX_PT], - MYCOMPLEX (*Fpt)[SRC_M_NUM][INTEG_NUM][PTAM_MAX_PT], MYINT (*Ipt)[SRC_M_NUM][INTEG_NUM], MYINT (*Gpt)[SRC_M_NUM][INTEG_NUM], bool *iendk0) + MYCOMPLEX (*J3)[GRT_PTAM_WINDOW_SIZE][GRT_SRC_M_NUM][GRT_INTEG_NUM], MYREAL (*Kpt)[GRT_SRC_M_NUM][GRT_INTEG_NUM][GRT_PTAM_PT_MAX], + MYCOMPLEX (*Fpt)[GRT_SRC_M_NUM][GRT_INTEG_NUM][GRT_PTAM_PT_MAX], MYINT (*Ipt)[GRT_SRC_M_NUM][GRT_INTEG_NUM], MYINT (*Gpt)[GRT_SRC_M_NUM][GRT_INTEG_NUM], bool *iendk0) { MYCOMPLEX tmp0; - if (Gpt[ir][im][v] >= PTAM_WINDOW_SIZE-1 && Ipt[ir][im][v] < PTAM_MAX_PT) { + if (Gpt[ir][im][v] >= GRT_PTAM_WINDOW_SIZE-1 && Ipt[ir][im][v] < GRT_PTAM_PT_MAX) { if (grt_cplx_peak_or_trough(im, v, J3[ir], k, dk, &Kpt[ir][im][v][Ipt[ir][im][v]], &tmp0) != 0) { Fpt[ir][im][v][Ipt[ir][im][v]++] = tmp0; Gpt[ir][im][v] = 0; - } else if (Gpt[ir][im][v] >= PTAM_MAX_WAITS) { // 不再等待,直接取中点作为波峰波谷 + } else if (Gpt[ir][im][v] >= GRT_PTAM_WAITS_MAX) { // 不再等待,直接取中点作为波峰波谷 Kpt[ir][im][v][Ipt[ir][im][v]] = k - dk; Fpt[ir][im][v][Ipt[ir][im][v]++] = J3[ir][1][im][v]; Gpt[ir][im][v] = 0; } } - *iendk0 = *iendk0 && (Ipt[ir][im][v] == PTAM_MAX_PT); + *iendk0 = *iendk0 && (Ipt[ir][im][v] == GRT_PTAM_PT_MAX); } @@ -82,30 +82,30 @@ static void process_peak_or_trough( */ static void ptam_once( const MYINT ir, const MYINT nr, const MYREAL precoef, MYREAL k, MYREAL dk, - MYCOMPLEX SUM3[nr][PTAM_WINDOW_SIZE][SRC_M_NUM][INTEG_NUM], - MYCOMPLEX sum_J[nr][SRC_M_NUM][INTEG_NUM], - MYREAL Kpt[nr][SRC_M_NUM][INTEG_NUM][PTAM_MAX_PT], - MYCOMPLEX Fpt[nr][SRC_M_NUM][INTEG_NUM][PTAM_MAX_PT], - MYINT Ipt[nr][SRC_M_NUM][INTEG_NUM], - MYINT Gpt[nr][SRC_M_NUM][INTEG_NUM], + MYCOMPLEX SUM3[nr][GRT_PTAM_WINDOW_SIZE][GRT_SRC_M_NUM][GRT_INTEG_NUM], + MYCOMPLEX sum_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], + MYREAL Kpt[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM][GRT_PTAM_PT_MAX], + MYCOMPLEX Fpt[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM][GRT_PTAM_PT_MAX], + MYINT Ipt[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], + MYINT Gpt[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], bool *iendk0) { *iendk0 = true; - for(MYINT i=0; ik3[2] - item_pt->k3[0]; - const MYCOMPLEX (*F3)[SRC_M_NUM][QWV_NUM] = (isuiz)? item_pt->F3_uiz : item_pt->F3; + const MYCOMPLEX (*F3)[GRT_SRC_M_NUM][GRT_QWV_NUM] = (isuiz)? item_pt->F3_uiz : item_pt->F3; // 使用F(k)*sqrt(k)来衡量积分值,这可以平衡后续计算F(k)*Jm(kr)k积分时的系数 MYREAL sk[3]; @@ -136,13 +136,13 @@ static MYCOMPLEX simpson(const KInterval *item_pt, MYINT im, MYINT iqwv, bool is } /** 比较QWV的最大绝对值 */ -static void get_maxabsQWV(const MYCOMPLEX F[SRC_M_NUM][QWV_NUM], MYREAL maxabsF[GTYPES_MAX]){ +static void get_maxabsQWV(const MYCOMPLEX F[GRT_SRC_M_NUM][GRT_QWV_NUM], MYREAL maxabsF[GRT_GTYPES_MAX]){ MYREAL tmp; - for(MYINT i=0; i maxabsF[SRC_M_GTYPES[i]]){ - maxabsF[SRC_M_GTYPES[i]] = tmp; + if(tmp > maxabsF[GRT_SRC_M_GTYPES[i]]){ + maxabsF[GRT_SRC_M_GTYPES[i]] = tmp; } } } @@ -152,14 +152,14 @@ static void get_maxabsQWV(const MYCOMPLEX F[SRC_M_NUM][QWV_NUM], MYREAL maxabsF[ /** 检查区间是否符合要求,返回True表示通过 */ static bool check_fit( const KInterval *ptKitv, const KInterval *ptKitvL, const KInterval *ptKitvR, MYREAL kref, - bool isuiz, MYREAL maxabsQWV[GTYPES_MAX], MYREAL tol) + bool isuiz, MYREAL maxabsQWV[GRT_GTYPES_MAX], MYREAL tol) { // 计算积分差异 MYCOMPLEX S11, S12, S21, S22; // 核函数 - const MYCOMPLEX (*F3L)[SRC_M_NUM][QWV_NUM] = (isuiz)? ptKitvL->F3_uiz : ptKitvL->F3; - const MYCOMPLEX (*F3R)[SRC_M_NUM][QWV_NUM] = (isuiz)? ptKitvR->F3_uiz : ptKitvR->F3; + const MYCOMPLEX (*F3L)[GRT_SRC_M_NUM][GRT_QWV_NUM] = (isuiz)? ptKitvL->F3_uiz : ptKitvL->F3; + const MYCOMPLEX (*F3R)[GRT_SRC_M_NUM][GRT_QWV_NUM] = (isuiz)? ptKitvR->F3_uiz : ptKitvR->F3; // 取近似积分 \int_k1^k2 k^0.5 dk MYREAL kcoef13 = RTWOTHIRD*( ptKitv->k3[2]*sqrt(ptKitv->k3[2]) - ptKitv->k3[0]*sqrt(ptKitv->k3[0]) ); @@ -168,16 +168,16 @@ static bool check_fit( MYREAL S_dif, S_ref; bool badtol = false; - for(MYINT im=0; imk3[0] > kref && qwvchs[c]!='v') continue; + if(ptKitv->k3[0] > kref && GRT_QWV_CODES[c]!='v') continue; S11 = simpson(ptKitv, im, c, isuiz, SIMPSON_A_ApH); S12 = simpson(ptKitv, im, c, isuiz, SIMPSON_ApH_Ap2H); @@ -239,26 +239,26 @@ static bool check_fit( static void interv_integ( const KInterval *ptKitv, MYINT nr, MYREAL *rs, - MYCOMPLEX sum_J[nr][SRC_M_NUM][INTEG_NUM], + MYCOMPLEX sum_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], bool calc_upar, - MYCOMPLEX sum_uiz_J[nr][SRC_M_NUM][INTEG_NUM], - MYCOMPLEX sum_uir_J[nr][SRC_M_NUM][INTEG_NUM]) + MYCOMPLEX sum_uiz_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], + MYCOMPLEX sum_uir_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM]) { - MYCOMPLEX SUM[SRC_M_NUM][INTEG_NUM]={0}; + MYCOMPLEX SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM]={0}; // 震中距rs循环 for(MYINT ir=0; irk3, rs[ir], ptKitv->F3, false, SUM); - for(MYINT i=0; ik3, rs[ir], ptKitv->F3_uiz, false, SUM); - for(MYINT i=0; ik3, rs[ir], ptKitv->F3, true, SUM); - for(MYINT i=0; it1 = grt_compute_travt1d(mod1d->Thk, mod1d->Vb, mod1d->n, mod1d->isrc, mod1d->ircv, dist); strcpy(pt_hd->kt1, "S"); - for(int im=0; im=3) continue; - int modr = SRC_M_ORDERS[im]; + int modr = GRT_SRC_M_ORDERS[im]; int sgn=1; // 用于反转Z分量 - for(int c=0; cdt, pt_fh->nt, delayT0); diff --git a/pygrt/C_extension/src/dynamic/grn.c b/pygrt/C_extension/src/dynamic/grn.c index b9a8a0c2..303fe509 100755 --- a/pygrt/C_extension/src/dynamic/grn.c +++ b/pygrt/C_extension/src/dynamic/grn.c @@ -44,20 +44,20 @@ * @param[out] grn 三分量频谱 */ static void recordin_GRN( - MYINT iw, MYINT nr, MYCOMPLEX coef, MYCOMPLEX sum_J[nr][SRC_M_NUM][INTEG_NUM], - MYCOMPLEX *grn[nr][SRC_M_NUM][CHANNEL_NUM] + MYINT iw, MYINT nr, MYCOMPLEX coef, MYCOMPLEX sum_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], + MYCOMPLEX *grn[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM] ) { // 局部变量,将某个频点的格林函数谱临时存放 - MYCOMPLEX (*tmp_grn)[SRC_M_NUM][CHANNEL_NUM] = (MYCOMPLEX(*)[SRC_M_NUM][CHANNEL_NUM])calloc(nr, sizeof(*tmp_grn)); + MYCOMPLEX (*tmp_grn)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM] = (MYCOMPLEX(*)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM])calloc(nr, sizeof(*tmp_grn)); for(MYINT ir=0; irRho[mod1d->isrc]; // 震源区密度 - const MYREAL fac = RONE/(RFOUR*PI*Rho); - const MYREAL hs = GRT_MAX(fabs(mod1d->depsrc - mod1d->deprcv), MIN_DEPTH_GAP_SRC_RCV); // hs=max(震源和台站深度差,1.0) + const MYREAL fac = 1.0/(4.0*PI*Rho); + const MYREAL hs = GRT_MAX(fabs(mod1d->depsrc - mod1d->deprcv), GRT_MIN_DEPTH_GAP_SRC_RCV); // hs=max(震源和台站深度差,1.0) // 乘相应系数 k0 *= PI/hs; const MYREAL k02 = k0*k0; const MYREAL ampk2 = ampk*ampk; - if(vmin_ref < RZERO) keps = RZERO; // 若使用峰谷平均法,则不使用keps进行收敛判断 + if(vmin_ref < 0.0) keps = 0.0; // 若使用峰谷平均法,则不使用keps进行收敛判断 - bool useFIM = (filonLength > RZERO) || (safilonTol > RZERO) ; // 是否使用Filon积分(包括自适应Filon) + bool useFIM = (filonLength > 0.0) || (safilonTol > 0.0) ; // 是否使用Filon积分(包括自适应Filon) const MYREAL dk=PI2/(Length*rmax); // 波数积分间隔 - const MYREAL filondk = (filonLength > RZERO) ? PI2/(filonLength*rmax) : RZERO; // Filon积分间隔 + const MYREAL filondk = (filonLength > 0.0) ? PI2/(filonLength*rmax) : 0.0; // Filon积分间隔 const MYREAL filonK = filonCut/rmax; // 波数积分和Filon积分的分割点 // PTAM的积分中间结果, 每个震中距两个文件,因为PTAM对不同震中距使用不同的dk // 在文件名后加后缀,区分不同震中距 char **ptam_fstatsdir = (char**)calloc(nr, sizeof(char*)); - if(statsstr!=NULL && nstatsidxs > 0 && vmin_ref < RZERO){ + if(statsstr!=NULL && nstatsidxs > 0 && vmin_ref < 0.0){ for(MYINT ir=0; ir RZERO){ + if(useFIM && freq_invstats[iw]==GRT_INVERSE_SUCCESS){ + if(filondk > 0.0){ // 基于线性插值的Filon积分,固定采样间隔 k = grt_linear_filon_integ( local_mod1d, k, dk, filondk, kmax, keps, omega, nr, rs, sum_J, calc_upar, sum_uiz_J, sum_uir_J, fstats, grt_kernel, &freq_invstats[iw]); } - else if(safilonTol > RZERO){ + else if(safilonTol > 0.0){ // 基于自适应采样的Filon积分 k = grt_sa_filon_integ( local_mod1d, fabs(vmin_ref)/ampk, k, dk, safilonTol, kmax, omega, nr, rs, @@ -232,7 +232,7 @@ void grt_integ_grn_spec( } // k之后的部分使用峰谷平均法进行显式收敛,建议在浅源地震的时候使用 - if(vmin_ref < RZERO && freq_invstats[iw]==INVERSE_SUCCESS){ + if(vmin_ref < 0.0 && freq_invstats[iw]==GRT_INVERSE_SUCCESS){ grt_PTA_method( local_mod1d, k, dk, omega, nr, rs, sum_J, calc_upar, sum_uiz_J, sum_uir_J, @@ -244,7 +244,7 @@ void grt_integ_grn_spec( // 记录到格林函数结构体内 // 如果计算核函数过程中存在除零错误,则放弃该频率【通常在大震中距的低频段】 - if(freq_invstats[iw]==INVERSE_SUCCESS){ + if(freq_invstats[iw]==GRT_INVERSE_SUCCESS){ recordin_GRN(iw, nr, coef, sum_J, grn); if(calc_upar){ recordin_GRN(iw, nr, coef, sum_uiz_J, grn_uiz); @@ -287,7 +287,7 @@ void grt_integ_grn_spec( // 打印freq_invstats for(MYINT iw=nf1; iw<=nf2; ++iw){ - if(freq_invstats[iw]==INVERSE_FAILURE){ + if(freq_invstats[iw]==GRT_INVERSE_FAILURE){ fprintf(stderr, "iw=%d, freq=%e(Hz), meet Zero Divison Error, results are filled with 0.\n", iw, freqs[iw]); } } diff --git a/pygrt/C_extension/src/dynamic/grt_greenfn.c b/pygrt/C_extension/src/dynamic/grt_greenfn.c index 9e2e49ff..7fcef708 100644 --- a/pygrt/C_extension/src/dynamic/grt_greenfn.c +++ b/pygrt/C_extension/src/dynamic/grt_greenfn.c @@ -356,7 +356,7 @@ printf("\n" " will be the minimum velocity\n" " of model, but limited to %.1f. and if the \n", GRT_GREENFN_V_VMIN_REF); printf( " depth gap between source and receiver is \n" -" thinner than %.1f km, PTAM will be appled\n", MIN_DEPTH_GAP_SRC_RCV); printf( +" thinner than %.1f km, PTAM will be appled\n", GRT_MIN_DEPTH_GAP_SRC_RCV); printf( " automatically.\n" " + manually set POSITIVE value. \n" " + manually set NEGATIVE value, \n" @@ -380,7 +380,7 @@ printf("\n" " : designed to give residual k at\n" " 0 frequency, default is %.1f, and \n", GRT_GREENFN_K_K0); printf( " multiply PI/hs in program, \n" -" where hs = max(fabs(depsrc-deprcv), %.1f).\n", MIN_DEPTH_GAP_SRC_RCV); printf( +" where hs = max(fabs(depsrc-deprcv), %.1f).\n", GRT_MIN_DEPTH_GAP_SRC_RCV); printf( " : amplification factor, default is %.2f.\n", GRT_GREENFN_K_AMPK); printf( " : a threshold for break wavenumber \n" " integration in advance. See \n" @@ -725,7 +725,7 @@ int greenfn_main(int argc, char **argv) { } // 如果没有主动设置vmin_ref,则判断是否要自动使用PTAM - if( !Ctrl->V.active && fabs(Ctrl->D.deprcv - Ctrl->D.depsrc) <= MIN_DEPTH_GAP_SRC_RCV) { + if( !Ctrl->V.active && fabs(Ctrl->D.deprcv - Ctrl->D.depsrc) <= GRT_MIN_DEPTH_GAP_SRC_RCV) { Ctrl->V.vmin_ref = - fabs(Ctrl->V.vmin_ref); } @@ -782,13 +782,13 @@ int greenfn_main(int argc, char **argv) { } // 建立格林函数的complex数组 - MYCOMPLEX *(*grn)[SRC_M_NUM][CHANNEL_NUM] = (MYCOMPLEX*(*)[SRC_M_NUM][CHANNEL_NUM]) calloc(Ctrl->R.nr, sizeof(*grn)); - MYCOMPLEX *(*grn_uiz)[SRC_M_NUM][CHANNEL_NUM] = (Ctrl->e.active)? (MYCOMPLEX*(*)[SRC_M_NUM][CHANNEL_NUM]) calloc(Ctrl->R.nr, sizeof(*grn_uiz)) : NULL; - MYCOMPLEX *(*grn_uir)[SRC_M_NUM][CHANNEL_NUM] = (Ctrl->e.active)? (MYCOMPLEX*(*)[SRC_M_NUM][CHANNEL_NUM]) calloc(Ctrl->R.nr, sizeof(*grn_uir)) : NULL; + MYCOMPLEX *(*grn)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM] = (MYCOMPLEX*(*)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]) calloc(Ctrl->R.nr, sizeof(*grn)); + MYCOMPLEX *(*grn_uiz)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM] = (Ctrl->e.active)? (MYCOMPLEX*(*)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]) calloc(Ctrl->R.nr, sizeof(*grn_uiz)) : NULL; + MYCOMPLEX *(*grn_uir)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM] = (Ctrl->e.active)? (MYCOMPLEX*(*)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]) calloc(Ctrl->R.nr, sizeof(*grn_uir)) : NULL; for(int ir=0; irR.nr; ++ir){ - for(int i=0; iN.nf, sizeof(MYCOMPLEX)); if(grn_uiz) grn_uiz[ir][i][c] = (MYCOMPLEX*)calloc(Ctrl->N.nf, sizeof(MYCOMPLEX)); if(grn_uir) grn_uir[ir][i][c] = (MYCOMPLEX*)calloc(Ctrl->N.nf, sizeof(MYCOMPLEX)); @@ -826,7 +826,7 @@ int greenfn_main(int argc, char **argv) { Ctrl->R.nr, Ctrl->R.s_rs, Ctrl->R.rs, travtPS, Ctrl->D.depsrc, Ctrl->D.deprcv, Ctrl->E.delayT0, Ctrl->E.delayV0, Ctrl->e.active, Ctrl->G.doEX, Ctrl->G.doVF, Ctrl->G.doHF, Ctrl->G.doDC, - ZRTchs, grn, grn_uiz, grn_uir); + GRT_ZRT_CODES, grn, grn_uiz, grn_uir); // 打印走时 @@ -843,8 +843,8 @@ int greenfn_main(int argc, char **argv) { // 释放内存 for(int ir=0; irR.nr; ++ir){ - for(int i=0; iN.active; // 根据参数设置,选择分量名 - const char *chs = (rot2ZNE)? ZNEchs : ZRTchs; + const char *chs = (rot2ZNE)? GRT_ZNE_CODES : GRT_ZRT_CODES; float **ptarrout=NULL, *arrout=NULL; float *arrsyn[3] = {NULL, NULL, NULL}; @@ -686,8 +686,8 @@ int syn_main(int argc, char **argv){ // 重新计算方向因子 grt_set_source_radiation(Ctrl->srcRadi, Ctrl->computeType, (ityp==3), Ctrl->S.M0, upar_scale, Ctrl->A.azrad, Ctrl->mchn); - for(int c=0; csrcRadi[k][c]; if(coef == 0.0) continue; char *buffer = NULL; if(ityp==0 || ityp==3){ - GRT_SAFE_ASPRINTF(&buffer, "%s/%s%c.sac", Ctrl->G.s_grnpath, SRC_M_NAME_ABBR[k], ch); + GRT_SAFE_ASPRINTF(&buffer, "%s/%s%c.sac", Ctrl->G.s_grnpath, GRT_SRC_M_NAME_ABBR[k], ch); } else { - GRT_SAFE_ASPRINTF(&buffer, "%s/%c%s%c.sac", Ctrl->G.s_grnpath, tolower(ZRTchs[ityp-1]), SRC_M_NAME_ABBR[k], ch); + GRT_SAFE_ASPRINTF(&buffer, "%s/%c%s%c.sac", Ctrl->G.s_grnpath, tolower(GRT_ZRT_CODES[ityp-1]), GRT_SRC_M_NAME_ABBR[k], ch); } float *arr = grt_read_SAC(command, buffer, pthd, NULL); @@ -795,11 +795,11 @@ int syn_main(int argc, char **argv){ } // 保存到SAC文件 - for(int i1=0; i1e.active){ - for(int i2=0; i2mu[iy]; mod1d_kaka1 = mod1d->kaka[iy]; mod1d_kbkb1 = mod1d->kbkb[iy]; - mod1d_xa1 = sqrt(RONE - mod1d_kaka1/(k*k)); - mod1d_xb1 = sqrt(RONE - mod1d_kbkb1/(k*k)); + mod1d_xa1 = sqrt(1.0 - mod1d_kaka1/(k*k)); + mod1d_xb1 = sqrt(1.0 - mod1d_kbkb1/(k*k)); // mod1d_xa1 = sqrt(k*k - mod1d_kaka1)/k; // mod1d_xb1 = sqrt(k*k - mod1d_kbkb1)/k; @@ -173,7 +173,7 @@ void grt_kernel( // omega, k, // RD, pRDL, RU, pRUL, // TD, pTDL, TU, pTUL, stats); - if(*stats==INVERSE_FAILURE) goto BEFORE_RETURN; + if(*stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; } #if Print_GRTCOEF == 1 @@ -214,7 +214,7 @@ void grt_kernel( TD, TDL, TU, TUL, RD_FA, &RDL_FA, RU_FA, &RUL_FA, TD_FA, &TDL_FA, TU_FA, &TUL_FA, stats); - if(*stats==INVERSE_FAILURE) goto BEFORE_RETURN; + if(*stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; } } else if(iy==imin){ // 虚拟层位,可对递推公式简化 @@ -236,7 +236,7 @@ void grt_kernel( TD, TDL, TU, TUL, RD_RS, &RDL_RS, RU_RS, &RUL_RS, TD_RS, &TDL_RS, TU_RS, &TUL_RS, stats); // 写入原地址 - if(*stats==INVERSE_FAILURE) goto BEFORE_RETURN; + if(*stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; } } else if(iy==imax){ // 虚拟层位,可对递推公式简化 @@ -259,7 +259,7 @@ void grt_kernel( TD, TDL, TU, TUL, RD_BL, &RDL_BL, RU_BL, &RUL_BL, TD_BL, &TDL_BL, TU_BL, &TUL_BL, stats); // 写入原地址 - if(*stats==INVERSE_FAILURE) goto BEFORE_RETURN; + if(*stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; } } // END if @@ -271,8 +271,8 @@ void grt_kernel( // 计算震源系数 - MYCOMPLEX src_coef_PSV[SRC_M_NUM][QWV_NUM-1][2] = {0}; - MYCOMPLEX src_coef_SH[SRC_M_NUM][2] = {0}; + MYCOMPLEX src_coef_PSV[GRT_SRC_M_NUM][GRT_QWV_NUM-1][2] = {0}; + MYCOMPLEX src_coef_SH[GRT_SRC_M_NUM][2] = {0}; grt_source_coef_PSV(src_xa, src_xb, src_kaka, src_kbkb, k, src_coef_PSV); grt_source_coef_SH(src_xb, src_kbkb, k, src_coef_SH); @@ -282,15 +282,15 @@ void grt_kernel( // 递推RU_FA grt_calc_R_tilt_PSV(top_xa, top_xb, top_kbkb, k, R_tilt, stats); - if(*stats==INVERSE_FAILURE) goto BEFORE_RETURN; + if(*stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; grt_recursion_RU( - R_tilt, RONE, + R_tilt, 1.0, RD_FA, RDL_FA, RU_FA, RUL_FA, TD_FA, TDL_FA, TU_FA, TUL_FA, RU_FA, &RUL_FA, NULL, NULL, stats); - if(*stats==INVERSE_FAILURE) goto BEFORE_RETURN; + if(*stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; // 根据震源和台站相对位置,计算最终的系数 if(ircvup){ // A接收 B震源 @@ -307,7 +307,7 @@ void grt_kernel( TD_RS, TDL_RS, TU_RS, TUL_RS, RU_FB, &RUL_FB, inv_2x2T, &invT, stats); - if(*stats==INVERSE_FAILURE) goto BEFORE_RETURN; + if(*stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; #if Print_GRTCOEF == 1 // TEST------------------------------------------------------------- @@ -347,16 +347,16 @@ void grt_kernel( grt_cmat2x2_mul(RD_BL, RU_FB, tmpR2); grt_cmat2x2_one_sub(tmpR2); grt_cmat2x2_inv(tmpR2, tmpR2, stats);// (I - xx)^-1 - if(*stats==INVERSE_FAILURE) goto BEFORE_RETURN; + if(*stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; grt_cmat2x2_mul(inv_2x2T, tmpR2, tmp2x2); if(calc_uiz) grt_cmat2x2_assign(tmp2x2, tmp2x2_uiz); // 为后续计算空间导数备份 grt_cmat2x2_mul(R_EV, tmp2x2, tmp2x2); - tmpRL = invT / (RONE - RDL_BL * RUL_FB); + tmpRL = invT / (1.0 - RDL_BL * RUL_FB); tmpRL2 = R_EVL * tmpRL; - for(MYINT i=0; iDep[mod1d->isrc] == RZERO && mod1d->Dep[mod1d->ircv] == RZERO) + if(mod1d->Dep[mod1d->isrc] == 0.0 && mod1d->Dep[mod1d->ircv] == 0.0) { - for(MYINT c=0; c will be the minimum velocity\n" " of model, but limited to %.1f. and if the \n", GRT_GREENFN_V_VMIN_REF); printf( " depth gap between source and receiver is \n" -" thinner than %.1f km, PTAM will be appled\n", MIN_DEPTH_GAP_SRC_RCV); printf( +" thinner than %.1f km, PTAM will be appled\n", GRT_MIN_DEPTH_GAP_SRC_RCV); printf( " automatically.\n" " + manually set POSITIVE value. \n" " + manually set NEGATIVE value, \n" @@ -197,7 +197,7 @@ printf("\n" " bound is k0,\n" " : default is %.1f, and \n", GRT_GREENFN_K_K0); printf( " multiply PI/hs in program, \n" -" where hs = max(fabs(depsrc-deprcv), %.1f).\n", MIN_DEPTH_GAP_SRC_RCV); printf( +" where hs = max(fabs(depsrc-deprcv), %.1f).\n", GRT_MIN_DEPTH_GAP_SRC_RCV); printf( " : a threshold for break wavenumber \n" " integration in advance. See \n" " (Yao and Harkrider, 1983) for details.\n" @@ -399,13 +399,13 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ * @param[in] prefix 前缀字符串 */ static void print_grn_title(const char *prefix){ - for(int i=0; iV.active && fabs(Ctrl->D.deprcv - Ctrl->D.depsrc) <= MIN_DEPTH_GAP_SRC_RCV) { + if( !Ctrl->V.active && fabs(Ctrl->D.deprcv - Ctrl->D.depsrc) <= GRT_MIN_DEPTH_GAP_SRC_RCV) { Ctrl->V.vmin_ref = - fabs(Ctrl->V.vmin_ref); } @@ -476,9 +476,9 @@ int static_greenfn_main(int argc, char **argv){ } // 建立格林函数的浮点数 - MYREAL (*grn)[SRC_M_NUM][CHANNEL_NUM] = (MYREAL (*)[SRC_M_NUM][CHANNEL_NUM]) calloc(Ctrl->nr, sizeof(*grn)); - MYREAL (*grn_uiz)[SRC_M_NUM][CHANNEL_NUM] = (Ctrl->e.active)? (MYREAL (*)[SRC_M_NUM][CHANNEL_NUM]) calloc(Ctrl->nr, sizeof(*grn_uiz)) : NULL; - MYREAL (*grn_uir)[SRC_M_NUM][CHANNEL_NUM] = (Ctrl->e.active)? (MYREAL (*)[SRC_M_NUM][CHANNEL_NUM]) calloc(Ctrl->nr, sizeof(*grn_uir)) : NULL; + MYREAL (*grn)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM] = (MYREAL (*)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]) calloc(Ctrl->nr, sizeof(*grn)); + MYREAL (*grn_uiz)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM] = (Ctrl->e.active)? (MYREAL (*)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]) calloc(Ctrl->nr, sizeof(*grn_uiz)) : NULL; + MYREAL (*grn_uir)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM] = (Ctrl->e.active)? (MYREAL (*)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]) calloc(Ctrl->nr, sizeof(*grn_uir)) : NULL; //============================================================================== diff --git a/pygrt/C_extension/src/static/grt_static_rotation.c b/pygrt/C_extension/src/static/grt_static_rotation.c index 2587c6a1..498f7b5d 100644 --- a/pygrt/C_extension/src/static/grt_static_rotation.c +++ b/pygrt/C_extension/src/static/grt_static_rotation.c @@ -128,7 +128,7 @@ int static_rotation_main(int argc, char **argv){ GRT_SAFE_FREE_PTR(copyline); // 指示特定的通道名 - chs = (rot2ZNE)? ZNEchs : ZRTchs; + chs = (rot2ZNE)? GRT_ZNE_CODES : GRT_ZRT_CODES; // 想合成位移空间导数但输入的格林函数没有 if(ncols < 14){ diff --git a/pygrt/C_extension/src/static/grt_static_strain.c b/pygrt/C_extension/src/static/grt_static_strain.c index 6d16ffe3..d68e8213 100644 --- a/pygrt/C_extension/src/static/grt_static_strain.c +++ b/pygrt/C_extension/src/static/grt_static_strain.c @@ -126,7 +126,7 @@ int static_strain_main(int argc, char **argv){ GRT_SAFE_FREE_PTR(copyline); // 指示特定的通道名 - chs = (rot2ZNE)? ZNEchs : ZRTchs; + chs = (rot2ZNE)? GRT_ZNE_CODES : GRT_ZRT_CODES; // 想合成位移空间导数但输入的格林函数没有 if(ncols < 14){ diff --git a/pygrt/C_extension/src/static/grt_static_stress.c b/pygrt/C_extension/src/static/grt_static_stress.c index 84256520..4f4c784e 100644 --- a/pygrt/C_extension/src/static/grt_static_stress.c +++ b/pygrt/C_extension/src/static/grt_static_stress.c @@ -137,7 +137,7 @@ int static_stress_main(int argc, char **argv){ GRT_SAFE_FREE_PTR(copyline); // 指示特定的通道名 - chs = (rot2ZNE)? ZNEchs : ZRTchs; + chs = (rot2ZNE)? GRT_ZNE_CODES : GRT_ZRT_CODES; // 想合成位移空间导数但输入的格林函数没有 if(ncols < 14){ diff --git a/pygrt/C_extension/src/static/grt_static_syn.c b/pygrt/C_extension/src/static/grt_static_syn.c index 840e7146..87b07206 100644 --- a/pygrt/C_extension/src/static/grt_static_syn.c +++ b/pygrt/C_extension/src/static/grt_static_syn.c @@ -46,10 +46,10 @@ typedef struct { } e; // 存储不同震源的震源机制相关参数的数组 - MYREAL mchn[MECHANISM_NUM]; + MYREAL mchn[GRT_MECHANISM_NUM]; // 方向因子数组 - MYREAL srcRadi[SRC_M_NUM][CHANNEL_NUM]; + MYREAL srcRadi[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]; // 最终要计算的震源类型 MYINT computeType; @@ -265,17 +265,17 @@ int static_syn_main(int argc, char **argv){ getopt_from_command(Ctrl, argc, argv); // 辐射因子 - // double srcRadi[SRC_M_NUM][CHANNEL_NUM]={0}; + // double srcRadi[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]={0}; // 从标准输入中读取静态格林函数表 - double x0, y0, grn[SRC_M_NUM][CHANNEL_NUM]={0}, syn[CHANNEL_NUM]={0}, syn_upar[CHANNEL_NUM][CHANNEL_NUM]={0}; - double grn_uiz[SRC_M_NUM][CHANNEL_NUM]={0}, grn_uir[SRC_M_NUM][CHANNEL_NUM]={0}; + double x0, y0, grn[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]={0}, syn[GRT_CHANNEL_NUM]={0}, syn_upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]={0}; + double grn_uiz[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]={0}, grn_uir[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]={0}; // 输出分量格式,即是否需要旋转到ZNE bool rot2ZNE = Ctrl->N.active; // 根据参数设置,选择分量名 - const char *chs = (rot2ZNE)? ZNEchs : ZRTchs; + const char *chs = (rot2ZNE)? GRT_ZNE_CODES : GRT_ZRT_CODES; // 建立一个指针数组,方便读取多列数据 @@ -287,9 +287,9 @@ int static_syn_main(int argc, char **argv){ *(pt++) = &x0; *(pt++) = &y0; for(int m=0; m<3; ++m){ - for(int k=0; ks_computeType, toupper(chs[i])); fprintf(stdout, GRT_STRING_FMT, s_channel); } if(Ctrl->e.active){ - for(int k=0; ks_computeType, toupper(chs[i])); fprintf(stdout, GRT_STRING_FMT, s_channel); } @@ -418,8 +418,8 @@ int static_syn_main(int argc, char **argv){ printHead = true; } - double (*grn3)[CHANNEL_NUM]; // 使用对应类型的格林函数 - double tmpsyn[CHANNEL_NUM]; + double (*grn3)[GRT_CHANNEL_NUM]; // 使用对应类型的格林函数 + double tmpsyn[GRT_CHANNEL_NUM]; for(int ityp=0; itypsrcRadi, Ctrl->computeType, ityp==3, Ctrl->S.M0, upar_scale, azrad, Ctrl->mchn); - for(int i=0; isrcRadi[k][i]; } } // 保存数据 - for(int i=0; ie.active){ - for(int i=0; idepsrc - mod1d->deprcv), MIN_DEPTH_GAP_SRC_RCV); // hs=max(震源和台站深度差,1.0) + const MYREAL hs = GRT_MAX(fabs(mod1d->depsrc - mod1d->deprcv), GRT_MIN_DEPTH_GAP_SRC_RCV); // hs=max(震源和台站深度差,1.0) // 乘相应系数 k0 *= PI/hs; - if(vmin_ref < RZERO) keps = -RONE; // 若使用峰谷平均法,则不使用keps进行收敛判断 + if(vmin_ref < 0.0) keps = -1.0; // 若使用峰谷平均法,则不使用keps进行收敛判断 MYREAL k=0.0; - bool useFIM = (filonLength > RZERO) || (safilonTol > RZERO) ; // 是否使用Filon积分(包括自适应Filon) + bool useFIM = (filonLength > 0.0) || (safilonTol > 0.0) ; // 是否使用Filon积分(包括自适应Filon) const MYREAL dk=fabs(PI2/(Length*rmax)); // 波数积分间隔 - const MYREAL filondk = (filonLength > RZERO) ? PI2/(filonLength*rmax) : RZERO; // Filon积分间隔 + const MYREAL filondk = (filonLength > 0.0) ? PI2/(filonLength*rmax) : 0.0; // Filon积分间隔 const MYREAL filonK = filonCut/rmax; // 波数积分和Filon积分的分割点 const MYREAL kmax = k0; // 求和 sum F(ki,w)Jm(ki*r)ki // 关于形状详见int_Pk()函数内的注释 - MYCOMPLEX (*sum_J)[SRC_M_NUM][INTEG_NUM] = (MYCOMPLEX(*)[SRC_M_NUM][INTEG_NUM])calloc(nr, sizeof(*sum_J)); - MYCOMPLEX (*sum_uiz_J)[SRC_M_NUM][INTEG_NUM] = (calc_upar)? (MYCOMPLEX(*)[SRC_M_NUM][INTEG_NUM])calloc(nr, sizeof(*sum_uiz_J)) : NULL; - MYCOMPLEX (*sum_uir_J)[SRC_M_NUM][INTEG_NUM] = (calc_upar)? (MYCOMPLEX(*)[SRC_M_NUM][INTEG_NUM])calloc(nr, sizeof(*sum_uir_J)) : NULL; + MYCOMPLEX (*sum_J)[GRT_SRC_M_NUM][GRT_INTEG_NUM] = (MYCOMPLEX(*)[GRT_SRC_M_NUM][GRT_INTEG_NUM])calloc(nr, sizeof(*sum_J)); + MYCOMPLEX (*sum_uiz_J)[GRT_SRC_M_NUM][GRT_INTEG_NUM] = (calc_upar)? (MYCOMPLEX(*)[GRT_SRC_M_NUM][GRT_INTEG_NUM])calloc(nr, sizeof(*sum_uiz_J)) : NULL; + MYCOMPLEX (*sum_uir_J)[GRT_SRC_M_NUM][GRT_INTEG_NUM] = (calc_upar)? (MYCOMPLEX(*)[GRT_SRC_M_NUM][GRT_INTEG_NUM])calloc(nr, sizeof(*sum_uir_J)) : NULL; // 是否要输出积分过程文件 bool needfstats = (statsstr!=NULL); @@ -104,7 +104,7 @@ void grt_integ_static_grn( // PTAM的积分中间结果, 每个震中距两个文件,因为PTAM对不同震中距使用不同的dk // 在文件名后加后缀,区分不同震中距 char **ptam_fstatsdir = (char**)calloc(nr, sizeof(char*)); - if(needfstats && vmin_ref < RZERO){ + if(needfstats && vmin_ref < 0.0){ for(MYINT ir=0; ir RZERO){ + if(filondk > 0.0){ // 基于线性插值的Filon积分,固定采样间隔 k = grt_linear_filon_integ( mod1d, k, dk, filondk, kmax, keps, 0.0, nr, rs, sum_J, calc_upar, sum_uiz_J, sum_uir_J, fstats, grt_static_kernel, &inv_stats); } - else if(safilonTol > RZERO){ + else if(safilonTol > 0.0){ // 基于自适应采样的Filon积分 k = grt_sa_filon_integ( mod1d, kmax, k, dk, safilonTol, kmax, 0.0, nr, rs, @@ -179,7 +179,7 @@ void grt_integ_static_grn( } // k之后的部分使用峰谷平均法进行显式收敛,建议在浅源地震的时候使用 - if(vmin_ref < RZERO){ + if(vmin_ref < 0.0){ grt_PTA_method( mod1d, k, dk, 0.0, nr, rs, sum_J, calc_upar, sum_uiz_J, sum_uir_J, @@ -189,7 +189,7 @@ void grt_integ_static_grn( MYCOMPLEX src_mu = mod1d->mu[mod1d->isrc]; - MYCOMPLEX fac = dk * RONE/(RFOUR*PI * src_mu); + MYCOMPLEX fac = dk * 1.0/(4.0*PI * src_mu); // 将积分结果记录到浮点数数组中 recordin_GRN(nr, fac, sum_J, grn); diff --git a/pygrt/C_extension/src/static/static_layer.c b/pygrt/C_extension/src/static/static_layer.c index 2091f76b..509cd2b7 100644 --- a/pygrt/C_extension/src/static/static_layer.c +++ b/pygrt/C_extension/src/static/static_layer.c @@ -21,15 +21,15 @@ void grt_calc_static_R_tilt_PSV(MYCOMPLEX delta1, MYCOMPLEX R_tilt[2][2]){ // 公式(6.3.12) - R_tilt[0][0] = R_tilt[1][1] = CZERO; + R_tilt[0][0] = R_tilt[1][1] = 0.0; R_tilt[0][1] = -delta1; - R_tilt[1][0] = -RONE/delta1; + R_tilt[1][0] = -1.0/delta1; } void grt_calc_static_R_EV_PSV(bool ircvup, const MYCOMPLEX R[2][2], MYCOMPLEX R_EV[2][2]) { - MYCOMPLEX D11[2][2] = {{RONE, -RONE}, {RONE, RONE}}; - MYCOMPLEX D12[2][2] = {{RONE, -RONE}, {-RONE, -RONE}}; + MYCOMPLEX D11[2][2] = {{1.0, -1.0}, {1.0, 1.0}}; + MYCOMPLEX D12[2][2] = {{1.0, -1.0}, {-1.0, -1.0}}; // 公式(6.3.35,37) if(ircvup){// 震源更深 @@ -43,7 +43,7 @@ void grt_calc_static_R_EV_PSV(bool ircvup, const MYCOMPLEX R[2][2], MYCOMPLEX R_ void grt_calc_static_R_EV_SH(MYCOMPLEX RL, MYCOMPLEX *R_EVL) { - *R_EVL = (RONE + (RL)); + *R_EVL = (1.0 + (RL)); } void grt_calc_static_uiz_R_EV_PSV( @@ -51,7 +51,7 @@ void grt_calc_static_uiz_R_EV_PSV( const MYCOMPLEX R[2][2], MYCOMPLEX R_EV[2][2]) { // 新推导公式 - MYCOMPLEX kd2 = RTWO*k*delta1; + MYCOMPLEX kd2 = 2.0*k*delta1; MYCOMPLEX D11[2][2] = {{k, -k-kd2}, {k, k-kd2}}; MYCOMPLEX D12[2][2] = {{-k, k+kd2}, {k, k-kd2}}; if(ircvup){// 震源更深 @@ -67,9 +67,9 @@ void grt_calc_static_uiz_R_EV_SH(bool ircvup, MYREAL k, MYCOMPLEX RL, MYCOMPLEX { // 新推导公式 if(ircvup){// 震源更深 - *R_EVL = (RONE - (RL))*k; + *R_EVL = (1.0 - (RL))*k; } else { // 接收点更深 - *R_EVL = (RL - RONE)*k; + *R_EVL = (RL - 1.0)*k; } } @@ -96,26 +96,26 @@ void grt_calc_static_RT_PSV( // REFELCTION //------------------ RD ----------------------------------- - RD[0][0] = -RTWO*delta1*k*thk*dmu/A112 * ex2; - RD[0][1] = - ( RFOUR*del11*k2*thk2*A221*dmu + A112*B ) / (A221*A112) * ex2; + RD[0][0] = -2.0*delta1*k*thk*dmu/A112 * ex2; + RD[0][1] = - ( 4.0*del11*k2*thk2*A221*dmu + A112*B ) / (A221*A112) * ex2; RD[1][0] = - dmu/A112 * ex2; RD[1][1] = RD[0][0]; //------------------ RU ----------------------------------- - RU[0][0] = RZERO; + RU[0][0] = 0.0; RU[0][1] = B/A112; RU[1][0] = dmu/A221; - RU[1][1] = RZERO; + RU[1][1] = 0.0; // Transmission //------------------ TD ----------------------------------- - TD[0][0] = mu1*(RONE+delta1)/(A112) * ex; - TD[0][1] = RTWO*mu1*delta1*k*thk*(RONE+delta1)/(A112) * ex; - TD[1][0] = RZERO; + TD[0][0] = mu1*(1.0+delta1)/(A112) * ex; + TD[0][1] = 2.0*mu1*delta1*k*thk*(1.0+delta1)/(A112) * ex; + TD[1][0] = 0.0; TD[1][1] = TD[0][0]*A112/A221; //------------------ TU ----------------------------------- - TU[0][0] = mu2*(RONE+delta2)/A221 * ex; - TU[0][1] = RTWO*delta1*k*thk*mu2*(RONE+delta2)/A112 * ex; - TU[1][0] = RZERO; + TU[0][0] = mu2*(1.0+delta2)/A221 * ex; + TU[0][1] = 2.0*delta1*k*thk*mu2*(1.0+delta2)/A112 * ex; + TU[1][0] = 0.0; TU[1][1] = TU[0][0]*A221/A112; // printf("delta1=%.5e%+.5ej, delta2=%.5e%+.5ej, mu1=%.5e%+.5ej, mu2=%.5e%+.5ej, thk=%e, k=%e\n", @@ -147,6 +147,6 @@ void grt_calc_static_RT_SH( *RUL = - dmu/amu; // Transmission - *TDL = RTWO*mu1/amu * ex; + *TDL = 2.0*mu1/amu * ex; *TUL = (*TDL)*mu2/mu1; } \ No newline at end of file diff --git a/pygrt/C_extension/src/static/static_propagate.c b/pygrt/C_extension/src/static/static_propagate.c index 65774176..49dfce64 100644 --- a/pygrt/C_extension/src/static/static_propagate.c +++ b/pygrt/C_extension/src/static/static_propagate.c @@ -25,14 +25,14 @@ void grt_static_kernel( - const GRT_MODEL1D *mod1d, MYCOMPLEX omega, MYREAL k, MYCOMPLEX QWV[SRC_M_NUM][QWV_NUM], - bool calc_uiz, MYCOMPLEX QWV_uiz[SRC_M_NUM][QWV_NUM], MYINT *stats) + const GRT_MODEL1D *mod1d, MYCOMPLEX omega, MYREAL k, MYCOMPLEX QWV[GRT_SRC_M_NUM][GRT_QWV_NUM], + bool calc_uiz, MYCOMPLEX QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM], MYINT *stats) { // 初始化qwv为0 - for(MYINT i=0; i Date: Wed, 3 Sep 2025 13:52:57 +0800 Subject: [PATCH 3/5] REFAC: use macro for comment head '#', and add related helper function --- pygrt/C_extension/include/grt/common/const.h | 2 + pygrt/C_extension/include/grt/common/model.h | 21 +++++++++ pygrt/C_extension/include/grt/common/sacio2.h | 14 +++--- pygrt/C_extension/include/grt/common/util.h | 26 ++++++++++- pygrt/C_extension/src/common/iostats.c | 4 +- pygrt/C_extension/src/common/model.c | 35 ++++++++++++++- pygrt/C_extension/src/common/util.c | 43 +++++++++++++++++++ pygrt/C_extension/src/dynamic/signals.c | 13 +++--- .../src/static/grt_static_greenfn.c | 2 +- .../src/static/grt_static_rotation.c | 6 ++- .../src/static/grt_static_strain.c | 6 ++- .../src/static/grt_static_stress.c | 6 ++- pygrt/C_extension/src/static/grt_static_syn.c | 6 ++- 13 files changed, 156 insertions(+), 28 deletions(-) diff --git a/pygrt/C_extension/include/grt/common/const.h b/pygrt/C_extension/include/grt/common/const.h index 6233c769..0e0d02f6 100755 --- a/pygrt/C_extension/include/grt/common/const.h +++ b/pygrt/C_extension/include/grt/common/const.h @@ -77,6 +77,8 @@ typedef int MYINT; ///< 整数 #define GRT_STR_CMPLX_LONG_FMT "%51s" ///< 与复数格式同长度的字符串输出格式(更长) #define GRT_CMPLX_SPLIT(x) creal(x), cimag(x) ///< 用于打印复数时将实部虚部分开 +#define GRT_COMMENT_HEAD '#' ///< # 号, 作为注释的字符 + #define GRT_MAX(a, b) ((a) > (b) ? (a) : (b)) ///< 求两者较大值 #define GRT_MIN(a, b) ((a) < (b) ? (a) : (b)) ///< 求两者较小值 diff --git a/pygrt/C_extension/include/grt/common/model.h b/pygrt/C_extension/include/grt/common/model.h index f6fbfcd0..e6034a59 100755 --- a/pygrt/C_extension/include/grt/common/model.h +++ b/pygrt/C_extension/include/grt/common/model.h @@ -99,6 +99,27 @@ void grt_update_mod1d_omega(GRT_MODEL1D *mod1d, MYCOMPLEX omega); GRT_MODEL1D * grt_read_mod1d_from_file(const char *command, const char *modelpath, double depsrc, double deprcv, bool allowLiquid); +/** + * 从模型文件中判断各个量的大致精度(字符串长度),以确定浮点数输出位数 + * + * @param[in] command 命令名称 + * @param[in] modelpath 模型文件路径 + * @param[out] diglen 每一列的最大字符串长度 + * + */ +void grt_get_model_diglen_from_file(const char *command, const char *modelpath, MYINT diglen[6]); + +/** + * 浮点数比较,检查模型中是否存在该速度(不论Vp,Vs) + * + * @param[in] mod1d 模型 + * @param[in] vel 输入速度 + * @param[in] tol 浮点数比较精度 + * + * @return 是否存在 + */ +bool grt_check_vel_in_mod(const GRT_MODEL1D *mod1d, const MYREAL vel, const MYREAL tol); + /** * 计算最大最小速度(非零值) * diff --git a/pygrt/C_extension/include/grt/common/sacio2.h b/pygrt/C_extension/include/grt/common/sacio2.h index 1ab43f40..abd1aae1 100644 --- a/pygrt/C_extension/include/grt/common/sacio2.h +++ b/pygrt/C_extension/include/grt/common/sacio2.h @@ -14,9 +14,9 @@ /** * 读取SAC头段变量 * - * @param command (in)当前程序命令名称 - * @param name (in)SAC文件路径 - * @param hd (out)SAC头段变量结构体 + * @param[in] command 当前程序命令名称 + * @param[in] name SAC文件路径 + * @param[out] hd SAC头段变量结构体 */ void grt_read_SAC_HEAD(const char *command, const char *name, SACHEAD *hd); @@ -24,10 +24,10 @@ void grt_read_SAC_HEAD(const char *command, const char *name, SACHEAD *hd); /** * 读取SAC文件 * - * @param command (in)当前程序命令名称 - * @param name (in)SAC文件路径 - * @param hd (out)SAC头段变量结构体 - * @param arrout (inout)预分配内存,不需要则设为NULL + * @param[in] command 当前程序命令名称 + * @param[in] name SAC文件路径 + * @param[out] hd SAC头段变量结构体 + * @param[in,out] arrout 预分配内存,不需要则设为NULL * * @return 浮点数指针 */ diff --git a/pygrt/C_extension/include/grt/common/util.h b/pygrt/C_extension/include/grt/common/util.h index 4016950a..bb701cad 100644 --- a/pygrt/C_extension/include/grt/common/util.h +++ b/pygrt/C_extension/include/grt/common/util.h @@ -26,17 +26,41 @@ */ char ** grt_string_split(const char *string, const char *delim, int *size); +/** + * 指定分隔符,获得字符串的分割出的子字符串数。 + * 相当于是 grt_string_split 函数的简化版本 + * + * @param[in] string 原字符串 + * @param[in] delim 分隔符 + * + * @return 子字符串数 + */ +int grt_string_ncols(const char *string, const char* delim); /** * 从路径字符串中找到用/或\\分隔的最后一项 * - * @param path 路径字符串指针 + * @param[in] path 路径字符串指针 * * @return 指向最后一项字符串的指针 */ const char* grt_get_basename(const char* path); +/** + * 去除字符串首尾空白 + * + * @param[in,out] str 字符串 + */ +void grt_trim_whitespace(char* str); + + +/** + * 检查是否为注释行或空行 + * + * @param[in] line 读入一行的字符串 + */ +bool grt_is_comment_or_empty(const char* line); /** * 处理单个震中距对应的数据逆变换和SAC保存 diff --git a/pygrt/C_extension/src/common/iostats.c b/pygrt/C_extension/src/common/iostats.c index fe455a23..4ae2ab54 100755 --- a/pygrt/C_extension/src/common/iostats.c +++ b/pygrt/C_extension/src/common/iostats.c @@ -37,7 +37,7 @@ MYINT grt_extract_stats(FILE *bf0, FILE *af0){ // 打印标题 if(bf0 == NULL){ char K[20]; - snprintf(K, sizeof(K), GRT_STRING_FMT, "k"); K[0]='#'; + snprintf(K, sizeof(K), GRT_STRING_FMT, "k"); K[0]=GRT_COMMENT_HEAD; fprintf(af0, "%s", K); for(MYINT im=0; imn; ++i){ + if(fabs(vel - mod1d->Va[i])Vb[i]) #include #include +#include #include "grt/common/util.h" #include "grt/common/model.h" @@ -41,6 +42,23 @@ char ** grt_string_split(const char *string, const char *delim, int *size) } +int grt_string_ncols(const char *string, const char* delim){ + int count = 0; + + const char *str = string; + while (*str) { + // 跳过所有分隔符 + while (*str && strchr(delim, *str)) str++; + // 如果还有非分隔符字符,增加计数 + if (*str) count++; + // 跳过所有非分隔符字符 + while (*str && !strchr(delim, *str)) str++; + } + + return count; +} + + const char* grt_get_basename(const char* path) { // 找到最后一个 '/' char* last_slash = strrchr(path, '/'); @@ -60,6 +78,31 @@ const char* grt_get_basename(const char* path) { } +void grt_trim_whitespace(char* str) { + char* end; + + // 去除首部空白 + while (isspace((unsigned char)*str)) str++; + + // 去除尾部空白 + end = str + strlen(str) - 1; + while (end > str && isspace((unsigned char)*end)) end--; + + // 写入终止符 + *(end + 1) = '\0'; +} + + +bool grt_is_comment_or_empty(const char* line) { + // 跳过前导空白 + while (isspace((unsigned char)*line)) line++; + + // 检查是否为空行或注释行 + return (*line == '\0' || *line == GRT_COMMENT_HEAD); +} + + + /** * 将一条数据反变换回时间域再进行处理,保存到SAC文件 diff --git a/pygrt/C_extension/src/dynamic/signals.c b/pygrt/C_extension/src/dynamic/signals.c index efd668e4..1edb9c68 100644 --- a/pygrt/C_extension/src/dynamic/signals.c +++ b/pygrt/C_extension/src/dynamic/signals.c @@ -22,6 +22,7 @@ #include "grt/dynamic/signals.h" #include "grt/common/const.h" #include "grt/common/colorstr.h" +#include "grt/common/util.h" bool grt_check_tftype_tfparams(const char tftype, const char *tfparams){ @@ -324,16 +325,10 @@ float * grt_get_custom_wave(int *Nt, const char *tfparams){ // 逐行读入 char line[1024]; - char first_char = '\0'; int nt = 0; while(fgets(line, sizeof(line), fp)) { - first_char = '\0'; - // 读入第一个非空字符,判断是否是注释行 - if(sscanf(line, " %c", &first_char) == 1){ - if(first_char == '#') continue; - } else { - continue; // 空行,跳过 - } + // 注释行 + if(grt_is_comment_or_empty(line)) continue; tfarr = (float*)realloc(tfarr, sizeof(float)*(nt+1)); if(sscanf(line, " %f", &tfarr[nt]) < 1){ @@ -348,6 +343,8 @@ float * grt_get_custom_wave(int *Nt, const char *tfparams){ return NULL; } + fclose(fp); + *Nt = nt; return tfarr; } diff --git a/pygrt/C_extension/src/static/grt_static_greenfn.c b/pygrt/C_extension/src/static/grt_static_greenfn.c index b4bcbcc5..417ed4b5 100644 --- a/pygrt/C_extension/src/static/grt_static_greenfn.c +++ b/pygrt/C_extension/src/static/grt_static_greenfn.c @@ -504,7 +504,7 @@ int static_greenfn_main(int argc, char **argv){ // 输出标题 char XX[20]; - sprintf(XX, GRT_STRING_FMT, "X(km)"); XX[0]='#'; + sprintf(XX, GRT_STRING_FMT, "X(km)"); XX[0]=GRT_COMMENT_HEAD; fprintf(stdout, "%s", XX); fprintf(stdout, GRT_STRING_FMT, "Y(km)"); print_grn_title(""); diff --git a/pygrt/C_extension/src/static/grt_static_rotation.c b/pygrt/C_extension/src/static/grt_static_rotation.c index 498f7b5d..30a65b77 100644 --- a/pygrt/C_extension/src/static/grt_static_rotation.c +++ b/pygrt/C_extension/src/static/grt_static_rotation.c @@ -9,6 +9,7 @@ #include "grt/common/const.h" +#include "grt/common/util.h" #include "grt.h" @@ -135,7 +136,8 @@ int static_rotation_main(int argc, char **argv){ GRTRaiseError("[%s] Error! The input has no spatial derivatives. \n", command); } } - if(line[0] == '#') continue; + // 注释行 + if(grt_is_comment_or_empty(line)) continue; // 读取该行数据 char *copyline = strdup(line); @@ -156,7 +158,7 @@ int static_rotation_main(int argc, char **argv){ fprintf(stdout, "# "GRT_REAL_FMT" "GRT_REAL_FMT" "GRT_REAL_FMT"\n", rcv_va, rcv_vb, rcv_rho); char XX[20]; - sprintf(XX, GRT_STRING_FMT, "X(km)"); XX[0]='#'; + sprintf(XX, GRT_STRING_FMT, "X(km)"); XX[0]=GRT_COMMENT_HEAD; fprintf(stdout, "%s", XX); fprintf(stdout, GRT_STRING_FMT, "Y(km)"); char s_channel[15]; diff --git a/pygrt/C_extension/src/static/grt_static_strain.c b/pygrt/C_extension/src/static/grt_static_strain.c index d68e8213..94b67f3e 100644 --- a/pygrt/C_extension/src/static/grt_static_strain.c +++ b/pygrt/C_extension/src/static/grt_static_strain.c @@ -8,6 +8,7 @@ */ #include "grt/common/const.h" +#include "grt/common/util.h" #include "grt.h" @@ -133,7 +134,8 @@ int static_strain_main(int argc, char **argv){ GRTRaiseError("[%s] Error! The input has no spatial derivatives. \n", command); } } - if(line[0] == '#') continue; + // 注释行 + if(grt_is_comment_or_empty(line)) continue; // 读取该行数据 char *copyline = strdup(line); @@ -154,7 +156,7 @@ int static_strain_main(int argc, char **argv){ fprintf(stdout, "# "GRT_REAL_FMT" "GRT_REAL_FMT" "GRT_REAL_FMT"\n", rcv_va, rcv_vb, rcv_rho); char XX[20]; - sprintf(XX, GRT_STRING_FMT, "X(km)"); XX[0]='#'; + sprintf(XX, GRT_STRING_FMT, "X(km)"); XX[0]=GRT_COMMENT_HEAD; fprintf(stdout, "%s", XX); fprintf(stdout, GRT_STRING_FMT, "Y(km)"); char s_channel[15]; diff --git a/pygrt/C_extension/src/static/grt_static_stress.c b/pygrt/C_extension/src/static/grt_static_stress.c index 4f4c784e..6fe285ef 100644 --- a/pygrt/C_extension/src/static/grt_static_stress.c +++ b/pygrt/C_extension/src/static/grt_static_stress.c @@ -8,6 +8,7 @@ */ #include "grt/common/const.h" +#include "grt/common/util.h" #include "grt.h" @@ -144,7 +145,8 @@ int static_stress_main(int argc, char **argv){ GRTRaiseError("[%s] Error! The input has no spatial derivatives. \n", command); } } - if(line[0] == '#') continue; + // 注释行 + if(grt_is_comment_or_empty(line)) continue; // 读取该行数据 char *copyline = strdup(line); @@ -169,7 +171,7 @@ int static_stress_main(int argc, char **argv){ fprintf(stdout, "# "GRT_REAL_FMT" "GRT_REAL_FMT" "GRT_REAL_FMT"\n", rcv_va, rcv_vb, rcv_rho); char XX[20]; - sprintf(XX, GRT_STRING_FMT, "X(km)"); XX[0]='#'; + sprintf(XX, GRT_STRING_FMT, "X(km)"); XX[0]=GRT_COMMENT_HEAD; fprintf(stdout, "%s", XX); fprintf(stdout, GRT_STRING_FMT, "Y(km)"); char s_channel[15]; diff --git a/pygrt/C_extension/src/static/grt_static_syn.c b/pygrt/C_extension/src/static/grt_static_syn.c index 87b07206..81232f7c 100644 --- a/pygrt/C_extension/src/static/grt_static_syn.c +++ b/pygrt/C_extension/src/static/grt_static_syn.c @@ -11,6 +11,7 @@ #include "grt/common/const.h" #include "grt/common/radiation.h" #include "grt/common/coord.h" +#include "grt/common/util.h" #include "grt.h" @@ -372,7 +373,8 @@ int static_syn_main(int argc, char **argv){ GRTRaiseError("[%s] Error! The input has no spatial derivatives. \n", command); } } - if(line[0] == '#') continue; + // 注释行 + if(grt_is_comment_or_empty(line)) continue; // 读取该行数据 char *copyline = strdup(line); @@ -396,7 +398,7 @@ int static_syn_main(int argc, char **argv){ fprintf(stdout, "# "GRT_REAL_FMT" "GRT_REAL_FMT" "GRT_REAL_FMT"\n", rcv_va, rcv_vb, rcv_rho); char XX[20]; - sprintf(XX, GRT_STRING_FMT, "X(km)"); XX[0]='#'; + sprintf(XX, GRT_STRING_FMT, "X(km)"); XX[0]=GRT_COMMENT_HEAD; fprintf(stdout, "%s", XX); fprintf(stdout, GRT_STRING_FMT, "Y(km)"); char s_channel[5]; From 5835f2d59dc81ff6c0984d9690eb48f5db1a14bb Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Wed, 3 Sep 2025 14:54:26 +0800 Subject: [PATCH 4/5] FEAT: add `grt_getline` function --- pygrt/C_extension/include/grt/common/util.h | 8 +++ pygrt/C_extension/src/common/model.c | 10 ++-- pygrt/C_extension/src/common/util.c | 56 +++++++++++++++++++ pygrt/C_extension/src/dynamic/signals.c | 7 ++- .../src/static/grt_static_rotation.c | 7 ++- .../src/static/grt_static_strain.c | 7 ++- .../src/static/grt_static_stress.c | 7 ++- pygrt/C_extension/src/static/grt_static_syn.c | 8 ++- 8 files changed, 96 insertions(+), 14 deletions(-) diff --git a/pygrt/C_extension/include/grt/common/util.h b/pygrt/C_extension/include/grt/common/util.h index bb701cad..7af79f5f 100644 --- a/pygrt/C_extension/include/grt/common/util.h +++ b/pygrt/C_extension/include/grt/common/util.h @@ -62,6 +62,14 @@ void grt_trim_whitespace(char* str); */ bool grt_is_comment_or_empty(const char* line); + +/** + * 由于 Windows MSYS2 环境没有 getline 函数(即使定义了 _GNU_SOURCE) + * 所以这里需要使用自定义的 getline 函数,参数与 POSIX 定义相同 + */ +ssize_t grt_getline(char **lineptr, size_t *n, FILE *stream); + + /** * 处理单个震中距对应的数据逆变换和SAC保存 * diff --git a/pygrt/C_extension/src/common/model.c b/pygrt/C_extension/src/common/model.c index 9feadb17..e36a0567 100755 --- a/pygrt/C_extension/src/common/model.c +++ b/pygrt/C_extension/src/common/model.c @@ -211,7 +211,6 @@ GRT_MODEL1D * grt_read_mod1d_from_file(const char *command, const char *modelpat const int ncols = 6; // 模型文件有6列,或除去qa qb有四列 const int ncols_noQ = 4; - char line[1024]; int iline = 0; double h, va, vb, rho, qa, qb; double (*modarr)[ncols] = NULL; @@ -219,7 +218,10 @@ GRT_MODEL1D * grt_read_mod1d_from_file(const char *command, const char *modelpat int nlay = 0; mod1d->io_depth = false; - while(fgets(line, sizeof(line), fp)) { + size_t len; + char *line = NULL; + + while(grt_getline(&line, &len, fp) != -1) { iline++; // 注释行 @@ -372,6 +374,7 @@ GRT_MODEL1D * grt_read_mod1d_from_file(const char *command, const char *modelpat fclose(fp); GRT_SAFE_FREE_PTR(modarr); + GRT_SAFE_FREE_PTR(line); return mod1d; } @@ -379,13 +382,12 @@ GRT_MODEL1D * grt_read_mod1d_from_file(const char *command, const char *modelpat void grt_get_model_diglen_from_file(const char *command, const char *modelpath, MYINT diglen[6]){ FILE *fp = GRTCheckOpenFile(command, modelpath, "r"); - ssize_t read; size_t len; char *line = NULL; memset(diglen, 0, sizeof(MYINT)*6); - while((read = getline(&line, &len, fp) != -1)){ + while(grt_getline(&line, &len, fp) != -1){ char *token = strtok(line, " \n"); for(MYINT i=0; i<6; ++i){ if(token == NULL) break; diff --git a/pygrt/C_extension/src/common/util.c b/pygrt/C_extension/src/common/util.c index 4ab1863d..fadaa9e3 100644 --- a/pygrt/C_extension/src/common/util.c +++ b/pygrt/C_extension/src/common/util.c @@ -102,6 +102,62 @@ bool grt_is_comment_or_empty(const char* line) { } +ssize_t grt_getline(char **lineptr, size_t *n, FILE *stream){ + if (!lineptr || !n || !stream) { + return -1; + } + + char *buf = *lineptr; + size_t size = *n; + size_t len = 0; + int c; + + // 如果缓冲区为空,分配初始缓冲区 + if (buf == NULL || size == 0) { + size = 128; + buf = malloc(size); + if (buf == NULL) { + return -1; + } + } + + // 逐字符读取直到换行符或EOF + while ((c = fgetc(stream)) != EOF) { + // 检查是否需要扩展缓冲区 + if (len + 1 >= size) { + size_t new_size = size * 2; + char *new_buf = realloc(buf, new_size); + if (new_buf == NULL) { + free(buf); + return -1; + } + buf = new_buf; + size = new_size; + } + + buf[len++] = c; + + // 遇到换行符停止读取 + if (c == '\n') { + break; + } + } + + // 如果没有读取到任何字符且遇到EOF + if (len == 0 && c == EOF) { + return -1; + } + + // 添加字符串终止符 + buf[len] = '\0'; + + *lineptr = buf; + *n = size; + + return len; +} + + /** diff --git a/pygrt/C_extension/src/dynamic/signals.c b/pygrt/C_extension/src/dynamic/signals.c index 1edb9c68..b481d21f 100644 --- a/pygrt/C_extension/src/dynamic/signals.c +++ b/pygrt/C_extension/src/dynamic/signals.c @@ -324,9 +324,11 @@ float * grt_get_custom_wave(int *Nt, const char *tfparams){ } // 逐行读入 - char line[1024]; + size_t len; + char *line = NULL; + int nt = 0; - while(fgets(line, sizeof(line), fp)) { + while(grt_getline(&line, &len, fp) != -1) { // 注释行 if(grt_is_comment_or_empty(line)) continue; @@ -344,6 +346,7 @@ float * grt_get_custom_wave(int *Nt, const char *tfparams){ } fclose(fp); + GRT_SAFE_FREE_PTR(line); *Nt = nt; return tfarr; diff --git a/pygrt/C_extension/src/static/grt_static_rotation.c b/pygrt/C_extension/src/static/grt_static_rotation.c index 30a65b77..04e7822f 100644 --- a/pygrt/C_extension/src/static/grt_static_rotation.c +++ b/pygrt/C_extension/src/static/grt_static_rotation.c @@ -100,9 +100,11 @@ int static_rotation_main(int argc, char **argv){ bool rot2ZNE = false; // 逐行读入 - char line[1024]; + size_t len; + char *line = NULL; + int iline = 0; - while(fgets(line, sizeof(line), stdin)){ + while(grt_getline(&line, &len, stdin) != -1){ iline++; if(iline == 1){ // 读取震源物性参数 @@ -203,6 +205,7 @@ int static_rotation_main(int argc, char **argv){ GRTRaiseError("[%s] Error! Empty input. \n", command); } + GRT_SAFE_FREE_PTR(line); free_Ctrl(Ctrl); return EXIT_SUCCESS; } \ No newline at end of file diff --git a/pygrt/C_extension/src/static/grt_static_strain.c b/pygrt/C_extension/src/static/grt_static_strain.c index 94b67f3e..b3457691 100644 --- a/pygrt/C_extension/src/static/grt_static_strain.c +++ b/pygrt/C_extension/src/static/grt_static_strain.c @@ -98,9 +98,11 @@ int static_strain_main(int argc, char **argv){ bool rot2ZNE = false; // 逐行读入 - char line[1024]; + size_t len; + char *line = NULL; + int iline = 0; - while(fgets(line, sizeof(line), stdin)){ + while(grt_getline(&line, &len, stdin) != -1){ iline++; if(iline == 1){ // 读取震源物性参数 @@ -203,6 +205,7 @@ int static_strain_main(int argc, char **argv){ GRTRaiseError("[%s] Error! Empty input. \n", command); } + GRT_SAFE_FREE_PTR(line); free_Ctrl(Ctrl); return EXIT_SUCCESS; } \ No newline at end of file diff --git a/pygrt/C_extension/src/static/grt_static_stress.c b/pygrt/C_extension/src/static/grt_static_stress.c index 6fe285ef..e4753fa6 100644 --- a/pygrt/C_extension/src/static/grt_static_stress.c +++ b/pygrt/C_extension/src/static/grt_static_stress.c @@ -101,9 +101,11 @@ int static_stress_main(int argc, char **argv){ double lam_ukk=0.0; // 逐行读入 - char line[1024]; + size_t len; + char *line = NULL; + int iline = 0; - while(fgets(line, sizeof(line), stdin)){ + while(grt_getline(&line, &len, stdin) != -1){ iline++; if(iline == 1){ // 读取震源物性参数 @@ -221,6 +223,7 @@ int static_stress_main(int argc, char **argv){ GRTRaiseError("[%s] Error! Empty input. \n", command); } + GRT_SAFE_FREE_PTR(line); free_Ctrl(Ctrl); return EXIT_SUCCESS; } \ No newline at end of file diff --git a/pygrt/C_extension/src/static/grt_static_syn.c b/pygrt/C_extension/src/static/grt_static_syn.c index 81232f7c..50ba1d56 100644 --- a/pygrt/C_extension/src/static/grt_static_syn.c +++ b/pygrt/C_extension/src/static/grt_static_syn.c @@ -327,9 +327,12 @@ int static_syn_main(int argc, char **argv){ // 震中距 double dist = 0.0; - char line[1024]; + // 逐行读入 + size_t len; + char *line = NULL; + int iline = 0; - while(fgets(line, sizeof(line), stdin)){ + while(grt_getline(&line, &len, stdin) != -1){ iline++; if(iline == 1){ // 读取震源物性参数 @@ -505,6 +508,7 @@ int static_syn_main(int argc, char **argv){ GRTRaiseError("[%s] Error! Empty input. \n", command); } + GRT_SAFE_FREE_PTR(line); free_Ctrl(Ctrl); return EXIT_SUCCESS; } \ No newline at end of file From 3972d5474e4b7be69fc869981d2e92759a26661d Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Wed, 3 Sep 2025 16:06:55 +0800 Subject: [PATCH 5/5] FEAT: rearange fnctions in `search.h/c` and `matrix.h`, with X macros --- pygrt/C_extension/include/grt/common/matrix.h | 72 +++++++- pygrt/C_extension/include/grt/common/search.h | 87 +++++++-- pygrt/C_extension/src/common/search.c | 173 +++++++++++++----- pygrt/C_extension/src/dynamic/grn.c | 2 +- pygrt/C_extension/src/dynamic/grt_greenfn.c | 2 +- pygrt/C_extension/src/static/static_grn.c | 2 +- 6 files changed, 279 insertions(+), 59 deletions(-) diff --git a/pygrt/C_extension/include/grt/common/matrix.h b/pygrt/C_extension/include/grt/common/matrix.h index abc11b2b..2aa72827 100755 --- a/pygrt/C_extension/include/grt/common/matrix.h +++ b/pygrt/C_extension/include/grt/common/matrix.h @@ -3,7 +3,7 @@ * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn) * @date 2024-07-24 * - * 2x2小矩阵的加、减、乘、除、求逆等操作,由于均为小型数组操作,所有函数均为内联函数 + * 小矩阵的加、减、乘、除、求逆等操作,由于均为小型数组操作,所有函数均为内联函数 */ #pragma once @@ -142,7 +142,7 @@ inline GCC_ALWAYS_INLINE void grt_cmat2x2_assign(const MYCOMPLEX M1[2][2], MYCOM } /** - * 计算nxn复矩阵的积(小矩阵)(最暴力的方式) + * 计算mxn复矩阵的积(小矩阵)(最暴力的方式) * * @param[in] m1 M1矩阵行数 * @param[in] n1 M1矩阵列数 @@ -171,6 +171,31 @@ inline GCC_ALWAYS_INLINE void grt_cmatmxn_mul(MYINT m1, MYINT n1, MYINT p1, cons } } +/** + * 计算mxn复矩阵的转置矩阵(不共轭) + * + * @param[in] m1 M1矩阵行数 + * @param[in] n1 M1矩阵列数 + * @param[in] M1 M1矩阵 + * @param[out] M2 M2矩阵 (M1^T) + */ +inline GCC_ALWAYS_INLINE void grt_cmatmxn_transpose(MYINT m1, MYINT n1, const MYCOMPLEX M1[m1][n1], MYCOMPLEX M2[n1][m1]){ + MYINT m, n; + MYCOMPLEX M0[n1][m1]; + for(m=0; m #include "grt/common/const.h" +// 定义 X 宏,为多个类型定义查找函数 +#define __FOR_EACH_TYPE \ + X(MYINT) X(MYREAL) X(float) X(double) + + /** - * 该函数对输入的整数数组进行线性搜索,找到目标值时返回其索引。 + * 该函数对输入数组进行线性搜索,找到目标值时返回其索引。 * 如果目标值在数组中未找到,则返回 -1。 * - * @param[in] array 要搜索的整数数组。 + * @param[in] array 要搜索的数组。 * @param[in] size 数组的大小(元素个数)。 * @param[in] target 要查找的目标值。 * @@ -24,13 +29,19 @@ * @note 如果数组中存在多个目标值,该函数返回第一个匹配的索引。 * */ -MYINT grt_findElement_MYINT(const MYINT array[], MYINT size, MYINT target); +#define X(T) \ +MYINT grt_findElement_##T(const T *array, MYINT size, T target); + +__FOR_EACH_TYPE +#undef X + + /** - * 搜索浮点数数组中最接近目标值且小于目标值的索引。 + * 搜索数组中最接近目标值且小于目标值的索引。 * 如果目标值在数组中未找到,则返回 -1。 * - * @param[in] array 要搜索的浮点数数组。 + * @param[in] array 要搜索的数组。 * @param[in] size 数组的大小(元素个数)。 * @param[in] target 要查找的目标值。 * @@ -39,12 +50,18 @@ MYINT grt_findElement_MYINT(const MYINT array[], MYINT size, MYINT target); * @note 如果数组中存在多个目标值,该函数返回第一个匹配的索引。 * */ -MYINT grt_findLessEqualClosest_MYREAL(const MYREAL array[], MYINT size, MYREAL target); +#define X(T) \ +MYINT grt_findLessEqualClosest_##T(const T *array, MYINT size, T target); + +__FOR_EACH_TYPE +#undef X + + /** - * 搜索浮点数数组中最接近目标值的索引。 + * 搜索数组中最接近目标值的索引。 * - * @param[in] array 要搜索的浮点数数组。 + * @param[in] array 要搜索的数组。 * @param[in] size 数组的大小(元素个数)。 * @param[in] target 要查找的目标值。 * @@ -53,18 +70,62 @@ MYINT grt_findLessEqualClosest_MYREAL(const MYREAL array[], MYINT size, MYREAL t * @note 如果数组中存在多个目标值,该函数返回第一个匹配的索引。 * */ -MYINT grt_findClosest_MYREAL(const MYREAL array[], MYINT size, MYREAL target); +#define X(T) \ +MYINT grt_findClosest_##T(const T *array, MYINT size, T target); + +__FOR_EACH_TYPE +#undef X + /** - * 搜索浮点数数组的最大或最小值,返回其索引。 + * 搜索数组的最值,返回其索引。 * - * @param[in] array 要搜索的浮点数数组。 + * @param[in] array 要搜索的数组。 * @param[in] size 数组的大小(元素个数)。 - * @param[in] isMax 是否要找最大值,否则找最小值。 * * @return idx 目标值的索引。 * * @note 如果数组中存在相同最值,该函数返回第一个匹配的索引。 * */ -MYINT grt_findMinMax_MYREAL(const MYREAL array[], MYINT size, bool isMax); \ No newline at end of file +#define X(T) \ +MYINT grt_findMin_##T(const T *array, MYINT size);\ +MYINT grt_findMax_##T(const T *array, MYINT size);\ + +__FOR_EACH_TYPE +#undef X + + +/** + * 比较函数 + * + * @param[in] a 元素 a 地址 + * @param[in] b 元素 b 地址 + * + * @return flag 比较结果,(1) a > b, (0) a == b, (-1) a < b + */ +#define X(T) \ +int grt_compare_##T(const void *a, const void *b); + +__FOR_EACH_TYPE +#undef X +#undef __FOR_EACH_TYPE + + + +/** + * 在有序数组中插入元素,元素类型和数组类型需匹配 + * + * @param[in,out] arr 有序数组地址 + * @param[in,out] size 数组大小地址 + * @param[in] capacity 数组最大容量 + * @param[in] target 元素地址 + * @param[in] elementSize 元素和数组内元素的字节长度 + * @param[in] ascending 升序(true) 或 降序(false) + * @param[in] compare 比较函数 + * + * @return pos 插入位置的索引 + */ +MYINT grt_insertOrdered( + void *arr, MYINT *size, MYINT capacity, const void *target, size_t elementSize, bool ascending, + int (*compare)(const void *, const void *)); diff --git a/pygrt/C_extension/src/common/search.c b/pygrt/C_extension/src/common/search.c index 21b6a7d8..dfa2047b 100644 --- a/pygrt/C_extension/src/common/search.c +++ b/pygrt/C_extension/src/common/search.c @@ -8,59 +8,150 @@ #include #include +#include #include "grt/common/search.h" #include "grt/common/const.h" -static bool _gt_(MYREAL a1, MYREAL a2) { return a1 > a2; } -static bool _lt_(MYREAL a1, MYREAL a2) { return a1 < a2; } +// 定义 X 宏,为多个类型定义查找函数 +#define __FOR_EACH_TYPE \ + X(MYINT) X(MYREAL) X(float) X(double) -MYINT grt_findElement_MYINT(const MYINT array[], MYINT size, MYINT target) { - for (MYINT i = 0; i < size; ++i) { - if (array[i] == target) { - return i; // 找到目标元素,返回索引 - } - } - return -1; // 未找到目标元素,返回-1 +#define X(T) \ +MYINT grt_findElement_##T(const T *array, MYINT size, T target) {\ + for (MYINT i = 0; i < size; ++i) {\ + if (array[i] == target) {\ + /** 找到目标元素,返回索引 */\ + return i;\ + }\ + }\ + /** 未找到目标元素,返回-1 */\ + return -1; \ } -MYINT grt_findLessEqualClosest_MYREAL(const MYREAL array[], MYINT size, MYREAL target) { - MYINT ires=-1; - MYREAL mindist=-1.0, dist=0.0; - for (MYINT i = 0; i < size; ++i) { - dist = target-array[i]; - if(dist >= 0.0 && (mindist < 0.0 || dist < mindist)){ - ires = i; - mindist = dist; - } - } - return ires; +__FOR_EACH_TYPE +#undef X + + + +#define X(T) \ +MYINT grt_findLessEqualClosest_##T(const T *array, MYINT size, T target) {\ + MYINT ires=-1;\ + T mindist=-1, dist=0;\ + for (MYINT i = 0; i < size; ++i) {\ + dist = target-array[i];\ + if(dist >= 0 && (mindist < 0 || dist < mindist)){\ + ires = i;\ + mindist = dist;\ + }\ + }\ + return ires;\ } -MYINT grt_findClosest_MYREAL(const MYREAL array[], MYINT size, MYREAL target) { - MYINT ires=0; - MYREAL mindist=-1.0, dist=0.0; - for (MYINT i = 0; i < size; ++i) { - dist = fabs(target-array[i]); - if(mindist < 0.0 || dist < mindist){ - ires = i; - mindist = dist; - } - } - return ires; +__FOR_EACH_TYPE +#undef X + + + +#define X(T) \ +MYINT grt_findClosest_##T(const T *array, MYINT size, T target) {\ + MYINT ires=0;\ + T mindist=-1, dist=0;\ + for (MYINT i = 0; i < size; ++i) {\ + dist = fabs(target-array[i]);\ + if(mindist < 0 || dist < mindist){\ + ires = i;\ + mindist = dist;\ + }\ + }\ + return ires;\ +} + +__FOR_EACH_TYPE +#undef X + + +#define X(T) \ +MYINT grt_findMin_##T(const T *array, MYINT size) {\ + T rval = array[0];\ + MYINT idx=0;\ + for(MYINT ir=0; ir rval){\ + rval = array[ir];\ + idx = ir;\ + }\ + }\ + return idx;\ } -MYINT grt_findMinMax_MYREAL(const MYREAL array[], MYINT size, bool isMax) { - MYREAL rmax = array[0]; - MYINT idx=0; - bool (*_func)(MYREAL, MYREAL) = (isMax)? _gt_ : _lt_; - for(MYINT ir=0; ir valb){\ + return 1;\ + } else if (vala < valb){\ + return -1;\ + } else {\ + return 0;\ + }\ +}\ + +__FOR_EACH_TYPE +#undef X +#undef __FOR_EACH_TYPE + + +MYINT grt_insertOrdered( + void *arr, MYINT *size, MYINT capacity, const void *target, size_t elementSize, bool ascending, + int (*compare)(const void *, const void *)) +{ + MYINT sgn = (ascending)? 1 : -1; + + // 数组满载情况下,只可能插入更小(升序)或更大(降序)的数值 + if(*size == capacity && sgn*compare(target, arr+(*size-1)*elementSize) >= 0) return -1; + + // 找到插入位置 + MYINT pos=*size; + for(MYINT i=0; i<*size; ++i){ + if(sgn*compare(target, arr+i*elementSize) < 0){ + pos = i; + break; } } - return idx; -} + // 截断式插入,防止越界 + MYINT lastpos = *size; + if(lastpos >= capacity){ + lastpos = capacity-1; + } else { + ++(*size); + } + pos = GRT_MIN(pos, lastpos); + + // 移动插入位置后的元素 + memmove(arr + (pos + 1) * elementSize, + arr + pos * elementSize, + (lastpos - pos) * elementSize); + + // 插入新元素 + memcpy(arr + pos * elementSize, target, elementSize); + + return pos; +} \ No newline at end of file diff --git a/pygrt/C_extension/src/dynamic/grn.c b/pygrt/C_extension/src/dynamic/grn.c index 303fe509..6f4a4a3e 100755 --- a/pygrt/C_extension/src/dynamic/grn.c +++ b/pygrt/C_extension/src/dynamic/grn.c @@ -92,7 +92,7 @@ void grt_integ_grn_spec( gettimeofday(&begin_t, NULL); // 最大震中距 - MYINT irmax = grt_findMinMax_MYREAL(rs, nr, true); + MYINT irmax = grt_findMax_MYREAL(rs, nr); MYREAL rmax=rs[irmax]; const MYREAL Rho = mod1d->Rho[mod1d->isrc]; // 震源区密度 diff --git a/pygrt/C_extension/src/dynamic/grt_greenfn.c b/pygrt/C_extension/src/dynamic/grt_greenfn.c index 7fcef708..03a7890b 100644 --- a/pygrt/C_extension/src/dynamic/grt_greenfn.c +++ b/pygrt/C_extension/src/dynamic/grt_greenfn.c @@ -733,7 +733,7 @@ int greenfn_main(int argc, char **argv) { Ctrl->N.winT = Ctrl->N.nt*Ctrl->N.dt; // 最大震中距 - MYREAL rmax = Ctrl->R.rs[grt_findMinMax_MYREAL(Ctrl->R.rs, Ctrl->R.nr, true)]; + MYREAL rmax = Ctrl->R.rs[grt_findMax_MYREAL(Ctrl->R.rs, Ctrl->R.nr)]; // 时窗最大截止时刻 MYREAL tmax = Ctrl->E.delayT0 + Ctrl->N.winT; diff --git a/pygrt/C_extension/src/static/static_grn.c b/pygrt/C_extension/src/static/static_grn.c index ca7f2f1b..4be6c325 100644 --- a/pygrt/C_extension/src/static/static_grn.c +++ b/pygrt/C_extension/src/static/static_grn.c @@ -76,7 +76,7 @@ void grt_integ_static_grn( const char *statsstr // 积分结果输出 ){ - MYREAL rmax=rs[grt_findMinMax_MYREAL(rs, nr, true)]; // 最大震中距 + MYREAL rmax=rs[grt_findMax_MYREAL(rs, nr)]; // 最大震中距 const MYREAL hs = GRT_MAX(fabs(mod1d->depsrc - mod1d->deprcv), GRT_MIN_DEPTH_GAP_SRC_RCV); // hs=max(震源和台站深度差,1.0)