diff --git a/pygrt/C_extension/include/grt/common/RT_matrix.h b/pygrt/C_extension/include/grt/common/RT_matrix.h index 27c1187a..7155de13 100644 --- a/pygrt/C_extension/include/grt/common/RT_matrix.h +++ b/pygrt/C_extension/include/grt/common/RT_matrix.h @@ -16,18 +16,18 @@ /** 应该把 PSV 和 SH 的小矩阵统一放在一个结构体中管理 */ typedef struct { - MYCOMPLEX RD[2][2]; ///< P-SV 下传反射系数矩阵 - MYCOMPLEX RU[2][2]; ///< P-SV 上传反射系数矩阵 - MYCOMPLEX TD[2][2]; ///< P-SV 下传透射系数矩阵 - MYCOMPLEX TU[2][2]; ///< P-SV 上传透射系数矩阵 - MYCOMPLEX RDL; ///< SH 下传反射系数 - MYCOMPLEX RUL; ///< SH 上传反射系数 - MYCOMPLEX TDL; ///< SH 下传透射系数 - MYCOMPLEX TUL; ///< SH 上传透射系数 + cplx_t RD[2][2]; ///< P-SV 下传反射系数矩阵 + cplx_t RU[2][2]; ///< P-SV 上传反射系数矩阵 + cplx_t TD[2][2]; ///< P-SV 下传透射系数矩阵 + cplx_t TU[2][2]; ///< P-SV 上传透射系数矩阵 + cplx_t RDL; ///< SH 下传反射系数 + cplx_t RUL; ///< SH 上传反射系数 + cplx_t TDL; ///< SH 下传透射系数 + cplx_t TUL; ///< SH 上传透射系数 /* 一些辅助变量 */ - MYCOMPLEX invT[2][2]; - MYCOMPLEX invTL; + cplx_t invT[2][2]; + cplx_t invTL; /* 状态变量 */ MYINT stats; ///< 是否有除零错误,非0为异常值 diff --git a/pygrt/C_extension/include/grt/common/attenuation.h b/pygrt/C_extension/include/grt/common/attenuation.h index 921b10a2..4254d5a2 100755 --- a/pygrt/C_extension/include/grt/common/attenuation.h +++ b/pygrt/C_extension/include/grt/common/attenuation.h @@ -23,9 +23,9 @@ * * @return atncoef 系数因子,作用在 \f$ k=\omega / c(\omega)\f$的计算 */ -MYCOMPLEX grt_attenuation_law(MYREAL Qinv, MYCOMPLEX omega); +cplx_t grt_attenuation_law(real_t Qinv, cplx_t omega); /** * attenuation_law函数在python中被调用的版本,长度2的数组分别表示复数的实部和虚部 */ -void grt_py_attenuation_law(MYREAL Qinv, MYREAL omg[2], MYREAL atte[2]); \ No newline at end of file +void grt_py_attenuation_law(real_t Qinv, real_t omg[2], real_t atte[2]); \ No newline at end of file diff --git a/pygrt/C_extension/include/grt/common/bessel.h b/pygrt/C_extension/include/grt/common/bessel.h index 5d17d248..c19bff53 100755 --- a/pygrt/C_extension/include/grt/common/bessel.h +++ b/pygrt/C_extension/include/grt/common/bessel.h @@ -19,7 +19,7 @@ * @param[out] bj2 \f$ J_2(x) \f$ * */ -void grt_bessel012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2); +void grt_bessel012(real_t x, real_t *bj0, real_t *bj1, real_t *bj2); /** @@ -31,4 +31,4 @@ void grt_bessel012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2); * @param[in,out] bj2 传入 \f$ J_2(x) \f$, 返回\f$ J_2^{'}(x) \f$ * */ -void grt_besselp012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2); \ No newline at end of file +void grt_besselp012(real_t x, real_t *bj0, real_t *bj1, real_t *bj2); \ No newline at end of file diff --git a/pygrt/C_extension/include/grt/common/const.h b/pygrt/C_extension/include/grt/common/const.h index fed2c284..ebdb6180 100755 --- a/pygrt/C_extension/include/grt/common/const.h +++ b/pygrt/C_extension/include/grt/common/const.h @@ -32,19 +32,17 @@ #endif +typedef int MYINT; ///< 整数 -// #define GRT_USE_FLOAT ///< 是否使用单精度浮点数 +typedef double real_t; +typedef double complex cplx_t; -typedef int MYINT; ///< 整数 +#define GRT_NC_MYINT NC_INT +#define GRT_NC_FUNC_MYINT(func) func##_int +#define NC_REAL NC_DOUBLE +#define NC_FUNC_REAL(func) func##_double -#ifdef GRT_USE_FLOAT - typedef float _Complex MYCOMPLEX; ///< 复数 - typedef float MYREAL; ///< 浮点数 -#else - typedef double _Complex MYCOMPLEX; - typedef double MYREAL; -#endif // 常数 #define RTWOTHIRD 0.6666666666666667 ///< 2/3 diff --git a/pygrt/C_extension/include/grt/common/dwm.h b/pygrt/C_extension/include/grt/common/dwm.h index f737ebab..4a93a95e 100644 --- a/pygrt/C_extension/include/grt/common/dwm.h +++ b/pygrt/C_extension/include/grt/common/dwm.h @@ -41,11 +41,11 @@ * * @return k 积分截至时的波数 */ -MYREAL grt_discrete_integ( - GRT_MODEL1D *mod1d, MYREAL dk, MYREAL kmax, MYREAL keps, - MYINT nr, MYREAL *rs, - MYCOMPLEX sum_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], +real_t grt_discrete_integ( + GRT_MODEL1D *mod1d, real_t dk, real_t kmax, real_t keps, + MYINT nr, real_t *rs, + cplx_t sum_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], bool calc_upar, - MYCOMPLEX sum_uiz_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], - MYCOMPLEX sum_uir_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], + cplx_t sum_uiz_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], + cplx_t sum_uir_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], FILE *fstats, GRT_KernelFunc kerfunc); diff --git a/pygrt/C_extension/include/grt/common/fim.h b/pygrt/C_extension/include/grt/common/fim.h index 5c2e167a..896f5618 100755 --- a/pygrt/C_extension/include/grt/common/fim.h +++ b/pygrt/C_extension/include/grt/common/fim.h @@ -48,13 +48,13 @@ * * @return k 积分截至时的波数 */ -MYREAL grt_linear_filon_integ( - GRT_MODEL1D *mod1d, MYREAL k0, MYREAL dk0, MYREAL filondk, MYREAL kmax, MYREAL keps, - MYINT nr, MYREAL *rs, - MYCOMPLEX sum_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], +real_t grt_linear_filon_integ( + GRT_MODEL1D *mod1d, real_t k0, real_t dk0, real_t filondk, real_t kmax, real_t keps, + MYINT nr, real_t *rs, + cplx_t sum_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], bool calc_upar, - MYCOMPLEX sum_uiz_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], - MYCOMPLEX sum_uir_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], + cplx_t sum_uiz_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], + cplx_t sum_uir_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], FILE *fstats, GRT_KernelFunc kerfunc); diff --git a/pygrt/C_extension/include/grt/common/integral.h b/pygrt/C_extension/include/grt/common/integral.h index 0d368799..2db3524a 100644 --- a/pygrt/C_extension/include/grt/common/integral.h +++ b/pygrt/C_extension/include/grt/common/integral.h @@ -25,10 +25,10 @@ * */ void grt_int_Pk( - MYREAL k, MYREAL r, - const MYCOMPLEX QWV[GRT_SRC_M_NUM][GRT_QWV_NUM], + real_t k, real_t r, + const cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM], bool calc_uir, - MYCOMPLEX SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM]); + cplx_t SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM]); @@ -40,7 +40,7 @@ void grt_int_Pk( * @param[out] tol Z、R、T分量结果 */ void grt_merge_Pk( - const MYCOMPLEX sum_J[GRT_SRC_M_NUM][GRT_INTEG_NUM], MYCOMPLEX tol[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]); + const cplx_t sum_J[GRT_SRC_M_NUM][GRT_INTEG_NUM], cplx_t tol[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]); @@ -59,10 +59,10 @@ void grt_merge_Pk( * */ void grt_int_Pk_filon( - MYREAL k, MYREAL r, bool iscos, - const MYCOMPLEX QWV[GRT_SRC_M_NUM][GRT_QWV_NUM], + real_t k, real_t r, bool iscos, + const cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM], bool calc_uir, - MYCOMPLEX SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM]); + cplx_t SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM]); /** @@ -78,7 +78,7 @@ void grt_int_Pk_filon( * */ void grt_int_Pk_sa_filon( - const MYREAL k3[3], MYREAL r, - const MYCOMPLEX QWV3[3][GRT_SRC_M_NUM][GRT_QWV_NUM], + const real_t k3[3], real_t r, + const cplx_t QWV3[3][GRT_SRC_M_NUM][GRT_QWV_NUM], bool calc_uir, - MYCOMPLEX SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM]); \ No newline at end of file + cplx_t SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM]); \ No newline at end of file diff --git a/pygrt/C_extension/include/grt/common/iostats.h b/pygrt/C_extension/include/grt/common/iostats.h index 36ced6e5..77c7faf6 100755 --- a/pygrt/C_extension/include/grt/common/iostats.h +++ b/pygrt/C_extension/include/grt/common/iostats.h @@ -27,7 +27,7 @@ * 记录其值主要用于参考其变化趋势。 */ void grt_write_stats( - FILE *f0, MYREAL k, const MYCOMPLEX QWV[GRT_SRC_M_NUM][GRT_QWV_NUM]); + FILE *f0, real_t k, const cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM]); /** @@ -55,8 +55,8 @@ MYINT grt_extract_stats(FILE *bf0, FILE *af0); */ void grt_write_stats_ptam( FILE *f0, - MYREAL Kpt[GRT_SRC_M_NUM][GRT_INTEG_NUM][GRT_PTAM_PT_MAX], - MYCOMPLEX Fpt[GRT_SRC_M_NUM][GRT_INTEG_NUM][GRT_PTAM_PT_MAX]); + real_t Kpt[GRT_SRC_M_NUM][GRT_INTEG_NUM][GRT_PTAM_PT_MAX], + cplx_t Fpt[GRT_SRC_M_NUM][GRT_INTEG_NUM][GRT_PTAM_PT_MAX]); /** diff --git a/pygrt/C_extension/include/grt/common/kernel.h b/pygrt/C_extension/include/grt/common/kernel.h index f710fcd1..4c5efa84 100644 --- a/pygrt/C_extension/include/grt/common/kernel.h +++ b/pygrt/C_extension/include/grt/common/kernel.h @@ -17,5 +17,5 @@ * 计算核函数的函数指针,动态与静态的接口一致 */ typedef void (*GRT_KernelFunc) ( - const GRT_MODEL1D *mod1d, MYCOMPLEX QWV[GRT_SRC_M_NUM][GRT_QWV_NUM], - bool calc_uiz, MYCOMPLEX QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]); \ No newline at end of file + const GRT_MODEL1D *mod1d, cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM], + bool calc_uiz, cplx_t QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]); \ No newline at end of file diff --git a/pygrt/C_extension/include/grt/common/matrix.h b/pygrt/C_extension/include/grt/common/matrix.h index 2dd031f6..0732fa56 100755 --- a/pygrt/C_extension/include/grt/common/matrix.h +++ b/pygrt/C_extension/include/grt/common/matrix.h @@ -20,12 +20,12 @@ * @param[out] invM 逆矩阵 * @param[out] stats 状态代码,是否有除零错误,非0为异常值 */ -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]; - MYCOMPLEX M11 = M[1][1]; - MYCOMPLEX det = M00*M11 - M01*M10; +inline GCC_ALWAYS_INLINE void grt_cmat2x2_inv(const cplx_t M[2][2], cplx_t invM[2][2], MYINT *stats) { + cplx_t M00 = M[0][0]; + cplx_t M01 = M[0][1]; + cplx_t M10 = M[1][0]; + cplx_t M11 = M[1][1]; + cplx_t det = M00*M11 - M01*M10; if ( det == 0.0 ){ // fprintf(stderr, "%.5e+%.5ej %.5e+%.5ej \n", creal(M[0][0]), cimag(M[0][0]), creal(M[0][1]), cimag(M[0][1])); // fprintf(stderr, "%.5e+%.5ej %.5e+%.5ej \n", creal(M[1][0]), cimag(M[1][0]), creal(M[1][1]), cimag(M[1][1])); @@ -50,7 +50,7 @@ inline GCC_ALWAYS_INLINE void grt_cmat2x2_inv(const MYCOMPLEX M[2][2], MYCOMPLEX * @param[in] M2 矩阵2 * @param[out] M 和矩阵 */ -inline GCC_ALWAYS_INLINE void grt_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 cplx_t M1[2][2], const cplx_t M2[2][2], cplx_t 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]; @@ -64,7 +64,7 @@ inline GCC_ALWAYS_INLINE void grt_cmat2x2_add(const MYCOMPLEX M1[2][2], const MY * @param[in] M2 矩阵2 * @param[out] M 差矩阵, M1 - M2 */ -inline GCC_ALWAYS_INLINE void grt_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 cplx_t M1[2][2], const cplx_t M2[2][2], cplx_t 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]; @@ -76,7 +76,7 @@ inline GCC_ALWAYS_INLINE void grt_cmat2x2_sub(const MYCOMPLEX M1[2][2], const MY * * @param[in,out] M 差矩阵 I-M2 */ -inline GCC_ALWAYS_INLINE void grt_cmat2x2_one_sub(MYCOMPLEX M[2][2]){ +inline GCC_ALWAYS_INLINE void grt_cmat2x2_one_sub(cplx_t M[2][2]){ M[0][0] = 1.0 - M[0][0]; M[0][1] = - M[0][1]; M[1][0] = - M[1][0]; @@ -90,9 +90,9 @@ inline GCC_ALWAYS_INLINE void grt_cmat2x2_one_sub(MYCOMPLEX M[2][2]){ * @param[in] M2 矩阵2 * @param[out] M 积矩阵, M1 * M2 */ -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; +inline GCC_ALWAYS_INLINE void grt_cmat2x2_mul(const cplx_t M1[2][2], const cplx_t M2[2][2], cplx_t M[2][2]){ + cplx_t M011, M012, M021, M022; + cplx_t M111, M112, M121, M122; M011 = M1[0][0]; M012 = M1[0][1]; M021 = M1[1][0]; M022 = M1[1][1]; M111 = M2[0][0]; M112 = M2[0][1]; @@ -110,7 +110,7 @@ inline GCC_ALWAYS_INLINE void grt_cmat2x2_mul(const MYCOMPLEX M1[2][2], const MY * @param[in] k 常数 * @param[out] M 积矩阵, k * M2 */ -inline GCC_ALWAYS_INLINE void grt_cmat2x2_k(const MYCOMPLEX M1[2][2], MYCOMPLEX k0, MYCOMPLEX M[2][2]){ +inline GCC_ALWAYS_INLINE void grt_cmat2x2_k(const cplx_t M1[2][2], cplx_t k0, cplx_t M[2][2]){ M[0][0] = M1[0][0] * k0; M[0][1] = M1[0][1] * k0; M[1][0] = M1[1][0] * k0; @@ -124,8 +124,8 @@ inline GCC_ALWAYS_INLINE void grt_cmat2x2_k(const MYCOMPLEX M1[2][2], MYCOMPLEX * @param[in] M2 向量 * @param[out] M 积矩阵, M1 * M2 */ -inline GCC_ALWAYS_INLINE void grt_cmat2x1_mul(const MYCOMPLEX M1[2][2], const MYCOMPLEX M2[2], MYCOMPLEX M[2]){ - MYCOMPLEX M00, M10; +inline GCC_ALWAYS_INLINE void grt_cmat2x1_mul(const cplx_t M1[2][2], const cplx_t M2[2], cplx_t M[2]){ + cplx_t M00, M10; M00 = M1[0][0]*M2[0] + M1[0][1]*M2[1]; M10 = M1[1][0]*M2[0] + M1[1][1]*M2[1]; M[0] = M00; @@ -138,7 +138,7 @@ inline GCC_ALWAYS_INLINE void grt_cmat2x1_mul(const MYCOMPLEX M1[2][2], const MY * @param[in] M1 源矩阵 * @param[out] M2 目标矩阵 */ -inline GCC_ALWAYS_INLINE void grt_cmat2x2_assign(const MYCOMPLEX M1[2][2], MYCOMPLEX M2[2][2]){ +inline GCC_ALWAYS_INLINE void grt_cmat2x2_assign(const cplx_t M1[2][2], cplx_t M2[2][2]){ M2[0][0] = M1[0][0]; M2[0][1] = M1[0][1]; M2[1][0] = M1[1][0]; @@ -155,9 +155,9 @@ inline GCC_ALWAYS_INLINE void grt_cmat2x2_assign(const MYCOMPLEX M1[2][2], MYCOM * @param[in] M2 M2矩阵 * @param[out] M 积矩阵 M1 * M2 */ -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]){ +inline GCC_ALWAYS_INLINE void grt_cmatmxn_mul(MYINT m1, MYINT n1, MYINT p1, const cplx_t M1[m1][n1], const cplx_t M2[n1][p1], cplx_t M[m1][p1]){ MYINT m, n, k; - MYCOMPLEX M0[m1][p1]; + cplx_t M0[m1][p1]; for(m=0; m=3) - MYREAL depsrc; ///< 震源深度 km - MYREAL deprcv; ///< 接收点深度 km + real_t depsrc; ///< 震源深度 km + real_t deprcv; ///< 接收点深度 km MYINT isrc; ///< 震源所在虚拟层位, isrc>=1 MYINT ircv; ///< 接收点所在虚拟层位, ircv>=1, ircv != isrc bool ircvup; ///< 接收点位于浅层, ircv < isrc bool io_depth; ///< 读取的模型首列为每层顶界面深度 - MYCOMPLEX omega; ///< 圆频率 - MYREAL k; ///< 波数 - MYCOMPLEX c_phase; ///< 当前相速度 - - MYREAL *Thk; ///< Thk[n], 最后一层厚度不使用(当作正无穷), km - MYREAL *Dep; ///< Dep[n], 每一层顶界面深度,第一层必须为 0.0 - MYREAL *Va; ///< Va[n] P波速度 km/s - MYREAL *Vb; ///< Vb[n] S波速度 km/s - MYREAL *Rho; ///< Rho[n] 密度 g/cm^3 - MYREAL *Qa; ///< Qa[n] P波Q值 - MYREAL *Qb; ///< Qb[n] S波Q值 - MYREAL *Qainv; ///< 1/Q_p - MYREAL *Qbinv; ///< 1/Q_s - - MYCOMPLEX *mu; ///< mu[n] \f$ V_b^2 * \rho \f$ - MYCOMPLEX *lambda; ///< lambda[n] \f$ V_a^2 * \rho - 2*\mu \f$ - MYCOMPLEX *delta; ///< delta[n] \f$ (\lambda+\mu)/(\lambda+3*\mu) \f$ - MYCOMPLEX *atna; - MYCOMPLEX *atnb; - MYCOMPLEX *xa; - MYCOMPLEX *xb; - MYCOMPLEX *caca; - MYCOMPLEX *cbcb; + cplx_t omega; ///< 圆频率 + real_t k; ///< 波数 + cplx_t c_phase; ///< 当前相速度 + + real_t *Thk; ///< Thk[n], 最后一层厚度不使用(当作正无穷), km + real_t *Dep; ///< Dep[n], 每一层顶界面深度,第一层必须为 0.0 + real_t *Va; ///< Va[n] P波速度 km/s + real_t *Vb; ///< Vb[n] S波速度 km/s + real_t *Rho; ///< Rho[n] 密度 g/cm^3 + real_t *Qa; ///< Qa[n] P波Q值 + real_t *Qb; ///< Qb[n] S波Q值 + real_t *Qainv; ///< 1/Q_p + real_t *Qbinv; ///< 1/Q_s + + cplx_t *mu; ///< mu[n] \f$ V_b^2 * \rho \f$ + cplx_t *lambda; ///< lambda[n] \f$ V_a^2 * \rho - 2*\mu \f$ + cplx_t *delta; ///< delta[n] \f$ (\lambda+\mu)/(\lambda+3*\mu) \f$ + cplx_t *atna; + cplx_t *atnb; + cplx_t *xa; + cplx_t *xb; + cplx_t *caca; + cplx_t *cbcb; /* 状态变量,非 0 为异常值 */ MYINT stats; @@ -93,7 +93,7 @@ GRT_MODEL1D * grt_copy_mod1d(const GRT_MODEL1D *mod1d1); * @param[in,out] mod1d `MODEL1D` 结构体指针 * @param[in] omega 复数频率 */ -void grt_attenuate_mod1d(GRT_MODEL1D *mod1d, MYCOMPLEX omega); +void grt_attenuate_mod1d(GRT_MODEL1D *mod1d, cplx_t omega); /** * 根据记录好的圆频率和波数,计算相速度和每层的 xa, xb, caca, cbcb @@ -101,7 +101,7 @@ void grt_attenuate_mod1d(GRT_MODEL1D *mod1d, MYCOMPLEX omega); * @param[in,out] mod1d 模型结构体指针 * @param[in] k 波数 */ -void grt_mod1d_xa_xb(GRT_MODEL1D *mod1d, const MYREAL k); +void grt_mod1d_xa_xb(GRT_MODEL1D *mod1d, const real_t k); /** @@ -146,7 +146,7 @@ void grt_get_model_diglen_from_file(const char *command, const char *modelpath, * * @return 是否存在 */ -bool grt_check_vel_in_mod(const GRT_MODEL1D *mod1d, const MYREAL vel, const MYREAL tol); +bool grt_check_vel_in_mod(const GRT_MODEL1D *mod1d, const real_t vel, const real_t tol); /** * 计算最大最小速度(非零值) diff --git a/pygrt/C_extension/include/grt/common/myfftw.h b/pygrt/C_extension/include/grt/common/myfftw.h index af1c0dbd..c94cf28b 100644 --- a/pygrt/C_extension/include/grt/common/myfftw.h +++ b/pygrt/C_extension/include/grt/common/myfftw.h @@ -37,13 +37,13 @@ typedef struct {\ /*时间序列长度*/ \ MYINT nt; \ /* 时间间隔 */ \ - MYREAL dt; \ + real_t dt; \ /*有效频谱长度*/\ MYINT nf_valid; \ /*总频谱长度 nf = nt/2+1*/\ MYINT nf; \ /*频率间隔*/\ - MYREAL df; \ + real_t df; \ /*时间域实序列*/\ T *w_t; \ /*频率域复序列*/\ @@ -52,7 +52,7 @@ typedef struct {\ fftw##s##_plan plan; \ /*为一些特殊情况而准备的变量*/ \ /*发生频移*/\ - MYREAL f0; \ + real_t f0; \ /*最终不执行FFTW而使用最朴素的傅里叶逆变换*/\ bool naive_inv; \ } GRT_FFTW##S##_HOLDER;\ @@ -68,9 +68,9 @@ typedef struct {\ * @return FFTW_HOLDER 结构体指针 * */\ -GRT_FFTW##S##_HOLDER * grt_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 real_t dt, const MYINT nf_valid, const real_t df); \ /** 初始化 FFTW_HOLDER 结构体指针,进行实数序列到复数频谱的正变换 */\ -GRT_FFTW##S##_HOLDER * grt_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 real_t dt, const MYINT nf_valid, const real_t df); \ \ /** 将内部数据全部置零 */\ void grt_reset_fftw##s##_holder_zero(GRT_FFTW##S##_HOLDER *fh);\ diff --git a/pygrt/C_extension/include/grt/common/mynetcdf.h b/pygrt/C_extension/include/grt/common/mynetcdf.h index 9c40a9d3..6b7b19cd 100644 --- a/pygrt/C_extension/include/grt/common/mynetcdf.h +++ b/pygrt/C_extension/include/grt/common/mynetcdf.h @@ -13,7 +13,7 @@ #include "grt/common/const.h" -#define GRT_NC_CHECK(call) ({\ +#define NC_CHECK(call) ({\ int status = (call); \ if (status != NC_NOERR) { \ GRTRaiseError("NetCDF error at %s:%d: %s\n", \ @@ -21,15 +21,3 @@ } \ }) - -#define GRT_NC_MYINT NC_INT -#define GRT_NC_FUNC_MYINT(func) func##_int - -#ifdef GRT_USE_FLOAT - #define GRT_NC_MYREAL NC_FLOAT - #define GRT_NC_FUNC_MYREAL(func) func##_float -#else - #define GRT_NC_MYREAL NC_DOUBLE - #define GRT_NC_FUNC_MYREAL(func) func##_double -#endif - diff --git a/pygrt/C_extension/include/grt/common/ptam.h b/pygrt/C_extension/include/grt/common/ptam.h index a6ef3b51..51eea166 100755 --- a/pygrt/C_extension/include/grt/common/ptam.h +++ b/pygrt/C_extension/include/grt/common/ptam.h @@ -45,12 +45,12 @@ * */ void grt_PTA_method( - GRT_MODEL1D *mod1d, MYREAL k0, MYREAL predk, - MYINT nr, MYREAL *rs, - MYCOMPLEX sum_J0[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], + GRT_MODEL1D *mod1d, real_t k0, real_t predk, + MYINT nr, real_t *rs, + cplx_t sum_J0[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], bool calc_upar, - MYCOMPLEX sum_uiz_J0[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], - MYCOMPLEX sum_uir_J0[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], + cplx_t sum_uiz_J0[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], + cplx_t sum_uir_J0[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], FILE *ptam_fstatsnr[nr][2], GRT_KernelFunc kerfunc); @@ -72,8 +72,8 @@ void grt_PTA_method( * */ MYINT grt_cplx_peak_or_trough( - MYINT idx1, MYINT idx2, const MYCOMPLEX arr[GRT_PTAM_WINDOW_SIZE][GRT_SRC_M_NUM][GRT_INTEG_NUM], - MYREAL k, MYREAL dk, MYREAL *pk, MYCOMPLEX *value); + MYINT idx1, MYINT idx2, const cplx_t arr[GRT_PTAM_WINDOW_SIZE][GRT_SRC_M_NUM][GRT_INTEG_NUM], + real_t k, real_t dk, real_t *pk, cplx_t *value); /** @@ -86,4 +86,4 @@ MYINT grt_cplx_peak_or_trough( * @param[in,out] arr 振荡的数组,最终收敛值在第一个,arr[0] * */ -void grt_cplx_shrink(MYINT n1, MYCOMPLEX *arr); \ No newline at end of file +void grt_cplx_shrink(MYINT n1, cplx_t *arr); \ No newline at end of file diff --git a/pygrt/C_extension/include/grt/common/quadratic.h b/pygrt/C_extension/include/grt/common/quadratic.h index 5531185d..c5b7c24f 100755 --- a/pygrt/C_extension/include/grt/common/quadratic.h +++ b/pygrt/C_extension/include/grt/common/quadratic.h @@ -25,7 +25,7 @@ * @param[out] pb 拟合b值 * @param[out] pc 拟合c值 */ -void grt_quad_term(const MYREAL x[3], const MYCOMPLEX f[3], MYCOMPLEX *pa, MYCOMPLEX *pb, MYCOMPLEX *pc); +void grt_quad_term(const real_t x[3], const cplx_t f[3], cplx_t *pa, cplx_t *pb, cplx_t *pc); /** @@ -38,7 +38,7 @@ void grt_quad_term(const MYREAL x[3], const MYCOMPLEX f[3], MYCOMPLEX *pa, MYCOM * @return \f$ f(x) = ax^2 + bx + c \f$ * */ -MYCOMPLEX grt_quad_eval(MYREAL x, MYCOMPLEX a, MYCOMPLEX b, MYCOMPLEX c); +cplx_t grt_quad_eval(real_t x, cplx_t a, cplx_t b, cplx_t c); /** @@ -52,4 +52,4 @@ MYCOMPLEX grt_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 grt_quad_integral(MYREAL x1, MYREAL x2, MYCOMPLEX a, MYCOMPLEX b, MYCOMPLEX c); \ No newline at end of file +cplx_t grt_quad_integral(real_t x1, real_t x2, cplx_t a, cplx_t b, cplx_t c); \ No newline at end of file diff --git a/pygrt/C_extension/include/grt/common/safim.h b/pygrt/C_extension/include/grt/common/safim.h index 6614d9e6..4bd9d1c9 100755 --- a/pygrt/C_extension/include/grt/common/safim.h +++ b/pygrt/C_extension/include/grt/common/safim.h @@ -49,13 +49,13 @@ * * @return k 积分截至时的波数 */ -MYREAL grt_sa_filon_integ( - GRT_MODEL1D *mod1d, MYREAL k0, MYREAL dk0, MYREAL tol, MYREAL kmax, MYREAL kref, - MYINT nr, MYREAL *rs, - MYCOMPLEX sum_J0[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], +real_t grt_sa_filon_integ( + GRT_MODEL1D *mod1d, real_t k0, real_t dk0, real_t tol, real_t kmax, real_t kref, + MYINT nr, real_t *rs, + cplx_t sum_J0[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], bool calc_upar, - MYCOMPLEX sum_uiz_J0[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], - MYCOMPLEX sum_uir_J0[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], + cplx_t sum_uiz_J0[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], + cplx_t sum_uir_J0[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], FILE *fstats, GRT_KernelFunc kerfunc); diff --git a/pygrt/C_extension/include/grt/common/search.h b/pygrt/C_extension/include/grt/common/search.h index 8f2df742..dadba160 100644 --- a/pygrt/C_extension/include/grt/common/search.h +++ b/pygrt/C_extension/include/grt/common/search.h @@ -13,7 +13,7 @@ // 定义 X 宏,为多个类型定义查找函数 #define __FOR_EACH_TYPE \ - X(MYINT) X(MYREAL) X(float) X(double) + X(MYINT) X(real_t) X(float) X(double) /** diff --git a/pygrt/C_extension/include/grt/common/util.h b/pygrt/C_extension/include/grt/common/util.h index fc63671e..38c04d44 100644 --- a/pygrt/C_extension/include/grt/common/util.h +++ b/pygrt/C_extension/include/grt/common/util.h @@ -125,12 +125,12 @@ ssize_t grt_getline(char **lineptr, size_t *n, FILE *stream); void grt_GF_freq2time_write_to_file( const char *command, const GRT_MODEL1D *mod1d, const char *s_output_dir, const char *s_modelname, const char *s_depsrc, const char *s_deprcv, - 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, + const real_t wI, GRT_FFTW_HOLDER *pt_fh, + const MYINT nr, char *s_dists[nr], const real_t dists[nr], real_t travtPS[nr][2], + const real_t depsrc, const real_t deprcv, + const real_t delayT0, const real_t delayV0, const bool calc_upar, const bool doEX, const bool doVF, const bool doHF, const bool doDC, const char *chalst, - MYCOMPLEX *grn[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM], - MYCOMPLEX *grn_uiz[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM], - MYCOMPLEX *grn_uir[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM]); \ No newline at end of file + cplx_t *grn[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM], + cplx_t *grn_uiz[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM], + cplx_t *grn_uir[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM]); \ No newline at end of file diff --git a/pygrt/C_extension/include/grt/dynamic/grn.h b/pygrt/C_extension/include/grt/dynamic/grn.h index df870eb3..328885a5 100755 --- a/pygrt/C_extension/include/grt/dynamic/grn.h +++ b/pygrt/C_extension/include/grt/dynamic/grn.h @@ -49,17 +49,17 @@ * */ void grt_integ_grn_spec( - GRT_MODEL1D *mod1d, 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, + GRT_MODEL1D *mod1d, MYINT nf1, MYINT nf2, real_t *freqs, + MYINT nr, real_t *rs, real_t wI, + real_t vmin_ref, real_t keps, real_t ampk, real_t k0, real_t Length, real_t filonLength, real_t safilonTol, real_t filonCut, bool print_progressbar, // 返回值,代表Z、R、T分量 - MYCOMPLEX *grn[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM], + cplx_t *grn[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM], bool calc_upar, - MYCOMPLEX *grn_uiz[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM], - MYCOMPLEX *grn_uir[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM], + cplx_t *grn_uiz[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM], + cplx_t *grn_uir[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM], const char *statsstr, // 积分过程输出 MYINT nstatsidxs, // 仅输出特定频点 diff --git a/pygrt/C_extension/include/grt/dynamic/layer.h b/pygrt/C_extension/include/grt/dynamic/layer.h index 80806d02..6574bdd4 100755 --- a/pygrt/C_extension/include/grt/dynamic/layer.h +++ b/pygrt/C_extension/include/grt/dynamic/layer.h @@ -35,7 +35,7 @@ void grt_topfree_RU(const GRT_MODEL1D *mod1d, RT_MATRIX *M); * @param[out] R_EV P-SV接收函数矩阵 * */ -void grt_wave2qwv_REV_PSV(const GRT_MODEL1D *mod1d, MYCOMPLEX R[2][2], MYCOMPLEX R_EV[2][2]); +void grt_wave2qwv_REV_PSV(const GRT_MODEL1D *mod1d, cplx_t R[2][2], cplx_t R_EV[2][2]); /** * 计算接收点位置的 SH 波接收矩阵,将波场转为位移,公式(5.2.19) + (5.7.7,25) @@ -45,7 +45,7 @@ void grt_wave2qwv_REV_PSV(const GRT_MODEL1D *mod1d, MYCOMPLEX R[2][2], MYCOMPLEX * @param[out] R_EVL SH接收函数值 * */ -void grt_wave2qwv_REV_SH(const GRT_MODEL1D *mod1d, MYCOMPLEX RL, MYCOMPLEX *R_EVL); +void grt_wave2qwv_REV_SH(const GRT_MODEL1D *mod1d, cplx_t RL, cplx_t *R_EVL); /** @@ -57,7 +57,7 @@ void grt_wave2qwv_REV_SH(const GRT_MODEL1D *mod1d, MYCOMPLEX RL, MYCOMPLEX *R_EV * @param[out] R_EV P-SV接收函数矩阵 * */ -void grt_wave2qwv_z_REV_PSV(const GRT_MODEL1D *mod1d, const MYCOMPLEX R[2][2], MYCOMPLEX R_EV[2][2]); +void grt_wave2qwv_z_REV_PSV(const GRT_MODEL1D *mod1d, const cplx_t R[2][2], cplx_t R_EV[2][2]); /** @@ -69,7 +69,7 @@ void grt_wave2qwv_z_REV_PSV(const GRT_MODEL1D *mod1d, const MYCOMPLEX R[2][2], M * @param[out] R_EVL SH接收函数值 * */ -void grt_wave2qwv_z_REV_SH(const GRT_MODEL1D *mod1d, MYCOMPLEX RL, MYCOMPLEX *R_EVL); +void grt_wave2qwv_z_REV_SH(const GRT_MODEL1D *mod1d, cplx_t RL, cplx_t *R_EVL); /** @@ -147,34 +147,34 @@ void grt_delay_RT_matrix(const GRT_MODEL1D *mod1d, const MYINT iy, RT_MATRIX *M) * */ void grt_get_layer_D( - MYCOMPLEX xa, MYCOMPLEX xb, MYCOMPLEX kbkb, MYCOMPLEX mu, - MYCOMPLEX omega, MYREAL rho, MYREAL k, MYCOMPLEX D[4][4], bool inverse, MYINT liquid_invtype); + cplx_t xa, cplx_t xb, cplx_t kbkb, cplx_t mu, + cplx_t omega, real_t rho, real_t k, cplx_t D[4][4], bool inverse, MYINT liquid_invtype); /** 子矩阵 D11,函数参数见 get_layer_D 函数 */ void grt_get_layer_D11( - MYCOMPLEX xa, MYCOMPLEX xb, MYREAL k, MYCOMPLEX D[2][2]); + cplx_t xa, cplx_t xb, real_t k, cplx_t D[2][2]); /** 子矩阵 D12,函数参数见 get_layer_D 函数 */ void grt_get_layer_D12( - MYCOMPLEX xa, MYCOMPLEX xb, MYREAL k, MYCOMPLEX D[2][2]); + cplx_t xa, cplx_t xb, real_t k, cplx_t D[2][2]); /** 子矩阵 D21,函数参数见 get_layer_D 函数 */ void grt_get_layer_D21( - MYCOMPLEX xa, MYCOMPLEX xb, MYCOMPLEX kbkb, MYCOMPLEX mu, - MYCOMPLEX omega, MYREAL rho, MYREAL k, MYCOMPLEX D[2][2]); + cplx_t xa, cplx_t xb, cplx_t kbkb, cplx_t mu, + cplx_t omega, real_t rho, real_t k, cplx_t D[2][2]); /** 子矩阵 D22,函数参数见 get_layer_D 函数 */ void grt_get_layer_D22( - MYCOMPLEX xa, MYCOMPLEX xb, MYCOMPLEX kbkb, MYCOMPLEX mu, - MYCOMPLEX omega, MYREAL rho, MYREAL k, MYCOMPLEX D[2][2]); + cplx_t xa, cplx_t xb, cplx_t kbkb, cplx_t mu, + cplx_t omega, real_t rho, real_t k, cplx_t D[2][2]); /** 子矩阵 D11_uiz,后缀uiz表示连接位移对z的偏导和垂直波函数,函数参数见 get_layer_D 函数 */ void grt_get_layer_D11_uiz( - MYCOMPLEX xa, MYCOMPLEX xb, MYREAL k, MYCOMPLEX D[2][2]); + cplx_t xa, cplx_t xb, real_t k, cplx_t D[2][2]); /** 子矩阵 D12_uiz,函数参数见 get_layer_D 函数 */ void grt_get_layer_D12_uiz( - MYCOMPLEX xa, MYCOMPLEX xb, MYREAL k, MYCOMPLEX D[2][2]); + cplx_t xa, cplx_t xb, real_t k, cplx_t D[2][2]); /** @@ -190,14 +190,14 @@ void grt_get_layer_D12_uiz( * */ void grt_get_layer_T( - MYCOMPLEX xb, MYCOMPLEX mu, - MYCOMPLEX omega, MYREAL k, MYCOMPLEX T[2][2], bool inverse); + cplx_t xb, cplx_t mu, + cplx_t omega, real_t k, cplx_t T[2][2], bool inverse); /** 计算 P-SV 型垂直波函数的时间延迟矩阵,公式(5.2.27) */ -void grt_get_layer_E_Rayl(MYCOMPLEX xa1, MYCOMPLEX xb1, MYREAL thk, MYREAL k, MYCOMPLEX E[4][4], bool inverse); +void grt_get_layer_E_Rayl(cplx_t xa1, cplx_t xb1, real_t thk, real_t k, cplx_t E[4][4], bool inverse); /** 计算 SH 型垂直波函数的时间延迟矩阵,公式(5.2.28) */ -void grt_get_layer_E_Love(MYCOMPLEX xb1, MYREAL thk, MYREAL k, MYCOMPLEX E[2][2], bool inverse); +void grt_get_layer_E_Love(cplx_t xb1, real_t thk, real_t k, cplx_t E[2][2], bool inverse); @@ -207,9 +207,9 @@ void grt_get_layer_E_Love(MYCOMPLEX xb1, MYREAL thk, MYREAL k, MYCOMPLEX E[2][2] * 函数接口也类似 */ void grt_RT_matrix_from_4x4( - MYCOMPLEX xa1, MYCOMPLEX xb1, MYCOMPLEX kbkb1, MYCOMPLEX mu1, MYREAL rho1, - MYCOMPLEX xa2, MYCOMPLEX xb2, MYCOMPLEX kbkb2, MYCOMPLEX mu2, MYREAL rho2, - MYCOMPLEX omega, MYREAL thk, - MYREAL k, - 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); \ No newline at end of file + cplx_t xa1, cplx_t xb1, cplx_t kbkb1, cplx_t mu1, real_t rho1, + cplx_t xa2, cplx_t xb2, cplx_t kbkb2, cplx_t mu2, real_t rho2, + cplx_t omega, real_t thk, + real_t k, + cplx_t RD[2][2], cplx_t *RDL, cplx_t RU[2][2], cplx_t *RUL, + cplx_t TD[2][2], cplx_t *TDL, cplx_t TU[2][2], cplx_t *TUL, MYINT *stats); \ No newline at end of file diff --git a/pygrt/C_extension/include/grt/dynamic/propagate.h b/pygrt/C_extension/include/grt/dynamic/propagate.h index 20a528f0..6df9cdd9 100755 --- a/pygrt/C_extension/include/grt/dynamic/propagate.h +++ b/pygrt/C_extension/include/grt/dynamic/propagate.h @@ -86,8 +86,8 @@ * */ void grt_kernel( - const GRT_MODEL1D *mod1d, MYCOMPLEX QWV[GRT_SRC_M_NUM][GRT_QWV_NUM], - bool calc_uiz, MYCOMPLEX QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]); + const GRT_MODEL1D *mod1d, cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM], + bool calc_uiz, cplx_t QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]); @@ -126,7 +126,7 @@ void grt_kernel( */ void grt_construct_qwv( bool ircvup, - const MYCOMPLEX R1[2][2], MYCOMPLEX RL1, - const MYCOMPLEX R2[2][2], MYCOMPLEX RL2, - const MYCOMPLEX coef_PSV[GRT_QWV_NUM-1][2], const MYCOMPLEX coef_SH[2], - MYCOMPLEX qwv[GRT_QWV_NUM]); \ No newline at end of file + const cplx_t R1[2][2], cplx_t RL1, + const cplx_t R2[2][2], cplx_t RL2, + const cplx_t coef_PSV[GRT_QWV_NUM-1][2], const cplx_t coef_SH[2], + cplx_t qwv[GRT_QWV_NUM]); \ No newline at end of file diff --git a/pygrt/C_extension/include/grt/dynamic/source.h b/pygrt/C_extension/include/grt/dynamic/source.h index 5280404d..e04c8377 100755 --- a/pygrt/C_extension/include/grt/dynamic/source.h +++ b/pygrt/C_extension/include/grt/dynamic/source.h @@ -24,7 +24,7 @@ * @param[out] coef 震源系数 \f$ P_m, SV_m, SH_m \f$ * */ -void grt_source_coef_PSV(const GRT_MODEL1D *mod1d, MYCOMPLEX coef[GRT_SRC_M_NUM][GRT_QWV_NUM-1][2]); +void grt_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM-1][2]); /* SH 波的震源系数 */ -void grt_source_coef_SH(const GRT_MODEL1D *mod1d, MYCOMPLEX coef[GRT_SRC_M_NUM][2]); +void grt_source_coef_SH(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][2]); diff --git a/pygrt/C_extension/include/grt/lamb/elliptic.h b/pygrt/C_extension/include/grt/lamb/elliptic.h index 89493fbc..3f59d1b0 100644 --- a/pygrt/C_extension/include/grt/lamb/elliptic.h +++ b/pygrt/C_extension/include/grt/lamb/elliptic.h @@ -24,7 +24,7 @@ * @return 积分值 * */ -MYREAL grt_ellipticK(const MYREAL m); +real_t grt_ellipticK(const real_t m); /** @@ -37,7 +37,7 @@ MYREAL grt_ellipticK(const MYREAL m); * @return 积分值 * */ -MYREAL grt_ellipticE(const MYREAL m); +real_t grt_ellipticE(const real_t m); /** @@ -51,4 +51,4 @@ MYREAL grt_ellipticE(const MYREAL m); * @return 积分值 * */ -MYCOMPLEX grt_ellipticPi(const MYCOMPLEX n, const MYREAL m); +cplx_t grt_ellipticPi(const cplx_t n, const real_t m); diff --git a/pygrt/C_extension/include/grt/lamb/lamb1.h b/pygrt/C_extension/include/grt/lamb/lamb1.h index cbd2f578..cc7f2b1d 100644 --- a/pygrt/C_extension/include/grt/lamb/lamb1.h +++ b/pygrt/C_extension/include/grt/lamb/lamb1.h @@ -25,4 +25,4 @@ * */ void grt_solve_lamb1( - const MYREAL nu, const MYREAL *ts, const int nt, const MYREAL azimuth, MYREAL (*u)[3][3]); \ No newline at end of file + const real_t nu, const real_t *ts, const int nt, const real_t azimuth, real_t (*u)[3][3]); \ No newline at end of file diff --git a/pygrt/C_extension/include/grt/lamb/lamb_util.h b/pygrt/C_extension/include/grt/lamb/lamb_util.h index eca4762a..62e9c0ca 100644 --- a/pygrt/C_extension/include/grt/lamb/lamb_util.h +++ b/pygrt/C_extension/include/grt/lamb/lamb_util.h @@ -18,7 +18,7 @@ * @param[in] c 系数 c * @param[out] y3 三个根,其中 y3[2] 为正根 */ -void grt_roots3(const MYREAL a, const MYREAL b, const MYREAL c, MYCOMPLEX y3[3]); +void grt_roots3(const real_t a, const real_t b, const real_t c, cplx_t y3[3]); /** * 做如下多项式求值, \f$ \sum_{m=0}^n C_{2m+o} y^m \f$ @@ -31,4 +31,4 @@ void grt_roots3(const MYREAL a, const MYREAL b, const MYREAL c, MYCOMPLEX y3[3]) * @return 多项式结果 * */ -MYCOMPLEX grt_evalpoly2(const MYCOMPLEX *C, const int n, const MYCOMPLEX y, const int offset); \ No newline at end of file +cplx_t grt_evalpoly2(const cplx_t *C, const int n, const cplx_t y, const int offset); \ No newline at end of file diff --git a/pygrt/C_extension/include/grt/static/static_grn.h b/pygrt/C_extension/include/grt/static/static_grn.h index 363e7f29..b11398b6 100644 --- a/pygrt/C_extension/include/grt/static/static_grn.h +++ b/pygrt/C_extension/include/grt/static/static_grn.h @@ -42,15 +42,15 @@ * */ void grt_integ_static_grn( - GRT_MODEL1D *mod1d, MYINT nr, MYREAL *rs, MYREAL vmin_ref, MYREAL keps, MYREAL k0, MYREAL Length, - MYREAL filonLength, MYREAL safilonTol, MYREAL filonCut, + GRT_MODEL1D *mod1d, MYINT nr, real_t *rs, real_t vmin_ref, real_t keps, real_t k0, real_t Length, + real_t filonLength, real_t safilonTol, real_t filonCut, // 返回值,代表Z、R、T分量 - MYREAL grn[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM], + real_t grn[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM], bool calc_upar, - MYREAL grn_uiz[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM], - MYREAL grn_uir[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM], + real_t grn_uiz[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM], + real_t grn_uir[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM], const char *statsstr // 积分结果输出 ); \ No newline at end of file diff --git a/pygrt/C_extension/include/grt/static/static_layer.h b/pygrt/C_extension/include/grt/static/static_layer.h index a7e029f1..2ca49745 100644 --- a/pygrt/C_extension/include/grt/static/static_layer.h +++ b/pygrt/C_extension/include/grt/static/static_layer.h @@ -37,7 +37,7 @@ void grt_static_topfree_RU(const GRT_MODEL1D *mod1d, RT_MATRIX *M); * @param[out] R_EV P-SV接收函数矩阵 * */ -void grt_static_wave2qwv_REV_PSV(const GRT_MODEL1D *mod1d, const MYCOMPLEX R[2][2], MYCOMPLEX R_EV[2][2]); +void grt_static_wave2qwv_REV_PSV(const GRT_MODEL1D *mod1d, const cplx_t R[2][2], cplx_t R_EV[2][2]); /** * 计算接收点位置的 SH 静态接收矩阵,将波场转为位移,公式(6.3.37) @@ -46,7 +46,7 @@ void grt_static_wave2qwv_REV_PSV(const GRT_MODEL1D *mod1d, const MYCOMPLEX R[2][ * @param[out] R_EVL SH接收函数值 * */ -void grt_static_wave2qwv_REV_SH(MYCOMPLEX RL, MYCOMPLEX *R_EVL); +void grt_static_wave2qwv_REV_SH(cplx_t RL, cplx_t *R_EVL); /** * 计算接收点位置的ui_z的 P-SV 静态接收矩阵,即将波场转为ui_z。 @@ -57,7 +57,7 @@ void grt_static_wave2qwv_REV_SH(MYCOMPLEX RL, MYCOMPLEX *R_EVL); * @param[out] R_EV P-SV接收函数矩阵 * */ -void grt_static_wave2qwv_z_REV_PSV(const GRT_MODEL1D *mod1d, const MYCOMPLEX R[2][2], MYCOMPLEX R_EV[2][2]); +void grt_static_wave2qwv_z_REV_PSV(const GRT_MODEL1D *mod1d, const cplx_t R[2][2], cplx_t R_EV[2][2]); /** * 计算接收点位置的ui_z的 SH 静态接收矩阵,即将波场转为ui_z。 @@ -68,7 +68,7 @@ void grt_static_wave2qwv_z_REV_PSV(const GRT_MODEL1D *mod1d, const MYCOMPLEX R[2 * @param[out] R_EVL SH接收函数值 * */ -void grt_static_wave2qwv_z_REV_SH(const GRT_MODEL1D *mod1d, MYCOMPLEX RL, MYCOMPLEX *R_EVL); +void grt_static_wave2qwv_z_REV_SH(const GRT_MODEL1D *mod1d, cplx_t RL, cplx_t *R_EVL); /** diff --git a/pygrt/C_extension/include/grt/static/static_propagate.h b/pygrt/C_extension/include/grt/static/static_propagate.h index 4e5e56ce..f8adcde1 100644 --- a/pygrt/C_extension/include/grt/static/static_propagate.h +++ b/pygrt/C_extension/include/grt/static/static_propagate.h @@ -25,5 +25,5 @@ * */ void grt_static_kernel( - const GRT_MODEL1D *mod1d, MYCOMPLEX QWV[GRT_SRC_M_NUM][GRT_QWV_NUM], - bool calc_uiz, MYCOMPLEX QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]); \ No newline at end of file + const GRT_MODEL1D *mod1d, cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM], + bool calc_uiz, cplx_t QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]); \ No newline at end of file diff --git a/pygrt/C_extension/include/grt/static/static_source.h b/pygrt/C_extension/include/grt/static/static_source.h index 0dd1588d..cc11e54c 100644 --- a/pygrt/C_extension/include/grt/static/static_source.h +++ b/pygrt/C_extension/include/grt/static/static_source.h @@ -22,7 +22,7 @@ * @param[in] mod1d 模型结构体指针 * @param[out] coef 震源系数 \f$ P_m, SV_m, SH_m \f$ */ -void grt_static_source_coef_PSV(const GRT_MODEL1D *mod1d, MYCOMPLEX coef[GRT_SRC_M_NUM][GRT_QWV_NUM-1][2]); +void grt_static_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM-1][2]); /* SH 波的静态震源系数 */ -void grt_static_source_coef_SH(const GRT_MODEL1D *mod1d, MYCOMPLEX coef[GRT_SRC_M_NUM][2]); \ No newline at end of file +void grt_static_source_coef_SH(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][2]); \ No newline at end of file diff --git a/pygrt/C_extension/include/grt/travt/travt.h b/pygrt/C_extension/include/grt/travt/travt.h index dd23518f..1ce45636 100644 --- a/pygrt/C_extension/include/grt/travt/travt.h +++ b/pygrt/C_extension/include/grt/travt/travt.h @@ -23,6 +23,6 @@ * @param[in] ircv 场点所在层位 * @param[in] dist 震中距 */ -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 +real_t grt_compute_travt1d( + const real_t *Thk, const real_t *Vel0, const int nlay, + const int isrc, const int ircv, const real_t dist); \ No newline at end of file diff --git a/pygrt/C_extension/src/common/RT_matrix.c b/pygrt/C_extension/src/common/RT_matrix.c index 652309e3..6c6c0630 100644 --- a/pygrt/C_extension/src/common/RT_matrix.c +++ b/pygrt/C_extension/src/common/RT_matrix.c @@ -27,7 +27,7 @@ void grt_recursion_RD(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MATRIX *M) void grt_recursion_RD_PSV(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MATRIX *M) { - MYCOMPLEX tmp1[2][2], tmp2[2][2]; + cplx_t tmp1[2][2], tmp2[2][2]; // RD, RDL grt_cmat2x2_mul(M1->RU, M2->RD, tmp1); @@ -43,7 +43,7 @@ void grt_recursion_RD_PSV(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MATRIX *M void grt_recursion_RD_SH(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MATRIX *M) { - MYCOMPLEX inv1; + cplx_t inv1; inv1 = 1.0 - M1->RUL * M2->RDL; if(inv1 == 0.0){ @@ -64,7 +64,7 @@ void grt_recursion_TD(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MATRIX *M) void grt_recursion_TD_PSV(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MATRIX *M) { - MYCOMPLEX tmp1[2][2], tmp2[2][2]; + cplx_t tmp1[2][2], tmp2[2][2]; // TD, TDL grt_cmat2x2_mul(M1->RU, M2->RD, tmp2); @@ -77,7 +77,7 @@ void grt_recursion_TD_PSV(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MATRIX *M void grt_recursion_TD_SH(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MATRIX *M) { - MYCOMPLEX inv1; + cplx_t inv1; inv1 = 1.0 - M1->RUL * M2->RDL; if(inv1 == 0.0){ @@ -97,7 +97,7 @@ void grt_recursion_RU(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MATRIX *M) void grt_recursion_RU_PSV(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MATRIX *M) { - MYCOMPLEX tmp1[2][2], tmp2[2][2]; + cplx_t tmp1[2][2], tmp2[2][2]; // RU, RUL grt_cmat2x2_mul(M2->RD, M1->RU, tmp2); @@ -114,7 +114,7 @@ void grt_recursion_RU_PSV(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MATRIX *M void grt_recursion_RU_SH(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MATRIX *M) { - MYCOMPLEX inv1; + cplx_t inv1; inv1 = 1.0 - M1->RUL * M2->RDL; if(inv1 == 0.0){ @@ -136,7 +136,7 @@ void grt_recursion_TU(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MATRIX *M) void grt_recursion_TU_PSV(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MATRIX *M) { - MYCOMPLEX tmp1[2][2], tmp2[2][2]; + cplx_t tmp1[2][2], tmp2[2][2]; // TU, TUL grt_cmat2x2_mul(M2->RD, M1->RU, tmp2); @@ -151,7 +151,7 @@ void grt_recursion_TU_PSV(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MATRIX *M void grt_recursion_TU_SH(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MATRIX *M) { - MYCOMPLEX inv1; + cplx_t inv1; inv1 = 1.0 - M1->RUL * M2->RDL; if(inv1 == 0.0){ @@ -174,7 +174,7 @@ void grt_recursion_RT_matrix(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MATRIX void grt_recursion_RT_matrix_PSV(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MATRIX *M) { // 临时矩阵 - MYCOMPLEX tmp1[2][2], tmp2[2][2]; + cplx_t tmp1[2][2], tmp2[2][2]; grt_cmat2x2_mul(M1->RU, M2->RD, tmp1); grt_cmat2x2_one_sub(tmp1); @@ -207,7 +207,7 @@ void grt_recursion_RT_matrix_PSV(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MA void grt_recursion_RT_matrix_SH(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MATRIX *M) { // 临时 - MYCOMPLEX inv0, inv1T; + cplx_t inv0, inv1T; inv0 = 1.0 - M1->RUL * M2->RDL; if(inv0 == 0.0){ diff --git a/pygrt/C_extension/src/common/attenuation.c b/pygrt/C_extension/src/common/attenuation.c index 3dd07d36..db6ee62e 100755 --- a/pygrt/C_extension/src/common/attenuation.c +++ b/pygrt/C_extension/src/common/attenuation.c @@ -12,15 +12,15 @@ -MYCOMPLEX grt_attenuation_law(MYREAL Qinv, MYCOMPLEX omega){ +cplx_t grt_attenuation_law(real_t Qinv, cplx_t omega){ return 1.0 + Qinv/PI * log(omega/PI2) + 0.5*Qinv*I; // return 1.0; } -void grt_py_attenuation_law(MYREAL Qinv, MYREAL omg[2], MYREAL atte[2]){ +void grt_py_attenuation_law(real_t Qinv, real_t omg[2], real_t atte[2]){ // 用于在python中调用attenuation_law - MYCOMPLEX omega = omg[0] + I*omg[1]; - MYCOMPLEX atte0 = grt_attenuation_law(Qinv, omega); + cplx_t omega = omg[0] + I*omg[1]; + cplx_t 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 05cd7a8b..a0554418 100755 --- a/pygrt/C_extension/src/common/bessel.c +++ b/pygrt/C_extension/src/common/bessel.c @@ -10,16 +10,16 @@ #include "grt/common/bessel.h" #include "grt/common/const.h" -void grt_bessel012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2){ +void grt_bessel012(real_t x, real_t *bj0, real_t *bj1, real_t *bj2){ *bj0 = j0(x); *bj1 = j1(x); *bj2 = jn(2, x); } -void grt_besselp012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2){ - MYREAL j0=*bj0; - MYREAL j1=*bj1; - MYREAL j2=*bj2; +void grt_besselp012(real_t x, real_t *bj0, real_t *bj1, real_t *bj2){ + real_t j0=*bj0; + real_t j1=*bj1; + real_t j2=*bj2; *bj0 = -j1; *bj1 = j0 - 1.0/x * j1; *bj2 = j1 - 2.0/x * j2; diff --git a/pygrt/C_extension/src/common/dwm.c b/pygrt/C_extension/src/common/dwm.c index 34069c0b..6c0ff0b9 100644 --- a/pygrt/C_extension/src/common/dwm.c +++ b/pygrt/C_extension/src/common/dwm.c @@ -23,22 +23,22 @@ #include "grt/common/const.h" -MYREAL grt_discrete_integ( - GRT_MODEL1D *mod1d, MYREAL dk, MYREAL kmax, MYREAL keps, - MYINT nr, MYREAL *rs, - MYCOMPLEX sum_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], +real_t grt_discrete_integ( + GRT_MODEL1D *mod1d, real_t dk, real_t kmax, real_t keps, + MYINT nr, real_t *rs, + cplx_t sum_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], bool calc_upar, - MYCOMPLEX sum_uiz_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], - MYCOMPLEX sum_uir_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], + cplx_t sum_uiz_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], + cplx_t sum_uir_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], FILE *fstats, GRT_KernelFunc kerfunc) { - MYCOMPLEX SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM]; + cplx_t SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM]; // 不同震源不同阶数的核函数 F(k, w) - MYCOMPLEX QWV[GRT_SRC_M_NUM][GRT_QWV_NUM]; - MYCOMPLEX QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]; + cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM]; + cplx_t QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]; - MYREAL k = 0.0; + real_t k = 0.0; MYINT ik = 0; // 所有震中距的k循环是否结束 diff --git a/pygrt/C_extension/src/common/fim.c b/pygrt/C_extension/src/common/fim.c index b288afc4..3d1d0440 100755 --- a/pygrt/C_extension/src/common/fim.c +++ b/pygrt/C_extension/src/common/fim.c @@ -22,27 +22,27 @@ -MYREAL grt_linear_filon_integ( - GRT_MODEL1D *mod1d, MYREAL k0, MYREAL dk0, MYREAL dk, MYREAL kmax, MYREAL keps, - MYINT nr, MYREAL *rs, - MYCOMPLEX sum_J0[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], +real_t grt_linear_filon_integ( + GRT_MODEL1D *mod1d, real_t k0, real_t dk0, real_t dk, real_t kmax, real_t keps, + MYINT nr, real_t *rs, + cplx_t sum_J0[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], bool calc_upar, - MYCOMPLEX sum_uiz_J0[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], - MYCOMPLEX sum_uir_J0[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], + cplx_t sum_uiz_J0[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], + cplx_t sum_uir_J0[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], FILE *fstats, GRT_KernelFunc kerfunc) { // 从0开始,存储第二部分Filon积分的结果 - MYCOMPLEX (*sum_J)[GRT_SRC_M_NUM][GRT_INTEG_NUM] = (MYCOMPLEX(*)[GRT_SRC_M_NUM][GRT_INTEG_NUM])calloc(nr, sizeof(*sum_J)); - MYCOMPLEX (*sum_uiz_J)[GRT_SRC_M_NUM][GRT_INTEG_NUM] = (calc_upar)? (MYCOMPLEX(*)[GRT_SRC_M_NUM][GRT_INTEG_NUM])calloc(nr, sizeof(*sum_uiz_J)) : NULL; - MYCOMPLEX (*sum_uir_J)[GRT_SRC_M_NUM][GRT_INTEG_NUM] = (calc_upar)? (MYCOMPLEX(*)[GRT_SRC_M_NUM][GRT_INTEG_NUM])calloc(nr, sizeof(*sum_uir_J)) : NULL; + cplx_t (*sum_J)[GRT_SRC_M_NUM][GRT_INTEG_NUM] = (cplx_t(*)[GRT_SRC_M_NUM][GRT_INTEG_NUM])calloc(nr, sizeof(*sum_J)); + cplx_t (*sum_uiz_J)[GRT_SRC_M_NUM][GRT_INTEG_NUM] = (calc_upar)? (cplx_t(*)[GRT_SRC_M_NUM][GRT_INTEG_NUM])calloc(nr, sizeof(*sum_uiz_J)) : NULL; + cplx_t (*sum_uir_J)[GRT_SRC_M_NUM][GRT_INTEG_NUM] = (calc_upar)? (cplx_t(*)[GRT_SRC_M_NUM][GRT_INTEG_NUM])calloc(nr, sizeof(*sum_uir_J)) : NULL; - MYCOMPLEX SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM]; + cplx_t SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM]; // 不同震源不同阶数的核函数 F(k, w) - MYCOMPLEX QWV[GRT_SRC_M_NUM][GRT_QWV_NUM]; - MYCOMPLEX QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]; + cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM]; + cplx_t QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]; - MYREAL k=k0; + real_t k=k0; MYINT ik=0; bool iendk, iendk0; @@ -141,7 +141,7 @@ MYREAL grt_linear_filon_integ( // ------------------------------------------------------------------------------ // 为累计项乘系数 for(MYINT ir=0; irn; ++i){ Va0 = mod1d->Va[i]; Vb0 = mod1d->Vb[i]; @@ -181,7 +181,7 @@ void grt_attenuate_mod1d(GRT_MODEL1D *mod1d, MYCOMPLEX omega){ } -void grt_mod1d_xa_xb(GRT_MODEL1D *mod1d, const MYREAL k) +void grt_mod1d_xa_xb(GRT_MODEL1D *mod1d, const real_t k) { mod1d->k = k; // 完全为0,没有虚频率,这只可能是在计算静态解,此时不需要xa, xb等物理量 @@ -190,14 +190,14 @@ void grt_mod1d_xa_xb(GRT_MODEL1D *mod1d, const MYREAL k) mod1d->c_phase = mod1d->omega/k; for(MYINT i=0; in; ++i){ - MYREAL va, vb; + real_t va, vb; va = mod1d->Va[i]; vb = mod1d->Vb[i]; - MYCOMPLEX atna, atnb; + cplx_t atna, atnb; atna = mod1d->atna[i]; atnb = mod1d->atnb[i]; - MYCOMPLEX caca, cbcb; + cplx_t caca, cbcb; caca = mod1d->c_phase / (va*atna); caca *= caca; mod1d->caca[i] = caca; @@ -437,7 +437,7 @@ void grt_get_model_diglen_from_file(const char *command, const char *modelpath, } -bool grt_check_vel_in_mod(const GRT_MODEL1D *mod1d, const MYREAL vel, const MYREAL tol){ +bool grt_check_vel_in_mod(const GRT_MODEL1D *mod1d, const real_t vel, const real_t tol){ // 浮点数比较,检查是否存在该速度值 for(MYINT i=0; in; ++i){ if(fabs(vel - mod1d->Va[i])Vb[i])Va; - const MYREAL *Vb = mod1d->Vb; + const real_t *Va = mod1d->Va; + const real_t *Vb = mod1d->Vb; for(MYINT i=0; in; ++i){ if(Va[i] < *vmin) *vmin = Va[i]; if(Va[i] > *vmax) *vmax = Va[i]; diff --git a/pygrt/C_extension/src/common/myfftw.c b/pygrt/C_extension/src/common/myfftw.c index 92003f40..5cf32f01 100644 --- a/pygrt/C_extension/src/common/myfftw.c +++ b/pygrt/C_extension/src/common/myfftw.c @@ -12,7 +12,7 @@ #include "grt/common/checkerror.h" #define X(T, s, S) \ -static GRT_FFTW##S##_HOLDER * grt_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 real_t dt, const MYINT nf_valid, const real_t df)\ {\ GRT_FFTW##S##_HOLDER *fh = (GRT_FFTW##S##_HOLDER*)calloc(1, sizeof(GRT_FFTW##S##_HOLDER));\ if (!fh) {\ @@ -41,7 +41,7 @@ void grt_reset_fftw##s##_holder_zero(GRT_FFTW##S##_HOLDER *fh)\ }\ \ \ -GRT_FFTW##S##_HOLDER * grt_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 real_t dt, const MYINT nf_valid, const real_t 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);\ @@ -49,7 +49,7 @@ GRT_FFTW##S##_HOLDER * grt_create_fftw##s##_holder_C2R_1D(const MYINT nt, const }\ \ \ -GRT_FFTW##S##_HOLDER * grt_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 real_t dt, const MYINT nf_valid, const real_t 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);\ diff --git a/pygrt/C_extension/src/common/ptam.c b/pygrt/C_extension/src/common/ptam.c index aa996d3f..1eafc64e 100755 --- a/pygrt/C_extension/src/common/ptam.c +++ b/pygrt/C_extension/src/common/ptam.c @@ -42,11 +42,11 @@ * @param[in,out] iendk0 一个布尔指针,用于指示是否满足结束条件 */ static void process_peak_or_trough( - MYINT ir, MYINT im, MYINT v, MYREAL k, MYREAL dk, - MYCOMPLEX (*J3)[GRT_PTAM_WINDOW_SIZE][GRT_SRC_M_NUM][GRT_INTEG_NUM], MYREAL (*Kpt)[GRT_SRC_M_NUM][GRT_INTEG_NUM][GRT_PTAM_PT_MAX], - MYCOMPLEX (*Fpt)[GRT_SRC_M_NUM][GRT_INTEG_NUM][GRT_PTAM_PT_MAX], MYINT (*Ipt)[GRT_SRC_M_NUM][GRT_INTEG_NUM], MYINT (*Gpt)[GRT_SRC_M_NUM][GRT_INTEG_NUM], bool *iendk0) + MYINT ir, MYINT im, MYINT v, real_t k, real_t dk, + cplx_t (*J3)[GRT_PTAM_WINDOW_SIZE][GRT_SRC_M_NUM][GRT_INTEG_NUM], real_t (*Kpt)[GRT_SRC_M_NUM][GRT_INTEG_NUM][GRT_PTAM_PT_MAX], + cplx_t (*Fpt)[GRT_SRC_M_NUM][GRT_INTEG_NUM][GRT_PTAM_PT_MAX], MYINT (*Ipt)[GRT_SRC_M_NUM][GRT_INTEG_NUM], MYINT (*Gpt)[GRT_SRC_M_NUM][GRT_INTEG_NUM], bool *iendk0) { - MYCOMPLEX tmp0; + cplx_t tmp0; if (Gpt[ir][im][v] >= GRT_PTAM_WINDOW_SIZE-1 && Ipt[ir][im][v] < GRT_PTAM_PT_MAX) { if (grt_cplx_peak_or_trough(im, v, J3[ir], k, dk, &Kpt[ir][im][v][Ipt[ir][im][v]], &tmp0) != 0) { Fpt[ir][im][v][Ipt[ir][im][v]++] = tmp0; @@ -81,11 +81,11 @@ static void process_peak_or_trough( * */ static void ptam_once( - const MYINT ir, const MYINT nr, const MYREAL precoef, MYREAL k, MYREAL dk, - MYCOMPLEX SUM3[nr][GRT_PTAM_WINDOW_SIZE][GRT_SRC_M_NUM][GRT_INTEG_NUM], - MYCOMPLEX sum_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], - MYREAL Kpt[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM][GRT_PTAM_PT_MAX], - MYCOMPLEX Fpt[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM][GRT_PTAM_PT_MAX], + const MYINT ir, const MYINT nr, const real_t precoef, real_t k, real_t dk, + cplx_t SUM3[nr][GRT_PTAM_WINDOW_SIZE][GRT_SRC_M_NUM][GRT_INTEG_NUM], + cplx_t sum_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], + real_t Kpt[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM][GRT_PTAM_PT_MAX], + cplx_t Fpt[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM][GRT_PTAM_PT_MAX], MYINT Ipt[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], MYINT Gpt[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], bool *iendk0) @@ -117,23 +117,23 @@ static void ptam_once( void grt_PTA_method( - GRT_MODEL1D *mod1d, MYREAL k0, MYREAL predk, - MYINT nr, MYREAL *rs, - MYCOMPLEX sum_J0[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], + GRT_MODEL1D *mod1d, real_t k0, real_t predk, + MYINT nr, real_t *rs, + cplx_t sum_J0[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], bool calc_upar, - MYCOMPLEX sum_uiz_J0[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], - MYCOMPLEX sum_uir_J0[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], + cplx_t sum_uiz_J0[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], + cplx_t sum_uir_J0[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], FILE *ptam_fstatsnr[nr][2], GRT_KernelFunc kerfunc) { // 需要兼容对正常收敛而不具有规律波峰波谷的序列 // 有时序列收敛比较好,不表现为规律的波峰波谷, // 此时设置最大等待次数,超过直接设置为中间值 - MYREAL k=0.0; + real_t k=0.0; // 不同震源不同阶数的核函数 F(k, w) - MYCOMPLEX QWV[GRT_SRC_M_NUM][GRT_QWV_NUM]; - MYCOMPLEX QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]; + cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM]; + cplx_t QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]; // 使用宏函数,方便定义 #define __CALLOC_ARRAY(VAR, TYP, __ARR) \ @@ -143,28 +143,28 @@ void grt_PTA_method( // 存储采样的值,维度3表示通过连续3个点来判断波峰或波谷 // 既用于存储被积函数,也最后用于存储求和的结果 #define __ARR [GRT_PTAM_WINDOW_SIZE][GRT_SRC_M_NUM][GRT_INTEG_NUM] - __CALLOC_ARRAY(SUM3, MYCOMPLEX, __ARR); - __CALLOC_ARRAY(SUM3_uiz, MYCOMPLEX, __ARR); - __CALLOC_ARRAY(SUM3_uir, MYCOMPLEX, __ARR); + __CALLOC_ARRAY(SUM3, cplx_t, __ARR); + __CALLOC_ARRAY(SUM3_uiz, cplx_t, __ARR); + __CALLOC_ARRAY(SUM3_uir, cplx_t, __ARR); #undef __ARR // 之前求和的值 #define __ARR [GRT_SRC_M_NUM][GRT_INTEG_NUM] - __CALLOC_ARRAY(sum_J, MYCOMPLEX, __ARR); - __CALLOC_ARRAY(sum_uiz_J, MYCOMPLEX, __ARR); - __CALLOC_ARRAY(sum_uir_J, MYCOMPLEX, __ARR); + __CALLOC_ARRAY(sum_J, cplx_t, __ARR); + __CALLOC_ARRAY(sum_uiz_J, cplx_t, __ARR); + __CALLOC_ARRAY(sum_uir_J, cplx_t, __ARR); #undef __ARR // 存储波峰波谷的位置和值 #define __ARR [GRT_SRC_M_NUM][GRT_INTEG_NUM][GRT_PTAM_PT_MAX] - __CALLOC_ARRAY(Kpt, MYREAL, __ARR); - __CALLOC_ARRAY(Fpt, MYCOMPLEX, __ARR); + __CALLOC_ARRAY(Kpt, real_t, __ARR); + __CALLOC_ARRAY(Fpt, cplx_t, __ARR); - __CALLOC_ARRAY(Kpt_uiz, MYREAL, __ARR); - __CALLOC_ARRAY(Fpt_uiz, MYCOMPLEX, __ARR); + __CALLOC_ARRAY(Kpt_uiz, real_t, __ARR); + __CALLOC_ARRAY(Fpt_uiz, cplx_t, __ARR); - __CALLOC_ARRAY(Kpt_uir, MYREAL, __ARR); - __CALLOC_ARRAY(Fpt_uir, MYCOMPLEX, __ARR); + __CALLOC_ARRAY(Kpt_uir, real_t, __ARR); + __CALLOC_ARRAY(Fpt_uir, cplx_t, __ARR); #undef __ARR #define __ARR [GRT_SRC_M_NUM][GRT_INTEG_NUM] @@ -201,10 +201,10 @@ void grt_PTA_method( // 对于PTAM,不同震中距使用不同dk for(MYINT ir=0; ir1; --n){ for(MYINT i=0; ik3[2] - item_pt->k3[0]; - const MYCOMPLEX (*F3)[GRT_SRC_M_NUM][GRT_QWV_NUM] = (isuiz)? item_pt->F3_uiz : item_pt->F3; +static cplx_t simpson(const KInterval *item_pt, MYINT im, MYINT iqwv, bool isuiz, SIMPSON_INTV stats) { + cplx_t Fint = 0.0; + real_t klen = item_pt->k3[2] - item_pt->k3[0]; + const cplx_t (*F3)[GRT_SRC_M_NUM][GRT_QWV_NUM] = (isuiz)? item_pt->F3_uiz : item_pt->F3; // 使用F(k)*sqrt(k)来衡量积分值,这可以平衡后续计算F(k)*Jm(kr)k积分时的系数 - MYREAL sk[3]; + real_t sk[3]; for(MYINT i=0; i<3; ++i){ sk[i] = sqrt(item_pt->k3[i]); } @@ -136,8 +136,8 @@ static MYCOMPLEX simpson(const KInterval *item_pt, MYINT im, MYINT iqwv, bool is } /** 比较QWV的最大绝对值 */ -static void get_maxabsQWV(const MYCOMPLEX F[GRT_SRC_M_NUM][GRT_QWV_NUM], MYREAL maxabsF[GRT_GTYPES_MAX]){ - MYREAL tmp; +static void get_maxabsQWV(const cplx_t F[GRT_SRC_M_NUM][GRT_QWV_NUM], real_t maxabsF[GRT_GTYPES_MAX]){ + real_t tmp; for(MYINT i=0; iF3_uiz : ptKitvL->F3; - const MYCOMPLEX (*F3R)[GRT_SRC_M_NUM][GRT_QWV_NUM] = (isuiz)? ptKitvR->F3_uiz : ptKitvR->F3; + const cplx_t (*F3L)[GRT_SRC_M_NUM][GRT_QWV_NUM] = (isuiz)? ptKitvL->F3_uiz : ptKitvL->F3; + const cplx_t (*F3R)[GRT_SRC_M_NUM][GRT_QWV_NUM] = (isuiz)? ptKitvR->F3_uiz : ptKitvR->F3; // 取近似积分 \int_k1^k2 k^0.5 dk - MYREAL kcoef13 = RTWOTHIRD*( ptKitv->k3[2]*sqrt(ptKitv->k3[2]) - ptKitv->k3[0]*sqrt(ptKitv->k3[0]) ); - MYREAL kcoef12 = RTWOTHIRD*( ptKitvL->k3[2]*sqrt(ptKitvL->k3[2]) - ptKitvL->k3[0]*sqrt(ptKitvL->k3[0]) ); - MYREAL kcoef23 = RTWOTHIRD*( ptKitvR->k3[2]*sqrt(ptKitvR->k3[2]) - ptKitvR->k3[0]*sqrt(ptKitvR->k3[0]) ); + real_t kcoef13 = RTWOTHIRD*( ptKitv->k3[2]*sqrt(ptKitv->k3[2]) - ptKitv->k3[0]*sqrt(ptKitv->k3[0]) ); + real_t kcoef12 = RTWOTHIRD*( ptKitvL->k3[2]*sqrt(ptKitvL->k3[2]) - ptKitvL->k3[0]*sqrt(ptKitvL->k3[0]) ); + real_t kcoef23 = RTWOTHIRD*( ptKitvR->k3[2]*sqrt(ptKitvR->k3[2]) - ptKitvR->k3[0]*sqrt(ptKitvR->k3[0]) ); - MYREAL S_dif, S_ref; + real_t S_dif, S_ref; bool badtol = false; for(MYINT im=0; imkcmpnm, sizeof(pt_hd->kcmpnm), "%s%s%c", s_prefix, srcname, ch); @@ -228,7 +228,7 @@ static void write_one_to_sac( grt_reset_fftw_holder_zero(pt_fh); // 赋值复数,包括时移 - MYCOMPLEX cfac, ccoef; + cplx_t cfac, ccoef; cfac = exp(I*PI2*pt_fh->df*delayT); ccoef = sgn; // 只赋值有效长度,其余均为0 @@ -294,15 +294,15 @@ static void write_one_to_sac( */ static void single_freq2time_write_to_file( const char *command, const GRT_MODEL1D *mod1d, 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, + const real_t wI, GRT_FFTW_HOLDER *pt_fh, + const char *s_dist, const real_t dist, + const real_t depsrc, const real_t deprcv, + const real_t delayT0, const real_t delayV0, const bool calc_upar, const bool doEX, const bool doVF, const bool doHF, const bool doDC, SACHEAD *pt_hd, const char *chalst, - MYCOMPLEX *grn[GRT_SRC_M_NUM][GRT_CHANNEL_NUM], - MYCOMPLEX *grn_uiz[GRT_SRC_M_NUM][GRT_CHANNEL_NUM], - MYCOMPLEX *grn_uir[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]) + cplx_t *grn[GRT_SRC_M_NUM][GRT_CHANNEL_NUM], + cplx_t *grn_uiz[GRT_SRC_M_NUM][GRT_CHANNEL_NUM], + cplx_t *grn_uir[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]) { // 文件保存子目录 char *s_output_subdir = NULL; @@ -311,7 +311,7 @@ static void single_freq2time_write_to_file( GRTCheckMakeDir(command, s_output_subdir); // 时间延迟 - MYREAL delayT = delayT0; + real_t delayT = delayT0; if(delayV0 > 0.0) delayT += sqrt( GRT_SQUARE(dist) + GRT_SQUARE(deprcv - depsrc) ) / delayV0; // 修改SAC头段时间变量 pt_hd->b = delayT; @@ -371,15 +371,15 @@ static void single_freq2time_write_to_file( void grt_GF_freq2time_write_to_file( const char *command, const GRT_MODEL1D *mod1d, const char *s_output_dir, const char *s_modelname, const char *s_depsrc, const char *s_deprcv, - 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, + const real_t wI, GRT_FFTW_HOLDER *pt_fh, + const MYINT nr, char *s_dists[nr], const real_t dists[nr], real_t travtPS[nr][2], + const real_t depsrc, const real_t deprcv, + const real_t delayT0, const real_t delayV0, const bool calc_upar, const bool doEX, const bool doVF, const bool doHF, const bool doDC, const char *chalst, - MYCOMPLEX *grn[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM], - MYCOMPLEX *grn_uiz[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM], - MYCOMPLEX *grn_uir[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM]) + cplx_t *grn[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM], + cplx_t *grn_uiz[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM], + cplx_t *grn_uir[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM]) { // 建立SAC头文件,包含必要的头变量 SACHEAD hd = new_sac_head(pt_fh->dt, pt_fh->nt, delayT0); diff --git a/pygrt/C_extension/src/dynamic/grn.c b/pygrt/C_extension/src/dynamic/grn.c index f36e27c0..7276659d 100755 --- a/pygrt/C_extension/src/dynamic/grn.c +++ b/pygrt/C_extension/src/dynamic/grn.c @@ -44,12 +44,12 @@ * @param[out] grn 三分量频谱 */ static void recordin_GRN( - MYINT iw, MYINT nr, MYCOMPLEX coef, MYCOMPLEX sum_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], - MYCOMPLEX *grn[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM] + MYINT iw, MYINT nr, cplx_t coef, cplx_t sum_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM], + cplx_t *grn[nr][GRT_SRC_M_NUM][GRT_CHANNEL_NUM] ) { // 局部变量,将某个频点的格林函数谱临时存放 - MYCOMPLEX (*tmp_grn)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM] = (MYCOMPLEX(*)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM])calloc(nr, sizeof(*tmp_grn)); + cplx_t (*tmp_grn)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM] = (cplx_t(*)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM])calloc(nr, sizeof(*tmp_grn)); for(MYINT ir=0; irRho[mod1d->isrc]; // 震源区密度 - const MYREAL fac = 1.0/(4.0*PI*Rho); - const MYREAL hs = GRT_MAX(fabs(mod1d->depsrc - mod1d->deprcv), GRT_MIN_DEPTH_GAP_SRC_RCV); // hs=max(震源和台站深度差,1.0) + const real_t Rho = mod1d->Rho[mod1d->isrc]; // 震源区密度 + const real_t fac = 1.0/(4.0*PI*Rho); + const real_t hs = GRT_MAX(fabs(mod1d->depsrc - mod1d->deprcv), GRT_MIN_DEPTH_GAP_SRC_RCV); // hs=max(震源和台站深度差,1.0) // 乘相应系数 k0 *= PI/hs; - const MYREAL k02 = k0*k0; - const MYREAL ampk2 = ampk*ampk; + const real_t k02 = k0*k0; + const real_t ampk2 = ampk*ampk; if(vmin_ref < 0.0) keps = 0.0; // 若使用峰谷平均法,则不使用keps进行收敛判断 bool useFIM = (filonLength > 0.0) || (safilonTol > 0.0) ; // 是否使用Filon积分(包括自适应Filon) - const MYREAL dk=PI2/(Length*rmax); // 波数积分间隔 - const MYREAL filondk = (filonLength > 0.0) ? PI2/(filonLength*rmax) : 0.0; // Filon积分间隔 - const MYREAL filonK = filonCut/rmax; // 波数积分和Filon积分的分割点 + const real_t dk=PI2/(Length*rmax); // 波数积分间隔 + const real_t filondk = (filonLength > 0.0) ? PI2/(filonLength*rmax) : 0.0; // Filon积分间隔 + const real_t filonK = filonCut/rmax; // 波数积分和Filon积分的分割点 // PTAM的积分中间结果, 每个震中距两个文件,因为PTAM对不同震中距使用不同的dk @@ -139,18 +139,18 @@ void grt_integ_grn_spec( // schedule语句可以动态调度任务,最大程度地使用计算资源 #pragma omp parallel for schedule(guided) default(shared) for(MYINT iw=nf1; iw<=nf2; ++iw){ - MYREAL k=0.0; // 波数 - MYREAL w = freqs[iw]*PI2; // 实频率 - MYCOMPLEX omega = w - wI*I; // 复数频率 omega = w - i*wI - MYCOMPLEX omega2_inv = 1.0/omega; // 1/omega^2 + real_t k=0.0; // 波数 + real_t w = freqs[iw]*PI2; // 实频率 + cplx_t omega = w - wI*I; // 复数频率 omega = w - i*wI + cplx_t omega2_inv = 1.0/omega; // 1/omega^2 omega2_inv = omega2_inv*omega2_inv; - MYCOMPLEX coef = -dk*fac*omega2_inv; // 最终要乘上的系数 + cplx_t coef = -dk*fac*omega2_inv; // 最终要乘上的系数 // 局部变量,用于求和 sum F(ki,w)Jm(ki*r)ki // 关于形状详见int_Pk()函数内的注释 - MYCOMPLEX (*sum_J)[GRT_SRC_M_NUM][GRT_INTEG_NUM] = (MYCOMPLEX(*)[GRT_SRC_M_NUM][GRT_INTEG_NUM])calloc(nr, sizeof(*sum_J)); - MYCOMPLEX (*sum_uiz_J)[GRT_SRC_M_NUM][GRT_INTEG_NUM] = (calc_upar)? (MYCOMPLEX(*)[GRT_SRC_M_NUM][GRT_INTEG_NUM])calloc(nr, sizeof(*sum_uiz_J)) : NULL; - MYCOMPLEX (*sum_uir_J)[GRT_SRC_M_NUM][GRT_INTEG_NUM] = (calc_upar)? (MYCOMPLEX(*)[GRT_SRC_M_NUM][GRT_INTEG_NUM])calloc(nr, sizeof(*sum_uir_J)) : NULL; + cplx_t (*sum_J)[GRT_SRC_M_NUM][GRT_INTEG_NUM] = (cplx_t(*)[GRT_SRC_M_NUM][GRT_INTEG_NUM])calloc(nr, sizeof(*sum_J)); + cplx_t (*sum_uiz_J)[GRT_SRC_M_NUM][GRT_INTEG_NUM] = (calc_upar)? (cplx_t(*)[GRT_SRC_M_NUM][GRT_INTEG_NUM])calloc(nr, sizeof(*sum_uiz_J)) : NULL; + cplx_t (*sum_uir_J)[GRT_SRC_M_NUM][GRT_INTEG_NUM] = (calc_upar)? (cplx_t(*)[GRT_SRC_M_NUM][GRT_INTEG_NUM])calloc(nr, sizeof(*sum_uir_J)) : NULL; GRT_MODEL1D *local_mod1d = NULL; #ifdef _OPENMP @@ -204,7 +204,7 @@ void grt_integ_grn_spec( - MYREAL kmax; + real_t kmax; // vmin_ref的正负性在这里不影响 kmax = sqrt(k02 + ampk2*(w/vmin_ref)*(w/vmin_ref)); diff --git a/pygrt/C_extension/src/dynamic/grt_greenfn.c b/pygrt/C_extension/src/dynamic/grt_greenfn.c index 5ec1ae98..bafbce07 100644 --- a/pygrt/C_extension/src/dynamic/grt_greenfn.c +++ b/pygrt/C_extension/src/dynamic/grt_greenfn.c @@ -49,8 +49,8 @@ typedef struct { /** 震源和接收器深度 */ struct { bool active; - MYREAL depsrc; - MYREAL deprcv; + real_t depsrc; + real_t deprcv; char *s_depsrc; char *s_deprcv; } D; @@ -59,12 +59,12 @@ typedef struct { bool active; MYINT nt; MYINT nf; - MYREAL dt; - MYREAL df; - MYREAL winT; ///< 时窗长度 - MYREAL zeta; ///< 虚频率系数, w <- w - zeta*PI/r* 1j - MYREAL wI; ///< 虚频率 zeta*PI/r - MYREAL *freqs; + real_t dt; + real_t df; + real_t winT; ///< 时窗长度 + real_t zeta; ///< 虚频率系数, w <- w - zeta*PI/r* 1j + real_t wI; ///< 虚频率 zeta*PI/r + real_t *freqs; MYINT upsample_n; ///< 升采样倍数 } N; /** 输出目录 */ @@ -75,8 +75,8 @@ typedef struct { /** 频段 */ struct { bool active; - MYREAL freq1; - MYREAL freq2; + real_t freq1; + real_t freq2; MYINT nf1; MYINT nf2; } H; @@ -84,25 +84,25 @@ typedef struct { struct { bool active; MYINT method; - MYREAL Length; - MYREAL filonLength; - MYREAL safilonTol; - MYREAL filonCut; + real_t Length; + real_t filonLength; + real_t safilonTol; + real_t filonCut; } L; /** 波数积分上限 */ struct { bool active; - MYREAL keps; - MYREAL ampk; - MYREAL k0; - MYREAL vmin; + real_t keps; + real_t ampk; + real_t k0; + real_t vmin; bool v_active; } K; /** 时间延迟 */ struct { bool active; - MYREAL delayT0; - MYREAL delayV0; + real_t delayT0; + real_t delayV0; } E; /** 波数积分过程的核函数文件 */ struct { @@ -118,7 +118,7 @@ typedef struct { bool active; char *s_raw; char **s_rs; - MYREAL *rs; + real_t *rs; MYINT nr; } R; /** 多线程 */ @@ -663,7 +663,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ fclose(fp); } // 转为浮点数 - Ctrl->R.rs = (MYREAL*)realloc(Ctrl->R.rs, sizeof(MYREAL)*(Ctrl->R.nr)); + Ctrl->R.rs = (real_t*)realloc(Ctrl->R.rs, sizeof(real_t)*(Ctrl->R.nr)); for(MYINT i=0; iR.nr; ++i){ Ctrl->R.rs[i] = atof(Ctrl->R.s_rs[i]); if(Ctrl->R.rs[i] < 0.0){ @@ -785,7 +785,7 @@ int greenfn_main(int argc, char **argv) { } // 最大最小速度 - MYREAL vmin, vmax; + real_t vmin, vmax; grt_get_mod1d_vmin_vmax(mod1d, &vmin, &vmax); // 参考最小速度 @@ -802,10 +802,10 @@ int greenfn_main(int argc, char **argv) { Ctrl->N.winT = Ctrl->N.nt*Ctrl->N.dt; // 最大震中距 - MYREAL rmax = Ctrl->R.rs[grt_findMax_MYREAL(Ctrl->R.rs, Ctrl->R.nr)]; + real_t rmax = Ctrl->R.rs[grt_findMax_real_t(Ctrl->R.rs, Ctrl->R.nr)]; // 时窗最大截止时刻 - MYREAL tmax = Ctrl->E.delayT0 + Ctrl->N.winT; + real_t tmax = Ctrl->E.delayT0 + Ctrl->N.winT; if(Ctrl->E.delayV0 > 0.0) tmax += rmax/Ctrl->E.delayV0; // 自动选择积分间隔,默认使用传统离散波数积分 @@ -824,7 +824,7 @@ int greenfn_main(int argc, char **argv) { // 定义要计算的频率、时窗等 Ctrl->N.nf = Ctrl->N.nt/2 + 1; Ctrl->N.df = 1.0/Ctrl->N.winT; - Ctrl->N.freqs = (MYREAL*)malloc(Ctrl->N.nf*sizeof(MYREAL)); + Ctrl->N.freqs = (real_t*)malloc(Ctrl->N.nf*sizeof(real_t)); for(int i=0; iN.nf; ++i){ Ctrl->N.freqs[i] = i*Ctrl->N.df; } @@ -851,16 +851,16 @@ int greenfn_main(int argc, char **argv) { } // 建立格林函数的complex数组 - MYCOMPLEX *(*grn)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM] = (MYCOMPLEX*(*)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]) calloc(Ctrl->R.nr, sizeof(*grn)); - MYCOMPLEX *(*grn_uiz)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM] = (Ctrl->e.active)? (MYCOMPLEX*(*)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]) calloc(Ctrl->R.nr, sizeof(*grn_uiz)) : NULL; - MYCOMPLEX *(*grn_uir)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM] = (Ctrl->e.active)? (MYCOMPLEX*(*)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]) calloc(Ctrl->R.nr, sizeof(*grn_uir)) : NULL; + cplx_t *(*grn)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM] = (cplx_t*(*)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]) calloc(Ctrl->R.nr, sizeof(*grn)); + cplx_t *(*grn_uiz)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM] = (Ctrl->e.active)? (cplx_t*(*)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]) calloc(Ctrl->R.nr, sizeof(*grn_uiz)) : NULL; + cplx_t *(*grn_uir)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM] = (Ctrl->e.active)? (cplx_t*(*)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]) calloc(Ctrl->R.nr, sizeof(*grn_uir)) : NULL; for(int ir=0; irR.nr; ++ir){ for(int i=0; iN.nf, sizeof(MYCOMPLEX)); - if(grn_uiz) grn_uiz[ir][i][c] = (MYCOMPLEX*)calloc(Ctrl->N.nf, sizeof(MYCOMPLEX)); - if(grn_uir) grn_uir[ir][i][c] = (MYCOMPLEX*)calloc(Ctrl->N.nf, sizeof(MYCOMPLEX)); + grn[ir][i][c] = (cplx_t*)calloc(Ctrl->N.nf, sizeof(cplx_t)); + if(grn_uiz) grn_uiz[ir][i][c] = (cplx_t*)calloc(Ctrl->N.nf, sizeof(cplx_t)); + if(grn_uir) grn_uir[ir][i][c] = (cplx_t*)calloc(Ctrl->N.nf, sizeof(cplx_t)); } } } @@ -887,7 +887,7 @@ int greenfn_main(int argc, char **argv) { 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); + real_t (* travtPS)[2] = (real_t (*)[2])calloc(Ctrl->R.nr, sizeof(real_t)*2); grt_GF_freq2time_write_to_file( command, mod1d, Ctrl->O.s_output_dir, Ctrl->M.s_modelname, Ctrl->D.s_depsrc, Ctrl->D.s_deprcv, diff --git a/pygrt/C_extension/src/dynamic/grt_syn.c b/pygrt/C_extension/src/dynamic/grt_syn.c index 29087044..0c1f9211 100644 --- a/pygrt/C_extension/src/dynamic/grt_syn.c +++ b/pygrt/C_extension/src/dynamic/grt_syn.c @@ -37,9 +37,9 @@ typedef struct { /** 方位角 */ struct { bool active; - MYREAL azimuth; - MYREAL azrad; - MYREAL backazimuth; + real_t azimuth; + real_t azrad; + real_t backazimuth; } A; /** 旋转到 Z, N, E */ struct { @@ -49,8 +49,8 @@ typedef struct { struct { bool active; bool mult_src_mu; - MYREAL M0; - MYREAL src_mu; + real_t M0; + real_t src_mu; } S; /** 剪切源 */ struct { @@ -92,10 +92,10 @@ typedef struct { } e; // 存储不同震源的震源机制相关参数的数组 - MYREAL mchn[GRT_MECHANISM_NUM]; + real_t mchn[GRT_MECHANISM_NUM]; // 方向因子数组 - MYREAL srcRadi[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]; + real_t srcRadi[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]; // 最终要计算的震源类型 MYINT computeType; diff --git a/pygrt/C_extension/src/dynamic/layer.c b/pygrt/C_extension/src/dynamic/layer.c index 7b5e90b6..1841ea45 100644 --- a/pygrt/C_extension/src/dynamic/layer.c +++ b/pygrt/C_extension/src/dynamic/layer.c @@ -31,15 +31,15 @@ void grt_topfree_RU(const GRT_MODEL1D *mod1d, RT_MATRIX *M) { - MYCOMPLEX cbcb0 = mod1d->cbcb[0]; + cplx_t cbcb0 = mod1d->cbcb[0]; if(cbcb0 != 0.0){ // 固体表面 // 公式(5.3.10-14) - MYCOMPLEX Delta = 0.0; - MYCOMPLEX cbcb02 = 0.25*cbcb0*cbcb0; + cplx_t Delta = 0.0; + cplx_t cbcb02 = 0.25*cbcb0*cbcb0; - MYCOMPLEX xa0 = mod1d->xa[0]; - MYCOMPLEX xb0 = mod1d->xb[0]; + cplx_t xa0 = mod1d->xa[0]; + cplx_t xb0 = mod1d->xb[0]; // 对公式(5.3.10-14)进行重新整理,对浮点数友好一些 Delta = -1.0 + xa0*xb0 + cbcb0 - cbcb02; @@ -63,14 +63,14 @@ void grt_topfree_RU(const GRT_MODEL1D *mod1d, RT_MATRIX *M) } -void grt_wave2qwv_REV_PSV(const GRT_MODEL1D *mod1d, MYCOMPLEX R[2][2], MYCOMPLEX R_EV[2][2]) +void grt_wave2qwv_REV_PSV(const GRT_MODEL1D *mod1d, cplx_t R[2][2], cplx_t R_EV[2][2]) { MYINT ircv = mod1d->ircv; - MYCOMPLEX xa = mod1d->xa[ircv]; - MYCOMPLEX xb = mod1d->xb[ircv]; - MYREAL k = mod1d->k; + cplx_t xa = mod1d->xa[ircv]; + cplx_t xb = mod1d->xb[ircv]; + real_t k = mod1d->k; - MYCOMPLEX D11[2][2], D12[2][2]; + cplx_t D11[2][2], D12[2][2]; if(xb != 1.0){ // 位于固体层 // 公式(5.2.19) @@ -97,11 +97,11 @@ void grt_wave2qwv_REV_PSV(const GRT_MODEL1D *mod1d, MYCOMPLEX R[2][2], MYCOMPLEX } -void grt_wave2qwv_REV_SH(const GRT_MODEL1D *mod1d, MYCOMPLEX RL, MYCOMPLEX *R_EVL) +void grt_wave2qwv_REV_SH(const GRT_MODEL1D *mod1d, cplx_t RL, cplx_t *R_EVL) { MYINT ircv = mod1d->ircv; - MYCOMPLEX xb = mod1d->xb[ircv]; - MYREAL k = mod1d->k; + cplx_t xb = mod1d->xb[ircv]; + real_t k = mod1d->k; if(xb != 1.0){ // 位于固体层 @@ -114,22 +114,22 @@ void grt_wave2qwv_REV_SH(const GRT_MODEL1D *mod1d, MYCOMPLEX RL, MYCOMPLEX *R_EV } -void grt_wave2qwv_z_REV_PSV(const GRT_MODEL1D *mod1d, const MYCOMPLEX R[2][2], MYCOMPLEX R_EV[2][2]) +void grt_wave2qwv_z_REV_PSV(const GRT_MODEL1D *mod1d, const cplx_t R[2][2], cplx_t R_EV[2][2]) { MYINT ircv = mod1d->ircv; - MYCOMPLEX xa = mod1d->xa[ircv]; - MYCOMPLEX xb = mod1d->xb[ircv]; - MYREAL k = mod1d->k; + cplx_t xa = mod1d->xa[ircv]; + cplx_t xb = mod1d->xb[ircv]; + real_t k = mod1d->k; // 将垂直波函数转为ui,z在(B_m, P_m, C_m)系下的分量 // 新推导的公式 - MYCOMPLEX ak = k*k*xa; - MYCOMPLEX bk = k*k*xb; - MYCOMPLEX bb = xb*bk; - MYCOMPLEX aa = xa*ak; - MYCOMPLEX D11[2][2] = {{ak, bb}, {aa, bk}}; - MYCOMPLEX D12[2][2] = {{-ak, bb}, {aa, -bk}}; + cplx_t ak = k*k*xa; + cplx_t bk = k*k*xb; + cplx_t bb = xb*bk; + cplx_t aa = xa*ak; + cplx_t D11[2][2] = {{ak, bb}, {aa, bk}}; + cplx_t D12[2][2] = {{-ak, bb}, {aa, -bk}}; // 位于液体层 if(xb == 1.0){ @@ -147,15 +147,15 @@ void grt_wave2qwv_z_REV_PSV(const GRT_MODEL1D *mod1d, const MYCOMPLEX R[2][2], M } -void grt_wave2qwv_z_REV_SH(const GRT_MODEL1D *mod1d, MYCOMPLEX RL, MYCOMPLEX *R_EVL) +void grt_wave2qwv_z_REV_SH(const GRT_MODEL1D *mod1d, cplx_t RL, cplx_t *R_EVL) { MYINT ircv = mod1d->ircv; - MYCOMPLEX xb = mod1d->xb[ircv]; - MYREAL k = mod1d->k; + cplx_t xb = mod1d->xb[ircv]; + real_t k = mod1d->k; // 将垂直波函数转为ui,z在(B_m, P_m, C_m)系下的分量 // 新推导的公式 - MYCOMPLEX bk = k*k*xb; + cplx_t bk = k*k*xb; if(xb != 1.0){ // 位于固体层 @@ -173,10 +173,10 @@ void grt_wave2qwv_z_REV_SH(const GRT_MODEL1D *mod1d, MYCOMPLEX RL, MYCOMPLEX *R_ void grt_RT_matrix_ll_PSV(const GRT_MODEL1D *mod1d, MYINT iy, RT_MATRIX *M) { - MODEL_2LAYS_ATTRIB(MYCOMPLEX, xa); - MODEL_2LAYS_ATTRIB(MYREAL, Rho); + MODEL_2LAYS_ATTRIB(cplx_t, xa); + MODEL_2LAYS_ATTRIB(real_t, Rho); - MYCOMPLEX A = xa1*Rho2 + xa2*Rho1; + cplx_t A = xa1*Rho2 + xa2*Rho1; if(A==0.0){ M->stats = GRT_INVERSE_FAILURE; return; @@ -208,11 +208,11 @@ void grt_RT_matrix_ll_SH(RT_MATRIX *M) void grt_RT_matrix_ls_PSV(const GRT_MODEL1D *mod1d, const MYINT iy, RT_MATRIX *M) { - MODEL_2LAYS_ATTRIB(MYCOMPLEX, xa); - MODEL_2LAYS_ATTRIB(MYCOMPLEX, xb); - MODEL_2LAYS_ATTRIB(MYCOMPLEX, mu); - MODEL_2LAYS_ATTRIB(MYCOMPLEX, cbcb); - MODEL_2LAYS_ATTRIB(MYREAL, Rho); + MODEL_2LAYS_ATTRIB(cplx_t, xa); + MODEL_2LAYS_ATTRIB(cplx_t, xb); + MODEL_2LAYS_ATTRIB(cplx_t, mu); + MODEL_2LAYS_ATTRIB(cplx_t, cbcb); + MODEL_2LAYS_ATTRIB(real_t, Rho); // 后缀1表示上层的液体的物理参数,后缀2表示下层的固体的物理参数 // 若mu2==0, 则下层为液体,参数需相互交换 @@ -225,32 +225,32 @@ void grt_RT_matrix_ls_PSV(const GRT_MODEL1D *mod1d, const MYINT iy, RT_MATRIX *M } // 使用指针 - MYCOMPLEX (*pRD)[2], (*pRU)[2]; - MYCOMPLEX (*pTD)[2], (*pTU)[2]; + cplx_t (*pRD)[2], (*pRU)[2]; + cplx_t (*pTD)[2], (*pTU)[2]; if(isfluidUp){ pRD = M->RD; pRU = M->RU; pTD = M->TD; pTU = M->TU; } else { pRD = M->RU; pRU = M->RD; pTD = M->TU; pTU = M->TD; - GRT_SWAP(MYREAL, Rho1, Rho2); - GRT_SWAP(MYCOMPLEX, xa1, xa2); - GRT_SWAP(MYCOMPLEX, xb1, xb2); - GRT_SWAP(MYCOMPLEX, cbcb1, cbcb2); - GRT_SWAP(MYCOMPLEX, mu1, mu2); + GRT_SWAP(real_t, Rho1, Rho2); + GRT_SWAP(cplx_t, xa1, xa2); + GRT_SWAP(cplx_t, xb1, xb2); + GRT_SWAP(cplx_t, cbcb1, cbcb2); + GRT_SWAP(cplx_t, mu1, mu2); sgn = -1; } // 定义一些中间变量来简化运算和书写 - MYCOMPLEX lamka1k = Rho1*GRT_SQUARE(mod1d->c_phase); - MYCOMPLEX kb2k = cbcb2; - MYCOMPLEX Og2k = 1.0 - 0.5*kb2k; - MYCOMPLEX Og2k2 = Og2k*Og2k; - MYCOMPLEX A = 2.0*Og2k2*xa1*mu2 + 0.5*lamka1k*kb2k*xa2 - 2.0*mu2*xa1*xa2*xb2; - MYCOMPLEX B = 2.0*Og2k2*xa1*mu2 - 0.5*lamka1k*kb2k*xa2 + 2.0*mu2*xa1*xa2*xb2; - MYCOMPLEX C = 2.0*Og2k2*xa1*mu2 + 0.5*lamka1k*kb2k*xa2 + 2.0*mu2*xa1*xa2*xb2; - MYCOMPLEX D = 2.0*Og2k2*xa1*mu2 - 0.5*lamka1k*kb2k*xa2 - 2.0*mu2*xa1*xa2*xb2; + cplx_t lamka1k = Rho1*GRT_SQUARE(mod1d->c_phase); + cplx_t kb2k = cbcb2; + cplx_t Og2k = 1.0 - 0.5*kb2k; + cplx_t Og2k2 = Og2k*Og2k; + cplx_t A = 2.0*Og2k2*xa1*mu2 + 0.5*lamka1k*kb2k*xa2 - 2.0*mu2*xa1*xa2*xb2; + cplx_t B = 2.0*Og2k2*xa1*mu2 - 0.5*lamka1k*kb2k*xa2 + 2.0*mu2*xa1*xa2*xb2; + cplx_t C = 2.0*Og2k2*xa1*mu2 + 0.5*lamka1k*kb2k*xa2 + 2.0*mu2*xa1*xa2*xb2; + cplx_t D = 2.0*Og2k2*xa1*mu2 - 0.5*lamka1k*kb2k*xa2 - 2.0*mu2*xa1*xa2*xb2; if(A == 0.0){ M->stats = GRT_INVERSE_FAILURE; @@ -279,7 +279,7 @@ void grt_RT_matrix_ls_SH(const GRT_MODEL1D *mod1d, const MYINT iy, RT_MATRIX *M) { // TEMPORARY!!!!!! // 之后不再使用mu或其它变量来判断是否为液体,直接定义一个新的数组 - MODEL_2LAYS_ATTRIB(MYCOMPLEX, mu); + MODEL_2LAYS_ATTRIB(cplx_t, mu); // 后缀1表示上层的液体的物理参数,后缀2表示下层的固体的物理参数 // 若mu2==0, 则下层为液体,参数需相互交换 @@ -291,8 +291,8 @@ void grt_RT_matrix_ls_SH(const GRT_MODEL1D *mod1d, const MYINT iy, RT_MATRIX *M) } // 使用指针 - MYCOMPLEX *pRDL, *pRUL; - MYCOMPLEX *pTDL, *pTUL; + cplx_t *pRDL, *pRUL; + cplx_t *pTDL, *pTUL; if(isfluidUp){ pRDL = &M->RDL; pRUL = &M->RUL; pTDL = &M->TDL; pTUL = &M->TUL; @@ -311,23 +311,23 @@ void grt_RT_matrix_ls_SH(const GRT_MODEL1D *mod1d, const MYINT iy, RT_MATRIX *M) void grt_RT_matrix_ss_PSV(const GRT_MODEL1D *mod1d, const MYINT iy, RT_MATRIX *M) { - MODEL_2LAYS_ATTRIB(MYCOMPLEX, xa); - MODEL_2LAYS_ATTRIB(MYCOMPLEX, xb); - MODEL_2LAYS_ATTRIB(MYCOMPLEX, mu); - MODEL_2LAYS_ATTRIB(MYCOMPLEX, cbcb); - MODEL_2LAYS_ATTRIB(MYREAL, Rho); + MODEL_2LAYS_ATTRIB(cplx_t, xa); + MODEL_2LAYS_ATTRIB(cplx_t, xb); + MODEL_2LAYS_ATTRIB(cplx_t, mu); + MODEL_2LAYS_ATTRIB(cplx_t, cbcb); + MODEL_2LAYS_ATTRIB(real_t, Rho); // 定义一些中间变量来简化运算和书写 - // MYREAL kk = k*k; - MYCOMPLEX rmu = mu1/mu2; - MYCOMPLEX dmu = rmu - 1.0; // mu1 - mu2; 分子分母同除mu2 - MYCOMPLEX dmu2 = dmu*dmu; + // real_t kk = k*k; + cplx_t rmu = mu1/mu2; + cplx_t dmu = rmu - 1.0; // mu1 - mu2; 分子分母同除mu2 + cplx_t dmu2 = dmu*dmu; - MYCOMPLEX mu1cbcb1 = rmu*cbcb1;// mu1*kb1_k2; - MYCOMPLEX mu2cbcb2 = cbcb2; // mu2*kb2_k2; + cplx_t mu1cbcb1 = rmu*cbcb1;// mu1*kb1_k2; + cplx_t mu2cbcb2 = cbcb2; // mu2*kb2_k2; - MYREAL rho12 = Rho1 / Rho2; - MYREAL rho21 = Rho2 / Rho1; + real_t rho12 = Rho1 / Rho2; + real_t rho21 = Rho2 / Rho1; // 从原公式上,分母包含5项,但前四项会随着k的增大迅速超过最后一项 // 最后一项要小前几项10余个数量级,但计算结果还是保持在最后一项的量级, @@ -335,7 +335,7 @@ void grt_RT_matrix_ss_PSV(const GRT_MODEL1D *mod1d, const MYINT iy, RT_MATRIX *M // 故会发生严重精度损失的情况。目前只在实部上观察到这个现象,虚部基本都在相近量级(或许是相对不明显) // // 以下对公式重新整理,提出k的高阶项,以避免上述问题 - MYCOMPLEX Delta; + cplx_t Delta; Delta = dmu2*(1.0-xa1*xb1)*(1.0-xa2*xb2) + mu1cbcb1*dmu*(rho21*(1.0-xa1*xb1) - (1.0-xa2*xb2)) + 0.25*mu1cbcb1*mu2cbcb2*(rho12*(1.0-xa2*xb2) + rho21*(1.0-xa1*xb1) - 2.0 - (xa1*xb2+xa2*xb1)); @@ -374,7 +374,7 @@ void grt_RT_matrix_ss_PSV(const GRT_MODEL1D *mod1d, const MYINT iy, RT_MATRIX *M - 0.25*mu1cbcb1*mu2cbcb2*(rho12*(1.0+xa2*xb2) + rho21*(1.0-xa1*xb1) - 2.0 + (xa1*xb2-xa2*xb1))) * Delta; // REFRACTION - MYCOMPLEX tmp; + cplx_t tmp; tmp = mu1cbcb1*xa1*(dmu*(xb2-xb1) - 0.5*mu1cbcb1*(rho21*xb1+xb2)) * Delta; M->TD[0][0] = tmp; M->TU[0][0] = (rho21*xa2/xa1) * tmp; tmp = mu1cbcb1*xb1*(dmu*(1.0-xa1*xb2) - 0.5*mu1cbcb1*(1.0-rho21)) * Delta; @@ -388,15 +388,15 @@ void grt_RT_matrix_ss_PSV(const GRT_MODEL1D *mod1d, const MYINT iy, RT_MATRIX *M void grt_RT_matrix_ss_SH(const GRT_MODEL1D *mod1d, const MYINT iy, RT_MATRIX *M) { - MODEL_2LAYS_ATTRIB(MYCOMPLEX, xb); - MODEL_2LAYS_ATTRIB(MYCOMPLEX, mu); + MODEL_2LAYS_ATTRIB(cplx_t, xb); + MODEL_2LAYS_ATTRIB(cplx_t, mu); // REFELCTION M->RUL = (mu2*xb2 - mu1*xb1) / (mu2*xb2 + mu1*xb1) ; M->RDL = - (M->RUL); // REFRACTION - MYCOMPLEX tmp; + cplx_t tmp; tmp = 2.0 / (mu2*xb2 + mu1*xb1); M->TDL = mu1*xb1 * tmp; M->TUL = mu2*xb2 * tmp; @@ -408,7 +408,7 @@ void grt_RT_matrix_PSV(const GRT_MODEL1D *mod1d, const MYINT iy, RT_MATRIX *M) { // TEMPORARY!!!!!! // 之后不再使用mu或其它变量来判断是否为液体,直接定义一个新的数组 - MODEL_2LAYS_ATTRIB(MYCOMPLEX, mu); + MODEL_2LAYS_ATTRIB(cplx_t, mu); // 根据界面两侧的具体情况选择函数 if(mu1 != 0.0 && mu2 != 0.0){ @@ -427,7 +427,7 @@ void grt_RT_matrix_SH(const GRT_MODEL1D *mod1d, const MYINT iy, RT_MATRIX *M) { // TEMPORARY!!!!!! // 之后不再使用mu或其它变量来判断是否为液体,直接定义一个新的数组 - MODEL_2LAYS_ATTRIB(MYCOMPLEX, mu); + MODEL_2LAYS_ATTRIB(cplx_t, mu); // 根据界面两侧的具体情况选择函数 if(mu1 != 0.0 && mu2 != 0.0){ @@ -444,12 +444,12 @@ void grt_RT_matrix_SH(const GRT_MODEL1D *mod1d, const MYINT iy, RT_MATRIX *M) void grt_delay_RT_matrix(const GRT_MODEL1D *mod1d, const MYINT iy, RT_MATRIX *M) { - MYREAL thk = mod1d->Thk[iy-1]; - MYCOMPLEX xa1 = mod1d->xa[iy-1]; - MYCOMPLEX xb1 = mod1d->xb[iy-1]; - MYREAL k = mod1d->k; + real_t thk = mod1d->Thk[iy-1]; + cplx_t xa1 = mod1d->xa[iy-1]; + cplx_t xb1 = mod1d->xb[iy-1]; + real_t k = mod1d->k; - MYCOMPLEX exa, exb, ex2a, ex2b, exab; + cplx_t exa, exb, ex2a, ex2b, exab; exa = exp(- k*thk*xa1); exb = exp(- k*thk*xb1); // exa = 0.9; @@ -474,11 +474,11 @@ void grt_delay_RT_matrix(const GRT_MODEL1D *mod1d, const MYINT iy, RT_MATRIX *M) void grt_get_layer_D( - MYCOMPLEX xa, MYCOMPLEX xb, MYCOMPLEX kbkb, MYCOMPLEX mu, - MYCOMPLEX omega, MYREAL rho, MYREAL k, MYCOMPLEX D[4][4], bool inverse, MYINT liquid_invtype) + cplx_t xa, cplx_t xb, cplx_t kbkb, cplx_t mu, + cplx_t omega, real_t rho, real_t k, cplx_t D[4][4], bool inverse, MYINT liquid_invtype) { // 第iy层物理量 - MYCOMPLEX Omg; + cplx_t Omg; if(xb != 1.0){ Omg = k*k - 0.5*kbkb; if( ! inverse ){ @@ -539,7 +539,7 @@ void grt_get_layer_D( } void grt_get_layer_D11( - MYCOMPLEX xa, MYCOMPLEX xb, MYREAL k, MYCOMPLEX D[2][2]) + cplx_t xa, cplx_t xb, real_t k, cplx_t D[2][2]) { // 第iy层物理量 if(xb != 1.0){ @@ -553,7 +553,7 @@ void grt_get_layer_D11( } void grt_get_layer_D12( - MYCOMPLEX xa, MYCOMPLEX xb, MYREAL k, MYCOMPLEX D[2][2]) + cplx_t xa, cplx_t xb, real_t k, cplx_t D[2][2]) { // 第iy层物理量 if(xb != 1.0){ @@ -567,11 +567,11 @@ void grt_get_layer_D12( } void grt_get_layer_D11_uiz( - MYCOMPLEX xa, MYCOMPLEX xb, MYREAL k, MYCOMPLEX D[2][2]) + cplx_t xa, cplx_t xb, real_t k, cplx_t D[2][2]) { // 第iy层物理量 - MYCOMPLEX a = k*xa; - MYCOMPLEX b = k*xb; + cplx_t a = k*xa; + cplx_t b = k*xb; if(xb != 1.0){ D[0][0] = a*k; D[0][1] = b*b; @@ -583,11 +583,11 @@ void grt_get_layer_D11_uiz( } void grt_get_layer_D12_uiz( - MYCOMPLEX xa, MYCOMPLEX xb, MYREAL k, MYCOMPLEX D[2][2]) + cplx_t xa, cplx_t xb, real_t k, cplx_t D[2][2]) { // 第iy层物理量 - MYCOMPLEX a = k*xa; - MYCOMPLEX b = k*xb; + cplx_t a = k*xa; + cplx_t b = k*xb; if(xb != 1.0){ D[0][0] = - a*k; D[0][1] = b*b; @@ -599,11 +599,11 @@ void grt_get_layer_D12_uiz( } void grt_get_layer_D21( - MYCOMPLEX xa, MYCOMPLEX xb, MYCOMPLEX kbkb, MYCOMPLEX mu, - MYCOMPLEX omega, MYREAL rho, MYREAL k, MYCOMPLEX D[2][2]) + cplx_t xa, cplx_t xb, cplx_t kbkb, cplx_t mu, + cplx_t omega, real_t rho, real_t k, cplx_t D[2][2]) { // 第iy层物理量 - MYCOMPLEX Omg; + cplx_t Omg; if(xb != 1.0){ Omg = k*k - 0.5*kbkb; D[0][0] = 2*mu*Omg; D[0][1] = 2*k*mu*k*xb; @@ -616,11 +616,11 @@ void grt_get_layer_D21( } void grt_get_layer_D22( - MYCOMPLEX xa, MYCOMPLEX xb, MYCOMPLEX kbkb, MYCOMPLEX mu, - MYCOMPLEX omega, MYREAL rho, MYREAL k, MYCOMPLEX D[2][2]) + cplx_t xa, cplx_t xb, cplx_t kbkb, cplx_t mu, + cplx_t omega, real_t rho, real_t k, cplx_t D[2][2]) { // 第iy层物理量 - MYCOMPLEX Omg; + cplx_t Omg; if(xb != 1.0){ Omg = k*k - 0.5*kbkb; D[0][0] = 2*mu*Omg; D[0][1] = -2*k*mu*k*xb; @@ -632,8 +632,8 @@ void grt_get_layer_D22( } void grt_get_layer_T( - MYCOMPLEX xb, MYCOMPLEX mu, - MYCOMPLEX omega, MYREAL k, MYCOMPLEX T[2][2], bool inverse) + cplx_t xb, cplx_t mu, + cplx_t omega, real_t k, cplx_t T[2][2], bool inverse) { // 液体层不应该使用该函数 if(xb == 1.0){ @@ -654,11 +654,11 @@ void grt_get_layer_T( } } -void grt_get_layer_E_Love(MYCOMPLEX xb1, MYREAL thk, MYREAL k, MYCOMPLEX E[2][2], bool inverse) +void grt_get_layer_E_Love(cplx_t xb1, real_t thk, real_t k, cplx_t E[2][2], bool inverse) { - MYCOMPLEX exb = exp(k*thk*xb1); + cplx_t exb = exp(k*thk*xb1); - memset(E, 0, sizeof(MYCOMPLEX) * 4); + memset(E, 0, sizeof(cplx_t) * 4); if(! inverse){ E[0][0] = exb; E[1][1] = 1.0/exb; @@ -670,14 +670,14 @@ void grt_get_layer_E_Love(MYCOMPLEX xb1, MYREAL thk, MYREAL k, MYCOMPLEX E[2][2] } void grt_get_layer_E_Rayl( - MYCOMPLEX xa1, MYCOMPLEX xb1, MYREAL thk, MYREAL k, MYCOMPLEX E[4][4], bool inverse) + cplx_t xa1, cplx_t xb1, real_t thk, real_t k, cplx_t E[4][4], bool inverse) { - MYCOMPLEX exa, exb; + cplx_t exa, exb; exa = exp(k*thk*xa1); exb = exp(k*thk*xb1); - memset(E, 0, sizeof(MYCOMPLEX) * 16); + memset(E, 0, sizeof(cplx_t) * 16); if( ! inverse){ E[0][0] = exa; @@ -693,26 +693,26 @@ void grt_get_layer_E_Rayl( } void grt_RT_matrix_from_4x4( - MYCOMPLEX xa1, MYCOMPLEX xb1, MYCOMPLEX kbkb1, MYCOMPLEX mu1, MYREAL rho1, - MYCOMPLEX xa2, MYCOMPLEX xb2, MYCOMPLEX kbkb2, MYCOMPLEX mu2, MYREAL rho2, - MYCOMPLEX omega, MYREAL thk, - MYREAL k, - 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) + cplx_t xa1, cplx_t xb1, cplx_t kbkb1, cplx_t mu1, real_t rho1, + cplx_t xa2, cplx_t xb2, cplx_t kbkb2, cplx_t mu2, real_t rho2, + cplx_t omega, real_t thk, + real_t k, + cplx_t RD[2][2], cplx_t *RDL, cplx_t RU[2][2], cplx_t *RUL, + cplx_t TD[2][2], cplx_t *TDL, cplx_t TU[2][2], cplx_t *TUL, MYINT *stats) { - MYCOMPLEX D1_inv[4][4], D2[4][4], Q[4][4]; + cplx_t D1_inv[4][4], D2[4][4], Q[4][4]; grt_get_layer_D(xa1, xb1, kbkb1, mu1, omega, rho1, k, D1_inv, true, 2); grt_get_layer_D(xa2, xb2, kbkb2, mu2, omega, rho2, k, D2, false, 2); grt_cmatmxn_mul(4, 4, 4, D1_inv, D2, Q); - MYCOMPLEX exa, exb; + cplx_t exa, exb; exa = exp(-k*thk*xa1); exb = exp(-k*thk*xb1); - MYCOMPLEX E[4][4] = {0}; + cplx_t E[4][4] = {0}; E[0][0] = exa; E[1][1] = exb; E[2][2] = 1/exa; @@ -720,7 +720,7 @@ void grt_RT_matrix_from_4x4( grt_cmatmxn_mul(4, 4, 4, E, Q, Q); // 对Q矩阵划分子矩阵 - MYCOMPLEX Q11[2][2], Q12[2][2], Q21[2][2], Q22[2][2]; + cplx_t Q11[2][2], Q12[2][2], Q21[2][2], Q22[2][2]; 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); diff --git a/pygrt/C_extension/src/dynamic/propagate.c b/pygrt/C_extension/src/dynamic/propagate.c index 8fa00adf..fa1540a3 100755 --- a/pygrt/C_extension/src/dynamic/propagate.c +++ b/pygrt/C_extension/src/dynamic/propagate.c @@ -28,8 +28,8 @@ void grt_kernel( - const GRT_MODEL1D *mod1d, MYCOMPLEX QWV[GRT_SRC_M_NUM][GRT_QWV_NUM], - bool calc_uiz, MYCOMPLEX QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]) + const GRT_MODEL1D *mod1d, cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM], + bool calc_uiz, cplx_t QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]) { // 初始化qwv为0 for(MYINT i=0; in; ++iy){ // 因为n>=3, 故一定会进入该循环 @@ -120,13 +120,13 @@ void grt_kernel( // 计算震源系数 - MYCOMPLEX src_coef_PSV[GRT_SRC_M_NUM][GRT_QWV_NUM-1][2] = {0}; - MYCOMPLEX src_coef_SH[GRT_SRC_M_NUM][2] = {0}; + cplx_t src_coef_PSV[GRT_SRC_M_NUM][GRT_QWV_NUM-1][2] = {0}; + cplx_t src_coef_SH[GRT_SRC_M_NUM][2] = {0}; grt_source_coef_PSV(mod1d, src_coef_PSV); grt_source_coef_SH(mod1d, src_coef_SH); // 临时中转矩阵 (temperary) - MYCOMPLEX tmpR2[2][2], tmp2x2[2][2], tmpRL, tmp2x2_uiz[2][2], tmpRL2; + cplx_t tmpR2[2][2], tmp2x2[2][2], tmpRL, tmp2x2_uiz[2][2], tmpRL2; // 递推RU_FA grt_topfree_RU(mod1d, M_top); @@ -238,14 +238,14 @@ void grt_kernel( void grt_construct_qwv( bool ircvup, - const MYCOMPLEX R1[2][2], MYCOMPLEX RL1, - const MYCOMPLEX R2[2][2], MYCOMPLEX RL2, - const MYCOMPLEX coef_PSV[GRT_QWV_NUM-1][2], const MYCOMPLEX coef_SH[2], - MYCOMPLEX qwv[GRT_QWV_NUM]) + const cplx_t R1[2][2], cplx_t RL1, + const cplx_t R2[2][2], cplx_t RL2, + const cplx_t coef_PSV[GRT_QWV_NUM-1][2], const cplx_t coef_SH[2], + cplx_t qwv[GRT_QWV_NUM]) { - MYCOMPLEX qw0[2], qw1[2], v0; - MYCOMPLEX coefD[2] = {coef_PSV[0][0], coef_PSV[1][0]}; - MYCOMPLEX coefU[2] = {coef_PSV[0][1], coef_PSV[1][1]}; + cplx_t qw0[2], qw1[2], v0; + cplx_t coefD[2] = {coef_PSV[0][0], coef_PSV[1][0]}; + cplx_t coefU[2] = {coef_PSV[0][1], coef_PSV[1][1]}; if(ircvup){ grt_cmat2x1_mul(R2, coefD, qw0); qw0[0] += coefU[0]; qw0[1] += coefU[1]; diff --git a/pygrt/C_extension/src/dynamic/source.c b/pygrt/C_extension/src/dynamic/source.c index 53cd4035..c0e311db 100755 --- a/pygrt/C_extension/src/dynamic/source.c +++ b/pygrt/C_extension/src/dynamic/source.c @@ -18,7 +18,7 @@ #include "grt/common/prtdbg.h" -void grt_source_coef_PSV(const GRT_MODEL1D *mod1d, MYCOMPLEX coef[GRT_SRC_M_NUM][GRT_QWV_NUM-1][2]) +void grt_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM-1][2]) { // 先全部赋0 for(MYINT i=0; iisrc; - MYCOMPLEX xa = mod1d->xa[isrc]; - MYCOMPLEX caca = mod1d->caca[isrc]; - MYCOMPLEX xb = mod1d->xb[isrc]; - MYCOMPLEX cbcb = mod1d->cbcb[isrc]; - MYREAL k = mod1d->k; + cplx_t xa = mod1d->xa[isrc]; + cplx_t caca = mod1d->caca[isrc]; + cplx_t xb = mod1d->xb[isrc]; + cplx_t cbcb = mod1d->cbcb[isrc]; + real_t k = mod1d->k; - MYCOMPLEX tmp; + cplx_t tmp; // 爆炸源, 通过(4.9.8)的矩张量源公式,提取各向同性的量(M11+M22+M33),-a+k^2/a -> ka^2/a @@ -67,7 +67,7 @@ void grt_source_coef_PSV(const GRT_MODEL1D *mod1d, MYCOMPLEX coef[GRT_SRC_M_NUM] } -void grt_source_coef_SH(const GRT_MODEL1D *mod1d, MYCOMPLEX coef[GRT_SRC_M_NUM][2]) +void grt_source_coef_SH(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][2]) { // 先全部赋0 for(MYINT i=0; iisrc; - MYCOMPLEX xb = mod1d->xb[isrc]; - MYCOMPLEX cbcb = mod1d->cbcb[isrc]; - MYREAL k = mod1d->k; + cplx_t xb = mod1d->xb[isrc]; + cplx_t cbcb = mod1d->cbcb[isrc]; + real_t k = mod1d->k; - // MYCOMPLEX b = k*src_xb; - MYCOMPLEX tmp; + // cplx_t b = k*src_xb; + cplx_t tmp; // 水平力源 (4.6.21,26), 这里可以把x1,x2方向的力转到r,theta方向 // 推导可发现,r方向的力形成P,SV波, theta方向的力形成SH波 diff --git a/pygrt/C_extension/src/lamb/elliptic.c b/pygrt/C_extension/src/lamb/elliptic.c index 518e40c6..5641afbb 100644 --- a/pygrt/C_extension/src/lamb/elliptic.c +++ b/pygrt/C_extension/src/lamb/elliptic.c @@ -20,7 +20,7 @@ #include "grt/common/checkerror.h" -MYREAL grt_ellipticK(const MYREAL m) +real_t grt_ellipticK(const real_t m) { if(m <= 0.0 || m>= 1.0){ GRTRaiseError("For the first complete elliptic integral, m should be in (0,1), but get %f.\n", m); @@ -28,9 +28,9 @@ MYREAL grt_ellipticK(const MYREAL m) #define N 5 - MYREAL m0 = 1 - m; + real_t m0 = 1 - m; - static const MYREAL a[N] = { + static const real_t a[N] = { 1.38629436112, 0.09666344259, 0.03590092383, @@ -38,7 +38,7 @@ MYREAL grt_ellipticK(const MYREAL m) 0.01451196212, }; - static const MYREAL b[N] = { + static const real_t b[N] = { 0.50000000000, 0.12498593597, 0.06880248576, @@ -46,8 +46,8 @@ MYREAL grt_ellipticK(const MYREAL m) 0.00441787012, }; - MYREAL K1 = 0.0, K2 = 0.0; - MYREAL p = 1.0; + real_t K1 = 0.0, K2 = 0.0; + real_t p = 1.0; K1 += a[0]; K2 += b[0]; for(int i = 1; i < N; ++i){ @@ -60,7 +60,7 @@ MYREAL grt_ellipticK(const MYREAL m) #undef N } -MYREAL grt_ellipticE(const MYREAL m) +real_t grt_ellipticE(const real_t m) { if(m <= 0.0 || m>= 1.0){ GRTRaiseError("For the second complete elliptic integral, m should be in (0,1), but get %f.\n", m); @@ -68,9 +68,9 @@ MYREAL grt_ellipticE(const MYREAL m) #define N 5 - MYREAL m0 = 1 - m; + real_t m0 = 1 - m; - static const MYREAL a[N] = { + static const real_t a[N] = { 1.00000000000, 0.44325141463, 0.06260601220, @@ -78,7 +78,7 @@ MYREAL grt_ellipticE(const MYREAL m) 0.01736506451, }; - static const MYREAL b[N] = { + static const real_t b[N] = { 0.00000000000, 0.24998368310, 0.09200180037, @@ -86,8 +86,8 @@ MYREAL grt_ellipticE(const MYREAL m) 0.00526449639 }; - MYREAL K1 = 0.0, K2 = 0.0; - MYREAL p = 1.0; + real_t K1 = 0.0, K2 = 0.0; + real_t p = 1.0; K1 += a[0]; K2 += b[0]; for(int i = 1; i < N; ++i){ @@ -100,29 +100,29 @@ MYREAL grt_ellipticE(const MYREAL m) #undef N } -static MYREAL MAX3(const MYREAL x, const MYREAL y, const MYREAL z) +static real_t MAX3(const real_t x, const real_t y, const real_t z) { - MYREAL xy = ((x) > (y) ? (x) : (y)); + real_t xy = ((x) > (y) ? (x) : (y)); return ((xy) > (z) ? (xy) : (z)); } -static MYREAL MAX4(const MYREAL x, const MYREAL y, const MYREAL z, const MYREAL p) +static real_t MAX4(const real_t x, const real_t y, const real_t z, const real_t p) { - MYREAL xy = ((x) > (y) ? (x) : (y)); - MYREAL xyz = ((xy) > (z) ? (xy) : (z)); + real_t xy = ((x) > (y) ? (x) : (y)); + real_t xyz = ((xy) > (z) ? (xy) : (z)); return ((xyz) > (p) ? (xyz) : (p)); } -static MYCOMPLEX grt_ellipticRF(const MYCOMPLEX x, const MYCOMPLEX y, const MYCOMPLEX z) +static cplx_t grt_ellipticRF(const cplx_t x, const cplx_t y, const cplx_t z) { - MYREAL errtol = 0.001; + real_t errtol = 0.001; int nmax = 10000; - MYREAL eps; - MYCOMPLEX lam; - MYCOMPLEX X, Y, Z, An; - MYCOMPLEX xn, yn, zn; - MYCOMPLEX xr, yr, zr; + real_t eps; + cplx_t lam; + cplx_t X, Y, Z, An; + cplx_t xn, yn, zn; + cplx_t xr, yr, zr; int n = 0; xn = x; yn = y; @@ -146,11 +146,11 @@ static MYCOMPLEX grt_ellipticRF(const MYCOMPLEX x, const MYCOMPLEX y, const MYCO if(n >= nmax) break; } - MYCOMPLEX E2, E3; + cplx_t E2, E3; E2 = X*Y - Z*Z; E3 = X*Y*Z; - MYCOMPLEX R; + cplx_t R; R = 1.0 - 0.1*E2 + 1/14.0 * E3 + 1/24.0 * E2*E2 - 3/44.0 *E2*E3; R /= sqrt(An); @@ -158,18 +158,18 @@ static MYCOMPLEX grt_ellipticRF(const MYCOMPLEX x, const MYCOMPLEX y, const MYCO } -static MYCOMPLEX grt_ellipticRJ(const MYCOMPLEX x, const MYCOMPLEX y, const MYCOMPLEX z, const MYCOMPLEX p) +static cplx_t grt_ellipticRJ(const cplx_t x, const cplx_t y, const cplx_t z, const cplx_t p) { - MYREAL errtol = 0.001; + real_t errtol = 0.001; int nmax = 10000; - MYREAL eps; - MYCOMPLEX lam, dm, em; - MYCOMPLEX X, Y, Z, P, An; - MYCOMPLEX xn, yn, zn, pn; - MYCOMPLEX xr, yr, zr, pr; - MYCOMPLEX sum1 = 0.0; - MYREAL pow4 = 1.0; + real_t eps; + cplx_t lam, dm, em; + cplx_t X, Y, Z, P, An; + cplx_t xn, yn, zn, pn; + cplx_t xr, yr, zr, pr; + cplx_t sum1 = 0.0; + real_t pow4 = 1.0; int n = 0; xn = x; yn = y; @@ -203,14 +203,14 @@ static MYCOMPLEX grt_ellipticRJ(const MYCOMPLEX x, const MYCOMPLEX y, const MYCO if(n >= nmax) break; } - MYCOMPLEX E2, E3, E4, E5; - MYCOMPLEX PPP = P*P*P; + cplx_t E2, E3, E4, E5; + cplx_t PPP = P*P*P; E2 = X*Y + Y*Z + X*Z - 3.0*P*P; E3 = X*Y*Z + 3.0*E2*P + 4.0*PPP; E4 = (2.0*X*Y*Z + E2*P + 3.0*PPP)*P; E5 = X*Y*Z*P*P; - MYCOMPLEX res = 0.0; + cplx_t res = 0.0; res = 1.0 - 3.0/14 * E2 + 1.0/6 * E3 + 9.0/88*E2*E2 - 3.0/22*E4 - 9.0/52*E2*E3 + 3.0/26*E5; res /= An * sqrt(An); res *= pow4; @@ -221,7 +221,7 @@ static MYCOMPLEX grt_ellipticRJ(const MYCOMPLEX x, const MYCOMPLEX y, const MYCO } -MYCOMPLEX grt_ellipticPi(const MYCOMPLEX n, const MYREAL m) +cplx_t grt_ellipticPi(const cplx_t n, const real_t m) { if(m <= 0.0 || m>= 1.0){ GRTRaiseError("For the third complete elliptic integral, m should be in (0,1), but get %f.\n", m); diff --git a/pygrt/C_extension/src/lamb/grt_lamb1.c b/pygrt/C_extension/src/lamb/grt_lamb1.c index a49c8465..c05b1f25 100644 --- a/pygrt/C_extension/src/lamb/grt_lamb1.c +++ b/pygrt/C_extension/src/lamb/grt_lamb1.c @@ -21,20 +21,20 @@ typedef struct { /** 模型参数 */ struct { bool active; - MYREAL nu; ///< 泊松比 + real_t nu; ///< 泊松比 } P; /** 归一化时间序列 */ struct { bool active; - MYREAL *ts; + real_t *ts; int nt; } T; /** 方位角 */ struct { bool active; - MYREAL azimuth; ///< 方位角,单位为度 + real_t azimuth; ///< 方位角,单位为度 } A; } GRT_MODULE_CTRL; @@ -107,7 +107,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ case 'T': Ctrl->T.active = true; { - MYREAL a1, a2, delta; + real_t a1, a2, delta; if(3 != sscanf(optarg, "%lf/%lf/%lf", &a1, &a2, &delta)){ GRTBadOptionError(command, T, ""); }; @@ -122,7 +122,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ } Ctrl->T.nt = floor((a2-a1)/delta) + 1; - Ctrl->T.ts = (MYREAL*)calloc(Ctrl->T.nt, sizeof(MYREAL)); + Ctrl->T.ts = (real_t*)calloc(Ctrl->T.nt, sizeof(real_t)); for(int i=0; iT.nt; ++i){ Ctrl->T.ts[i] = a1 + delta*i; } diff --git a/pygrt/C_extension/src/lamb/lamb1.c b/pygrt/C_extension/src/lamb/lamb1.c index 2dff2f07..2991640c 100644 --- a/pygrt/C_extension/src/lamb/lamb1.c +++ b/pygrt/C_extension/src/lamb/lamb1.c @@ -19,49 +19,49 @@ typedef struct { // 常量,不随时间变化 - MYREAL nu; - MYREAL k; ///< k = beta/alpha; - MYREAL kk; ///< kk = k*k; - MYREAL kp; ///< kp = sqrt(1 - k^2); - MYREAL kpkp; ///< kpkp = kp*kp; - MYREAL h1; ///< h1 = 2 * kk - 1; - MYREAL h2; ///< h2 = kk * kpkp; - MYREAL h3; ///< h3 = kpkp * h1^2; - MYREAL b3; ///< b3 = 3 * kp^2 - 1; - MYREAL b6; ///< b6 = 6 * kp^2 - 1; - MYREAL c; ///< c = 6 * h2 - 1; - MYREAL sf; ///< sin(phi) - MYREAL cf; ///< cos(phi) + real_t nu; + real_t k; ///< k = beta/alpha; + real_t kk; ///< kk = k*k; + real_t kp; ///< kp = sqrt(1 - k^2); + real_t kpkp; ///< kpkp = kp*kp; + real_t h1; ///< h1 = 2 * kk - 1; + real_t h2; ///< h2 = kk * kpkp; + real_t h3; ///< h3 = kpkp * h1^2; + real_t b3; ///< b3 = 3 * kp^2 - 1; + real_t b6; ///< b6 = 6 * kp^2 - 1; + real_t c; ///< c = 6 * h2 - 1; + real_t sf; ///< sin(phi) + real_t cf; ///< cos(phi) // 一元三次方程的三个根,其中 y30 为正根 - MYCOMPLEX ys[3]; - MYCOMPLEX ysp[3]; + cplx_t ys[3]; + cplx_t ysp[3]; // Rayleigh 波项 - MYREAL kpa; - MYREAL kpakpa; - MYREAL RaylR; - MYREAL RaylQ[3][3]; + real_t kpa; + real_t kpakpa; + real_t RaylR; + real_t RaylQ[3][3]; } VARS; -static void ckim_P(MYREAL tbar, VARS *V, MYCOMPLEX ckim[3][6][10]) +static void ckim_P(real_t tbar, VARS *V, cplx_t ckim[3][6][10]) { - MYREAL kk = V->kk; - MYREAL h1 = V->h1; - MYREAL h2 = V->h2; - MYREAL h3 = V->h3; - MYREAL b3 = V->b3; - MYREAL b6 = V->b6; - MYREAL c = V->c; + real_t kk = V->kk; + real_t h1 = V->h1; + real_t h2 = V->h2; + real_t h3 = V->h3; + real_t b3 = V->b3; + real_t b6 = V->b6; + real_t c = V->c; - MYREAL tbtb = tbar*tbar; - MYREAL A = 5.0 + 8.0*tbtb - 12.0*kk*tbtb - 24.0*h2; - MYREAL B = h1 * (1.0 + 5.0*tbtb - 6.0*kk*tbtb - 8.0*h2); - MYREAL d = 2.0*h1 - tbtb; + real_t tbtb = tbar*tbar; + real_t A = 5.0 + 8.0*tbtb - 12.0*kk*tbtb - 24.0*h2; + real_t B = h1 * (1.0 + 5.0*tbtb - 6.0*kk*tbtb - 8.0*h2); + real_t d = 2.0*h1 - tbtb; - MYREAL aa = tbtb - kk; + real_t aa = tbtb - kk; // 表 6.6.1 // (ξ, i, m), ξ=1,2, i=1...5, m=0...9 @@ -117,19 +117,19 @@ static void ckim_P(MYREAL tbar, VARS *V, MYCOMPLEX ckim[3][6][10]) ckim[2][5][3] = -4.0*I*h1*h2; } -static void ckim_S1(MYREAL tbar, VARS *V, MYCOMPLEX ckim[3][6][10]) +static void ckim_S1(real_t tbar, VARS *V, cplx_t ckim[3][6][10]) { - MYREAL kk = V->kk; - MYREAL kpkp = V->kpkp; + real_t kk = V->kk; + real_t kpkp = V->kpkp; - MYREAL tbtb = tbar*tbar; - MYREAL bb = tbtb - 1.0; - MYREAL A = 15.0 - 10.0*tbtb + 8.0*kk*tbtb - 16.0*kk; - MYREAL B = 9.0 - 10.0*tbtb + 16.0*kk*bb; - MYREAL C0 = 4.0 - 3.0*tbtb + 2.0*kk*tbtb - 3.0*kk; + real_t tbtb = tbar*tbar; + real_t bb = tbtb - 1.0; + real_t A = 15.0 - 10.0*tbtb + 8.0*kk*tbtb - 16.0*kk; + real_t B = 9.0 - 10.0*tbtb + 16.0*kk*bb; + real_t C0 = 4.0 - 3.0*tbtb + 2.0*kk*tbtb - 3.0*kk; - MYREAL h[6][10] = {0}; + real_t h[6][10] = {0}; for(int i=0; i<6; ++i){ for(int j=0; j<10; ++j){ h[i][j] = i + j*kpkp; @@ -201,7 +201,7 @@ static void ckim_S1(MYREAL tbar, VARS *V, MYCOMPLEX ckim[3][6][10]) ckim[2][5][1] = -2.0*kpkp; } -static void ckim_S2(MYREAL tbar, VARS *V, MYCOMPLEX ckim[3][6][10]) +static void ckim_S2(real_t tbar, VARS *V, cplx_t ckim[3][6][10]) { ckim_S1(tbar, V, ckim); @@ -220,10 +220,10 @@ static void ckim_S2(MYREAL tbar, VARS *V, MYCOMPLEX ckim[3][6][10]) } } -static void Cmat(VARS *V, MYCOMPLEX ckim[3][6][10], MYCOMPLEX C[3][3][3][10]) +static void Cmat(VARS *V, cplx_t ckim[3][6][10], cplx_t C[3][3][3][10]) { - MYREAL sf = V->sf; - MYREAL cf = V->cf; + real_t sf = V->sf; + real_t cf = V->cf; // 构建矩阵 C (3x3) for(int k=1; k<=2; ++k){ for(int m=0; m<=9; ++m){ @@ -242,16 +242,16 @@ static void Cmat(VARS *V, MYCOMPLEX ckim[3][6][10], MYCOMPLEX C[3][3][3][10]) } } -static void uv(int sgn, MYREAL kpkp, MYCOMPLEX ys[3], MYCOMPLEX C[3][3][3][10], MYCOMPLEX u[3][3][10], MYCOMPLEX v[3][3][11]) +static void uv(int sgn, real_t kpkp, cplx_t ys[3], cplx_t C[3][3][3][10], cplx_t u[3][3][10], cplx_t v[3][3][11]) { - MYCOMPLEX y1, y2, y3; + cplx_t y1, y2, y3; y1 = ys[0]*sgn; y2 = ys[1]*sgn; y3 = ys[2]*sgn; // 构建系数 u_ij,m 和 v_ij,n // i,j = 1,2,3 m=1...9 n=1...10 - MYCOMPLEX delta[4] = {0}; + cplx_t delta[4] = {0}; delta[1] = -16.0 * sgn * kpkp *(y1 - y2)*(y1 - y3); delta[2] = -16.0 * sgn * kpkp *(y2 - y1)*(y2 - y3); delta[3] = -16.0 * sgn * kpkp *(y3 - y1)*(y3 - y2); @@ -283,14 +283,14 @@ static void uv(int sgn, MYREAL kpkp, MYCOMPLEX ys[3], MYCOMPLEX C[3][3][3][10], } } -static MYCOMPLEX U_P(int n, MYREAL tbar, MYCOMPLEX c, VARS *V) +static cplx_t U_P(int n, real_t tbar, cplx_t c, VARS *V) { - MYREAL kk = V->kk; - MYREAL tbtb = tbar*tbar; - MYREAL aa = tbtb - kk; - MYREAL a = sqrt(aa); // 要求 tbar >= k + real_t kk = V->kk; + real_t tbtb = tbar*tbar; + real_t aa = tbtb - kk; + real_t a = sqrt(aa); // 要求 tbar >= k - MYCOMPLEX xi = sqrt(aa - c + 0*I); + cplx_t xi = sqrt(aa - c + 0*I); if(n==1){ if(cimag(c) == 0.0 && creal(c) < aa){ @@ -318,23 +318,23 @@ static MYCOMPLEX U_P(int n, MYREAL tbar, MYCOMPLEX c, VARS *V) GRTRaiseError("Wrong execution in function %s.", __func__); } -static MYCOMPLEX V_P(int n, MYREAL tbar, MYCOMPLEX c, VARS *V) +static cplx_t V_P(int n, real_t tbar, cplx_t c, VARS *V) { - MYREAL kp = V->kp; - MYREAL kk = V->kk; - MYREAL kpkp = V->kpkp; - MYREAL tbtb = tbar*tbar; - MYREAL aa = tbtb - kk; - MYREAL a = sqrt(aa); // 要求 tbar >= k - MYREAL bb = tbtb - 1.0; - - MYCOMPLEX xi = sqrt(aa - c + 0*I); - MYCOMPLEX et = sqrt(kpkp - c + 0*I); + real_t kp = V->kp; + real_t kk = V->kk; + real_t kpkp = V->kpkp; + real_t tbtb = tbar*tbar; + real_t aa = tbtb - kk; + real_t a = sqrt(aa); // 要求 tbar >= k + real_t bb = tbtb - 1.0; + + cplx_t xi = sqrt(aa - c + 0*I); + cplx_t et = sqrt(kpkp - c + 0*I); - MYREAL m1 = aa/kpkp; - MYREAL m1p = kpkp/aa; - MYREAL m2 = bb/aa; - MYCOMPLEX n1 = bb/(xi*xi); + real_t m1 = aa/kpkp; + real_t m1p = kpkp/aa; + real_t m2 = bb/aa; + cplx_t n1 = bb/(xi*xi); if(tbar > 1.0){ if(n == 1){ @@ -388,13 +388,13 @@ static MYCOMPLEX V_P(int n, MYREAL tbar, MYCOMPLEX c, VARS *V) GRTRaiseError("Wrong execution in function %s.", __func__); } -static MYCOMPLEX U_S1(int n, MYREAL tbar, MYCOMPLEX c, VARS *V) +static cplx_t U_S1(int n, real_t tbar, cplx_t c, VARS *V) { - MYREAL tbtb = tbar*tbar; - MYREAL bb = tbtb - 1.0; - MYREAL b = sqrt(bb); // 要求 tbar >= 1.0 + real_t tbtb = tbar*tbar; + real_t bb = tbtb - 1.0; + real_t b = sqrt(bb); // 要求 tbar >= 1.0 - MYCOMPLEX xip = sqrt(bb - c + 0*I); + cplx_t xip = sqrt(bb - c + 0*I); if(n==1){ if(cimag(c) == 0.0 && creal(c) < bb){ @@ -422,22 +422,22 @@ static MYCOMPLEX U_S1(int n, MYREAL tbar, MYCOMPLEX c, VARS *V) GRTRaiseError("Wrong execution in function %s.", __func__); } -static MYCOMPLEX V_S1(int n, MYREAL tbar, MYCOMPLEX c, VARS *V) +static cplx_t V_S1(int n, real_t tbar, cplx_t c, VARS *V) { - MYREAL kp = V->kp; - MYREAL kk = V->kk; - MYREAL kpkp = V->kpkp; - MYREAL tbtb = tbar*tbar; - MYREAL aa = tbtb - kk; - MYREAL a = sqrt(aa); // 要求 tbar >= k - MYREAL bb = tbtb - 1.0; - MYREAL b = sqrt(bb); // 要求 tbar >= 1.0 - - MYCOMPLEX xip = sqrt(bb - c + 0*I); - MYCOMPLEX etp = sqrt(kpkp + c + 0*I); + real_t kp = V->kp; + real_t kk = V->kk; + real_t kpkp = V->kpkp; + real_t tbtb = tbar*tbar; + real_t aa = tbtb - kk; + real_t a = sqrt(aa); // 要求 tbar >= k + real_t bb = tbtb - 1.0; + real_t b = sqrt(bb); // 要求 tbar >= 1.0 + + cplx_t xip = sqrt(bb - c + 0*I); + cplx_t etp = sqrt(kpkp + c + 0*I); - MYREAL m2 = bb/aa; - MYCOMPLEX n2 = bb/(xip*xip); + real_t m2 = bb/aa; + cplx_t n2 = bb/(xip*xip); if(n==1){ if(cimag(c) == 0.0 && creal(c) > - kpkp && creal(c) < bb){ @@ -464,19 +464,19 @@ static MYCOMPLEX V_S1(int n, MYREAL tbar, MYCOMPLEX c, VARS *V) GRTRaiseError("Wrong execution in function %s.", __func__); } -static MYCOMPLEX U_S2(int n, MYREAL tbar, MYCOMPLEX c, VARS *V) +static cplx_t U_S2(int n, real_t tbar, cplx_t c, VARS *V) { - MYREAL kp = V->kp; - MYREAL kk = V->kk; - MYREAL tbtb = tbar*tbar; - MYREAL aa = tbtb - kk; - MYREAL a = sqrt(aa); // 要求 tbar >= k - MYREAL bb = tbtb - 1.0; - MYREAL b = sqrt(bb); // 要求 tbar >= 1.0 - MYREAL f = kp + a; - MYCOMPLEX d = sqrt(c*(bb - c)); - - MYCOMPLEX xip = sqrt(bb - c + 0*I); + real_t kp = V->kp; + real_t kk = V->kk; + real_t tbtb = tbar*tbar; + real_t aa = tbtb - kk; + real_t a = sqrt(aa); // 要求 tbar >= k + real_t bb = tbtb - 1.0; + real_t b = sqrt(bb); // 要求 tbar >= 1.0 + real_t f = kp + a; + cplx_t d = sqrt(c*(bb - c)); + + cplx_t xip = sqrt(bb - c + 0*I); if(n == 1){ return 1.0/(I*xip) * atan((a - b)*I*xip/(c + b*(a-b))); @@ -491,29 +491,29 @@ static MYCOMPLEX U_S2(int n, MYREAL tbar, MYCOMPLEX c, VARS *V) return 0.5*(f-b)*(f-b)/f; } else if(n == 5){ - MYREAL f4 = f*f*f*f; - MYREAL b4 = bb*bb; + real_t f4 = f*f*f*f; + real_t b4 = bb*bb; return 0.125*(f4 - b4)/(f*f) - 0.5*bb*log(f/b); } GRTRaiseError("Wrong execution in function %s.", __func__); } -static MYCOMPLEX V_S2(int n, MYREAL tbar, MYCOMPLEX c, VARS *V) +static cplx_t V_S2(int n, real_t tbar, cplx_t c, VARS *V) { - MYREAL kp = V->kp; - MYREAL kk = V->kk; - MYREAL kpkp = V->kpkp; - MYREAL tbtb = tbar*tbar; - MYREAL aa = tbtb - kk; - MYREAL a = sqrt(aa); // 要求 tbar >= k - MYREAL bb = tbtb - 1.0; - MYREAL b = sqrt(bb); // 要求 tbar >= 1.0 - - MYCOMPLEX xip = sqrt(bb - c + 0*I); - MYCOMPLEX etp = sqrt(kpkp + c + 0*I); + real_t kp = V->kp; + real_t kk = V->kk; + real_t kpkp = V->kpkp; + real_t tbtb = tbar*tbar; + real_t aa = tbtb - kk; + real_t a = sqrt(aa); // 要求 tbar >= k + real_t bb = tbtb - 1.0; + real_t b = sqrt(bb); // 要求 tbar >= 1.0 + + cplx_t xip = sqrt(bb - c + 0*I); + cplx_t etp = sqrt(kpkp + c + 0*I); - MYREAL m1p = kpkp/aa; - MYCOMPLEX n3 = kpkp/(etp*etp); + real_t m1p = kpkp/aa; + cplx_t n3 = kpkp/(etp*etp); if(n==1){ if(cimag(c) == 0.0 && creal(c) > - kpkp && creal(c) < bb){ @@ -540,19 +540,19 @@ static MYCOMPLEX V_S2(int n, MYREAL tbar, MYCOMPLEX c, VARS *V) GRTRaiseError("Wrong execution in function %s.", __func__); } -static MYCOMPLEX U_SP(int n, MYREAL tbar, MYCOMPLEX c, VARS *V) +static cplx_t U_SP(int n, real_t tbar, cplx_t c, VARS *V) { - MYREAL kk = V->kk; - MYREAL kp = V->kp; - MYREAL tbtb = tbar*tbar; - MYREAL aa = tbtb - kk; - MYREAL a = sqrt(aa); // 要求 tbar >= k - MYREAL bb = tbtb - 1.0; - MYREAL bp = sqrt(- bb); // 要求 tbar <= 1.0 - MYREAL f = kp + a; - MYCOMPLEX d = sqrt(c*(bb - c)); - - MYCOMPLEX xip = sqrt(bb - c + 0*I); + real_t kk = V->kk; + real_t kp = V->kp; + real_t tbtb = tbar*tbar; + real_t aa = tbtb - kk; + real_t a = sqrt(aa); // 要求 tbar >= k + real_t bb = tbtb - 1.0; + real_t bp = sqrt(- bb); // 要求 tbar <= 1.0 + real_t f = kp + a; + cplx_t d = sqrt(c*(bb - c)); + + cplx_t xip = sqrt(bb - c + 0*I); if(n == 1){ return - 1.0/(I*xip) * atan(a/(I*xip)); @@ -567,26 +567,26 @@ static MYCOMPLEX U_SP(int n, MYREAL tbar, MYCOMPLEX c, VARS *V) return 0.5*(f*f + bb)/f; } else if(n == 5){ - MYREAL f4 = f*f*f*f; - MYREAL b4 = bb*bb; + real_t f4 = f*f*f*f; + real_t b4 = bb*bb; return 0.125*(f4 - b4)/(f*f) - 0.5*bb*log(f/bp); } GRTRaiseError("Wrong execution in function %s.", __func__); } -static MYCOMPLEX V_SP(int n, MYREAL tbar, MYCOMPLEX c, VARS *V) +static cplx_t V_SP(int n, real_t tbar, cplx_t c, VARS *V) { - MYREAL kk = V->kk; - MYREAL kp = V->kp; - MYREAL kpkp = V->kpkp; - MYREAL tbtb = tbar*tbar; - MYREAL aa = tbtb - kk; - MYREAL bb = tbtb - 1.0; - - MYCOMPLEX etp = sqrt(kpkp + c + 0*I); + real_t kk = V->kk; + real_t kp = V->kp; + real_t kpkp = V->kpkp; + real_t tbtb = tbar*tbar; + real_t aa = tbtb - kk; + real_t bb = tbtb - 1.0; + + cplx_t etp = sqrt(kpkp + c + 0*I); - MYREAL m1 = aa/kpkp; - MYCOMPLEX n4 = aa/(etp*etp); + real_t m1 = aa/kpkp; + cplx_t n4 = aa/(etp*etp); if(n==1){ if(cimag(c) == 0.0 && creal(c) > - kpkp && creal(c) < bb){ @@ -614,24 +614,24 @@ static MYCOMPLEX V_SP(int n, MYREAL tbar, MYCOMPLEX c, VARS *V) GRTRaiseError("Wrong execution in function %s.", __func__); } -static void build_raw(MYREAL tbar, VARS *V, MYCOMPLEX u[3][3], - MYCOMPLEX C[3][3][3][10], - MYCOMPLEX (*Ufunc)(int, MYREAL, MYCOMPLEX, VARS *), - MYCOMPLEX (*Vfunc)(int, MYREAL, MYCOMPLEX, VARS *), - int sgn, MYCOMPLEX ys[3]) +static void build_raw(real_t tbar, VARS *V, cplx_t u[3][3], + cplx_t C[3][3][3][10], + cplx_t (*Ufunc)(int, real_t, cplx_t, VARS *), + cplx_t (*Vfunc)(int, real_t, cplx_t, VARS *), + int sgn, cplx_t ys[3]) { // 计算 u, v 系数 - MYCOMPLEX uc[3][3][10] = {0}; - MYCOMPLEX vc[3][3][11] = {0}; + cplx_t uc[3][3][10] = {0}; + cplx_t vc[3][3][11] = {0}; uv(sgn, V->kpkp, ys, C, uc, vc); - MYCOMPLEX y1, y2, y3; + cplx_t y1, y2, y3; y1 = ys[0]; y2 = ys[1]; y3 = ys[2]; // 计算 UP, VP 函数 - MYCOMPLEX UF[6][4] = {0}; + cplx_t UF[6][4] = {0}; UF[1][1] = Ufunc(1, tbar, y1, V); UF[1][2] = Ufunc(1, tbar, y2, V); UF[1][3] = Ufunc(1, tbar, y3, V); @@ -642,7 +642,7 @@ static void build_raw(MYREAL tbar, VARS *V, MYCOMPLEX u[3][3], UF[i][1] = Ufunc(i, tbar, y3, V); } - MYCOMPLEX VF[7][4] = {0}; + cplx_t VF[7][4] = {0}; VF[1][1] = Vfunc(1, tbar, y1, V); VF[1][2] = Vfunc(1, tbar, y2, V); VF[1][3] = Vfunc(1, tbar, y3, V); @@ -656,7 +656,7 @@ static void build_raw(MYREAL tbar, VARS *V, MYCOMPLEX u[3][3], // 组合成 P 波项 for(int i1=0; i1<3; ++i1){ for(int i2=0; i2<3; ++i2){ - MYCOMPLEX tmp = 0.0; + cplx_t tmp = 0.0; tmp += uc[i1][i2][1]*UF[1][1] + uc[i1][i2][2]*UF[2][1] + uc[i1][i2][3]*UF[1][2] + uc[i1][i2][4]*UF[2][2] + uc[i1][i2][5]*UF[1][3] + uc[i1][i2][6]*UF[2][3]; @@ -674,14 +674,14 @@ static void build_raw(MYREAL tbar, VARS *V, MYCOMPLEX u[3][3], } } -static void build_P(MYREAL tbar, VARS *V, MYREAL u[3][3]) +static void build_P(real_t tbar, VARS *V, real_t u[3][3]) { - MYREAL k = V->k; + real_t k = V->k; if(tbar < k) return; - MYCOMPLEX cu[3][3]; - MYCOMPLEX ckim[3][6][10] = {0}; - MYCOMPLEX C[3][3][3][10] = {0}; + cplx_t cu[3][3]; + cplx_t ckim[3][6][10] = {0}; + cplx_t C[3][3][3][10] = {0}; ckim_P(tbar, V, ckim); Cmat(V, ckim, C); build_raw(tbar, V, cu, C, U_P, V_P, 1, V->ys); @@ -692,13 +692,13 @@ static void build_P(MYREAL tbar, VARS *V, MYREAL u[3][3]) } } -static void build_S1(MYREAL tbar, VARS *V, MYREAL u[3][3]) +static void build_S1(real_t tbar, VARS *V, real_t u[3][3]) { if(tbar < 1.0) return; - MYCOMPLEX cu[3][3]; - MYCOMPLEX ckim[3][6][10] = {0}; - MYCOMPLEX C[3][3][3][10] = {0}; + cplx_t cu[3][3]; + cplx_t ckim[3][6][10] = {0}; + cplx_t C[3][3][3][10] = {0}; ckim_S1(tbar, V, ckim); Cmat(V, ckim, C); build_raw(tbar, V, cu, C, U_S1, V_S1, 1, V->ysp); @@ -709,14 +709,14 @@ static void build_S1(MYREAL tbar, VARS *V, MYREAL u[3][3]) } } -static void build_S2_SP(MYREAL tbar, VARS *V, MYREAL us2[3][3], MYREAL usp[3][3]) +static void build_S2_SP(real_t tbar, VARS *V, real_t us2[3][3], real_t usp[3][3]) { - MYREAL k = V->k; + real_t k = V->k; if(tbar < k) return; - MYCOMPLEX cu[3][3]; - MYCOMPLEX ckim[3][6][10] = {0}; - MYCOMPLEX C[3][3][3][10] = {0}; + cplx_t cu[3][3]; + cplx_t ckim[3][6][10] = {0}; + cplx_t C[3][3][3][10] = {0}; ckim_S2(tbar, V, ckim); Cmat(V, ckim, C); @@ -739,13 +739,13 @@ static void build_S2_SP(MYREAL tbar, VARS *V, MYREAL us2[3][3], MYREAL usp[3][3] } -static void build_R(MYREAL tbar, VARS *V, MYREAL u[3][3]) +static void build_R(real_t tbar, VARS *V, real_t u[3][3]) { - MYREAL kpakpa = V->kpakpa; - MYREAL kpa = V->kpa; + real_t kpakpa = V->kpakpa; + real_t kpa = V->kpa; if(tbar < kpa) return; - MYREAL coef = M_PI_4 * tbar * V->RaylR / sqrt(tbar*tbar - kpakpa); + real_t coef = M_PI_4 * tbar * V->RaylR / sqrt(tbar*tbar - kpakpa); for(int i1=0; i1<3; ++i1){ for(int i2=0; i2<3; ++i2){ @@ -759,7 +759,7 @@ static void build_R(MYREAL tbar, VARS *V, MYREAL u[3][3]) void grt_solve_lamb1( - const MYREAL nu, const MYREAL *ts, const int nt, const MYREAL azimuth, MYREAL (*u)[3][3]) + const real_t nu, const real_t *ts, const int nt, const real_t azimuth, real_t (*u)[3][3]) { // 检查泊松比范围 if(nu <= 0.0 || nu >= 0.5){ @@ -783,9 +783,9 @@ void grt_solve_lamb1( printf("\n"); } - MYREAL phi = azimuth * DEG1; + real_t phi = azimuth * DEG1; - MYREAL tbar_eps = GRT_MIN(1e-8, (ts[1]-ts[0]) * 1e-5); + real_t tbar_eps = GRT_MIN(1e-8, (ts[1]-ts[0]) * 1e-5); // 初始化相关变量 VARS V0 = {0}; @@ -806,14 +806,14 @@ void grt_solve_lamb1( // 求一元三次方程的根 { - MYREAL a, b, c; - MYREAL nu2, nu3, nu4; + real_t a, b, c; + real_t nu2, nu3, nu4; nu2 = V0.nu*V0.nu; nu3 = nu2*V0.nu; nu4 = nu3*V0.nu; - MYREAL snu = 1.0 - V0.nu; - MYREAL snu2 = snu*snu; - MYREAL snu3 = snu2*snu; + real_t snu = 1.0 - V0.nu; + real_t snu2 = snu*snu; + real_t snu3 = snu2*snu; a = -0.5 * (2.0*nu2 + 1.0)/snu; b = 0.25 * (4.0*nu3 - 4.0*nu2 + 4.0*V0.nu - 1.0)/snu2; c = -0.125*nu4/snu3; @@ -825,9 +825,9 @@ void grt_solve_lamb1( // 另一种形式的Rayleigh波函数 { - MYCOMPLEX y3[3]; - MYREAL a, b, c; - MYREAL m = V0.kk; + cplx_t y3[3]; + real_t a, b, c; + real_t m = V0.kk; a = 0.5*(2.0*m - 3.0)/(1 - m); b = 0.5/(1 - m); c = - 0.0625/(1 - m); @@ -841,22 +841,22 @@ void grt_solve_lamb1( V0.kpakpa = creal(y3[2]); V0.kpa = sqrt(V0.kpakpa); - MYREAL u0 = sqrt(V0.kpakpa - V0.kk); - MYREAL v0 = sqrt(V0.kpakpa - 1.0); + real_t u0 = sqrt(V0.kpakpa - V0.kk); + real_t v0 = sqrt(V0.kpakpa - 1.0); - MYREAL R1, R2; + real_t R1, R2; R1 = (1.0 - 2.0*V0.kpakpa)*u0*v0 + 2.0*u0*u0*v0*v0; R2 = 2.0*(1.0 - 2.0*V0.kpakpa)*u0*v0 + 2.0*u0*u0*v0*v0 + V0.kpakpa*(u0*u0 + v0*v0); V0.RaylR = R1 / R2; } for(int i=0; i < nt; ++i){ - MYREAL up[3][3] = {0}; - MYREAL us1[3][3] = {0}; - MYREAL us2[3][3] = {0}; - MYREAL usp[3][3] = {0}; - MYREAL uR[3][3] = {0}; - MYREAL tbar = ts[i]; + real_t up[3][3] = {0}; + real_t us1[3][3] = {0}; + real_t us2[3][3] = {0}; + real_t usp[3][3] = {0}; + real_t uR[3][3] = {0}; + real_t tbar = ts[i]; // 跳过一些震相到时处的奇点 if(tbar == 1.0 || tbar == V0.k || tbar == V0.kpa) tbar += tbar_eps; diff --git a/pygrt/C_extension/src/lamb/lamb_util.c b/pygrt/C_extension/src/lamb/lamb_util.c index 0d32bd7d..3520610c 100644 --- a/pygrt/C_extension/src/lamb/lamb_util.c +++ b/pygrt/C_extension/src/lamb/lamb_util.c @@ -8,7 +8,7 @@ #include "grt/lamb/lamb_util.h" -void grt_roots3(double a, double b, double c, MYCOMPLEX y3[3]) +void grt_roots3(double a, double b, double c, cplx_t y3[3]) { double Q, R; Q = (a*a - 3.0*b) / 9.0; @@ -37,10 +37,10 @@ void grt_roots3(double a, double b, double c, MYCOMPLEX y3[3]) } } -MYCOMPLEX grt_evalpoly2(const MYCOMPLEX *C, const int n, const MYCOMPLEX y, const int offset) +cplx_t grt_evalpoly2(const cplx_t *C, const int n, const cplx_t y, const int offset) { - MYCOMPLEX res = 0.0; - MYCOMPLEX p = 1.0; + cplx_t res = 0.0; + cplx_t p = 1.0; for(int i=0; i<=n; ++i){ res += C[2*i+offset] * p; p *= y; diff --git a/pygrt/C_extension/src/static/grt_static_greenfn.c b/pygrt/C_extension/src/static/grt_static_greenfn.c index 88720fc8..e1f080d9 100644 --- a/pygrt/C_extension/src/static/grt_static_greenfn.c +++ b/pygrt/C_extension/src/static/grt_static_greenfn.c @@ -38,8 +38,8 @@ typedef struct { /** 震源和接收器深度 */ struct { bool active; - MYREAL depsrc; - MYREAL deprcv; + real_t depsrc; + real_t deprcv; char *s_depsrc; char *s_deprcv; } D; @@ -47,17 +47,17 @@ typedef struct { struct { bool active; MYINT method; - MYREAL Length; - MYREAL filonLength; - MYREAL safilonTol; - MYREAL filonCut; + real_t Length; + real_t filonLength; + real_t safilonTol; + real_t filonCut; } L; /** 波数积分上限 */ struct { bool active; - MYREAL keps; - MYREAL k0; - MYREAL vmin; + real_t keps; + real_t k0; + real_t vmin; bool v_active; } K; /** 波数积分过程的核函数文件 */ @@ -69,13 +69,13 @@ typedef struct { struct { bool active; MYINT nx; - MYREAL *xs; + real_t *xs; } X; /** Y 坐标 */ struct { bool active; MYINT ny; - MYREAL *ys; + real_t *ys; } Y; /** 输出 nc 文件名 */ struct { @@ -88,7 +88,7 @@ typedef struct { } e; MYINT nr; - MYREAL *rs; + real_t *rs; } GRT_MODULE_CTRL; @@ -369,7 +369,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ case 'X': Ctrl->X.active = true; { - MYREAL a1, a2, delta; + real_t a1, a2, delta; if(3 != sscanf(optarg, "%lf/%lf/%lf", &a1, &a2, &delta)){ GRTBadOptionError(command, X, ""); }; @@ -381,7 +381,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ } Ctrl->X.nx = floor((a2-a1)/delta) + 1; - Ctrl->X.xs = (MYREAL*)calloc(Ctrl->X.nx, sizeof(MYREAL)); + Ctrl->X.xs = (real_t*)calloc(Ctrl->X.nx, sizeof(real_t)); for(int i=0; iX.nx; ++i){ Ctrl->X.xs[i] = a1 + delta*i; } @@ -392,7 +392,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ case 'Y': Ctrl->Y.active = true; { - MYREAL a1, a2, delta; + real_t a1, a2, delta; if(3 != sscanf(optarg, "%lf/%lf/%lf", &a1, &a2, &delta)){ GRTBadOptionError(command, Y, ""); }; @@ -404,7 +404,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ } Ctrl->Y.ny = floor((a2-a1)/delta) + 1; - Ctrl->Y.ys = (MYREAL*)calloc(Ctrl->Y.ny, sizeof(MYREAL)); + Ctrl->Y.ys = (real_t*)calloc(Ctrl->Y.ny, sizeof(real_t)); for(int i=0; iY.ny; ++i){ Ctrl->Y.ys[i] = a1 + delta*i; } @@ -441,7 +441,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ // 设置震中距数组 Ctrl->nr = Ctrl->X.nx*Ctrl->Y.ny; - Ctrl->rs = (MYREAL*)calloc(Ctrl->nr, sizeof(MYREAL)); + Ctrl->rs = (real_t*)calloc(Ctrl->nr, sizeof(real_t)); for(int ix=0; ixX.nx; ++ix){ for(int iy=0; iyY.ny; ++iy){ Ctrl->rs[iy + ix*Ctrl->Y.ny] = GRT_MAX(sqrt(GRT_SQUARE(Ctrl->X.xs[ix]) + GRT_SQUARE(Ctrl->Y.ys[iy])), GRT_MIN_DISTANCE); // 避免0震中距 @@ -469,7 +469,7 @@ int static_greenfn_main(int argc, char **argv){ GRT_MODEL1D *mod1d = Ctrl->M.mod1d; // 最大最小速度 - MYREAL vmin, vmax; + real_t vmin, vmax; grt_get_mod1d_vmin_vmax(mod1d, &vmin, &vmax); // 参考最小速度 @@ -496,9 +496,9 @@ int static_greenfn_main(int argc, char **argv){ } // 建立格林函数的浮点数 - MYREAL (*grn)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM] = (MYREAL (*)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]) calloc(Ctrl->nr, sizeof(*grn)); - MYREAL (*grn_uiz)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM] = (Ctrl->e.active)? (MYREAL (*)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]) calloc(Ctrl->nr, sizeof(*grn_uiz)) : NULL; - MYREAL (*grn_uir)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM] = (Ctrl->e.active)? (MYREAL (*)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]) calloc(Ctrl->nr, sizeof(*grn_uir)) : NULL; + real_t (*grn)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM] = (real_t (*)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]) calloc(Ctrl->nr, sizeof(*grn)); + real_t (*grn_uiz)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM] = (Ctrl->e.active)? (real_t (*)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]) calloc(Ctrl->nr, sizeof(*grn_uiz)) : NULL; + real_t (*grn_uir)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM] = (Ctrl->e.active)? (real_t (*)[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]) calloc(Ctrl->nr, sizeof(*grn_uir)) : NULL; //============================================================================== @@ -510,12 +510,12 @@ int static_greenfn_main(int argc, char **argv){ ); //============================================================================== - MYREAL src_va = mod1d->Va[mod1d->isrc]; - MYREAL src_vb = mod1d->Vb[mod1d->isrc]; - MYREAL src_rho = mod1d->Rho[mod1d->isrc]; - MYREAL rcv_va = mod1d->Va[mod1d->ircv]; - MYREAL rcv_vb = mod1d->Vb[mod1d->ircv]; - MYREAL rcv_rho = mod1d->Rho[mod1d->ircv]; + real_t src_va = mod1d->Va[mod1d->isrc]; + real_t src_vb = mod1d->Vb[mod1d->isrc]; + real_t src_rho = mod1d->Rho[mod1d->isrc]; + real_t rcv_va = mod1d->Va[mod1d->ircv]; + real_t rcv_vb = mod1d->Vb[mod1d->ircv]; + real_t rcv_rho = mod1d->Rho[mod1d->ircv]; // ================================================================================== @@ -530,30 +530,30 @@ int static_greenfn_main(int argc, char **argv){ int uir_varids[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]; // 创建 NC 文件 - GRT_NC_CHECK(nc_create(Ctrl->O.s_outgrid, NC_CLOBBER, &ncid)); + NC_CHECK(nc_create(Ctrl->O.s_outgrid, NC_CLOBBER, &ncid)); // 写入全局属性 - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_put_att) (ncid, NC_GLOBAL, "src_va", GRT_NC_MYREAL, 1, &src_va)); - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_put_att) (ncid, NC_GLOBAL, "src_vb", GRT_NC_MYREAL, 1, &src_vb)); - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_put_att) (ncid, NC_GLOBAL, "src_rho", GRT_NC_MYREAL, 1, &src_rho)); - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_put_att) (ncid, NC_GLOBAL, "rcv_va", GRT_NC_MYREAL, 1, &rcv_va)); - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_put_att) (ncid, NC_GLOBAL, "rcv_vb", GRT_NC_MYREAL, 1, &rcv_vb)); - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_put_att) (ncid, NC_GLOBAL, "rcv_rho", GRT_NC_MYREAL, 1, &rcv_rho)); + NC_CHECK(NC_FUNC_REAL(nc_put_att) (ncid, NC_GLOBAL, "src_va", NC_REAL, 1, &src_va)); + NC_CHECK(NC_FUNC_REAL(nc_put_att) (ncid, NC_GLOBAL, "src_vb", NC_REAL, 1, &src_vb)); + NC_CHECK(NC_FUNC_REAL(nc_put_att) (ncid, NC_GLOBAL, "src_rho", NC_REAL, 1, &src_rho)); + NC_CHECK(NC_FUNC_REAL(nc_put_att) (ncid, NC_GLOBAL, "rcv_va", NC_REAL, 1, &rcv_va)); + NC_CHECK(NC_FUNC_REAL(nc_put_att) (ncid, NC_GLOBAL, "rcv_vb", NC_REAL, 1, &rcv_vb)); + NC_CHECK(NC_FUNC_REAL(nc_put_att) (ncid, NC_GLOBAL, "rcv_rho", NC_REAL, 1, &rcv_rho)); // 是否计算了位移偏导也直接写到全局属性 { MYINT tmp = Ctrl->e.active; - GRT_NC_CHECK(GRT_NC_FUNC_MYINT(nc_put_att) (ncid, NC_GLOBAL, "calc_upar", GRT_NC_MYINT, 1, &tmp)); + NC_CHECK(GRT_NC_FUNC_MYINT(nc_put_att) (ncid, NC_GLOBAL, "calc_upar", GRT_NC_MYINT, 1, &tmp)); } // 定义维度 - GRT_NC_CHECK(nc_def_dim(ncid, "north", Ctrl->X.nx, &x_dimid)); - GRT_NC_CHECK(nc_def_dim(ncid, "east", Ctrl->Y.ny, &y_dimid)); + NC_CHECK(nc_def_dim(ncid, "north", Ctrl->X.nx, &x_dimid)); + NC_CHECK(nc_def_dim(ncid, "east", Ctrl->Y.ny, &y_dimid)); dimids[0] = x_dimid; dimids[1] = y_dimid; // 定义维度数组 - GRT_NC_CHECK(nc_def_var(ncid, "north", GRT_NC_MYREAL, 1, &x_dimid, &x_varid)); - GRT_NC_CHECK(nc_def_var(ncid, "east", GRT_NC_MYREAL, 1, &y_dimid, &y_varid)); + NC_CHECK(nc_def_var(ncid, "north", NC_REAL, 1, &x_dimid, &x_varid)); + NC_CHECK(nc_def_var(ncid, "east", NC_REAL, 1, &y_dimid, &y_varid)); // 定义不同震源不同分量的格林函数数组 for(int i=0; ie.active){ GRT_SAFE_ASPRINTF(&s_title, "z%s%c", GRT_SRC_M_NAME_ABBR[i], GRT_ZRT_CODES[c]); - GRT_NC_CHECK(nc_def_var(ncid, s_title, GRT_NC_MYREAL, ndims, dimids, &uiz_varids[i][c])); + NC_CHECK(nc_def_var(ncid, s_title, NC_REAL, ndims, dimids, &uiz_varids[i][c])); GRT_SAFE_ASPRINTF(&s_title, "r%s%c", GRT_SRC_M_NAME_ABBR[i], GRT_ZRT_CODES[c]); - GRT_NC_CHECK(nc_def_var(ncid, s_title, GRT_NC_MYREAL, ndims, dimids, &uir_varids[i][c])); + NC_CHECK(nc_def_var(ncid, s_title, NC_REAL, ndims, dimids, &uir_varids[i][c])); } } GRT_SAFE_FREE_PTR(s_title); } // 结束定义模式 - GRT_NC_CHECK(nc_enddef(ncid)); + NC_CHECK(nc_enddef(ncid)); // 写入数据 - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_put_var) (ncid, x_varid, Ctrl->X.xs)); - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_put_var) (ncid, y_varid, Ctrl->Y.ys)); - MYREAL *tmpdata = (MYREAL *)calloc(Ctrl->nr, sizeof(MYREAL)); + NC_CHECK(NC_FUNC_REAL(nc_put_var) (ncid, x_varid, Ctrl->X.xs)); + NC_CHECK(NC_FUNC_REAL(nc_put_var) (ncid, y_varid, Ctrl->Y.ys)); + real_t *tmpdata = (real_t *)calloc(Ctrl->nr, sizeof(real_t)); for(int i=0; ie.active){ for(int ir=0; ir < Ctrl->nr; ++ir){ tmpdata[ir] = (-1) * sgn0 * grn_uiz[ir][i][c]; // 这里多乘的(-1)是因为对z的偏导,z需反向 } - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_put_var) (ncid, uiz_varids[i][c], tmpdata)); + NC_CHECK(NC_FUNC_REAL(nc_put_var) (ncid, uiz_varids[i][c], tmpdata)); for(int ir=0; ir < Ctrl->nr; ++ir){ tmpdata[ir] = sgn0 * grn_uir[ir][i][c]; // 这里多乘的(-1)是因为对z的偏导,z需反向 } - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_put_var) (ncid, uir_varids[i][c], tmpdata)); + NC_CHECK(NC_FUNC_REAL(nc_put_var) (ncid, uir_varids[i][c], tmpdata)); } } } GRT_SAFE_FREE_PTR(tmpdata); // 关闭文件 - GRT_NC_CHECK(nc_close(ncid)); + NC_CHECK(nc_close(ncid)); // 释放内存 diff --git a/pygrt/C_extension/src/static/grt_static_rotation.c b/pygrt/C_extension/src/static/grt_static_rotation.c index 2f4c0a8c..a1d1ea34 100644 --- a/pygrt/C_extension/src/static/grt_static_rotation.c +++ b/pygrt/C_extension/src/static/grt_static_rotation.c @@ -78,7 +78,7 @@ int static_rotation_main(int argc, char **argv){ int out_varids[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]; // 打开 nc 文件 - GRT_NC_CHECK(nc_open(s_ingrid, NC_WRITE, &in_ncid)); + NC_CHECK(nc_open(s_ingrid, NC_WRITE, &in_ncid)); // 输出分量格式,即是否需要旋转到ZNE bool rot2ZNE = false; @@ -86,7 +86,7 @@ int static_rotation_main(int argc, char **argv){ // 读入数据是否旋转到ZNE { MYINT rot2ZNE_int; - GRT_NC_CHECK(GRT_NC_FUNC_MYINT(nc_get_att) (in_ncid, NC_GLOBAL, "rot2ZNE", &rot2ZNE_int)); + NC_CHECK(GRT_NC_FUNC_MYINT(nc_get_att) (in_ncid, NC_GLOBAL, "rot2ZNE", &rot2ZNE_int)); rot2ZNE = !(rot2ZNE_int == 0); } @@ -95,43 +95,43 @@ int static_rotation_main(int argc, char **argv){ // 读入的数据是否有位移偏导 MYINT calc_upar; - GRT_NC_CHECK(GRT_NC_FUNC_MYINT(nc_get_att) (in_ncid, NC_GLOBAL, "calc_upar", &calc_upar)); + NC_CHECK(GRT_NC_FUNC_MYINT(nc_get_att) (in_ncid, NC_GLOBAL, "calc_upar", &calc_upar)); if(calc_upar == 0){ GRTRaiseError("[%s] Input grid didn't have displacement derivatives.", command); } // 读入坐标变量 dimid, varid size_t nx, ny; - GRT_NC_CHECK(nc_inq_dimid(in_ncid, "north", &in_x_dimid)); - GRT_NC_CHECK(nc_inq_dimlen(in_ncid, in_x_dimid, &nx)); - GRT_NC_CHECK(nc_inq_dimid(in_ncid, "east", &in_y_dimid)); - GRT_NC_CHECK(nc_inq_dimlen(in_ncid, in_y_dimid, &ny)); + NC_CHECK(nc_inq_dimid(in_ncid, "north", &in_x_dimid)); + NC_CHECK(nc_inq_dimlen(in_ncid, in_x_dimid, &nx)); + NC_CHECK(nc_inq_dimid(in_ncid, "east", &in_y_dimid)); + NC_CHECK(nc_inq_dimlen(in_ncid, in_y_dimid, &ny)); in_dimids[0] = in_x_dimid; in_dimids[1] = in_y_dimid; // 读取坐标变量 - MYREAL *xs = (MYREAL *)calloc(nx, sizeof(MYREAL)); - MYREAL *ys = (MYREAL *)calloc(ny, sizeof(MYREAL)); - GRT_NC_CHECK(nc_inq_varid(in_ncid, "north", &in_x_varid)); - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_get_var) (in_ncid, in_x_varid, xs)); - GRT_NC_CHECK(nc_inq_varid(in_ncid, "east", &in_y_varid)); - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_get_var) (in_ncid, in_y_varid, ys)); + real_t *xs = (real_t *)calloc(nx, sizeof(real_t)); + real_t *ys = (real_t *)calloc(ny, sizeof(real_t)); + NC_CHECK(nc_inq_varid(in_ncid, "north", &in_x_varid)); + NC_CHECK(NC_FUNC_REAL(nc_get_var) (in_ncid, in_x_varid, xs)); + NC_CHECK(nc_inq_varid(in_ncid, "east", &in_y_varid)); + NC_CHECK(NC_FUNC_REAL(nc_get_var) (in_ncid, in_y_varid, ys)); // 读入合成位移偏导 varid for(int c=0; c 0.5*(u_{z,r} - u_{r,z}) GRT_SAFE_ASPRINTF(&s_title, "rotation_%c%c", toupper(chs[c]), toupper(chs[c2])); - GRT_NC_CHECK(nc_def_var(in_ncid, s_title, GRT_NC_MYREAL, ndims, in_dimids, &out_varids[c2][c])); + NC_CHECK(nc_def_var(in_ncid, s_title, NC_REAL, ndims, in_dimids, &out_varids[c2][c])); } GRT_SAFE_FREE_PTR(s_title); } // 结束定义模式 - GRT_NC_CHECK(nc_enddef(in_ncid)); + NC_CHECK(nc_enddef(in_ncid)); // 总震中距数 size_t nr = nx * ny; // 先读入内存, - MYREAL *u[GRT_CHANNEL_NUM]; - MYREAL *upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]; + real_t *u[GRT_CHANNEL_NUM]; + real_t *upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]; // 计算结果 - MYREAL *res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]; + real_t *res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]; for(int c=0; ccm if(chs[c]=='R' && chs[c2]=='T'){ @@ -194,13 +194,13 @@ int static_rotation_main(int argc, char **argv){ // 写入 nc 文件 for(int c=0; c 0.5*(u_{z,r} + u_{r,z}) GRT_SAFE_ASPRINTF(&s_title, "strain_%c%c", toupper(chs[c]), toupper(chs[c2])); - GRT_NC_CHECK(nc_def_var(in_ncid, s_title, GRT_NC_MYREAL, ndims, in_dimids, &out_varids[c2][c])); + NC_CHECK(nc_def_var(in_ncid, s_title, NC_REAL, ndims, in_dimids, &out_varids[c2][c])); } GRT_SAFE_FREE_PTR(s_title); } // 结束定义模式 - GRT_NC_CHECK(nc_enddef(in_ncid)); + NC_CHECK(nc_enddef(in_ncid)); // 总震中距数 size_t nr = nx * ny; // 先读入内存, - MYREAL *u[GRT_CHANNEL_NUM]; - MYREAL *upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]; + real_t *u[GRT_CHANNEL_NUM]; + real_t *upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]; // 计算结果 - MYREAL *res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]; + real_t *res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]; for(int c=0; ccm if(chs[c]=='R' && chs[c2]=='T'){ @@ -194,13 +194,13 @@ int static_strain_main(int argc, char **argv){ // 写入 nc 文件 for(int c=0; c mu*(u_{z,r} + u_{r,z}) GRT_SAFE_ASPRINTF(&s_title, "stress_%c%c", toupper(chs[c]), toupper(chs[c2])); - GRT_NC_CHECK(nc_def_var(in_ncid, s_title, GRT_NC_MYREAL, ndims, in_dimids, &out_varids[c2][c])); + NC_CHECK(nc_def_var(in_ncid, s_title, NC_REAL, ndims, in_dimids, &out_varids[c2][c])); } GRT_SAFE_FREE_PTR(s_title); } // 结束定义模式 - GRT_NC_CHECK(nc_enddef(in_ncid)); + NC_CHECK(nc_enddef(in_ncid)); // 总震中距数 size_t nr = nx * ny; // 先读入内存, - MYREAL *u[GRT_CHANNEL_NUM]; - MYREAL *upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]; + real_t *u[GRT_CHANNEL_NUM]; + real_t *upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]; // 计算结果 - MYREAL *res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]; + real_t *res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]; for(int c=0; cG.s_ingrid, NC_NOWRITE, &in_ncid)); - GRT_NC_CHECK(nc_create(Ctrl->O.s_outgrid, NC_CLOBBER, &out_ncid)); + NC_CHECK(nc_open(Ctrl->G.s_ingrid, NC_NOWRITE, &in_ncid)); + NC_CHECK(nc_create(Ctrl->O.s_outgrid, NC_CLOBBER, &out_ncid)); // 读取全局属性,视情况计算 src_mu - MYREAL src_va=0.0, src_vb=0.0, src_rho=0.0, src_mu=0.0; - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_get_att) (in_ncid, NC_GLOBAL, "src_va", &src_va)); - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_get_att) (in_ncid, NC_GLOBAL, "src_vb", &src_vb)); - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_get_att) (in_ncid, NC_GLOBAL, "src_rho", &src_rho)); + real_t src_va=0.0, src_vb=0.0, src_rho=0.0, src_mu=0.0; + NC_CHECK(NC_FUNC_REAL(nc_get_att) (in_ncid, NC_GLOBAL, "src_va", &src_va)); + NC_CHECK(NC_FUNC_REAL(nc_get_att) (in_ncid, NC_GLOBAL, "src_vb", &src_vb)); + NC_CHECK(NC_FUNC_REAL(nc_get_att) (in_ncid, NC_GLOBAL, "src_rho", &src_rho)); src_mu = src_vb*src_vb*src_rho*1e10; if(Ctrl->S.mult_src_mu) Ctrl->S.M0 *= src_mu; // 读入的数据是否有位移偏导 MYINT calc_upar; - GRT_NC_CHECK(GRT_NC_FUNC_MYINT(nc_get_att) (in_ncid, NC_GLOBAL, "calc_upar", &calc_upar)); + NC_CHECK(GRT_NC_FUNC_MYINT(nc_get_att) (in_ncid, NC_GLOBAL, "calc_upar", &calc_upar)); if(Ctrl->e.active && calc_upar == 0){ GRTRaiseError("[%s] Input grid didn't have displacement derivatives, you can't set -e.", command); } // 复制属性 - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_put_att) (out_ncid, NC_GLOBAL, "src_va", GRT_NC_MYREAL, 1, &src_va)); - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_put_att) (out_ncid, NC_GLOBAL, "src_vb", GRT_NC_MYREAL, 1, &src_vb)); - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_put_att) (out_ncid, NC_GLOBAL, "src_rho", GRT_NC_MYREAL, 1, &src_rho)); - GRT_NC_CHECK(GRT_NC_FUNC_MYINT(nc_put_att) (out_ncid, NC_GLOBAL, "calc_upar", GRT_NC_MYINT, 1, &calc_upar)); + NC_CHECK(NC_FUNC_REAL(nc_put_att) (out_ncid, NC_GLOBAL, "src_va", NC_REAL, 1, &src_va)); + NC_CHECK(NC_FUNC_REAL(nc_put_att) (out_ncid, NC_GLOBAL, "src_vb", NC_REAL, 1, &src_vb)); + NC_CHECK(NC_FUNC_REAL(nc_put_att) (out_ncid, NC_GLOBAL, "src_rho", NC_REAL, 1, &src_rho)); + NC_CHECK(GRT_NC_FUNC_MYINT(nc_put_att) (out_ncid, NC_GLOBAL, "calc_upar", GRT_NC_MYINT, 1, &calc_upar)); { - MYREAL rcv_va=0.0, rcv_vb=0.0, rcv_rho=0.0; - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_get_att) (in_ncid, NC_GLOBAL, "rcv_va", &rcv_va)); - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_get_att) (in_ncid, NC_GLOBAL, "rcv_vb", &rcv_vb)); - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_get_att) (in_ncid, NC_GLOBAL, "rcv_rho", &rcv_rho)); - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_put_att) (out_ncid, NC_GLOBAL, "rcv_va", GRT_NC_MYREAL, 1, &rcv_va)); - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_put_att) (out_ncid, NC_GLOBAL, "rcv_vb", GRT_NC_MYREAL, 1, &rcv_vb)); - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_put_att) (out_ncid, NC_GLOBAL, "rcv_rho", GRT_NC_MYREAL, 1, &rcv_rho)); + real_t rcv_va=0.0, rcv_vb=0.0, rcv_rho=0.0; + NC_CHECK(NC_FUNC_REAL(nc_get_att) (in_ncid, NC_GLOBAL, "rcv_va", &rcv_va)); + NC_CHECK(NC_FUNC_REAL(nc_get_att) (in_ncid, NC_GLOBAL, "rcv_vb", &rcv_vb)); + NC_CHECK(NC_FUNC_REAL(nc_get_att) (in_ncid, NC_GLOBAL, "rcv_rho", &rcv_rho)); + NC_CHECK(NC_FUNC_REAL(nc_put_att) (out_ncid, NC_GLOBAL, "rcv_va", NC_REAL, 1, &rcv_va)); + NC_CHECK(NC_FUNC_REAL(nc_put_att) (out_ncid, NC_GLOBAL, "rcv_vb", NC_REAL, 1, &rcv_vb)); + NC_CHECK(NC_FUNC_REAL(nc_put_att) (out_ncid, NC_GLOBAL, "rcv_rho", NC_REAL, 1, &rcv_rho)); } // 是否旋转到ZNE记录到全局属性 { MYINT rot2ZNE_int = rot2ZNE; - GRT_NC_CHECK(GRT_NC_FUNC_MYINT(nc_put_att) (out_ncid, NC_GLOBAL, "rot2ZNE", GRT_NC_MYINT, 1, &rot2ZNE_int)); + NC_CHECK(GRT_NC_FUNC_MYINT(nc_put_att) (out_ncid, NC_GLOBAL, "rot2ZNE", GRT_NC_MYINT, 1, &rot2ZNE_int)); } // 震源类型写入全局属性 - GRT_NC_CHECK(nc_put_att_text(out_ncid, NC_GLOBAL, "computeType", strlen(Ctrl->s_computeType), Ctrl->s_computeType)); + NC_CHECK(nc_put_att_text(out_ncid, NC_GLOBAL, "computeType", strlen(Ctrl->s_computeType), Ctrl->s_computeType)); // 读入坐标变量 dimid, varid size_t nx, ny; - GRT_NC_CHECK(nc_inq_dimid(in_ncid, "north", &in_x_dimid)); - GRT_NC_CHECK(nc_inq_dimlen(in_ncid, in_x_dimid, &nx)); - GRT_NC_CHECK(nc_inq_dimid(in_ncid, "east", &in_y_dimid)); - GRT_NC_CHECK(nc_inq_dimlen(in_ncid, in_y_dimid, &ny)); + NC_CHECK(nc_inq_dimid(in_ncid, "north", &in_x_dimid)); + NC_CHECK(nc_inq_dimlen(in_ncid, in_x_dimid, &nx)); + NC_CHECK(nc_inq_dimid(in_ncid, "east", &in_y_dimid)); + NC_CHECK(nc_inq_dimlen(in_ncid, in_y_dimid, &ny)); // 写入坐标变量 dimid, varid - GRT_NC_CHECK(nc_def_dim(out_ncid, "north", nx, &out_x_dimid)); - GRT_NC_CHECK(nc_def_dim(out_ncid, "east", ny, &out_y_dimid)); - GRT_NC_CHECK(nc_def_var(out_ncid, "north", GRT_NC_MYREAL, 1, &out_x_dimid, &out_x_varid)); - GRT_NC_CHECK(nc_def_var(out_ncid, "east", GRT_NC_MYREAL, 1, &out_y_dimid, &out_y_varid)); + NC_CHECK(nc_def_dim(out_ncid, "north", nx, &out_x_dimid)); + NC_CHECK(nc_def_dim(out_ncid, "east", ny, &out_y_dimid)); + NC_CHECK(nc_def_var(out_ncid, "north", NC_REAL, 1, &out_x_dimid, &out_x_varid)); + NC_CHECK(nc_def_var(out_ncid, "east", NC_REAL, 1, &out_y_dimid, &out_y_varid)); out_dimids[0] = out_x_dimid; out_dimids[1] = out_y_dimid; @@ -387,14 +387,14 @@ int static_syn_main(int argc, char **argv){ if(modr==0 && GRT_ZRT_CODES[c]=='T') continue; GRT_SAFE_ASPRINTF(&s_title, "%s%c", GRT_SRC_M_NAME_ABBR[i], GRT_ZRT_CODES[c]); - GRT_NC_CHECK(nc_inq_varid(in_ncid, s_title, &in_u_varids[i][c])); + NC_CHECK(nc_inq_varid(in_ncid, s_title, &in_u_varids[i][c])); // 位移偏导 if(Ctrl->e.active){ GRT_SAFE_ASPRINTF(&s_title, "z%s%c", GRT_SRC_M_NAME_ABBR[i], GRT_ZRT_CODES[c]); - GRT_NC_CHECK(nc_inq_varid(in_ncid, s_title, &in_uiz_varids[i][c])); + NC_CHECK(nc_inq_varid(in_ncid, s_title, &in_uiz_varids[i][c])); GRT_SAFE_ASPRINTF(&s_title, "r%s%c", GRT_SRC_M_NAME_ABBR[i], GRT_ZRT_CODES[c]); - GRT_NC_CHECK(nc_inq_varid(in_ncid, s_title, &in_uir_varids[i][c])); + NC_CHECK(nc_inq_varid(in_ncid, s_title, &in_uir_varids[i][c])); } } GRT_SAFE_FREE_PTR(s_title); @@ -404,83 +404,83 @@ int static_syn_main(int argc, char **argv){ for(int c=0; ce.active){ for(int c2=0; c2e.active){ - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_get_var) (in_ncid, in_uiz_varids[i][c], uiz[i][c])); - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_get_var) (in_ncid, in_uir_varids[i][c], uir[i][c])); + NC_CHECK(NC_FUNC_REAL(nc_get_var) (in_ncid, in_uiz_varids[i][c], uiz[i][c])); + NC_CHECK(NC_FUNC_REAL(nc_get_var) (in_ncid, in_uir_varids[i][c], uir[i][c])); } } } // 最终计算的结果 - MYREAL (*syn)[GRT_CHANNEL_NUM] = (MYREAL (*)[GRT_CHANNEL_NUM])calloc(nr, sizeof(MYREAL)*GRT_CHANNEL_NUM); - MYREAL (*syn_upar)[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM] = (MYREAL (*)[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM])calloc(nr, sizeof(MYREAL)*GRT_CHANNEL_NUM*GRT_CHANNEL_NUM); + real_t (*syn)[GRT_CHANNEL_NUM] = (real_t (*)[GRT_CHANNEL_NUM])calloc(nr, sizeof(real_t)*GRT_CHANNEL_NUM); + real_t (*syn_upar)[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM] = (real_t (*)[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM])calloc(nr, sizeof(real_t)*GRT_CHANNEL_NUM*GRT_CHANNEL_NUM); // 每个点逐个处理 for(size_t ix=0; ix < nx; ++ix){ - MYREAL x = xs[ix]; + real_t x = xs[ix]; for(size_t iy=0; iy < ny; ++iy){ - MYREAL y = ys[iy]; + real_t y = ys[iy]; size_t ir = iy + ix*ny; // 方位角 - MYREAL azrad = atan2(y, x); + real_t azrad = atan2(y, x); // 震中距 - MYREAL dist = GRT_MAX(sqrt(x*x + y*y), GRT_MIN_DISTANCE); + real_t dist = GRT_MAX(sqrt(x*x + y*y), GRT_MIN_DISTANCE); // 计算和位移相关量的种类(1-位移,2-ui_z,3-ui_r,4-ui_t) int calcUTypes = (Ctrl->e.active)? 4 : 1; - MYREAL upar_scale = 1.0; + real_t upar_scale = 1.0; - MYREAL *(*up)[GRT_CHANNEL_NUM]; // 使用对应类型的格林函数 - MYREAL tmpsyn[GRT_CHANNEL_NUM]; + real_t *(*up)[GRT_CHANNEL_NUM]; // 使用对应类型的格林函数 + real_t tmpsyn[GRT_CHANNEL_NUM]; for(int ityp=0; itypsrcRadi, Ctrl->computeType, (ityp > 0) && GRT_ZRT_CODES[ityp-1]=='T', Ctrl->S.M0, upar_scale, azrad, Ctrl->mchn); @@ -545,12 +545,12 @@ int static_syn_main(int argc, char **argv){ } // 写入 nc 文件 - MYREAL *tmpdata = (MYREAL *)calloc(nr, sizeof(MYREAL)); + real_t *tmpdata = (real_t *)calloc(nr, sizeof(real_t)); for(int c=0; ce.active){ @@ -558,7 +558,7 @@ int static_syn_main(int argc, char **argv){ for(size_t ir=0; ir < nr; ++ir){ tmpdata[ir] = syn_upar[ir][c2][c]; } - GRT_NC_CHECK(GRT_NC_FUNC_MYREAL(nc_put_var) (out_ncid, out_syn_upar_varids[c2][c], tmpdata)); + NC_CHECK(NC_FUNC_REAL(nc_put_var) (out_ncid, out_syn_upar_varids[c2][c], tmpdata)); } } } @@ -566,8 +566,8 @@ int static_syn_main(int argc, char **argv){ // 关闭文件 - GRT_NC_CHECK(nc_close(in_ncid)); - GRT_NC_CHECK(nc_close(out_ncid)); + NC_CHECK(nc_close(in_ncid)); + NC_CHECK(nc_close(out_ncid)); // 释放内存 for(int i=0; idepsrc - mod1d->deprcv), GRT_MIN_DEPTH_GAP_SRC_RCV); // hs=max(震源和台站深度差,1.0) + const real_t hs = GRT_MAX(fabs(mod1d->depsrc - mod1d->deprcv), GRT_MIN_DEPTH_GAP_SRC_RCV); // hs=max(震源和台站深度差,1.0) // 乘相应系数 k0 *= PI/hs; if(vmin_ref < 0.0) keps = -1.0; // 若使用峰谷平均法,则不使用keps进行收敛判断 - MYREAL k=0.0; + real_t k=0.0; bool useFIM = (filonLength > 0.0) || (safilonTol > 0.0) ; // 是否使用Filon积分(包括自适应Filon) - const MYREAL dk=fabs(PI2/(Length*rmax)); // 波数积分间隔 - const MYREAL filondk = (filonLength > 0.0) ? PI2/(filonLength*rmax) : 0.0; // Filon积分间隔 - const MYREAL filonK = filonCut/rmax; // 波数积分和Filon积分的分割点 + const real_t dk=fabs(PI2/(Length*rmax)); // 波数积分间隔 + const real_t filondk = (filonLength > 0.0) ? PI2/(filonLength*rmax) : 0.0; // Filon积分间隔 + const real_t filonK = filonCut/rmax; // 波数积分和Filon积分的分割点 - const MYREAL kmax = k0; + const real_t kmax = k0; // 求和 sum F(ki,w)Jm(ki*r)ki // 关于形状详见int_Pk()函数内的注释 - MYCOMPLEX (*sum_J)[GRT_SRC_M_NUM][GRT_INTEG_NUM] = (MYCOMPLEX(*)[GRT_SRC_M_NUM][GRT_INTEG_NUM])calloc(nr, sizeof(*sum_J)); - MYCOMPLEX (*sum_uiz_J)[GRT_SRC_M_NUM][GRT_INTEG_NUM] = (calc_upar)? (MYCOMPLEX(*)[GRT_SRC_M_NUM][GRT_INTEG_NUM])calloc(nr, sizeof(*sum_uiz_J)) : NULL; - MYCOMPLEX (*sum_uir_J)[GRT_SRC_M_NUM][GRT_INTEG_NUM] = (calc_upar)? (MYCOMPLEX(*)[GRT_SRC_M_NUM][GRT_INTEG_NUM])calloc(nr, sizeof(*sum_uir_J)) : NULL; + cplx_t (*sum_J)[GRT_SRC_M_NUM][GRT_INTEG_NUM] = (cplx_t(*)[GRT_SRC_M_NUM][GRT_INTEG_NUM])calloc(nr, sizeof(*sum_J)); + cplx_t (*sum_uiz_J)[GRT_SRC_M_NUM][GRT_INTEG_NUM] = (calc_upar)? (cplx_t(*)[GRT_SRC_M_NUM][GRT_INTEG_NUM])calloc(nr, sizeof(*sum_uiz_J)) : NULL; + cplx_t (*sum_uir_J)[GRT_SRC_M_NUM][GRT_INTEG_NUM] = (calc_upar)? (cplx_t(*)[GRT_SRC_M_NUM][GRT_INTEG_NUM])calloc(nr, sizeof(*sum_uir_J)) : NULL; // 是否要输出积分过程文件 bool needfstats = (statsstr!=NULL); @@ -185,8 +185,8 @@ void grt_integ_static_grn( - MYCOMPLEX src_mu = mod1d->mu[mod1d->isrc]; - MYCOMPLEX fac = dk * 1.0/(4.0*PI * src_mu); + cplx_t src_mu = mod1d->mu[mod1d->isrc]; + cplx_t fac = dk * 1.0/(4.0*PI * src_mu); // 将积分结果记录到浮点数数组中 recordin_GRN(nr, fac, sum_J, grn); diff --git a/pygrt/C_extension/src/static/static_layer.c b/pygrt/C_extension/src/static/static_layer.c index 3f84cafb..a727ecc5 100644 --- a/pygrt/C_extension/src/static/static_layer.c +++ b/pygrt/C_extension/src/static/static_layer.c @@ -27,7 +27,7 @@ void grt_static_topfree_RU(const GRT_MODEL1D *mod1d, RT_MATRIX *M) { - MYCOMPLEX delta1 = mod1d->delta[0]; + cplx_t delta1 = mod1d->delta[0]; // 公式(6.3.12) M->RU[0][0] = M->RU[1][1] = 0.0; M->RU[0][1] = -delta1; @@ -35,10 +35,10 @@ void grt_static_topfree_RU(const GRT_MODEL1D *mod1d, RT_MATRIX *M) M->RUL = 1.0; } -void grt_static_wave2qwv_REV_PSV(const GRT_MODEL1D *mod1d, const MYCOMPLEX R[2][2], MYCOMPLEX R_EV[2][2]) +void grt_static_wave2qwv_REV_PSV(const GRT_MODEL1D *mod1d, const cplx_t R[2][2], cplx_t R_EV[2][2]) { - MYCOMPLEX D11[2][2] = {{1.0, -1.0}, {1.0, 1.0}}; - MYCOMPLEX D12[2][2] = {{1.0, -1.0}, {-1.0, -1.0}}; + cplx_t D11[2][2] = {{1.0, -1.0}, {1.0, 1.0}}; + cplx_t D12[2][2] = {{1.0, -1.0}, {-1.0, -1.0}}; // 公式(6.3.35,37) if(mod1d->ircvup){// 震源更深 @@ -50,23 +50,23 @@ void grt_static_wave2qwv_REV_PSV(const GRT_MODEL1D *mod1d, const MYCOMPLEX R[2][ } } -void grt_static_wave2qwv_REV_SH(MYCOMPLEX RL, MYCOMPLEX *R_EVL) +void grt_static_wave2qwv_REV_SH(cplx_t RL, cplx_t *R_EVL) { *R_EVL = (1.0 + (RL)); } void grt_static_wave2qwv_z_REV_PSV( const GRT_MODEL1D *mod1d, - const MYCOMPLEX R[2][2], MYCOMPLEX R_EV[2][2]) + const cplx_t R[2][2], cplx_t R_EV[2][2]) { - MYREAL k = mod1d->k; + real_t k = mod1d->k; MYINT ircv = mod1d->ircv; - MYCOMPLEX delta1 = mod1d->delta[ircv]; + cplx_t delta1 = mod1d->delta[ircv]; // 新推导公式 - MYCOMPLEX kd2 = 2.0*k*delta1; - MYCOMPLEX D11[2][2] = {{k, -k-kd2}, {k, k-kd2}}; - MYCOMPLEX D12[2][2] = {{-k, k+kd2}, {k, k-kd2}}; + cplx_t kd2 = 2.0*k*delta1; + cplx_t D11[2][2] = {{k, -k-kd2}, {k, k-kd2}}; + cplx_t D12[2][2] = {{-k, k+kd2}, {k, k-kd2}}; if(mod1d->ircvup){// 震源更深 grt_cmat2x2_mul(D12, R, R_EV); grt_cmat2x2_add(D11, R_EV, R_EV); @@ -76,9 +76,9 @@ void grt_static_wave2qwv_z_REV_PSV( } } -void grt_static_wave2qwv_z_REV_SH(const GRT_MODEL1D *mod1d, MYCOMPLEX RL, MYCOMPLEX *R_EVL) +void grt_static_wave2qwv_z_REV_SH(const GRT_MODEL1D *mod1d, cplx_t RL, cplx_t *R_EVL) { - MYREAL k = mod1d->k; + real_t k = mod1d->k; // 新推导公式 if(mod1d->ircvup){// 震源更深 *R_EVL = (1.0 - (RL))*k; @@ -90,19 +90,19 @@ void grt_static_wave2qwv_z_REV_SH(const GRT_MODEL1D *mod1d, MYCOMPLEX RL, MYCOMP void grt_static_RT_matrix_PSV(const GRT_MODEL1D *mod1d, const MYINT iy, RT_MATRIX *M) { - MODEL_2LAYS_ATTRIB(MYCOMPLEX, mu); - MODEL_2LAYS_ATTRIB(MYCOMPLEX, delta); - MYREAL thk = mod1d->Thk[iy-1]; - MYREAL k = mod1d->k; + MODEL_2LAYS_ATTRIB(cplx_t, mu); + MODEL_2LAYS_ATTRIB(cplx_t, delta); + real_t thk = mod1d->Thk[iy-1]; + real_t k = mod1d->k; // 公式(6.3.18) - MYCOMPLEX dmu = mu1 - mu2; - MYCOMPLEX A112 = mu1*delta1 + mu2; - MYCOMPLEX A221 = mu2*delta2 + mu1; - MYCOMPLEX B = mu1*delta1 - mu2*delta2; - MYCOMPLEX del11 = delta1*delta1; - MYREAL k2 = k*k; - MYREAL thk2 = thk*thk; + cplx_t dmu = mu1 - mu2; + cplx_t A112 = mu1*delta1 + mu2; + cplx_t A221 = mu2*delta2 + mu1; + cplx_t B = mu1*delta1 - mu2*delta2; + cplx_t del11 = delta1*delta1; + real_t k2 = k*k; + real_t thk2 = thk*thk; // Reflection //------------------ RD ----------------------------------- @@ -132,11 +132,11 @@ void grt_static_RT_matrix_PSV(const GRT_MODEL1D *mod1d, const MYINT iy, RT_MATRI void grt_static_RT_matrix_SH(const GRT_MODEL1D *mod1d, const MYINT iy, RT_MATRIX *M) { - MODEL_2LAYS_ATTRIB(MYCOMPLEX, mu); + MODEL_2LAYS_ATTRIB(cplx_t, mu); // 公式(6.3.18) - MYCOMPLEX dmu = mu1 - mu2; - MYCOMPLEX amu = mu1 + mu2; + cplx_t dmu = mu1 - mu2; + cplx_t amu = mu1 + mu2; // Reflection M->RDL = dmu/amu; @@ -150,10 +150,10 @@ void grt_static_RT_matrix_SH(const GRT_MODEL1D *mod1d, const MYINT iy, RT_MATRIX void grt_static_delay_RT_matrix(const GRT_MODEL1D *mod1d, const MYINT iy, RT_MATRIX *M) { - MYREAL thk = mod1d->Thk[iy-1]; - MYREAL k = mod1d->k; + real_t thk = mod1d->Thk[iy-1]; + real_t k = mod1d->k; - MYCOMPLEX ex, ex2; + cplx_t ex, ex2; ex = exp(- k*thk); ex2 = ex * ex; diff --git a/pygrt/C_extension/src/static/static_propagate.c b/pygrt/C_extension/src/static/static_propagate.c index a9a854c9..197fc2ca 100644 --- a/pygrt/C_extension/src/static/static_propagate.c +++ b/pygrt/C_extension/src/static/static_propagate.c @@ -26,8 +26,8 @@ void grt_static_kernel( - const GRT_MODEL1D *mod1d, MYCOMPLEX QWV[GRT_SRC_M_NUM][GRT_QWV_NUM], - bool calc_uiz, MYCOMPLEX QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]) + const GRT_MODEL1D *mod1d, cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM], + bool calc_uiz, cplx_t QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]) { // 初始化qwv为0 for(MYINT i=0; in; ++iy){ // 因为n>=3, 故一定会进入该循环 @@ -111,13 +111,13 @@ void grt_static_kernel( // 计算震源系数 - MYCOMPLEX src_coef_PSV[GRT_SRC_M_NUM][GRT_QWV_NUM-1][2] = {0}; - MYCOMPLEX src_coef_SH[GRT_SRC_M_NUM][2] = {0}; + cplx_t src_coef_PSV[GRT_SRC_M_NUM][GRT_QWV_NUM-1][2] = {0}; + cplx_t src_coef_SH[GRT_SRC_M_NUM][2] = {0}; grt_static_source_coef_PSV(mod1d, src_coef_PSV); grt_static_source_coef_SH(mod1d, src_coef_SH); // 临时中转矩阵 (temperary) - MYCOMPLEX tmpR2[2][2], tmp2x2[2][2], tmpRL, tmp2x2_uiz[2][2], tmpRL_uiz; + cplx_t tmpR2[2][2], tmp2x2[2][2], tmpRL, tmp2x2_uiz[2][2], tmpRL_uiz; // 递推RU_FA grt_static_topfree_RU(mod1d, M_top); diff --git a/pygrt/C_extension/src/static/static_source.c b/pygrt/C_extension/src/static/static_source.c index 1fe59c25..b44804fe 100644 --- a/pygrt/C_extension/src/static/static_source.c +++ b/pygrt/C_extension/src/static/static_source.c @@ -16,7 +16,7 @@ #include "grt/static/static_source.h" #include "grt/common/const.h" -void grt_static_source_coef_PSV(const GRT_MODEL1D *mod1d, MYCOMPLEX coef[GRT_SRC_M_NUM][GRT_QWV_NUM-1][2]) +void grt_static_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM-1][2]) { // 先全部赋0 for(MYINT i=0; iisrc; - MYCOMPLEX delta = mod1d->delta[isrc]; - MYREAL k = mod1d->k; + cplx_t delta = mod1d->delta[isrc]; + real_t k = mod1d->k; - MYCOMPLEX tmp; - MYCOMPLEX A = 1.0+delta; + cplx_t tmp; + cplx_t A = 1.0+delta; // 爆炸源 coef[0][0][0] = tmp = (delta-1.0)/A; coef[0][0][1] = tmp; @@ -58,9 +58,9 @@ void grt_static_source_coef_PSV(const GRT_MODEL1D *mod1d, MYCOMPLEX coef[GRT_SRC } -void grt_static_source_coef_SH(const GRT_MODEL1D *mod1d, MYCOMPLEX coef[GRT_SRC_M_NUM][2]) +void grt_static_source_coef_SH(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][2]) { - MYREAL k = mod1d->k; + real_t k = mod1d->k; // 先全部赋0 for(MYINT i=0; i 1.0/Vel[i]) pmax0 = 1.0/Vel[i]; } // 初始化一些迭代变量 - MYREAL pmin=0.0, pmax=pmax0; - MYREAL p; - MYREAL s, c, v, h; - MYREAL x = 0.0; - MYREAL t = 0.0; - MYREAL tint = 0.0; - MYREAL dxdp = 0.0; + real_t pmin=0.0, pmax=pmax0; + real_t p; + real_t s, c, v, h; + real_t x = 0.0; + real_t t = 0.0; + real_t tint = 0.0; + real_t dxdp = 0.0; for(int iter=0; iter 0.0){ // 存在射线向上的基本条件 - MYREAL v, p, h, c; - MYREAL sumt, sumx; + real_t v, p, h, c; + real_t sumt, sumx; bool badrefrac = false; // 找到透射位置 for(int m=imin-1; m>=0; --m){ @@ -328,8 +328,8 @@ MYREAL grt_compute_travt1d( //------------------- 向下出射的射线,考虑透射 ----------------- // 找到透射位置 for(int m=imax+1; mR.rs = (MYREAL*)realloc(Ctrl->R.rs, sizeof(MYREAL)*(Ctrl->R.nr)); + Ctrl->R.rs = (real_t*)realloc(Ctrl->R.rs, sizeof(real_t)*(Ctrl->R.nr)); for(MYINT i=0; iR.nr; ++i){ Ctrl->R.rs[i] = atof(Ctrl->R.s_rs[i]); if(Ctrl->R.rs[i] < 0.0){ diff --git a/pygrt/c_structures.py b/pygrt/c_structures.py index edf7b9c3..da200a18 100755 --- a/pygrt/c_structures.py +++ b/pygrt/c_structures.py @@ -13,7 +13,6 @@ from ctypes import * __all__ = [ - "USE_FLOAT", "CHANNEL_NUM", "QWV_NUM", "INTEG_NUM", @@ -35,7 +34,6 @@ ] -USE_FLOAT = False CHANNEL_NUM = 3 QWV_NUM = 3 INTEG_NUM = 4 @@ -47,12 +45,12 @@ qwvchs = ['q', 'w', 'v'] -NPCT_REAL_TYPE = 'f4' if USE_FLOAT else 'f8' -NPCT_CMPLX_TYPE = f'c{int(NPCT_REAL_TYPE[1:])*2}' +NPCT_REAL_TYPE = 'f8' +NPCT_CMPLX_TYPE = 'c16' -REAL = c_float if USE_FLOAT else c_double +REAL = c_double CPLX = REAL*2 PREAL = POINTER(REAL) PCPLX = POINTER(CPLX)