diff --git a/docs/create_api_rst.sh b/docs/create_api_rst.sh index 98e92de..ba9961d 100755 --- a/docs/create_api_rst.sh +++ b/docs/create_api_rst.sh @@ -26,49 +26,76 @@ C API EOF -for item in $(ls $SDIR); do - echo $item - if [ -d $SDIR/$item ]; then - mkdir -p $PDIR/$item - content=( - "$item" - "=========================================" - "" - ".. toctree::" - " :maxdepth: 2" - "" - ) - printf "%s\n" "${content[@]}" > $PDIR/$item.rst - - for hnm in $(ls $SDIR/$item); do +# 递归处理目录和文件 +process_directory() { + local source_dir=$1 + local target_dir=$2 + local relative_path=$3 + + echo $source_dir $target_dir $relative_path + + for item in $(ls $source_dir); do + local source_item="$source_dir/$item" + local target_item="$target_dir/$item" + local relative_item="$relative_path/$item" + + if [ -d $source_item ]; then + # 创建对应的目标目录 + mkdir -p "$target_item" + + # 为目录创建.rst文件 + local rst_file="${target_dir}/${item}.rst" + content=( + "$item" + "=========================================" + "" + ".. toctree::" + " :maxdepth: 2" + "" + ) + printf "%s\n" "${content[@]}" > "$rst_file" + + # 将目录添加到父目录的toctree中 + if [ "$relative_path" != "." ]; then + local parent_rst="${target_dir}.rst" + echo " $(basename $source_dir)/${item}" >> "$parent_rst" + else + # 如果是顶级目录下的目录,添加到主toctree + echo " C_extension/${item}" >> "$c_api_rst" + fi + + # 递归处理子目录 + process_directory "$source_item" "$target_item" "$relative_item" + + elif [ -f $source_item ] && [[ $item == *.h ]]; then + # 处理.h文件 + # local base_name="${item%%.*}" + local base_name="${item}" content=( - "$hnm" + "$item" "--------------------------------------" "" - ".. doxygenfile:: $hnm" + ".. doxygenfile:: ${relative_item#./}" " :project: h_PyGRT" "" ) - printf "%s\n" "${content[@]}" > $PDIR/$item/${hnm%%.*}.rst - echo " $item/${hnm%%.*}" >> $PDIR/$item.rst - done - - elif [ -f $SDIR/$item ]; then - content=( - "$item" - "--------------------------------------" - "" - ".. doxygenfile:: $item" - " :project: h_PyGRT" - "" - ) - printf "%s\n" "${content[@]}" > $PDIR/${item%%.*}.rst - fi - - echo " C_extension/${item%%.*}" >> $c_api_rst + printf "%s\n" "${content[@]}" > "${target_dir}/${base_name}.rst" + + # 如果是顶级目录下的文件,添加到主toctree + if [ "$relative_path" == "." ]; then + echo " C_extension/${base_name}" >> "$c_api_rst" + else + # 将文件添加到父目录的toctree中 + local parent_rst="${target_dir}.rst" + echo " $(basename ${target_dir})/${base_name}" >> "$parent_rst" + fi + fi + done +} + +# 从SDIR开始处理 +process_directory "$SDIR" "$PDIR" "." - -done diff --git a/pygrt/C_extension/include/grt.h b/pygrt/C_extension/include/grt.h index 21161b3..8335e50 100644 --- a/pygrt/C_extension/include/grt.h +++ b/pygrt/C_extension/include/grt.h @@ -18,10 +18,10 @@ #include #include -#include "common/const.h" -#include "common/logo.h" +#include "grt/common/const.h" +#include "grt/common/logo.h" +#include "grt/common/checkerror.h" -#include "grt_error.h" #define GRT_MAIN_COMMAND "grt" ///< 主程序名 diff --git a/pygrt/C_extension/include/common/attenuation.h b/pygrt/C_extension/include/grt/common/attenuation.h similarity index 83% rename from pygrt/C_extension/include/common/attenuation.h rename to pygrt/C_extension/include/grt/common/attenuation.h index ed09c73..921b10a 100755 --- a/pygrt/C_extension/include/common/attenuation.h +++ b/pygrt/C_extension/include/grt/common/attenuation.h @@ -9,7 +9,7 @@ #pragma once -#include "common/const.h" +#include "grt/common/const.h" /** * 品质因子Q 对 波速的影响, Futterman causal Q, 参考Aki&Richards, 2009, Chapter 5.5.1 @@ -23,9 +23,9 @@ * * @return atncoef 系数因子,作用在 \f$ k=\omega / c(\omega)\f$的计算 */ -MYCOMPLEX attenuation_law(MYREAL Qinv, MYCOMPLEX omega); +MYCOMPLEX grt_attenuation_law(MYREAL Qinv, MYCOMPLEX omega); /** * attenuation_law函数在python中被调用的版本,长度2的数组分别表示复数的实部和虚部 */ -void py_attenuation_law(MYREAL Qinv, MYREAL omg[2], MYREAL atte[2]); \ No newline at end of file +void grt_py_attenuation_law(MYREAL Qinv, MYREAL omg[2], MYREAL atte[2]); \ No newline at end of file diff --git a/pygrt/C_extension/include/common/bessel.h b/pygrt/C_extension/include/grt/common/bessel.h similarity index 81% rename from pygrt/C_extension/include/common/bessel.h rename to pygrt/C_extension/include/grt/common/bessel.h index 361d61a..5d17d24 100755 --- a/pygrt/C_extension/include/common/bessel.h +++ b/pygrt/C_extension/include/grt/common/bessel.h @@ -8,7 +8,7 @@ #pragma once -#include "common/const.h" +#include "grt/common/const.h" /** * 计算Bessel函数 \f$ J_m(x), m=0,1,2 \f$ @@ -19,7 +19,7 @@ * @param[out] bj2 \f$ J_2(x) \f$ * */ -void bessel012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2); +void grt_bessel012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2); /** @@ -31,4 +31,4 @@ void bessel012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2); * @param[in,out] bj2 传入 \f$ J_2(x) \f$, 返回\f$ J_2^{'}(x) \f$ * */ -void besselp012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2); \ No newline at end of file +void grt_besselp012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2); \ No newline at end of file diff --git a/pygrt/C_extension/include/grt_error.h b/pygrt/C_extension/include/grt/common/checkerror.h similarity index 97% rename from pygrt/C_extension/include/grt_error.h rename to pygrt/C_extension/include/grt/common/checkerror.h index 23848b9..b22c5f8 100644 --- a/pygrt/C_extension/include/grt_error.h +++ b/pygrt/C_extension/include/grt/common/checkerror.h @@ -1,9 +1,9 @@ /** - * @file grt_error.h + * @file checkerror.h * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn) * @date 2025-08 * - * 一些用于报错的宏 + * 一些检查和报错的宏 * */ @@ -14,7 +14,7 @@ #include #include -#include "common/colorstr.h" +#include "grt/common/colorstr.h" // GRT自定义报错信息 #define GRTRaiseError(ErrorMessage, ...) ({\ diff --git a/pygrt/C_extension/include/common/colorstr.h b/pygrt/C_extension/include/grt/common/colorstr.h similarity index 100% rename from pygrt/C_extension/include/common/colorstr.h rename to pygrt/C_extension/include/grt/common/colorstr.h diff --git a/pygrt/C_extension/include/common/const.h b/pygrt/C_extension/include/grt/common/const.h similarity index 98% rename from pygrt/C_extension/include/common/const.h rename to pygrt/C_extension/include/grt/common/const.h index 8a636f3..8d42966 100755 --- a/pygrt/C_extension/include/common/const.h +++ b/pygrt/C_extension/include/grt/common/const.h @@ -11,7 +11,7 @@ #include #include -#include "grt_error.h" +#include "grt/common/checkerror.h" // CMPLX macro not exist on MacOS #ifndef CMPLX @@ -109,8 +109,6 @@ typedef int MYINT; ///< 整数 };\ }) -#define GRT_SAFE_ - // ----------------------------------------------------------------------------- #define CHANNEL_NUM 3 ///< 3, 代码中分量个数(ZRT,ZNE) @@ -163,4 +161,4 @@ extern const char ZNEchs[]; * * @param[in] num_threads 线程数 */ -void set_num_threads(int num_threads); \ No newline at end of file +void grt_set_num_threads(int num_threads); \ No newline at end of file diff --git a/pygrt/C_extension/include/common/coord.h b/pygrt/C_extension/include/grt/common/coord.h similarity index 88% rename from pygrt/C_extension/include/common/coord.h rename to pygrt/C_extension/include/grt/common/coord.h index 9cc4482..2e2b461 100644 --- a/pygrt/C_extension/include/common/coord.h +++ b/pygrt/C_extension/include/grt/common/coord.h @@ -17,7 +17,7 @@ * @param[in] theta r轴相对x轴的旋转弧度(负数表示逆变换,即zrt->zxy) * @param[out] A 待旋转的矢量(s1, s2, s3) */ -void rot_zxy2zrt_vec(double theta, double A[3]); +void grt_rot_zxy2zrt_vec(double theta, double A[3]); @@ -27,7 +27,7 @@ void rot_zxy2zrt_vec(double theta, double A[3]); * @param[in] theta r轴相对x轴的旋转弧度(负数表示逆变换,即zrt->zxy) * @param[out] A 待旋转的二阶对称张量(s11, s12, s13, s22, s23, s33) */ -void rot_zxy2zrt_symtensor2odr(double theta, double A[6]); +void grt_rot_zxy2zrt_symtensor2odr(double theta, double A[6]); /** @@ -53,4 +53,4 @@ void rot_zxy2zrt_symtensor2odr(double theta, double A[6]); * @param[in,out] upar 柱坐标下的位移空间偏导 * @param[in] r r轴坐标 */ -void rot_zrt2zxy_upar(const double theta, double u[3], double upar[3][3], const double r); \ No newline at end of file +void grt_rot_zrt2zxy_upar(const double theta, double u[3], double upar[3][3], const double r); \ No newline at end of file diff --git a/pygrt/C_extension/include/common/dwm.h b/pygrt/C_extension/include/grt/common/dwm.h similarity index 89% rename from pygrt/C_extension/include/common/dwm.h rename to pygrt/C_extension/include/grt/common/dwm.h index d43c117..f40f992 100644 --- a/pygrt/C_extension/include/common/dwm.h +++ b/pygrt/C_extension/include/grt/common/dwm.h @@ -15,8 +15,8 @@ #include -#include "common/model.h" -#include "common/kernel.h" +#include "grt/common/model.h" +#include "grt/common/kernel.h" /** @@ -43,11 +43,11 @@ * * @return k 积分截至时的波数 */ -MYREAL discrete_integ( - const MODEL1D *mod1d, MYREAL dk, MYREAL kmax, MYREAL keps, MYCOMPLEX omega, +MYREAL grt_discrete_integ( + const GRT_MODEL1D *mod1d, MYREAL dk, MYREAL kmax, MYREAL keps, MYCOMPLEX omega, MYINT nr, MYREAL *rs, MYCOMPLEX sum_J[nr][SRC_M_NUM][INTEG_NUM], bool calc_upar, MYCOMPLEX sum_uiz_J[nr][SRC_M_NUM][INTEG_NUM], MYCOMPLEX sum_uir_J[nr][SRC_M_NUM][INTEG_NUM], - FILE *fstats, KernelFunc kerfunc, MYINT *stats); + FILE *fstats, GRT_KernelFunc kerfunc, MYINT *stats); diff --git a/pygrt/C_extension/include/common/fim.h b/pygrt/C_extension/include/grt/common/fim.h similarity index 87% rename from pygrt/C_extension/include/common/fim.h rename to pygrt/C_extension/include/grt/common/fim.h index 83b5a52..10e744c 100755 --- a/pygrt/C_extension/include/common/fim.h +++ b/pygrt/C_extension/include/grt/common/fim.h @@ -14,9 +14,9 @@ #include -#include "common/const.h" -#include "common/model.h" -#include "common/kernel.h" +#include "grt/common/const.h" +#include "grt/common/model.h" +#include "grt/common/kernel.h" @@ -50,13 +50,13 @@ * * @return k 积分截至时的波数 */ -MYREAL linear_filon_integ( - const MODEL1D *mod1d, MYREAL k0, MYREAL dk0, MYREAL filondk, MYREAL kmax, MYREAL keps, MYCOMPLEX omega, +MYREAL grt_linear_filon_integ( + const GRT_MODEL1D *mod1d, MYREAL k0, MYREAL dk0, MYREAL filondk, MYREAL kmax, MYREAL keps, MYCOMPLEX omega, MYINT nr, MYREAL *rs, MYCOMPLEX sum_J[nr][SRC_M_NUM][INTEG_NUM], bool calc_upar, MYCOMPLEX sum_uiz_J[nr][SRC_M_NUM][INTEG_NUM], MYCOMPLEX sum_uir_J[nr][SRC_M_NUM][INTEG_NUM], - FILE *fstats, KernelFunc kerfunc, MYINT *stats); + FILE *fstats, GRT_KernelFunc kerfunc, MYINT *stats); diff --git a/pygrt/C_extension/include/common/integral.h b/pygrt/C_extension/include/grt/common/integral.h similarity index 96% rename from pygrt/C_extension/include/common/integral.h rename to pygrt/C_extension/include/grt/common/integral.h index a15f971..e8e1388 100644 --- a/pygrt/C_extension/include/common/integral.h +++ b/pygrt/C_extension/include/grt/common/integral.h @@ -9,7 +9,7 @@ #pragma once -#include "common/const.h" +#include "grt/common/const.h" /** * 计算核函数和Bessel函数的乘积,相当于计算了一个小积分区间内的值。参数中涉及两种数组形状: @@ -24,7 +24,7 @@ * @param[out] SUM 该dk区间内的积分值 * */ -void int_Pk( +void grt_int_Pk( MYREAL k, MYREAL r, const MYCOMPLEX QWV[SRC_M_NUM][QWV_NUM], bool calc_uir, @@ -39,7 +39,7 @@ void int_Pk( * @param[in] sum_J 积分结果 * @param[out] tol Z、R、T分量结果 */ -void merge_Pk( +void grt_merge_Pk( const MYCOMPLEX sum_J[SRC_M_NUM][INTEG_NUM], MYCOMPLEX tol[SRC_M_NUM][CHANNEL_NUM]); @@ -58,7 +58,7 @@ void merge_Pk( * @param[out] SUM 该dk区间内的积分值 * */ -void int_Pk_filon( +void grt_int_Pk_filon( MYREAL k, MYREAL r, bool iscos, const MYCOMPLEX QWV[SRC_M_NUM][QWV_NUM], bool calc_uir, @@ -77,7 +77,7 @@ void int_Pk_filon( * @param[out] SUM 该三点区间内的积分值 * */ -void int_Pk_sa_filon( +void grt_int_Pk_sa_filon( const MYREAL k3[3], MYREAL r, const MYCOMPLEX QWV3[3][SRC_M_NUM][QWV_NUM], bool calc_uir, diff --git a/pygrt/C_extension/include/common/iostats.h b/pygrt/C_extension/include/grt/common/iostats.h similarity index 91% rename from pygrt/C_extension/include/common/iostats.h rename to pygrt/C_extension/include/grt/common/iostats.h index 20289db..56de1e8 100755 --- a/pygrt/C_extension/include/common/iostats.h +++ b/pygrt/C_extension/include/grt/common/iostats.h @@ -11,7 +11,7 @@ #include #include -#include "common/const.h" +#include "grt/common/const.h" @@ -26,7 +26,7 @@ * @note 文件记录的值均为波数积分的中间结果,与最终的结果还差一系列的系数, * 记录其值主要用于参考其变化趋势。 */ -void write_stats( +void grt_write_stats( FILE *f0, MYREAL k, const MYCOMPLEX QWV[SRC_M_NUM][QWV_NUM]); @@ -38,7 +38,7 @@ void write_stats( * * @return 0表示读取成功,-1表示读取结果/失败 */ -MYINT extract_stats(FILE *bf0, FILE *af0); +MYINT grt_extract_stats(FILE *bf0, FILE *af0); @@ -53,7 +53,7 @@ MYINT extract_stats(FILE *bf0, FILE *af0); * 记录其值主要用于参考其变化趋势。 * */ -void write_stats_ptam( +void grt_write_stats_ptam( FILE *f0, MYREAL Kpt[SRC_M_NUM][INTEG_NUM][PTAM_MAX_PT], MYCOMPLEX Fpt[SRC_M_NUM][INTEG_NUM][PTAM_MAX_PT]); @@ -67,5 +67,5 @@ void write_stats_ptam( * * @return 0表示读取成功,-1表示读取结果/失败 */ -MYINT extract_stats_ptam(FILE *bf0, FILE *af0); +MYINT grt_extract_stats_ptam(FILE *bf0, FILE *af0); diff --git a/pygrt/C_extension/include/common/kernel.h b/pygrt/C_extension/include/grt/common/kernel.h similarity index 68% rename from pygrt/C_extension/include/common/kernel.h rename to pygrt/C_extension/include/grt/common/kernel.h index 0ebac2f..1182919 100644 --- a/pygrt/C_extension/include/common/kernel.h +++ b/pygrt/C_extension/include/grt/common/kernel.h @@ -11,11 +11,11 @@ #pragma once -#include "common/model.h" +#include "grt/common/model.h" /** * 计算核函数的函数指针,动态与静态的接口一致 */ -typedef void (*KernelFunc) ( - const MODEL1D *mod1d, MYCOMPLEX omega, MYREAL k, MYCOMPLEX QWV[SRC_M_NUM][QWV_NUM], +typedef void (*GRT_KernelFunc) ( + const GRT_MODEL1D *mod1d, MYCOMPLEX omega, MYREAL k, MYCOMPLEX QWV[SRC_M_NUM][QWV_NUM], bool calc_uiz, MYCOMPLEX QWV_uiz[SRC_M_NUM][QWV_NUM], MYINT *stats); \ No newline at end of file diff --git a/pygrt/C_extension/include/common/logo.h b/pygrt/C_extension/include/grt/common/logo.h similarity index 94% rename from pygrt/C_extension/include/common/logo.h rename to pygrt/C_extension/include/grt/common/logo.h index ba223fd..5f14806 100644 --- a/pygrt/C_extension/include/common/logo.h +++ b/pygrt/C_extension/include/grt/common/logo.h @@ -12,16 +12,16 @@ #include -#include "common/const.h" -#include "common/version.h" -#include "common/colorstr.h" +#include "grt/common/const.h" +#include "grt/common/version.h" +#include "grt/common/colorstr.h" #if _TEST_WHETHER_WIN32_ #include #endif -inline GCC_ALWAYS_INLINE void print_logo(){ +inline GCC_ALWAYS_INLINE void grt_print_logo(){ #if _TEST_WHETHER_WIN32_ // 在Windows上设置控制台代码页为UTF-8,以显示下方由特殊字符组成的logo diff --git a/pygrt/C_extension/include/common/matrix.h b/pygrt/C_extension/include/grt/common/matrix.h similarity index 79% rename from pygrt/C_extension/include/common/matrix.h rename to pygrt/C_extension/include/grt/common/matrix.h index 963b9b8..5596b11 100755 --- a/pygrt/C_extension/include/common/matrix.h +++ b/pygrt/C_extension/include/grt/common/matrix.h @@ -8,7 +8,7 @@ #pragma once -#include "common/const.h" +#include "grt/common/const.h" /** * 计算2x2复矩阵的逆 @@ -17,7 +17,7 @@ * @param[out] invM 逆矩阵 * @param[out] stats 状态代码,是否有除零错误,非0为异常值 */ -inline GCC_ALWAYS_INLINE void cmat2x2_inv(const MYCOMPLEX M[2][2], MYCOMPLEX invM[2][2], MYINT *stats) { +inline GCC_ALWAYS_INLINE void grt_cmat2x2_inv(const MYCOMPLEX M[2][2], MYCOMPLEX invM[2][2], MYINT *stats) { MYCOMPLEX M00 = M[0][0]; MYCOMPLEX M01 = M[0][1]; MYCOMPLEX M10 = M[1][0]; @@ -46,7 +46,7 @@ inline GCC_ALWAYS_INLINE void cmat2x2_inv(const MYCOMPLEX M[2][2], MYCOMPLEX inv * @param[in] M2 矩阵2 * @param[out] M 和矩阵 */ -inline GCC_ALWAYS_INLINE void cmat2x2_add(const MYCOMPLEX M1[2][2], const MYCOMPLEX M2[2][2], MYCOMPLEX M[2][2]){ +inline GCC_ALWAYS_INLINE void grt_cmat2x2_add(const MYCOMPLEX M1[2][2], const MYCOMPLEX M2[2][2], MYCOMPLEX M[2][2]){ M[0][0] = M1[0][0] + M2[0][0]; M[0][1] = M1[0][1] + M2[0][1]; M[1][0] = M1[1][0] + M2[1][0]; @@ -60,7 +60,7 @@ inline GCC_ALWAYS_INLINE void cmat2x2_add(const MYCOMPLEX M1[2][2], const MYCOMP * @param[in] M2 矩阵2 * @param[out] M 差矩阵, M1 - M2 */ -inline GCC_ALWAYS_INLINE void cmat2x2_sub(const MYCOMPLEX M1[2][2], const MYCOMPLEX M2[2][2], MYCOMPLEX M[2][2]){ +inline GCC_ALWAYS_INLINE void grt_cmat2x2_sub(const MYCOMPLEX M1[2][2], const MYCOMPLEX M2[2][2], MYCOMPLEX M[2][2]){ M[0][0] = M1[0][0] - M2[0][0]; M[0][1] = M1[0][1] - M2[0][1]; M[1][0] = M1[1][0] - M2[1][0]; @@ -72,7 +72,7 @@ inline GCC_ALWAYS_INLINE void cmat2x2_sub(const MYCOMPLEX M1[2][2], const MYCOMP * * @param[in,out] M 差矩阵 I-M2 */ -inline GCC_ALWAYS_INLINE void cmat2x2_one_sub(MYCOMPLEX M[2][2]){ +inline GCC_ALWAYS_INLINE void grt_cmat2x2_one_sub(MYCOMPLEX M[2][2]){ M[0][0] = RONE - M[0][0]; M[0][1] = - M[0][1]; M[1][0] = - M[1][0]; @@ -86,7 +86,7 @@ inline GCC_ALWAYS_INLINE void cmat2x2_one_sub(MYCOMPLEX M[2][2]){ * @param[in] M2 矩阵2 * @param[out] M 积矩阵, M1 * M2 */ -inline GCC_ALWAYS_INLINE void cmat2x2_mul(const MYCOMPLEX M1[2][2], const MYCOMPLEX M2[2][2], MYCOMPLEX M[2][2]){ +inline GCC_ALWAYS_INLINE void grt_cmat2x2_mul(const MYCOMPLEX M1[2][2], const MYCOMPLEX M2[2][2], MYCOMPLEX M[2][2]){ MYCOMPLEX M011, M012, M021, M022; MYCOMPLEX M111, M112, M121, M122; M011 = M1[0][0]; M012 = M1[0][1]; @@ -106,7 +106,7 @@ inline GCC_ALWAYS_INLINE void cmat2x2_mul(const MYCOMPLEX M1[2][2], const MYCOMP * @param[in] k 常数 * @param[out] M 积矩阵, k * M2 */ -inline GCC_ALWAYS_INLINE void cmat2x2_k(const MYCOMPLEX M1[2][2], MYCOMPLEX k0, MYCOMPLEX M[2][2]){ +inline GCC_ALWAYS_INLINE void grt_cmat2x2_k(const MYCOMPLEX M1[2][2], MYCOMPLEX k0, MYCOMPLEX M[2][2]){ M[0][0] = M1[0][0] * k0; M[0][1] = M1[0][1] * k0; M[1][0] = M1[1][0] * k0; @@ -120,7 +120,7 @@ inline GCC_ALWAYS_INLINE void cmat2x2_k(const MYCOMPLEX M1[2][2], MYCOMPLEX k0, * @param[in] M2 向量 * @param[out] M 积矩阵, M1 * M2 */ -inline GCC_ALWAYS_INLINE void cmat2x1_mul(const MYCOMPLEX M1[2][2], const MYCOMPLEX M2[2], MYCOMPLEX M[2]){ +inline GCC_ALWAYS_INLINE void grt_cmat2x1_mul(const MYCOMPLEX M1[2][2], const MYCOMPLEX M2[2], MYCOMPLEX M[2]){ MYCOMPLEX M00, M10; M00 = M1[0][0]*M2[0] + M1[0][1]*M2[1]; M10 = M1[1][0]*M2[0] + M1[1][1]*M2[1]; @@ -134,7 +134,7 @@ inline GCC_ALWAYS_INLINE void cmat2x1_mul(const MYCOMPLEX M1[2][2], const MYCOMP * @param[in] M1 源矩阵 * @param[out] M2 目标矩阵 */ -inline GCC_ALWAYS_INLINE void cmat2x2_assign(const MYCOMPLEX M1[2][2], MYCOMPLEX M2[2][2]){ +inline GCC_ALWAYS_INLINE void grt_cmat2x2_assign(const MYCOMPLEX M1[2][2], MYCOMPLEX M2[2][2]){ M2[0][0] = M1[0][0]; M2[0][1] = M1[0][1]; M2[1][0] = M1[1][0]; @@ -151,7 +151,7 @@ inline GCC_ALWAYS_INLINE void cmat2x2_assign(const MYCOMPLEX M1[2][2], MYCOMPLEX * @param[in] M2 M2矩阵 * @param[out] M 积矩阵 M1 * M2 */ -inline GCC_ALWAYS_INLINE void cmatmxn_mul(MYINT m1, MYINT n1, MYINT p1, const MYCOMPLEX M1[m1][n1], const MYCOMPLEX M2[n1][p1], MYCOMPLEX M[m1][p1]){ +inline GCC_ALWAYS_INLINE void grt_cmatmxn_mul(MYINT m1, MYINT n1, MYINT p1, const MYCOMPLEX M1[m1][n1], const MYCOMPLEX M2[n1][p1], MYCOMPLEX M[m1][p1]){ MYINT m, n, k; MYCOMPLEX M0[m1][p1]; for(m=0; m #include -#include "common/const.h" +#include "grt/common/const.h" /** 单个水平层的结构体 */ typedef struct _LAYER { @@ -38,7 +38,7 @@ typedef struct _MODEL1D { MYINT imax; ///< 较深的虚拟层位,max(isrc, ircv) bool ircvup; ///< 接收点是否浅于震源层, ircv < isrc -} MODEL1D; +} GRT_MODEL1D; /** 用于从Python传递数据的中间结构体 */ @@ -57,7 +57,7 @@ typedef struct _PYMODEL1D { MYREAL *Rho; ///< Rho[n] 密度 g/cm^3 MYREAL *Qa; ///< Qa[n] P波Q值 MYREAL *Qb; ///< Qb[n] S波Q值 -} PYMODEL1D; +} GRT_PYMODEL1D; /** @@ -66,7 +66,7 @@ typedef struct _PYMODEL1D { * @param[in] mod1d `MODEL1D` 结构体指针 * */ -void print_mod1d(const MODEL1D *mod1d); +void grt_print_mod1d(const GRT_MODEL1D *mod1d); /** * 打印PYMODEL1D模型参数信息,主要用于调试程序 @@ -74,14 +74,14 @@ void print_mod1d(const MODEL1D *mod1d); * @param[in] pymod `PYMODEL1D` 结构体指针 * */ -void print_pymod(const PYMODEL1D *pymod); +void grt_print_pymod(const GRT_PYMODEL1D *pymod); /** * 释放 `PYMODEL1D` 结构体指针 * * @param[out] pymod `PYMODEL1D` 结构体指针 */ -void free_pymod(PYMODEL1D *pymod); +void grt_free_pymod(GRT_PYMODEL1D *pymod); /** * 初始化模型内存空间 @@ -91,7 +91,7 @@ void free_pymod(PYMODEL1D *pymod); * @return `MODEL1D` 结构体指针 * */ -MODEL1D * init_mod1d(MYINT n); +GRT_MODEL1D * grt_init_mod1d(MYINT n); /** * 将 `PYMODEL1D` 结构体信息转到 `MODEL1D` 结构体中 @@ -100,7 +100,7 @@ MODEL1D * init_mod1d(MYINT n); * @param[out] mod1d `MODEL1D` 结构体指针 * */ -void get_mod1d(const PYMODEL1D *pymod1d, MODEL1D *mod1d); +void grt_get_mod1d(const GRT_PYMODEL1D *pymod1d, GRT_MODEL1D *mod1d); /** * 复制 `MODEL1D` 结构体,要求都以初始化 @@ -109,14 +109,14 @@ void get_mod1d(const PYMODEL1D *pymod1d, MODEL1D *mod1d); * @param[out] mod1d2 `MODEL1D` 目标结构体指针 * */ -void copy_mod1d(const MODEL1D *mod1d1, MODEL1D *mod1d2); +void grt_copy_mod1d(const GRT_MODEL1D *mod1d1, GRT_MODEL1D *mod1d2); /** * 释放 `MODEL1D` 结构体指针 * * @param[out] mod1d `MODEL1D` 结构体指针 */ -void free_mod1d(MODEL1D *mod1d); +void grt_free_mod1d(GRT_MODEL1D *mod1d); /** @@ -125,7 +125,7 @@ void free_mod1d(MODEL1D *mod1d); * @param[in,out] mod1d `MODEL1D` 结构体指针 * @param[in] omega 复数频率 */ -void update_mod1d_omega(MODEL1D *mod1d, MYCOMPLEX omega); +void grt_update_mod1d_omega(GRT_MODEL1D *mod1d, MYCOMPLEX omega); /** * 初始化PYMODEL1D模型内存空间 @@ -135,7 +135,7 @@ void update_mod1d_omega(MODEL1D *mod1d, MYCOMPLEX omega); * @return `PYMODEL1D` 结构体指针 * */ -PYMODEL1D * init_pymod(MYINT n); +GRT_PYMODEL1D * grt_init_pymod(MYINT n); /** * 从文件中读取模型文件,以PYMODEL1D结构体形式 @@ -149,7 +149,7 @@ PYMODEL1D * init_pymod(MYINT n); * @return `PYMODEL1D` 结构体指针 * */ -PYMODEL1D * read_pymod_from_file(const char *command, const char *modelpath, double depsrc, double deprcv, bool allowLiquid); +GRT_PYMODEL1D * grt_read_pymod_from_file(const char *command, const char *modelpath, double depsrc, double deprcv, bool allowLiquid); /** @@ -160,4 +160,4 @@ PYMODEL1D * read_pymod_from_file(const char *command, const char *modelpath, dou * @param vmax (out)最大速度 * */ -void get_pymod_vmin_vmax(const PYMODEL1D *pymod, double *vmin, double *vmax); \ No newline at end of file +void grt_get_pymod_vmin_vmax(const GRT_PYMODEL1D *pymod, double *vmin, double *vmax); \ No newline at end of file diff --git a/pygrt/C_extension/include/common/myfftw.h b/pygrt/C_extension/include/grt/common/myfftw.h similarity index 78% rename from pygrt/C_extension/include/common/myfftw.h rename to pygrt/C_extension/include/grt/common/myfftw.h index 26b52d0..af1c0db 100644 --- a/pygrt/C_extension/include/common/myfftw.h +++ b/pygrt/C_extension/include/grt/common/myfftw.h @@ -17,7 +17,7 @@ #include #include -#include "common/const.h" +#include "grt/common/const.h" #define FOR_EACH_FFTW_TYPE \ @@ -55,7 +55,7 @@ typedef struct {\ MYREAL f0; \ /*最终不执行FFTW而使用最朴素的傅里叶逆变换*/\ bool naive_inv; \ -} FFTW##S##_HOLDER;\ +} GRT_FFTW##S##_HOLDER;\ \ /** * 初始化 FFTW_HOLDER 结构体指针,进行复数频谱到实数序列的逆变换 @@ -68,22 +68,22 @@ typedef struct {\ * @return FFTW_HOLDER 结构体指针 * */\ -FFTW##S##_HOLDER * create_fftw##s##_holder_C2R_1D(const MYINT nt, const MYREAL dt, const MYINT nf_valid, const MYREAL df); \ +GRT_FFTW##S##_HOLDER * grt_create_fftw##s##_holder_C2R_1D(const MYINT nt, const MYREAL dt, const MYINT nf_valid, const MYREAL df); \ /** 初始化 FFTW_HOLDER 结构体指针,进行实数序列到复数频谱的正变换 */\ -FFTW##S##_HOLDER * create_fftw##s##_holder_R2C_1D(const MYINT nt, const MYREAL dt, const MYINT nf_valid, const MYREAL df); \ +GRT_FFTW##S##_HOLDER * grt_create_fftw##s##_holder_R2C_1D(const MYINT nt, const MYREAL dt, const MYINT nf_valid, const MYREAL df); \ \ /** 将内部数据全部置零 */\ -void reset_fftw##s##_holder_zero(FFTW##S##_HOLDER *fh);\ +void grt_reset_fftw##s##_holder_zero(GRT_FFTW##S##_HOLDER *fh);\ \ /** 清理函数:释放结构体内存,防止内存泄漏 */\ -void destroy_fftw##s##_holder(FFTW##S##_HOLDER *fh);\ +void grt_destroy_fftw##s##_holder(GRT_FFTW##S##_HOLDER *fh);\ \ \ /** * 最朴素的非均匀频域采样逆变换 * 严格按照傅里叶逆变换定义:g(t) = \int G(f) * e^(i2πft) df */ \ -void naive_inverse_transform_##T(FFTW##S##_HOLDER *fh);\ +void grt_naive_inverse_transform_##T(GRT_FFTW##S##_HOLDER *fh);\ FOR_EACH_FFTW_TYPE #undef X diff --git a/pygrt/C_extension/include/common/progressbar.h b/pygrt/C_extension/include/grt/common/progressbar.h similarity index 80% rename from pygrt/C_extension/include/common/progressbar.h rename to pygrt/C_extension/include/grt/common/progressbar.h index 45274f8..85dee83 100755 --- a/pygrt/C_extension/include/common/progressbar.h +++ b/pygrt/C_extension/include/grt/common/progressbar.h @@ -7,7 +7,7 @@ * */ -#include "common/const.h" +#include "grt/common/const.h" #define _PROGRESSBAR_WIDTH_ 45 ///< 定义进度条的长度 @@ -19,4 +19,4 @@ * @param[in] prefix 进度条前缀字符串 * @param[in] percentage 百分比(整数) */ -void printprogressBar(const char *prefix, MYINT percentage); \ No newline at end of file +void grt_printprogressBar(const char *prefix, MYINT percentage); \ No newline at end of file diff --git a/pygrt/C_extension/include/common/prtdbg.h b/pygrt/C_extension/include/grt/common/prtdbg.h similarity index 100% rename from pygrt/C_extension/include/common/prtdbg.h rename to pygrt/C_extension/include/grt/common/prtdbg.h diff --git a/pygrt/C_extension/include/common/ptam.h b/pygrt/C_extension/include/grt/common/ptam.h similarity index 90% rename from pygrt/C_extension/include/common/ptam.h rename to pygrt/C_extension/include/grt/common/ptam.h index 7d47a4c..30d48b2 100755 --- a/pygrt/C_extension/include/common/ptam.h +++ b/pygrt/C_extension/include/grt/common/ptam.h @@ -18,8 +18,8 @@ #include -#include "common/model.h" -#include "common/kernel.h" +#include "grt/common/model.h" +#include "grt/common/kernel.h" @@ -46,14 +46,14 @@ * * */ -void PTA_method( - const MODEL1D *mod1d, MYREAL k0, MYREAL predk, MYCOMPLEX omega, +void grt_PTA_method( + const GRT_MODEL1D *mod1d, MYREAL k0, MYREAL predk, MYCOMPLEX omega, MYINT nr, MYREAL *rs, MYCOMPLEX sum_J0[nr][SRC_M_NUM][INTEG_NUM], bool calc_upar, MYCOMPLEX sum_uiz_J0[nr][SRC_M_NUM][INTEG_NUM], MYCOMPLEX sum_uir_J0[nr][SRC_M_NUM][INTEG_NUM], - FILE *ptam_fstatsnr[nr][2], KernelFunc kerfunc, MYINT *stats); + FILE *ptam_fstatsnr[nr][2], GRT_KernelFunc kerfunc, MYINT *stats); @@ -73,7 +73,7 @@ void PTA_method( * @return 波峰(1),波谷(-1),其它(0) * */ -MYINT cplx_peak_or_trough( +MYINT grt_cplx_peak_or_trough( MYINT idx1, MYINT idx2, const MYCOMPLEX arr[PTAM_WINDOW_SIZE][SRC_M_NUM][INTEG_NUM], MYREAL k, MYREAL dk, MYREAL *pk, MYCOMPLEX *value); @@ -88,4 +88,4 @@ MYINT cplx_peak_or_trough( * @param[in,out] arr 振荡的数组,最终收敛值在第一个,arr[0] * */ -void cplx_shrink(MYINT n1, MYCOMPLEX *arr); \ No newline at end of file +void grt_cplx_shrink(MYINT n1, MYCOMPLEX *arr); \ No newline at end of file diff --git a/pygrt/C_extension/include/common/quadratic.h b/pygrt/C_extension/include/grt/common/quadratic.h similarity index 80% rename from pygrt/C_extension/include/common/quadratic.h rename to pygrt/C_extension/include/grt/common/quadratic.h index d1c325a..5531185 100755 --- a/pygrt/C_extension/include/common/quadratic.h +++ b/pygrt/C_extension/include/grt/common/quadratic.h @@ -12,7 +12,7 @@ #pragma once -#include "common/const.h" +#include "grt/common/const.h" /** @@ -25,7 +25,7 @@ * @param[out] pb 拟合b值 * @param[out] pc 拟合c值 */ -void quad_term(const MYREAL x[3], const MYCOMPLEX f[3], MYCOMPLEX *pa, MYCOMPLEX *pb, MYCOMPLEX *pc); +void grt_quad_term(const MYREAL x[3], const MYCOMPLEX f[3], MYCOMPLEX *pa, MYCOMPLEX *pb, MYCOMPLEX *pc); /** @@ -38,7 +38,7 @@ void quad_term(const MYREAL x[3], const MYCOMPLEX f[3], MYCOMPLEX *pa, MYCOMPLEX * @return \f$ f(x) = ax^2 + bx + c \f$ * */ -MYCOMPLEX quad_eval(MYREAL x, MYCOMPLEX a, MYCOMPLEX b, MYCOMPLEX c); +MYCOMPLEX grt_quad_eval(MYREAL x, MYCOMPLEX a, MYCOMPLEX b, MYCOMPLEX c); /** @@ -52,4 +52,4 @@ MYCOMPLEX quad_eval(MYREAL x, MYCOMPLEX a, MYCOMPLEX b, MYCOMPLEX c); * @return \f$ \frac{1}{3}a(x_2^3-x_1^3) + \frac{1}{2}b(x_2^2-x_1^2) + c(x_2-x_1) \f$ * */ -MYCOMPLEX quad_integral(MYREAL x1, MYREAL x2, MYCOMPLEX a, MYCOMPLEX b, MYCOMPLEX c); \ No newline at end of file +MYCOMPLEX grt_quad_integral(MYREAL x1, MYREAL x2, MYCOMPLEX a, MYCOMPLEX b, MYCOMPLEX c); \ No newline at end of file diff --git a/pygrt/C_extension/include/common/radiation.h b/pygrt/C_extension/include/grt/common/radiation.h similarity index 96% rename from pygrt/C_extension/include/common/radiation.h rename to pygrt/C_extension/include/grt/common/radiation.h index abce6d2..0d0dc0a 100644 --- a/pygrt/C_extension/include/common/radiation.h +++ b/pygrt/C_extension/include/grt/common/radiation.h @@ -11,7 +11,7 @@ #include -#include "common/const.h" +#include "grt/common/const.h" #define GRT_SYN_COMPUTE_EX 0 ///< 计算爆炸源 #define GRT_SYN_COMPUTE_SF 1 ///< 计算单力源 @@ -33,7 +33,7 @@ * 对于剪切源,mchn={strike, dip, rake}, * 对于张量源,mchn={Mxx, Mxy, Mxz, Myy, Myz, Mzz} */ -void set_source_radiation( +void grt_set_source_radiation( double srcRadi[SRC_M_NUM][CHANNEL_NUM], const int computeType, const bool par_theta, const double M0, const double coef, const double azrad, const double mchn[MECHANISM_NUM] ); \ No newline at end of file diff --git a/pygrt/C_extension/include/common/recursion.h b/pygrt/C_extension/include/grt/common/recursion.h similarity index 96% rename from pygrt/C_extension/include/common/recursion.h rename to pygrt/C_extension/include/grt/common/recursion.h index 8cf62da..698c1df 100644 --- a/pygrt/C_extension/include/common/recursion.h +++ b/pygrt/C_extension/include/grt/common/recursion.h @@ -13,7 +13,7 @@ #pragma once -#include "common/const.h" +#include "grt/common/const.h" // 对4个矩阵赋值 @@ -46,7 +46,7 @@ /** 合并 recursion_RD_PSV(SH) */ -void recursion_RD( +void grt_recursion_RD( const MYCOMPLEX RD1[2][2], MYCOMPLEX RDL1, const MYCOMPLEX RU1[2][2], MYCOMPLEX RUL1, const MYCOMPLEX TD1[2][2], MYCOMPLEX TDL1, const MYCOMPLEX TU1[2][2], MYCOMPLEX TUL1, const MYCOMPLEX RD2[2][2], MYCOMPLEX RDL2, @@ -65,7 +65,7 @@ void recursion_RD( * @param[out] stats 状态代码,是否有除零错误,非0为异常值 * */ -void recursion_RD_PSV( +void grt_recursion_RD_PSV( const MYCOMPLEX RD1[2][2], const MYCOMPLEX RU1[2][2], const MYCOMPLEX TD1[2][2], const MYCOMPLEX TU1[2][2], const MYCOMPLEX RD2[2][2], @@ -84,14 +84,14 @@ void recursion_RD_PSV( * @param[out] stats 状态代码,是否有除零错误,非0为异常值 * */ -void recursion_RD_SH( +void grt_recursion_RD_SH( MYCOMPLEX RDL1, MYCOMPLEX RUL1, MYCOMPLEX TDL1, MYCOMPLEX TUL1, MYCOMPLEX RDL2, MYCOMPLEX *RDL, MYCOMPLEX *invT, MYINT *stats); /** 合并 recursion_TD_PSV(SH) */ -void recursion_TD( +void grt_recursion_TD( const MYCOMPLEX RU1[2][2], MYCOMPLEX RUL1, const MYCOMPLEX TD1[2][2], MYCOMPLEX TDL1, const MYCOMPLEX RD2[2][2], MYCOMPLEX RDL2, @@ -110,7 +110,7 @@ void recursion_TD( * @param[out] stats 状态代码,是否有除零错误,非0为异常值 * */ -void recursion_TD_PSV( +void grt_recursion_TD_PSV( const MYCOMPLEX RU1[2][2], const MYCOMPLEX TD1[2][2], const MYCOMPLEX RD2[2][2], const MYCOMPLEX TD2[2][2], MYCOMPLEX TD[2][2], MYCOMPLEX inv_2x2T[2][2], MYINT *stats); @@ -127,14 +127,14 @@ void recursion_TD_PSV( * @param[out] stats 状态代码,是否有除零错误,非0为异常值 * */ -void recursion_TD_SH( +void grt_recursion_TD_SH( MYCOMPLEX RUL1, MYCOMPLEX TDL1, MYCOMPLEX RDL2, MYCOMPLEX TDL2, MYCOMPLEX *TDL, MYCOMPLEX *invT, MYINT *stats); /** 合并 recursion_RU_PSV(SH) */ -void recursion_RU( +void grt_recursion_RU( const MYCOMPLEX RU1[2][2], MYCOMPLEX RUL1, const MYCOMPLEX RD2[2][2], MYCOMPLEX RDL2, const MYCOMPLEX RU2[2][2], MYCOMPLEX RUL2, const MYCOMPLEX TD2[2][2], MYCOMPLEX TDL2, const MYCOMPLEX TU2[2][2], MYCOMPLEX TUL2, @@ -153,7 +153,7 @@ void recursion_RU( * @param[out] stats 状态代码,是否有除零错误,非0为异常值 * */ -void recursion_RU_PSV( +void grt_recursion_RU_PSV( const MYCOMPLEX RU1[2][2], const MYCOMPLEX RD2[2][2], const MYCOMPLEX RU2[2][2], const MYCOMPLEX TD2[2][2], const MYCOMPLEX TU2[2][2], @@ -172,7 +172,7 @@ void recursion_RU_PSV( * @param[out] stats 状态代码,是否有除零错误,非0为异常值 * */ -void recursion_RU_SH( +void grt_recursion_RU_SH( MYCOMPLEX RUL1, MYCOMPLEX RDL2, MYCOMPLEX RUL2, MYCOMPLEX TDL2, MYCOMPLEX TUL2, @@ -180,7 +180,7 @@ void recursion_RU_SH( /** 合并 recursion_TU_PSV(SH) */ -void recursion_TU( +void grt_recursion_TU( const MYCOMPLEX RU1[2][2], MYCOMPLEX RUL1, const MYCOMPLEX TU1[2][2], MYCOMPLEX TUL1, const MYCOMPLEX RD2[2][2], MYCOMPLEX RDL2, @@ -199,7 +199,7 @@ void recursion_TU( * @param[out] stats 状态代码,是否有除零错误,非0为异常值 * */ -void recursion_TU_PSV( +void grt_recursion_TU_PSV( const MYCOMPLEX RU1[2][2], const MYCOMPLEX TU1[2][2], const MYCOMPLEX RD2[2][2], const MYCOMPLEX TU2[2][2], MYCOMPLEX TU[2][2], MYCOMPLEX inv_2x2T[2][2], MYINT *stats); @@ -216,7 +216,7 @@ void recursion_TU_PSV( * @param[out] stats 状态代码,是否有除零错误,非0为异常值 * */ -void recursion_TU_SH( +void grt_recursion_TU_SH( MYCOMPLEX RUL1, MYCOMPLEX TUL1, MYCOMPLEX RDL2, MYCOMPLEX TUL2, MYCOMPLEX *TUL, MYCOMPLEX *invT, MYINT *stats); @@ -225,7 +225,7 @@ void recursion_TU_SH( /** 合并 recursion_RT_PSV(SH) */ -void recursion_RT( +void grt_recursion_RT( const MYCOMPLEX RD1[2][2], MYCOMPLEX RDL1, const MYCOMPLEX RU1[2][2], MYCOMPLEX RUL1, const MYCOMPLEX TD1[2][2], MYCOMPLEX TDL1, const MYCOMPLEX TU1[2][2], MYCOMPLEX TUL1, const MYCOMPLEX RD2[2][2], MYCOMPLEX RDL2, const MYCOMPLEX RU2[2][2], MYCOMPLEX RUL2, @@ -252,7 +252,7 @@ void recursion_RT( * @param[out] stats 状态代码,是否有除零错误,非0为异常值 * */ -void recursion_RT_PSV( +void grt_recursion_RT_PSV( const MYCOMPLEX RD1[2][2], const MYCOMPLEX RU1[2][2], const MYCOMPLEX TD1[2][2], const MYCOMPLEX TU1[2][2], const MYCOMPLEX RD2[2][2], const MYCOMPLEX RU2[2][2], @@ -279,7 +279,7 @@ void recursion_RT_PSV( * @param[out] stats 状态代码,是否有除零错误,非0为异常值 * */ -void recursion_RT_SH( +void grt_recursion_RT_SH( MYCOMPLEX RDL1, MYCOMPLEX RUL1, MYCOMPLEX TDL1, MYCOMPLEX TUL1, MYCOMPLEX RDL2, MYCOMPLEX RUL2, @@ -290,7 +290,7 @@ void recursion_RT_SH( /** 合并 recursion_RT_PSV(SH)_imaginary */ -void recursion_RT_imaginary( +void grt_recursion_RT_imaginary( MYCOMPLEX xa1, MYCOMPLEX xb1, MYREAL thk, MYREAL k, // 使用上层的厚度 MYCOMPLEX RU[2][2], MYCOMPLEX *RUL, MYCOMPLEX TD[2][2], MYCOMPLEX *TDL, MYCOMPLEX TU[2][2], MYCOMPLEX *TUL); @@ -306,7 +306,7 @@ void recursion_RT_imaginary( * @param[in,out] TD 上层 P-SV 下传透射系数矩阵 * @param[in,out] TU 上层 P-SV 上传透射系数矩阵 */ -void recursion_RT_PSV_imaginary( +void grt_recursion_RT_PSV_imaginary( MYCOMPLEX xa1, MYCOMPLEX xb1, MYREAL thk, MYREAL k, // 使用上层的厚度 MYCOMPLEX RU[2][2], MYCOMPLEX TD[2][2], MYCOMPLEX TU[2][2]); @@ -320,7 +320,7 @@ void recursion_RT_PSV_imaginary( * @param[in,out] TDL 上层 SH 下传透射系数 * @param[in,out] TUL 上层 SH 上传透射系数 */ -void recursion_RT_SH_imaginary( +void grt_recursion_RT_SH_imaginary( MYCOMPLEX xb1, MYREAL thk, MYREAL k, // 使用上层的厚度 MYCOMPLEX *RUL, MYCOMPLEX *TDL, MYCOMPLEX *TUL); @@ -358,7 +358,7 @@ void recursion_RT_SH_imaginary( * @param[in] coef_SH SH 波震源系数,\f$ SH_m \f$ ,维度2表示下行波(p=0)和上行波(p=1) * @param[out] qwv 最终通过矩阵传播计算出的在台站位置的\f$ q_m,w_m,v_m\f$ */ -void get_qwv( +void grt_get_qwv( bool ircvup, const MYCOMPLEX R1[2][2], MYCOMPLEX RL1, const MYCOMPLEX R2[2][2], MYCOMPLEX RL2, diff --git a/pygrt/C_extension/include/common/sacio.h b/pygrt/C_extension/include/grt/common/sacio.h similarity index 100% rename from pygrt/C_extension/include/common/sacio.h rename to pygrt/C_extension/include/grt/common/sacio.h diff --git a/pygrt/C_extension/include/common/sacio2.h b/pygrt/C_extension/include/grt/common/sacio2.h similarity index 78% rename from pygrt/C_extension/include/common/sacio2.h rename to pygrt/C_extension/include/grt/common/sacio2.h index a5789b8..1ab43f4 100644 --- a/pygrt/C_extension/include/common/sacio2.h +++ b/pygrt/C_extension/include/grt/common/sacio2.h @@ -9,7 +9,7 @@ #pragma once -#include "sacio.h" +#include "grt/common/sacio.h" /** * 读取SAC头段变量 @@ -18,7 +18,7 @@ * @param name (in)SAC文件路径 * @param hd (out)SAC头段变量结构体 */ -void read_SAC_HEAD(const char *command, const char *name, SACHEAD *hd); +void grt_read_SAC_HEAD(const char *command, const char *name, SACHEAD *hd); /** @@ -31,4 +31,4 @@ void read_SAC_HEAD(const char *command, const char *name, SACHEAD *hd); * * @return 浮点数指针 */ -float * read_SAC(const char *command, const char *name, SACHEAD *hd, float *arrout); \ No newline at end of file +float * grt_read_SAC(const char *command, const char *name, SACHEAD *hd, float *arrout); \ No newline at end of file diff --git a/pygrt/C_extension/include/common/safim.h b/pygrt/C_extension/include/grt/common/safim.h similarity index 87% rename from pygrt/C_extension/include/common/safim.h rename to pygrt/C_extension/include/grt/common/safim.h index 93ddbb6..38b788d 100755 --- a/pygrt/C_extension/include/common/safim.h +++ b/pygrt/C_extension/include/grt/common/safim.h @@ -15,9 +15,9 @@ #include -#include "common/const.h" -#include "common/model.h" -#include "common/kernel.h" +#include "grt/common/const.h" +#include "grt/common/model.h" +#include "grt/common/kernel.h" @@ -51,13 +51,13 @@ * * @return k 积分截至时的波数 */ -MYREAL sa_filon_integ( - const MODEL1D *mod1d, MYREAL vmin, MYREAL k0, MYREAL dk0, MYREAL tol, MYREAL kmax, MYCOMPLEX omega, +MYREAL grt_sa_filon_integ( + const GRT_MODEL1D *mod1d, MYREAL vmin, MYREAL k0, MYREAL dk0, MYREAL tol, MYREAL kmax, MYCOMPLEX omega, MYINT nr, MYREAL *rs, MYCOMPLEX sum_J0[nr][SRC_M_NUM][INTEG_NUM], bool calc_upar, MYCOMPLEX sum_uiz_J0[nr][SRC_M_NUM][INTEG_NUM], MYCOMPLEX sum_uir_J0[nr][SRC_M_NUM][INTEG_NUM], - FILE *fstats, KernelFunc kerfunc, MYINT *stats); + FILE *fstats, GRT_KernelFunc kerfunc, MYINT *stats); diff --git a/pygrt/C_extension/include/common/search.h b/pygrt/C_extension/include/grt/common/search.h similarity index 84% rename from pygrt/C_extension/include/common/search.h rename to pygrt/C_extension/include/grt/common/search.h index e0d02d2..2143eb2 100644 --- a/pygrt/C_extension/include/common/search.h +++ b/pygrt/C_extension/include/grt/common/search.h @@ -9,7 +9,7 @@ #pragma once #include -#include "common/const.h" +#include "grt/common/const.h" /** * 该函数对输入的整数数组进行线性搜索,找到目标值时返回其索引。 @@ -24,7 +24,7 @@ * @note 如果数组中存在多个目标值,该函数返回第一个匹配的索引。 * */ -MYINT findElement_MYINT(const MYINT array[], MYINT size, MYINT target); +MYINT grt_findElement_MYINT(const MYINT array[], MYINT size, MYINT target); /** * 搜索浮点数数组中最接近目标值且小于目标值的索引。 @@ -39,7 +39,7 @@ MYINT findElement_MYINT(const MYINT array[], MYINT size, MYINT target); * @note 如果数组中存在多个目标值,该函数返回第一个匹配的索引。 * */ -MYINT findLessEqualClosest_MYREAL(const MYREAL array[], MYINT size, MYREAL target); +MYINT grt_findLessEqualClosest_MYREAL(const MYREAL array[], MYINT size, MYREAL target); /** * 搜索浮点数数组中最接近目标值的索引。 @@ -53,7 +53,7 @@ MYINT findLessEqualClosest_MYREAL(const MYREAL array[], MYINT size, MYREAL targe * @note 如果数组中存在多个目标值,该函数返回第一个匹配的索引。 * */ -MYINT findClosest_MYREAL(const MYREAL array[], MYINT size, MYREAL target); +MYINT grt_findClosest_MYREAL(const MYREAL array[], MYINT size, MYREAL target); /** * 搜索浮点数数组的最大或最小值,返回其索引。 @@ -67,4 +67,4 @@ MYINT findClosest_MYREAL(const MYREAL array[], MYINT size, MYREAL target); * @note 如果数组中存在相同最值,该函数返回第一个匹配的索引。 * */ -MYINT findMinMax_MYREAL(const MYREAL array[], MYINT size, bool isMax); \ No newline at end of file +MYINT grt_findMinMax_MYREAL(const MYREAL array[], MYINT size, bool isMax); \ No newline at end of file diff --git a/pygrt/C_extension/include/common/util.h b/pygrt/C_extension/include/grt/common/util.h similarity index 88% rename from pygrt/C_extension/include/common/util.h rename to pygrt/C_extension/include/grt/common/util.h index bb9de7b..dc3b9b7 100644 --- a/pygrt/C_extension/include/common/util.h +++ b/pygrt/C_extension/include/grt/common/util.h @@ -11,9 +11,9 @@ #include -#include "common/const.h" -#include "common/model.h" -#include "common/myfftw.h" +#include "grt/common/const.h" +#include "grt/common/model.h" +#include "grt/common/myfftw.h" /** * 指定分隔符,从一串字符串中分割出子字符串数组 @@ -24,7 +24,7 @@ * * @return split 子字符串数组 */ -char ** string_split(const char *string, const char *delim, int *size); +char ** grt_string_split(const char *string, const char *delim, int *size); /** @@ -34,7 +34,7 @@ char ** string_split(const char *string, const char *delim, int *size); * * @return 指向最后一项字符串的指针 */ -const char* get_basename(const char* path); +const char* grt_get_basename(const char* path); @@ -69,10 +69,10 @@ const char* get_basename(const char* path); * @param[in] grn_uir 格林函数对r偏导频谱结果 * */ -void GF_freq2time_write_to_file( - const char *command, const PYMODEL1D *pymod, +void grt_GF_freq2time_write_to_file( + const char *command, const GRT_PYMODEL1D *pymod, const char *s_output_dir, const char *s_modelname, const char *s_depsrc, const char *s_deprcv, - const MYREAL wI, FFTW_HOLDER *pt_fh, + const MYREAL wI, GRT_FFTW_HOLDER *pt_fh, const MYINT nr, char *s_dists[nr], const MYREAL dists[nr], MYREAL travtPS[nr][2], const MYREAL depsrc, const MYREAL deprcv, const MYREAL delayT0, const MYREAL delayV0, const bool calc_upar, diff --git a/pygrt/C_extension/include/common/version.h b/pygrt/C_extension/include/grt/common/version.h similarity index 100% rename from pygrt/C_extension/include/common/version.h rename to pygrt/C_extension/include/grt/common/version.h diff --git a/pygrt/C_extension/include/dynamic/grn.h b/pygrt/C_extension/include/grt/dynamic/grn.h similarity index 96% rename from pygrt/C_extension/include/dynamic/grn.h rename to pygrt/C_extension/include/grt/dynamic/grn.h index b9a6610..414315d 100755 --- a/pygrt/C_extension/include/dynamic/grn.h +++ b/pygrt/C_extension/include/grt/dynamic/grn.h @@ -14,7 +14,7 @@ #pragma once -#include "common/model.h" +#include "grt/common/model.h" /** @@ -48,8 +48,8 @@ * @param[in] statsidxs 特定频点的索引值 * */ -void integ_grn_spec( - PYMODEL1D *pymod1d, MYINT nf1, MYINT nf2, MYREAL *freqs, +void grt_integ_grn_spec( + GRT_PYMODEL1D *pymod1d, MYINT nf1, MYINT nf2, MYREAL *freqs, MYINT nr, MYREAL *rs, MYREAL wI, MYREAL vmin_ref, MYREAL keps, MYREAL ampk, MYREAL k0, MYREAL Length, MYREAL filonLength, MYREAL safilonTol, MYREAL filonCut, bool print_progressbar, diff --git a/pygrt/C_extension/include/dynamic/layer.h b/pygrt/C_extension/include/grt/dynamic/layer.h similarity index 92% rename from pygrt/C_extension/include/dynamic/layer.h rename to pygrt/C_extension/include/grt/dynamic/layer.h index 48d22b2..b988e1d 100755 --- a/pygrt/C_extension/include/dynamic/layer.h +++ b/pygrt/C_extension/include/grt/dynamic/layer.h @@ -13,8 +13,8 @@ #pragma once -#include "common/model.h" -#include "common/const.h" +#include "grt/common/model.h" +#include "grt/common/const.h" /** * 计算自由表面的 P-SV 波反射系数,公式(5.3.10-14) @@ -27,7 +27,7 @@ * @param[out] stats 状态代码,是否有除零错误,非0为异常值 * */ -void calc_R_tilt_PSV(MYCOMPLEX xa0, MYCOMPLEX xb0, MYCOMPLEX kbkb0, MYREAL k, MYCOMPLEX R_tilt[2][2], MYINT *stats); +void grt_calc_R_tilt_PSV(MYCOMPLEX xa0, MYCOMPLEX xb0, MYCOMPLEX kbkb0, MYREAL k, MYCOMPLEX R_tilt[2][2], MYINT *stats); /** @@ -41,7 +41,7 @@ void calc_R_tilt_PSV(MYCOMPLEX xa0, MYCOMPLEX xb0, MYCOMPLEX kbkb0, MYREAL k, MY * @param[out] R_EV P-SV接收函数矩阵 * */ -void calc_R_EV_PSV( +void grt_calc_R_EV_PSV( MYCOMPLEX xa_rcv, MYCOMPLEX xb_rcv, bool ircvup, MYREAL k, const MYCOMPLEX R[2][2], MYCOMPLEX R_EV[2][2]); @@ -55,7 +55,7 @@ void calc_R_EV_PSV( * @param[out] R_EVL SH接收函数值 * */ -void calc_R_EV_SH( +void grt_calc_R_EV_SH( MYCOMPLEX xb_rcv, MYREAL k, MYCOMPLEX RL, MYCOMPLEX *R_EVL); @@ -73,7 +73,7 @@ void calc_R_EV_SH( * @param[out] R_EV P-SV接收函数矩阵 * */ -void calc_uiz_R_EV_PSV( +void grt_calc_uiz_R_EV_PSV( MYCOMPLEX xa_rcv, MYCOMPLEX xb_rcv, bool ircvup, MYREAL k, const MYCOMPLEX R[2][2], MYCOMPLEX R_EV[2][2]); @@ -90,7 +90,7 @@ void calc_uiz_R_EV_PSV( * @param[out] R_EVL SH接收函数值 * */ -void calc_uiz_R_EV_SH( +void grt_calc_uiz_R_EV_SH( MYCOMPLEX xb_rcv, bool ircvup, MYREAL k, MYCOMPLEX RL, MYCOMPLEX *R_EVL); @@ -122,7 +122,7 @@ void calc_uiz_R_EV_SH( * @param[out] stats 状态代码,是否有除零错误,非0为异常值 * */ -void calc_RT_PSV( +void grt_calc_RT_PSV( MYREAL Rho1, MYCOMPLEX xa1, MYCOMPLEX xb1, MYCOMPLEX kbkb1, MYCOMPLEX mu1, MYREAL Rho2, MYCOMPLEX xa2, MYCOMPLEX xb2, MYCOMPLEX kbkb2, MYCOMPLEX mu2, MYREAL thk, // 使用上层的厚度 @@ -149,7 +149,7 @@ void calc_RT_PSV( * @param[out] stats 状态代码,是否有除零错误,非0为异常值 * */ -void calc_RT_SH( +void grt_calc_RT_SH( MYCOMPLEX xb1, MYCOMPLEX mu1, MYCOMPLEX xb2, MYCOMPLEX mu2, MYREAL thk, // 使用上层的厚度 @@ -158,7 +158,7 @@ void calc_RT_SH( MYCOMPLEX *TDL, MYCOMPLEX *TUL); /** 液-液 界面,函数参数见 calc_RT_PSV 函数 */ -void calc_RT_ll_PSV( +void grt_calc_RT_ll_PSV( MYREAL Rho1, MYCOMPLEX xa1, MYREAL Rho2, MYCOMPLEX xa2, MYREAL thk, // 使用上层的厚度 @@ -167,12 +167,12 @@ void calc_RT_ll_PSV( MYCOMPLEX TD[2][2], MYCOMPLEX TU[2][2], MYINT *stats); /** 液-液 界面,函数参数见 calc_RT_SH 函数 */ -void calc_RT_ll_SH( +void grt_calc_RT_ll_SH( MYCOMPLEX *RDL, MYCOMPLEX *RUL, MYCOMPLEX *TDL, MYCOMPLEX *TUL); /** 液-固 界面,函数参数见 calc_RT_PSV 函数 */ -void calc_RT_ls_PSV( +void grt_calc_RT_ls_PSV( MYREAL Rho1, MYCOMPLEX xa1, MYCOMPLEX xb1, MYCOMPLEX kbkb1, MYCOMPLEX mu1, MYREAL Rho2, MYCOMPLEX xa2, MYCOMPLEX xb2, MYCOMPLEX kbkb2, MYCOMPLEX mu2, MYREAL thk, // 使用上层的厚度 @@ -181,7 +181,7 @@ void calc_RT_ls_PSV( MYCOMPLEX TD[2][2], MYCOMPLEX TU[2][2], MYINT *stats); /** 液-固 界面,函数参数见 calc_RT_SH 函数 */ -void calc_RT_ls_SH( +void grt_calc_RT_ls_SH( MYCOMPLEX xb1, MYCOMPLEX mu1, MYCOMPLEX mu2, MYREAL thk, // 使用上层的厚度 MYCOMPLEX omega, MYREAL k, @@ -189,7 +189,7 @@ void calc_RT_ls_SH( MYCOMPLEX *TDL, MYCOMPLEX *TUL); /** 固-固 界面,函数参数见 calc_RT_PSV 函数 */ -void calc_RT_ss_PSV( +void grt_calc_RT_ss_PSV( MYREAL Rho1, MYCOMPLEX xa1, MYCOMPLEX xb1, MYCOMPLEX kbkb1, MYCOMPLEX mu1, MYREAL Rho2, MYCOMPLEX xa2, MYCOMPLEX xb2, MYCOMPLEX kbkb2, MYCOMPLEX mu2, MYREAL thk, // 使用上层的厚度 @@ -198,7 +198,7 @@ void calc_RT_ss_PSV( MYCOMPLEX TD[2][2], MYCOMPLEX TU[2][2], MYINT *stats); /** 固-固 界面,函数参数见 calc_RT_SH 函数 */ -void calc_RT_ss_SH( +void grt_calc_RT_ss_SH( MYCOMPLEX xb1, MYCOMPLEX mu1, MYCOMPLEX xb2, MYCOMPLEX mu2, MYREAL thk, // 使用上层的厚度 @@ -221,34 +221,34 @@ void calc_RT_ss_SH( * @param[in] inverse 是否生成逆矩阵 * */ -void get_layer_D( +void grt_get_layer_D( MYCOMPLEX xa, MYCOMPLEX xb, MYCOMPLEX kbkb, MYCOMPLEX mu, MYCOMPLEX omega, MYREAL k, MYCOMPLEX D[4][4], bool inverse); /** 子矩阵 D11,函数参数见 get_layer_D 函数 */ -void get_layer_D11( +void grt_get_layer_D11( MYCOMPLEX xa, MYCOMPLEX xb, MYREAL k, MYCOMPLEX D[2][2]); /** 子矩阵 D12,函数参数见 get_layer_D 函数 */ -void get_layer_D12( +void grt_get_layer_D12( MYCOMPLEX xa, MYCOMPLEX xb, MYREAL k, MYCOMPLEX D[2][2]); /** 子矩阵 D21,函数参数见 get_layer_D 函数 */ -void get_layer_D21( +void grt_get_layer_D21( MYCOMPLEX xa, MYCOMPLEX xb, MYCOMPLEX kbkb, MYCOMPLEX mu, MYCOMPLEX omega, MYREAL k, MYCOMPLEX D[2][2]); /** 子矩阵 D22,函数参数见 get_layer_D 函数 */ -void get_layer_D22( +void grt_get_layer_D22( MYCOMPLEX xa, MYCOMPLEX xb, MYCOMPLEX kbkb, MYCOMPLEX mu, MYCOMPLEX omega, MYREAL k, MYCOMPLEX D[2][2]); /** 子矩阵 D11_uiz,后缀uiz表示连接位移对z的偏导和垂直波函数,函数参数见 get_layer_D 函数 */ -void get_layer_D11_uiz( +void grt_get_layer_D11_uiz( MYCOMPLEX xa, MYCOMPLEX xb, MYREAL k, MYCOMPLEX D[2][2]); /** 子矩阵 D12_uiz,函数参数见 get_layer_D 函数 */ -void get_layer_D12_uiz( +void grt_get_layer_D12_uiz( MYCOMPLEX xa, MYCOMPLEX xb, MYREAL k, MYCOMPLEX D[2][2]); @@ -264,15 +264,15 @@ void get_layer_D12_uiz( * @param[in] inverse 是否生成逆矩阵 * */ -void get_layer_T( +void grt_get_layer_T( MYCOMPLEX xb, MYCOMPLEX mu, MYCOMPLEX omega, MYREAL k, MYCOMPLEX T[2][2], bool inverse); /** 计算 P-SV 型垂直波函数的时间延迟矩阵,公式(5.2.27) */ -void get_layer_E_Rayl(MYCOMPLEX xa1, MYCOMPLEX xb1, MYREAL thk, MYREAL k, MYCOMPLEX E[4][4], bool inverse); +void grt_get_layer_E_Rayl(MYCOMPLEX xa1, MYCOMPLEX xb1, MYREAL thk, MYREAL k, MYCOMPLEX E[4][4], bool inverse); /** 计算 SH 型垂直波函数的时间延迟矩阵,公式(5.2.28) */ -void get_layer_E_Love(MYCOMPLEX xb1, MYREAL thk, MYREAL k, MYCOMPLEX E[2][2], bool inverse); +void grt_get_layer_E_Love(MYCOMPLEX xb1, MYREAL thk, MYREAL k, MYCOMPLEX E[2][2], bool inverse); @@ -281,7 +281,7 @@ void get_layer_E_Love(MYCOMPLEX xb1, MYREAL thk, MYREAL k, MYCOMPLEX E[2][2], bo * 和 calc_RT_PSV(SH) 函数解决相同问题,但没有使用显式推导的公式,而是直接做矩阵运算, * 函数接口也类似 */ -void calc_RT_from_4x4( +void grt_calc_RT_from_4x4( MYCOMPLEX xa1, MYCOMPLEX xb1, MYCOMPLEX kbkb1, MYCOMPLEX mu1, MYCOMPLEX xa2, MYCOMPLEX xb2, MYCOMPLEX kbkb2, MYCOMPLEX mu2, MYCOMPLEX omega, MYREAL thk, diff --git a/pygrt/C_extension/include/dynamic/propagate.h b/pygrt/C_extension/include/grt/dynamic/propagate.h similarity index 95% rename from pygrt/C_extension/include/dynamic/propagate.h rename to pygrt/C_extension/include/grt/dynamic/propagate.h index 6b95d84..c6422df 100755 --- a/pygrt/C_extension/include/dynamic/propagate.h +++ b/pygrt/C_extension/include/grt/dynamic/propagate.h @@ -14,8 +14,8 @@ #pragma once -#include "common/const.h" -#include "common/model.h" +#include "grt/common/const.h" +#include "grt/common/model.h" /** @@ -87,8 +87,8 @@ * @param[out] stats 状态代码,是否有除零错误,非0为异常值 * */ -void kernel( - const MODEL1D *mod1d, MYCOMPLEX omega, MYREAL k, MYCOMPLEX QWV[SRC_M_NUM][QWV_NUM], +void grt_kernel( + const GRT_MODEL1D *mod1d, MYCOMPLEX omega, MYREAL k, MYCOMPLEX QWV[SRC_M_NUM][QWV_NUM], bool calc_uiz, MYCOMPLEX QWV_uiz[SRC_M_NUM][QWV_NUM], MYINT *stats); diff --git a/pygrt/C_extension/include/dynamic/signals.h b/pygrt/C_extension/include/grt/dynamic/signals.h similarity index 82% rename from pygrt/C_extension/include/dynamic/signals.h rename to pygrt/C_extension/include/grt/dynamic/signals.h index bc3e33f..edf4bb7 100644 --- a/pygrt/C_extension/include/dynamic/signals.h +++ b/pygrt/C_extension/include/grt/dynamic/signals.h @@ -26,7 +26,7 @@ * * @return 检查是否通过 */ -bool check_tftype_tfparams(const char tftype, const char *tfparams); +bool grt_check_tftype_tfparams(const char tftype, const char *tfparams); /** * 获得时间函数,要求提前运行check_tftype_tfparams函数以检查参数 @@ -38,7 +38,7 @@ bool check_tftype_tfparams(const char tftype, const char *tfparams); * * @return 时间函数指针 */ -float * get_time_function(int *TFnt, float dt, const char tftype, const char *tfparams); +float * grt_get_time_function(int *TFnt, float dt, const char tftype, const char *tfparams); /** @@ -53,7 +53,7 @@ float * get_time_function(int *TFnt, float dt, const char tftype, const char *tf * @param[out] TFarr 指向时间函数的指针的指针 * @param[out] TFnt 返回的时间函数点数 */ -void linear_convolve_time_function(float *arr, int nt, float dt, const char tftype, const char *tfparams, float **TFarr, int *TFnt); +void grt_linear_convolve_time_function(float *arr, int nt, float dt, const char tftype, const char *tfparams, float **TFarr, int *TFnt); /** @@ -67,7 +67,7 @@ void linear_convolve_time_function(float *arr, int nt, float dt, const char tfty * @param[in] ny 输出数组点数 * @param[in] iscircular 是否使用循环卷积 */ -void oaconvolve(float *x, int nx, float *h, int nh, float *y, int ny, bool iscircular); +void grt_oaconvolve(float *x, int nx, float *h, int nh, float *y, int ny, bool iscircular); /** @@ -79,7 +79,7 @@ void oaconvolve(float *x, int nx, float *h, int nh, float *y, int ny, bool iscir * * @return 积分结果 */ -float trap_area(const float *x, int nx, float dt); +float grt_trap_area(const float *x, int nx, float dt); /** @@ -89,7 +89,7 @@ float trap_area(const float *x, int nx, float dt); * @param[in] nx 数组长度 * @param[in] dt 时间间隔 */ -void trap_integral(float *x, int nx, float dt); +void grt_trap_integral(float *x, int nx, float dt); /** * 对时间序列做中心一阶差分 @@ -98,7 +98,7 @@ void trap_integral(float *x, int nx, float dt); * @param[in] nx 数组长度 * @param[in] dt 时间间隔 */ -void differential(float *x, int nx, float dt); +void grt_differential(float *x, int nx, float dt); @@ -112,7 +112,7 @@ void differential(float *x, int nx, float dt); * * @return float指针 */ -float * get_parabola_wave(float dt, float *Tlen, int *Nt); +float * grt_get_parabola_wave(float dt, float *Tlen, int *Nt); @@ -144,7 +144,7 @@ float * get_parabola_wave(float dt, float *Tlen, int *Nt); * * @return float指针 */ -float * get_trap_wave(float dt, float *T1, float *T2, float *T3, int *Nt); +float * grt_get_trap_wave(float dt, float *T1, float *T2, float *T3, int *Nt); @@ -159,7 +159,7 @@ float * get_trap_wave(float dt, float *T1, float *T2, float *T3, int *Nt); * * @return float指针 */ -float * get_ricker_wave(float dt, float f0, int *Nt); +float * grt_get_ricker_wave(float dt, float f0, int *Nt); /** @@ -170,11 +170,11 @@ float * get_ricker_wave(float dt, float f0, int *Nt); * * @return float指针 */ -float * get_custom_wave(int *Nt, const char *tfparams); +float * grt_get_custom_wave(int *Nt, const char *tfparams); /** * 专用于在Python端释放C中申请的内存 * * @param[out] pt 指针 */ -void free1d(void *pt); \ No newline at end of file +void grt_free1d(void *pt); \ No newline at end of file diff --git a/pygrt/C_extension/include/dynamic/source.h b/pygrt/C_extension/include/grt/dynamic/source.h similarity index 92% rename from pygrt/C_extension/include/dynamic/source.h rename to pygrt/C_extension/include/grt/dynamic/source.h index fc9fbae..6ca736c 100755 --- a/pygrt/C_extension/include/dynamic/source.h +++ b/pygrt/C_extension/include/grt/dynamic/source.h @@ -11,8 +11,8 @@ #pragma once -#include "common/const.h" -#include "common/model.h" +#include "grt/common/const.h" +#include "grt/common/model.h" /** @@ -28,13 +28,13 @@ * @param[out] coef 震源系数 \f$ P_m, SV_m, SH_m \f$ * */ -void source_coef_PSV( +void grt_source_coef_PSV( MYCOMPLEX src_xa, MYCOMPLEX src_xb, MYCOMPLEX src_kaka, MYCOMPLEX src_kbkb, MYREAL k, MYCOMPLEX coef[SRC_M_NUM][QWV_NUM-1][2]); /* SH 波的震源系数,参数见 source_coef_PSV */ -void source_coef_SH( +void grt_source_coef_SH( MYCOMPLEX src_xb, MYCOMPLEX src_kbkb, MYREAL k, MYCOMPLEX coef[SRC_M_NUM][2]); diff --git a/pygrt/C_extension/include/static/static_grn.h b/pygrt/C_extension/include/grt/static/static_grn.h similarity index 92% rename from pygrt/C_extension/include/static/static_grn.h rename to pygrt/C_extension/include/grt/static/static_grn.h index f90c70a..415bf67 100644 --- a/pygrt/C_extension/include/static/static_grn.h +++ b/pygrt/C_extension/include/grt/static/static_grn.h @@ -16,7 +16,7 @@ #pragma once -#include "common/model.h" +#include "grt/common/model.h" /** @@ -41,8 +41,8 @@ * @param[in] statsstr 积分过程输出目录 * */ -void integ_static_grn( - PYMODEL1D *pymod1d, MYINT nr, MYREAL *rs, MYREAL vmin_ref, MYREAL keps, MYREAL k0, MYREAL Length, +void grt_integ_static_grn( + GRT_PYMODEL1D *pymod1d, MYINT nr, MYREAL *rs, MYREAL vmin_ref, MYREAL keps, MYREAL k0, MYREAL Length, MYREAL filonLength, MYREAL safilonTol, MYREAL filonCut, // 返回值,代表Z、R、T分量 diff --git a/pygrt/C_extension/include/static/static_layer.h b/pygrt/C_extension/include/grt/static/static_layer.h similarity index 89% rename from pygrt/C_extension/include/static/static_layer.h rename to pygrt/C_extension/include/grt/static/static_layer.h index ffe7654..2862273 100644 --- a/pygrt/C_extension/include/static/static_layer.h +++ b/pygrt/C_extension/include/grt/static/static_layer.h @@ -15,7 +15,7 @@ #include -#include "common/const.h" +#include "grt/common/const.h" /** * 计算自由表面的 P-SV 静态反射系数,公式(6.3.12) @@ -24,7 +24,7 @@ * @param[out] R_tilt P-SV系数矩阵,SH系数为1 * */ -void calc_static_R_tilt_PSV(MYCOMPLEX delta1, MYCOMPLEX R_tilt[2][2]); +void grt_calc_static_R_tilt_PSV(MYCOMPLEX delta1, MYCOMPLEX R_tilt[2][2]); /** @@ -35,7 +35,7 @@ void calc_static_R_tilt_PSV(MYCOMPLEX delta1, MYCOMPLEX R_tilt[2][2]); * @param[out] R_EV P-SV接收函数矩阵 * */ -void calc_static_R_EV_PSV(bool ircvup, const MYCOMPLEX R[2][2], MYCOMPLEX R_EV[2][2]); +void grt_calc_static_R_EV_PSV(bool ircvup, const MYCOMPLEX R[2][2], MYCOMPLEX R_EV[2][2]); /** * 计算接收点位置的 SH 静态接收矩阵,将波场转为位移,公式(6.3.37) @@ -44,7 +44,7 @@ void calc_static_R_EV_PSV(bool ircvup, const MYCOMPLEX R[2][2], MYCOMPLEX R_EV[2 * @param[out] R_EVL SH接收函数值 * */ -void calc_static_R_EV_SH(MYCOMPLEX RL, MYCOMPLEX *R_EVL); +void grt_calc_static_R_EV_SH(MYCOMPLEX RL, MYCOMPLEX *R_EVL); /** * 计算接收点位置的ui_z的 P-SV 静态接收矩阵,即将波场转为ui_z。 @@ -57,7 +57,7 @@ void calc_static_R_EV_SH(MYCOMPLEX RL, MYCOMPLEX *R_EVL); * @param[out] R_EV P-SV接收函数矩阵 * */ -void calc_static_uiz_R_EV_PSV( +void grt_calc_static_uiz_R_EV_PSV( MYCOMPLEX delta1, bool ircvup, MYREAL k, const MYCOMPLEX R[2][2], MYCOMPLEX R_EV[2][2]); @@ -71,7 +71,7 @@ void calc_static_uiz_R_EV_PSV( * @param[out] R_EVL SH接收函数值 * */ -void calc_static_uiz_R_EV_SH(bool ircvup, MYREAL k, MYCOMPLEX RL, MYCOMPLEX *R_EVL); +void grt_calc_static_uiz_R_EV_SH(bool ircvup, MYREAL k, MYCOMPLEX RL, MYCOMPLEX *R_EVL); /** @@ -90,7 +90,7 @@ void calc_static_uiz_R_EV_SH(bool ircvup, MYREAL k, MYCOMPLEX RL, MYCOMPLEX *R_E * @param[out] TU P-SV 上传透射系数矩阵 * */ -void calc_static_RT_PSV( +void grt_calc_static_RT_PSV( MYCOMPLEX delta1, MYCOMPLEX mu1, MYCOMPLEX delta2, MYCOMPLEX mu2, MYREAL thk, MYREAL k, @@ -110,7 +110,7 @@ void calc_static_RT_PSV( * @param[out] TUL SH 上传透射系数 * */ -void calc_static_RT_SH( +void grt_calc_static_RT_SH( MYCOMPLEX mu1, MYCOMPLEX mu2, MYREAL thk, MYREAL k, MYCOMPLEX *RDL, MYCOMPLEX *RUL, MYCOMPLEX *TDL, MYCOMPLEX *TUL); \ No newline at end of file diff --git a/pygrt/C_extension/include/static/static_propagate.h b/pygrt/C_extension/include/grt/static/static_propagate.h similarity index 82% rename from pygrt/C_extension/include/static/static_propagate.h rename to pygrt/C_extension/include/grt/static/static_propagate.h index ef4e74f..7badf1e 100644 --- a/pygrt/C_extension/include/static/static_propagate.h +++ b/pygrt/C_extension/include/grt/static/static_propagate.h @@ -13,8 +13,8 @@ #pragma once -#include "common/const.h" -#include "common/model.h" +#include "grt/common/const.h" +#include "grt/common/model.h" /** @@ -22,6 +22,6 @@ * 函数参数与动态kernel函数保持一致,具体说明详见`dynamic/propagate.h`。 * */ -void static_kernel( - const MODEL1D *mod1d, MYCOMPLEX omega, MYREAL k, MYCOMPLEX QWV[SRC_M_NUM][QWV_NUM], +void grt_static_kernel( + const GRT_MODEL1D *mod1d, MYCOMPLEX omega, MYREAL k, MYCOMPLEX QWV[SRC_M_NUM][QWV_NUM], bool calc_uiz, MYCOMPLEX QWV_uiz[SRC_M_NUM][QWV_NUM], MYINT *stats); \ No newline at end of file diff --git a/pygrt/C_extension/include/static/static_source.h b/pygrt/C_extension/include/grt/static/static_source.h similarity index 82% rename from pygrt/C_extension/include/static/static_source.h rename to pygrt/C_extension/include/grt/static/static_source.h index f24cdcf..5bf2644 100644 --- a/pygrt/C_extension/include/static/static_source.h +++ b/pygrt/C_extension/include/grt/static/static_source.h @@ -10,7 +10,7 @@ */ #pragma once -#include "common/const.h" +#include "grt/common/const.h" /** * 计算不同震源的静态震源系数,文献/书中仅提供剪切源的震源系数,其它震源系数重新推导 @@ -22,7 +22,7 @@ * @param[in] k 波数 * @param[out] coef 震源系数 \f$ P_m, SV_m, SH_m \f$ */ -void static_source_coef_PSV(MYCOMPLEX delta, MYREAL k, MYCOMPLEX coef[SRC_M_NUM][QWV_NUM-1][2]); +void grt_static_source_coef_PSV(MYCOMPLEX delta, MYREAL k, MYCOMPLEX coef[SRC_M_NUM][QWV_NUM-1][2]); /* SH 波的静态震源系数,参数见 static_source_coef_PSV */ -void static_source_coef_SH(MYREAL k, MYCOMPLEX coef[SRC_M_NUM][2]); \ No newline at end of file +void grt_static_source_coef_SH(MYREAL k, MYCOMPLEX coef[SRC_M_NUM][2]); \ No newline at end of file diff --git a/pygrt/C_extension/include/travt/travt.h b/pygrt/C_extension/include/grt/travt/travt.h similarity index 93% rename from pygrt/C_extension/include/travt/travt.h rename to pygrt/C_extension/include/grt/travt/travt.h index 5db4ced..dd23518 100644 --- a/pygrt/C_extension/include/travt/travt.h +++ b/pygrt/C_extension/include/grt/travt/travt.h @@ -9,7 +9,7 @@ #pragma once -#include "common/const.h" +#include "grt/common/const.h" /** * 已知每层的厚度和速度,且震源和场点位于(虚拟)界面上, @@ -23,6 +23,6 @@ * @param[in] ircv 场点所在层位 * @param[in] dist 震中距 */ -MYREAL compute_travt1d( +MYREAL grt_compute_travt1d( const MYREAL *Thk, const MYREAL *Vel0, const int nlay, const int isrc, const int ircv, const MYREAL dist); \ No newline at end of file diff --git a/pygrt/C_extension/src/common/attenuation.c b/pygrt/C_extension/src/common/attenuation.c index 83507ac..69c6697 100755 --- a/pygrt/C_extension/src/common/attenuation.c +++ b/pygrt/C_extension/src/common/attenuation.c @@ -7,20 +7,20 @@ */ -#include "common/attenuation.h" -#include "common/const.h" +#include "grt/common/attenuation.h" +#include "grt/common/const.h" -MYCOMPLEX attenuation_law(MYREAL Qinv, MYCOMPLEX omega){ +MYCOMPLEX grt_attenuation_law(MYREAL Qinv, MYCOMPLEX omega){ return RONE + Qinv/PI * log(omega/PI2) + RHALF*Qinv*I; // return RONE; } -void py_attenuation_law(MYREAL Qinv, MYREAL omg[2], MYREAL atte[2]){ +void grt_py_attenuation_law(MYREAL Qinv, MYREAL omg[2], MYREAL atte[2]){ // 用于在python中调用attenuation_law MYCOMPLEX omega = omg[0] + I*omg[1]; - MYCOMPLEX atte0 = attenuation_law(Qinv, omega); + MYCOMPLEX atte0 = grt_attenuation_law(Qinv, omega); atte[0] = creal(atte0); atte[1] = cimag(atte0); } \ No newline at end of file diff --git a/pygrt/C_extension/src/common/bessel.c b/pygrt/C_extension/src/common/bessel.c index 3028408..d3f3d5b 100755 --- a/pygrt/C_extension/src/common/bessel.c +++ b/pygrt/C_extension/src/common/bessel.c @@ -7,16 +7,16 @@ */ -#include "common/bessel.h" -#include "common/const.h" +#include "grt/common/bessel.h" +#include "grt/common/const.h" -void bessel012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2){ +void grt_bessel012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2){ *bj0 = j0(x); *bj1 = j1(x); *bj2 = jn(2, x); } -void besselp012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2){ +void grt_besselp012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2){ MYREAL j0=*bj0; MYREAL j1=*bj1; MYREAL j2=*bj2; diff --git a/pygrt/C_extension/src/common/const.c b/pygrt/C_extension/src/common/const.c index ccd3625..1141ac9 100644 --- a/pygrt/C_extension/src/common/const.c +++ b/pygrt/C_extension/src/common/const.c @@ -6,7 +6,7 @@ * 将全局变量放在该文件中 */ -#include "common/const.h" +#include "grt/common/const.h" #include @@ -30,7 +30,7 @@ const char ZNEchs[] = {'Z', 'N', 'E'}; -void set_num_threads(int num_threads){ +void grt_set_num_threads(int num_threads){ #ifdef _OPENMP omp_set_num_threads(num_threads); #endif diff --git a/pygrt/C_extension/src/common/coord.c b/pygrt/C_extension/src/common/coord.c index 66392e3..63ce187 100644 --- a/pygrt/C_extension/src/common/coord.c +++ b/pygrt/C_extension/src/common/coord.c @@ -10,9 +10,9 @@ #include #include +#include "grt/common/coord.h" - -void rot_zxy2zrt_vec(double theta, double A[3]){ +void grt_rot_zxy2zrt_vec(double theta, double A[3]){ double s1, s2, s3; s1 = A[0]; s2 = A[1]; s3 = A[2]; double st = sin(theta); @@ -24,7 +24,7 @@ void rot_zxy2zrt_vec(double theta, double A[3]){ -void rot_zxy2zrt_symtensor2odr(double theta, double A[6]) { +void grt_rot_zxy2zrt_symtensor2odr(double theta, double A[6]) { double s11, s12, s13, s22, s23, s33; s11 = A[0]; s12 = A[1]; s13 = A[2]; s22 = A[3]; s23 = A[4]; @@ -45,7 +45,7 @@ void rot_zxy2zrt_symtensor2odr(double theta, double A[6]) { -void rot_zrt2zxy_upar(const double theta, double u[3], double upar[3][3], const double r){ +void grt_rot_zrt2zxy_upar(const double theta, double u[3], double upar[3][3], const double r){ double s00, s01, s02; double s10, s11, s12; double s20, s21, s22; diff --git a/pygrt/C_extension/src/common/dwm.c b/pygrt/C_extension/src/common/dwm.c index c176fb5..d151d55 100644 --- a/pygrt/C_extension/src/common/dwm.c +++ b/pygrt/C_extension/src/common/dwm.c @@ -15,22 +15,22 @@ #include #include -#include "common/dwm.h" -#include "common/kernel.h" -#include "common/integral.h" -#include "common/iostats.h" -#include "common/model.h" -#include "common/const.h" +#include "grt/common/dwm.h" +#include "grt/common/kernel.h" +#include "grt/common/integral.h" +#include "grt/common/iostats.h" +#include "grt/common/model.h" +#include "grt/common/const.h" -MYREAL discrete_integ( - const MODEL1D *mod1d, MYREAL dk, MYREAL kmax, MYREAL keps, MYCOMPLEX omega, +MYREAL grt_discrete_integ( + const GRT_MODEL1D *mod1d, MYREAL dk, MYREAL kmax, MYREAL keps, MYCOMPLEX omega, MYINT nr, MYREAL *rs, MYCOMPLEX sum_J[nr][SRC_M_NUM][INTEG_NUM], bool calc_upar, MYCOMPLEX sum_uiz_J[nr][SRC_M_NUM][INTEG_NUM], MYCOMPLEX sum_uir_J[nr][SRC_M_NUM][INTEG_NUM], - FILE *fstats, KernelFunc kerfunc, MYINT *stats) + FILE *fstats, GRT_KernelFunc kerfunc, MYINT *stats) { MYCOMPLEX SUM[SRC_M_NUM][INTEG_NUM]; @@ -62,7 +62,7 @@ MYREAL discrete_integ( if(*stats==INVERSE_FAILURE) goto BEFORE_RETURN; // 记录积分核函数 - if(fstats!=NULL) write_stats(fstats, k, QWV); + if(fstats!=NULL) grt_write_stats(fstats, k, QWV); // 震中距rs循环 iendk = true; @@ -76,7 +76,7 @@ MYREAL discrete_integ( } // 计算被积函数一项 F(k,w)Jm(kr)k - int_Pk(k, rs[ir], QWV, false, SUM); + grt_int_Pk(k, rs[ir], QWV, false, SUM); iendk0 = true; for(MYINT i=0; i #include -#include "common/fim.h" -#include "common/integral.h" -#include "common/iostats.h" -#include "common/const.h" -#include "common/model.h" +#include "grt/common/fim.h" +#include "grt/common/integral.h" +#include "grt/common/iostats.h" +#include "grt/common/const.h" +#include "grt/common/model.h" -MYREAL linear_filon_integ( - const MODEL1D *mod1d, MYREAL k0, MYREAL dk0, MYREAL dk, MYREAL kmax, MYREAL keps, MYCOMPLEX omega, +MYREAL grt_linear_filon_integ( + const GRT_MODEL1D *mod1d, MYREAL k0, MYREAL dk0, MYREAL dk, MYREAL kmax, MYREAL keps, MYCOMPLEX omega, MYINT nr, MYREAL *rs, MYCOMPLEX sum_J0[nr][SRC_M_NUM][INTEG_NUM], bool calc_upar, MYCOMPLEX sum_uiz_J0[nr][SRC_M_NUM][INTEG_NUM], MYCOMPLEX sum_uir_J0[nr][SRC_M_NUM][INTEG_NUM], - FILE *fstats, KernelFunc kerfunc, MYINT *stats) + FILE *fstats, GRT_KernelFunc kerfunc, MYINT *stats) { // 从0开始,存储第二部分Filon积分的结果 MYCOMPLEX (*sum_J)[SRC_M_NUM][INTEG_NUM] = (MYCOMPLEX(*)[SRC_M_NUM][INTEG_NUM])calloc(nr, sizeof(*sum_J)); @@ -63,7 +63,7 @@ MYREAL linear_filon_integ( if(*stats==INVERSE_FAILURE) goto BEFORE_RETURN; // 记录积分结果 - if(fstats!=NULL) write_stats(fstats, k, QWV); + if(fstats!=NULL) grt_write_stats(fstats, k, QWV); // 震中距rs循环 iendk = true; @@ -77,7 +77,7 @@ MYREAL linear_filon_integ( } // F(k, w)*Jm(kr)k 的近似公式, sqrt(k) * F(k,w) * cos - int_Pk_filon(k, rs[ir], true, QWV, false, SUM); + grt_int_Pk_filon(k, rs[ir], true, QWV, false, SUM); iendk0 = true; for(MYINT i=0; iname, Ctrl->s_filepath, &hd, NULL); + float *arr = grt_read_SAC(Ctrl->name, Ctrl->s_filepath, &hd, NULL); // 将波形写入标准输出,第一列时间,第二列振幅 float begt = hd.b; diff --git a/pygrt/C_extension/src/common/grt_k2a.c b/pygrt/C_extension/src/common/grt_k2a.c index a78d617..1ab6663 100644 --- a/pygrt/C_extension/src/common/grt_k2a.c +++ b/pygrt/C_extension/src/common/grt_k2a.c @@ -8,9 +8,9 @@ * */ -#include "common/const.h" -#include "common/iostats.h" -#include "common/util.h" +#include "grt/common/const.h" +#include "grt/common/iostats.h" +#include "grt/common/util.h" #include "grt.h" @@ -66,12 +66,12 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ */ static void print_K(FILE *fp){ // 打印标题 - extract_stats(NULL, stdout); + grt_extract_stats(NULL, stdout); fprintf(stdout, "\n"); // 读取数据 while (true) { - if(0 != extract_stats(fp, stdout)) break; + if(0 != grt_extract_stats(fp, stdout)) break; fprintf(stdout, "\n"); } @@ -84,12 +84,12 @@ static void print_K(FILE *fp){ */ static void print_PTAM(FILE *fp){ // 打印标题 - extract_stats_ptam(NULL, stdout); + grt_extract_stats(NULL, stdout); fprintf(stdout, "\n"); // 读取数据 while (true) { - if(0 != extract_stats_ptam(fp, stdout)) break; + if(0 != grt_extract_stats_ptam(fp, stdout)) break; fprintf(stdout, "\n"); } @@ -113,7 +113,7 @@ int k2a_main(int argc, char **argv){ FILE *fp = GRTCheckOpenFile(Ctrl->name, Ctrl->s_filepath, "rb"); // 根据文件名确定函数 - const char *basename = get_basename(Ctrl->s_filepath); + const char *basename = grt_get_basename(Ctrl->s_filepath); if(strncmp(basename, "PTAM", 4) == 0) { print_PTAM(fp); } else if(strncmp(basename, "K", 1) == 0) { diff --git a/pygrt/C_extension/src/common/integral.c b/pygrt/C_extension/src/common/integral.c index 18e2e54..6f869c9 100644 --- a/pygrt/C_extension/src/common/integral.c +++ b/pygrt/C_extension/src/common/integral.c @@ -12,14 +12,14 @@ #include #include -#include "common/integral.h" -#include "common/const.h" -#include "common/bessel.h" -#include "common/quadratic.h" +#include "grt/common/integral.h" +#include "grt/common/const.h" +#include "grt/common/bessel.h" +#include "grt/common/quadratic.h" -void int_Pk( +void grt_int_Pk( MYREAL k, MYREAL r, // F(ki,w), 第一个维度表示不同震源,不同阶数,第二个维度3代表三类系数qm,wm,vm const MYCOMPLEX QWV[SRC_M_NUM][QWV_NUM], @@ -34,12 +34,12 @@ void int_Pk( MYREAL Jmcoef[MORDER_MAX+1] = {0}; - bessel012(kr, &bjmk[0], &bjmk[1], &bjmk[2]); + grt_bessel012(kr, &bjmk[0], &bjmk[1], &bjmk[2]); if(calc_uir){ MYREAL bjmk0[MORDER_MAX+1] = {0}; for(MYINT i=0; i<=MORDER_MAX; ++i) bjmk0[i] = bjmk[i]; - besselp012(kr, &bjmk[0], &bjmk[1], &bjmk[2]); + grt_besselp012(kr, &bjmk[0], &bjmk[1], &bjmk[2]); kcoef = k*k; for(MYINT i=1; i<=MORDER_MAX; ++i) Jmcoef[i] = kr_inv * (-kr_inv * bjmk0[i] + bjmk[i]); @@ -70,7 +70,7 @@ void int_Pk( } -void int_Pk_filon( +void grt_int_Pk_filon( MYREAL k, MYREAL r, bool iscos, const MYCOMPLEX QWV[SRC_M_NUM][QWV_NUM], bool calc_uir, @@ -159,7 +159,7 @@ static MYCOMPLEX interg_quad_cos( -void int_Pk_sa_filon( +void grt_int_Pk_sa_filon( const MYREAL k3[3], MYREAL r, const MYCOMPLEX QWV3[3][SRC_M_NUM][QWV_NUM], bool calc_uir, @@ -185,7 +185,7 @@ void int_Pk_sa_filon( for(MYINT d=0; d<3; ++d) F3[d] = QWV3[d][im][c] * sqrt(k3[d]) * sgn; // 拟合参数 - quad_term(k3, F3, &quad_a[im][c], &quad_b[im][c], &quad_c[im][c]); + grt_quad_term(k3, F3, &quad_a[im][c], &quad_b[im][c], &quad_c[im][c]); } } @@ -208,7 +208,7 @@ void int_Pk_sa_filon( -void merge_Pk( +void grt_merge_Pk( // F(ki,w)Jm(ki*r)ki, const MYCOMPLEX sum_J[SRC_M_NUM][INTEG_NUM], // 累积求和,Z、R、T分量 diff --git a/pygrt/C_extension/src/common/iostats.c b/pygrt/C_extension/src/common/iostats.c index 5705b43..e612ae6 100755 --- a/pygrt/C_extension/src/common/iostats.c +++ b/pygrt/C_extension/src/common/iostats.c @@ -12,12 +12,12 @@ #include #include -#include "common/iostats.h" -#include "common/const.h" +#include "grt/common/iostats.h" +#include "grt/common/const.h" -void write_stats( +void grt_write_stats( FILE *f0, MYREAL k, const MYCOMPLEX QWV[SRC_M_NUM][QWV_NUM]) { fwrite(&k, sizeof(MYREAL), 1, f0); @@ -33,7 +33,7 @@ void write_stats( } -MYINT extract_stats(FILE *bf0, FILE *af0){ +MYINT grt_extract_stats(FILE *bf0, FILE *af0){ // 打印标题 if(bf0 == NULL){ char K[20]; @@ -73,7 +73,7 @@ MYINT extract_stats(FILE *bf0, FILE *af0){ } -void write_stats_ptam( +void grt_write_stats_ptam( FILE *f0, MYREAL Kpt[SRC_M_NUM][INTEG_NUM][PTAM_MAX_PT], MYCOMPLEX Fpt[SRC_M_NUM][INTEG_NUM][PTAM_MAX_PT] @@ -94,7 +94,7 @@ void write_stats_ptam( } -MYINT extract_stats_ptam(FILE *bf0, FILE *af0){ +MYINT grt_extract_stats_ptam(FILE *bf0, FILE *af0){ // 打印标题 if(bf0 == NULL){ char K[20], K2[20]; diff --git a/pygrt/C_extension/src/common/model.c b/pygrt/C_extension/src/common/model.c index 743ad80..49ce9a3 100755 --- a/pygrt/C_extension/src/common/model.c +++ b/pygrt/C_extension/src/common/model.c @@ -12,14 +12,14 @@ #include #include -#include "common/model.h" -#include "common/prtdbg.h" -#include "common/attenuation.h" -#include "common/colorstr.h" +#include "grt/common/model.h" +#include "grt/common/prtdbg.h" +#include "grt/common/attenuation.h" +#include "grt/common/colorstr.h" -#include "grt_error.h" +#include "grt/common/checkerror.h" -void print_mod1d(const MODEL1D *mod1d){ +void grt_print_mod1d(const GRT_MODEL1D *mod1d){ LAYER *lay; for(MYINT u=0; u<50; ++u){printf("---"); } printf("\n"); for(MYINT i=0; in; ++i){ @@ -35,7 +35,7 @@ void print_mod1d(const MODEL1D *mod1d){ } } -void print_pymod(const PYMODEL1D *pymod){ +void grt_print_pymod(const GRT_PYMODEL1D *pymod){ // 模拟表格,打印速度 // 每列字符宽度 // [isrc/ircv] [h(km)] [Vp(km/s)] [Vs(km/s)] [Rho(g/cm^3)] [Qp] [Qs] @@ -102,7 +102,7 @@ void print_pymod(const PYMODEL1D *pymod){ printf("\n"); } -void free_pymod(PYMODEL1D *pymod){ +void grt_free_pymod(GRT_PYMODEL1D *pymod){ GRT_SAFE_FREE_PTR(pymod->Thk); GRT_SAFE_FREE_PTR(pymod->Va); GRT_SAFE_FREE_PTR(pymod->Vb); @@ -113,8 +113,8 @@ void free_pymod(PYMODEL1D *pymod){ } -MODEL1D * init_mod1d(MYINT n){ - MODEL1D *mod1d = (MODEL1D *)malloc(sizeof(MODEL1D)); +GRT_MODEL1D * grt_init_mod1d(MYINT n){ + GRT_MODEL1D *mod1d = (GRT_MODEL1D *)malloc(sizeof(GRT_MODEL1D)); mod1d->n = n; mod1d->lays = (LAYER *)malloc(sizeof(LAYER) * n); return mod1d; @@ -122,7 +122,7 @@ MODEL1D * init_mod1d(MYINT n){ -void get_mod1d(const PYMODEL1D *pymod1d, MODEL1D *mod1d){ +void grt_get_mod1d(const GRT_PYMODEL1D *pymod1d, GRT_MODEL1D *mod1d){ MYINT n = pymod1d->n; mod1d->n = n; mod1d->isrc = pymod1d->isrc; @@ -161,7 +161,7 @@ void get_mod1d(const PYMODEL1D *pymod1d, MODEL1D *mod1d){ } -void copy_mod1d(const MODEL1D *mod1d1, MODEL1D *mod1d2){ +void grt_copy_mod1d(const GRT_MODEL1D *mod1d1, GRT_MODEL1D *mod1d2){ MYINT n = mod1d1->n; mod1d2->n = mod1d1->n; mod1d2->isrc = mod1d1->isrc; @@ -195,14 +195,14 @@ void copy_mod1d(const MODEL1D *mod1d1, MODEL1D *mod1d2){ } -void free_mod1d(MODEL1D *mod1d){ +void grt_free_mod1d(GRT_MODEL1D *mod1d){ GRT_SAFE_FREE_PTR(mod1d->lays); GRT_SAFE_FREE_PTR(mod1d); } -void update_mod1d_omega(MODEL1D *mod1d, MYCOMPLEX omega){ +void grt_update_mod1d_omega(GRT_MODEL1D *mod1d, MYCOMPLEX omega){ MYREAL Va0, Vb0; MYCOMPLEX ka0, kb0; MYCOMPLEX atna, atnb; @@ -212,8 +212,8 @@ void update_mod1d_omega(MODEL1D *mod1d, MYCOMPLEX omega){ Va0 = lay->Va; Vb0 = lay->Vb; - atna = (lay->Qainv > 0.0)? attenuation_law(lay->Qainv, omega) : 1.0; - atnb = (lay->Qbinv > 0.0)? attenuation_law(lay->Qbinv, omega) : 1.0; + atna = (lay->Qainv > 0.0)? grt_attenuation_law(lay->Qainv, omega) : 1.0; + atnb = (lay->Qbinv > 0.0)? grt_attenuation_law(lay->Qbinv, omega) : 1.0; ka0 = omega/(Va0*atna); kb0 = (Vb0>RZERO)? omega/(Vb0*atnb) : CZERO; @@ -231,8 +231,8 @@ void update_mod1d_omega(MODEL1D *mod1d, MYCOMPLEX omega){ } -PYMODEL1D * init_pymod(MYINT n){ - PYMODEL1D *pymod = (PYMODEL1D *)calloc(1, sizeof(PYMODEL1D)); +GRT_PYMODEL1D * grt_init_pymod(MYINT n){ + GRT_PYMODEL1D *pymod = (GRT_PYMODEL1D *)calloc(1, sizeof(GRT_PYMODEL1D)); pymod->n = n; pymod->Thk = (MYREAL*)calloc(n, sizeof(MYREAL)); @@ -245,7 +245,7 @@ PYMODEL1D * init_pymod(MYINT n){ return pymod; } -void realloc_pymod(PYMODEL1D *pymod, MYINT n){ +void grt_realloc_pymod(GRT_PYMODEL1D *pymod, MYINT n){ pymod->n = n; pymod->Thk = (MYREAL*)realloc(pymod->Thk, n*sizeof(MYREAL)); @@ -257,7 +257,7 @@ void realloc_pymod(PYMODEL1D *pymod, MYINT n){ } -PYMODEL1D * read_pymod_from_file(const char *command, const char *modelpath, double depsrc, double deprcv, bool allowLiquid){ +GRT_PYMODEL1D * grt_read_pymod_from_file(const char *command, const char *modelpath, double depsrc, double deprcv, bool allowLiquid){ GRTCheckFileExist(command, modelpath); FILE *fp = GRTCheckOpenFile(command, modelpath, "r"); @@ -281,7 +281,7 @@ PYMODEL1D * read_pymod_from_file(const char *command, const char *modelpath, dou pimg_idx = pmin_idx; // 初始化 - PYMODEL1D *pymod = init_pymod(1); + GRT_PYMODEL1D *pymod = grt_init_pymod(1); const int ncols = 6; // 模型文件有6列,或除去qa qb有四列 const int ncols_noQ = 4; @@ -374,7 +374,7 @@ PYMODEL1D * read_pymod_from_file(const char *command, const char *modelpath, dou for(int k=0; k<2; ++k){ // printf("%d, %d, %lf, %lf, %e ", i, k, depth+h, depimg, depth+h- depimg); if(*pimg_idx < 0 && depth+h >= depimg && depsrc >= 0.0 && deprcv >= 0.0){ - realloc_pymod(pymod, nlay+1); + grt_realloc_pymod(pymod, nlay+1); pymod->Thk[nlay] = depimg - depth; pymod->Va[nlay] = va; pymod->Vb[nlay] = vb; @@ -393,7 +393,7 @@ PYMODEL1D * read_pymod_from_file(const char *command, const char *modelpath, dou } - realloc_pymod(pymod, nlay+1); + grt_realloc_pymod(pymod, nlay+1); pymod->Thk[nlay] = h; pymod->Va[nlay] = va; pymod->Vb[nlay] = vb; @@ -440,7 +440,7 @@ PYMODEL1D * read_pymod_from_file(const char *command, const char *modelpath, dou } -void get_pymod_vmin_vmax(const PYMODEL1D *pymod, double *vmin, double *vmax){ +void grt_get_pymod_vmin_vmax(const GRT_PYMODEL1D *pymod, double *vmin, double *vmax){ *vmin = __DBL_MAX__; *vmax = RZERO; const MYREAL *Va = pymod->Va; diff --git a/pygrt/C_extension/src/common/myfftw.c b/pygrt/C_extension/src/common/myfftw.c index 7b5cd2e..92003f4 100644 --- a/pygrt/C_extension/src/common/myfftw.c +++ b/pygrt/C_extension/src/common/myfftw.c @@ -8,14 +8,13 @@ */ #include -#include "common/myfftw.h" - -#include "grt_error.h" +#include "grt/common/myfftw.h" +#include "grt/common/checkerror.h" #define X(T, s, S) \ -static FFTW##S##_HOLDER * init_fftw##s##_holder(const MYINT nt, const MYREAL dt, const MYINT nf_valid, const MYREAL df)\ +static GRT_FFTW##S##_HOLDER * grt_init_fftw##s##_holder(const MYINT nt, const MYREAL dt, const MYINT nf_valid, const MYREAL df)\ {\ - FFTW##S##_HOLDER *fh = (FFTW##S##_HOLDER*)calloc(1, sizeof(FFTW##S##_HOLDER));\ + GRT_FFTW##S##_HOLDER *fh = (GRT_FFTW##S##_HOLDER*)calloc(1, sizeof(GRT_FFTW##S##_HOLDER));\ if (!fh) {\ GRTRaiseError("Failed to allocate memory in function %s\n.", __func__);\ }\ @@ -35,30 +34,30 @@ static FFTW##S##_HOLDER * init_fftw##s##_holder(const MYINT nt, const MYREAL dt }\ \ \ -void reset_fftw##s##_holder_zero(FFTW##S##_HOLDER *fh)\ +void grt_reset_fftw##s##_holder_zero(GRT_FFTW##S##_HOLDER *fh)\ {\ memset(fh->w_t, 0, sizeof(T)*fh->nt);\ memset(fh->W_f, 0, sizeof(fftw##s##_complex)*fh->nf);\ }\ \ \ -FFTW##S##_HOLDER * create_fftw##s##_holder_C2R_1D(const MYINT nt, const MYREAL dt, const MYINT nf_valid, const MYREAL df)\ +GRT_FFTW##S##_HOLDER * grt_create_fftw##s##_holder_C2R_1D(const MYINT nt, const MYREAL dt, const MYINT nf_valid, const MYREAL df)\ {\ - FFTW##S##_HOLDER * fh = init_fftw##s##_holder(nt, dt, nf_valid, df);\ + GRT_FFTW##S##_HOLDER * fh = grt_init_fftw##s##_holder(nt, dt, nf_valid, df);\ fh->plan = fftw##s##_plan_dft_c2r_1d(nt, fh->W_f, fh->w_t, FFTW_ESTIMATE);\ return fh;\ }\ \ \ -FFTW##S##_HOLDER * create_fftw##s##_holder_R2C_1D(const MYINT nt, const MYREAL dt, const MYINT nf_valid, const MYREAL df)\ +GRT_FFTW##S##_HOLDER * grt_create_fftw##s##_holder_R2C_1D(const MYINT nt, const MYREAL dt, const MYINT nf_valid, const MYREAL df)\ {\ - FFTW##S##_HOLDER * fh = init_fftw##s##_holder(nt, dt, nf_valid, df);\ + GRT_FFTW##S##_HOLDER * fh = grt_init_fftw##s##_holder(nt, dt, nf_valid, df);\ fh->plan = fftw##s##_plan_dft_r2c_1d(nt, fh->w_t, fh->W_f, FFTW_ESTIMATE);\ return fh;\ }\ \ \ -void destroy_fftw##s##_holder(FFTW##S##_HOLDER *fh)\ +void grt_destroy_fftw##s##_holder(GRT_FFTW##S##_HOLDER *fh)\ {\ if (fh) {\ fftw##s##_destroy_plan(fh->plan);\ @@ -69,7 +68,7 @@ void destroy_fftw##s##_holder(FFTW##S##_HOLDER *fh)\ }\ \ \ -void naive_inverse_transform_##T(FFTW##S##_HOLDER *fh)\ +void grt_naive_inverse_transform_##T(GRT_FFTW##S##_HOLDER *fh)\ {\ MYINT nt = fh->nt;\ MYINT nf = fh->nf;\ diff --git a/pygrt/C_extension/src/common/progressbar.c b/pygrt/C_extension/src/common/progressbar.c index 2614f92..2ea65fb 100755 --- a/pygrt/C_extension/src/common/progressbar.c +++ b/pygrt/C_extension/src/common/progressbar.c @@ -9,10 +9,10 @@ #include -#include "common/progressbar.h" +#include "grt/common/progressbar.h" -void printprogressBar(const char *prefix, MYINT percentage) { +void grt_printprogressBar(const char *prefix, MYINT percentage) { printf("\r\033[K"); // 移动到行首并清空行 if(prefix!=NULL) printf("%s", prefix); printf("["); diff --git a/pygrt/C_extension/src/common/ptam.c b/pygrt/C_extension/src/common/ptam.c index dd5b243..e7cb1d5 100755 --- a/pygrt/C_extension/src/common/ptam.c +++ b/pygrt/C_extension/src/common/ptam.c @@ -18,12 +18,12 @@ #include #include -#include "common/ptam.h" -#include "common/quadratic.h" -#include "common/integral.h" -#include "common/iostats.h" -#include "common/const.h" -#include "common/model.h" +#include "grt/common/ptam.h" +#include "grt/common/quadratic.h" +#include "grt/common/integral.h" +#include "grt/common/iostats.h" +#include "grt/common/const.h" +#include "grt/common/model.h" /** @@ -48,7 +48,7 @@ static void process_peak_or_trough( { MYCOMPLEX tmp0; if (Gpt[ir][im][v] >= PTAM_WINDOW_SIZE-1 && Ipt[ir][im][v] < PTAM_MAX_PT) { - if (cplx_peak_or_trough(im, v, J3[ir], k, dk, &Kpt[ir][im][v][Ipt[ir][im][v]], &tmp0) != 0) { + if (grt_cplx_peak_or_trough(im, v, J3[ir], k, dk, &Kpt[ir][im][v][Ipt[ir][im][v]], &tmp0) != 0) { Fpt[ir][im][v][Ipt[ir][im][v]++] = tmp0; Gpt[ir][im][v] = 0; } else if (Gpt[ir][im][v] >= PTAM_MAX_WAITS) { // 不再等待,直接取中点作为波峰波谷 @@ -116,14 +116,14 @@ static void ptam_once( } -void PTA_method( - const MODEL1D *mod1d, MYREAL k0, MYREAL predk, MYCOMPLEX omega, +void grt_PTA_method( + const GRT_MODEL1D *mod1d, MYREAL k0, MYREAL predk, MYCOMPLEX omega, MYINT nr, MYREAL *rs, MYCOMPLEX sum_J0[nr][SRC_M_NUM][INTEG_NUM], bool calc_upar, MYCOMPLEX sum_uiz_J0[nr][SRC_M_NUM][INTEG_NUM], MYCOMPLEX sum_uir_J0[nr][SRC_M_NUM][INTEG_NUM], - FILE *ptam_fstatsnr[nr][2], KernelFunc kerfunc, MYINT *stats) + FILE *ptam_fstatsnr[nr][2], GRT_KernelFunc kerfunc, MYINT *stats) { // 需要兼容对正常收敛而不具有规律波峰波谷的序列 // 有时序列收敛比较好,不表现为规律的波峰波谷, @@ -205,10 +205,10 @@ void PTA_method( if(*stats==INVERSE_FAILURE) goto BEFORE_RETURN; // 记录核函数 - if(fstatsK!=NULL) write_stats(fstatsK, k, QWV); + if(fstatsK!=NULL) grt_write_stats(fstatsK, k, QWV); // 计算被积函数一项 F(k,w)Jm(kr)k - int_Pk(k, rs[ir], QWV, false, SUM3[ir][PTAM_WINDOW_SIZE-1]); // [PTAM_WINDOW_SIZE-1]表示把新点值放在最后 + grt_int_Pk(k, rs[ir], QWV, false, SUM3[ir][PTAM_WINDOW_SIZE-1]); // [PTAM_WINDOW_SIZE-1]表示把新点值放在最后 // 判断和记录波峰波谷 ptam_once( ir, nr, precoef, k, dk, SUM3, sum_J, Kpt, Fpt, Ipt, Gpt, &iendk0); @@ -216,13 +216,13 @@ void PTA_method( if(calc_upar){ // ------------------------------- ui_z ----------------------------------- // 计算被积函数一项 F(k,w)Jm(kr)k - int_Pk(k, rs[ir], QWV_uiz, false, SUM3_uiz[ir][PTAM_WINDOW_SIZE-1]); // [PTAM_WINDOW_SIZE-1]表示把新点值放在最后 + grt_int_Pk(k, rs[ir], QWV_uiz, false, SUM3_uiz[ir][PTAM_WINDOW_SIZE-1]); // [PTAM_WINDOW_SIZE-1]表示把新点值放在最后 // 判断和记录波峰波谷 ptam_once(ir, nr, precoef, k, dk, SUM3_uiz, sum_uiz_J, Kpt_uiz, Fpt_uiz, Ipt_uiz, Gpt_uiz, &iendk0); // ------------------------------- ui_r ----------------------------------- // 计算被积函数一项 F(k,w)Jm(kr)k - int_Pk(k, rs[ir], QWV, true, SUM3_uir[ir][PTAM_WINDOW_SIZE-1]); // [PTAM_WINDOW_SIZE-1]表示把新点值放在最后 + grt_int_Pk(k, rs[ir], QWV, true, SUM3_uir[ir][PTAM_WINDOW_SIZE-1]); // [PTAM_WINDOW_SIZE-1]表示把新点值放在最后 // 判断和记录波峰波谷 ptam_once(ir, nr, precoef, k, dk, SUM3_uir, sum_uir_J, Kpt_uir, Fpt_uir, Ipt_uir, Gpt_uir, &iendk0); @@ -240,18 +240,18 @@ void PTA_method( for(MYINT ir=0; ir1; --n){ for(MYINT i=0; i -#include "common/quadratic.h" -#include "common/const.h" +#include "grt/common/quadratic.h" +#include "grt/common/const.h" -void quad_term(const MYREAL x[3], const MYCOMPLEX f[3], MYCOMPLEX *pa, MYCOMPLEX *pb, MYCOMPLEX *pc) +void grt_quad_term(const MYREAL x[3], const MYCOMPLEX f[3], MYCOMPLEX *pa, MYCOMPLEX *pb, MYCOMPLEX *pc) { MYREAL x1, x2, x3, w1, w2, w3; x1 = x[0]; @@ -42,13 +42,13 @@ void quad_term(const MYREAL x[3], const MYCOMPLEX f[3], MYCOMPLEX *pa, MYCOMPLEX -MYCOMPLEX quad_eval(MYREAL x, MYCOMPLEX a, MYCOMPLEX b, MYCOMPLEX c) +MYCOMPLEX grt_quad_eval(MYREAL x, MYCOMPLEX a, MYCOMPLEX b, MYCOMPLEX c) { return a*x*x + b*x + c; } -MYCOMPLEX quad_integral(MYREAL x1, MYREAL x2, MYCOMPLEX a, MYCOMPLEX b, MYCOMPLEX c) +MYCOMPLEX grt_quad_integral(MYREAL x1, MYREAL x2, MYCOMPLEX a, MYCOMPLEX b, MYCOMPLEX c) { MYREAL xx1, xx2, xxx1, xxx2; xx1 = x1*x1; xx2 = x2*x2; diff --git a/pygrt/C_extension/src/common/radiation.c b/pygrt/C_extension/src/common/radiation.c index 3e970d0..1a06d40 100644 --- a/pygrt/C_extension/src/common/radiation.c +++ b/pygrt/C_extension/src/common/radiation.c @@ -9,10 +9,10 @@ #include -#include "common/radiation.h" -#include "common/const.h" +#include "grt/common/radiation.h" +#include "grt/common/const.h" -void set_source_radiation( +void grt_set_source_radiation( double srcRadi[SRC_M_NUM][CHANNEL_NUM], const int computeType, const bool par_theta, const double M0, const double coef, const double azrad, const double mchn[MECHANISM_NUM] ){ diff --git a/pygrt/C_extension/src/common/recursion.c b/pygrt/C_extension/src/common/recursion.c index 676ca5b..3d228fb 100644 --- a/pygrt/C_extension/src/common/recursion.c +++ b/pygrt/C_extension/src/common/recursion.c @@ -14,28 +14,28 @@ #include #include -#include "common/recursion.h" -#include "common/const.h" -#include "common/matrix.h" +#include "grt/common/recursion.h" +#include "grt/common/const.h" +#include "grt/common/matrix.h" -void recursion_RD( +void grt_recursion_RD( const MYCOMPLEX RD1[2][2], MYCOMPLEX RDL1, const MYCOMPLEX RU1[2][2], MYCOMPLEX RUL1, const MYCOMPLEX TD1[2][2], MYCOMPLEX TDL1, const MYCOMPLEX TU1[2][2], MYCOMPLEX TUL1, const MYCOMPLEX RD2[2][2], MYCOMPLEX RDL2, MYCOMPLEX RD[2][2], MYCOMPLEX *RDL, MYCOMPLEX inv_2x2T[2][2], MYCOMPLEX *invT, MYINT *stats) { - recursion_RD_PSV( + grt_recursion_RD_PSV( RD1, RU1, TD1, TU1, RD2, RD, inv_2x2T, stats); - recursion_RD_SH( + grt_recursion_RD_SH( RDL1, RUL1, TDL1, TUL1, RDL2, RDL, invT, stats); } -void recursion_RD_PSV( +void grt_recursion_RD_PSV( const MYCOMPLEX RD1[2][2], const MYCOMPLEX RU1[2][2], const MYCOMPLEX TD1[2][2], const MYCOMPLEX TU1[2][2], const MYCOMPLEX RD2[2][2], @@ -44,18 +44,18 @@ void recursion_RD_PSV( MYCOMPLEX tmp1[2][2], tmp2[2][2]; // RD, RDL - cmat2x2_mul(RU1, RD2, tmp1); - cmat2x2_one_sub(tmp1); - cmat2x2_inv(tmp1, tmp1, stats); if(*stats==INVERSE_FAILURE) return; - cmat2x2_mul(tmp1, TD1, tmp2); - if(inv_2x2T!=NULL) cmat2x2_assign(tmp2, inv_2x2T); - - cmat2x2_mul(RD2, tmp2, tmp1); - cmat2x2_mul(TU1, tmp1, tmp2); - cmat2x2_add(RD1, tmp2, RD); + grt_cmat2x2_mul(RU1, RD2, tmp1); + grt_cmat2x2_one_sub(tmp1); + grt_cmat2x2_inv(tmp1, tmp1, stats); if(*stats==INVERSE_FAILURE) return; + grt_cmat2x2_mul(tmp1, TD1, tmp2); + if(inv_2x2T!=NULL) grt_cmat2x2_assign(tmp2, inv_2x2T); + + grt_cmat2x2_mul(RD2, tmp2, tmp1); + grt_cmat2x2_mul(TU1, tmp1, tmp2); + grt_cmat2x2_add(RD1, tmp2, RD); } -void recursion_RD_SH( +void grt_recursion_RD_SH( MYCOMPLEX RDL1, MYCOMPLEX RUL1, MYCOMPLEX TDL1, MYCOMPLEX TUL1, MYCOMPLEX RDL2, MYCOMPLEX *RDL, MYCOMPLEX *invT, MYINT *stats) @@ -73,22 +73,22 @@ void recursion_RD_SH( } -void recursion_TD( +void grt_recursion_TD( const MYCOMPLEX RU1[2][2], MYCOMPLEX RUL1, const MYCOMPLEX TD1[2][2], MYCOMPLEX TDL1, const MYCOMPLEX RD2[2][2], MYCOMPLEX RDL2, const MYCOMPLEX TD2[2][2], MYCOMPLEX TDL2, MYCOMPLEX TD[2][2], MYCOMPLEX *TDL, MYCOMPLEX inv_2x2T[2][2], MYCOMPLEX *invT, MYINT *stats) { - recursion_TD_PSV( + grt_recursion_TD_PSV( RU1, TD1, RD2, TD2, TD, inv_2x2T, stats); - recursion_TD_SH( + grt_recursion_TD_SH( RUL1, TDL1, RDL2, TDL2, TDL, invT, stats); } -void recursion_TD_PSV( +void grt_recursion_TD_PSV( const MYCOMPLEX RU1[2][2], const MYCOMPLEX TD1[2][2], const MYCOMPLEX RD2[2][2], const MYCOMPLEX TD2[2][2], MYCOMPLEX TD[2][2], MYCOMPLEX inv_2x2T[2][2], MYINT *stats) @@ -96,15 +96,15 @@ void recursion_TD_PSV( MYCOMPLEX tmp1[2][2], tmp2[2][2]; // TD, TDL - cmat2x2_mul(RU1, RD2, tmp2); - cmat2x2_one_sub(tmp2); - cmat2x2_inv(tmp2, tmp1, stats); if(*stats==INVERSE_FAILURE) return; - cmat2x2_mul(tmp1, TD1, tmp2); - if(inv_2x2T!=NULL) cmat2x2_assign(tmp2, inv_2x2T); - cmat2x2_mul(TD2, tmp2, TD); + grt_cmat2x2_mul(RU1, RD2, tmp2); + grt_cmat2x2_one_sub(tmp2); + grt_cmat2x2_inv(tmp2, tmp1, stats); if(*stats==INVERSE_FAILURE) return; + grt_cmat2x2_mul(tmp1, TD1, tmp2); + if(inv_2x2T!=NULL) grt_cmat2x2_assign(tmp2, inv_2x2T); + grt_cmat2x2_mul(TD2, tmp2, TD); } -void recursion_TD_SH( +void grt_recursion_TD_SH( MYCOMPLEX RUL1, MYCOMPLEX TDL1, MYCOMPLEX RDL2, MYCOMPLEX TDL2, MYCOMPLEX *TDL, MYCOMPLEX *invT, MYINT *stats) @@ -121,23 +121,23 @@ void recursion_TD_SH( if(invT!=NULL) *invT = inv1; } -void recursion_RU( +void grt_recursion_RU( const MYCOMPLEX RU1[2][2], MYCOMPLEX RUL1, const MYCOMPLEX RD2[2][2], MYCOMPLEX RDL2, const MYCOMPLEX RU2[2][2], MYCOMPLEX RUL2, const MYCOMPLEX TD2[2][2], MYCOMPLEX TDL2, const MYCOMPLEX TU2[2][2], MYCOMPLEX TUL2, MYCOMPLEX RU[2][2], MYCOMPLEX *RUL, MYCOMPLEX inv_2x2T[2][2], MYCOMPLEX *invT, MYINT *stats) { - recursion_RU_PSV( + grt_recursion_RU_PSV( RU1, RD2, RU2, TD2, TU2, RU, inv_2x2T, stats); - recursion_RU_SH( + grt_recursion_RU_SH( RUL1, RDL2, RUL2, TDL2, TUL2, RUL, invT, stats); } -void recursion_RU_PSV( +void grt_recursion_RU_PSV( const MYCOMPLEX RU1[2][2], const MYCOMPLEX RD2[2][2], const MYCOMPLEX RU2[2][2], const MYCOMPLEX TD2[2][2], const MYCOMPLEX TU2[2][2], @@ -146,19 +146,19 @@ void recursion_RU_PSV( MYCOMPLEX tmp1[2][2], tmp2[2][2]; // RU, RUL - cmat2x2_mul(RD2, RU1, tmp2); - cmat2x2_one_sub(tmp2); - cmat2x2_inv(tmp2, tmp1, stats); if(*stats==INVERSE_FAILURE) return; - cmat2x2_mul(tmp1, TU2, tmp2); - if(inv_2x2T!=NULL) cmat2x2_assign(tmp2, inv_2x2T); - - cmat2x2_mul(RU1, tmp2, tmp1); - cmat2x2_mul(TD2, tmp1, tmp2); - cmat2x2_add(RU2, tmp2, RU); + grt_cmat2x2_mul(RD2, RU1, tmp2); + grt_cmat2x2_one_sub(tmp2); + grt_cmat2x2_inv(tmp2, tmp1, stats); if(*stats==INVERSE_FAILURE) return; + grt_cmat2x2_mul(tmp1, TU2, tmp2); + if(inv_2x2T!=NULL) grt_cmat2x2_assign(tmp2, inv_2x2T); + + grt_cmat2x2_mul(RU1, tmp2, tmp1); + grt_cmat2x2_mul(TD2, tmp1, tmp2); + grt_cmat2x2_add(RU2, tmp2, RU); } -void recursion_RU_SH( +void grt_recursion_RU_SH( MYCOMPLEX RUL1, MYCOMPLEX RDL2, MYCOMPLEX RUL2, MYCOMPLEX TDL2, MYCOMPLEX TUL2, @@ -177,23 +177,23 @@ void recursion_RU_SH( } -void recursion_TU( +void grt_recursion_TU( const MYCOMPLEX RU1[2][2], MYCOMPLEX RUL1, const MYCOMPLEX TU1[2][2], MYCOMPLEX TUL1, const MYCOMPLEX RD2[2][2], MYCOMPLEX RDL2, const MYCOMPLEX TU2[2][2], MYCOMPLEX TUL2, MYCOMPLEX TU[2][2], MYCOMPLEX *TUL, MYCOMPLEX inv_2x2T[2][2], MYCOMPLEX *invT, MYINT *stats) { - recursion_TU_PSV( + grt_recursion_TU_PSV( RU1, TU1, RD2, TU2, TU, inv_2x2T, stats); - recursion_TU_SH( + grt_recursion_TU_SH( RUL1, TUL1, RDL2, TUL2, TUL, invT, stats); } -void recursion_TU_PSV( +void grt_recursion_TU_PSV( const MYCOMPLEX RU1[2][2], const MYCOMPLEX TU1[2][2], const MYCOMPLEX RD2[2][2], const MYCOMPLEX TU2[2][2], MYCOMPLEX TU[2][2], MYCOMPLEX inv_2x2T[2][2], MYINT *stats) @@ -201,17 +201,17 @@ void recursion_TU_PSV( MYCOMPLEX tmp1[2][2], tmp2[2][2]; // TU, TUL - cmat2x2_mul(RD2, RU1, tmp2); - cmat2x2_one_sub(tmp2); - cmat2x2_inv(tmp2, tmp1, stats); if(*stats==INVERSE_FAILURE) return; - cmat2x2_mul(tmp1, TU2, tmp2); - if(inv_2x2T!=NULL) cmat2x2_assign(tmp2, inv_2x2T); - cmat2x2_mul(TU1, tmp2, TU); + grt_cmat2x2_mul(RD2, RU1, tmp2); + grt_cmat2x2_one_sub(tmp2); + grt_cmat2x2_inv(tmp2, tmp1, stats); if(*stats==INVERSE_FAILURE) return; + grt_cmat2x2_mul(tmp1, TU2, tmp2); + if(inv_2x2T!=NULL) grt_cmat2x2_assign(tmp2, inv_2x2T); + grt_cmat2x2_mul(TU1, tmp2, TU); } -void recursion_TU_SH( +void grt_recursion_TU_SH( MYCOMPLEX RUL1, MYCOMPLEX TUL1, MYCOMPLEX RDL2, MYCOMPLEX TUL2, MYCOMPLEX *TUL, MYCOMPLEX *invT, MYINT *stats) @@ -229,7 +229,7 @@ void recursion_TU_SH( } -void recursion_RT( +void grt_recursion_RT( const MYCOMPLEX RD1[2][2], MYCOMPLEX RDL1, const MYCOMPLEX RU1[2][2], MYCOMPLEX RUL1, const MYCOMPLEX TD1[2][2], MYCOMPLEX TDL1, const MYCOMPLEX TU1[2][2], MYCOMPLEX TUL1, const MYCOMPLEX RD2[2][2], MYCOMPLEX RDL2, const MYCOMPLEX RU2[2][2], MYCOMPLEX RUL2, @@ -237,18 +237,18 @@ void recursion_RT( MYCOMPLEX RD[2][2], MYCOMPLEX *RDL, MYCOMPLEX RU[2][2], MYCOMPLEX *RUL, MYCOMPLEX TD[2][2], MYCOMPLEX *TDL, MYCOMPLEX TU[2][2], MYCOMPLEX *TUL, MYINT *stats) { - recursion_RT_PSV( + grt_recursion_RT_PSV( RD1, RU1, TD1, TU1, RD2, RU2, TD2, TU2, RD, RU, TD, TU, stats); - recursion_RT_SH( + grt_recursion_RT_SH( RDL1, RUL1, TDL1, TUL1, RDL2, RUL2, TDL2, TUL2, RDL, RUL, TDL, TUL, stats); } -void recursion_RT_PSV( +void grt_recursion_RT_PSV( const MYCOMPLEX RD1[2][2], const MYCOMPLEX RU1[2][2], const MYCOMPLEX TD1[2][2], const MYCOMPLEX TU1[2][2], const MYCOMPLEX RD2[2][2], const MYCOMPLEX RU2[2][2], @@ -259,35 +259,35 @@ void recursion_RT_PSV( // 临时矩阵 MYCOMPLEX tmp1[2][2], tmp2[2][2]; - cmat2x2_mul(RU1, RD2, tmp1); - cmat2x2_one_sub(tmp1); - cmat2x2_inv(tmp1, tmp1, stats); if(*stats==INVERSE_FAILURE) return; - cmat2x2_mul(tmp1, TD1, tmp2); + grt_cmat2x2_mul(RU1, RD2, tmp1); + grt_cmat2x2_one_sub(tmp1); + grt_cmat2x2_inv(tmp1, tmp1, stats); if(*stats==INVERSE_FAILURE) return; + grt_cmat2x2_mul(tmp1, TD1, tmp2); // TD - cmat2x2_mul(TD2, tmp2, TD); // 相同的逆阵,节省计算量 + grt_cmat2x2_mul(TD2, tmp2, TD); // 相同的逆阵,节省计算量 // RD - cmat2x2_mul(RD2, tmp2, tmp1); - cmat2x2_mul(TU1, tmp1, tmp2); - cmat2x2_add(RD1, tmp2, RD); + grt_cmat2x2_mul(RD2, tmp2, tmp1); + grt_cmat2x2_mul(TU1, tmp1, tmp2); + grt_cmat2x2_add(RD1, tmp2, RD); - cmat2x2_mul(RD2, RU1, tmp1); - cmat2x2_one_sub(tmp1); - cmat2x2_inv(tmp1, tmp1, stats); if(*stats==INVERSE_FAILURE) return; - cmat2x2_mul(tmp1, TU2, tmp2); + grt_cmat2x2_mul(RD2, RU1, tmp1); + grt_cmat2x2_one_sub(tmp1); + grt_cmat2x2_inv(tmp1, tmp1, stats); if(*stats==INVERSE_FAILURE) return; + grt_cmat2x2_mul(tmp1, TU2, tmp2); // TU - cmat2x2_mul(TU1, tmp2, TU); + grt_cmat2x2_mul(TU1, tmp2, TU); // RU - cmat2x2_mul(RU1, tmp2, tmp1); - cmat2x2_mul(TD2, tmp1, tmp2); - cmat2x2_add(RU2, tmp2, RU); + grt_cmat2x2_mul(RU1, tmp2, tmp1); + grt_cmat2x2_mul(TD2, tmp1, tmp2); + grt_cmat2x2_add(RU2, tmp2, RU); } -void recursion_RT_SH( +void grt_recursion_RT_SH( MYCOMPLEX RDL1, MYCOMPLEX RUL1, MYCOMPLEX TDL1, MYCOMPLEX TUL1, MYCOMPLEX RDL2, MYCOMPLEX RUL2, @@ -320,21 +320,21 @@ void recursion_RT_SH( } -void recursion_RT_imaginary( +void grt_recursion_RT_imaginary( MYCOMPLEX xa1, MYCOMPLEX xb1, MYREAL thk, MYREAL k, // 使用上层的厚度 MYCOMPLEX RU[2][2], MYCOMPLEX *RUL, MYCOMPLEX TD[2][2], MYCOMPLEX *TDL, MYCOMPLEX TU[2][2], MYCOMPLEX *TUL) { - recursion_RT_PSV_imaginary( + grt_recursion_RT_PSV_imaginary( xa1, xb1, thk, k, RU, TD, TU); - recursion_RT_SH_imaginary( + grt_recursion_RT_SH_imaginary( xb1, thk, k, RUL, TDL, TUL); } -void recursion_RT_PSV_imaginary( +void grt_recursion_RT_PSV_imaginary( MYCOMPLEX xa1, MYCOMPLEX xb1, MYREAL thk, MYREAL k, // 使用上层的厚度 MYCOMPLEX RU[2][2], MYCOMPLEX TD[2][2], MYCOMPLEX TU[2][2]) { @@ -359,7 +359,7 @@ void recursion_RT_PSV_imaginary( } -void recursion_RT_SH_imaginary( +void grt_recursion_RT_SH_imaginary( MYCOMPLEX xb1, MYREAL thk, MYREAL k, // 使用上层的厚度 MYCOMPLEX *RUL, MYCOMPLEX *TDL, MYCOMPLEX *TUL) { @@ -376,7 +376,7 @@ void recursion_RT_SH_imaginary( -void get_qwv( +void grt_get_qwv( bool ircvup, const MYCOMPLEX R1[2][2], MYCOMPLEX RL1, const MYCOMPLEX R2[2][2], MYCOMPLEX RL2, @@ -387,15 +387,15 @@ void get_qwv( MYCOMPLEX coefD[2] = {coef_PSV[0][0], coef_PSV[1][0]}; MYCOMPLEX coefU[2] = {coef_PSV[0][1], coef_PSV[1][1]}; if(ircvup){ - cmat2x1_mul(R2, coefD, qw0); + grt_cmat2x1_mul(R2, coefD, qw0); qw0[0] += coefU[0]; qw0[1] += coefU[1]; v0 = RL1 * (RL2*coef_SH[0] + coef_SH[1]); } else { - cmat2x1_mul(R2, coefU, qw0); + grt_cmat2x1_mul(R2, coefU, qw0); qw0[0] += coefD[0]; qw0[1] += coefD[1]; v0 = RL1 * (coef_SH[0] + RL2*coef_SH[1]); } - cmat2x1_mul(R1, qw0, qw1); + grt_cmat2x1_mul(R1, qw0, qw1); qwv[0] = qw1[0]; qwv[1] = qw1[1]; diff --git a/pygrt/C_extension/src/common/sacio.c b/pygrt/C_extension/src/common/sacio.c index 1210c5e..e6dbd1b 100644 --- a/pygrt/C_extension/src/common/sacio.c +++ b/pygrt/C_extension/src/common/sacio.c @@ -29,7 +29,7 @@ #include #include #include -#include "common/sacio.h" +#include "grt/common/sacio.h" /* function prototype for local use */ static void byte_swap (char *pt, size_t n); diff --git a/pygrt/C_extension/src/common/sacio2.c b/pygrt/C_extension/src/common/sacio2.c index 9b5e64e..d371aed 100644 --- a/pygrt/C_extension/src/common/sacio2.c +++ b/pygrt/C_extension/src/common/sacio2.c @@ -10,12 +10,12 @@ #include #include -#include "common/sacio2.h" -#include "common/sacio.h" -#include "common/colorstr.h" +#include "grt/common/sacio2.h" +#include "grt/common/sacio.h" +#include "grt/common/colorstr.h" -void read_SAC_HEAD(const char *command, const char *name, SACHEAD *hd){ +void grt_read_SAC_HEAD(const char *command, const char *name, SACHEAD *hd){ int lswap = read_sac_head(name, hd); if(lswap == -1){ fprintf(stderr, "[%s] " BOLD_RED "read %s head failed.\n" DEFAULT_RESTORE, command, name); @@ -24,7 +24,7 @@ void read_SAC_HEAD(const char *command, const char *name, SACHEAD *hd){ } -float * read_SAC(const char *command, const char *name, SACHEAD *hd, float *arrout){ +float * grt_read_SAC(const char *command, const char *name, SACHEAD *hd, float *arrout){ float *arrin=NULL; if((arrin = read_sac(name, hd)) == NULL){ fprintf(stderr, "[%s] " BOLD_RED "read %s failed.\n" DEFAULT_RESTORE, command, name); diff --git a/pygrt/C_extension/src/common/safim.c b/pygrt/C_extension/src/common/safim.c index a1d3d70..a2594aa 100644 --- a/pygrt/C_extension/src/common/safim.c +++ b/pygrt/C_extension/src/common/safim.c @@ -17,11 +17,11 @@ #include #include -#include "common/safim.h" -#include "common/integral.h" -#include "common/iostats.h" -#include "common/const.h" -#include "common/model.h" +#include "grt/common/safim.h" +#include "grt/common/integral.h" +#include "grt/common/iostats.h" +#include "grt/common/const.h" +#include "grt/common/model.h" /** @@ -255,7 +255,7 @@ static void interv_integ( } // 该分段内的积分 - int_Pk_sa_filon(ptKitv->k3, rs[ir], ptKitv->F3, false, SUM); + grt_int_Pk_sa_filon(ptKitv->k3, rs[ir], ptKitv->F3, false, SUM); for(MYINT i=0; ik3, rs[ir], ptKitv->F3_uiz, false, SUM); + grt_int_Pk_sa_filon(ptKitv->k3, rs[ir], ptKitv->F3_uiz, false, SUM); for(MYINT i=0; ik3, rs[ir], ptKitv->F3, true, SUM); + grt_int_Pk_sa_filon(ptKitv->k3, rs[ir], ptKitv->F3, true, SUM); for(MYINT i=0; i 0) { @@ -407,10 +407,10 @@ MYREAL sa_filon_integ( // 记录后四个采样值 if(fstats!=NULL){ for(MYINT i=1; i<3; ++i){ - write_stats(fstats, Kitv_left.k3[i], Kitv_left.F3[i]); + grt_write_stats(fstats, Kitv_left.k3[i], Kitv_left.F3[i]); } for(MYINT i=1; i<3; ++i){ - write_stats(fstats, Kitv_right.k3[i], Kitv_right.F3[i]); + grt_write_stats(fstats, Kitv_right.k3[i], Kitv_right.F3[i]); } } // 计算积分 diff --git a/pygrt/C_extension/src/common/search.c b/pygrt/C_extension/src/common/search.c index 72bb7f8..21b6a7d 100644 --- a/pygrt/C_extension/src/common/search.c +++ b/pygrt/C_extension/src/common/search.c @@ -9,14 +9,14 @@ #include #include -#include "common/search.h" -#include "common/const.h" +#include "grt/common/search.h" +#include "grt/common/const.h" static bool _gt_(MYREAL a1, MYREAL a2) { return a1 > a2; } static bool _lt_(MYREAL a1, MYREAL a2) { return a1 < a2; } -MYINT findElement_MYINT(const MYINT array[], MYINT size, MYINT target) { +MYINT grt_findElement_MYINT(const MYINT array[], MYINT size, MYINT target) { for (MYINT i = 0; i < size; ++i) { if (array[i] == target) { return i; // 找到目标元素,返回索引 @@ -25,7 +25,7 @@ MYINT findElement_MYINT(const MYINT array[], MYINT size, MYINT target) { return -1; // 未找到目标元素,返回-1 } -MYINT findLessEqualClosest_MYREAL(const MYREAL array[], MYINT size, MYREAL target) { +MYINT grt_findLessEqualClosest_MYREAL(const MYREAL array[], MYINT size, MYREAL target) { MYINT ires=-1; MYREAL mindist=-1.0, dist=0.0; for (MYINT i = 0; i < size; ++i) { @@ -38,7 +38,7 @@ MYINT findLessEqualClosest_MYREAL(const MYREAL array[], MYINT size, MYREAL targe return ires; } -MYINT findClosest_MYREAL(const MYREAL array[], MYINT size, MYREAL target) { +MYINT grt_findClosest_MYREAL(const MYREAL array[], MYINT size, MYREAL target) { MYINT ires=0; MYREAL mindist=-1.0, dist=0.0; for (MYINT i = 0; i < size; ++i) { @@ -51,7 +51,7 @@ MYINT findClosest_MYREAL(const MYREAL array[], MYINT size, MYREAL target) { return ires; } -MYINT findMinMax_MYREAL(const MYREAL array[], MYINT size, bool isMax) { +MYINT grt_findMinMax_MYREAL(const MYREAL array[], MYINT size, bool isMax) { MYREAL rmax = array[0]; MYINT idx=0; bool (*_func)(MYREAL, MYREAL) = (isMax)? _gt_ : _lt_; diff --git a/pygrt/C_extension/src/common/util.c b/pygrt/C_extension/src/common/util.c index 1234cf7..7fa8ade 100644 --- a/pygrt/C_extension/src/common/util.c +++ b/pygrt/C_extension/src/common/util.c @@ -11,14 +11,14 @@ #include #include -#include "common/util.h" -#include "common/model.h" -#include "common/const.h" -#include "common/sacio.h" -#include "common/myfftw.h" -#include "travt/travt.h" - -char ** string_split(const char *string, const char *delim, int *size) +#include "grt/common/util.h" +#include "grt/common/model.h" +#include "grt/common/const.h" +#include "grt/common/sacio.h" +#include "grt/common/myfftw.h" +#include "grt/travt/travt.h" + +char ** grt_string_split(const char *string, const char *delim, int *size) { char *str_copy = strdup(string); // 创建字符串副本,以免修改原始字符串 char *token = strtok(str_copy, delim); @@ -41,7 +41,7 @@ char ** string_split(const char *string, const char *delim, int *size) } -const char* get_basename(const char* path) { +const char* grt_get_basename(const char* path) { // 找到最后一个 '/' char* last_slash = strrchr(path, '/'); @@ -78,7 +78,7 @@ const char* get_basename(const char* path) { */ static void write_one_to_sac( const char *srcname, const char ch, const MYREAL delayT, - const MYREAL wI, FFTW_HOLDER *pt_fh, + const MYREAL wI, GRT_FFTW_HOLDER *pt_fh, SACHEAD *pt_hd, const char *s_output_subdir, const char *s_prefix, const int sgn, const MYCOMPLEX *grncplx) { @@ -88,7 +88,7 @@ static void write_one_to_sac( GRT_SAFE_ASPRINTF(&s_outpath, "%s/%s.sac", s_output_subdir, pt_hd->kcmpnm); // 执行fft任务会修改数组,需重新置零 - reset_fftw_holder_zero(pt_fh); + grt_reset_fftw_holder_zero(pt_fh); // 赋值复数,包括时移 MYCOMPLEX cfac, ccoef; @@ -105,7 +105,7 @@ static void write_one_to_sac( // 发起fft任务 fftw_execute(pt_fh->plan); } else { - naive_inverse_transform_double(pt_fh); + grt_naive_inverse_transform_double(pt_fh); } @@ -156,8 +156,8 @@ static void write_one_to_sac( * */ static void single_freq2time_write_to_file( - const char *command, const PYMODEL1D *pymod, const char *s_prefix, - const MYREAL wI, FFTW_HOLDER *pt_fh, + const char *command, const GRT_PYMODEL1D *pymod, const char *s_prefix, + const MYREAL wI, GRT_FFTW_HOLDER *pt_fh, const char *s_dist, const MYREAL dist, const MYREAL depsrc, const MYREAL deprcv, const MYREAL delayT0, const MYREAL delayV0, const bool calc_upar, @@ -180,9 +180,9 @@ static void single_freq2time_write_to_file( pt_hd->b = delayT; // 计算理论走时 - pt_hd->t0 = compute_travt1d(pymod->Thk, pymod->Va, pymod->n, pymod->isrc, pymod->ircv, dist); + pt_hd->t0 = grt_compute_travt1d(pymod->Thk, pymod->Va, pymod->n, pymod->isrc, pymod->ircv, dist); strcpy(pt_hd->kt0, "P"); - pt_hd->t1 = compute_travt1d(pymod->Thk, pymod->Vb, pymod->n, pymod->isrc, pymod->ircv, dist); + pt_hd->t1 = grt_compute_travt1d(pymod->Thk, pymod->Vb, pymod->n, pymod->isrc, pymod->ircv, dist); strcpy(pt_hd->kt1, "S"); for(int im=0; im #include -#include "dynamic/grn.h" -#include "dynamic/propagate.h" -#include "common/ptam.h" -#include "common/fim.h" -#include "common/safim.h" -#include "common/dwm.h" -#include "common/integral.h" -#include "common/iostats.h" -#include "common/const.h" -#include "common/model.h" -#include "common/prtdbg.h" -#include "common/search.h" -#include "common/progressbar.h" +#include "grt/dynamic/grn.h" +#include "grt/dynamic/propagate.h" +#include "grt/common/ptam.h" +#include "grt/common/fim.h" +#include "grt/common/safim.h" +#include "grt/common/dwm.h" +#include "grt/common/integral.h" +#include "grt/common/iostats.h" +#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" /** @@ -52,7 +52,7 @@ static void recordin_GRN( MYCOMPLEX (*tmp_grn)[SRC_M_NUM][CHANNEL_NUM] = (MYCOMPLEX(*)[SRC_M_NUM][CHANNEL_NUM])calloc(nr, sizeof(*tmp_grn)); for(MYINT ir=0; ir mod1d - MODEL1D *main_mod1d = init_mod1d(pymod1d->n); - get_mod1d(pymod1d, main_mod1d); + GRT_MODEL1D *main_mod1d = grt_init_mod1d(pymod1d->n); + grt_get_mod1d(pymod1d, main_mod1d); const LAYER *src_lay = main_mod1d->lays + main_mod1d->isrc; const MYREAL Rho = src_lay->Rho; // 震源区密度 @@ -158,18 +158,18 @@ void integ_grn_spec( MYCOMPLEX (*sum_uiz_J)[SRC_M_NUM][INTEG_NUM] = (calc_upar)? (MYCOMPLEX(*)[SRC_M_NUM][INTEG_NUM])calloc(nr, sizeof(*sum_uiz_J)) : NULL; MYCOMPLEX (*sum_uir_J)[SRC_M_NUM][INTEG_NUM] = (calc_upar)? (MYCOMPLEX(*)[SRC_M_NUM][INTEG_NUM])calloc(nr, sizeof(*sum_uir_J)) : NULL; - MODEL1D *local_mod1d = NULL; + GRT_MODEL1D *local_mod1d = NULL; #ifdef _OPENMP // 定义局部模型对象 - local_mod1d = init_mod1d(main_mod1d->n); - copy_mod1d(main_mod1d, local_mod1d); + local_mod1d = grt_init_mod1d(main_mod1d->n); + grt_copy_mod1d(main_mod1d, local_mod1d); #else local_mod1d = main_mod1d; #endif - update_mod1d_omega(local_mod1d, omega); + grt_update_mod1d_omega(local_mod1d, omega); // 是否要输出积分过程文件 - bool needfstats = (statsstr!=NULL && ((findElement_MYINT(statsidxs, nstatsidxs, iw) >= 0) || (findElement_MYINT(statsidxs, nstatsidxs, -1) >= 0))); + bool needfstats = (statsstr!=NULL && ((grt_findElement_MYINT(statsidxs, nstatsidxs, iw) >= 0) || (grt_findElement_MYINT(statsidxs, nstatsidxs, -1) >= 0))); // 为当前频率创建波数积分记录文件 FILE *fstats = NULL; @@ -215,35 +215,35 @@ void integ_grn_spec( freq_invstats[iw]=INVERSE_SUCCESS; // 常规的波数积分 - k = discrete_integ( + k = grt_discrete_integ( local_mod1d, dk, (useFIM)? filonK : kmax, keps, omega, nr, rs, sum_J, calc_upar, sum_uiz_J, sum_uir_J, - fstats, kernel, &freq_invstats[iw]); + fstats, grt_kernel, &freq_invstats[iw]); // 使用Filon积分 if(useFIM && freq_invstats[iw]==INVERSE_SUCCESS){ if(filondk > RZERO){ // 基于线性插值的Filon积分,固定采样间隔 - k = linear_filon_integ( + k = grt_linear_filon_integ( local_mod1d, k, dk, filondk, kmax, keps, omega, nr, rs, sum_J, calc_upar, sum_uiz_J, sum_uir_J, - fstats, kernel, &freq_invstats[iw]); + fstats, grt_kernel, &freq_invstats[iw]); } else if(safilonTol > RZERO){ // 基于自适应采样的Filon积分 - k = sa_filon_integ( + k = grt_sa_filon_integ( local_mod1d, fabs(vmin_ref)/ampk, k, dk, safilonTol, kmax, omega, nr, rs, sum_J, calc_upar, sum_uiz_J, sum_uir_J, - fstats, kernel, &freq_invstats[iw]); + fstats, grt_kernel, &freq_invstats[iw]); } } // k之后的部分使用峰谷平均法进行显式收敛,建议在浅源地震的时候使用 if(vmin_ref < RZERO && freq_invstats[iw]==INVERSE_SUCCESS){ - PTA_method( + grt_PTA_method( local_mod1d, k, dk, omega, nr, rs, sum_J, calc_upar, sum_uiz_J, sum_uir_J, - ptam_fstatsnr, kernel, &freq_invstats[iw]); + ptam_fstatsnr, grt_kernel, &freq_invstats[iw]); } // fprintf(stderr, "iw=%d, w=%.5e, k=%.5e, dk=%.5e, nk=%d\n", iw, w, k, dk, (int)(k/dk)); @@ -271,14 +271,14 @@ void integ_grn_spec( GRT_SAFE_FREE_PTR(ptam_fstatsnr); #ifdef _OPENMP - free_mod1d(local_mod1d); + grt_free_mod1d(local_mod1d); #endif // 记录进度条变量 #pragma omp critical { progress++; - if(print_progressbar) printprogressBar("Computing Green Functions: ", progress*100/(nf2-nf1+1)); + if(print_progressbar) grt_printprogressBar("Computing Green Functions: ", progress*100/(nf2-nf1+1)); } @@ -291,7 +291,7 @@ void integ_grn_spec( - free_mod1d(main_mod1d); + grt_free_mod1d(main_mod1d); GRT_SAFE_FREE_PTR_ARRAY(ptam_fstatsdir, nr); diff --git a/pygrt/C_extension/src/dynamic/grt_greenfn.c b/pygrt/C_extension/src/dynamic/grt_greenfn.c index bebb331..05dd020 100644 --- a/pygrt/C_extension/src/dynamic/grt_greenfn.c +++ b/pygrt/C_extension/src/dynamic/grt_greenfn.c @@ -9,14 +9,14 @@ */ -#include "dynamic/grn.h" -#include "dynamic/signals.h" -#include "common/const.h" -#include "common/model.h" -#include "common/search.h" -#include "common/sacio.h" -#include "common/util.h" -#include "common/myfftw.h" +#include "grt/dynamic/grn.h" +#include "grt/dynamic/signals.h" +#include "grt/common/const.h" +#include "grt/common/model.h" +#include "grt/common/search.h" +#include "grt/common/sacio.h" +#include "grt/common/util.h" +#include "grt/common/myfftw.h" #include "grt.h" @@ -44,7 +44,7 @@ typedef struct { bool active; char *s_modelpath; ///< 模型路径 const char *s_modelname; ///< 模型名称 - PYMODEL1D *pymod; ///< 模型PYMODEL1D结构体指针 + GRT_PYMODEL1D *pymod; ///< 模型PYMODEL1D结构体指针 } M; /** 震源和接收器深度 */ struct { @@ -153,7 +153,7 @@ static void free_Ctrl(GRT_MODULE_CTRL *Ctrl){ // M GRT_SAFE_FREE_PTR(Ctrl->M.s_modelpath); - free_pymod(Ctrl->M.pymod); + grt_free_pymod(Ctrl->M.pymod); // D GRT_SAFE_FREE_PTR(Ctrl->D.s_depsrc); @@ -184,7 +184,7 @@ static void free_Ctrl(GRT_MODULE_CTRL *Ctrl){ /** 打印结构体中的参数 */ static void print_Ctrl(const GRT_MODULE_CTRL *Ctrl){ - print_pymod(Ctrl->M.pymod); + grt_print_pymod(Ctrl->M.pymod); const char format[] = " \%-20s \%s\n"; const char format_real[] = " \%-20s \%.3f\n"; @@ -455,7 +455,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ case 'M': Ctrl->M.active = true; Ctrl->M.s_modelpath = strdup(optarg); - Ctrl->M.s_modelname = get_basename(Ctrl->M.s_modelpath); + Ctrl->M.s_modelname = grt_get_basename(Ctrl->M.s_modelpath); break; // 震源和场点深度, -Ddepsrc/deprcv @@ -600,7 +600,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ case 'R': Ctrl->R.active = true; Ctrl->R.s_raw = strdup(optarg); - Ctrl->R.s_rs = string_split(optarg, ",", &Ctrl->R.nr); + Ctrl->R.s_rs = grt_string_split(optarg, ",", &Ctrl->R.nr); // 转为浮点数 Ctrl->R.rs = (MYREAL*)realloc(Ctrl->R.rs, sizeof(MYREAL)*(Ctrl->R.nr)); for(MYINT i=0; iR.nr; ++i){ @@ -620,7 +620,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ if(Ctrl->P.nthreads <= 0){ GRTBadOptionError(command, P, "Nonpositive value is not supported."); } - set_num_threads(Ctrl->P.nthreads); + grt_set_num_threads(Ctrl->P.nthreads); break; // 选择要计算的格林函数 -G1/1/1/1 @@ -649,7 +649,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ case 'S': Ctrl->S.active = true; Ctrl->S.s_raw = strdup(optarg); - Ctrl->S.s_statsidxs = string_split(optarg, ",", &Ctrl->S.nstatsidxs); + Ctrl->S.s_statsidxs = grt_string_split(optarg, ",", &Ctrl->S.nstatsidxs); // 转为浮点数 Ctrl->S.statsidxs = (MYINT*)realloc(Ctrl->S.statsidxs, sizeof(MYINT)*(Ctrl->S.nstatsidxs)); for(MYINT i=0; iS.nstatsidxs; ++i){ @@ -704,10 +704,10 @@ int greenfn_main(int argc, char **argv) { getopt_from_command(Ctrl, argc, argv); // 读入模型文件 - if((Ctrl->M.pymod = read_pymod_from_file(command, Ctrl->M.s_modelpath, Ctrl->D.depsrc, Ctrl->D.deprcv, true)) == NULL){ + if((Ctrl->M.pymod = grt_read_pymod_from_file(command, Ctrl->M.s_modelpath, Ctrl->D.depsrc, Ctrl->D.deprcv, true)) == NULL){ exit(EXIT_FAILURE); } - PYMODEL1D *pymod = Ctrl->M.pymod; + GRT_PYMODEL1D *pymod = Ctrl->M.pymod; // 当震源位于液体层中时,仅允许计算爆炸源对应的格林函数 // 程序结束前会输出对应警告 @@ -717,7 +717,7 @@ int greenfn_main(int argc, char **argv) { // 最大最小速度 MYREAL vmin, vmax; - get_pymod_vmin_vmax(pymod, &vmin, &vmax); + grt_get_pymod_vmin_vmax(pymod, &vmin, &vmax); // 参考最小速度 if(!Ctrl->V.active){ @@ -733,7 +733,7 @@ int greenfn_main(int argc, char **argv) { Ctrl->N.winT = Ctrl->N.nt*Ctrl->N.dt; // 最大震中距 - MYREAL rmax = Ctrl->R.rs[findMinMax_MYREAL(Ctrl->R.rs, Ctrl->R.nr, true)]; + MYREAL rmax = Ctrl->R.rs[grt_findMinMax_MYREAL(Ctrl->R.rs, Ctrl->R.nr, true)]; // 时窗最大截止时刻 MYREAL tmax = Ctrl->E.delayT0 + Ctrl->N.winT; @@ -805,7 +805,7 @@ int greenfn_main(int argc, char **argv) { //============================================================================== // 计算格林函数 - integ_grn_spec( + grt_integ_grn_spec( pymod, Ctrl->H.nf1, Ctrl->H.nf2, Ctrl->N.freqs, Ctrl->R.nr, Ctrl->R.rs, Ctrl->N.wI, Ctrl->V.vmin_ref, Ctrl->K.keps, Ctrl->K.ampk, Ctrl->K.k0, Ctrl->L.Length, Ctrl->L.filonLength, Ctrl->L.safilonTol, Ctrl->L.filonCut, !Ctrl->s.active, grn, Ctrl->e.active, grn_uiz, grn_uir, @@ -816,11 +816,11 @@ int greenfn_main(int argc, char **argv) { // 使用fftw3做反傅里叶变换,并保存到 SAC // 其中考虑了升采样倍数 - FFTW_HOLDER *fftw_holder = create_fftw_holder_C2R_1D( + GRT_FFTW_HOLDER *fftw_holder = grt_create_fftw_holder_C2R_1D( Ctrl->N.nt*Ctrl->N.upsample_n, Ctrl->N.dt/Ctrl->N.upsample_n, Ctrl->N.nf, Ctrl->N.df); MYREAL (* travtPS)[2] = (MYREAL (*)[2])calloc(Ctrl->R.nr, sizeof(MYREAL)*2); - GF_freq2time_write_to_file( + grt_GF_freq2time_write_to_file( command, pymod, Ctrl->O.s_output_dir, Ctrl->M.s_modelname, Ctrl->D.s_depsrc, Ctrl->D.s_deprcv, Ctrl->N.wI, fftw_holder, @@ -857,7 +857,7 @@ int greenfn_main(int argc, char **argv) { GRT_SAFE_FREE_PTR(grn_uir); GRT_SAFE_FREE_PTR(travtPS); - destroy_fftw_holder(fftw_holder); + grt_destroy_fftw_holder(fftw_holder); free_Ctrl(Ctrl); return EXIT_SUCCESS; diff --git a/pygrt/C_extension/src/dynamic/grt_rotation.c b/pygrt/C_extension/src/dynamic/grt_rotation.c index 25bfd43..0957972 100644 --- a/pygrt/C_extension/src/dynamic/grt_rotation.c +++ b/pygrt/C_extension/src/dynamic/grt_rotation.c @@ -8,8 +8,8 @@ */ -#include "common/sacio2.h" -#include "common/const.h" +#include "grt/common/sacio2.h" +#include "grt/common/const.h" #include "grt.h" @@ -106,7 +106,7 @@ int rotation_main(int argc, char **argv){ // 读取一个头段变量,获得基本参数,分配数组内存 SACHEAD hd; GRT_SAFE_ASPRINTF(&s_filepath, "%s/%c%s%c.sac", Ctrl->s_synpath, tolower(chs[0]), Ctrl->s_prefix, chs[0]); - read_SAC_HEAD(command, s_filepath, &hd); + grt_read_SAC_HEAD(command, s_filepath, &hd); int npts=hd.npts; float dist=hd.dist; float *arrout = (float*)calloc(npts, sizeof(float)); @@ -120,14 +120,14 @@ int rotation_main(int argc, char **argv){ // 读取数据 u_{i,j} GRT_SAFE_ASPRINTF(&s_filepath, "%s/%c%s%c.sac", Ctrl->s_synpath, tolower(c2), Ctrl->s_prefix, c1); - arrin = read_SAC(command, s_filepath, &hd, arrin); + arrin = grt_read_SAC(command, s_filepath, &hd, arrin); // 累加 for(int i=0; is_synpath, tolower(c1), Ctrl->s_prefix, c2); - arrin = read_SAC(command, s_filepath, &hd, arrin); + arrin = grt_read_SAC(command, s_filepath, &hd, arrin); // 累加 for(int i=0; is_synpath, Ctrl->s_prefix); - arrin = read_SAC(command, s_filepath, &hd, arrin); + arrin = grt_read_SAC(command, s_filepath, &hd, arrin); for(int i=0; is_synpath, tolower(chs[0]), Ctrl->s_prefix, chs[0]); - read_SAC_HEAD(command, s_filepath, &hd); + grt_read_SAC_HEAD(command, s_filepath, &hd); int npts=hd.npts; float dist=hd.dist; float *arrout = (float*)calloc(npts, sizeof(float)); @@ -119,14 +119,14 @@ int strain_main(int argc, char **argv){ // 读取数据 u_{i,j} GRT_SAFE_ASPRINTF(&s_filepath, "%s/%c%s%c.sac", Ctrl->s_synpath, tolower(c2), Ctrl->s_prefix, c1); - arrin = read_SAC(command, s_filepath, &hd, arrin); + arrin = grt_read_SAC(command, s_filepath, &hd, arrin); // 累加 for(int i=0; is_synpath, tolower(c1), Ctrl->s_prefix, c2); - arrin = read_SAC(command, s_filepath, &hd, arrin); + arrin = grt_read_SAC(command, s_filepath, &hd, arrin); // 累加 for(int i=0; is_synpath, Ctrl->s_prefix); - arrin = read_SAC(command, s_filepath, &hd, arrin); + arrin = grt_read_SAC(command, s_filepath, &hd, arrin); for(int i=0; is_synpath, Ctrl->s_prefix); - arrin = read_SAC(command, s_filepath, &hd, arrin); + arrin = grt_read_SAC(command, s_filepath, &hd, arrin); for(int i=0; is_synpath, tolower(chs[0]), Ctrl->s_prefix, chs[0]); - read_SAC_HEAD(command, s_filepath, &hd); + grt_read_SAC_HEAD(command, s_filepath, &hd); int npts=hd.npts; float dt=hd.delta; float dist=hd.dist; @@ -129,8 +129,8 @@ int stress_main(int argc, char **argv){ fftwf_complex *lams = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex)*nf); fftwf_complex *mus = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex)*nf); // 分配FFTW - FFTWF_HOLDER *fwd_fftw_holder = create_fftwf_holder_R2C_1D(npts, dt, nf, df); - FFTWF_HOLDER *inv_fftw_holder = create_fftwf_holder_C2R_1D(npts, dt, nf, df); + GRT_FFTWF_HOLDER *fwd_fftw_holder = grt_create_fftwf_holder_R2C_1D(npts, dt, nf, df); + GRT_FFTWF_HOLDER *inv_fftw_holder = grt_create_fftwf_holder_C2R_1D(npts, dt, nf, df); // 初始化 memset(lam_ukk, 0, sizeof(fftwf_complex)*nf); memset(lams, 0, sizeof(fftwf_complex)*nf); @@ -141,8 +141,8 @@ int stress_main(int argc, char **argv){ freq = (i==0) ? 0.01f : df*i; // 计算衰减因子不能为0频 w = PI2 * freq; fftwf_complex atta, attb; - atta = attenuation_law(Qainv, w); - attb = attenuation_law(Qbinv, w); + atta = grt_attenuation_law(Qainv, w); + attb = grt_attenuation_law(Qbinv, w); // 乘上1e10,转为dyne/(cm^2) mus[i] = vb*vb*attb*attb*rho*1e10; lams[i] = va*va*atta*atta*rho*1e10 - 2.0*mus[i]; @@ -155,7 +155,7 @@ int stress_main(int argc, char **argv){ // 读取数据 u_{k,k} GRT_SAFE_ASPRINTF(&s_filepath, "%s/%c%s%c.sac", Ctrl->s_synpath, tolower(c1), Ctrl->s_prefix, c1); - read_SAC(command, s_filepath, &hd, fwd_fftw_holder->w_t); + grt_read_SAC(command, s_filepath, &hd, fwd_fftw_holder->w_t); // 累加 fftwf_execute(fwd_fftw_holder->plan); @@ -164,7 +164,7 @@ int stress_main(int argc, char **argv){ // 加上协变导数 if(!rot2ZNE){ GRT_SAFE_ASPRINTF(&s_filepath, "%s/%sR.sac", Ctrl->s_synpath, Ctrl->s_prefix); - read_SAC(command, s_filepath, &hd, fwd_fftw_holder->w_t); + grt_read_SAC(command, s_filepath, &hd, fwd_fftw_holder->w_t); fftwf_execute(fwd_fftw_holder->plan); for(int i=0; iW_f[i]/dist*1e-5; } @@ -173,8 +173,8 @@ int stress_main(int argc, char **argv){ for(int i=0; is_synpath, tolower(c2), Ctrl->s_prefix, c1); - read_SAC(command, s_filepath, &hd, fwd_fftw_holder->w_t); + grt_read_SAC(command, s_filepath, &hd, fwd_fftw_holder->w_t); // 累加 fftwf_execute(fwd_fftw_holder->plan); @@ -193,7 +193,7 @@ int stress_main(int argc, char **argv){ // 读取数据 u_{j,i} GRT_SAFE_ASPRINTF(&s_filepath, "%s/%c%s%c.sac", Ctrl->s_synpath, tolower(c1), Ctrl->s_prefix, c2); - read_SAC(command, s_filepath, &hd, fwd_fftw_holder->w_t); + grt_read_SAC(command, s_filepath, &hd, fwd_fftw_holder->w_t); // 累加 fftwf_execute(fwd_fftw_holder->plan); @@ -208,14 +208,14 @@ int stress_main(int argc, char **argv){ if(c1=='R' && c2=='T'){ // 读取数据 u_T GRT_SAFE_ASPRINTF(&s_filepath, "%s/%sT.sac", Ctrl->s_synpath, Ctrl->s_prefix); - read_SAC(command, s_filepath, &hd, fwd_fftw_holder->w_t); + grt_read_SAC(command, s_filepath, &hd, fwd_fftw_holder->w_t); fftwf_execute(fwd_fftw_holder->plan); for(int i=0; iW_f[i] -= mus[i] * fwd_fftw_holder->W_f[i] / dist * 1e-5; } else if(c1=='T' && c2=='T'){ // 读取数据 u_R GRT_SAFE_ASPRINTF(&s_filepath, "%s/%sR.sac", Ctrl->s_synpath, Ctrl->s_prefix); - read_SAC(command, s_filepath, &hd, fwd_fftw_holder->w_t); + grt_read_SAC(command, s_filepath, &hd, fwd_fftw_holder->w_t); fftwf_execute(fwd_fftw_holder->plan); for(int i=0; iW_f[i] += 2.0f * mus[i] * fwd_fftw_holder->W_f[i] / dist * 1e-5; } @@ -228,14 +228,14 @@ int stress_main(int argc, char **argv){ write_sac(s_filepath, hd, inv_fftw_holder->w_t); // 置零 - reset_fftwf_holder_zero(inv_fftw_holder); + grt_reset_fftwf_holder_zero(inv_fftw_holder); } } - destroy_fftwf_holder(fwd_fftw_holder); - destroy_fftwf_holder(inv_fftw_holder); + grt_destroy_fftwf_holder(fwd_fftw_holder); + grt_destroy_fftwf_holder(inv_fftw_holder); GRT_SAFE_FFTW_FREE_PTR(lam_ukk, f); GRT_SAFE_FFTW_FREE_PTR(lams, f); diff --git a/pygrt/C_extension/src/dynamic/grt_syn.c b/pygrt/C_extension/src/dynamic/grt_syn.c index e879516..c465887 100644 --- a/pygrt/C_extension/src/dynamic/grt_syn.c +++ b/pygrt/C_extension/src/dynamic/grt_syn.c @@ -7,11 +7,11 @@ * */ -#include "dynamic/signals.h" -#include "common/sacio2.h" -#include "common/const.h" -#include "common/radiation.h" -#include "common/coord.h" +#include "grt/dynamic/signals.h" +#include "grt/common/sacio2.h" +#include "grt/common/const.h" +#include "grt/common/radiation.h" +#include "grt/common/coord.h" #include "grt.h" @@ -269,7 +269,7 @@ static void check_grn_exist(GRT_MODULE_CTRL *Ctrl, const char *name){ // 检查文件的同时将src_mu计算出来 if(Ctrl->S.src_mu == 0.0 && Ctrl->S.mult_src_mu){ SACHEAD hd; - read_SAC_HEAD(command, buffer, &hd); + grt_read_SAC_HEAD(command, buffer, &hd); double va, vb, rho; va = hd.user6; vb = hd.user7; @@ -417,7 +417,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ GRTBadOptionError(command, D, ""); } // 检查测试 - if(! check_tftype_tfparams(Ctrl->D.tftype, Ctrl->D.tfparams)){ + if(! grt_check_tftype_tfparams(Ctrl->D.tftype, Ctrl->D.tfparams)){ GRTBadOptionError(command, D, ""); } break; @@ -610,9 +610,9 @@ static void data_zrt2zne(float *syn[3], float *syn_upar[3][3], int nt, double az } if(doupar) { - rot_zrt2zxy_upar(azrad, dblsyn, dblupar, dist*1e5); + grt_rot_zrt2zxy_upar(azrad, dblsyn, dblupar, dist*1e5); } else { - rot_zxy2zrt_vec(-azrad, dblsyn); + grt_rot_zxy2zrt_vec(-azrad, dblsyn); } @@ -684,7 +684,7 @@ int syn_main(int argc, char **argv){ } // 重新计算方向因子 - set_source_radiation(Ctrl->srcRadi, Ctrl->computeType, (ityp==3), Ctrl->S.M0, upar_scale, Ctrl->A.azrad, Ctrl->mchn); + grt_set_source_radiation(Ctrl->srcRadi, Ctrl->computeType, (ityp==3), Ctrl->S.M0, upar_scale, Ctrl->A.azrad, Ctrl->mchn); for(int c=0; cG.s_grnpath, tolower(ZRTchs[ityp-1]), SRC_M_NAME_ABBR[k], ch); } - float *arr = read_SAC(command, buffer, pthd, NULL); + float *arr = grt_read_SAC(command, buffer, pthd, NULL); hd0 = *pthd; // 备份一份 nt = pthd->npts; @@ -746,7 +746,7 @@ int syn_main(int argc, char **argv){ if(Ctrl->D.active && Ctrl->D.tfarr==NULL){ // 获得时间函数 - Ctrl->D.tfarr = get_time_function(&Ctrl->D.tfnt, dt, Ctrl->D.tftype, Ctrl->D.tfparams); + Ctrl->D.tfarr = grt_get_time_function(&Ctrl->D.tfnt, dt, Ctrl->D.tftype, Ctrl->D.tfparams); if(Ctrl->D.tfarr==NULL){ GRTRaiseError("[%s] get time function error.\n", command); } @@ -762,7 +762,7 @@ int syn_main(int argc, char **argv){ // 时域循环卷积 if(Ctrl->D.tfarr!=NULL){ float *convarr = (float*)calloc(nt, sizeof(float)); - oaconvolve(arrout, nt, Ctrl->D.tfarr, Ctrl->D.tfnt, convarr, nt, true); + grt_oaconvolve(arrout, nt, Ctrl->D.tfarr, Ctrl->D.tfnt, convarr, nt, true); for(int i=0; iI.int_times; ++i){ - trap_integral(arrout, nt, dt); + grt_trap_integral(arrout, nt, dt); } for(int i=0; iJ.dif_times; ++i){ - differential(arrout, nt, dt); + grt_differential(arrout, nt, dt); } } // ENDFOR 三分量 diff --git a/pygrt/C_extension/src/dynamic/layer.c b/pygrt/C_extension/src/dynamic/layer.c index e16f6b5..97ca0a1 100644 --- a/pygrt/C_extension/src/dynamic/layer.c +++ b/pygrt/C_extension/src/dynamic/layer.c @@ -16,14 +16,14 @@ #include #include -#include "dynamic/layer.h" -#include "common/model.h" -#include "common/prtdbg.h" -#include "common/matrix.h" -#include "common/colorstr.h" +#include "grt/dynamic/layer.h" +#include "grt/common/model.h" +#include "grt/common/prtdbg.h" +#include "grt/common/matrix.h" +#include "grt/common/colorstr.h" -void calc_R_tilt_PSV(MYCOMPLEX xa0, MYCOMPLEX xb0, MYCOMPLEX kbkb0, MYREAL k, MYCOMPLEX R_tilt[2][2], MYINT *stats) +void grt_calc_R_tilt_PSV(MYCOMPLEX xa0, MYCOMPLEX xb0, MYCOMPLEX kbkb0, MYREAL k, MYCOMPLEX R_tilt[2][2], MYINT *stats) { if(kbkb0 != CZERO){ // 固体表面 @@ -52,7 +52,7 @@ void calc_R_tilt_PSV(MYCOMPLEX xa0, MYCOMPLEX xb0, MYCOMPLEX kbkb0, MYREAL k, MY } -void calc_R_EV_PSV( +void grt_calc_R_EV_PSV( MYCOMPLEX xa_rcv, MYCOMPLEX xb_rcv, bool ircvup, MYREAL k, const MYCOMPLEX R[2][2], MYCOMPLEX R_EV[2][2]) { @@ -74,16 +74,16 @@ void calc_R_EV_PSV( // 公式(5.7.7,25) if(ircvup){// 震源更深 - cmat2x2_mul(D12, R, R_EV); - cmat2x2_add(D11, R_EV, R_EV); + grt_cmat2x2_mul(D12, R, R_EV); + grt_cmat2x2_add(D11, R_EV, R_EV); } else { // 接收点更深 - cmat2x2_mul(D11, R, R_EV); - cmat2x2_add(D12, R_EV, R_EV); + grt_cmat2x2_mul(D11, R, R_EV); + grt_cmat2x2_add(D12, R_EV, R_EV); } } -void calc_R_EV_SH(MYCOMPLEX xb_rcv, MYREAL k, MYCOMPLEX RL, MYCOMPLEX *R_EVL) +void grt_calc_R_EV_SH(MYCOMPLEX xb_rcv, MYREAL k, MYCOMPLEX RL, MYCOMPLEX *R_EVL) { if(xb_rcv != CONE){ // 位于固体层 @@ -96,7 +96,7 @@ void calc_R_EV_SH(MYCOMPLEX xb_rcv, MYREAL k, MYCOMPLEX RL, MYCOMPLEX *R_EVL) } -void calc_uiz_R_EV_PSV( +void grt_calc_uiz_R_EV_PSV( MYCOMPLEX xa_rcv, MYCOMPLEX xb_rcv, bool ircvup, MYREAL k, const MYCOMPLEX R[2][2], MYCOMPLEX R_EV[2][2]) @@ -117,16 +117,16 @@ void calc_uiz_R_EV_PSV( // 公式(5.7.7,25) if(ircvup){// 震源更深 - cmat2x2_mul(D12, R, R_EV); - cmat2x2_add(D11, R_EV, R_EV); + grt_cmat2x2_mul(D12, R, R_EV); + grt_cmat2x2_add(D11, R_EV, R_EV); } else { // 接收点更深 - cmat2x2_mul(D11, R, R_EV); - cmat2x2_add(D12, R_EV, R_EV); + grt_cmat2x2_mul(D11, R, R_EV); + grt_cmat2x2_add(D12, R_EV, R_EV); } } -void calc_uiz_R_EV_SH(MYCOMPLEX xb_rcv, bool ircvup, MYREAL k, MYCOMPLEX RL, MYCOMPLEX *R_EVL) +void grt_calc_uiz_R_EV_SH(MYCOMPLEX xb_rcv, bool ircvup, MYREAL k, MYCOMPLEX RL, MYCOMPLEX *R_EVL) { // 将势函数转为ui,z在(B_m, P_m, C_m)系下的分量 // 新推导的公式 @@ -146,7 +146,7 @@ void calc_uiz_R_EV_SH(MYCOMPLEX xb_rcv, bool ircvup, MYREAL k, MYCOMPLEX RL, MYC } -void calc_RT_ll_PSV( +void grt_calc_RT_ll_PSV( MYREAL Rho1, MYCOMPLEX xa1, MYREAL Rho2, MYCOMPLEX xa2, MYREAL thk, // 使用上层的厚度 @@ -179,7 +179,7 @@ void calc_RT_ll_PSV( } -void calc_RT_ll_SH(MYCOMPLEX *RDL, MYCOMPLEX *RUL, MYCOMPLEX *TDL, MYCOMPLEX *TUL) +void grt_calc_RT_ll_SH(MYCOMPLEX *RDL, MYCOMPLEX *RUL, MYCOMPLEX *TDL, MYCOMPLEX *TUL) { *RDL = RZERO; *RUL = RZERO; @@ -189,7 +189,7 @@ void calc_RT_ll_SH(MYCOMPLEX *RDL, MYCOMPLEX *RUL, MYCOMPLEX *TDL, MYCOMPLEX *TU -void calc_RT_ls_PSV( +void grt_calc_RT_ls_PSV( MYREAL Rho1, MYCOMPLEX xa1, MYCOMPLEX xb1, MYCOMPLEX kbkb1, MYCOMPLEX mu1, MYREAL Rho2, MYCOMPLEX xa2, MYCOMPLEX xb2, MYCOMPLEX kbkb2, MYCOMPLEX mu2, MYREAL thk, // 使用上层的厚度 @@ -280,7 +280,7 @@ void calc_RT_ls_PSV( } -void calc_RT_ls_SH( +void grt_calc_RT_ls_SH( MYCOMPLEX xb1, MYCOMPLEX mu1, MYCOMPLEX mu2, MYREAL thk, // 使用上层的厚度 MYCOMPLEX omega, MYREAL k, @@ -326,7 +326,7 @@ void calc_RT_ls_SH( -void calc_RT_ss_PSV( +void grt_calc_RT_ss_PSV( MYREAL Rho1, MYCOMPLEX xa1, MYCOMPLEX xb1, MYCOMPLEX kbkb1, MYCOMPLEX mu1, MYREAL Rho2, MYCOMPLEX xa2, MYCOMPLEX xb2, MYCOMPLEX kbkb2, MYCOMPLEX mu2, MYREAL thk, // 使用上层的厚度 @@ -411,7 +411,7 @@ void calc_RT_ss_PSV( } -void calc_RT_ss_SH( +void grt_calc_RT_ss_SH( MYCOMPLEX xb1, MYCOMPLEX mu1, MYCOMPLEX xb2, MYCOMPLEX mu2, MYREAL thk, // 使用上层的厚度 @@ -437,7 +437,7 @@ void calc_RT_ss_SH( -void calc_RT_PSV( +void grt_calc_RT_PSV( MYREAL Rho1, MYCOMPLEX xa1, MYCOMPLEX xb1, MYCOMPLEX kbkb1, MYCOMPLEX mu1, MYREAL Rho2, MYCOMPLEX xa2, MYCOMPLEX xb2, MYCOMPLEX kbkb2, MYCOMPLEX mu2, MYREAL thk, // 使用上层的厚度 @@ -447,21 +447,21 @@ void calc_RT_PSV( { // 根据界面两侧的具体情况选择函数 if(mu1 != CZERO && mu2 != CZERO){ - calc_RT_ss_PSV( + grt_calc_RT_ss_PSV( Rho1, xa1, xb1, kbkb1, mu1, Rho2, xa2, xb2, kbkb2, mu2, thk, omega, k, RD, RU, TD, TU, stats); } else if(mu1 == CZERO && mu2 == CZERO){ - calc_RT_ll_PSV( + grt_calc_RT_ll_PSV( Rho1, xa1, Rho2, xa2, thk, omega, k, RD, RU, TD, TU, stats); } else{ - calc_RT_ls_PSV( + grt_calc_RT_ls_PSV( Rho1, xa1, xb1, kbkb1, mu1, Rho2, xa2, xb2, kbkb2, mu2, thk, omega, k, @@ -470,7 +470,7 @@ void calc_RT_PSV( } -void calc_RT_SH( +void grt_calc_RT_SH( MYCOMPLEX xb1, MYCOMPLEX mu1, MYCOMPLEX xb2, MYCOMPLEX mu2, MYREAL thk, // 使用上层的厚度 @@ -480,18 +480,18 @@ void calc_RT_SH( { // 根据界面两侧的具体情况选择函数 if(mu1 != CZERO && mu2 != CZERO){ - calc_RT_ss_SH( + grt_calc_RT_ss_SH( xb1, mu1, xb2, mu2, thk, omega, k, RDL, RUL, TDL, TUL); } else if(mu1 == CZERO && mu2 == CZERO){ - calc_RT_ll_SH( + grt_calc_RT_ll_SH( RDL, RUL, TDL, TUL); } else{ - calc_RT_ls_SH( + grt_calc_RT_ls_SH( xb1, mu1, mu2, thk, omega, k, RDL, RUL, TDL, TUL); @@ -501,7 +501,7 @@ void calc_RT_SH( -void get_layer_D( +void grt_get_layer_D( MYCOMPLEX xa, MYCOMPLEX xb, MYCOMPLEX kbkb, MYCOMPLEX mu, MYCOMPLEX omega, MYREAL k, MYCOMPLEX D[4][4], bool inverse) { @@ -527,7 +527,7 @@ void get_layer_D( } } -void get_layer_D11( +void grt_get_layer_D11( MYCOMPLEX xa, MYCOMPLEX xb, MYREAL k, MYCOMPLEX D[2][2]) { // 第iy层物理量 @@ -535,7 +535,7 @@ void get_layer_D11( D[1][0] = k*xa; D[1][1] = k; } -void get_layer_D12( +void grt_get_layer_D12( MYCOMPLEX xa, MYCOMPLEX xb, MYREAL k, MYCOMPLEX D[2][2]) { // 第iy层物理量 @@ -543,7 +543,7 @@ void get_layer_D12( D[1][0] = -k*xa; D[1][1] = k; } -void get_layer_D11_uiz( +void grt_get_layer_D11_uiz( MYCOMPLEX xa, MYCOMPLEX xb, MYREAL k, MYCOMPLEX D[2][2]) { // 第iy层物理量 @@ -554,7 +554,7 @@ void get_layer_D11_uiz( D[1][0] = a*a; D[1][1] = b*k; } -void get_layer_D12_uiz( +void grt_get_layer_D12_uiz( MYCOMPLEX xa, MYCOMPLEX xb, MYREAL k, MYCOMPLEX D[2][2]) { // 第iy层物理量 @@ -565,7 +565,7 @@ void get_layer_D12_uiz( D[1][0] = a*a; D[1][1] = - b*k; } -void get_layer_D21( +void grt_get_layer_D21( MYCOMPLEX xa, MYCOMPLEX xb, MYCOMPLEX kbkb, MYCOMPLEX mu, MYCOMPLEX omega, MYREAL k, MYCOMPLEX D[2][2]) { @@ -577,7 +577,7 @@ void get_layer_D21( D[1][0] = 2*k*mu*k*xa; D[1][1] = 2*mu*Omg; } -void get_layer_D22( +void grt_get_layer_D22( MYCOMPLEX xa, MYCOMPLEX xb, MYCOMPLEX kbkb, MYCOMPLEX mu, MYCOMPLEX omega, MYREAL k, MYCOMPLEX D[2][2]) { @@ -589,7 +589,7 @@ void get_layer_D22( D[1][0] = -2*k*mu*k*xa; D[1][1] = 2*mu*Omg; } -void get_layer_T( +void grt_get_layer_T( MYCOMPLEX xb, MYCOMPLEX mu, MYCOMPLEX omega, MYREAL k, MYCOMPLEX T[2][2], bool inverse) { @@ -607,7 +607,7 @@ void get_layer_T( } } -void get_layer_E_Love(MYCOMPLEX xb1, MYREAL thk, MYREAL k, MYCOMPLEX E[2][2], bool inverse) +void grt_get_layer_E_Love(MYCOMPLEX xb1, MYREAL thk, MYREAL k, MYCOMPLEX E[2][2], bool inverse) { MYCOMPLEX exb = exp(k*thk*xb1); @@ -622,7 +622,7 @@ void get_layer_E_Love(MYCOMPLEX xb1, MYREAL thk, MYREAL k, MYCOMPLEX E[2][2], bo } -void get_layer_E_Rayl( +void grt_get_layer_E_Rayl( MYCOMPLEX xa1, MYCOMPLEX xb1, MYREAL thk, MYREAL k, MYCOMPLEX E[4][4], bool inverse) { MYCOMPLEX exa, exb; @@ -645,7 +645,7 @@ void get_layer_E_Rayl( } } -void calc_RT_from_4x4( +void grt_calc_RT_from_4x4( MYCOMPLEX xa1, MYCOMPLEX xb1, MYCOMPLEX kbkb1, MYCOMPLEX mu1, MYCOMPLEX xa2, MYCOMPLEX xb2, MYCOMPLEX kbkb2, MYCOMPLEX mu2, MYCOMPLEX omega, MYREAL thk, @@ -656,10 +656,10 @@ void calc_RT_from_4x4( MYCOMPLEX D1_inv[4][4], D2[4][4], Q[4][4]; - get_layer_D(xa1, xb1, kbkb1, mu1, omega, k, D1_inv, true); - get_layer_D(xa2, xb2, kbkb2, mu2, omega, k, D2, false); + grt_get_layer_D(xa1, xb1, kbkb1, mu1, omega, k, D1_inv, true); + grt_get_layer_D(xa2, xb2, kbkb2, mu2, omega, k, D2, false); - cmatmxn_mul(4, 4, 4, D1_inv, D2, Q); + grt_cmatmxn_mul(4, 4, 4, D1_inv, D2, Q); MYCOMPLEX exa, exb; @@ -671,26 +671,26 @@ void calc_RT_from_4x4( E[1][1] = exb; E[2][2] = 1/exa; E[3][3] = 1/exb; - cmatmxn_mul(4, 4, 4, E, Q, Q); + grt_cmatmxn_mul(4, 4, 4, E, Q, Q); // 对Q矩阵划分子矩阵 MYCOMPLEX Q11[2][2], Q12[2][2], Q21[2][2], Q22[2][2]; - cmatmxn_block(4, 4, Q, 0, 0, 2, 2, Q11); - cmatmxn_block(4, 4, Q, 0, 2, 2, 2, Q12); - cmatmxn_block(4, 4, Q, 2, 0, 2, 2, Q21); - cmatmxn_block(4, 4, Q, 2, 2, 2, 2, Q22); + grt_cmatmxn_block(4, 4, Q, 0, 0, 2, 2, Q11); + grt_cmatmxn_block(4, 4, Q, 0, 2, 2, 2, Q12); + grt_cmatmxn_block(4, 4, Q, 2, 0, 2, 2, Q21); + grt_cmatmxn_block(4, 4, Q, 2, 2, 2, 2, Q22); // 计算反射透射系数 // TD - cmat2x2_inv(Q22, TD, stats); + grt_cmat2x2_inv(Q22, TD, stats); // RD - cmat2x2_mul(Q12, TD, RD); + grt_cmat2x2_mul(Q12, TD, RD); // RU - cmat2x2_mul(TD, Q21, RU); - cmat2x2_k(RU, -1, RU); + grt_cmat2x2_mul(TD, Q21, RU); + grt_cmat2x2_k(RU, -1, RU); // TU - cmat2x2_mul(Q12, RU, TU); - cmat2x2_add(Q11, TU, TU); + grt_cmat2x2_mul(Q12, RU, TU); + grt_cmat2x2_add(Q11, TU, TU); *RDL = (mu1*xb1 - mu2*xb2) / (mu1*xb1 + mu2*xb2) * exa*exa; *RUL = - (*RDL); diff --git a/pygrt/C_extension/src/dynamic/propagate.c b/pygrt/C_extension/src/dynamic/propagate.c index 12b065d..77a7cd0 100755 --- a/pygrt/C_extension/src/dynamic/propagate.c +++ b/pygrt/C_extension/src/dynamic/propagate.c @@ -15,17 +15,17 @@ #include #include -#include "dynamic/propagate.h" -#include "dynamic/layer.h" -#include "dynamic/source.h" -#include "common/recursion.h" -#include "common/model.h" -#include "common/matrix.h" -#include "common/prtdbg.h" +#include "grt/dynamic/propagate.h" +#include "grt/dynamic/layer.h" +#include "grt/dynamic/source.h" +#include "grt/common/recursion.h" +#include "grt/common/model.h" +#include "grt/common/matrix.h" +#include "grt/common/prtdbg.h" -void kernel( - const MODEL1D *mod1d, MYCOMPLEX omega, MYREAL k, MYCOMPLEX QWV[SRC_M_NUM][QWV_NUM], +void grt_kernel( + const GRT_MODEL1D *mod1d, MYCOMPLEX omega, MYREAL k, MYCOMPLEX QWV[SRC_M_NUM][QWV_NUM], bool calc_uiz, MYCOMPLEX QWV_uiz[SRC_M_NUM][QWV_NUM], MYINT *stats) { @@ -156,13 +156,13 @@ void kernel( src_kbkb = mod1d_kbkb1; } else { // 对第iy层的系数矩阵赋值,加入时间延迟因子(第iy-1界面与第iy界面之间) - calc_RT_PSV( + grt_calc_RT_PSV( mod1d_Rho0, mod1d_xa0, mod1d_xb0, mod1d_kbkb0, mod1d_mu0, mod1d_Rho1, mod1d_xa1, mod1d_xb1, mod1d_kbkb1, mod1d_mu1, mod1d_thk0, // 使用iy-1层的厚度 omega, k, RD, RU, TD, TU, stats); - calc_RT_SH( + grt_calc_RT_SH( mod1d_xb0, mod1d_mu0, mod1d_xb1, mod1d_mu1, mod1d_thk0, // 使用iy-1层的厚度 @@ -209,7 +209,7 @@ void kernel( GRT_RT_PSV_ASSIGN(FA); GRT_RT_SH_ASSIGN(FA); } else { // 递推FA - recursion_RT( + grt_recursion_RT( RD_FA, RDL_FA, RU_FA, RUL_FA, TD_FA, TDL_FA, TU_FA, TUL_FA, RD, RDL, RU, RUL, @@ -220,7 +220,7 @@ void kernel( } } else if(iy==imin){ // 虚拟层位,可对递推公式简化 - recursion_RT_imaginary( + grt_recursion_RT_imaginary( mod1d_xa0, mod1d_xb0, mod1d_thk0, k, RU_FA, &RUL_FA, TD_FA, &TDL_FA, TU_FA, &TUL_FA); @@ -231,7 +231,7 @@ void kernel( GRT_RT_PSV_ASSIGN(RS); GRT_RT_SH_ASSIGN(RS); } else { // 递推RS - recursion_RT( + grt_recursion_RT( RD_RS, RDL_RS, RU_RS, RUL_RS, TD_RS, TDL_RS, TU_RS, TUL_RS, RD, RDL, RU, RUL, @@ -242,7 +242,7 @@ void kernel( } } else if(iy==imax){ // 虚拟层位,可对递推公式简化 - recursion_RT_imaginary( + grt_recursion_RT_imaginary( mod1d_xa0, mod1d_xb0, mod1d_thk0, k, RU_RS, &RUL_RS, TD_RS, &TDL_RS, TU_RS, &TUL_RS); @@ -254,7 +254,7 @@ void kernel( GRT_RT_SH_ASSIGN(BL); } else { // 递推BL // 只有 RD 矩阵最终会被使用到 - recursion_RT( + grt_recursion_RT( RD_BL, RDL_BL, RU_BL, RUL_BL, TD_BL, TDL_BL, TU_BL, TUL_BL, RD, RDL, RU, RUL, @@ -275,17 +275,17 @@ void kernel( // 计算震源系数 MYCOMPLEX src_coef_PSV[SRC_M_NUM][QWV_NUM-1][2] = {0}; MYCOMPLEX src_coef_SH[SRC_M_NUM][2] = {0}; - source_coef_PSV(src_xa, src_xb, src_kaka, src_kbkb, k, src_coef_PSV); - source_coef_SH(src_xb, src_kbkb, k, src_coef_SH); + grt_source_coef_PSV(src_xa, src_xb, src_kaka, src_kbkb, k, src_coef_PSV); + grt_source_coef_SH(src_xb, src_kbkb, k, src_coef_SH); // 临时中转矩阵 (temperary) MYCOMPLEX tmpR2[2][2], tmp2x2[2][2], tmpRL, tmp2x2_uiz[2][2], tmpRL2; MYCOMPLEX inv_2x2T[2][2], invT; // 递推RU_FA - calc_R_tilt_PSV(top_xa, top_xb, top_kbkb, k, R_tilt, stats); + grt_calc_R_tilt_PSV(top_xa, top_xb, top_kbkb, k, R_tilt, stats); if(*stats==INVERSE_FAILURE) goto BEFORE_RETURN; - recursion_RU( + grt_recursion_RU( R_tilt, RONE, RD_FA, RDL_FA, RU_FA, RUL_FA, @@ -298,11 +298,11 @@ void kernel( if(ircvup){ // A接收 B震源 // 计算R_EV - calc_R_EV_PSV(rcv_xa, rcv_xb, ircvup, k, RU_FA, R_EV); - calc_R_EV_SH(rcv_xb, k, RUL_FA, &R_EVL); + grt_calc_R_EV_PSV(rcv_xa, rcv_xb, ircvup, k, RU_FA, R_EV); + grt_calc_R_EV_SH(rcv_xb, k, RUL_FA, &R_EVL); // 递推RU_FS - recursion_RU( + grt_recursion_RU( RU_FA, RUL_FA, // 已从ZR变为FR,加入了自由表面的效应 RD_RS, RDL_RS, RU_RS, RUL_RS, @@ -346,42 +346,42 @@ void kernel( #endif // 公式(5.7.12-14) - cmat2x2_mul(RD_BL, RU_FB, tmpR2); - cmat2x2_one_sub(tmpR2); - cmat2x2_inv(tmpR2, tmpR2, stats);// (I - xx)^-1 + grt_cmat2x2_mul(RD_BL, RU_FB, tmpR2); + grt_cmat2x2_one_sub(tmpR2); + grt_cmat2x2_inv(tmpR2, tmpR2, stats);// (I - xx)^-1 if(*stats==INVERSE_FAILURE) goto BEFORE_RETURN; - cmat2x2_mul(inv_2x2T, tmpR2, tmp2x2); + grt_cmat2x2_mul(inv_2x2T, tmpR2, tmp2x2); - if(calc_uiz) cmat2x2_assign(tmp2x2, tmp2x2_uiz); // 为后续计算空间导数备份 + if(calc_uiz) grt_cmat2x2_assign(tmp2x2, tmp2x2_uiz); // 为后续计算空间导数备份 - cmat2x2_mul(R_EV, tmp2x2, tmp2x2); + grt_cmat2x2_mul(R_EV, tmp2x2, tmp2x2); tmpRL = invT / (RONE - RDL_BL * RUL_FB); tmpRL2 = R_EVL * tmpRL; for(MYINT i=0; i #include -#include "dynamic/signals.h" -#include "common/const.h" -#include "common/colorstr.h" +#include "grt/dynamic/signals.h" +#include "grt/common/const.h" +#include "grt/common/colorstr.h" -bool check_tftype_tfparams(const char tftype, const char *tfparams){ +bool grt_check_tftype_tfparams(const char tftype, const char *tfparams){ // 抛物波 if(GRT_SIG_PARABOLA == tftype){ @@ -76,7 +76,7 @@ bool check_tftype_tfparams(const char tftype, const char *tfparams){ } -float * get_time_function(int *TFnt, float dt, const char tftype, const char *tfparams){ +float * grt_get_time_function(int *TFnt, float dt, const char tftype, const char *tfparams){ // 获得时间函数 float *tfarr=NULL; int tfnt=0; @@ -84,23 +84,23 @@ float * get_time_function(int *TFnt, float dt, const char tftype, const char *tf if(GRT_SIG_PARABOLA == tftype){ float t0=0.0; sscanf(tfparams, "%f", &t0); - tfarr = get_parabola_wave(dt, &t0, &tfnt); + tfarr = grt_get_parabola_wave(dt, &t0, &tfnt); } // 梯形波 else if(GRT_SIG_TRAPEZOID == tftype){ float t1=0.0, t2=0.0, t3=0.0; sscanf(tfparams, "%f/%f/%f", &t1, &t2, &t3); - tfarr = get_trap_wave(dt, &t1, &t2, &t3, &tfnt); + tfarr = grt_get_trap_wave(dt, &t1, &t2, &t3, &tfnt); } // 雷克子波 else if(GRT_SIG_RICKER == tftype){ float f0=0.0; sscanf(tfparams, "%f", &f0); - tfarr = get_ricker_wave(dt, f0, &tfnt); + tfarr = grt_get_ricker_wave(dt, f0, &tfnt); } // 自定义时间函数 else if(GRT_SIG_CUSTOM == tftype){ - tfarr = get_custom_wave(&tfnt, tfparams); + tfarr = grt_get_custom_wave(&tfnt, tfparams); } *TFnt = tfnt; @@ -113,11 +113,11 @@ void linear_convolve_time_function(float *arr, int nt, float dt, const char tfty // 获得时间函数 float *tfarr=NULL; int tfnt=0; - tfarr = get_time_function(&tfnt, dt, tftype, tfparams); + tfarr = grt_get_time_function(&tfnt, dt, tftype, tfparams); float *yarr = (float*)calloc(nt, sizeof(float)); // 线性卷积 - oaconvolve(arr, nt, tfarr, tfnt, yarr, nt, false); + grt_oaconvolve(arr, nt, tfarr, tfnt, yarr, nt, false); // 原地更改 for(int i=0; i #include -#include "dynamic/source.h" -#include "common/model.h" -#include "common/matrix.h" -#include "common/prtdbg.h" +#include "grt/dynamic/source.h" +#include "grt/common/model.h" +#include "grt/common/matrix.h" +#include "grt/common/prtdbg.h" -void source_coef_PSV( +void grt_source_coef_PSV( MYCOMPLEX src_xa, MYCOMPLEX src_xb, MYCOMPLEX src_kaka, MYCOMPLEX src_kbkb, MYREAL k, MYCOMPLEX coef[SRC_M_NUM][QWV_NUM-1][2]) @@ -66,7 +66,7 @@ void source_coef_PSV( } -void source_coef_SH( +void grt_source_coef_SH( MYCOMPLEX src_xb, MYCOMPLEX src_kbkb, MYREAL k, MYCOMPLEX coef[SRC_M_NUM][2]) diff --git a/pygrt/C_extension/src/grt.c b/pygrt/C_extension/src/grt.c index 3ce9db0..397dc23 100644 --- a/pygrt/C_extension/src/grt.c +++ b/pygrt/C_extension/src/grt.c @@ -42,7 +42,7 @@ static void free_Ctrl(GRT_MAIN_CTRL *Ctrl){ /** 打印使用说明 */ static void print_help(){ -print_logo(); +grt_print_logo(); printf("\n" "Usage: \n" "----------------------------------------------------------------\n" diff --git a/pygrt/C_extension/src/static/grt_static_greenfn.c b/pygrt/C_extension/src/static/grt_static_greenfn.c index 54ae7ad..c58edf2 100644 --- a/pygrt/C_extension/src/static/grt_static_greenfn.c +++ b/pygrt/C_extension/src/static/grt_static_greenfn.c @@ -8,13 +8,13 @@ */ -#include "static/static_grn.h" -#include "common/const.h" -#include "common/model.h" -#include "common/integral.h" -#include "common/iostats.h" -#include "common/search.h" -#include "common/util.h" +#include "grt/static/static_grn.h" +#include "grt/common/const.h" +#include "grt/common/model.h" +#include "grt/common/integral.h" +#include "grt/common/iostats.h" +#include "grt/common/search.h" +#include "grt/common/util.h" #include "grt.h" @@ -32,7 +32,7 @@ typedef struct { bool active; char *s_modelpath; ///< 模型路径 const char *s_modelname; ///< 模型名称 - PYMODEL1D *pymod; ///< 模型PYMODEL1D结构体指针 + GRT_PYMODEL1D *pymod; ///< 模型PYMODEL1D结构体指针 } M; /** 震源和接收器深度 */ struct { @@ -94,7 +94,7 @@ static void free_Ctrl(GRT_MODULE_CTRL *Ctrl){ // M GRT_SAFE_FREE_PTR(Ctrl->M.s_modelpath); - free_pymod(Ctrl->M.pymod); + grt_free_pymod(Ctrl->M.pymod); // D GRT_SAFE_FREE_PTR(Ctrl->D.s_depsrc); @@ -244,7 +244,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ case 'M': Ctrl->M.active = true; Ctrl->M.s_modelpath = strdup(optarg); - Ctrl->M.s_modelname = get_basename(Ctrl->M.s_modelpath); + Ctrl->M.s_modelname = grt_get_basename(Ctrl->M.s_modelpath); break; // 震源和场点深度, -Ddepsrc/deprcv @@ -443,14 +443,14 @@ int static_greenfn_main(int argc, char **argv){ getopt_from_command(Ctrl, argc, argv); // 读入模型文件(暂先不考虑液体层) - if((Ctrl->M.pymod = read_pymod_from_file(command, Ctrl->M.s_modelpath, Ctrl->D.depsrc, Ctrl->D.deprcv, false)) == NULL){ + if((Ctrl->M.pymod = grt_read_pymod_from_file(command, Ctrl->M.s_modelpath, Ctrl->D.depsrc, Ctrl->D.deprcv, false)) == NULL){ exit(EXIT_FAILURE); } - PYMODEL1D *pymod = Ctrl->M.pymod; + GRT_PYMODEL1D *pymod = Ctrl->M.pymod; // 最大最小速度 MYREAL vmin, vmax; - get_pymod_vmin_vmax(pymod, &vmin, &vmax); + grt_get_pymod_vmin_vmax(pymod, &vmin, &vmax); // 参考最小速度 if(!Ctrl->V.active){ @@ -483,7 +483,7 @@ int static_greenfn_main(int argc, char **argv){ //============================================================================== // 计算静态格林函数 - integ_static_grn( + grt_integ_static_grn( pymod, Ctrl->nr, Ctrl->rs, Ctrl->V.vmin_ref, Ctrl->K.keps, Ctrl->K.k0, Ctrl->L.Length, Ctrl->L.filonLength, Ctrl->L.safilonTol, Ctrl->L.filonCut, grn, Ctrl->e.active, grn_uiz, grn_uir, Ctrl->S.s_statsdir diff --git a/pygrt/C_extension/src/static/grt_static_rotation.c b/pygrt/C_extension/src/static/grt_static_rotation.c index 20bf647..2587c6a 100644 --- a/pygrt/C_extension/src/static/grt_static_rotation.c +++ b/pygrt/C_extension/src/static/grt_static_rotation.c @@ -8,7 +8,7 @@ */ -#include "common/const.h" +#include "grt/common/const.h" #include "grt.h" diff --git a/pygrt/C_extension/src/static/grt_static_strain.c b/pygrt/C_extension/src/static/grt_static_strain.c index f227fe5..6d16ffe 100644 --- a/pygrt/C_extension/src/static/grt_static_strain.c +++ b/pygrt/C_extension/src/static/grt_static_strain.c @@ -7,7 +7,7 @@ * */ -#include "common/const.h" +#include "grt/common/const.h" #include "grt.h" diff --git a/pygrt/C_extension/src/static/grt_static_stress.c b/pygrt/C_extension/src/static/grt_static_stress.c index da644d7..8425652 100644 --- a/pygrt/C_extension/src/static/grt_static_stress.c +++ b/pygrt/C_extension/src/static/grt_static_stress.c @@ -7,7 +7,7 @@ * */ -#include "common/const.h" +#include "grt/common/const.h" #include "grt.h" diff --git a/pygrt/C_extension/src/static/grt_static_syn.c b/pygrt/C_extension/src/static/grt_static_syn.c index e22bf2b..840e714 100644 --- a/pygrt/C_extension/src/static/grt_static_syn.c +++ b/pygrt/C_extension/src/static/grt_static_syn.c @@ -8,9 +8,9 @@ */ -#include "common/const.h" -#include "common/radiation.h" -#include "common/coord.h" +#include "grt/common/const.h" +#include "grt/common/radiation.h" +#include "grt/common/coord.h" #include "grt.h" @@ -454,7 +454,7 @@ int static_syn_main(int argc, char **argv){ tmpsyn[0] = tmpsyn[1] = tmpsyn[2] = 0.0; // 计算震源辐射因子 - set_source_radiation(Ctrl->srcRadi, Ctrl->computeType, ityp==3, Ctrl->S.M0, upar_scale, azrad, Ctrl->mchn); + grt_set_source_radiation(Ctrl->srcRadi, Ctrl->computeType, ityp==3, Ctrl->S.M0, upar_scale, azrad, Ctrl->mchn); for(int i=0; ie.active){ - rot_zrt2zxy_upar(azrad, syn, syn_upar, dist*1e5); + grt_rot_zrt2zxy_upar(azrad, syn, syn_upar, dist*1e5); } else { - rot_zxy2zrt_vec(-azrad, syn); + grt_rot_zxy2zrt_vec(-azrad, syn); } } diff --git a/pygrt/C_extension/src/static/static_grn.c b/pygrt/C_extension/src/static/static_grn.c index e1853d3..80894ab 100644 --- a/pygrt/C_extension/src/static/static_grn.c +++ b/pygrt/C_extension/src/static/static_grn.c @@ -19,16 +19,16 @@ #include #include -#include "static/static_grn.h" -#include "static/static_propagate.h" -#include "common/dwm.h" -#include "common/ptam.h" -#include "common/fim.h" -#include "common/safim.h" -#include "common/const.h" -#include "common/model.h" -#include "common/integral.h" -#include "common/search.h" +#include "grt/static/static_grn.h" +#include "grt/static/static_propagate.h" +#include "grt/common/dwm.h" +#include "grt/common/ptam.h" +#include "grt/common/fim.h" +#include "grt/common/safim.h" +#include "grt/common/const.h" +#include "grt/common/model.h" +#include "grt/common/integral.h" +#include "grt/common/search.h" @@ -48,7 +48,7 @@ static void recordin_GRN( MYCOMPLEX (*tmp_grn)[SRC_M_NUM][CHANNEL_NUM] = (MYCOMPLEX(*)[SRC_M_NUM][CHANNEL_NUM])calloc(nr, sizeof(*tmp_grn)); for(MYINT ir=0; ir mod1d - MODEL1D *mod1d = init_mod1d(pymod1d->n); - get_mod1d(pymod1d, mod1d); + GRT_MODEL1D *mod1d = grt_init_mod1d(pymod1d->n); + grt_get_mod1d(pymod1d, mod1d); const MYREAL hs = (fabs(pymod1d->depsrc - pymod1d->deprcv) < MIN_DEPTH_GAP_SRC_RCV)? MIN_DEPTH_GAP_SRC_RCV : fabs(pymod1d->depsrc - pymod1d->deprcv); // hs=max(震源和台站深度差,1.0) @@ -159,35 +159,35 @@ void integ_static_grn( MYINT inv_stats=INVERSE_SUCCESS; // 常规的波数积分 - k = discrete_integ( + k = grt_discrete_integ( mod1d, dk, (useFIM)? filonK : kmax, keps, 0.0, nr, rs, sum_J, calc_upar, sum_uiz_J, sum_uir_J, - fstats, static_kernel, &inv_stats); + fstats, grt_static_kernel, &inv_stats); // 基于线性插值的Filon积分 if(useFIM){ if(filondk > RZERO){ // 基于线性插值的Filon积分,固定采样间隔 - k = linear_filon_integ( + k = grt_linear_filon_integ( mod1d, k, dk, filondk, kmax, keps, 0.0, nr, rs, sum_J, calc_upar, sum_uiz_J, sum_uir_J, - fstats, static_kernel, &inv_stats); + fstats, grt_static_kernel, &inv_stats); } else if(safilonTol > RZERO){ // 基于自适应采样的Filon积分 - k = sa_filon_integ( + k = grt_sa_filon_integ( mod1d, kmax, k, dk, safilonTol, kmax, 0.0, nr, rs, sum_J, calc_upar, sum_uiz_J, sum_uir_J, - fstats, static_kernel, &inv_stats); + fstats, grt_static_kernel, &inv_stats); } } // k之后的部分使用峰谷平均法进行显式收敛,建议在浅源地震的时候使用 if(vmin_ref < RZERO){ - PTA_method( + grt_PTA_method( mod1d, k, dk, 0.0, nr, rs, sum_J, calc_upar, sum_uiz_J, sum_uir_J, - ptam_fstatsnr, static_kernel, &inv_stats); + ptam_fstatsnr, grt_static_kernel, &inv_stats); } @@ -208,7 +208,7 @@ void integ_static_grn( GRT_SAFE_FREE_PTR(sum_uiz_J); GRT_SAFE_FREE_PTR(sum_uir_J); - free_mod1d(mod1d); + grt_free_mod1d(mod1d); GRT_SAFE_FREE_PTR_ARRAY(ptam_fstatsdir, nr); diff --git a/pygrt/C_extension/src/static/static_layer.c b/pygrt/C_extension/src/static/static_layer.c index a226cd3..2091f76 100644 --- a/pygrt/C_extension/src/static/static_layer.c +++ b/pygrt/C_extension/src/static/static_layer.c @@ -15,38 +15,38 @@ #include #include -#include "static/static_layer.h" -#include "common/model.h" -#include "common/matrix.h" +#include "grt/static/static_layer.h" +#include "grt/common/model.h" +#include "grt/common/matrix.h" -void calc_static_R_tilt_PSV(MYCOMPLEX delta1, MYCOMPLEX R_tilt[2][2]){ +void grt_calc_static_R_tilt_PSV(MYCOMPLEX delta1, MYCOMPLEX R_tilt[2][2]){ // 公式(6.3.12) R_tilt[0][0] = R_tilt[1][1] = CZERO; R_tilt[0][1] = -delta1; R_tilt[1][0] = -RONE/delta1; } -void calc_static_R_EV_PSV(bool ircvup, const MYCOMPLEX R[2][2], MYCOMPLEX R_EV[2][2]) +void grt_calc_static_R_EV_PSV(bool ircvup, const MYCOMPLEX R[2][2], MYCOMPLEX R_EV[2][2]) { MYCOMPLEX D11[2][2] = {{RONE, -RONE}, {RONE, RONE}}; MYCOMPLEX D12[2][2] = {{RONE, -RONE}, {-RONE, -RONE}}; // 公式(6.3.35,37) if(ircvup){// 震源更深 - cmat2x2_mul(D12, R, R_EV); - cmat2x2_add(D11, R_EV, R_EV); + grt_cmat2x2_mul(D12, R, R_EV); + grt_cmat2x2_add(D11, R_EV, R_EV); } else { // 接收点更深 - cmat2x2_mul(D11, R, R_EV); - cmat2x2_add(D12, R_EV, R_EV); + grt_cmat2x2_mul(D11, R, R_EV); + grt_cmat2x2_add(D12, R_EV, R_EV); } } -void calc_static_R_EV_SH(MYCOMPLEX RL, MYCOMPLEX *R_EVL) +void grt_calc_static_R_EV_SH(MYCOMPLEX RL, MYCOMPLEX *R_EVL) { *R_EVL = (RONE + (RL)); } -void calc_static_uiz_R_EV_PSV( +void grt_calc_static_uiz_R_EV_PSV( MYCOMPLEX delta1, bool ircvup, MYREAL k, const MYCOMPLEX R[2][2], MYCOMPLEX R_EV[2][2]) { @@ -55,15 +55,15 @@ void calc_static_uiz_R_EV_PSV( MYCOMPLEX D11[2][2] = {{k, -k-kd2}, {k, k-kd2}}; MYCOMPLEX D12[2][2] = {{-k, k+kd2}, {k, k-kd2}}; if(ircvup){// 震源更深 - cmat2x2_mul(D12, R, R_EV); - cmat2x2_add(D11, R_EV, R_EV); + grt_cmat2x2_mul(D12, R, R_EV); + grt_cmat2x2_add(D11, R_EV, R_EV); } else { // 接收点更深 - cmat2x2_mul(D11, R, R_EV); - cmat2x2_add(D12, R_EV, R_EV); + grt_cmat2x2_mul(D11, R, R_EV); + grt_cmat2x2_add(D12, R_EV, R_EV); } } -void calc_static_uiz_R_EV_SH(bool ircvup, MYREAL k, MYCOMPLEX RL, MYCOMPLEX *R_EVL) +void grt_calc_static_uiz_R_EV_SH(bool ircvup, MYREAL k, MYCOMPLEX RL, MYCOMPLEX *R_EVL) { // 新推导公式 if(ircvup){// 震源更深 @@ -74,7 +74,7 @@ void calc_static_uiz_R_EV_SH(bool ircvup, MYREAL k, MYCOMPLEX RL, MYCOMPLEX *R_E } -void calc_static_RT_PSV( +void grt_calc_static_RT_PSV( MYCOMPLEX delta1, MYCOMPLEX mu1, MYCOMPLEX delta2, MYCOMPLEX mu2, MYREAL thk, MYREAL k, @@ -128,7 +128,7 @@ void calc_static_RT_PSV( // printf("-----------------------------\n"); } -void calc_static_RT_SH( +void grt_calc_static_RT_SH( MYCOMPLEX mu1, MYCOMPLEX mu2, MYREAL thk, MYREAL k, MYCOMPLEX *RDL, MYCOMPLEX *RUL, MYCOMPLEX *TDL, MYCOMPLEX *TUL) diff --git a/pygrt/C_extension/src/static/static_propagate.c b/pygrt/C_extension/src/static/static_propagate.c index b56724e..4e8c06d 100644 --- a/pygrt/C_extension/src/static/static_propagate.c +++ b/pygrt/C_extension/src/static/static_propagate.c @@ -15,17 +15,17 @@ #include #include -#include "static/static_propagate.h" -#include "static/static_layer.h" -#include "static/static_source.h" -#include "common/recursion.h" -#include "common/model.h" -#include "common/const.h" -#include "common/matrix.h" +#include "grt/static/static_propagate.h" +#include "grt/static/static_layer.h" +#include "grt/static/static_source.h" +#include "grt/common/recursion.h" +#include "grt/common/model.h" +#include "grt/common/const.h" +#include "grt/common/matrix.h" -void static_kernel( - const MODEL1D *mod1d, MYCOMPLEX omega, MYREAL k, MYCOMPLEX QWV[SRC_M_NUM][QWV_NUM], +void grt_static_kernel( + const GRT_MODEL1D *mod1d, MYCOMPLEX omega, MYREAL k, MYCOMPLEX QWV[SRC_M_NUM][QWV_NUM], bool calc_uiz, MYCOMPLEX QWV_uiz[SRC_M_NUM][QWV_NUM], MYINT *stats) { // 初始化qwv为0 @@ -130,12 +130,12 @@ void static_kernel( } // 对第iy层的系数矩阵赋值,加入时间延迟因子(第iy-1界面与第iy界面之间) - calc_static_RT_PSV( + grt_calc_static_RT_PSV( mod1d_delta0, mod1d_mu0, mod1d_delta1, mod1d_mu1, mod1d_thk0, k, // 使用iy-1层的厚度 RD, RU, TD, TU); - calc_static_RT_SH( + grt_calc_static_RT_SH( mod1d_mu0, mod1d_mu1, mod1d_thk0, k, // 使用iy-1层的厚度 &RDL, &RUL, &TDL, &TUL); @@ -146,7 +146,7 @@ void static_kernel( GRT_RT_PSV_ASSIGN(FA); GRT_RT_SH_ASSIGN(FA); } else { // 递推FA - recursion_RT( + grt_recursion_RT( RD_FA, RDL_FA, RU_FA, RUL_FA, TD_FA, TDL_FA, TU_FA, TUL_FA, RD, RDL, RU, RUL, @@ -161,7 +161,7 @@ void static_kernel( GRT_RT_PSV_ASSIGN(RS); GRT_RT_SH_ASSIGN(RS); } else { // 递推RS - recursion_RT( + grt_recursion_RT( RD_RS, RDL_RS, RU_RS, RUL_RS, TD_RS, TDL_RS, TU_RS, TUL_RS, RD, RDL, RU, RUL, @@ -177,7 +177,7 @@ void static_kernel( GRT_RT_SH_ASSIGN(BL); } else { // 递推BL // 只有 RD 矩阵最终会被使用到 - recursion_RT( + grt_recursion_RT( RD_BL, RDL_BL, RU_BL, RUL_BL, TD_BL, TDL_BL, TU_BL, TUL_BL, RD, RDL, RU, RUL, @@ -195,16 +195,16 @@ void static_kernel( // 计算震源系数 MYCOMPLEX src_coef_PSV[SRC_M_NUM][QWV_NUM-1][2] = {0}; MYCOMPLEX src_coef_SH[SRC_M_NUM][2] = {0}; - static_source_coef_PSV(src_delta, k, src_coef_PSV); - static_source_coef_SH(k, src_coef_SH); + grt_static_source_coef_PSV(src_delta, k, src_coef_PSV); + grt_static_source_coef_SH(k, src_coef_SH); // 临时中转矩阵 (temperary) MYCOMPLEX tmpR2[2][2], tmp2x2[2][2], tmpRL, tmp2x2_uiz[2][2], tmpRL_uiz; MYCOMPLEX inv_2x2T[2][2], invT; // 递推RU_FA - calc_static_R_tilt_PSV(top_delta, R_tilt); - recursion_RU( + grt_calc_static_R_tilt_PSV(top_delta, R_tilt); + grt_recursion_RU( R_tilt, RONE, RD_FA, RDL_FA, RU_FA, RUL_FA, @@ -215,11 +215,11 @@ void static_kernel( // 根据震源和台站相对位置,计算最终的系数 if(ircvup){ // A接收 B震源 // 计算R_EV - calc_static_R_EV_PSV(ircvup, RU_FA, R_EV); - calc_static_R_EV_SH(RUL_FA, &R_EVL); + grt_calc_static_R_EV_PSV(ircvup, RU_FA, R_EV); + grt_calc_static_R_EV_SH(RUL_FA, &R_EVL); // 递推RU_FS - recursion_RU( + grt_recursion_RU( RU_FA, RUL_FA, // 已从ZR变为FR,加入了自由表面的效应 RD_RS, RDL_RS, RU_RS, RUL_RS, @@ -228,40 +228,40 @@ void static_kernel( RU_FB, &RUL_FB, inv_2x2T, &invT, stats); // 公式(5.7.12-14) - cmat2x2_mul(RD_BL, RU_FB, tmpR2); - cmat2x2_one_sub(tmpR2); - cmat2x2_inv(tmpR2, tmpR2, stats);// (I - xx)^-1 - cmat2x2_mul(inv_2x2T, tmpR2, tmp2x2); + grt_cmat2x2_mul(RD_BL, RU_FB, tmpR2); + grt_cmat2x2_one_sub(tmpR2); + grt_cmat2x2_inv(tmpR2, tmpR2, stats);// (I - xx)^-1 + grt_cmat2x2_mul(inv_2x2T, tmpR2, tmp2x2); - if(calc_uiz) cmat2x2_assign(tmp2x2, tmp2x2_uiz); // 为后续计算空间导数备份 + if(calc_uiz) grt_cmat2x2_assign(tmp2x2, tmp2x2_uiz); // 为后续计算空间导数备份 - cmat2x2_mul(R_EV, tmp2x2, tmp2x2); + grt_cmat2x2_mul(R_EV, tmp2x2, tmp2x2); tmpRL = R_EVL * invT / (RONE - RDL_BL * RUL_FB); for(MYINT i=0; i #include -#include "static/static_source.h" -#include "common/const.h" +#include "grt/static/static_source.h" +#include "grt/common/const.h" -void static_source_coef_PSV(MYCOMPLEX delta, MYREAL k, MYCOMPLEX coef[SRC_M_NUM][QWV_NUM-1][2]) +void grt_static_source_coef_PSV(MYCOMPLEX delta, MYREAL k, MYCOMPLEX coef[SRC_M_NUM][QWV_NUM-1][2]) { // 先全部赋0 for(MYINT i=0; iR.active = true; - Ctrl->R.s_rs = string_split(optarg, ",", &Ctrl->R.nr); + Ctrl->R.s_rs = grt_string_split(optarg, ",", &Ctrl->R.nr); // 转为浮点数 Ctrl->R.rs = (MYREAL*)realloc(Ctrl->R.rs, sizeof(MYREAL)*(Ctrl->R.nr)); for(MYINT i=0; iR.nr; ++i){ @@ -521,10 +521,10 @@ int travt_main(int argc, char **argv){ getopt_from_command(Ctrl, argc, argv); - PYMODEL1D *pymod; + GRT_PYMODEL1D *pymod; // 读入模型文件 - if((pymod = read_pymod_from_file(Ctrl->name, Ctrl->M.s_modelpath, Ctrl->D.depsrc, Ctrl->D.deprcv, true)) == NULL){ + if((pymod = grt_read_pymod_from_file(Ctrl->name, Ctrl->M.s_modelpath, Ctrl->D.depsrc, Ctrl->D.deprcv, true)) == NULL){ exit(EXIT_FAILURE); } // print_pymod(pymod); @@ -533,9 +533,9 @@ int travt_main(int argc, char **argv){ printf(" Distance(km) Tp(secs) Ts(secs) \n"); double travtP=-1, travtS=-1; for(int i=0; iR.nr; ++i){ - travtP = compute_travt1d( + travtP = grt_compute_travt1d( pymod->Thk, pymod->Va, pymod->n, pymod->isrc, pymod->ircv, Ctrl->R.rs[i]); - travtS = compute_travt1d( + travtS = grt_compute_travt1d( pymod->Thk, pymod->Vb, pymod->n, pymod->isrc, pymod->ircv, Ctrl->R.rs[i]); printf(" %-15s %-15.3f %-15.3f\n", Ctrl->R.s_rs[i], travtP, travtS); diff --git a/pygrt/c_interfaces.py b/pygrt/c_interfaces.py index 9391746..e773f99 100755 --- a/pygrt/c_interfaces.py +++ b/pygrt/c_interfaces.py @@ -27,9 +27,9 @@ """libgrt库""" -C_integ_grn_spec = libgrt.integ_grn_spec +C_grt_integ_grn_spec = libgrt.grt_integ_grn_spec """C库中计算格林函数的主函数 integ_grn_spec, 详见C API同名函数""" -C_integ_grn_spec.argtypes = [ +C_grt_integ_grn_spec.argtypes = [ POINTER(c_PyModel1D), c_int, c_int, PREAL, c_int, PREAL, REAL, REAL, REAL, REAL, REAL, REAL, REAL, REAL, REAL, @@ -49,10 +49,10 @@ ] -C_integ_static_grn = libgrt.integ_static_grn +C_grt_integ_static_grn = libgrt.grt_integ_static_grn """计算静态格林函数""" -C_integ_static_grn.restype = None -C_integ_static_grn.argtypes = [ +C_grt_integ_static_grn.restype = None +C_grt_integ_static_grn.argtypes = [ POINTER(c_PyModel1D), c_int, PREAL, REAL, REAL, REAL, REAL, REAL, REAL, REAL, POINTER((REAL*CHANNEL_NUM)*SRC_M_NUM), @@ -63,10 +63,10 @@ ] -C_set_num_threads = libgrt.set_num_threads +C_grt_set_num_threads = libgrt.grt_set_num_threads """设置多线程数""" -C_set_num_threads.restype = None -C_set_num_threads.argtypes = [c_int] +C_grt_set_num_threads.restype = None +C_grt_set_num_threads.argtypes = [c_int] def set_num_threads(n): @@ -75,76 +75,76 @@ def set_num_threads(n): :param n: 线程数 ''' - C_set_num_threads(n) + C_grt_set_num_threads(n) -C_compute_travt1d = libgrt.compute_travt1d +C_grt_compute_travt1d = libgrt.grt_compute_travt1d """计算1D层状半空间的初至波走时""" -C_compute_travt1d.restype = REAL -C_compute_travt1d.argtypes = [ +C_grt_compute_travt1d.restype = REAL +C_grt_compute_travt1d.argtypes = [ PREAL, PREAL, c_int, c_int, c_int, REAL ] -C_read_pymod_from_file = libgrt.read_pymod_from_file +C_grt_read_pymod_from_file = libgrt.grt_read_pymod_from_file """读取模型文件并进行预处理""" -C_read_pymod_from_file.restype = POINTER(c_PyModel1D) -C_read_pymod_from_file.argtypes = [c_char_p, c_char_p, c_double, c_double, c_bool] +C_grt_read_pymod_from_file.restype = POINTER(c_PyModel1D) +C_grt_read_pymod_from_file.argtypes = [c_char_p, c_char_p, c_double, c_double, c_bool] -C_free_pymod = libgrt.free_pymod +C_grt_free_pymod = libgrt.grt_free_pymod """释放C程序中申请的PYMODEL1D结构体内存""" -C_free_pymod.restype = None -C_free_pymod.argtypes = [POINTER(c_PyModel1D)] +C_grt_free_pymod.restype = None +C_grt_free_pymod.argtypes = [POINTER(c_PyModel1D)] # ------------------------------------------------------------------- # C函数定义的时间函数 # ------------------------------------------------------------------- -C_free = libgrt.free1d +C_grt_free = libgrt.grt_free1d """释放在C中申请的内存""" -C_free.restype = None -C_free.argtypes = [c_void_p] +C_grt_free.restype = None +C_grt_free.argtypes = [c_void_p] -C_get_trap_wave = libgrt.get_trap_wave +C_grt_get_trap_wave = libgrt.grt_get_trap_wave """梯形波""" -C_get_trap_wave.restype = FPOINTER -C_get_trap_wave.argtypes = [c_float, FPOINTER, FPOINTER, FPOINTER, IPOINTER] +C_grt_get_trap_wave.restype = FPOINTER +C_grt_get_trap_wave.argtypes = [c_float, FPOINTER, FPOINTER, FPOINTER, IPOINTER] -C_get_parabola_wave = libgrt.get_parabola_wave +C_grt_get_parabola_wave = libgrt.grt_get_parabola_wave """抛物波""" -C_get_parabola_wave.restype = FPOINTER -C_get_parabola_wave.argtypes = [c_float, FPOINTER, IPOINTER] +C_grt_get_parabola_wave.restype = FPOINTER +C_grt_get_parabola_wave.argtypes = [c_float, FPOINTER, IPOINTER] -C_get_ricker_wave = libgrt.get_ricker_wave +C_grt_get_ricker_wave = libgrt.grt_get_ricker_wave """雷克子波""" -C_get_ricker_wave.restype = FPOINTER -C_get_ricker_wave.argtypes = [c_float, c_float, IPOINTER] +C_grt_get_ricker_wave.restype = FPOINTER +C_grt_get_ricker_wave.argtypes = [c_float, c_float, IPOINTER] # ------------------------------------------------------------------- # C函数定义的旋转函数 # ------------------------------------------------------------------- -C_rot_zxy2zrt_vec = libgrt.rot_zxy2zrt_vec +C_grt_rot_zxy2zrt_vec = libgrt.grt_rot_zxy2zrt_vec """直角坐标zxy到柱坐标zrt的矢量旋转""" -C_rot_zxy2zrt_vec.restype = None -C_rot_zxy2zrt_vec.argtypes = [c_double, DPOINTER] # double, double[3] +C_grt_rot_zxy2zrt_vec.restype = None +C_grt_rot_zxy2zrt_vec.argtypes = [c_double, DPOINTER] # double, double[3] -C_rot_zxy2zrt_symtensor2odr = libgrt.rot_zxy2zrt_symtensor2odr +C_grt_rot_zxy2zrt_symtensor2odr = libgrt.grt_rot_zxy2zrt_symtensor2odr """直角坐标zxy到柱坐标zrt的二阶对称张量旋转""" -C_rot_zxy2zrt_symtensor2odr.restype = None -C_rot_zxy2zrt_symtensor2odr.argtypes = [c_double, DPOINTER] # double, double[6] +C_grt_rot_zxy2zrt_symtensor2odr.restype = None +C_grt_rot_zxy2zrt_symtensor2odr.argtypes = [c_double, DPOINTER] # double, double[6] -C_rot_zrt2zxy_upar = libgrt.rot_zrt2zxy_upar +C_grt_rot_zrt2zxy_upar = libgrt.grt_rot_zrt2zxy_upar """柱坐标下的位移偏导 ∂u(z,r,t)/∂(z,r,t) 转到 直角坐标 ∂u(z,x,y)/∂(z,x,y)""" -C_rot_zrt2zxy_upar.restype = None -C_rot_zrt2zxy_upar.argtypes = [c_double, DPOINTER, DPOINTER, c_double] # double, double[3], double[3][3], double +C_grt_rot_zrt2zxy_upar.restype = None +C_grt_rot_zrt2zxy_upar.argtypes = [c_double, DPOINTER, DPOINTER, c_double] # double, double[3], double[3][3], double # ------------------------------------------------------------------- # C函数定义的衰减函数 # ------------------------------------------------------------------- -C_py_attenuation_law = libgrt.py_attenuation_law +C_grt_py_attenuation_law = libgrt.grt_py_attenuation_law """品质因子Q 对 波速的影响""" -C_py_attenuation_law.restype = None -C_py_attenuation_law.argtypes = [REAL, DPOINTER, DPOINTER] # double, double[2], double[2] +C_grt_py_attenuation_law.restype = None +C_grt_py_attenuation_law.argtypes = [REAL, DPOINTER, DPOINTER] # double, double[2], double[2] diff --git a/pygrt/pymod.py b/pygrt/pymod.py index f9efd19..8861a57 100755 --- a/pygrt/pymod.py +++ b/pygrt/pymod.py @@ -54,7 +54,7 @@ def __init__(self, modarr0:np.ndarray, depsrc:float, deprcv:float): tmp_path = tmpfile.name # 获取临时文件路径 try: - c_pymod_ptr = C_read_pymod_from_file("pygrt".encode("utf-8"), tmp_path.encode("utf-8"), depsrc, deprcv, True) + c_pymod_ptr = C_grt_read_pymod_from_file("pygrt".encode("utf-8"), tmp_path.encode("utf-8"), depsrc, deprcv, True) self.c_pymod1d = c_pymod_ptr.contents # 这部分内存在C中申请,需由C函数释放。占用不多,这里跳过 finally: if os.path.exists(tmp_path): @@ -85,7 +85,7 @@ def compute_travt1d(self, dist:float): - **travtP** - 初至P波走时(s) - **travtS** - 初至S波走时(s) """ - travtP = C_compute_travt1d( + travtP = C_grt_compute_travt1d( self.c_pymod1d.Thk, self.c_pymod1d.Va, self.c_pymod1d.n, @@ -93,7 +93,7 @@ def compute_travt1d(self, dist:float): self.c_pymod1d.ircv, dist ) - travtS = C_compute_travt1d( + travtS = C_grt_compute_travt1d( self.c_pymod1d.Thk, self.c_pymod1d.Vb, self.c_pymod1d.n, @@ -374,7 +374,7 @@ def compute_grn( # 爆炸源 EX[ZR] 1e-20 cm/(dyne*cm) # 剪切源 DD[ZR],DS[ZRT],SS[ZRT] 1e-20 cm/(dyne*cm) #================================================================================= - C_integ_grn_spec( + C_grt_integ_grn_spec( self.c_pymod1d, nf1, nf2, c_freqs, nrs, c_rs, wI, vmin_ref, keps, ampk, k0, Length, filonLength, safilonTol, filonCut, print_runtime, c_grnArr, calc_upar, c_grnArr_uiz, c_grnArr_uir, @@ -555,7 +555,7 @@ def compute_static_grn( # 爆炸源 EX[ZR] 1e-20 cm/(dyne*cm) # 剪切源 DD[ZR],DS[ZRT],SS[ZRT] 1e-20 cm/(dyne*cm) #================================================================================= - C_integ_static_grn( + C_grt_integ_static_grn( self.c_pymod1d, nr, c_rs, vmin_ref, keps, k0, Length, filonLength, safilonTol, filonCut, c_pygrn, calc_upar, c_pygrn_uiz, c_pygrn_uir, c_statsfile diff --git a/pygrt/signals.py b/pygrt/signals.py index 16abef9..25fa968 100755 --- a/pygrt/signals.py +++ b/pygrt/signals.py @@ -48,10 +48,10 @@ def gen_parabola_wave(vlen, dt): ct1 = c_float(vlen) cnt = c_int(0) - carr = C_get_parabola_wave(dt, byref(ct1), byref(cnt)) + carr = C_grt_get_parabola_wave(dt, byref(ct1), byref(cnt)) arr = npct.as_array(carr, shape=(cnt.value,)).copy() - C_free(carr) + C_grt_free(carr) return arr @@ -72,10 +72,10 @@ def gen_trap_wave(t1, t2, t3, dt): ct3 = c_float(t3) cnt = c_int(0) - carr = C_get_trap_wave(dt, byref(ct1), byref(ct2), byref(ct3), byref(cnt)) + carr = C_grt_get_trap_wave(dt, byref(ct1), byref(ct2), byref(ct3), byref(cnt)) arr = npct.as_array(carr, shape=(cnt.value,)).copy() - C_free(carr) + C_grt_free(carr) return arr @@ -92,11 +92,11 @@ def gen_ricker_wave(f0:float, dt:float): ''' cnt = c_int(0) - carr = C_get_ricker_wave(dt, f0, byref(cnt)) + carr = C_grt_get_ricker_wave(dt, f0, byref(cnt)) if cast(carr, c_void_p).value is None: raise ValueError("NULL pointer") arr = npct.as_array(carr, shape=(cnt.value,)).copy() - C_free(carr) + C_grt_free(carr) return arr \ No newline at end of file diff --git a/pygrt/utils.py b/pygrt/utils.py index b10dcab..efabdbd 100755 --- a/pygrt/utils.py +++ b/pygrt/utils.py @@ -215,9 +215,9 @@ def _gen_syn_from_static_gf(grn:dict, calc_upar:bool, compute_type:str, M0:float for k in range(3): dblupar[k + i*3] = XX[i+1, k, ix, iy] if calc_upar: - C_rot_zrt2zxy_upar(azrad, dblsyn, dblupar, dist*1e5) + C_grt_rot_zrt2zxy_upar(azrad, dblsyn, dblupar, dist*1e5) else: - C_rot_zxy2zrt_vec(-azrad, dblsyn) + C_grt_rot_zxy2zrt_vec(-azrad, dblsyn) for i in range(3): XX[0, i, ix, iy] = dblsyn[i] @@ -291,9 +291,9 @@ def _data_zrt2zne(stall:Stream): dbleupar[i2 + i1*3] = uparLst[i2 + i1*3].data[n] if doupar: - C_rot_zrt2zxy_upar(azrad, dblsyn, dbleupar, dist*1e5) + C_grt_rot_zrt2zxy_upar(azrad, dblsyn, dbleupar, dist*1e5) else: - C_rot_zxy2zrt_vec(-azrad, dblsyn) + C_grt_rot_zxy2zrt_vec(-azrad, dblsyn) # 将结果写入原数组 for i1 in range(3): @@ -804,10 +804,10 @@ def _compute_stress(st_syn:Stream): freq = 0.01 if i==0 else df*i w = 2.0*np.pi*freq omega[0] = w - C_py_attenuation_law(Qbinv, omega, atte) + C_grt_py_attenuation_law(Qbinv, omega, atte) attb = atte[0] + atte[1]*1j mus[i] = vb*vb*attb*attb*rho*1e10 - C_py_attenuation_law(Qainv, omega, atte) + C_grt_py_attenuation_law(Qainv, omega, atte) atta = atte[0] + atte[1]*1j lams[i] = va*va*atta*atta*rho*1e10 - 2.0*mus[i]