From 0bc151ca2332fcc00143d35f5ab584c15f84c8a2 Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Tue, 31 Mar 2026 19:08:28 +0800 Subject: [PATCH 1/3] rename GRT_MODEL1D to MODEL1D --- pygrt/C_extension/include/grt/common/model.h | 50 ++++++++--------- pygrt/C_extension/include/grt/dynamic/grn.h | 4 +- pygrt/C_extension/include/grt/dynamic/layer.h | 54 +++++++++---------- .../C_extension/include/grt/dynamic/source.h | 6 +-- pygrt/C_extension/include/grt/integral/dwm.h | 4 +- pygrt/C_extension/include/grt/integral/fim.h | 4 +- .../include/grt/integral/integ_method.h | 4 +- .../C_extension/include/grt/integral/kernel.h | 14 ++--- pygrt/C_extension/include/grt/integral/ptam.h | 4 +- .../C_extension/include/grt/integral/safim.h | 4 +- .../include/grt/static/static_grn.h | 4 +- .../include/grt/static/static_layer.h | 22 ++++---- .../include/grt/static/static_source.h | 6 +-- pygrt/C_extension/src/common/model.c | 30 +++++------ pygrt/C_extension/src/dynamic/grn.c | 4 +- pygrt/C_extension/src/dynamic/grt_greenfn.c | 4 +- pygrt/C_extension/src/dynamic/grt_kernel.c | 6 +-- pygrt/C_extension/src/dynamic/layer.c | 54 +++++++++---------- pygrt/C_extension/src/dynamic/source.c | 6 +-- pygrt/C_extension/src/integral/dwm.c | 2 +- pygrt/C_extension/src/integral/fim.c | 2 +- pygrt/C_extension/src/integral/integ_method.c | 2 +- pygrt/C_extension/src/integral/ptam.c | 2 +- pygrt/C_extension/src/integral/safim.c | 2 +- .../src/static/grt_static_greenfn.c | 4 +- pygrt/C_extension/src/static/static_grn.c | 2 +- pygrt/C_extension/src/static/static_layer.c | 22 ++++---- pygrt/C_extension/src/static/static_source.c | 6 +-- pygrt/C_extension/src/tools/grt_travt.c | 4 +- 29 files changed, 166 insertions(+), 166 deletions(-) diff --git a/pygrt/C_extension/include/grt/common/model.h b/pygrt/C_extension/include/grt/common/model.h index 2b16abfe..fd801d52 100755 --- a/pygrt/C_extension/include/grt/common/model.h +++ b/pygrt/C_extension/include/grt/common/model.h @@ -3,7 +3,7 @@ * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn) * @date 2024-07-24 * - * `GRT_MODEL1D` 结构体相关操作函数 + * `MODEL1D` 结构体相关操作函数 */ #pragma once @@ -75,42 +75,42 @@ typedef struct { cplx_t uiz_R_EV[2][2]; cplx_t uiz_R_EVL; -} GRT_MODEL1D; +} MODEL1D; /** - * 打印 GRT_MODEL1D 模型参数信息，主要用于调试程序 + * 打印 MODEL1D 模型参数信息，主要用于调试程序 * - * @param[in] mod1d `GRT_MODEL1D` 结构体指针 + * @param[in] mod1d `MODEL1D` 结构体指针 * */ -void grt_print_mod1d(const GRT_MODEL1D *mod1d); +void grt_print_mod1d(const MODEL1D *mod1d); /** - * 释放 `GRT_MODEL1D` 结构体指针 + * 释放 `MODEL1D` 结构体指针 * - * @param[out] mod1d `GRT_MODEL1D` 结构体指针 + * @param[out] mod1d `MODEL1D` 结构体指针 */ -void grt_free_mod1d(GRT_MODEL1D *mod1d); +void grt_free_mod1d(MODEL1D *mod1d); /** - * 初始化 GRT_MODEL1D 模型内存空间 + * 初始化 MODEL1D 模型内存空间 * * @param[in] n 模型层数 * - * @return `GRT_MODEL1D` 结构体指针 + * @return `MODEL1D` 结构体指针 * */ -GRT_MODEL1D * grt_init_mod1d(size_t n); +MODEL1D * grt_init_mod1d(size_t n); /** - * 复制 `GRT_MODEL1D` 结构体 + * 复制 `MODEL1D` 结构体 * - * @param[in] mod1d1 `GRT_MODEL1D` 源结构体指针 - * @return 复制好的 `GRT_MODEL1D` 结构体指针 + * @param[in] mod1d1 `MODEL1D` 源结构体指针 + * @return 复制好的 `MODEL1D` 结构体指针 * */ -GRT_MODEL1D * grt_copy_mod1d(const GRT_MODEL1D *mod1d1); +MODEL1D * grt_copy_mod1d(const MODEL1D *mod1d1); /** * 根据不同的 omega， 计算衰减系数，更新弹性模量 @@ -118,7 +118,7 @@ GRT_MODEL1D * grt_copy_mod1d(const GRT_MODEL1D *mod1d1); * @param[in,out] mod1d `MODEL1D` 结构体指针 * @param[in] omega 复数频率 */ -void grt_attenuate_mod1d(GRT_MODEL1D *mod1d, cplx_t omega); +void grt_attenuate_mod1d(MODEL1D *mod1d, cplx_t omega); /** * 根据记录好的圆频率和波数，计算相速度和每层的 xa, xb, caca, cbcb @@ -126,16 +126,16 @@ void grt_attenuate_mod1d(GRT_MODEL1D *mod1d, cplx_t omega); * @param[in,out] mod1d 模型结构体指针 * @param[in] k 波数 */ -void grt_mod1d_xa_xb(GRT_MODEL1D *mod1d, const real_t k); +void grt_mod1d_xa_xb(MODEL1D *mod1d, const real_t k); /** - * 扩容 `GRT_MODEL1D` 结构体 + * 扩容 `MODEL1D` 结构体 * * @param[in,out] mod1d `MODEL1D` 结构体指针 * @param[in] n 新层数 */ -void grt_realloc_mod1d(GRT_MODEL1D *mod1d, size_t n); +void grt_realloc_mod1d(MODEL1D *mod1d, size_t n); /** * 从文件中读取模型文件 @@ -145,10 +145,10 @@ void grt_realloc_mod1d(GRT_MODEL1D *mod1d, size_t n); * @param[in] deprcv 接收深度 * @param[in] allowLiquid 是否允许液体层 * - * @return `GRT_MODEL1D` 结构体指针 + * @return `MODEL1D` 结构体指针 * */ -GRT_MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, real_t deprcv, bool allowLiquid); +MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, real_t deprcv, bool allowLiquid); /** * 设置模型的边界条件，并对底界面做检查 @@ -157,7 +157,7 @@ GRT_MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, rea * @param[in] topbound 顶层边界条件 * @param[in] botbound 底层边界条件 */ -void grt_set_mod1d_boundary(GRT_MODEL1D *mod1d, GRT_BOUND_TYPE topbound, GRT_BOUND_TYPE botbound); +void grt_set_mod1d_boundary(MODEL1D *mod1d, GRT_BOUND_TYPE topbound, GRT_BOUND_TYPE botbound); /** * 从模型文件中判断各个量的大致精度（字符串长度），以确定浮点数输出位数 @@ -177,14 +177,14 @@ void grt_get_model_diglen_from_file(const char *modelpath, size_t diglen[6]); * * @return 是否存在 */ -bool grt_check_vel_in_mod(const GRT_MODEL1D *mod1d, const real_t vel, const real_t tol); +bool grt_check_vel_in_mod(const MODEL1D *mod1d, const real_t vel, const real_t tol); /** * 计算最大最小速度（非零值） * - * @param mod1d (in)`GRT_MODEL1D` 结构体指针 + * @param mod1d (in)`MODEL1D` 结构体指针 * @param vmin (out)最小速度 * @param vmax (out)最大速度 * */ -void grt_get_mod1d_vmin_vmax(const GRT_MODEL1D *mod1d, real_t *vmin, real_t *vmax); \ No newline at end of file +void grt_get_mod1d_vmin_vmax(const MODEL1D *mod1d, real_t *vmin, real_t *vmax); \ No newline at end of file diff --git a/pygrt/C_extension/include/grt/dynamic/grn.h b/pygrt/C_extension/include/grt/dynamic/grn.h index 53769c04..f6ded88c 100755 --- a/pygrt/C_extension/include/grt/dynamic/grn.h +++ b/pygrt/C_extension/include/grt/dynamic/grn.h @@ -20,7 +20,7 @@ /** * 积分计算Z, R, T三个分量格林函数的频谱的核心函数（被Python调用） * - * @param[in,out] mod1d `GRT_MODEL1D` 结构体指针 + * @param[in,out] mod1d `MODEL1D` 结构体指针 * @param[in] nf1 开始计算频谱的频率索引值, 总范围在[nf1, nf2] * @param[in] nf2 结束计算频谱的频率索引值, 总范围在[nf1, nf2] * @param[in] freqs 所有频点的频率值（包括未计算的） @@ -41,7 +41,7 @@ * */ void grt_integ_grn_spec( - GRT_MODEL1D *mod1d, size_t nf1, size_t nf2, real_t *freqs, + MODEL1D *mod1d, size_t nf1, size_t nf2, real_t *freqs, size_t nr, real_t *rs, real_t wI, bool keepAllFreq, K_INTEG_METHOD *Kmet, bool print_progressbar, bool calc_upar, pcplxChnlGrid grn[nr], diff --git a/pygrt/C_extension/include/grt/dynamic/layer.h b/pygrt/C_extension/include/grt/dynamic/layer.h index 69e91a5a..7ab5b685 100755 --- a/pygrt/C_extension/include/grt/dynamic/layer.h +++ b/pygrt/C_extension/include/grt/dynamic/layer.h @@ -24,8 +24,8 @@ * @param[in,out] mod1d 模型结构体指针，结果保存在 M_top * */ -void grt_topbound_RU_PSV(GRT_MODEL1D *mod1d); -void grt_topbound_RU_SH(GRT_MODEL1D *mod1d); +void grt_topbound_RU_PSV(MODEL1D *mod1d); +void grt_topbound_RU_SH(MODEL1D *mod1d); /** * 计算不同边界条件下底界面的反射系数RD，其中自由表面的公式见(5.3.10-14) @@ -34,8 +34,8 @@ void grt_topbound_RU_SH(GRT_MODEL1D *mod1d); * @param[in,out] mod1d 模型结构体指针，结果保存在 M_bot * */ -void grt_botbound_RD_PSV(GRT_MODEL1D *mod1d); -void grt_botbound_RD_SH(GRT_MODEL1D *mod1d); +void grt_botbound_RD_PSV(MODEL1D *mod1d); +void grt_botbound_RD_SH(MODEL1D *mod1d); /** @@ -44,7 +44,7 @@ void grt_botbound_RD_SH(GRT_MODEL1D *mod1d); * @param[in,out] mod1d 模型结构体指针，结果保存在 R_EV * */ -void grt_wave2qwv_REV_PSV(GRT_MODEL1D *mod1d); +void grt_wave2qwv_REV_PSV(MODEL1D *mod1d); /** * 计算接收点位置的 SH 波接收矩阵，将波场转为位移，公式(5.2.19) + (5.7.7,25) @@ -52,7 +52,7 @@ void grt_wave2qwv_REV_PSV(GRT_MODEL1D *mod1d); * @param[in,out] mod1d 模型结构体指针，结果保存在 R_EVL * */ -void grt_wave2qwv_REV_SH(GRT_MODEL1D *mod1d); +void grt_wave2qwv_REV_SH(MODEL1D *mod1d); /** @@ -62,7 +62,7 @@ void grt_wave2qwv_REV_SH(GRT_MODEL1D *mod1d); * @param[in,out] mod1d 模型结构体指针，结果保存在 uiz_R_EV * */ -void grt_wave2qwv_z_REV_PSV(GRT_MODEL1D *mod1d); +void grt_wave2qwv_z_REV_PSV(MODEL1D *mod1d); /** @@ -72,7 +72,7 @@ void grt_wave2qwv_z_REV_PSV(GRT_MODEL1D *mod1d); * @param[in,out] mod1d 模型结构体指针，结果保存在 uiz_R_EVL * */ -void grt_wave2qwv_z_REV_SH(GRT_MODEL1D *mod1d); +void grt_wave2qwv_z_REV_SH(MODEL1D *mod1d); /** @@ -86,7 +86,7 @@ void grt_wave2qwv_z_REV_SH(GRT_MODEL1D *mod1d); * @param[out] M R/T矩阵 * */ -void grt_RT_matrix_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); +void grt_RT_matrix_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); /** @@ -98,25 +98,25 @@ void grt_RT_matrix_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); * @param[out] M R/T矩阵 * */ -void grt_RT_matrix_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); +void grt_RT_matrix_SH(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); /** 液-液 界面 */ -void grt_RT_matrix_ll_PSV(const GRT_MODEL1D *mod1d, size_t iy, RT_MATRIX *M); +void grt_RT_matrix_ll_PSV(const MODEL1D *mod1d, size_t iy, RT_MATRIX *M); /** 液-液 界面 */ void grt_RT_matrix_ll_SH(RT_MATRIX *M); /** 液-固 界面 */ -void grt_RT_matrix_ls_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); +void grt_RT_matrix_ls_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); /** 液-固 界面 */ -void grt_RT_matrix_ls_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); +void grt_RT_matrix_ls_SH(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); /** 固-固 界面 */ -void grt_RT_matrix_ss_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); +void grt_RT_matrix_ss_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); /** 固-固 界面 */ -void grt_RT_matrix_ss_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); +void grt_RT_matrix_ss_SH(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); /** * 为 R/T 矩阵添加时间延迟因子 @@ -126,9 +126,9 @@ void grt_RT_matrix_ss_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M * @param[out] M R/T矩阵 * */ -void grt_delay_RT_matrix(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); -void grt_delay_RT_matrix_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); -void grt_delay_RT_matrix_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); +void grt_delay_RT_matrix(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); +void grt_delay_RT_matrix_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); +void grt_delay_RT_matrix_SH(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); /** * 为虚拟层的广义 R/T 矩阵添加时间延迟因子 @@ -138,7 +138,7 @@ void grt_delay_RT_matrix_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX * @param[out] M R/T矩阵 * */ -void grt_delay_GRT_matrix(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); +void grt_delay_GRT_matrix(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); /** * 计算该层的连接 P-SV 应力位移矢量与垂直波函数的D矩阵(或其逆矩阵)， @@ -155,25 +155,25 @@ void grt_delay_GRT_matrix(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX * * @param[out] D D矩阵(或其逆矩阵) * */ -void grt_get_layer_D(const GRT_MODEL1D *mod1d, const size_t iy, bool inverse, int liquid_invtype, cplx_t D[4][4]); +void grt_get_layer_D(const MODEL1D *mod1d, const size_t iy, bool inverse, int liquid_invtype, cplx_t D[4][4]); /** 子矩阵 D11，函数参数见 get_layer_D 函数 */ -void grt_get_layer_D11(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]); +void grt_get_layer_D11(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]); /** 子矩阵 D12，函数参数见 get_layer_D 函数 */ -void grt_get_layer_D12(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]); +void grt_get_layer_D12(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]); /** 子矩阵 D21，函数参数见 get_layer_D 函数 */ -void grt_get_layer_D21(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]); +void grt_get_layer_D21(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]); /** 子矩阵 D22，函数参数见 get_layer_D 函数 */ -void grt_get_layer_D22(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]); +void grt_get_layer_D22(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]); /** 子矩阵 D11_uiz，后缀uiz表示连接位移对z的偏导和垂直波函数，函数参数见 get_layer_D 函数 */ -void grt_get_layer_D11_uiz(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]); +void grt_get_layer_D11_uiz(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]); /** 子矩阵 D12_uiz，函数参数见 get_layer_D 函数 */ -void grt_get_layer_D12_uiz(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]); +void grt_get_layer_D12_uiz(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]); /** @@ -186,7 +186,7 @@ void grt_get_layer_D12_uiz(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2 * @param[out] T T矩阵(或其逆矩阵) * */ -void grt_get_layer_T(const GRT_MODEL1D *mod1d, const size_t iy, bool inverse, cplx_t T[2][2]); +void grt_get_layer_T(const MODEL1D *mod1d, const size_t iy, bool inverse, cplx_t T[2][2]); // /** diff --git a/pygrt/C_extension/include/grt/dynamic/source.h b/pygrt/C_extension/include/grt/dynamic/source.h index 2f9ad013..545f042e 100755 --- a/pygrt/C_extension/include/grt/dynamic/source.h +++ b/pygrt/C_extension/include/grt/dynamic/source.h @@ -24,10 +24,10 @@ * @param[out] src_coefU 上行震源系数 * */ -void grt_source_coef(const GRT_MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU); +void grt_source_coef(const MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU); /* P-SV 波的震源系数 */ -void grt_source_coef_PSV(const GRT_MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU); +void grt_source_coef_PSV(const MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU); /* SH 波的震源系数 */ -void grt_source_coef_SH(const GRT_MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU); +void grt_source_coef_SH(const MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU); diff --git a/pygrt/C_extension/include/grt/integral/dwm.h b/pygrt/C_extension/include/grt/integral/dwm.h index 56ab0f87..b1411c29 100644 --- a/pygrt/C_extension/include/grt/integral/dwm.h +++ b/pygrt/C_extension/include/grt/integral/dwm.h @@ -24,7 +24,7 @@ * 传统的离散波数积分，结果以三维数组的形式返回，形状分别代表震中距、不同震源不同阶数 * 和4种积分类型(p=0,1,2,3) * - * @param[in,out] mod1d `GRT_MODEL1D` 结构体指针 + * @param[in,out] mod1d `MODEL1D` 结构体指针 * @param[in] dk 波数积分间隔 * @param[in] kmax 波数积分的上限 * @param[in] keps 波数积分的收敛条件，要求在某震中距下所有格林函数都收敛，为负数代表不提前判断收敛，按照波数积分上限进行积分 @@ -39,5 +39,5 @@ * @return k 积分截至时的波数 */ real_t grt_discrete_integ( - GRT_MODEL1D *mod1d, real_t dk, real_t kmax, real_t keps, + MODEL1D *mod1d, real_t dk, real_t kmax, real_t keps, size_t nr, real_t *rs, K_INTEG *K, FILE *fstats, GRT_KernelFunc kerfunc); diff --git a/pygrt/C_extension/include/grt/integral/fim.h b/pygrt/C_extension/include/grt/integral/fim.h index 8791b41f..0440564b 100755 --- a/pygrt/C_extension/include/grt/integral/fim.h +++ b/pygrt/C_extension/include/grt/integral/fim.h @@ -28,7 +28,7 @@ * 其中\f$x=kr\f$. * * - * @param[in,out] mod1d `GRT_MODEL1D` 结构体指针 + * @param[in,out] mod1d `MODEL1D` 结构体指针 * @param[in] k0 前一部分的波数积分结束点k值 * @param[in] dk0 前一部分的波数积分间隔 * @param[in] filondk filon积分间隔 @@ -45,7 +45,7 @@ * @return k 积分截至时的波数 */ real_t grt_linear_filon_integ( - GRT_MODEL1D *mod1d, real_t k0, real_t dk0, real_t dk, real_t kmax, real_t keps, + MODEL1D *mod1d, real_t k0, real_t dk0, real_t dk, real_t kmax, real_t keps, size_t nr, real_t *rs, K_INTEG *K, FILE *fstats, GRT_KernelFunc kerfunc); diff --git a/pygrt/C_extension/include/grt/integral/integ_method.h b/pygrt/C_extension/include/grt/integral/integ_method.h index ed0f632c..cc26ffd6 100644 --- a/pygrt/C_extension/include/grt/integral/integ_method.h +++ b/pygrt/C_extension/include/grt/integral/integ_method.h @@ -74,7 +74,7 @@ void grt_KMET_destroy_fstats(const size_t nr, K_INTEG_METHOD *Kmet); /** * 发起波数积分的总函数 * - * @param[in,out] mod1d `GRT_MODEL1D` 结构体指针 + * @param[in,out] mod1d `MODEL1D` 结构体指针 * @param[in] nr 震中距数量 * @param[in] rs 震中距数组 * @param[in,out] Kmet K_INTEG_METHOD 结构体指针 @@ -85,4 +85,4 @@ void grt_KMET_destroy_fstats(const size_t nr, K_INTEG_METHOD *Kmet); * */ K_INTEG * grt_wavenumber_integral( - GRT_MODEL1D *mod1d, size_t nr, real_t *rs, K_INTEG_METHOD *Kmet, bool calc_upar, GRT_KernelFunc kerfunc); + MODEL1D *mod1d, size_t nr, real_t *rs, K_INTEG_METHOD *Kmet, bool calc_upar, GRT_KernelFunc kerfunc); diff --git a/pygrt/C_extension/include/grt/integral/kernel.h b/pygrt/C_extension/include/grt/integral/kernel.h index f9df7fa6..7419adc2 100644 --- a/pygrt/C_extension/include/grt/integral/kernel.h +++ b/pygrt/C_extension/include/grt/integral/kernel.h @@ -16,7 +16,7 @@ * 计算核函数的函数指针，动态与静态的接口一致 */ typedef void (*GRT_KernelFunc) ( - GRT_MODEL1D *mod1d, const real_t k, cplxChnlGrid QWV, + MODEL1D *mod1d, const real_t k, cplxChnlGrid QWV, bool calc_uiz, cplxChnlGrid QWV_uiz); @@ -89,26 +89,26 @@ typedef void (*GRT_KernelFunc) ( * */ void grt_kernel( - GRT_MODEL1D *mod1d, const real_t k, cplxChnlGrid QWV, + MODEL1D *mod1d, const real_t k, cplxChnlGrid QWV, bool calc_uiz, cplxChnlGrid QWV_uiz); /** 构建广义反射透射系数矩阵。作为 kernel 函数中的第一部分 */ -void grt_GRT_matrix(GRT_MODEL1D *mod1d, const real_t k); +void grt_GRT_matrix(MODEL1D *mod1d, const real_t k); /** 从广义 R/T 矩阵出发，计算每个震源对应的核函数 QWV。 作为 kernel 函数中的第二部分 */ void grt_GRT_build_QWV( - GRT_MODEL1D *mod1d, cplxChnlGrid QWV, + MODEL1D *mod1d, cplxChnlGrid QWV, bool calc_uiz, cplxChnlGrid QWV_uiz); /** 静态解的核函数 */ void grt_static_kernel( - GRT_MODEL1D *mod1d, const real_t k, cplxChnlGrid QWV, + MODEL1D *mod1d, const real_t k, cplxChnlGrid QWV, bool calc_uiz, cplxChnlGrid QWV_uiz); /** 静态广义反射透射系数矩阵 */ -void grt_static_GRT_matrix(GRT_MODEL1D *mod1d, const real_t k); +void grt_static_GRT_matrix(MODEL1D *mod1d, const real_t k); /** 静态 QWV */ void grt_static_GRT_build_QWV( - GRT_MODEL1D *mod1d, cplxChnlGrid QWV, + MODEL1D *mod1d, cplxChnlGrid QWV, bool calc_uiz, cplxChnlGrid QWV_uiz); \ No newline at end of file diff --git a/pygrt/C_extension/include/grt/integral/ptam.h b/pygrt/C_extension/include/grt/integral/ptam.h index c4418dc8..4b081daf 100755 --- a/pygrt/C_extension/include/grt/integral/ptam.h +++ b/pygrt/C_extension/include/grt/integral/ptam.h @@ -26,7 +26,7 @@ /** * 峰谷平均法 Peak-Trough Averaging Method，最后收敛的积分结果以三维数组的形式返回， * - * @param[in,out] mod1d `GRT_MODEL1D` 结构体指针 + * @param[in,out] mod1d `MODEL1D` 结构体指针 * @param[in] k0 先前的积分已经进行到了波数k0 * @param[in] predk 先前的积分使用的积分间隔dk，因为峰谷平均法使用的 * 积分间隔会和之前的不一致，这里传入该系数以做预先调整 @@ -41,6 +41,6 @@ * */ void grt_PTA_method( - GRT_MODEL1D *mod1d, real_t k0, real_t predk, + MODEL1D *mod1d, real_t k0, real_t predk, size_t nr, real_t *rs, K_INTEG *K, FILE *ptam_fstatsnr[nr][2], GRT_KernelFunc kerfunc); diff --git a/pygrt/C_extension/include/grt/integral/safim.h b/pygrt/C_extension/include/grt/integral/safim.h index af5ec56d..f71afafc 100755 --- a/pygrt/C_extension/include/grt/integral/safim.h +++ b/pygrt/C_extension/include/grt/integral/safim.h @@ -30,7 +30,7 @@ * 其中\f$x=kr\f$. * * - * @param[in,out] mod1d `GRT_MODEL1D` 结构体指针 + * @param[in,out] mod1d `MODEL1D` 结构体指针 * @param[in] k0 前一部分的波数积分结束点k值 * @param[in] dk0 前一部分的波数积分间隔 * @param[in] tol 自适应Filon积分的采样精度 @@ -47,7 +47,7 @@ * @return k 积分截至时的波数 */ real_t grt_sa_filon_integ( - GRT_MODEL1D *mod1d, real_t k0, real_t dk0, real_t tol, real_t kmax, real_t kref, + MODEL1D *mod1d, real_t k0, real_t dk0, real_t tol, real_t kmax, real_t kref, size_t nr, real_t *rs, K_INTEG *K, FILE *fstats, GRT_KernelFunc kerfunc); diff --git a/pygrt/C_extension/include/grt/static/static_grn.h b/pygrt/C_extension/include/grt/static/static_grn.h index 9ec6d882..96514c76 100644 --- a/pygrt/C_extension/include/grt/static/static_grn.h +++ b/pygrt/C_extension/include/grt/static/static_grn.h @@ -23,7 +23,7 @@ /** * 积分计算Z, R, T三个分量静态格林函数的核心函数 * - * @param[in,out] mod1d `GRT_MODEL1D` 结构体指针 + * @param[in,out] mod1d `MODEL1D` 结构体指针 * @param[in] nr 震中距数量 * @param[in] rs 震中距数组 * @param[in,out] Kmet 波数积分相关参数的结构体指针 @@ -36,7 +36,7 @@ * */ void grt_integ_static_grn( - GRT_MODEL1D *mod1d, size_t nr, real_t *rs, K_INTEG_METHOD *Kmet, + MODEL1D *mod1d, size_t nr, real_t *rs, K_INTEG_METHOD *Kmet, bool calc_upar, realChnlGrid grn[nr], realChnlGrid grn_uiz[nr], diff --git a/pygrt/C_extension/include/grt/static/static_layer.h b/pygrt/C_extension/include/grt/static/static_layer.h index b1479e5f..e1a18eea 100644 --- a/pygrt/C_extension/include/grt/static/static_layer.h +++ b/pygrt/C_extension/include/grt/static/static_layer.h @@ -25,8 +25,8 @@ * @param[in,out] mod1d 模型结构体指针，结果保存在 M_top * */ -void grt_static_topbound_RU_PSV(GRT_MODEL1D *mod1d); -void grt_static_topbound_RU_SH(GRT_MODEL1D *mod1d); +void grt_static_topbound_RU_PSV(MODEL1D *mod1d); +void grt_static_topbound_RU_SH(MODEL1D *mod1d); /** * 计算不同边界条件下底界面的静态反射系数RD，其中自由表面的公式见(6.3.12) @@ -34,8 +34,8 @@ void grt_static_topbound_RU_SH(GRT_MODEL1D *mod1d); * @param[in,out] mod1d 模型结构体指针，结果保存在 M_bot * */ -void grt_static_botbound_RD_PSV(GRT_MODEL1D *mod1d); -void grt_static_botbound_RD_SH(GRT_MODEL1D *mod1d); +void grt_static_botbound_RD_PSV(MODEL1D *mod1d); +void grt_static_botbound_RD_SH(MODEL1D *mod1d); /** * 计算接收点位置的 P-SV 静态接收矩阵，将波场转为位移，公式(6.3.35) @@ -43,7 +43,7 @@ void grt_static_botbound_RD_SH(GRT_MODEL1D *mod1d); * @param[in,out] mod1d 模型结构体指针，结果保存在 R_EV * */ -void grt_static_wave2qwv_REV_PSV(GRT_MODEL1D *mod1d); +void grt_static_wave2qwv_REV_PSV(MODEL1D *mod1d); /** * 计算接收点位置的 SH 静态接收矩阵，将波场转为位移，公式(6.3.37) @@ -51,7 +51,7 @@ void grt_static_wave2qwv_REV_PSV(GRT_MODEL1D *mod1d); * @param[in,out] mod1d 模型结构体指针，结果保存在 R_EVL * */ -void grt_static_wave2qwv_REV_SH(GRT_MODEL1D *mod1d); +void grt_static_wave2qwv_REV_SH(MODEL1D *mod1d); /** * 计算接收点位置的ui_z的 P-SV 静态接收矩阵，即将波场转为ui_z。 @@ -60,7 +60,7 @@ void grt_static_wave2qwv_REV_SH(GRT_MODEL1D *mod1d); * @param[in,out] mod1d 模型结构体指针，结果保存在 uiz_R_EV * */ -void grt_static_wave2qwv_z_REV_PSV(GRT_MODEL1D *mod1d); +void grt_static_wave2qwv_z_REV_PSV(MODEL1D *mod1d); /** * 计算接收点位置的ui_z的 SH 静态接收矩阵，即将波场转为ui_z。 @@ -69,7 +69,7 @@ void grt_static_wave2qwv_z_REV_PSV(GRT_MODEL1D *mod1d); * @param[in,out] mod1d 模型结构体指针，结果保存在 uiz_R_EVL * */ -void grt_static_wave2qwv_z_REV_SH(GRT_MODEL1D *mod1d); +void grt_static_wave2qwv_z_REV_SH(MODEL1D *mod1d); /** @@ -81,7 +81,7 @@ void grt_static_wave2qwv_z_REV_SH(GRT_MODEL1D *mod1d); * @param[out] M R/T矩阵 * */ -void grt_static_RT_matrix_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); +void grt_static_RT_matrix_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); /** * 计算界面的 SH 波静态反射透射系数RDL/RUL/TDL/TUL @@ -92,7 +92,7 @@ void grt_static_RT_matrix_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATR * @param[out] M R/T矩阵 * */ -void grt_static_RT_matrix_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); +void grt_static_RT_matrix_SH(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); /** @@ -103,4 +103,4 @@ void grt_static_RT_matrix_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRI * @param[out] M R/T矩阵 * */ -void grt_static_delay_RT_matrix(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); \ No newline at end of file +void grt_static_delay_RT_matrix(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); \ No newline at end of file diff --git a/pygrt/C_extension/include/grt/static/static_source.h b/pygrt/C_extension/include/grt/static/static_source.h index 9119b5ef..821fc23f 100644 --- a/pygrt/C_extension/include/grt/static/static_source.h +++ b/pygrt/C_extension/include/grt/static/static_source.h @@ -23,10 +23,10 @@ * @param[out] src_coefU 上行震源系数 * */ -void grt_static_source_coef(const GRT_MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU); +void grt_static_source_coef(const MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU); /* P-SV 波的静态震源系数 */ -void grt_static_source_coef_PSV(const GRT_MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU); +void grt_static_source_coef_PSV(const MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU); /* SH 波的静态震源系数 */ -void grt_static_source_coef_SH(const GRT_MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU); +void grt_static_source_coef_SH(const MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU); diff --git a/pygrt/C_extension/src/common/model.c b/pygrt/C_extension/src/common/model.c index 043a897b..c128cdf8 100755 --- a/pygrt/C_extension/src/common/model.c +++ b/pygrt/C_extension/src/common/model.c @@ -3,7 +3,7 @@ * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn) * @date 2024-07-24 * - * GRT_MODEL1D 结构体的相关操作函数 + * MODEL1D 结构体的相关操作函数 * */ @@ -43,7 +43,7 @@ X(cbcb, cplx_t)\ -void grt_print_mod1d(const GRT_MODEL1D *mod1d){ +void grt_print_mod1d(const MODEL1D *mod1d){ // 模拟表格，打印速度 // 每列字符宽度 // [isrc/ircv] [h(km)] [Vp(km/s)] [Vs(km/s)] [Rho(g/cm^3)] [Qp] [Qs] @@ -110,7 +110,7 @@ void grt_print_mod1d(const GRT_MODEL1D *mod1d){ printf("\n"); } -void grt_free_mod1d(GRT_MODEL1D *mod1d){ +void grt_free_mod1d(MODEL1D *mod1d){ #define X(P, T) GRT_SAFE_FREE_PTR(mod1d->P); GRT_FOR_EACH_MODEL_QUANTITY_ARRAY #undef X @@ -119,8 +119,8 @@ void grt_free_mod1d(GRT_MODEL1D *mod1d){ } -GRT_MODEL1D * grt_init_mod1d(size_t n){ - GRT_MODEL1D *mod1d = (GRT_MODEL1D *)calloc(1, sizeof(GRT_MODEL1D)); +MODEL1D * grt_init_mod1d(size_t n){ + MODEL1D *mod1d = (MODEL1D *)calloc(1, sizeof(MODEL1D)); mod1d->n = n; mod1d->omgref = PI2*1.0; @@ -132,8 +132,8 @@ GRT_MODEL1D * grt_init_mod1d(size_t n){ } -GRT_MODEL1D * grt_copy_mod1d(const GRT_MODEL1D *mod1d1){ - GRT_MODEL1D *mod1d2 = (GRT_MODEL1D *)calloc(1, sizeof(GRT_MODEL1D)); +MODEL1D * grt_copy_mod1d(const MODEL1D *mod1d1){ + MODEL1D *mod1d2 = (MODEL1D *)calloc(1, sizeof(MODEL1D)); // 先直接赋值，实现浅拷贝 *mod1d2 = *mod1d1; @@ -151,7 +151,7 @@ GRT_MODEL1D * grt_copy_mod1d(const GRT_MODEL1D *mod1d1){ } -void grt_attenuate_mod1d(GRT_MODEL1D *mod1d, cplx_t omega){ +void grt_attenuate_mod1d(MODEL1D *mod1d, cplx_t omega){ real_t Va0, Vb0; cplx_t atna, atnb; for(size_t i=0; in; ++i){ @@ -177,7 +177,7 @@ void grt_attenuate_mod1d(GRT_MODEL1D *mod1d, cplx_t omega){ } -void grt_mod1d_xa_xb(GRT_MODEL1D *mod1d, const real_t k) +void grt_mod1d_xa_xb(MODEL1D *mod1d, const real_t k) { mod1d->k = k; // 不合理的频率值，只可能是在计算静态解，此时不需要xa, xb等物理量 @@ -218,7 +218,7 @@ void grt_mod1d_xa_xb(GRT_MODEL1D *mod1d, const real_t k) } -void grt_realloc_mod1d(GRT_MODEL1D *mod1d, size_t n){ +void grt_realloc_mod1d(MODEL1D *mod1d, size_t n){ mod1d->n = n; #define X(P, T) mod1d->P = (T*)realloc(mod1d->P, n*sizeof(T)); @@ -228,7 +228,7 @@ void grt_realloc_mod1d(GRT_MODEL1D *mod1d, size_t n){ -GRT_MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, real_t deprcv, bool allowLiquid){ +MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, real_t deprcv, bool allowLiquid){ GRTCheckFileExist(modelpath); if(depsrc * deprcv < 0.0){ @@ -239,7 +239,7 @@ GRT_MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, rea // 初始化 - GRT_MODEL1D *mod1d = grt_init_mod1d(1); + MODEL1D *mod1d = grt_init_mod1d(1); const int ncols = 6; // 模型文件有6列，或除去qa qb有四列 const int ncols_noQ = 4; @@ -441,7 +441,7 @@ GRT_MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, rea } -void grt_set_mod1d_boundary(GRT_MODEL1D *mod1d, GRT_BOUND_TYPE topbound, GRT_BOUND_TYPE botbound) +void grt_set_mod1d_boundary(MODEL1D *mod1d, GRT_BOUND_TYPE topbound, GRT_BOUND_TYPE botbound) { mod1d->topbound = topbound; mod1d->botbound = botbound; @@ -478,7 +478,7 @@ void grt_get_model_diglen_from_file(const char *modelpath, size_t diglen[6]){ } -bool grt_check_vel_in_mod(const GRT_MODEL1D *mod1d, const real_t vel, const real_t tol){ +bool grt_check_vel_in_mod(const MODEL1D *mod1d, const real_t vel, const real_t tol){ // 浮点数比较，检查是否存在该速度值 for(size_t i=0; in; ++i){ if(fabs(vel - mod1d->Va[i])Vb[i])Va; diff --git a/pygrt/C_extension/src/dynamic/grn.c b/pygrt/C_extension/src/dynamic/grn.c index 93e5d294..945359fe 100755 --- a/pygrt/C_extension/src/dynamic/grn.c +++ b/pygrt/C_extension/src/dynamic/grn.c @@ -62,7 +62,7 @@ static void recordin_GRN( void grt_integ_grn_spec( - GRT_MODEL1D *mod1d, size_t nf1, size_t nf2, real_t *freqs, + MODEL1D *mod1d, size_t nf1, size_t nf2, real_t *freqs, size_t nr, real_t *rs, real_t wI, bool keepAllFreq, K_INTEG_METHOD *Kmet, bool print_progressbar, bool calc_upar, pcplxChnlGrid grn[nr], @@ -118,7 +118,7 @@ void grt_integ_grn_spec( cplx_t coef = - dk*fac / GRT_SQUARE(omega); // 最终要乘上的系数 - GRT_MODEL1D *local_mod1d = NULL; + MODEL1D *local_mod1d = NULL; K_INTEG_METHOD *local_Kmet = NULL; #ifdef _OPENMP // 定义局部模型对象 diff --git a/pygrt/C_extension/src/dynamic/grt_greenfn.c b/pygrt/C_extension/src/dynamic/grt_greenfn.c index e3f80731..5ad070ef 100644 --- a/pygrt/C_extension/src/dynamic/grt_greenfn.c +++ b/pygrt/C_extension/src/dynamic/grt_greenfn.c @@ -44,7 +44,7 @@ typedef struct { bool active; char *s_modelpath; ///< 模型路径 const char *s_modelname; ///< 模型名称 - GRT_MODEL1D *mod1d; ///< 模型结构体指针 + MODEL1D *mod1d; ///< 模型结构体指针 } M; /** 震源和接收器深度 */ struct { @@ -978,7 +978,7 @@ int greenfn_main(int argc, char **argv) { if((Ctrl->M.mod1d = grt_read_mod1d_from_file(Ctrl->M.s_modelpath, Ctrl->D.depsrc, Ctrl->D.deprcv, true)) == NULL){ exit(EXIT_FAILURE); } - GRT_MODEL1D *mod1d = Ctrl->M.mod1d; + MODEL1D *mod1d = Ctrl->M.mod1d; // 边界条件 grt_set_mod1d_boundary(mod1d, Ctrl->B.topbound, Ctrl->B.botbound); diff --git a/pygrt/C_extension/src/dynamic/grt_kernel.c b/pygrt/C_extension/src/dynamic/grt_kernel.c index 87942a1a..75f4aa76 100644 --- a/pygrt/C_extension/src/dynamic/grt_kernel.c +++ b/pygrt/C_extension/src/dynamic/grt_kernel.c @@ -26,7 +26,7 @@ typedef struct { bool active; char *s_modelpath; ///< 模型路径 const char *s_modelname; ///< 模型名称 - GRT_MODEL1D *mod1d; ///< 模型结构体指针 + MODEL1D *mod1d; ///< 模型结构体指针 } M; /** 震源和接收器深度 */ struct { @@ -390,7 +390,7 @@ int kernel_main(int argc, char **argv) if((Ctrl->M.mod1d = grt_read_mod1d_from_file(Ctrl->M.s_modelpath, Ctrl->D.depsrc, Ctrl->D.deprcv, true)) == NULL){ exit(EXIT_FAILURE); } - GRT_MODEL1D *mod1d = Ctrl->M.mod1d; + MODEL1D *mod1d = Ctrl->M.mod1d; // 边界条件 grt_set_mod1d_boundary(mod1d, Ctrl->B.topbound, Ctrl->B.botbound); @@ -430,7 +430,7 @@ int kernel_main(int argc, char **argv) real_t w = PI2*freq; cplx_t omega = w - I*Ctrl->F.wI; - GRT_MODEL1D *local_mod1d = NULL; + MODEL1D *local_mod1d = NULL; #ifdef _OPENMP // 定义局部模型对象 local_mod1d = grt_copy_mod1d(mod1d); diff --git a/pygrt/C_extension/src/dynamic/layer.c b/pygrt/C_extension/src/dynamic/layer.c index bc45b641..5bc90211 100644 --- a/pygrt/C_extension/src/dynamic/layer.c +++ b/pygrt/C_extension/src/dynamic/layer.c @@ -106,7 +106,7 @@ static int __rigidBound_R_SH(bool isLiquid, cplx_t *ML) return GRT_INVERSE_SUCCESS; } -void grt_topbound_RU_PSV(GRT_MODEL1D *mod1d) +void grt_topbound_RU_PSV(MODEL1D *mod1d) { cplx_t cbcb = mod1d->cbcb[0]; cplx_t xa = mod1d->xa[0]; @@ -128,7 +128,7 @@ void grt_topbound_RU_PSV(GRT_MODEL1D *mod1d) // RU 不需要时延 } -void grt_topbound_RU_SH(GRT_MODEL1D *mod1d) +void grt_topbound_RU_SH(MODEL1D *mod1d) { bool isLiquid = mod1d->isLiquid[0]; if(mod1d->topbound == GRT_BOUND_FREE){ @@ -147,7 +147,7 @@ void grt_topbound_RU_SH(GRT_MODEL1D *mod1d) // RU 不需要时延 } -void grt_botbound_RD_PSV(GRT_MODEL1D *mod1d) +void grt_botbound_RD_PSV(MODEL1D *mod1d) { size_t nlay = mod1d->n; cplx_t cbcb = mod1d->cbcb[nlay-2]; @@ -182,7 +182,7 @@ void grt_botbound_RD_PSV(GRT_MODEL1D *mod1d) mod1d->M_bot.RD[1][0] *= exab; mod1d->M_bot.RD[1][1] *= ex2b; } -void grt_botbound_RD_SH(GRT_MODEL1D *mod1d) +void grt_botbound_RD_SH(MODEL1D *mod1d) { size_t nlay = mod1d->n; cplx_t xb = mod1d->xb[nlay-2]; @@ -212,7 +212,7 @@ void grt_botbound_RD_SH(GRT_MODEL1D *mod1d) } -void grt_wave2qwv_REV_PSV(GRT_MODEL1D *mod1d) +void grt_wave2qwv_REV_PSV(MODEL1D *mod1d) { size_t ircv = mod1d->ircv; cplx_t xa = mod1d->xa[ircv]; @@ -247,7 +247,7 @@ void grt_wave2qwv_REV_PSV(GRT_MODEL1D *mod1d) } -void grt_wave2qwv_REV_SH(GRT_MODEL1D *mod1d) +void grt_wave2qwv_REV_SH(MODEL1D *mod1d) { size_t ircv = mod1d->ircv; bool isLiquid = mod1d->isLiquid[ircv]; @@ -268,7 +268,7 @@ void grt_wave2qwv_REV_SH(GRT_MODEL1D *mod1d) } -void grt_wave2qwv_z_REV_PSV(GRT_MODEL1D *mod1d) +void grt_wave2qwv_z_REV_PSV(MODEL1D *mod1d) { size_t ircv = mod1d->ircv; cplx_t xa = mod1d->xa[ircv]; @@ -302,7 +302,7 @@ void grt_wave2qwv_z_REV_PSV(GRT_MODEL1D *mod1d) } -void grt_wave2qwv_z_REV_SH(GRT_MODEL1D *mod1d) +void grt_wave2qwv_z_REV_SH(MODEL1D *mod1d) { size_t ircv = mod1d->ircv; cplx_t xb = mod1d->xb[ircv]; @@ -327,7 +327,7 @@ void grt_wave2qwv_z_REV_SH(GRT_MODEL1D *mod1d) } -void grt_RT_matrix_ll_PSV(const GRT_MODEL1D *mod1d, size_t iy, RT_MATRIX *M) +void grt_RT_matrix_ll_PSV(const MODEL1D *mod1d, size_t iy, RT_MATRIX *M) { MODEL_2LAYS_ATTRIB(cplx_t, xa); MODEL_2LAYS_ATTRIB(real_t, Rho); @@ -362,7 +362,7 @@ void grt_RT_matrix_ll_SH(RT_MATRIX *M) -void grt_RT_matrix_ls_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_RT_matrix_ls_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) { MODEL_2LAYS_ATTRIB(cplx_t, xa); MODEL_2LAYS_ATTRIB(cplx_t, xb); @@ -434,7 +434,7 @@ void grt_RT_matrix_ls_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX * } -void grt_RT_matrix_ls_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_RT_matrix_ls_SH(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) { MODEL_2LAYS_ATTRIB(bool, isLiquid); @@ -465,7 +465,7 @@ void grt_RT_matrix_ls_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M -void grt_RT_matrix_ss_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_RT_matrix_ss_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) { MODEL_2LAYS_ATTRIB(cplx_t, xa); MODEL_2LAYS_ATTRIB(cplx_t, xb); @@ -551,7 +551,7 @@ void grt_RT_matrix_ss_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX * } -void grt_RT_matrix_ss_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_RT_matrix_ss_SH(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) { MODEL_2LAYS_ATTRIB(cplx_t, xb); MODEL_2LAYS_ATTRIB(cplx_t, mu); @@ -569,7 +569,7 @@ void grt_RT_matrix_ss_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M -void grt_RT_matrix_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_RT_matrix_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) { MODEL_2LAYS_ATTRIB(bool, isLiquid); @@ -586,7 +586,7 @@ void grt_RT_matrix_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) } -void grt_RT_matrix_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_RT_matrix_SH(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) { MODEL_2LAYS_ATTRIB(bool, isLiquid); @@ -603,7 +603,7 @@ void grt_RT_matrix_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) } -void grt_delay_RT_matrix(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_delay_RT_matrix(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) { real_t thk = mod1d->Thk[iy-1]; cplx_t xa1 = mod1d->xa[iy-1]; @@ -633,7 +633,7 @@ void grt_delay_RT_matrix(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M M->TUL *= exb; } -void grt_delay_RT_matrix_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_delay_RT_matrix_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) { real_t thk = mod1d->Thk[iy-1]; cplx_t xa1 = mod1d->xa[iy-1]; @@ -657,7 +657,7 @@ void grt_delay_RT_matrix_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRI M->TU[1][0] *= exb; M->TU[1][1] *= exb; } -void grt_delay_RT_matrix_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_delay_RT_matrix_SH(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) { real_t thk = mod1d->Thk[iy-1]; cplx_t xb1 = mod1d->xb[iy-1]; @@ -674,7 +674,7 @@ void grt_delay_RT_matrix_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX -void grt_delay_GRT_matrix(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_delay_GRT_matrix(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) { real_t thk = mod1d->Thk[iy-1]; cplx_t xa1 = mod1d->xa[iy-1]; @@ -703,7 +703,7 @@ void grt_delay_GRT_matrix(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX * } -void grt_get_layer_D(const GRT_MODEL1D *mod1d, const size_t iy, bool inverse, int liquid_invtype, cplx_t D[4][4]) +void grt_get_layer_D(const MODEL1D *mod1d, const size_t iy, bool inverse, int liquid_invtype, cplx_t D[4][4]) { // 第iy层物理量 real_t k = mod1d->k; @@ -774,7 +774,7 @@ void grt_get_layer_D(const GRT_MODEL1D *mod1d, const size_t iy, bool inverse, in } -void grt_get_layer_D11(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]) +void grt_get_layer_D11(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]) { real_t k = mod1d->k; cplx_t xa = mod1d->xa[iy]; @@ -789,7 +789,7 @@ void grt_get_layer_D11(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2][2] } -void grt_get_layer_D12(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]) +void grt_get_layer_D12(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]) { real_t k = mod1d->k; cplx_t xa = mod1d->xa[iy]; @@ -803,7 +803,7 @@ void grt_get_layer_D12(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2][2] } } -void grt_get_layer_D11_uiz(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]) +void grt_get_layer_D11_uiz(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]) { real_t k = mod1d->k; cplx_t xa = mod1d->xa[iy]; @@ -820,7 +820,7 @@ void grt_get_layer_D11_uiz(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2 } } -void grt_get_layer_D12_uiz(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]) +void grt_get_layer_D12_uiz(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]) { real_t k = mod1d->k; cplx_t xa = mod1d->xa[iy]; @@ -837,7 +837,7 @@ void grt_get_layer_D12_uiz(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2 } } -void grt_get_layer_D21(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]) +void grt_get_layer_D21(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]) { real_t k = mod1d->k; cplx_t xa = mod1d->xa[iy]; @@ -857,7 +857,7 @@ void grt_get_layer_D21(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2][2] } -void grt_get_layer_D22(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]) +void grt_get_layer_D22(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]) { real_t k = mod1d->k; cplx_t xa = mod1d->xa[iy]; @@ -876,7 +876,7 @@ void grt_get_layer_D22(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2][2] } } -void grt_get_layer_T(const GRT_MODEL1D *mod1d, const size_t iy, bool inverse, cplx_t T[2][2]) +void grt_get_layer_T(const MODEL1D *mod1d, const size_t iy, bool inverse, cplx_t T[2][2]) { // 液体层不应该使用该函数 if(mod1d->isLiquid[iy]){ diff --git a/pygrt/C_extension/src/dynamic/source.c b/pygrt/C_extension/src/dynamic/source.c index 33756640..a3026f52 100755 --- a/pygrt/C_extension/src/dynamic/source.c +++ b/pygrt/C_extension/src/dynamic/source.c @@ -63,7 +63,7 @@ inline GCC_ALWAYS_INLINE void _source_SH(const cplx_t xb, const cplx_t cbcb, con } -void grt_source_coef(const GRT_MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU) +void grt_source_coef(const MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU) { // 先全部赋0 memset(src_coefD, 0, sizeof(cplxChnlGrid)); @@ -74,7 +74,7 @@ void grt_source_coef(const GRT_MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlG } -void grt_source_coef_PSV(const GRT_MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU) +void grt_source_coef_PSV(const MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU) { size_t isrc = mod1d->isrc; cplx_t xa = mod1d->xa[isrc]; @@ -87,7 +87,7 @@ void grt_source_coef_PSV(const GRT_MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxC } -void grt_source_coef_SH(const GRT_MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU) +void grt_source_coef_SH(const MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU) { size_t isrc = mod1d->isrc; cplx_t xb = mod1d->xb[isrc]; diff --git a/pygrt/C_extension/src/integral/dwm.c b/pygrt/C_extension/src/integral/dwm.c index 69fb21d6..7d2511bb 100644 --- a/pygrt/C_extension/src/integral/dwm.c +++ b/pygrt/C_extension/src/integral/dwm.c @@ -25,7 +25,7 @@ real_t grt_discrete_integ( - GRT_MODEL1D *mod1d, real_t dk, real_t kmax, real_t keps, + MODEL1D *mod1d, real_t dk, real_t kmax, real_t keps, size_t nr, real_t *rs, K_INTEG *K, FILE *fstats, GRT_KernelFunc kerfunc) { if(kmax == 0.0) return 0.0; diff --git a/pygrt/C_extension/src/integral/fim.c b/pygrt/C_extension/src/integral/fim.c index 7f5a5d46..8b5039a6 100755 --- a/pygrt/C_extension/src/integral/fim.c +++ b/pygrt/C_extension/src/integral/fim.c @@ -24,7 +24,7 @@ real_t grt_linear_filon_integ( - GRT_MODEL1D *mod1d, real_t k0, real_t dk0, real_t dk, real_t kmax, real_t keps, + MODEL1D *mod1d, real_t k0, real_t dk0, real_t dk, real_t kmax, real_t keps, size_t nr, real_t *rs, K_INTEG *K, FILE *fstats, GRT_KernelFunc kerfunc) { if(k0 + dk0 >= kmax) return k0; diff --git a/pygrt/C_extension/src/integral/integ_method.c b/pygrt/C_extension/src/integral/integ_method.c index 1f331755..38b5f7ca 100644 --- a/pygrt/C_extension/src/integral/integ_method.c +++ b/pygrt/C_extension/src/integral/integ_method.c @@ -78,7 +78,7 @@ void grt_KMET_destroy_fstats(const size_t nr, K_INTEG_METHOD *Kmet) K_INTEG * grt_wavenumber_integral( - GRT_MODEL1D *mod1d, size_t nr, real_t *rs, K_INTEG_METHOD *Kmet, bool calc_upar, GRT_KernelFunc kerfunc) + MODEL1D *mod1d, size_t nr, real_t *rs, K_INTEG_METHOD *Kmet, bool calc_upar, GRT_KernelFunc kerfunc) { real_t k = 0.0; diff --git a/pygrt/C_extension/src/integral/ptam.c b/pygrt/C_extension/src/integral/ptam.c index d5e0e5cf..985a8dc8 100755 --- a/pygrt/C_extension/src/integral/ptam.c +++ b/pygrt/C_extension/src/integral/ptam.c @@ -212,7 +212,7 @@ static void _cplx_shrink(size_t n1, size_t ir, int im, int v, cplxIntegGrid (*F void grt_PTA_method( - GRT_MODEL1D *mod1d, real_t k0, real_t predk, + MODEL1D *mod1d, real_t k0, real_t predk, size_t nr, real_t *rs, K_INTEG *K, FILE *ptam_fstatsnr[nr][2], GRT_KernelFunc kerfunc) { // 需要兼容对正常收敛而不具有规律波峰波谷的序列 diff --git a/pygrt/C_extension/src/integral/safim.c b/pygrt/C_extension/src/integral/safim.c index c0f0caa0..27c02e93 100644 --- a/pygrt/C_extension/src/integral/safim.c +++ b/pygrt/C_extension/src/integral/safim.c @@ -276,7 +276,7 @@ static void interv_integ(const KInterval *ptKitv, size_t nr, real_t *rs, K_INTEG real_t grt_sa_filon_integ( - GRT_MODEL1D *mod1d, real_t k0, real_t dk0, real_t tol, real_t kmax, real_t kref, + MODEL1D *mod1d, real_t k0, real_t dk0, real_t tol, real_t kmax, real_t kref, size_t nr, real_t *rs, K_INTEG *K, FILE *fstats, GRT_KernelFunc kerfunc) { real_t kmin = k0 + dk0; diff --git a/pygrt/C_extension/src/static/grt_static_greenfn.c b/pygrt/C_extension/src/static/grt_static_greenfn.c index e4fd564a..f10bd142 100644 --- a/pygrt/C_extension/src/static/grt_static_greenfn.c +++ b/pygrt/C_extension/src/static/grt_static_greenfn.c @@ -31,7 +31,7 @@ typedef struct { bool active; char *s_modelpath; ///< 模型路径 const char *s_modelname; ///< 模型名称 - GRT_MODEL1D *mod1d; ///< 模型结构体指针 + MODEL1D *mod1d; ///< 模型结构体指针 } M; /** 震源和接收器深度 */ struct { @@ -581,7 +581,7 @@ int static_greenfn_main(int argc, char **argv){ if((Ctrl->M.mod1d = grt_read_mod1d_from_file(Ctrl->M.s_modelpath, Ctrl->D.depsrc, Ctrl->D.deprcv, false)) == NULL){ exit(EXIT_FAILURE); } - GRT_MODEL1D *mod1d = Ctrl->M.mod1d; + MODEL1D *mod1d = Ctrl->M.mod1d; // 边界条件 grt_set_mod1d_boundary(mod1d, Ctrl->B.topbound, Ctrl->B.botbound); diff --git a/pygrt/C_extension/src/static/static_grn.c b/pygrt/C_extension/src/static/static_grn.c index 6198b4e5..03ac00b2 100644 --- a/pygrt/C_extension/src/static/static_grn.c +++ b/pygrt/C_extension/src/static/static_grn.c @@ -59,7 +59,7 @@ static void recordin_GRN( void grt_integ_static_grn( - GRT_MODEL1D *mod1d, size_t nr, real_t *rs, K_INTEG_METHOD *Kmet, + MODEL1D *mod1d, size_t nr, real_t *rs, K_INTEG_METHOD *Kmet, bool calc_upar, realChnlGrid grn[nr], realChnlGrid grn_uiz[nr], diff --git a/pygrt/C_extension/src/static/static_layer.c b/pygrt/C_extension/src/static/static_layer.c index 83924998..17d5907c 100644 --- a/pygrt/C_extension/src/static/static_layer.c +++ b/pygrt/C_extension/src/static/static_layer.c @@ -67,7 +67,7 @@ static int __rigidBound_R_SH(cplx_t *ML) return GRT_INVERSE_SUCCESS; } -void grt_static_topbound_RU_PSV(GRT_MODEL1D *mod1d) +void grt_static_topbound_RU_PSV(MODEL1D *mod1d) { cplx_t delta = mod1d->delta[0]; real_t k = mod1d->k; @@ -87,7 +87,7 @@ void grt_static_topbound_RU_PSV(GRT_MODEL1D *mod1d) // RU 不需要时延 } -void grt_static_topbound_RU_SH(GRT_MODEL1D *mod1d) +void grt_static_topbound_RU_SH(MODEL1D *mod1d) { if(mod1d->topbound == GRT_BOUND_FREE){ mod1d->M_top.stats = __freeBound_R_SH(&mod1d->M_top.RUL); @@ -105,7 +105,7 @@ void grt_static_topbound_RU_SH(GRT_MODEL1D *mod1d) // RU 不需要时延 } -void grt_static_botbound_RD_PSV(GRT_MODEL1D *mod1d) +void grt_static_botbound_RD_PSV(MODEL1D *mod1d) { size_t nlay = mod1d->n; cplx_t delta = mod1d->delta[nlay-2]; @@ -133,7 +133,7 @@ void grt_static_botbound_RD_PSV(GRT_MODEL1D *mod1d) mod1d->M_bot.RD[1][0] *= ex2; mod1d->M_bot.RD[1][1] *= ex2; } -void grt_static_botbound_RD_SH(GRT_MODEL1D *mod1d) +void grt_static_botbound_RD_SH(MODEL1D *mod1d) { size_t nlay = mod1d->n; real_t thk = mod1d->Thk[nlay-2]; @@ -159,7 +159,7 @@ void grt_static_botbound_RD_SH(GRT_MODEL1D *mod1d) mod1d->M_bot.RDL *= ex2; } -void grt_static_wave2qwv_REV_PSV(GRT_MODEL1D *mod1d) +void grt_static_wave2qwv_REV_PSV(MODEL1D *mod1d) { cplx_t D11[2][2] = {{1.0, -1.0}, {1.0, 1.0}}; cplx_t D12[2][2] = {{1.0, -1.0}, {-1.0, -1.0}}; @@ -174,7 +174,7 @@ void grt_static_wave2qwv_REV_PSV(GRT_MODEL1D *mod1d) } } -void grt_static_wave2qwv_REV_SH(GRT_MODEL1D *mod1d) +void grt_static_wave2qwv_REV_SH(MODEL1D *mod1d) { if(mod1d->ircvup){// 震源更深 mod1d->R_EVL = 1.0 + mod1d->M_FA.RUL; @@ -183,7 +183,7 @@ void grt_static_wave2qwv_REV_SH(GRT_MODEL1D *mod1d) } } -void grt_static_wave2qwv_z_REV_PSV(GRT_MODEL1D *mod1d) +void grt_static_wave2qwv_z_REV_PSV(MODEL1D *mod1d) { real_t k = mod1d->k; size_t ircv = mod1d->ircv; @@ -202,7 +202,7 @@ void grt_static_wave2qwv_z_REV_PSV(GRT_MODEL1D *mod1d) } } -void grt_static_wave2qwv_z_REV_SH(GRT_MODEL1D *mod1d) +void grt_static_wave2qwv_z_REV_SH(MODEL1D *mod1d) { real_t k = mod1d->k; // 新推导公式 @@ -214,7 +214,7 @@ void grt_static_wave2qwv_z_REV_SH(GRT_MODEL1D *mod1d) } -void grt_static_RT_matrix_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_static_RT_matrix_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) { MODEL_2LAYS_ATTRIB(cplx_t, mu); MODEL_2LAYS_ATTRIB(cplx_t, delta); @@ -256,7 +256,7 @@ void grt_static_RT_matrix_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATR } -void grt_static_RT_matrix_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_static_RT_matrix_SH(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) { MODEL_2LAYS_ATTRIB(cplx_t, mu); @@ -274,7 +274,7 @@ void grt_static_RT_matrix_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRI } -void grt_static_delay_RT_matrix(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_static_delay_RT_matrix(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) { real_t thk = mod1d->Thk[iy-1]; real_t k = mod1d->k; diff --git a/pygrt/C_extension/src/static/static_source.c b/pygrt/C_extension/src/static/static_source.c index 76eae716..3fb0044e 100644 --- a/pygrt/C_extension/src/static/static_source.c +++ b/pygrt/C_extension/src/static/static_source.c @@ -61,7 +61,7 @@ inline GCC_ALWAYS_INLINE void _source_SH(const real_t k, cplxChnlGrid coefD, cpl } -void grt_static_source_coef(const GRT_MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU) +void grt_static_source_coef(const MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU) { // 先全部赋0 memset(src_coefD, 0, sizeof(cplxChnlGrid)); @@ -72,7 +72,7 @@ void grt_static_source_coef(const GRT_MODEL1D *mod1d, cplxChnlGrid src_coefD, cp } -void grt_static_source_coef_PSV(const GRT_MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU) +void grt_static_source_coef_PSV(const MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU) { size_t isrc = mod1d->isrc; cplx_t delta = mod1d->delta[isrc]; @@ -82,7 +82,7 @@ void grt_static_source_coef_PSV(const GRT_MODEL1D *mod1d, cplxChnlGrid src_coefD } -void grt_static_source_coef_SH(const GRT_MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU) +void grt_static_source_coef_SH(const MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU) { real_t k = mod1d->k; diff --git a/pygrt/C_extension/src/tools/grt_travt.c b/pygrt/C_extension/src/tools/grt_travt.c index ca3e0d4d..d8ca03cb 100644 --- a/pygrt/C_extension/src/tools/grt_travt.c +++ b/pygrt/C_extension/src/tools/grt_travt.c @@ -20,7 +20,7 @@ typedef struct { struct { bool active; char *s_modelpath; - GRT_MODEL1D *mod1d; ///< 模型结构体指针 + MODEL1D *mod1d; ///< 模型结构体指针 } M; /** 震源和接收器深度 */ struct { @@ -223,7 +223,7 @@ int travt_main(int argc, char **argv){ if((Ctrl->M.mod1d = grt_read_mod1d_from_file(Ctrl->M.s_modelpath, Ctrl->D.depsrc, Ctrl->D.deprcv, true)) == NULL){ exit(EXIT_FAILURE); } - GRT_MODEL1D *mod1d = Ctrl->M.mod1d; + MODEL1D *mod1d = Ctrl->M.mod1d; printf("------------------------------------------------\n"); printf(" Distance(km) Tp(secs) Ts(secs) \n"); From a2d319716ecab7b1a608bbad2d7649ee3ca57548 Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Tue, 31 Mar 2026 22:10:04 +0800 Subject: [PATCH 2/3] remove prtdbg.h --- pygrt/C_extension/include/grt/common/prtdbg.h | 13 ------------- 1 file changed, 13 deletions(-) delete mode 100755 pygrt/C_extension/include/grt/common/prtdbg.h diff --git a/pygrt/C_extension/include/grt/common/prtdbg.h b/pygrt/C_extension/include/grt/common/prtdbg.h deleted file mode 100755 index a8adb9d5..00000000 --- a/pygrt/C_extension/include/grt/common/prtdbg.h +++ /dev/null @@ -1,13 +0,0 @@ -/** - * @file prtdbg.h - * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn) - * @date 2024-07-24 - * - * 仅用于代码初期debug使用 - */ - -#pragma once - -#define Print_GRTCOEF 0 -#define Print_GRTCOEF_2 0 -#define Print_SOURCECOEF 0 From 0919fbf97d1663ad11fa79af5065c3f83c197964 Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Wed, 1 Apr 2026 14:44:27 +0800 Subject: [PATCH 3/3] REFAC: split the model-related struct into two structs: MODEL1D and MODEL1D_STATE --- pygrt/C_extension/include/grt/common/model.h | 139 ++++--- pygrt/C_extension/include/grt/dynamic/layer.h | 83 ++-- .../C_extension/include/grt/dynamic/source.h | 8 +- pygrt/C_extension/include/grt/integral/dwm.h | 4 +- pygrt/C_extension/include/grt/integral/fim.h | 4 +- .../include/grt/integral/integ_method.h | 4 +- .../C_extension/include/grt/integral/kernel.h | 26 +- pygrt/C_extension/include/grt/integral/ptam.h | 4 +- .../C_extension/include/grt/integral/safim.h | 4 +- .../include/grt/static/static_layer.h | 40 +- .../include/grt/static/static_source.h | 8 +- pygrt/C_extension/src/common/model.c | 386 ++++++++++-------- pygrt/C_extension/src/dynamic/grn.c | 20 +- pygrt/C_extension/src/dynamic/grt_kernel.c | 18 +- pygrt/C_extension/src/dynamic/layer.c | 329 +++++++-------- pygrt/C_extension/src/dynamic/source.c | 30 +- pygrt/C_extension/src/integral/dwm.c | 6 +- pygrt/C_extension/src/integral/fim.c | 10 +- pygrt/C_extension/src/integral/integ_method.c | 22 +- .../src/integral/kernel_template.c_ | 101 +++-- pygrt/C_extension/src/integral/ptam.c | 6 +- pygrt/C_extension/src/integral/safim.c | 14 +- pygrt/C_extension/src/static/static_grn.c | 8 +- pygrt/C_extension/src/static/static_layer.c | 127 +++--- pygrt/C_extension/src/static/static_source.c | 18 +- pygrt/c_interfaces.py | 10 +- pygrt/c_structures.py | 56 +-- pygrt/pymod.py | 2 +- 28 files changed, 748 insertions(+), 739 deletions(-) diff --git a/pygrt/C_extension/include/grt/common/model.h b/pygrt/C_extension/include/grt/common/model.h index fd801d52..59e374b2 100755 --- a/pygrt/C_extension/include/grt/common/model.h +++ b/pygrt/C_extension/include/grt/common/model.h @@ -24,11 +24,7 @@ typedef struct { bool ircvup; ///< 接收点位于浅层, ircv < isrc bool io_depth; ///< 读取的模型首列为每层顶界面深度 bool srcrcv_isInserted; ///< 震源和台站是否已经以虚拟层的形式插入 - cplx_t omgref; ///< 参考圆频率, 用于计算衰减 - cplx_t omega; ///< 圆频率 - real_t k; ///< 波数 - cplx_t c_phase; ///< 当前相速度 real_t *Thk; ///< Thk[n], 最后一层厚度不使用(当作正无穷), km real_t *Dep; ///< Dep[n], 每一层顶界面深度，第一层必须为 0.0 @@ -41,6 +37,18 @@ typedef struct { real_t *Qbinv; ///< 1/Q_s bool *isLiquid; ///< 每层是否为液体层 + GRT_BOUND_TYPE topbound; ///< 顶界面的边界条件 + GRT_BOUND_TYPE botbound; ///< 底界面的边界条件 +} MODEL1D; + + +typedef struct { + MODEL1D *mod1d; + + cplx_t omega; ///< 圆频率 + real_t k; ///< 波数 + cplx_t c_phase; ///< 当前相速度 + cplx_t *mu; ///< mu[n] \f$ V_b^2 * \rho \f$ cplx_t *lambda; ///< lambda[n] \f$ V_a^2 * \rho - 2*\mu \f$ cplx_t *delta; ///< delta[n] \f$ (\lambda+\mu)/(\lambda+3*\mu) \f$ @@ -60,13 +68,7 @@ typedef struct { RT_MATRIX M_RS; RT_MATRIX M_FA; RT_MATRIX M_FB; - - - /* 顶界面的反射系数矩阵，仅 RU, RUL 有用 */ - GRT_BOUND_TYPE topbound; RT_MATRIX M_top; - /* 底界面的反射系数矩阵，仅 RD, RDL 有用 */ - GRT_BOUND_TYPE botbound; RT_MATRIX M_bot; /* 接收点处的接收矩阵 (转为位移u和位移导数uiz的(B_m, C_m, P_m)系分量) */ @@ -75,29 +77,12 @@ typedef struct { cplx_t uiz_R_EV[2][2]; cplx_t uiz_R_EVL; -} MODEL1D; - - -/** - * 打印 MODEL1D 模型参数信息，主要用于调试程序 - * - * @param[in] mod1d `MODEL1D` 结构体指针 - * - */ -void grt_print_mod1d(const MODEL1D *mod1d); - -/** - * 释放 `MODEL1D` 结构体指针 - * - * @param[out] mod1d `MODEL1D` 结构体指针 - */ -void grt_free_mod1d(MODEL1D *mod1d); +} MODEL1D_STATE; /** - * 初始化 MODEL1D 模型内存空间 + * 初始化 `MODEL1D` 内存空间 * * @param[in] n 模型层数 - * * @return `MODEL1D` 结构体指针 * */ @@ -113,29 +98,21 @@ MODEL1D * grt_init_mod1d(size_t n); MODEL1D * grt_copy_mod1d(const MODEL1D *mod1d1); /** - * 根据不同的 omega， 计算衰减系数，更新弹性模量 + * 扩容 `MODEL1D` 结构体 * * @param[in,out] mod1d `MODEL1D` 结构体指针 - * @param[in] omega 复数频率 - */ -void grt_attenuate_mod1d(MODEL1D *mod1d, cplx_t omega); - -/** - * 根据记录好的圆频率和波数，计算相速度和每层的 xa, xb, caca, cbcb - * - * @param[in,out] mod1d 模型结构体指针 - * @param[in] k 波数 + * @param[in] n 新层数 */ -void grt_mod1d_xa_xb(MODEL1D *mod1d, const real_t k); +void grt_realloc_mod1d(MODEL1D *mod1d, size_t n); /** - * 扩容 `MODEL1D` 结构体 + * 释放 `MODEL1D` 结构体指针 * - * @param[in,out] mod1d `MODEL1D` 结构体指针 - * @param[in] n 新层数 + * @param[out] mod1d `MODEL1D` 结构体指针 */ -void grt_realloc_mod1d(MODEL1D *mod1d, size_t n); +void grt_free_mod1d(MODEL1D *mod1d); + /** * 从文件中读取模型文件 @@ -144,12 +121,12 @@ void grt_realloc_mod1d(MODEL1D *mod1d, size_t n); * @param[in] depsrc 震源深度 * @param[in] deprcv 接收深度 * @param[in] allowLiquid 是否允许液体层 - * * @return `MODEL1D` 结构体指针 * */ MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, real_t deprcv, bool allowLiquid); + /** * 设置模型的边界条件，并对底界面做检查 * @@ -159,14 +136,68 @@ MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, real_t */ void grt_set_mod1d_boundary(MODEL1D *mod1d, GRT_BOUND_TYPE topbound, GRT_BOUND_TYPE botbound); + + +// ===================================================================================================================== +// ===================================================================================================================== + + /** - * 从模型文件中判断各个量的大致精度（字符串长度），以确定浮点数输出位数 + * 初始化 `MODEL1D_STATE` 内存空间 * - * @param[in] modelpath 模型文件路径 - * @param[out] diglen 每一列的最大字符串长度 + * @param[in] mod1d `MODEL1D` 结构体指针 + * @return `MODEL1D_STATE` 结构体指针 + */ +MODEL1D_STATE * grt_init_mod1d_state(MODEL1D *mod1d); + + +/** + * 复制 `MODEL1D_STATE` 结构体，但其中 mod1d 指针不变 + * + * @param[in] mstat1 `MODEL1D_STATE` 结构体指针 + * @return `MODEL1D_STATE` 结构体指针 + */ +MODEL1D_STATE * grt_copy_mod1d_state(const MODEL1D_STATE *mstat1); + + +/** + * 根据给定频率，设置衰减后的弹性参数；若omega实部小于0则视为弹性介质 + * + * @param[in,out] mstat `MODEL1D_STATE` 结构体指针 + * @param[in] omega 圆频率 + * + */ +void grt_update_mod1d_state_omega(MODEL1D_STATE *mstat, const cplx_t omega); + +/** + * 根据记录好的圆频率，给定波数，计算相速度和每层的 xa, xb, caca, cbcb + * + * @param[in,out] mstat `MODEL1D_STATE` 结构体指针 + * @param[in] k 波数 + */ +void grt_update_mod1d_state_k(MODEL1D_STATE *mstat, const real_t k); + + +/** + * 释放 `MODEL1D_STATE` 结构体指针 + * + * @param[out] mstat `MODEL1D_STATE` 结构体指针 + */ +void grt_free_mod1d_state(MODEL1D_STATE *mstat); + + +// ===================================================================================================================== +// ===================================================================================================================== + + +/** + * 打印 MODEL1D 模型参数信息，主要用于调试程序 + * + * @param[in] mod1d `MODEL1D` 结构体指针 * */ -void grt_get_model_diglen_from_file(const char *modelpath, size_t diglen[6]); +void grt_print_mod1d(const MODEL1D *mod); + /** * 浮点数比较，检查模型中是否存在该速度（不论Vp,Vs） @@ -174,7 +205,6 @@ void grt_get_model_diglen_from_file(const char *modelpath, size_t diglen[6]); * @param[in] mod1d 模型 * @param[in] vel 输入速度 * @param[in] tol 浮点数比较精度 - * * @return 是否存在 */ bool grt_check_vel_in_mod(const MODEL1D *mod1d, const real_t vel, const real_t tol); @@ -182,9 +212,10 @@ bool grt_check_vel_in_mod(const MODEL1D *mod1d, const real_t vel, const real_t t /** * 计算最大最小速度（非零值） * - * @param mod1d (in)`MODEL1D` 结构体指针 - * @param vmin (out)最小速度 - * @param vmax (out)最大速度 + * @param[in] mod1d `MODEL1D` 结构体指针 + * @param[out] vmin 最小速度 + * @param[out] vmax 最大速度 * */ -void grt_get_mod1d_vmin_vmax(const MODEL1D *mod1d, real_t *vmin, real_t *vmax); \ No newline at end of file +void grt_get_mod1d_vmin_vmax(const MODEL1D *mod1d, real_t *vmin, real_t *vmax); + diff --git a/pygrt/C_extension/include/grt/dynamic/layer.h b/pygrt/C_extension/include/grt/dynamic/layer.h index 7ab5b685..50db1555 100755 --- a/pygrt/C_extension/include/grt/dynamic/layer.h +++ b/pygrt/C_extension/include/grt/dynamic/layer.h @@ -17,62 +17,59 @@ #include "grt/common/const.h" #include "grt/common/RT_matrix.h" -/** - * 计算不同边界条件下顶界面的反射系数RU，其中自由表面的公式见(5.3.10-14) - * +/** 计算不同边界条件下顶界面的反射系数RU，其中自由表面的公式见(5.3.10-14) * - * @param[in,out] mod1d 模型结构体指针，结果保存在 M_top + * @param[in,out] mstat 模型结构体指针，结果保存在 M_bot * */ -void grt_topbound_RU_PSV(MODEL1D *mod1d); -void grt_topbound_RU_SH(MODEL1D *mod1d); +void grt_topbound_RU_PSV(MODEL1D_STATE *mstat); +void grt_topbound_RU_SH(MODEL1D_STATE *mstat); /** * 计算不同边界条件下底界面的反射系数RD，其中自由表面的公式见(5.3.10-14) * - * - * @param[in,out] mod1d 模型结构体指针，结果保存在 M_bot + * @param[in,out] mstat 模型结构体指针，结果保存在 M_bot * */ -void grt_botbound_RD_PSV(MODEL1D *mod1d); -void grt_botbound_RD_SH(MODEL1D *mod1d); +void grt_botbound_RD_PSV(MODEL1D_STATE *mstat); +void grt_botbound_RD_SH(MODEL1D_STATE *mstat); /** * 计算接收点位置的 P-SV 波接收矩阵，将波场转为位移，公式(5.2.19) + (5.7.7,25) * - * @param[in,out] mod1d 模型结构体指针，结果保存在 R_EV + * @param[in,out] mstat 模型结构体指针，结果保存在 R_EV * */ -void grt_wave2qwv_REV_PSV(MODEL1D *mod1d); +void grt_wave2qwv_REV_PSV(MODEL1D_STATE *mstat); /** * 计算接收点位置的 SH 波接收矩阵，将波场转为位移，公式(5.2.19) + (5.7.7,25) * - * @param[in,out] mod1d 模型结构体指针，结果保存在 R_EVL + * @param[in,out] mstat 模型结构体指针，结果保存在 R_EVL * */ -void grt_wave2qwv_REV_SH(MODEL1D *mod1d); +void grt_wave2qwv_REV_SH(MODEL1D_STATE *mstat); /** * 计算接收点位置的ui_z的 P-SV 波接收矩阵，即将波场转为ui_z。 * 公式本质是推导ui_z关于q_m, w_m, v_m的连接矩阵（就是应力推导过程的一部分） * - * @param[in,out] mod1d 模型结构体指针，结果保存在 uiz_R_EV + * @param[in,out] mstat 模型结构体指针，结果保存在 uiz_R_EV * */ -void grt_wave2qwv_z_REV_PSV(MODEL1D *mod1d); +void grt_wave2qwv_z_REV_PSV(MODEL1D_STATE *mstat); /** * 计算接收点位置的ui_z的 SH 波接收矩阵，即将波场转为ui_z。 * 公式本质是推导ui_z关于q_m, w_m, v_m的连接矩阵（就是应力推导过程的一部分） * - * @param[in,out] mod1d 模型结构体指针，结果保存在 uiz_R_EVL + * @param[in,out] mstat 模型结构体指针，结果保存在 uiz_R_EVL * */ -void grt_wave2qwv_z_REV_SH(MODEL1D *mod1d); +void grt_wave2qwv_z_REV_SH(MODEL1D_STATE *mstat); /** @@ -81,70 +78,70 @@ void grt_wave2qwv_z_REV_SH(MODEL1D *mod1d); * * @note 对公式(5.4.14)进行了重新整理。原公式各项之间的数量级差别过大，浮点数计算损失精度严重。 * - * @param[in] mod1d 模型结构体指针 + * @param[in] mstat 模型结构体指针 * @param[in] iy 层位索引 * @param[out] M R/T矩阵 * */ -void grt_RT_matrix_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); +void grt_RT_matrix_PSV(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M); /** * 计算界面的 SH 波反射透射系数 RDL/RUL/TDL/TUL, * 根据公式(5.4.31)计算系数 * - * @param[in] mod1d 模型结构体指针 + * @param[in] mstat 模型结构体指针 * @param[in] iy 层位索引 * @param[out] M R/T矩阵 * */ -void grt_RT_matrix_SH(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); +void grt_RT_matrix_SH(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M); /** 液-液 界面 */ -void grt_RT_matrix_ll_PSV(const MODEL1D *mod1d, size_t iy, RT_MATRIX *M); +void grt_RT_matrix_ll_PSV(const MODEL1D_STATE *mstat, size_t iy, RT_MATRIX *M); /** 液-液 界面 */ void grt_RT_matrix_ll_SH(RT_MATRIX *M); /** 液-固 界面 */ -void grt_RT_matrix_ls_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); +void grt_RT_matrix_ls_PSV(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M); /** 液-固 界面 */ -void grt_RT_matrix_ls_SH(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); +void grt_RT_matrix_ls_SH(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M); /** 固-固 界面 */ -void grt_RT_matrix_ss_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); +void grt_RT_matrix_ss_PSV(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M); /** 固-固 界面 */ -void grt_RT_matrix_ss_SH(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); +void grt_RT_matrix_ss_SH(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M); /** * 为 R/T 矩阵添加时间延迟因子 * - * @param[in] mod1d 模型结构体指针 + * @param[in] mstat 模型结构体指针 * @param[in] iy 层位索引 * @param[out] M R/T矩阵 * */ -void grt_delay_RT_matrix(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); -void grt_delay_RT_matrix_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); -void grt_delay_RT_matrix_SH(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); +void grt_delay_RT_matrix(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M); +void grt_delay_RT_matrix_PSV(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M); +void grt_delay_RT_matrix_SH(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M); /** * 为虚拟层的广义 R/T 矩阵添加时间延迟因子 * - * @param[in] mod1d 模型结构体指针 + * @param[in] mstat 模型结构体指针 * @param[in] iy 层位索引 * @param[out] M R/T矩阵 * */ -void grt_delay_GRT_matrix(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); +void grt_delay_GRT_matrix(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M); /** * 计算该层的连接 P-SV 应力位移矢量与垂直波函数的D矩阵(或其逆矩阵)， * 见公式(5.2.19-20) * - * @param[in] mod1d 模型结构体指针 + * @param[in] mstat 模型结构体指针 * @param[in] iy 层位索引 * @param[in] inverse 是否生成逆矩阵 * @param[in] liquid_invtype 对于液体层，矩阵会有很多零，至少第二列、第四列和第四行均为零； @@ -155,38 +152,38 @@ void grt_delay_GRT_matrix(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); * @param[out] D D矩阵(或其逆矩阵) * */ -void grt_get_layer_D(const MODEL1D *mod1d, const size_t iy, bool inverse, int liquid_invtype, cplx_t D[4][4]); +void grt_get_layer_D(const MODEL1D_STATE *mstat, const size_t iy, bool inverse, int liquid_invtype, cplx_t D[4][4]); /** 子矩阵 D11，函数参数见 get_layer_D 函数 */ -void grt_get_layer_D11(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]); +void grt_get_layer_D11(const MODEL1D_STATE *mstat, const size_t iy, cplx_t D[2][2]); /** 子矩阵 D12，函数参数见 get_layer_D 函数 */ -void grt_get_layer_D12(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]); +void grt_get_layer_D12(const MODEL1D_STATE *mstat, const size_t iy, cplx_t D[2][2]); /** 子矩阵 D21，函数参数见 get_layer_D 函数 */ -void grt_get_layer_D21(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]); +void grt_get_layer_D21(const MODEL1D_STATE *mstat, const size_t iy, cplx_t D[2][2]); /** 子矩阵 D22，函数参数见 get_layer_D 函数 */ -void grt_get_layer_D22(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]); +void grt_get_layer_D22(const MODEL1D_STATE *mstat, const size_t iy, cplx_t D[2][2]); /** 子矩阵 D11_uiz，后缀uiz表示连接位移对z的偏导和垂直波函数，函数参数见 get_layer_D 函数 */ -void grt_get_layer_D11_uiz(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]); +void grt_get_layer_D11_uiz(const MODEL1D_STATE *mstat, const size_t iy, cplx_t D[2][2]); /** 子矩阵 D12_uiz，函数参数见 get_layer_D 函数 */ -void grt_get_layer_D12_uiz(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]); +void grt_get_layer_D12_uiz(const MODEL1D_STATE *mstat, const size_t iy, cplx_t D[2][2]); /** * 计算该层的连接 SH 应力位移矢量与垂直波函数的 T 矩阵(或其逆矩阵)， * 见公式(5.2.21-22) * - * @param[in] mod1d 模型结构体指针 + * @param[in] mstat 模型结构体指针 * @param[in] iy 层位索引 * @param[in] inverse 是否生成逆矩阵 * @param[out] T T矩阵(或其逆矩阵) * */ -void grt_get_layer_T(const MODEL1D *mod1d, const size_t iy, bool inverse, cplx_t T[2][2]); +void grt_get_layer_T(const MODEL1D_STATE *mstat, const size_t iy, bool inverse, cplx_t T[2][2]); // /** diff --git a/pygrt/C_extension/include/grt/dynamic/source.h b/pygrt/C_extension/include/grt/dynamic/source.h index 545f042e..d5cb0c95 100755 --- a/pygrt/C_extension/include/grt/dynamic/source.h +++ b/pygrt/C_extension/include/grt/dynamic/source.h @@ -19,15 +19,15 @@ * 根据公式(4.6.6)，(4.6.15)，(4.6.21,26)，(4.8.34)计算不同震源不同阶数的震源系数， * 数组形状代表在[i][j]时表示i类震源的P(j=0),SV(j=1),SH(j=2)的震源系数(分别对应q,w,v). * - * @param[in] mod1d 模型结构体指针 + * @param[in] mstat 模型结构体指针 * @param[out] src_coefD 下行震源系数 * @param[out] src_coefU 上行震源系数 * */ -void grt_source_coef(const MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU); +void grt_source_coef(const MODEL1D_STATE *mstat, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU); /* P-SV 波的震源系数 */ -void grt_source_coef_PSV(const MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU); +void grt_source_coef_PSV(const MODEL1D_STATE *mstat, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU); /* SH 波的震源系数 */ -void grt_source_coef_SH(const MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU); +void grt_source_coef_SH(const MODEL1D_STATE *mstat, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU); diff --git a/pygrt/C_extension/include/grt/integral/dwm.h b/pygrt/C_extension/include/grt/integral/dwm.h index b1411c29..678a7d8c 100644 --- a/pygrt/C_extension/include/grt/integral/dwm.h +++ b/pygrt/C_extension/include/grt/integral/dwm.h @@ -24,7 +24,7 @@ * 传统的离散波数积分，结果以三维数组的形式返回，形状分别代表震中距、不同震源不同阶数 * 和4种积分类型(p=0,1,2,3) * - * @param[in,out] mod1d `MODEL1D` 结构体指针 + * @param[in,out] mstat `MODEL1D_STATE` 结构体指针 * @param[in] dk 波数积分间隔 * @param[in] kmax 波数积分的上限 * @param[in] keps 波数积分的收敛条件，要求在某震中距下所有格林函数都收敛，为负数代表不提前判断收敛，按照波数积分上限进行积分 @@ -39,5 +39,5 @@ * @return k 积分截至时的波数 */ real_t grt_discrete_integ( - MODEL1D *mod1d, real_t dk, real_t kmax, real_t keps, + MODEL1D_STATE *mstat, real_t dk, real_t kmax, real_t keps, size_t nr, real_t *rs, K_INTEG *K, FILE *fstats, GRT_KernelFunc kerfunc); diff --git a/pygrt/C_extension/include/grt/integral/fim.h b/pygrt/C_extension/include/grt/integral/fim.h index 0440564b..2899277e 100755 --- a/pygrt/C_extension/include/grt/integral/fim.h +++ b/pygrt/C_extension/include/grt/integral/fim.h @@ -28,7 +28,7 @@ * 其中\f$x=kr\f$. * * - * @param[in,out] mod1d `MODEL1D` 结构体指针 + * @param[in,out] mstat `MODEL1D_STATE` 结构体指针 * @param[in] k0 前一部分的波数积分结束点k值 * @param[in] dk0 前一部分的波数积分间隔 * @param[in] filondk filon积分间隔 @@ -45,7 +45,7 @@ * @return k 积分截至时的波数 */ real_t grt_linear_filon_integ( - MODEL1D *mod1d, real_t k0, real_t dk0, real_t dk, real_t kmax, real_t keps, + MODEL1D_STATE *mstat, real_t k0, real_t dk0, real_t dk, real_t kmax, real_t keps, size_t nr, real_t *rs, K_INTEG *K, FILE *fstats, GRT_KernelFunc kerfunc); diff --git a/pygrt/C_extension/include/grt/integral/integ_method.h b/pygrt/C_extension/include/grt/integral/integ_method.h index cc26ffd6..b9abbe70 100644 --- a/pygrt/C_extension/include/grt/integral/integ_method.h +++ b/pygrt/C_extension/include/grt/integral/integ_method.h @@ -74,7 +74,7 @@ void grt_KMET_destroy_fstats(const size_t nr, K_INTEG_METHOD *Kmet); /** * 发起波数积分的总函数 * - * @param[in,out] mod1d `MODEL1D` 结构体指针 + * @param[in,out] mstat `MODEL1D_STATE` 结构体指针 * @param[in] nr 震中距数量 * @param[in] rs 震中距数组 * @param[in,out] Kmet K_INTEG_METHOD 结构体指针 @@ -85,4 +85,4 @@ void grt_KMET_destroy_fstats(const size_t nr, K_INTEG_METHOD *Kmet); * */ K_INTEG * grt_wavenumber_integral( - MODEL1D *mod1d, size_t nr, real_t *rs, K_INTEG_METHOD *Kmet, bool calc_upar, GRT_KernelFunc kerfunc); + MODEL1D_STATE *mstat, size_t nr, real_t *rs, K_INTEG_METHOD *Kmet, bool calc_upar, GRT_KernelFunc kerfunc); diff --git a/pygrt/C_extension/include/grt/integral/kernel.h b/pygrt/C_extension/include/grt/integral/kernel.h index 7419adc2..dacde939 100644 --- a/pygrt/C_extension/include/grt/integral/kernel.h +++ b/pygrt/C_extension/include/grt/integral/kernel.h @@ -15,9 +15,7 @@ /** * 计算核函数的函数指针，动态与静态的接口一致 */ -typedef void (*GRT_KernelFunc) ( - MODEL1D *mod1d, const real_t k, cplxChnlGrid QWV, - bool calc_uiz, cplxChnlGrid QWV_uiz); +typedef void (*GRT_KernelFunc) (MODEL1D_STATE *mstat, const real_t k, cplxChnlGrid QWV, bool calc_uiz, cplxChnlGrid QWV_uiz); /** @@ -81,34 +79,26 @@ typedef void (*GRT_KernelFunc) ( * 即空间划分为FA,AB,BL, 计算这三个广义层的系数矩阵，再讨论震源层和接收层的深浅， * 计算相应的矩阵。 * - * @param[in,out] mod1d `MODEL1D` 结构体指针 + * @param[in,out] mstat `MODEL1D_STATE` 结构体指针 * @param[in] k 波数 * @param[out] QWV 不同震源，不同阶数的核函数 \f$ q_m, w_m, v_m \f$ * @param[in] calc_uiz 是否计算ui_z（位移u对坐标z的偏导） * @param[out] QWV_uiz 不同震源，不同阶数的核函数对z的偏导 \f$ \frac{\partial q_m}{\partial z}, \frac{\partial w_m}{\partial z}, \frac{\partial v_m}{\partial z} \f$ * */ -void grt_kernel( - MODEL1D *mod1d, const real_t k, cplxChnlGrid QWV, - bool calc_uiz, cplxChnlGrid QWV_uiz); +void grt_kernel(MODEL1D_STATE *mstat, const real_t k, cplxChnlGrid QWV, bool calc_uiz, cplxChnlGrid QWV_uiz); /** 构建广义反射透射系数矩阵。作为 kernel 函数中的第一部分 */ -void grt_GRT_matrix(MODEL1D *mod1d, const real_t k); +void grt_GRT_matrix(MODEL1D_STATE *mstat, const real_t k); /** 从广义 R/T 矩阵出发，计算每个震源对应的核函数 QWV。 作为 kernel 函数中的第二部分 */ -void grt_GRT_build_QWV( - MODEL1D *mod1d, cplxChnlGrid QWV, - bool calc_uiz, cplxChnlGrid QWV_uiz); +void grt_GRT_build_QWV(MODEL1D_STATE *mstat, cplxChnlGrid QWV, bool calc_uiz, cplxChnlGrid QWV_uiz); /** 静态解的核函数 */ -void grt_static_kernel( - MODEL1D *mod1d, const real_t k, cplxChnlGrid QWV, - bool calc_uiz, cplxChnlGrid QWV_uiz); +void grt_static_kernel(MODEL1D_STATE *mstat, const real_t k, cplxChnlGrid QWV, bool calc_uiz, cplxChnlGrid QWV_uiz); /** 静态广义反射透射系数矩阵 */ -void grt_static_GRT_matrix(MODEL1D *mod1d, const real_t k); +void grt_static_GRT_matrix(MODEL1D_STATE *mstat, const real_t k); /** 静态 QWV */ -void grt_static_GRT_build_QWV( - MODEL1D *mod1d, cplxChnlGrid QWV, - bool calc_uiz, cplxChnlGrid QWV_uiz); \ No newline at end of file +void grt_static_GRT_build_QWV(MODEL1D_STATE *mstat, cplxChnlGrid QWV, bool calc_uiz, cplxChnlGrid QWV_uiz); \ No newline at end of file diff --git a/pygrt/C_extension/include/grt/integral/ptam.h b/pygrt/C_extension/include/grt/integral/ptam.h index 4b081daf..a8bca4ee 100755 --- a/pygrt/C_extension/include/grt/integral/ptam.h +++ b/pygrt/C_extension/include/grt/integral/ptam.h @@ -26,7 +26,7 @@ /** * 峰谷平均法 Peak-Trough Averaging Method，最后收敛的积分结果以三维数组的形式返回， * - * @param[in,out] mod1d `MODEL1D` 结构体指针 + * @param[in,out] mstat `MODEL1D_STATE` 结构体指针 * @param[in] k0 先前的积分已经进行到了波数k0 * @param[in] predk 先前的积分使用的积分间隔dk，因为峰谷平均法使用的 * 积分间隔会和之前的不一致，这里传入该系数以做预先调整 @@ -41,6 +41,6 @@ * */ void grt_PTA_method( - MODEL1D *mod1d, real_t k0, real_t predk, + MODEL1D_STATE *mstat, real_t k0, real_t predk, size_t nr, real_t *rs, K_INTEG *K, FILE *ptam_fstatsnr[nr][2], GRT_KernelFunc kerfunc); diff --git a/pygrt/C_extension/include/grt/integral/safim.h b/pygrt/C_extension/include/grt/integral/safim.h index f71afafc..b650c4d4 100755 --- a/pygrt/C_extension/include/grt/integral/safim.h +++ b/pygrt/C_extension/include/grt/integral/safim.h @@ -30,7 +30,7 @@ * 其中\f$x=kr\f$. * * - * @param[in,out] mod1d `MODEL1D` 结构体指针 + * @param[in,out] mstat `MODEL1D_STATE` 结构体指针 * @param[in] k0 前一部分的波数积分结束点k值 * @param[in] dk0 前一部分的波数积分间隔 * @param[in] tol 自适应Filon积分的采样精度 @@ -47,7 +47,7 @@ * @return k 积分截至时的波数 */ real_t grt_sa_filon_integ( - MODEL1D *mod1d, real_t k0, real_t dk0, real_t tol, real_t kmax, real_t kref, + MODEL1D_STATE *mstat, real_t k0, real_t dk0, real_t tol, real_t kmax, real_t kref, size_t nr, real_t *rs, K_INTEG *K, FILE *fstats, GRT_KernelFunc kerfunc); diff --git a/pygrt/C_extension/include/grt/static/static_layer.h b/pygrt/C_extension/include/grt/static/static_layer.h index e1a18eea..9f889013 100644 --- a/pygrt/C_extension/include/grt/static/static_layer.h +++ b/pygrt/C_extension/include/grt/static/static_layer.h @@ -22,85 +22,85 @@ /** * 计算不同边界条件下顶界面的静态反射系数RU，其中自由表面的公式见(6.3.12) * - * @param[in,out] mod1d 模型结构体指针，结果保存在 M_top + * @param[in,out] mstat 模型结构体指针，结果保存在 M_top * */ -void grt_static_topbound_RU_PSV(MODEL1D *mod1d); -void grt_static_topbound_RU_SH(MODEL1D *mod1d); +void grt_static_topbound_RU_PSV(MODEL1D_STATE *mstat); +void grt_static_topbound_RU_SH(MODEL1D_STATE *mstat); /** * 计算不同边界条件下底界面的静态反射系数RD，其中自由表面的公式见(6.3.12) * - * @param[in,out] mod1d 模型结构体指针，结果保存在 M_bot + * @param[in,out] mstat 模型结构体指针，结果保存在 M_bot * */ -void grt_static_botbound_RD_PSV(MODEL1D *mod1d); -void grt_static_botbound_RD_SH(MODEL1D *mod1d); +void grt_static_botbound_RD_PSV(MODEL1D_STATE *mstat); +void grt_static_botbound_RD_SH(MODEL1D_STATE *mstat); /** * 计算接收点位置的 P-SV 静态接收矩阵，将波场转为位移，公式(6.3.35) * - * @param[in,out] mod1d 模型结构体指针，结果保存在 R_EV + * @param[in,out] mstat 模型结构体指针，结果保存在 R_EV * */ -void grt_static_wave2qwv_REV_PSV(MODEL1D *mod1d); +void grt_static_wave2qwv_REV_PSV(MODEL1D_STATE *mstat); /** * 计算接收点位置的 SH 静态接收矩阵，将波场转为位移，公式(6.3.37) * - * @param[in,out] mod1d 模型结构体指针，结果保存在 R_EVL + * @param[in,out] mstat 模型结构体指针，结果保存在 R_EVL * */ -void grt_static_wave2qwv_REV_SH(MODEL1D *mod1d); +void grt_static_wave2qwv_REV_SH(MODEL1D_STATE *mstat); /** * 计算接收点位置的ui_z的 P-SV 静态接收矩阵，即将波场转为ui_z。 * 公式本质是推导ui_z关于q_m, w_m, v_m的连接矩阵（就是应力推导过程的一部分） * - * @param[in,out] mod1d 模型结构体指针，结果保存在 uiz_R_EV + * @param[in,out] mstat 模型结构体指针，结果保存在 uiz_R_EV * */ -void grt_static_wave2qwv_z_REV_PSV(MODEL1D *mod1d); +void grt_static_wave2qwv_z_REV_PSV(MODEL1D_STATE *mstat); /** * 计算接收点位置的ui_z的 SH 静态接收矩阵，即将波场转为ui_z。 * 公式本质是推导ui_z关于q_m, w_m, v_m的连接矩阵（就是应力推导过程的一部分） * - * @param[in,out] mod1d 模型结构体指针，结果保存在 uiz_R_EVL + * @param[in,out] mstat 模型结构体指针，结果保存在 uiz_R_EVL * */ -void grt_static_wave2qwv_z_REV_SH(MODEL1D *mod1d); +void grt_static_wave2qwv_z_REV_SH(MODEL1D_STATE *mstat); /** * 计算界面的 P-SV 波静态反射透射系数RD/RU/TD/TU * 根据公式(6.3.18) * - * @param[in] mod1d 模型结构体指针 + * @param[in] mstat 模型结构体指针 * @param[in] iy 层位索引 * @param[out] M R/T矩阵 * */ -void grt_static_RT_matrix_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); +void grt_static_RT_matrix_PSV(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M); /** * 计算界面的 SH 波静态反射透射系数RDL/RUL/TDL/TUL * 根据公式(6.3.18) * - * @param[in] mod1d 模型结构体指针 + * @param[in] mstat 模型结构体指针 * @param[in] iy 层位索引 * @param[out] M R/T矩阵 * */ -void grt_static_RT_matrix_SH(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); +void grt_static_RT_matrix_SH(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M); /** * 为静态 R/T 矩阵添加时间延迟因子 * - * @param[in] mod1d 模型结构体指针 + * @param[in] mstat 模型结构体指针 * @param[in] iy 层位索引 * @param[out] M R/T矩阵 * */ -void grt_static_delay_RT_matrix(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M); \ No newline at end of file +void grt_static_delay_RT_matrix(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M); \ No newline at end of file diff --git a/pygrt/C_extension/include/grt/static/static_source.h b/pygrt/C_extension/include/grt/static/static_source.h index 821fc23f..78555e10 100644 --- a/pygrt/C_extension/include/grt/static/static_source.h +++ b/pygrt/C_extension/include/grt/static/static_source.h @@ -18,15 +18,15 @@ * * 数组形状代表在[i][j]时表示i类震源的P(j=0),SV(j=1),SH(j=2)的震源系数(分别对应q,w,v). * - * @param[in] mod1d 模型结构体指针 + * @param[in] mstat 模型结构体指针 * @param[out] src_coefD 下行震源系数 * @param[out] src_coefU 上行震源系数 * */ -void grt_static_source_coef(const MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU); +void grt_static_source_coef(const MODEL1D_STATE *mstat, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU); /* P-SV 波的静态震源系数 */ -void grt_static_source_coef_PSV(const MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU); +void grt_static_source_coef_PSV(const MODEL1D_STATE *mstat, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU); /* SH 波的静态震源系数 */ -void grt_static_source_coef_SH(const MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU); +void grt_static_source_coef_SH(const MODEL1D_STATE *mstat, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU); diff --git a/pygrt/C_extension/src/common/model.c b/pygrt/C_extension/src/common/model.c index c128cdf8..4087c799 100755 --- a/pygrt/C_extension/src/common/model.c +++ b/pygrt/C_extension/src/common/model.c @@ -13,7 +13,6 @@ #include #include "grt/common/model.h" -#include "grt/common/prtdbg.h" #include "grt/common/attenuation.h" #include "grt/common/util.h" @@ -21,7 +20,7 @@ // 定义宏，方便写代码 -#define GRT_FOR_EACH_MODEL_QUANTITY_ARRAY \ +#define __MODEL1D_FOR_EACH_ARRAY \ X(Thk, real_t)\ X(Dep, real_t)\ X(Va, real_t)\ @@ -31,7 +30,9 @@ X(Qb, real_t)\ X(Qainv, real_t)\ X(Qbinv, real_t)\ - X(isLiquid, bool)\ + X(isLiquid, bool) + +#define __MODEL1D_STATE_FOR_EACH_ARRAY \ X(mu, cplx_t)\ X(lambda, cplx_t)\ X(delta, cplx_t)\ @@ -42,97 +43,21 @@ X(caca, cplx_t)\ X(cbcb, cplx_t)\ - -void grt_print_mod1d(const MODEL1D *mod1d){ - // 模拟表格，打印速度 - // 每列字符宽度 - // [isrc/ircv] [h(km)] [Vp(km/s)] [Vs(km/s)] [Rho(g/cm^3)] [Qp] [Qs] - const int ncols = 7; - const int nlens[] = {13, 12, 13, 13, 16, 13, 13}; - int Nlen=0; - for(int ic=0; icn; ++i){ - if(i==mod1d->isrc){ - snprintf(indexstr, sizeof(indexstr), "%zu [src]", i+1); - } else if(i==mod1d->ircv){ - snprintf(indexstr, sizeof(indexstr), "%zu [rcv]", i+1); - } else { - snprintf(indexstr, sizeof(indexstr), "%zu ", i+1); - } - - printf("| %*s ", nlens[0]-3, indexstr); - - if(i < mod1d->n-1){ - printf("| %-*.2f ", nlens[1]-3, mod1d->Thk[i]); - } else { - printf("| %-*s ", nlens[1]-3, "Inf"); - } - - printf("| %-*.2f ", nlens[2]-3, mod1d->Va[i]); - printf("| %-*.2f ", nlens[3]-3, mod1d->Vb[i]); - printf("| %-*.2f ", nlens[4]-3, mod1d->Rho[i]); - printf("| %-*.2e ", nlens[5]-3, mod1d->Qa[i]); - printf("| %-*.2e ", nlens[6]-3, mod1d->Qb[i]); - printf("|\n"); - } - printf("%s\n", splitline); - printf("\n"); -} - -void grt_free_mod1d(MODEL1D *mod1d){ - #define X(P, T) GRT_SAFE_FREE_PTR(mod1d->P); - GRT_FOR_EACH_MODEL_QUANTITY_ARRAY - #undef X - - GRT_SAFE_FREE_PTR(mod1d); -} - - -MODEL1D * grt_init_mod1d(size_t n){ + +MODEL1D * grt_init_mod1d(size_t n) +{ MODEL1D *mod1d = (MODEL1D *)calloc(1, sizeof(MODEL1D)); mod1d->n = n; - mod1d->omgref = PI2*1.0; #define X(P, T) mod1d->P = (T*)calloc(n, sizeof(T)); - GRT_FOR_EACH_MODEL_QUANTITY_ARRAY + __MODEL1D_FOR_EACH_ARRAY #undef X return mod1d; } - -MODEL1D * grt_copy_mod1d(const MODEL1D *mod1d1){ +MODEL1D * grt_copy_mod1d(const MODEL1D *mod1d1) +{ MODEL1D *mod1d2 = (MODEL1D *)calloc(1, sizeof(MODEL1D)); // 先直接赋值，实现浅拷贝 @@ -144,91 +69,33 @@ MODEL1D * grt_copy_mod1d(const MODEL1D *mod1d1){ mod1d2->P = (T*)calloc(n, sizeof(T));\ memcpy(mod1d2->P, mod1d1->P, sizeof(T)*n);\ - GRT_FOR_EACH_MODEL_QUANTITY_ARRAY + __MODEL1D_FOR_EACH_ARRAY #undef X return mod1d2; } - -void grt_attenuate_mod1d(MODEL1D *mod1d, cplx_t omega){ - real_t Va0, Vb0; - cplx_t atna, atnb; - for(size_t i=0; in; ++i){ - Va0 = mod1d->Va[i]; - Vb0 = mod1d->Vb[i]; - - // 圆频率实部为负数表明不考虑模型的 Q 值属性 - // 在读入模型后需要需要运行一次本函数以填充弹性模量，见 grt_read_mod1d_from_file 函数 - atna = (creal(omega) >= 0.0 && mod1d->Qainv[i] > 0.0)? grt_attenuation_law(mod1d->Qainv[i], mod1d->omgref, omega) : 1.0; - atnb = (creal(omega) >= 0.0 && mod1d->Qbinv[i] > 0.0)? grt_attenuation_law(mod1d->Qbinv[i], mod1d->omgref, omega) : 1.0; - - mod1d->atna[i] = atna; - mod1d->atnb[i] = atnb; - - mod1d->mu[i] = (Vb0*atnb)*(Vb0*atnb)*(mod1d->Rho[i]); - mod1d->lambda[i] = (Va0*atnb)*(Va0*atnb)*(mod1d->Rho[i]) - 2*mod1d->mu[i]; - mod1d->delta[i] = (mod1d->lambda[i] + mod1d->mu[i]) / (mod1d->lambda[i] + 3.0*mod1d->mu[i]); - } - -#if Print_GRTCOEF == 1 - print_mod1d(mod1d); -#endif -} - - -void grt_mod1d_xa_xb(MODEL1D *mod1d, const real_t k) +void grt_realloc_mod1d(MODEL1D *mod1d, size_t n) { - mod1d->k = k; - // 不合理的频率值，只可能是在计算静态解，此时不需要xa, xb等物理量 - if(creal(mod1d->omega) < 0.0) return; - - mod1d->c_phase = mod1d->omega/k; - - size_t isrc = mod1d->isrc; - size_t ircv = mod1d->ircv; - - for(size_t i=0; in; ++i){ - if( mod1d->srcrcv_isInserted && (i == isrc || i == ircv) ){ - mod1d->xa[i] = mod1d->xa[i-1]; - mod1d->caca[i] = mod1d->caca[i-1]; - mod1d->xb[i] = mod1d->xb[i-1]; - mod1d->cbcb[i] = mod1d->cbcb[i-1]; - continue; - } - - real_t va, vb; - va = mod1d->Va[i]; - vb = mod1d->Vb[i]; - cplx_t atna, atnb; - atna = mod1d->atna[i]; - atnb = mod1d->atnb[i]; - - cplx_t caca, cbcb; - caca = mod1d->c_phase / (va*atna); - caca *= caca; - mod1d->caca[i] = caca; - mod1d->xa[i] = sqrt(1.0 - caca); - - cbcb = (mod1d->isLiquid[i])? 0.0 : mod1d->c_phase / (vb*atnb); // 考虑液体层 - cbcb *= cbcb; - mod1d->cbcb[i] = cbcb; - mod1d->xb[i] = sqrt(1.0 - cbcb); - } -} - - -void grt_realloc_mod1d(MODEL1D *mod1d, size_t n){ mod1d->n = n; #define X(P, T) mod1d->P = (T*)realloc(mod1d->P, n*sizeof(T)); - GRT_FOR_EACH_MODEL_QUANTITY_ARRAY + __MODEL1D_FOR_EACH_ARRAY #undef X } +void grt_free_mod1d(MODEL1D *mod1d) +{ + #define X(P, T) GRT_SAFE_FREE_PTR(mod1d->P); + __MODEL1D_FOR_EACH_ARRAY + #undef X + GRT_SAFE_FREE_PTR(mod1d); +} -MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, real_t deprcv, bool allowLiquid){ + +MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, real_t deprcv, bool allowLiquid) +{ GRTCheckFileExist(modelpath); if(depsrc * deprcv < 0.0){ @@ -429,9 +296,6 @@ MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, real_t fclose(fp); GRT_SAFE_FREE_PTR(modarr); GRT_SAFE_FREE_PTR(line); - - // 先指定负频率，仅填充弹性模量 - grt_attenuate_mod1d(mod1d, -1); // 设置一个默认边界条件 mod1d->topbound = GRT_BOUND_FREE; @@ -457,28 +321,198 @@ void grt_set_mod1d_boundary(MODEL1D *mod1d, GRT_BOUND_TYPE topbound, GRT_BOUND_T } -void grt_get_model_diglen_from_file(const char *modelpath, size_t diglen[6]){ - FILE *fp = GRTCheckOpenFile(modelpath, "r"); - size_t len; - char *line = NULL; +// ===================================================================================================================== +// ===================================================================================================================== - memset(diglen, 0, sizeof(size_t[6])); +MODEL1D_STATE * grt_init_mod1d_state(MODEL1D *mod1d) +{ + MODEL1D_STATE *mstat = (MODEL1D_STATE *)calloc(1, sizeof(MODEL1D_STATE)); + size_t n = mod1d->n; + mstat->mod1d = mod1d; + + #define X(P, T) mstat->P = (T*)calloc(n, sizeof(T)); + __MODEL1D_STATE_FOR_EACH_ARRAY + #undef X - while(grt_getline(&line, &len, fp) != -1){ - char *token = strtok(line, " \n"); - for(int i=0; i<6; ++i){ - if(token == NULL) break; - diglen[i] = GRT_MAX(diglen[i], strlen(token)); - token = strtok(NULL, " \n"); + return mstat; +} + + +MODEL1D_STATE * grt_copy_mod1d_state(const MODEL1D_STATE *mstat1) +{ + MODEL1D_STATE *mstat2 = (MODEL1D_STATE *)calloc(1, sizeof(MODEL1D_STATE)); + + // 先直接赋值，实现浅拷贝 + *mstat2 = *mstat1; + + // 对指针部分再重新申请内存并赋值，实现深拷贝 + size_t n = mstat1->mod1d->n; + #define X(P, T) \ + mstat2->P = (T*)calloc(n, sizeof(T));\ + memcpy(mstat2->P, mstat1->P, sizeof(T)*n);\ + + __MODEL1D_STATE_FOR_EACH_ARRAY + #undef X + + // 注意！这里没有对 mstat1->mod1d 进行深拷贝，因为始终视作其中为常量，不做变动 + + return mstat2; +} + + +void grt_update_mod1d_state_omega(MODEL1D_STATE *mstat, const cplx_t omega) +{ + mstat->omega = omega; + + MODEL1D *mod1d = mstat->mod1d; + real_t Va0, Vb0; + cplx_t atna, atnb; + + for(size_t i = 0; i < mod1d->n; ++i){ + Va0 = mod1d->Va[i]; + Vb0 = mod1d->Vb[i]; + + // 圆频率实部为负数表明不考虑模型的 Q 值属性 + // 在读入模型后需要需要运行一次本函数以填充弹性模量，见 grt_read_mod1d_from_file 函数 + atna = (creal(omega) >= 0.0 && mod1d->Qainv[i] > 0.0)? grt_attenuation_law(mod1d->Qainv[i], mod1d->omgref, omega) : 1.0; + atnb = (creal(omega) >= 0.0 && mod1d->Qbinv[i] > 0.0)? grt_attenuation_law(mod1d->Qbinv[i], mod1d->omgref, omega) : 1.0; + + mstat->atna[i] = atna; + mstat->atnb[i] = atnb; + + mstat->mu[i] = (Vb0*atnb)*(Vb0*atnb)*(mod1d->Rho[i]); + mstat->lambda[i] = (Va0*atnb)*(Va0*atnb)*(mod1d->Rho[i]) - 2*mstat->mu[i]; + mstat->delta[i] = (mstat->lambda[i] + mstat->mu[i]) / (mstat->lambda[i] + 3.0*mstat->mu[i]); + } +} + + +void grt_update_mod1d_state_k(MODEL1D_STATE *mstat, const real_t k) +{ + MODEL1D *mod1d = mstat->mod1d; + mstat->k = k; + // 不合理的频率值，只可能是在计算静态解，此时不需要xa, xb等物理量 + if(creal(mstat->omega) < 0.0) return; + mstat->c_phase = mstat->omega/k; + + + size_t isrc = mod1d->isrc; + size_t ircv = mod1d->ircv; + + for(size_t i=0; in; ++i){ + if( mod1d->srcrcv_isInserted && (i == isrc || i == ircv) ){ + mstat->xa[i] = mstat->xa[i-1]; + mstat->caca[i] = mstat->caca[i-1]; + mstat->xb[i] = mstat->xb[i-1]; + mstat->cbcb[i] = mstat->cbcb[i-1]; + continue; } + + real_t va, vb; + va = mod1d->Va[i]; + vb = mod1d->Vb[i]; + cplx_t atna, atnb; + atna = mstat->atna[i]; + atnb = mstat->atnb[i]; + + cplx_t caca, cbcb; + caca = mstat->c_phase / (va*atna); + caca *= caca; + mstat->caca[i] = caca; + mstat->xa[i] = sqrt(1.0 - caca); + + cbcb = (mod1d->isLiquid[i])? 0.0 : mstat->c_phase / (vb*atnb); // 考虑液体层 + cbcb *= cbcb; + mstat->cbcb[i] = cbcb; + mstat->xb[i] = sqrt(1.0 - cbcb); } +} - GRT_SAFE_FREE_PTR(line); - fclose(fp); +void grt_free_mod1d_state(MODEL1D_STATE *mstat) +{ + #define X(P, T) GRT_SAFE_FREE_PTR(mstat->P); + __MODEL1D_STATE_FOR_EACH_ARRAY + #undef X + + GRT_SAFE_FREE_PTR(mstat); } -bool grt_check_vel_in_mod(const MODEL1D *mod1d, const real_t vel, const real_t tol){ +// ===================================================================================================================== +// ===================================================================================================================== + +void grt_print_mod1d(const MODEL1D *mod1d) +{ + // 模拟表格，打印速度 + // 每列字符宽度 + // [isrc/ircv] [h(km)] [Vp(km/s)] [Vs(km/s)] [Rho(g/cm^3)] [Qp] [Qs] + const int ncols = 7; + const int nlens[] = {13, 12, 13, 13, 16, 13, 13}; + int Nlen=0; + for(int ic=0; icn; ++i){ + if(i==mod1d->isrc){ + snprintf(indexstr, sizeof(indexstr), "%zu [src]", i+1); + } else if(i==mod1d->ircv){ + snprintf(indexstr, sizeof(indexstr), "%zu [rcv]", i+1); + } else { + snprintf(indexstr, sizeof(indexstr), "%zu ", i+1); + } + + printf("| %*s ", nlens[0]-3, indexstr); + + if(i < mod1d->n-1){ + printf("| %-*.2f ", nlens[1]-3, mod1d->Thk[i]); + } else { + printf("| %-*s ", nlens[1]-3, "Inf"); + } + + printf("| %-*.2f ", nlens[2]-3, mod1d->Va[i]); + printf("| %-*.2f ", nlens[3]-3, mod1d->Vb[i]); + printf("| %-*.2f ", nlens[4]-3, mod1d->Rho[i]); + printf("| %-*.2e ", nlens[5]-3, mod1d->Qa[i]); + printf("| %-*.2e ", nlens[6]-3, mod1d->Qb[i]); + printf("|\n"); + } + printf("%s\n", splitline); + printf("\n"); +} + + + +bool grt_check_vel_in_mod(const MODEL1D *mod1d, const real_t vel, const real_t tol) +{ // 浮点数比较，检查是否存在该速度值 for(size_t i=0; in; ++i){ if(fabs(vel - mod1d->Va[i])Vb[i])Va; @@ -499,4 +533,4 @@ void grt_get_mod1d_vmin_vmax(const MODEL1D *mod1d, real_t *vmin, real_t *vmax){ 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/dynamic/grn.c b/pygrt/C_extension/src/dynamic/grn.c index 945359fe..aaabbceb 100755 --- a/pygrt/C_extension/src/dynamic/grn.c +++ b/pygrt/C_extension/src/dynamic/grn.c @@ -24,7 +24,6 @@ #include "grt/common/const.h" #include "grt/common/model.h" -#include "grt/common/prtdbg.h" #include "grt/common/search.h" #include "grt/common/progressbar.h" @@ -118,22 +117,19 @@ void grt_integ_grn_spec( cplx_t coef = - dk*fac / GRT_SQUARE(omega); // 最终要乘上的系数 - MODEL1D *local_mod1d = NULL; K_INTEG_METHOD *local_Kmet = NULL; #ifdef _OPENMP - // 定义局部模型对象 - local_mod1d = grt_copy_mod1d(mod1d); + // 定义局部对象 K_INTEG_METHOD __KMET = *Kmet; // 只需浅拷贝 local_Kmet = &__KMET; #else - local_mod1d = mod1d; local_Kmet = Kmet; #endif - // 将 omega 计入模型结构体 - local_mod1d->omega = omega; + MODEL1D_STATE *local_mstat = grt_init_mod1d_state(mod1d); - grt_attenuate_mod1d(local_mod1d, omega); + // 将 omega 计入模型结构体 + grt_update_mod1d_state_omega(local_mstat, omega); // 是否要输出积分过程文件 bool needfstats = (statsstr!=NULL && (grt_findElement_size_t(statsidxs, nstatsidxs, iw) >= 0)); @@ -154,7 +150,7 @@ void grt_integ_grn_spec( // Wavenumber Integration // 波数积分上限 local_Kmet->kmax = hypot(local_Kmet->k0, local_Kmet->ampk * w / local_Kmet->vmin); - K_INTEG *Kint = grt_wavenumber_integral(local_mod1d, nr, rs, local_Kmet, calc_upar, grt_kernel); + K_INTEG *Kint = grt_wavenumber_integral(local_mstat, nr, rs, local_Kmet, calc_upar, grt_kernel); // 记录到格林函数结构体内 // 如果计算核函数过程中存在除零错误，则放弃该频率【通常在大震中距的低频段】 @@ -165,13 +161,11 @@ void grt_integ_grn_spec( } // =================================================================================== - freq_invstats[iw] = local_mod1d->stats; + freq_invstats[iw] = local_mstat->stats; if(needfstats) grt_KMET_destroy_fstats(nr, local_Kmet); - #ifdef _OPENMP - grt_free_mod1d(local_mod1d); - #endif + grt_free_mod1d_state(local_mstat); // 记录进度条变量 #pragma omp critical diff --git a/pygrt/C_extension/src/dynamic/grt_kernel.c b/pygrt/C_extension/src/dynamic/grt_kernel.c index 75f4aa76..923e3d98 100644 --- a/pygrt/C_extension/src/dynamic/grt_kernel.c +++ b/pygrt/C_extension/src/dynamic/grt_kernel.c @@ -430,18 +430,12 @@ int kernel_main(int argc, char **argv) real_t w = PI2*freq; cplx_t omega = w - I*Ctrl->F.wI; - MODEL1D *local_mod1d = NULL; - #ifdef _OPENMP - // 定义局部模型对象 - local_mod1d = grt_copy_mod1d(mod1d); - #else - local_mod1d = mod1d; - #endif + MODEL1D_STATE *local_mstat = grt_init_mod1d_state(mod1d); // 将 omega 计入模型结构体 - local_mod1d->omega = omega; + local_mstat->omega = omega; - grt_attenuate_mod1d(local_mod1d, omega); + grt_update_mod1d_state_omega(local_mstat, omega); // 为当前频率创建波数积分记录文件 FILE *fstats = NULL; @@ -459,7 +453,7 @@ int kernel_main(int argc, char **argv) { real_t cc = Ctrl->C.c_phases[ic]; real_t k = w/cc; - grt_kernel(local_mod1d, k, QWV, Ctrl->e.active, QWV_uiz); + grt_kernel(local_mstat, k, QWV, Ctrl->e.active, QWV_uiz); // 系数 GRT_LOOP_ChnlGrid(im, c){ @@ -476,9 +470,7 @@ int kernel_main(int argc, char **argv) fclose(fstats); - #ifdef _OPENMP - grt_free_mod1d(local_mod1d); - #endif + grt_free_mod1d_state(local_mstat); } diff --git a/pygrt/C_extension/src/dynamic/layer.c b/pygrt/C_extension/src/dynamic/layer.c index 5bc90211..41838f08 100644 --- a/pygrt/C_extension/src/dynamic/layer.c +++ b/pygrt/C_extension/src/dynamic/layer.c @@ -18,15 +18,14 @@ #include "grt/dynamic/layer.h" #include "grt/common/model.h" -#include "grt/common/prtdbg.h" #include "grt/common/matrix.h" #include "grt/common/checkerror.h" /* 定义用于提取相邻两层物性参数的宏 */ -#define MODEL_2LAYS_ATTRIB(T, N) \ - T N##1 = mod1d->N[iy-1];\ - T N##2 = mod1d->N[iy];\ +#define MODEL_2LAYS_ATTRIB(MM, T, N) \ + T N##1 = MM->N[iy-1];\ + T N##2 = MM->N[iy];\ static int __freeBound_R_PSV(bool isLiquid, cplx_t cbcb0, cplx_t xa0, cplx_t xb0, real_t adsgn, cplx_t M[2][2]) @@ -106,21 +105,22 @@ static int __rigidBound_R_SH(bool isLiquid, cplx_t *ML) return GRT_INVERSE_SUCCESS; } -void grt_topbound_RU_PSV(MODEL1D *mod1d) +void grt_topbound_RU_PSV(MODEL1D_STATE *mstat) { - cplx_t cbcb = mod1d->cbcb[0]; - cplx_t xa = mod1d->xa[0]; - cplx_t xb = mod1d->xb[0]; + cplx_t cbcb = mstat->cbcb[0]; + cplx_t xa = mstat->xa[0]; + cplx_t xb = mstat->xb[0]; + MODEL1D *mod1d = mstat->mod1d; bool isLiquid = mod1d->isLiquid[0]; if(mod1d->topbound == GRT_BOUND_FREE){ - mod1d->M_top.stats = __freeBound_R_PSV(isLiquid, cbcb, xa, xb, 1.0, mod1d->M_top.RU); + mstat->M_top.stats = __freeBound_R_PSV(isLiquid, cbcb, xa, xb, 1.0, mstat->M_top.RU); } else if(mod1d->topbound == GRT_BOUND_RIGID){ - mod1d->M_top.stats = __rigidBound_R_PSV(isLiquid, xa, xb, 1.0, mod1d->M_top.RU); + mstat->M_top.stats = __rigidBound_R_PSV(isLiquid, xa, xb, 1.0, mstat->M_top.RU); } else if(mod1d->topbound == GRT_BOUND_HALFSPACE){ - memset(mod1d->M_top.RU, 0, sizeof(cplx_t)*4); - mod1d->M_top.stats = GRT_INVERSE_SUCCESS; + memset(mstat->M_top.RU, 0, sizeof(cplx_t)*4); + mstat->M_top.stats = GRT_INVERSE_SUCCESS; } else{ GRTRaiseError("Wrong execution."); @@ -128,18 +128,19 @@ void grt_topbound_RU_PSV(MODEL1D *mod1d) // RU 不需要时延 } -void grt_topbound_RU_SH(MODEL1D *mod1d) +void grt_topbound_RU_SH(MODEL1D_STATE *mstat) { + MODEL1D *mod1d = mstat->mod1d; bool isLiquid = mod1d->isLiquid[0]; if(mod1d->topbound == GRT_BOUND_FREE){ - mod1d->M_top.stats = __freeBound_R_SH(isLiquid, &mod1d->M_top.RUL); + mstat->M_top.stats = __freeBound_R_SH(isLiquid, &mstat->M_top.RUL); } else if(mod1d->topbound == GRT_BOUND_RIGID){ - mod1d->M_top.stats = __rigidBound_R_SH(isLiquid, &mod1d->M_top.RUL); + mstat->M_top.stats = __rigidBound_R_SH(isLiquid, &mstat->M_top.RUL); } else if(mod1d->topbound == GRT_BOUND_HALFSPACE){ - mod1d->M_top.RUL = 0.0; - mod1d->M_top.stats = GRT_INVERSE_SUCCESS; + mstat->M_top.RUL = 0.0; + mstat->M_top.stats = GRT_INVERSE_SUCCESS; } else{ GRTRaiseError("Wrong execution."); @@ -147,28 +148,29 @@ void grt_topbound_RU_SH(MODEL1D *mod1d) // RU 不需要时延 } -void grt_botbound_RD_PSV(MODEL1D *mod1d) +void grt_botbound_RD_PSV(MODEL1D_STATE *mstat) { + MODEL1D *mod1d = mstat->mod1d; size_t nlay = mod1d->n; - cplx_t cbcb = mod1d->cbcb[nlay-2]; - cplx_t xa = mod1d->xa[nlay-2]; - cplx_t xb = mod1d->xb[nlay-2]; + cplx_t cbcb = mstat->cbcb[nlay-2]; + cplx_t xa = mstat->xa[nlay-2]; + cplx_t xb = mstat->xb[nlay-2]; bool isLiquid = mod1d->isLiquid[nlay-2]; if(mod1d->botbound == GRT_BOUND_FREE){ - mod1d->M_bot.stats = __freeBound_R_PSV(isLiquid, cbcb, xa, xb, -1.0, mod1d->M_bot.RD); + mstat->M_bot.stats = __freeBound_R_PSV(isLiquid, cbcb, xa, xb, -1.0, mstat->M_bot.RD); } else if(mod1d->botbound == GRT_BOUND_RIGID){ - mod1d->M_bot.stats = __rigidBound_R_PSV(isLiquid, xa, xb, -1.0, mod1d->M_bot.RD); + mstat->M_bot.stats = __rigidBound_R_PSV(isLiquid, xa, xb, -1.0, mstat->M_bot.RD); } else if(mod1d->botbound == GRT_BOUND_HALFSPACE){ - grt_RT_matrix_PSV(mod1d, nlay-1, &mod1d->M_bot); + grt_RT_matrix_PSV(mstat, nlay-1, &mstat->M_bot); } else{ GRTRaiseError("Wrong execution."); } // 时延 RD - real_t k = mod1d->k; + real_t k = mstat->k; real_t thk = mod1d->Thk[nlay-2]; cplx_t exa, exb, ex2a, ex2b, exab; @@ -178,47 +180,49 @@ void grt_botbound_RD_PSV(MODEL1D *mod1d) ex2b = exb * exb; exab = exa * exb; - mod1d->M_bot.RD[0][0] *= ex2a; mod1d->M_bot.RD[0][1] *= exab; - mod1d->M_bot.RD[1][0] *= exab; mod1d->M_bot.RD[1][1] *= ex2b; + mstat->M_bot.RD[0][0] *= ex2a; mstat->M_bot.RD[0][1] *= exab; + mstat->M_bot.RD[1][0] *= exab; mstat->M_bot.RD[1][1] *= ex2b; } -void grt_botbound_RD_SH(MODEL1D *mod1d) +void grt_botbound_RD_SH(MODEL1D_STATE *mstat) { + MODEL1D *mod1d = mstat->mod1d; size_t nlay = mod1d->n; - cplx_t xb = mod1d->xb[nlay-2]; + cplx_t xb = mstat->xb[nlay-2]; bool isLiquid = mod1d->isLiquid[nlay-2]; if(mod1d->botbound == GRT_BOUND_FREE){ - mod1d->M_bot.stats = __freeBound_R_SH(isLiquid, &mod1d->M_bot.RDL); + mstat->M_bot.stats = __freeBound_R_SH(isLiquid, &mstat->M_bot.RDL); } else if(mod1d->botbound == GRT_BOUND_RIGID){ - mod1d->M_bot.stats = __rigidBound_R_SH(isLiquid, &mod1d->M_bot.RDL); + mstat->M_bot.stats = __rigidBound_R_SH(isLiquid, &mstat->M_bot.RDL); } else if(mod1d->botbound == GRT_BOUND_HALFSPACE){ - grt_RT_matrix_SH(mod1d, nlay-1, &mod1d->M_bot); + grt_RT_matrix_SH(mstat, nlay-1, &mstat->M_bot); } else{ GRTRaiseError("Wrong execution."); } // 时延 RD - real_t k = mod1d->k; + real_t k = mstat->k; real_t thk = mod1d->Thk[nlay-2]; cplx_t exb, ex2b; exb = exp(- k*thk*xb); ex2b = exb * exb; - mod1d->M_bot.RDL *= ex2b; + mstat->M_bot.RDL *= ex2b; } -void grt_wave2qwv_REV_PSV(MODEL1D *mod1d) +void grt_wave2qwv_REV_PSV(MODEL1D_STATE *mstat) { + MODEL1D *mod1d = mstat->mod1d; size_t ircv = mod1d->ircv; - cplx_t xa = mod1d->xa[ircv]; - cplx_t xb = mod1d->xb[ircv]; + cplx_t xa = mstat->xa[ircv]; + cplx_t xb = mstat->xb[ircv]; bool isLiquid = mod1d->isLiquid[ircv]; - real_t k = mod1d->k; + real_t k = mstat->k; cplx_t D11[2][2], D12[2][2]; if( ! isLiquid){ @@ -238,43 +242,45 @@ void grt_wave2qwv_REV_PSV(MODEL1D *mod1d) // 公式(5.7.7,25) if(mod1d->ircvup){// 震源更深 - grt_cmat2x2_mul(D12, mod1d->M_FA.RU, mod1d->R_EV); - grt_cmat2x2_add(D11, mod1d->R_EV, mod1d->R_EV); + grt_cmat2x2_mul(D12, mstat->M_FA.RU, mstat->R_EV); + grt_cmat2x2_add(D11, mstat->R_EV, mstat->R_EV); } else { // 接收点更深 - grt_cmat2x2_mul(D11, mod1d->M_BL.RD, mod1d->R_EV); - grt_cmat2x2_add(D12, mod1d->R_EV, mod1d->R_EV); + grt_cmat2x2_mul(D11, mstat->M_BL.RD, mstat->R_EV); + grt_cmat2x2_add(D12, mstat->R_EV, mstat->R_EV); } } -void grt_wave2qwv_REV_SH(MODEL1D *mod1d) +void grt_wave2qwv_REV_SH(MODEL1D_STATE *mstat) { + MODEL1D *mod1d = mstat->mod1d; size_t ircv = mod1d->ircv; bool isLiquid = mod1d->isLiquid[ircv]; - real_t k = mod1d->k; + real_t k = mstat->k; if( ! isLiquid){ // 位于固体层 // 公式(5.2.19) if(mod1d->ircvup){// 震源更深 - mod1d->R_EVL = (1.0 + mod1d->M_FA.RUL)*k; + mstat->R_EVL = (1.0 + mstat->M_FA.RUL)*k; } else { - mod1d->R_EVL = (1.0 + mod1d->M_BL.RDL)*k; + mstat->R_EVL = (1.0 + mstat->M_BL.RDL)*k; } } else { // 位于液体层 - mod1d->R_EVL = 0.0; + mstat->R_EVL = 0.0; } } -void grt_wave2qwv_z_REV_PSV(MODEL1D *mod1d) +void grt_wave2qwv_z_REV_PSV(MODEL1D_STATE *mstat) { + MODEL1D *mod1d = mstat->mod1d; size_t ircv = mod1d->ircv; - cplx_t xa = mod1d->xa[ircv]; - cplx_t xb = mod1d->xb[ircv]; + cplx_t xa = mstat->xa[ircv]; + cplx_t xb = mstat->xb[ircv]; bool isLiquid = mod1d->isLiquid[ircv]; - real_t k = mod1d->k; + real_t k = mstat->k; // 将垂直波函数转为ui,z在(B_m, P_m, C_m)系下的分量 @@ -293,21 +299,22 @@ void grt_wave2qwv_z_REV_PSV(MODEL1D *mod1d) // 公式(5.7.7,25) if(mod1d->ircvup){// 震源更深 - grt_cmat2x2_mul(D12, mod1d->M_FA.RU, mod1d->uiz_R_EV); - grt_cmat2x2_add(D11, mod1d->uiz_R_EV, mod1d->uiz_R_EV); + grt_cmat2x2_mul(D12, mstat->M_FA.RU, mstat->uiz_R_EV); + grt_cmat2x2_add(D11, mstat->uiz_R_EV, mstat->uiz_R_EV); } else { // 接收点更深 - grt_cmat2x2_mul(D11, mod1d->M_BL.RD, mod1d->uiz_R_EV); - grt_cmat2x2_add(D12, mod1d->uiz_R_EV, mod1d->uiz_R_EV); + grt_cmat2x2_mul(D11, mstat->M_BL.RD, mstat->uiz_R_EV); + grt_cmat2x2_add(D12, mstat->uiz_R_EV, mstat->uiz_R_EV); } } -void grt_wave2qwv_z_REV_SH(MODEL1D *mod1d) +void grt_wave2qwv_z_REV_SH(MODEL1D_STATE *mstat) { + MODEL1D *mod1d = mstat->mod1d; size_t ircv = mod1d->ircv; - cplx_t xb = mod1d->xb[ircv]; + cplx_t xb = mstat->xb[ircv]; bool isLiquid = mod1d->isLiquid[ircv]; - real_t k = mod1d->k; + real_t k = mstat->k; // 将垂直波函数转为ui,z在(B_m, P_m, C_m)系下的分量 // 新推导的公式 @@ -316,21 +323,21 @@ void grt_wave2qwv_z_REV_SH(MODEL1D *mod1d) if( ! isLiquid){ // 位于固体层 if(mod1d->ircvup){// 震源更深 - mod1d->uiz_R_EVL = (1.0 - mod1d->M_FA.RUL)*bk; + mstat->uiz_R_EVL = (1.0 - mstat->M_FA.RUL)*bk; } else { // 接收点更深 - mod1d->uiz_R_EVL = (mod1d->M_BL.RDL - 1.0)*bk; + mstat->uiz_R_EVL = (mstat->M_BL.RDL - 1.0)*bk; } } else { // 位于液体层 - mod1d->uiz_R_EVL = 0.0; + mstat->uiz_R_EVL = 0.0; } } -void grt_RT_matrix_ll_PSV(const MODEL1D *mod1d, size_t iy, RT_MATRIX *M) +void grt_RT_matrix_ll_PSV(const MODEL1D_STATE *mstat, size_t iy, RT_MATRIX *M) { - MODEL_2LAYS_ATTRIB(cplx_t, xa); - MODEL_2LAYS_ATTRIB(real_t, Rho); + MODEL_2LAYS_ATTRIB(mstat, cplx_t, xa); + MODEL_2LAYS_ATTRIB(mstat->mod1d, real_t, Rho); cplx_t A = xa1*Rho2 + xa2*Rho1; if(A==0.0){ @@ -362,14 +369,14 @@ void grt_RT_matrix_ll_SH(RT_MATRIX *M) -void grt_RT_matrix_ls_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_RT_matrix_ls_PSV(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M) { - MODEL_2LAYS_ATTRIB(cplx_t, xa); - MODEL_2LAYS_ATTRIB(cplx_t, xb); - MODEL_2LAYS_ATTRIB(cplx_t, mu); - MODEL_2LAYS_ATTRIB(cplx_t, cbcb); - MODEL_2LAYS_ATTRIB(real_t, Rho); - MODEL_2LAYS_ATTRIB(bool, isLiquid); + MODEL_2LAYS_ATTRIB(mstat, cplx_t, xa); + MODEL_2LAYS_ATTRIB(mstat, cplx_t, xb); + MODEL_2LAYS_ATTRIB(mstat, cplx_t, mu); + MODEL_2LAYS_ATTRIB(mstat, cplx_t, cbcb); + MODEL_2LAYS_ATTRIB(mstat->mod1d, real_t, Rho); + MODEL_2LAYS_ATTRIB(mstat->mod1d, bool, isLiquid); // 后缀1表示上层的液体的物理参数，后缀2表示下层的固体的物理参数 // 若mu2==0, 则下层为液体，参数需相互交换 @@ -399,7 +406,7 @@ void grt_RT_matrix_ls_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) // 定义一些中间变量来简化运算和书写 - cplx_t lamka1k = Rho1*GRT_SQUARE(mod1d->c_phase); + cplx_t lamka1k = Rho1*GRT_SQUARE(mstat->c_phase); cplx_t kb2k = cbcb2; cplx_t Og2k = 1.0 - 0.5*kb2k; cplx_t Og2k2 = Og2k*Og2k; @@ -434,9 +441,9 @@ void grt_RT_matrix_ls_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) } -void grt_RT_matrix_ls_SH(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_RT_matrix_ls_SH(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M) { - MODEL_2LAYS_ATTRIB(bool, isLiquid); + MODEL_2LAYS_ATTRIB(mstat->mod1d, bool, isLiquid); // 后缀1表示上层的液体的物理参数，后缀2表示下层的固体的物理参数 // 若mu2==0, 则下层为液体，参数需相互交换 @@ -465,13 +472,13 @@ void grt_RT_matrix_ls_SH(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) -void grt_RT_matrix_ss_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_RT_matrix_ss_PSV(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M) { - MODEL_2LAYS_ATTRIB(cplx_t, xa); - MODEL_2LAYS_ATTRIB(cplx_t, xb); - MODEL_2LAYS_ATTRIB(cplx_t, mu); - MODEL_2LAYS_ATTRIB(cplx_t, cbcb); - MODEL_2LAYS_ATTRIB(real_t, Rho); + MODEL_2LAYS_ATTRIB(mstat, cplx_t, xa); + MODEL_2LAYS_ATTRIB(mstat, cplx_t, xb); + MODEL_2LAYS_ATTRIB(mstat, cplx_t, mu); + MODEL_2LAYS_ATTRIB(mstat, cplx_t, cbcb); + MODEL_2LAYS_ATTRIB(mstat->mod1d, real_t, Rho); // 定义一些中间变量来简化运算和书写 cplx_t rmu = mu1/mu2; @@ -551,10 +558,10 @@ void grt_RT_matrix_ss_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) } -void grt_RT_matrix_ss_SH(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_RT_matrix_ss_SH(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M) { - MODEL_2LAYS_ATTRIB(cplx_t, xb); - MODEL_2LAYS_ATTRIB(cplx_t, mu); + MODEL_2LAYS_ATTRIB(mstat, cplx_t, xb); + MODEL_2LAYS_ATTRIB(mstat, cplx_t, mu); // REFELCTION M->RUL = (mu2*xb2 - mu1*xb1) / (mu2*xb2 + mu1*xb1) ; @@ -569,46 +576,47 @@ void grt_RT_matrix_ss_SH(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) -void grt_RT_matrix_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_RT_matrix_PSV(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M) { - MODEL_2LAYS_ATTRIB(bool, isLiquid); + MODEL_2LAYS_ATTRIB(mstat->mod1d, bool, isLiquid); // 根据界面两侧的具体情况选择函数 if(!isLiquid1 && !isLiquid2){ - grt_RT_matrix_ss_PSV(mod1d, iy, M); + grt_RT_matrix_ss_PSV(mstat, iy, M); } else if(isLiquid1 && isLiquid2){ - grt_RT_matrix_ll_PSV(mod1d, iy, M); + grt_RT_matrix_ll_PSV(mstat, iy, M); } else{ - grt_RT_matrix_ls_PSV(mod1d, iy, M); + grt_RT_matrix_ls_PSV(mstat, iy, M); } } -void grt_RT_matrix_SH(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_RT_matrix_SH(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M) { - MODEL_2LAYS_ATTRIB(bool, isLiquid); + MODEL_2LAYS_ATTRIB(mstat->mod1d, bool, isLiquid); // 根据界面两侧的具体情况选择函数 if(!isLiquid1 && !isLiquid2){ - grt_RT_matrix_ss_SH(mod1d, iy, M); + grt_RT_matrix_ss_SH(mstat, iy, M); } else if(isLiquid1 && isLiquid2){ grt_RT_matrix_ll_SH(M); } else{ - grt_RT_matrix_ls_SH(mod1d, iy, M); + grt_RT_matrix_ls_SH(mstat, iy, M); } } -void grt_delay_RT_matrix(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_delay_RT_matrix(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M) { + MODEL1D *mod1d = mstat->mod1d; real_t thk = mod1d->Thk[iy-1]; - cplx_t xa1 = mod1d->xa[iy-1]; - cplx_t xb1 = mod1d->xb[iy-1]; - real_t k = mod1d->k; + cplx_t xa1 = mstat->xa[iy-1]; + cplx_t xb1 = mstat->xb[iy-1]; + real_t k = mstat->k; cplx_t exa, exb, ex2a, ex2b, exab; exa = exp(- k*thk*xa1); @@ -633,12 +641,13 @@ void grt_delay_RT_matrix(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) M->TUL *= exb; } -void grt_delay_RT_matrix_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_delay_RT_matrix_PSV(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M) { + MODEL1D *mod1d = mstat->mod1d; real_t thk = mod1d->Thk[iy-1]; - cplx_t xa1 = mod1d->xa[iy-1]; - cplx_t xb1 = mod1d->xb[iy-1]; - real_t k = mod1d->k; + cplx_t xa1 = mstat->xa[iy-1]; + cplx_t xb1 = mstat->xb[iy-1]; + real_t k = mstat->k; cplx_t exa, exb, ex2a, ex2b, exab; exa = exp(- k*thk*xa1); @@ -657,11 +666,12 @@ void grt_delay_RT_matrix_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M M->TU[1][0] *= exb; M->TU[1][1] *= exb; } -void grt_delay_RT_matrix_SH(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_delay_RT_matrix_SH(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M) { + MODEL1D *mod1d = mstat->mod1d; real_t thk = mod1d->Thk[iy-1]; - cplx_t xb1 = mod1d->xb[iy-1]; - real_t k = mod1d->k; + cplx_t xb1 = mstat->xb[iy-1]; + real_t k = mstat->k; cplx_t exb, ex2b; exb = exp(- k*thk*xb1); @@ -674,12 +684,13 @@ void grt_delay_RT_matrix_SH(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) -void grt_delay_GRT_matrix(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_delay_GRT_matrix(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M) { + MODEL1D *mod1d = mstat->mod1d; real_t thk = mod1d->Thk[iy-1]; - cplx_t xa1 = mod1d->xa[iy-1]; - cplx_t xb1 = mod1d->xb[iy-1]; - real_t k = mod1d->k; + cplx_t xa1 = mstat->xa[iy-1]; + cplx_t xb1 = mstat->xb[iy-1]; + real_t k = mstat->k; cplx_t exa, exb, ex2a, ex2b, exab; exa = exp(- k*thk*xa1); @@ -703,18 +714,18 @@ void grt_delay_GRT_matrix(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) } -void grt_get_layer_D(const MODEL1D *mod1d, const size_t iy, bool inverse, int liquid_invtype, cplx_t D[4][4]) +void grt_get_layer_D(const MODEL1D_STATE *mstat, const size_t iy, bool inverse, int liquid_invtype, cplx_t D[4][4]) { // 第iy层物理量 - real_t k = mod1d->k; + real_t k = mstat->k; real_t kk = k*k; - cplx_t xa = mod1d->xa[iy]; - cplx_t xb = mod1d->xb[iy]; - cplx_t mu = mod1d->mu[iy]; - cplx_t cbcb = mod1d->cbcb[iy]; - real_t rho = mod1d->Rho[iy]; + cplx_t xa = mstat->xa[iy]; + cplx_t xb = mstat->xb[iy]; + cplx_t mu = mstat->mu[iy]; + cplx_t cbcb = mstat->cbcb[iy]; + real_t rho = mstat->mod1d->Rho[iy]; cplx_t Omgn, Omg; - if(! mod1d->isLiquid[iy]){ + if(! mstat->mod1d->isLiquid[iy]){ Omgn = 1.0 - 0.5*cbcb; Omg = Omgn*kk; if( ! inverse ){ @@ -734,7 +745,7 @@ void grt_get_layer_D(const MODEL1D *mod1d, const size_t iy, bool inverse, int li } } } else { - Omg = rho * GRT_SQUARE(mod1d->omega) / k; + Omg = rho * GRT_SQUARE(mstat->omega) / k; if(liquid_invtype == 1){ if( ! inverse ){ D[0][0] = k; D[0][1] = 0.0; D[0][2] = k; D[0][3] = 0.0; @@ -774,12 +785,12 @@ void grt_get_layer_D(const MODEL1D *mod1d, const size_t iy, bool inverse, int li } -void grt_get_layer_D11(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]) +void grt_get_layer_D11(const MODEL1D_STATE *mstat, const size_t iy, cplx_t D[2][2]) { - real_t k = mod1d->k; - cplx_t xa = mod1d->xa[iy]; - cplx_t xb = mod1d->xb[iy]; - if(! mod1d->isLiquid[iy]){ + real_t k = mstat->k; + cplx_t xa = mstat->xa[iy]; + cplx_t xb = mstat->xb[iy]; + if(! mstat->mod1d->isLiquid[iy]){ D[0][0] = k; D[0][1] = k*xb; D[1][0] = k*xa; D[1][1] = k; } else { @@ -789,12 +800,12 @@ void grt_get_layer_D11(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]) } -void grt_get_layer_D12(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]) +void grt_get_layer_D12(const MODEL1D_STATE *mstat, const size_t iy, cplx_t D[2][2]) { - real_t k = mod1d->k; - cplx_t xa = mod1d->xa[iy]; - cplx_t xb = mod1d->xb[iy]; - if(! mod1d->isLiquid[iy]){ + real_t k = mstat->k; + cplx_t xa = mstat->xa[iy]; + cplx_t xb = mstat->xb[iy]; + if(! mstat->mod1d->isLiquid[iy]){ D[0][0] = k; D[0][1] = -k*xb; D[1][0] = -k*xa; D[1][1] = k; } else { @@ -803,15 +814,15 @@ void grt_get_layer_D12(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]) } } -void grt_get_layer_D11_uiz(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]) +void grt_get_layer_D11_uiz(const MODEL1D_STATE *mstat, const size_t iy, cplx_t D[2][2]) { - real_t k = mod1d->k; - cplx_t xa = mod1d->xa[iy]; - cplx_t xb = mod1d->xb[iy]; + real_t k = mstat->k; + cplx_t xa = mstat->xa[iy]; + cplx_t xb = mstat->xb[iy]; cplx_t a = k*xa; cplx_t b = k*xb; - if(! mod1d->isLiquid[iy]){ + if(! mstat->mod1d->isLiquid[iy]){ D[0][0] = a*k; D[0][1] = b*b; D[1][0] = a*a; D[1][1] = b*k; } else { @@ -820,15 +831,15 @@ void grt_get_layer_D11_uiz(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2] } } -void grt_get_layer_D12_uiz(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]) +void grt_get_layer_D12_uiz(const MODEL1D_STATE *mstat, const size_t iy, cplx_t D[2][2]) { - real_t k = mod1d->k; - cplx_t xa = mod1d->xa[iy]; - cplx_t xb = mod1d->xb[iy]; + real_t k = mstat->k; + cplx_t xa = mstat->xa[iy]; + cplx_t xb = mstat->xb[iy]; cplx_t a = k*xa; cplx_t b = k*xb; - if(! mod1d->isLiquid[iy]){ + if(! mstat->mod1d->isLiquid[iy]){ D[0][0] = - a*k; D[0][1] = b*b; D[1][0] = a*a; D[1][1] = - b*k; } else { @@ -837,55 +848,55 @@ void grt_get_layer_D12_uiz(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2] } } -void grt_get_layer_D21(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]) +void grt_get_layer_D21(const MODEL1D_STATE *mstat, const size_t iy, cplx_t D[2][2]) { - real_t k = mod1d->k; - cplx_t xa = mod1d->xa[iy]; - cplx_t xb = mod1d->xb[iy]; - cplx_t mu = mod1d->mu[iy]; - cplx_t cbcb = mod1d->cbcb[iy]; - real_t rho = mod1d->Rho[iy]; + real_t k = mstat->k; + cplx_t xa = mstat->xa[iy]; + cplx_t xb = mstat->xb[iy]; + cplx_t mu = mstat->mu[iy]; + cplx_t cbcb = mstat->cbcb[iy]; + real_t rho = mstat->mod1d->Rho[iy]; cplx_t Omg; - if(! mod1d->isLiquid[iy]){ + if(! mstat->mod1d->isLiquid[iy]){ Omg = k*k*(1.0 - 0.5*cbcb); D[0][0] = 2*mu*Omg; D[0][1] = 2*k*mu*k*xb; D[1][0] = 2*k*mu*k*xa; D[1][1] = 2*mu*Omg; } else { - D[0][0] = - rho * GRT_SQUARE(mod1d->omega); D[0][1] = 0.0; + D[0][0] = - rho * GRT_SQUARE(mstat->omega); D[0][1] = 0.0; D[1][0] = 0.0; D[1][1] = 0.0; } } -void grt_get_layer_D22(const MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]) +void grt_get_layer_D22(const MODEL1D_STATE *mstat, const size_t iy, cplx_t D[2][2]) { - real_t k = mod1d->k; - cplx_t xa = mod1d->xa[iy]; - cplx_t xb = mod1d->xb[iy]; - cplx_t mu = mod1d->mu[iy]; - cplx_t cbcb = mod1d->cbcb[iy]; - real_t rho = mod1d->Rho[iy]; + real_t k = mstat->k; + cplx_t xa = mstat->xa[iy]; + cplx_t xb = mstat->xb[iy]; + cplx_t mu = mstat->mu[iy]; + cplx_t cbcb = mstat->cbcb[iy]; + real_t rho = mstat->mod1d->Rho[iy]; cplx_t Omg; - if(! mod1d->isLiquid[iy]){ + if(! mstat->mod1d->isLiquid[iy]){ Omg = k*k*(1.0 - 0.5*cbcb); D[0][0] = 2*mu*Omg; D[0][1] = -2*k*mu*k*xb; D[1][0] = -2*k*mu*k*xa; D[1][1] = 2*mu*Omg; } else { - D[0][0] = - rho * GRT_SQUARE(mod1d->omega); D[0][1] = 0.0; + D[0][0] = - rho * GRT_SQUARE(mstat->omega); D[0][1] = 0.0; D[1][0] = 0.0; D[1][1] = 0.0; } } -void grt_get_layer_T(const MODEL1D *mod1d, const size_t iy, bool inverse, cplx_t T[2][2]) +void grt_get_layer_T(const MODEL1D_STATE *mstat, const size_t iy, bool inverse, cplx_t T[2][2]) { // 液体层不应该使用该函数 - if(mod1d->isLiquid[iy]){ + if(mstat->mod1d->isLiquid[iy]){ GRTRaiseError("Wrong execution."); } - real_t k = mod1d->k; - cplx_t xb = mod1d->xb[iy]; - cplx_t mu = mod1d->mu[iy]; + real_t k = mstat->k; + cplx_t xb = mstat->xb[iy]; + cplx_t mu = mstat->mu[iy]; if( ! inverse ){ T[0][0] = k; T[0][1] = k; diff --git a/pygrt/C_extension/src/dynamic/source.c b/pygrt/C_extension/src/dynamic/source.c index a3026f52..ec9eca48 100755 --- a/pygrt/C_extension/src/dynamic/source.c +++ b/pygrt/C_extension/src/dynamic/source.c @@ -63,36 +63,36 @@ inline GCC_ALWAYS_INLINE void _source_SH(const cplx_t xb, const cplx_t cbcb, con } -void grt_source_coef(const MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU) +void grt_source_coef(const MODEL1D_STATE *mstat, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU) { // 先全部赋0 memset(src_coefD, 0, sizeof(cplxChnlGrid)); memset(src_coefU, 0, sizeof(cplxChnlGrid)); - grt_source_coef_PSV(mod1d, src_coefD, src_coefU); - grt_source_coef_SH(mod1d, src_coefD, src_coefU); + grt_source_coef_PSV(mstat, src_coefD, src_coefU); + grt_source_coef_SH(mstat, src_coefD, src_coefU); } -void grt_source_coef_PSV(const MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU) +void grt_source_coef_PSV(const MODEL1D_STATE *mstat, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU) { - size_t isrc = mod1d->isrc; - cplx_t xa = mod1d->xa[isrc]; - cplx_t caca = mod1d->caca[isrc]; - cplx_t xb = mod1d->xb[isrc]; - cplx_t cbcb = mod1d->cbcb[isrc]; - real_t k = mod1d->k; + size_t isrc = mstat->mod1d->isrc; + cplx_t xa = mstat->xa[isrc]; + cplx_t caca = mstat->caca[isrc]; + cplx_t xb = mstat->xb[isrc]; + cplx_t cbcb = mstat->cbcb[isrc]; + real_t k = mstat->k; _source_PSV(xa, caca, xb, cbcb, k, src_coefD, src_coefU); } -void grt_source_coef_SH(const MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU) +void grt_source_coef_SH(const MODEL1D_STATE *mstat, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU) { - size_t isrc = mod1d->isrc; - cplx_t xb = mod1d->xb[isrc]; - cplx_t cbcb = mod1d->cbcb[isrc]; - real_t k = mod1d->k; + size_t isrc = mstat->mod1d->isrc; + cplx_t xb = mstat->xb[isrc]; + cplx_t cbcb = mstat->cbcb[isrc]; + real_t k = mstat->k; _source_SH(xb, cbcb, k, src_coefD, src_coefU); } diff --git a/pygrt/C_extension/src/integral/dwm.c b/pygrt/C_extension/src/integral/dwm.c index 7d2511bb..e74aa86b 100644 --- a/pygrt/C_extension/src/integral/dwm.c +++ b/pygrt/C_extension/src/integral/dwm.c @@ -25,7 +25,7 @@ real_t grt_discrete_integ( - MODEL1D *mod1d, real_t dk, real_t kmax, real_t keps, + MODEL1D_STATE *mstat, real_t dk, real_t kmax, real_t keps, size_t nr, real_t *rs, K_INTEG *K, FILE *fstats, GRT_KernelFunc kerfunc) { if(kmax == 0.0) return 0.0; @@ -47,8 +47,8 @@ real_t grt_discrete_integ( k += dk; // 计算核函数 F(k, w) - kerfunc(mod1d, k, K->QWV, K->calc_upar, K->QWVz); - if(mod1d->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; + kerfunc(mstat, k, K->QWV, K->calc_upar, K->QWVz); + if(mstat->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; if(K->applyDCM){ GRT_LOOP_ChnlGrid(im, c){ diff --git a/pygrt/C_extension/src/integral/fim.c b/pygrt/C_extension/src/integral/fim.c index 8b5039a6..c3eb4bb8 100755 --- a/pygrt/C_extension/src/integral/fim.c +++ b/pygrt/C_extension/src/integral/fim.c @@ -24,7 +24,7 @@ real_t grt_linear_filon_integ( - MODEL1D *mod1d, real_t k0, real_t dk0, real_t dk, real_t kmax, real_t keps, + MODEL1D_STATE *mstat, real_t k0, real_t dk0, real_t dk, real_t kmax, real_t keps, size_t nr, real_t *rs, K_INTEG *K, FILE *fstats, GRT_KernelFunc kerfunc) { if(k0 + dk0 >= kmax) return k0; @@ -49,8 +49,8 @@ real_t grt_linear_filon_integ( k += dk; // 计算核函数 F(k, w) - kerfunc(mod1d, k, K->QWV, K->calc_upar, K->QWVz); - if(mod1d->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; + kerfunc(mstat, k, K->QWV, K->calc_upar, K->QWVz); + if(mstat->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; if(K->applyDCM){ GRT_LOOP_ChnlGrid(im, c){ @@ -156,8 +156,8 @@ real_t grt_linear_filon_integ( } // 计算核函数 F(k, w) - kerfunc(mod1d, k0N, K->QWV, K->calc_upar, K->QWVz); - if(mod1d->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; + kerfunc(mstat, k0N, K->QWV, K->calc_upar, K->QWVz); + if(mstat->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; if(K->applyDCM){ GRT_LOOP_ChnlGrid(im, c){ diff --git a/pygrt/C_extension/src/integral/integ_method.c b/pygrt/C_extension/src/integral/integ_method.c index 38b5f7ca..35bf8de0 100644 --- a/pygrt/C_extension/src/integral/integ_method.c +++ b/pygrt/C_extension/src/integral/integ_method.c @@ -78,7 +78,7 @@ void grt_KMET_destroy_fstats(const size_t nr, K_INTEG_METHOD *Kmet) K_INTEG * grt_wavenumber_integral( - MODEL1D *mod1d, size_t nr, real_t *rs, K_INTEG_METHOD *Kmet, bool calc_upar, GRT_KernelFunc kerfunc) + MODEL1D_STATE *mstat, size_t nr, real_t *rs, K_INTEG_METHOD *Kmet, bool calc_upar, GRT_KernelFunc kerfunc) { real_t k = 0.0; @@ -97,37 +97,37 @@ K_INTEG * grt_wavenumber_integral( // 准备 DCM，计算波数上限处的核函数 if(Kint->applyDCM){ - kerfunc(mod1d, Kint->kmax, Kint->QWV_kmax, Kint->calc_upar, Kint->QWVz_kmax); + kerfunc(mstat, Kint->kmax, Kint->QWV_kmax, Kint->calc_upar, Kint->QWVz_kmax); } // DWM k = grt_discrete_integ( - mod1d, Kmet->dk, kcut, Kmet->keps, nr, rs, + mstat, Kmet->dk, kcut, Kmet->keps, nr, rs, Kint, Kmet->fstats, kerfunc); - if(mod1d->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; + if(mstat->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; // 基于线性插值的Filon积分，固定采样间隔 if(Kmet->applyFIM){ k = grt_linear_filon_integ( - mod1d, k, Kmet->dk, Kmet->filondk, Kmet->kmax, Kmet->keps, nr, rs, + mstat, k, Kmet->dk, Kmet->filondk, Kmet->kmax, Kmet->keps, nr, rs, Kint, Kmet->fstats, kerfunc); - if(mod1d->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; + if(mstat->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; } // 基于自适应采样的Filon积分 else if(Kmet->applySAFIM){ - real_t kref = (creal(mod1d->omega) > 0.0) ? creal(mod1d->omega) / Kmet->vmin * Kmet->ampk : 0.0; // 静态解中频率位置给了负值 + real_t kref = (creal(mstat->omega) > 0.0) ? creal(mstat->omega) / Kmet->vmin * Kmet->ampk : 0.0; // 静态解中频率位置给了负值 k = grt_sa_filon_integ( - mod1d, k, Kmet->dk, Kmet->sa_tol, Kmet->kmax, kref, nr, rs, + mstat, k, Kmet->dk, Kmet->sa_tol, Kmet->kmax, kref, nr, rs, Kint, Kmet->fstats, kerfunc); - if(mod1d->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; + if(mstat->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; } // 显式收敛 if(Kmet->applyPTAM){ grt_PTA_method( - mod1d, k, Kmet->dk, nr, rs, + mstat, k, Kmet->dk, nr, rs, Kint, Kmet->ptam_fstatsnr, kerfunc); - if(mod1d->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; + if(mstat->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; } else if(Kmet->applyDCM){ grt_dcm_correction(nr, rs, Kint, !isFilon); diff --git a/pygrt/C_extension/src/integral/kernel_template.c_ b/pygrt/C_extension/src/integral/kernel_template.c_ index 56146ab8..394d6f21 100644 --- a/pygrt/C_extension/src/integral/kernel_template.c_ +++ b/pygrt/C_extension/src/integral/kernel_template.c_ @@ -57,27 +57,27 @@ #define __grt_wave2qwv_z_REV_SH CONCAT_PREFIX(wave2qwv_z_REV_SH) -void __GRT_MATRIX_FUNC__(GRT_MODEL1D *mod1d, const real_t k) +void __GRT_MATRIX_FUNC__(MODEL1D_STATE *mstat, const real_t k) { - mod1d->k = k; + mstat->k = k; // 为动态解计算xa,xb,caca,cbcb #if defined(__DYNAMIC_KERNEL__) - grt_mod1d_xa_xb(mod1d, k); + grt_update_mod1d_state_k(mstat, k); #endif - size_t isrc = mod1d->isrc; // 震源所在虚拟层位, isrc>=1 - size_t ircv = mod1d->ircv; // 接收点所在虚拟层位, ircv>=1, ircv != isrc + size_t isrc = mstat->mod1d->isrc; // 震源所在虚拟层位, isrc>=1 + size_t ircv = mstat->mod1d->ircv; // 接收点所在虚拟层位, ircv>=1, ircv != isrc size_t imin, imax; // 相对浅层深层层位 imin = GRT_MIN(isrc, ircv); imax = GRT_MAX(isrc, ircv); // 广义层的反射透射系数 - RT_MATRIX *M_BL = &mod1d->M_BL; - RT_MATRIX *M_RS = &mod1d->M_RS; - RT_MATRIX *M_FA = &mod1d->M_FA; + RT_MATRIX *M_BL = &mstat->M_BL; + RT_MATRIX *M_RS = &mstat->M_RS; + RT_MATRIX *M_FA = &mstat->M_FA; // 顶底界面 - RT_MATRIX *M_top = &mod1d->M_top; - RT_MATRIX *M_bot = &mod1d->M_bot; + RT_MATRIX *M_top = &mstat->M_top; + RT_MATRIX *M_bot = &mstat->M_bot; grt_reset_RT_matrix(M_BL); grt_reset_RT_matrix(M_RS); @@ -85,8 +85,10 @@ void __GRT_MATRIX_FUNC__(GRT_MODEL1D *mod1d, const real_t k) grt_reset_RT_matrix(M_top); grt_reset_RT_matrix(M_bot); + size_t nlay = mstat->mod1d->n; + // 从顶到底进行矩阵递推, 公式(5.5.3) - size_t maxidx = (imax+1 == mod1d->n)? mod1d->n : mod1d->n - 1; + size_t maxidx = (imax+1 == nlay)? nlay : nlay - 1; for(size_t iy = 1; iy < maxidx; ++iy) // 因为n>=3, 故一定会进入该循环 { // 定义物理层内的反射透射系数矩阵 @@ -99,15 +101,15 @@ void __GRT_MATRIX_FUNC__(GRT_MODEL1D *mod1d, const real_t k) if(iy != isrc && iy != ircv){ #endif // 对第iy层的系数矩阵赋值 - __grt_RT_matrix_PSV(mod1d, iy, M); - __grt_RT_matrix_SH(mod1d, iy, M); + __grt_RT_matrix_PSV(mstat, iy, M); + __grt_RT_matrix_SH(mstat, iy, M); if(M->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; #ifdef __DYNAMIC_KERNEL__ } #endif // 加入时间延迟因子(第iy-1界面与第iy界面之间) - __grt_delay_RT_matrix(mod1d, iy, M); + __grt_delay_RT_matrix(mstat, iy, M); // FA if(iy <= imin){ @@ -116,7 +118,7 @@ void __GRT_MATRIX_FUNC__(GRT_MODEL1D *mod1d, const real_t k) } #ifdef __DYNAMIC_KERNEL__ else if(iy == imin){ - grt_delay_GRT_matrix(mod1d, iy, M_FA); + grt_delay_GRT_matrix(mstat, iy, M_FA); } #endif else { // 递推FA @@ -131,7 +133,7 @@ void __GRT_MATRIX_FUNC__(GRT_MODEL1D *mod1d, const real_t k) } #ifdef __DYNAMIC_KERNEL__ else if(iy == imax){ - grt_delay_GRT_matrix(mod1d, iy, M_RS); + grt_delay_GRT_matrix(mstat, iy, M_RS); } #endif else { // 递推RS @@ -155,16 +157,16 @@ void __GRT_MATRIX_FUNC__(GRT_MODEL1D *mod1d, const real_t k) //=================================================================================== // 递推RU_FA - __grt_topbound_RU_PSV(mod1d); - __grt_topbound_RU_SH(mod1d); + __grt_topbound_RU_PSV(mstat); + __grt_topbound_RU_SH(mstat); if(M_top->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; grt_recursion_RU(M_top, M_FA, M_FA); if(M_FA->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; // 递推RD_BL, 在模型定义阶段已检查震源/台站的虚拟层成为底层时，只能是半无限空间 - if(imax+1 < mod1d->n){ // 此 if 仅减少不必要的计算 - __grt_botbound_RD_PSV(mod1d); - __grt_botbound_RD_SH(mod1d); + if(imax+1 < nlay){ // 此 if 仅减少不必要的计算 + __grt_botbound_RD_PSV(mstat); + __grt_botbound_RD_SH(mstat); if(M_bot->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; grt_recursion_RD(M_BL, M_bot, M_BL); if(M_BL->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; @@ -176,26 +178,25 @@ void __GRT_MATRIX_FUNC__(GRT_MODEL1D *mod1d, const real_t k) } -void __BUILD_QWV__( - GRT_MODEL1D *mod1d, cplxChnlGrid QWV, - bool calc_uiz, cplxChnlGrid QWV_uiz) +void __BUILD_QWV__(MODEL1D_STATE *mstat, cplxChnlGrid QWV, bool calc_uiz, cplxChnlGrid QWV_uiz) { // 初始化qwv为0 memset(QWV, 0, sizeof(cplxChnlGrid)); if(calc_uiz) memset(QWV_uiz, 0, sizeof(cplxChnlGrid)); + MODEL1D *mod1d = mstat->mod1d; bool ircvup = mod1d->ircvup; // 广义层的反射透射系数 - RT_MATRIX *M_BL = &mod1d->M_BL; - RT_MATRIX *M_AL = &mod1d->M_AL; - RT_MATRIX *M_RS = &mod1d->M_RS; - RT_MATRIX *M_FA = &mod1d->M_FA; - RT_MATRIX *M_FB = &mod1d->M_FB; + RT_MATRIX *M_BL = &mstat->M_BL; + RT_MATRIX *M_AL = &mstat->M_AL; + RT_MATRIX *M_RS = &mstat->M_RS; + RT_MATRIX *M_FA = &mstat->M_FA; + RT_MATRIX *M_FB = &mstat->M_FB; // 计算震源系数 P_m, SV_m, SH_m cplxChnlGrid src_coefD; cplxChnlGrid src_coefU; - __grt_source_coef(mod1d, src_coefD, src_coefU); + __grt_source_coef(mstat, src_coefD, src_coefU); // 临时中转矩阵 (temperary) cplx_t tmpR2[2][2], tmp2x2[2][2], tmpRL, tmp2x2_uiz[2][2], tmpRL2; @@ -204,8 +205,8 @@ void __BUILD_QWV__( if(ircvup){ // A接收 B震源 // 计算R_EV - __grt_wave2qwv_REV_PSV(mod1d); - __grt_wave2qwv_REV_SH(mod1d); + __grt_wave2qwv_REV_PSV(mstat); + __grt_wave2qwv_REV_SH(mstat); // 递推RU_FS grt_recursion_RU(M_FA, M_RS, M_FB); // 已从ZR变为FR，加入了自由表面的效应 @@ -219,18 +220,18 @@ void __BUILD_QWV__( if(calc_uiz) grt_cmat2x2_assign(tmp2x2, tmp2x2_uiz); // 为后续计算空间导数备份 - grt_cmat2x2_mul(mod1d->R_EV, tmp2x2, tmp2x2); + grt_cmat2x2_mul(mstat->R_EV, tmp2x2, tmp2x2); tmpRL = M_FB->invTL / (1.0 - M_BL->RDL * M_FB->RUL); - tmpRL2 = mod1d->R_EVL * tmpRL; + tmpRL2 = mstat->R_EVL * tmpRL; grt_construct_qwv(ircvup, tmp2x2, tmpRL2, M_BL->RD, M_BL->RDL, src_coefD, src_coefU, QWV); if(calc_uiz){ - __grt_wave2qwv_z_REV_PSV(mod1d); - __grt_wave2qwv_z_REV_SH(mod1d); - grt_cmat2x2_mul(mod1d->uiz_R_EV, tmp2x2_uiz, tmp2x2_uiz); - tmpRL2 = mod1d->uiz_R_EVL * tmpRL; + __grt_wave2qwv_z_REV_PSV(mstat); + __grt_wave2qwv_z_REV_SH(mstat); + grt_cmat2x2_mul(mstat->uiz_R_EV, tmp2x2_uiz, tmp2x2_uiz); + tmpRL2 = mstat->uiz_R_EVL * tmpRL; grt_construct_qwv(ircvup, tmp2x2_uiz, tmpRL2, M_BL->RD, M_BL->RDL, src_coefD, src_coefU, QWV_uiz); } @@ -238,8 +239,8 @@ void __BUILD_QWV__( else { // A震源 B接收 // 计算R_EV - __grt_wave2qwv_REV_PSV(mod1d); - __grt_wave2qwv_REV_SH(mod1d); + __grt_wave2qwv_REV_PSV(mstat); + __grt_wave2qwv_REV_SH(mstat); // 递推RD_SL grt_recursion_RD(M_RS, M_BL, M_AL); @@ -253,18 +254,18 @@ void __BUILD_QWV__( if(calc_uiz) grt_cmat2x2_assign(tmp2x2, tmp2x2_uiz); // 为后续计算空间导数备份 - grt_cmat2x2_mul(mod1d->R_EV, tmp2x2, tmp2x2); + grt_cmat2x2_mul(mstat->R_EV, tmp2x2, tmp2x2); tmpRL = M_AL->invTL / (1.0 - M_FA->RUL * M_AL->RDL); - tmpRL2 = mod1d->R_EVL * tmpRL; + tmpRL2 = mstat->R_EVL * tmpRL; grt_construct_qwv(ircvup, tmp2x2, tmpRL2, M_FA->RU, M_FA->RUL, src_coefD, src_coefU, QWV); if(calc_uiz){ - __grt_wave2qwv_z_REV_PSV(mod1d); - __grt_wave2qwv_z_REV_SH(mod1d); - grt_cmat2x2_mul(mod1d->uiz_R_EV, tmp2x2_uiz, tmp2x2_uiz); - tmpRL2 = mod1d->uiz_R_EVL * tmpRL; + __grt_wave2qwv_z_REV_PSV(mstat); + __grt_wave2qwv_z_REV_SH(mstat); + grt_cmat2x2_mul(mstat->uiz_R_EV, tmp2x2_uiz, tmp2x2_uiz); + tmpRL2 = mstat->uiz_R_EVL * tmpRL; grt_construct_qwv(ircvup, tmp2x2_uiz, tmpRL2, M_FA->RU, M_FA->RUL, src_coefD, src_coefU, QWV_uiz); } @@ -287,12 +288,10 @@ void __BUILD_QWV__( } -void __KERNEL_FUNC__( - GRT_MODEL1D *mod1d, const real_t k, cplxChnlGrid QWV, - bool calc_uiz, cplxChnlGrid QWV_uiz) +void __KERNEL_FUNC__(MODEL1D_STATE *mstat, const real_t k, cplxChnlGrid QWV, bool calc_uiz, cplxChnlGrid QWV_uiz) { - __GRT_MATRIX_FUNC__(mod1d, k); - __BUILD_QWV__(mod1d, QWV, calc_uiz, QWV_uiz); + __GRT_MATRIX_FUNC__(mstat, k); + __BUILD_QWV__(mstat, QWV, calc_uiz, QWV_uiz); } diff --git a/pygrt/C_extension/src/integral/ptam.c b/pygrt/C_extension/src/integral/ptam.c index 985a8dc8..7e10341d 100755 --- a/pygrt/C_extension/src/integral/ptam.c +++ b/pygrt/C_extension/src/integral/ptam.c @@ -212,7 +212,7 @@ static void _cplx_shrink(size_t n1, size_t ir, int im, int v, cplxIntegGrid (*F void grt_PTA_method( - MODEL1D *mod1d, real_t k0, real_t predk, + MODEL1D_STATE *mstat, real_t k0, real_t predk, size_t nr, real_t *rs, K_INTEG *K, FILE *ptam_fstatsnr[nr][2], GRT_KernelFunc kerfunc) { // 需要兼容对正常收敛而不具有规律波峰波谷的序列 @@ -275,8 +275,8 @@ void grt_PTA_method( k += dk; // 计算核函数 F(k, w) - kerfunc(mod1d, k, K->QWV, K->calc_upar, K->QWVz); - if(mod1d->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; + kerfunc(mstat, k, K->QWV, K->calc_upar, K->QWVz); + if(mstat->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; // 记录核函数 if(ptam_fstatsnr != NULL) grt_write_stats(ptam_fstatsnr[ir][0], k, (K->calc_upar)? K->QWVz : K->QWV); diff --git a/pygrt/C_extension/src/integral/safim.c b/pygrt/C_extension/src/integral/safim.c index 27c02e93..4d8ffec2 100644 --- a/pygrt/C_extension/src/integral/safim.c +++ b/pygrt/C_extension/src/integral/safim.c @@ -276,7 +276,7 @@ static void interv_integ(const KInterval *ptKitv, size_t nr, real_t *rs, K_INTEG real_t grt_sa_filon_integ( - MODEL1D *mod1d, real_t k0, real_t dk0, real_t tol, real_t kmax, real_t kref, + MODEL1D_STATE *mstat, real_t k0, real_t dk0, real_t tol, real_t kmax, real_t kref, size_t nr, real_t *rs, K_INTEG *K, FILE *fstats, GRT_KernelFunc kerfunc) { real_t kmin = k0 + dk0; @@ -300,8 +300,8 @@ real_t grt_sa_filon_integ( .Fz3 = {{{0}}} }; for(int i=0; i<3; ++i) { - kerfunc(mod1d, Kitv.k3[i], Kitv.F3[i], K->calc_upar, Kitv.Fz3[i]); - if(mod1d->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; + kerfunc(mstat, Kitv.k3[i], Kitv.F3[i], K->calc_upar, Kitv.Fz3[i]); + if(mstat->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; if(K->applyDCM){ GRT_LOOP_ChnlGrid(im, c){ @@ -355,11 +355,11 @@ real_t grt_sa_filon_integ( memcpy(Kitv_right.Fz3[2], Kitv.Fz3[2], sizeof(cplxChnlGrid)); } - kerfunc(mod1d, Kitv_left.k3[1], Kitv_left.F3[1], K->calc_upar, Kitv_left.Fz3[1]); - if(mod1d->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; + kerfunc(mstat, Kitv_left.k3[1], Kitv_left.F3[1], K->calc_upar, Kitv_left.Fz3[1]); + if(mstat->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; - kerfunc(mod1d, Kitv_right.k3[1], Kitv_right.F3[1], K->calc_upar, Kitv_right.Fz3[1]); - if(mod1d->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; + kerfunc(mstat, Kitv_right.k3[1], Kitv_right.F3[1], K->calc_upar, Kitv_right.Fz3[1]); + if(mstat->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; if(K->applyDCM){ GRT_LOOP_ChnlGrid(im, c){ diff --git a/pygrt/C_extension/src/static/static_grn.c b/pygrt/C_extension/src/static/static_grn.c index 03ac00b2..691527f9 100644 --- a/pygrt/C_extension/src/static/static_grn.c +++ b/pygrt/C_extension/src/static/static_grn.c @@ -95,9 +95,12 @@ void grt_integ_static_grn( // Wavenumber Integration // 波数积分上限 Kmet->kmax = Kmet->k0; - K_INTEG *Kint = grt_wavenumber_integral(mod1d, uniq_nr, uniq_rs, Kmet, calc_upar, grt_static_kernel); + // 模型状态 + MODEL1D_STATE *mstat = grt_init_mod1d_state(mod1d); + grt_update_mod1d_state_omega(mstat, -1.0); + K_INTEG *Kint = grt_wavenumber_integral(mstat, uniq_nr, uniq_rs, Kmet, calc_upar, grt_static_kernel); - cplx_t src_mu = mod1d->mu[mod1d->isrc]; + cplx_t src_mu = mstat->mu[mod1d->isrc]; cplx_t fac = Kmet->dk * 1.0/(4.0*PI * src_mu); // 将积分结果记录到浮点数数组中 @@ -125,6 +128,7 @@ void grt_integ_static_grn( // Free allocated memory for temporary variables grt_free_K_INTEG(Kint); + grt_free_mod1d_state(mstat); if(needfstats) grt_KMET_destroy_fstats(uniq_nr, Kmet); diff --git a/pygrt/C_extension/src/static/static_layer.c b/pygrt/C_extension/src/static/static_layer.c index 17d5907c..79f3a304 100644 --- a/pygrt/C_extension/src/static/static_layer.c +++ b/pygrt/C_extension/src/static/static_layer.c @@ -21,9 +21,9 @@ #include "grt/common/matrix.h" /* 定义用于提取相邻两层物性参数的宏 */ -#define MODEL_2LAYS_ATTRIB(T, N) \ - T N##1 = mod1d->N[iy-1];\ - T N##2 = mod1d->N[iy];\ +#define MODEL_2LAYS_ATTRIB(MM, T, N) \ + T N##1 = MM->N[iy-1];\ + T N##2 = MM->N[iy];\ static int __freeBound_R_PSV(real_t d, real_t k, cplx_t delta1, cplx_t M[2][2]) @@ -67,19 +67,20 @@ static int __rigidBound_R_SH(cplx_t *ML) return GRT_INVERSE_SUCCESS; } -void grt_static_topbound_RU_PSV(MODEL1D *mod1d) +void grt_static_topbound_RU_PSV(MODEL1D_STATE *mstat) { - cplx_t delta = mod1d->delta[0]; - real_t k = mod1d->k; + MODEL1D *mod1d = mstat->mod1d; + cplx_t delta = mstat->delta[0]; + real_t k = mstat->k; if(mod1d->topbound == GRT_BOUND_FREE){ - mod1d->M_top.stats = __freeBound_R_PSV(0.0, k, delta, mod1d->M_top.RU); + mstat->M_top.stats = __freeBound_R_PSV(0.0, k, delta, mstat->M_top.RU); } else if(mod1d->topbound == GRT_BOUND_RIGID){ - mod1d->M_top.stats = __rigidBound_R_PSV(0.0, k, delta, mod1d->M_top.RU); + mstat->M_top.stats = __rigidBound_R_PSV(0.0, k, delta, mstat->M_top.RU); } else if(mod1d->topbound == GRT_BOUND_HALFSPACE){ - memset(mod1d->M_top.RU, 0, sizeof(cplx_t)*4); - mod1d->M_top.stats = GRT_INVERSE_SUCCESS; + memset(mstat->M_top.RU, 0, sizeof(cplx_t)*4); + mstat->M_top.stats = GRT_INVERSE_SUCCESS; } else{ GRTRaiseError("Wrong execution."); @@ -87,17 +88,18 @@ void grt_static_topbound_RU_PSV(MODEL1D *mod1d) // RU 不需要时延 } -void grt_static_topbound_RU_SH(MODEL1D *mod1d) +void grt_static_topbound_RU_SH(MODEL1D_STATE *mstat) { + MODEL1D *mod1d = mstat->mod1d; if(mod1d->topbound == GRT_BOUND_FREE){ - mod1d->M_top.stats = __freeBound_R_SH(&mod1d->M_top.RUL); + mstat->M_top.stats = __freeBound_R_SH(&mstat->M_top.RUL); } else if(mod1d->topbound == GRT_BOUND_RIGID){ - mod1d->M_top.stats = __rigidBound_R_SH(&mod1d->M_top.RUL); + mstat->M_top.stats = __rigidBound_R_SH(&mstat->M_top.RUL); } else if(mod1d->topbound == GRT_BOUND_HALFSPACE){ - mod1d->M_top.RUL = 0.0; - mod1d->M_top.stats = GRT_INVERSE_SUCCESS; + mstat->M_top.RUL = 0.0; + mstat->M_top.stats = GRT_INVERSE_SUCCESS; } else{ GRTRaiseError("Wrong execution."); @@ -105,20 +107,21 @@ void grt_static_topbound_RU_SH(MODEL1D *mod1d) // RU 不需要时延 } -void grt_static_botbound_RD_PSV(MODEL1D *mod1d) +void grt_static_botbound_RD_PSV(MODEL1D_STATE *mstat) { + MODEL1D *mod1d = mstat->mod1d; size_t nlay = mod1d->n; - cplx_t delta = mod1d->delta[nlay-2]; + cplx_t delta = mstat->delta[nlay-2]; real_t thk = mod1d->Thk[nlay-2]; - real_t k = mod1d->k; + real_t k = mstat->k; if(mod1d->botbound == GRT_BOUND_FREE){ - mod1d->M_bot.stats = __freeBound_R_PSV(thk, k, delta, mod1d->M_bot.RD); + mstat->M_bot.stats = __freeBound_R_PSV(thk, k, delta, mstat->M_bot.RD); } else if(mod1d->botbound == GRT_BOUND_RIGID){ - mod1d->M_bot.stats = __rigidBound_R_PSV(thk, k, delta, mod1d->M_bot.RD); + mstat->M_bot.stats = __rigidBound_R_PSV(thk, k, delta, mstat->M_bot.RD); } else if(mod1d->botbound == GRT_BOUND_HALFSPACE){ - grt_static_RT_matrix_PSV(mod1d, nlay-1, &mod1d->M_bot); + grt_static_RT_matrix_PSV(mstat, nlay-1, &mstat->M_bot); } else{ GRTRaiseError("Wrong execution."); @@ -129,23 +132,24 @@ void grt_static_botbound_RD_PSV(MODEL1D *mod1d) ex = exp(- k*thk); ex2 = ex * ex; - mod1d->M_bot.RD[0][0] *= ex2; mod1d->M_bot.RD[0][1] *= ex2; - mod1d->M_bot.RD[1][0] *= ex2; mod1d->M_bot.RD[1][1] *= ex2; + mstat->M_bot.RD[0][0] *= ex2; mstat->M_bot.RD[0][1] *= ex2; + mstat->M_bot.RD[1][0] *= ex2; mstat->M_bot.RD[1][1] *= ex2; } -void grt_static_botbound_RD_SH(MODEL1D *mod1d) +void grt_static_botbound_RD_SH(MODEL1D_STATE *mstat) { + MODEL1D *mod1d = mstat->mod1d; size_t nlay = mod1d->n; real_t thk = mod1d->Thk[nlay-2]; - real_t k = mod1d->k; + real_t k = mstat->k; if(mod1d->botbound == GRT_BOUND_FREE){ - mod1d->M_bot.stats = __freeBound_R_SH(&mod1d->M_bot.RDL); + mstat->M_bot.stats = __freeBound_R_SH(&mstat->M_bot.RDL); } else if(mod1d->botbound == GRT_BOUND_RIGID){ - mod1d->M_bot.stats = __rigidBound_R_SH(&mod1d->M_bot.RDL); + mstat->M_bot.stats = __rigidBound_R_SH(&mstat->M_bot.RDL); } else if(mod1d->botbound == GRT_BOUND_HALFSPACE){ - grt_static_RT_matrix_SH(mod1d, nlay-1, &mod1d->M_bot); + grt_static_RT_matrix_SH(mstat, nlay-1, &mstat->M_bot); } else{ GRTRaiseError("Wrong execution."); @@ -156,70 +160,71 @@ void grt_static_botbound_RD_SH(MODEL1D *mod1d) ex = exp(- k*thk); ex2 = ex * ex; - mod1d->M_bot.RDL *= ex2; + mstat->M_bot.RDL *= ex2; } -void grt_static_wave2qwv_REV_PSV(MODEL1D *mod1d) +void grt_static_wave2qwv_REV_PSV(MODEL1D_STATE *mstat) { cplx_t D11[2][2] = {{1.0, -1.0}, {1.0, 1.0}}; cplx_t D12[2][2] = {{1.0, -1.0}, {-1.0, -1.0}}; // 公式(6.3.35,37) - if(mod1d->ircvup){// 震源更深 - grt_cmat2x2_mul(D12, mod1d->M_FA.RU, mod1d->R_EV); - grt_cmat2x2_add(D11, mod1d->R_EV, mod1d->R_EV); + if(mstat->mod1d->ircvup){// 震源更深 + grt_cmat2x2_mul(D12, mstat->M_FA.RU, mstat->R_EV); + grt_cmat2x2_add(D11, mstat->R_EV, mstat->R_EV); } else { // 接收点更深 - grt_cmat2x2_mul(D11, mod1d->M_BL.RD, mod1d->R_EV); - grt_cmat2x2_add(D12, mod1d->R_EV, mod1d->R_EV); + grt_cmat2x2_mul(D11, mstat->M_BL.RD, mstat->R_EV); + grt_cmat2x2_add(D12, mstat->R_EV, mstat->R_EV); } } -void grt_static_wave2qwv_REV_SH(MODEL1D *mod1d) +void grt_static_wave2qwv_REV_SH(MODEL1D_STATE *mstat) { - if(mod1d->ircvup){// 震源更深 - mod1d->R_EVL = 1.0 + mod1d->M_FA.RUL; + if(mstat->mod1d->ircvup){// 震源更深 + mstat->R_EVL = 1.0 + mstat->M_FA.RUL; } else { - mod1d->R_EVL = 1.0 + mod1d->M_BL.RDL; + mstat->R_EVL = 1.0 + mstat->M_BL.RDL; } } -void grt_static_wave2qwv_z_REV_PSV(MODEL1D *mod1d) +void grt_static_wave2qwv_z_REV_PSV(MODEL1D_STATE *mstat) { - real_t k = mod1d->k; + MODEL1D *mod1d = mstat->mod1d; + real_t k = mstat->k; size_t ircv = mod1d->ircv; - cplx_t delta1 = mod1d->delta[ircv]; + cplx_t delta1 = mstat->delta[ircv]; // 新推导公式 cplx_t kd2 = 2.0*k*delta1; cplx_t D11[2][2] = {{k, -k-kd2}, {k, k-kd2}}; cplx_t D12[2][2] = {{-k, k+kd2}, {k, k-kd2}}; if(mod1d->ircvup){// 震源更深 - grt_cmat2x2_mul(D12, mod1d->M_FA.RU, mod1d->uiz_R_EV); - grt_cmat2x2_add(D11, mod1d->uiz_R_EV, mod1d->uiz_R_EV); + grt_cmat2x2_mul(D12, mstat->M_FA.RU, mstat->uiz_R_EV); + grt_cmat2x2_add(D11, mstat->uiz_R_EV, mstat->uiz_R_EV); } else { // 接收点更深 - grt_cmat2x2_mul(D11, mod1d->M_BL.RD, mod1d->uiz_R_EV); - grt_cmat2x2_add(D12, mod1d->uiz_R_EV, mod1d->uiz_R_EV); + grt_cmat2x2_mul(D11, mstat->M_BL.RD, mstat->uiz_R_EV); + grt_cmat2x2_add(D12, mstat->uiz_R_EV, mstat->uiz_R_EV); } } -void grt_static_wave2qwv_z_REV_SH(MODEL1D *mod1d) +void grt_static_wave2qwv_z_REV_SH(MODEL1D_STATE *mstat) { - real_t k = mod1d->k; + real_t k = mstat->k; // 新推导公式 - if(mod1d->ircvup){// 震源更深 - mod1d->uiz_R_EVL = (1.0 - mod1d->M_FA.RUL)*k; + if(mstat->mod1d->ircvup){// 震源更深 + mstat->uiz_R_EVL = (1.0 - mstat->M_FA.RUL)*k; } else { // 接收点更深 - mod1d->uiz_R_EVL = (mod1d->M_BL.RDL - 1.0)*k; + mstat->uiz_R_EVL = (mstat->M_BL.RDL - 1.0)*k; } } -void grt_static_RT_matrix_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_static_RT_matrix_PSV(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M) { - MODEL_2LAYS_ATTRIB(cplx_t, mu); - MODEL_2LAYS_ATTRIB(cplx_t, delta); - real_t thk = mod1d->Thk[iy-1]; - real_t k = mod1d->k; + MODEL_2LAYS_ATTRIB(mstat, cplx_t, mu); + MODEL_2LAYS_ATTRIB(mstat, cplx_t, delta); + real_t thk = mstat->mod1d->Thk[iy-1]; + real_t k = mstat->k; // 公式(6.3.18) cplx_t dmu = mu1 - mu2; @@ -256,9 +261,9 @@ void grt_static_RT_matrix_PSV(const MODEL1D *mod1d, const size_t iy, RT_MATRIX * } -void grt_static_RT_matrix_SH(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_static_RT_matrix_SH(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M) { - MODEL_2LAYS_ATTRIB(cplx_t, mu); + MODEL_2LAYS_ATTRIB(mstat, cplx_t, mu); // 公式(6.3.18) cplx_t dmu = mu1 - mu2; @@ -274,10 +279,10 @@ void grt_static_RT_matrix_SH(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M } -void grt_static_delay_RT_matrix(const MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) +void grt_static_delay_RT_matrix(const MODEL1D_STATE *mstat, const size_t iy, RT_MATRIX *M) { - real_t thk = mod1d->Thk[iy-1]; - real_t k = mod1d->k; + real_t thk = mstat->mod1d->Thk[iy-1]; + real_t k = mstat->k; cplx_t ex, ex2; ex = exp(- k*thk); diff --git a/pygrt/C_extension/src/static/static_source.c b/pygrt/C_extension/src/static/static_source.c index 3fb0044e..ed956197 100644 --- a/pygrt/C_extension/src/static/static_source.c +++ b/pygrt/C_extension/src/static/static_source.c @@ -61,30 +61,30 @@ inline GCC_ALWAYS_INLINE void _source_SH(const real_t k, cplxChnlGrid coefD, cpl } -void grt_static_source_coef(const MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU) +void grt_static_source_coef(const MODEL1D_STATE *mstat, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU) { // 先全部赋0 memset(src_coefD, 0, sizeof(cplxChnlGrid)); memset(src_coefU, 0, sizeof(cplxChnlGrid)); - grt_static_source_coef_PSV(mod1d, src_coefD, src_coefU); - grt_static_source_coef_SH(mod1d, src_coefD, src_coefU); + grt_static_source_coef_PSV(mstat, src_coefD, src_coefU); + grt_static_source_coef_SH(mstat, src_coefD, src_coefU); } -void grt_static_source_coef_PSV(const MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU) +void grt_static_source_coef_PSV(const MODEL1D_STATE *mstat, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU) { - size_t isrc = mod1d->isrc; - cplx_t delta = mod1d->delta[isrc]; - real_t k = mod1d->k; + size_t isrc = mstat->mod1d->isrc; + cplx_t delta = mstat->delta[isrc]; + real_t k = mstat->k; _source_PSV(delta, k, src_coefD, src_coefU); } -void grt_static_source_coef_SH(const MODEL1D *mod1d, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU) +void grt_static_source_coef_SH(const MODEL1D_STATE *mstat, cplxChnlGrid src_coefD, cplxChnlGrid src_coefU) { - real_t k = mod1d->k; + real_t k = mstat->k; _source_SH(k, src_coefD, src_coefU); } diff --git a/pygrt/c_interfaces.py b/pygrt/c_interfaces.py index a3cf2fa5..9a6c8bc8 100755 --- a/pygrt/c_interfaces.py +++ b/pygrt/c_interfaces.py @@ -30,7 +30,7 @@ C_grt_integ_grn_spec = libgrt.grt_integ_grn_spec """C库中计算格林函数的主函数 integ_grn_spec, 详见C API同名函数""" C_grt_integ_grn_spec.argtypes = [ - POINTER(c_GRT_MODEL1D), c_size_t, c_size_t, PREAL, + POINTER(c_MODEL1D), c_size_t, c_size_t, PREAL, c_size_t, PREAL, REAL, c_bool, POINTER(c_K_INTEG_METHOD), c_bool, c_bool, POINTER((PCPLX*CHANNEL_NUM)*SRC_M_NUM), @@ -47,7 +47,7 @@ """计算静态格林函数""" C_grt_integ_static_grn.restype = None C_grt_integ_static_grn.argtypes = [ - POINTER(c_GRT_MODEL1D), c_size_t, PREAL, POINTER(c_K_INTEG_METHOD), + POINTER(c_MODEL1D), c_size_t, PREAL, POINTER(c_K_INTEG_METHOD), c_bool, POINTER((REAL*CHANNEL_NUM)*SRC_M_NUM), POINTER((REAL*CHANNEL_NUM)*SRC_M_NUM), @@ -96,18 +96,18 @@ def set_num_threads(n): C_grt_read_mod1d_from_file = libgrt.grt_read_mod1d_from_file """读取模型文件并进行预处理""" -C_grt_read_mod1d_from_file.restype = POINTER(c_GRT_MODEL1D) +C_grt_read_mod1d_from_file.restype = POINTER(c_MODEL1D) C_grt_read_mod1d_from_file.argtypes = [c_char_p, c_double, c_double, c_bool] C_grt_set_mod1d_boundary = libgrt.grt_set_mod1d_boundary """设置模型边界条件并检查底界面""" C_grt_set_mod1d_boundary.restype = None -C_grt_set_mod1d_boundary.argtypes = [POINTER(c_GRT_MODEL1D), c_int, c_int] +C_grt_set_mod1d_boundary.argtypes = [POINTER(c_MODEL1D), c_int, c_int] C_grt_free_mod1d = libgrt.grt_free_mod1d """释放C程序中申请的 GRT_MODEL1D 结构体内存""" C_grt_free_mod1d.restype = None -C_grt_free_mod1d.argtypes = [POINTER(c_GRT_MODEL1D)] +C_grt_free_mod1d.argtypes = [POINTER(c_MODEL1D)] # ------------------------------------------------------------------- # C函数定义的时间函数 diff --git a/pygrt/c_structures.py b/pygrt/c_structures.py index bda99656..4d528a85 100755 --- a/pygrt/c_structures.py +++ b/pygrt/c_structures.py @@ -32,7 +32,7 @@ "CPLX", "PCPLX", - "c_GRT_MODEL1D", + "c_MODEL1D", "c_K_INTEG_METHOD" ] @@ -63,29 +63,10 @@ class CPLX(Structure): PREAL = POINTER(REAL) PCPLX = POINTER(CPLX) -class c_RT_MATRIX(Structure): - """ - 和C结构体 RT_MATRIX 匹配 - """ - - _fields_ = [ - ('RD', CPLX * 4), - ('RU', CPLX * 4), - ('TD', CPLX * 4), - ('TU', CPLX * 4), - ('RDL', CPLX), - ('RUL', CPLX), - ('TDL', CPLX), - ('TUL', CPLX), - ('invT', CPLX * 4), - ('invTL', CPLX), - ('stats', c_int) - ] - -class c_GRT_MODEL1D(Structure): +class c_MODEL1D(Structure): """ - 和C结构体 GRT_MODEL1D 作匹配 + 和C结构体 MODEL1D 作匹配 """ _fields_ = [ @@ -97,11 +78,7 @@ class c_GRT_MODEL1D(Structure): ('ircvup', c_bool), ('io_depth', c_bool), ('srcrcv_isInserted', c_bool), - ('omgref', CPLX), - ('omega', CPLX), - ('k', REAL), - ('c_phase', CPLX), ('Thk', PREAL), ('Dep', PREAL), @@ -112,35 +89,10 @@ class c_GRT_MODEL1D(Structure): ('Qb', PREAL), ('Qainv', PREAL), ('Qbinv', PREAL), - ('isLiquid', c_bool), - - ('mu', PCPLX), - ('lambda', PCPLX), - ('delta', PCPLX), - ('atna', PCPLX), - ('atnb', PCPLX), - ('xa', PCPLX), - ('xb', PCPLX), - ('caca', PCPLX), - ('cbcb', PCPLX), - - ('stats', c_int), - - ('M_AL', c_RT_MATRIX), - ('M_BL', c_RT_MATRIX), - ('M_RS', c_RT_MATRIX), - ('M_FA', c_RT_MATRIX), - ('M_FB', c_RT_MATRIX), + ('isLiquid', POINTER(c_bool)), ('topbound', c_int), - ('M_top', c_RT_MATRIX), ('botbound', c_int), - ('M_bot', c_RT_MATRIX), - - ('R_EV', CPLX * 4), - ('R_EVL', CPLX), - ('uiz_R_EV', CPLX * 4), - ('uiz_R_EVL', CPLX), ] diff --git a/pygrt/pymod.py b/pygrt/pymod.py index e9d56555..81809d78 100755 --- a/pygrt/pymod.py +++ b/pygrt/pymod.py @@ -49,7 +49,7 @@ def __init__(self, modarr0:np.ndarray, depsrc:float, deprcv:float, allowLiquid:b ''' self.depsrc:float = depsrc self.deprcv:float = deprcv - self.c_mod1d:c_GRT_MODEL1D + self.c_mod1d:c_MODEL1D self.topbound:str = topbound self.botbound:str = botbound self.hasLiquid = False