diff --git a/pygrt/C_extension/include/common/bessel.h b/pygrt/C_extension/include/common/bessel.h index aca4e03c..02dd4991 100755 --- a/pygrt/C_extension/include/common/bessel.h +++ b/pygrt/C_extension/include/common/bessel.h @@ -19,4 +19,16 @@ * @param bj2 (out) \f$ J_2(x) \f$ * */ -void bessel012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2); \ No newline at end of file +void bessel012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2); + + +/** + * 计算Bessel函数的一阶导数 \f$ J_m^'(x), m=0,1,2 \f$ + * + * @param x 自变量 + * @param bj0 (inout) 传入\f$ J_0(x) \f$, 返回\f$ J_0^'(x) \f$ + * @param bj1 (inout) 传入\f$ J_1(x) \f$, 返回\f$ J_1^'(x) \f$ + * @param bj2 (inout) 传入\f$ J_2(x) \f$, 返回\f$ J_2^'(x) \f$ + * + */ +void besselp012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2); \ No newline at end of file diff --git a/pygrt/C_extension/include/common/const.h b/pygrt/C_extension/include/common/const.h index 746d913a..bf0f75be 100755 --- a/pygrt/C_extension/include/common/const.h +++ b/pygrt/C_extension/include/common/const.h @@ -68,6 +68,9 @@ typedef int MYINT; ///< 整数 #define PI2 6.2831853f ///< \f$ 2\pi \f$ #define HALFPI 1.5707963f ///< \f$ \frac{\pi}{2} \f$ #define QUARTERPI 0.78539816f ///< \f$ \frac{\pi}{4} \f$ + #define THREEQUARTERPI 2.35619450f ///< \f$ \frac{3\pi}{4} \f$ + #define FIVEQUARTERPI 3.92699082f ///< \f$ \frac{5\pi}{4} \f$ + #define SEVENQUARTERPI 5.49778714f ///< \f$ \frac{7\pi}{4} \f$ #define INV_SQRT_TWO 0.70710678f ///< \f$ \frac{1}{\sqrt{2}} \f$ #else @@ -102,6 +105,9 @@ typedef int MYINT; ///< 整数 #define PI2 6.283185307179586 ///< \f$ 2\pi \f$ #define HALFPI 1.5707963267948966 ///< \f$ \frac{\pi}{2} \f$ #define QUARTERPI 0.7853981633974483 ///< \f$ \frac{\pi}{4} \f$ + #define THREEQUARTERPI 2.356194490192345 ///< \f$ \frac{3\pi}{4} \f$ + #define FIVEQUARTERPI 3.9269908169872414 ///< \f$ \frac{5\pi}{4} \f$ + #define SEVENQUARTERPI 5.497787143782138 ///< \f$ \frac{7\pi}{4} \f$ #define INV_SQRT_TWO 0.7071067811865475 ///< \f$ \frac{1}{\sqrt{2}} \f$ #endif diff --git a/pygrt/C_extension/include/dynamic/dwm.h b/pygrt/C_extension/include/dynamic/dwm.h index 33701c90..813e79e8 100644 --- a/pygrt/C_extension/include/dynamic/dwm.h +++ b/pygrt/C_extension/include/dynamic/dwm.h @@ -32,6 +32,17 @@ * @param sum_VF_J[nr][3][4] (out)垂直力源 * @param sum_HF_J[nr][3][4] (out)水平力源 * @param sum_DC_J[nr][3][4] (out)双力偶源 + * + * @param calc_upar (in)是否计算位移u的空间导数 + * @param sum_EXP_uiz_J[nr][3][4] (out)爆炸源产生的ui_z(位移u对坐标z的偏导)，下同 + * @param sum_VF_uiz_J[nr][3][4] (out)垂直力源 + * @param sum_HF_uiz_J[nr][3][4] (out)水平力源 + * @param sum_DC_uiz_J[nr][3][4] (out)双力偶源 + * @param sum_EXP_uir_J[nr][3][4] (out)爆炸源产生的ui_r(位移u对坐标r的偏导)，下同 + * @param sum_VF_uir_J[nr][3][4] (out)垂直力源 + * @param sum_HF_uir_J[nr][3][4] (out)水平力源 + * @param sum_DC_uir_J[nr][3][4] (out)双力偶源 + * * @param fstats[nr]) (out)不同震中距的格林函数积分过程文件 * * @return k 积分截至时的波数 @@ -41,4 +52,9 @@ MYREAL discrete_integ( MYINT nr, MYREAL *rs, MYCOMPLEX sum_EXP_J[nr][3][4], MYCOMPLEX sum_VF_J[nr][3][4], MYCOMPLEX sum_HF_J[nr][3][4], MYCOMPLEX sum_DC_J[nr][3][4], + bool calc_upar, + MYCOMPLEX sum_EXP_uiz_J[nr][3][4], MYCOMPLEX sum_VF_uiz_J[nr][3][4], + MYCOMPLEX sum_HF_uiz_J[nr][3][4], MYCOMPLEX sum_DC_uiz_J[nr][3][4], + MYCOMPLEX sum_EXP_uir_J[nr][3][4], MYCOMPLEX sum_VF_uir_J[nr][3][4], + MYCOMPLEX sum_HF_uir_J[nr][3][4], MYCOMPLEX sum_DC_uir_J[nr][3][4], FILE *(fstats[nr])); diff --git a/pygrt/C_extension/include/dynamic/fim.h b/pygrt/C_extension/include/dynamic/fim.h index d2a2f92d..9be2ebfd 100755 --- a/pygrt/C_extension/include/dynamic/fim.h +++ b/pygrt/C_extension/include/dynamic/fim.h @@ -38,6 +38,17 @@ * @param sum_VF_J[nr][3][4] (out)垂直力源 * @param sum_HF_J[nr][3][4] (out)水平力源 * @param sum_DC_J[nr][3][4] (out)双力偶源 + * + * @param calc_upar (in)是否计算位移u的空间导数 + * @param sum_EXP_uiz_J[nr][3][4] (out)爆炸源产生的ui_z(位移u对坐标z的偏导)，下同 + * @param sum_VF_uiz_J[nr][3][4] (out)垂直力源 + * @param sum_HF_uiz_J[nr][3][4] (out)水平力源 + * @param sum_DC_uiz_J[nr][3][4] (out)双力偶源 + * @param sum_EXP_uir_J[nr][3][4] (out)爆炸源产生的ui_r(位移u对坐标r的偏导)，下同 + * @param sum_VF_uir_J[nr][3][4] (out)垂直力源 + * @param sum_HF_uir_J[nr][3][4] (out)水平力源 + * @param sum_DC_uir_J[nr][3][4] (out)双力偶源 + * * @param fstats[nr] (out)不同震中距的格林函数积分过程文件 * * @return k 积分截至时的波数 @@ -47,6 +58,11 @@ MYREAL linear_filon_integ( MYINT nr, MYREAL *rs, MYCOMPLEX sum_EXP_J[nr][3][4], MYCOMPLEX sum_VF_J[nr][3][4], MYCOMPLEX sum_HF_J[nr][3][4], MYCOMPLEX sum_DC_J[nr][3][4], + bool calc_upar, + MYCOMPLEX sum_EXP_uiz_J[nr][3][4], MYCOMPLEX sum_VF_uiz_J[nr][3][4], + MYCOMPLEX sum_HF_uiz_J[nr][3][4], MYCOMPLEX sum_DC_uiz_J[nr][3][4], + MYCOMPLEX sum_EXP_uir_J[nr][3][4], MYCOMPLEX sum_VF_uir_J[nr][3][4], + MYCOMPLEX sum_HF_uir_J[nr][3][4], MYCOMPLEX sum_DC_uir_J[nr][3][4], FILE *(fstats[nr])); @@ -63,6 +79,7 @@ MYREAL linear_filon_integ( * @param VF_qwv[3][3] (in)垂直力源核函数 * @param HF_qwv[3][3] (in)水平力源核函数 * @param DC_qwv[3][3] (in)双力偶源核函数 + * @param calc_uir (in)是否计算ui_r（位移u对坐标r的偏导） * @param EXP_J[3][4] (out)爆炸源，该dk区间内的积分值，下同 * @param VF_J[3][4] (out)垂直力源 * @param HF_J[3][4] (out)水平力源 @@ -73,6 +90,7 @@ void int_Pk_filon( MYREAL k, MYREAL r, const MYCOMPLEX EXP_qwv[3][3], const MYCOMPLEX VF_qwv[3][3], const MYCOMPLEX HF_qwv[3][3], const MYCOMPLEX DC_qwv[3][3], + bool calc_uir, MYCOMPLEX EXP_J[3][4], MYCOMPLEX VF_J[3][4], MYCOMPLEX HF_J[3][4], MYCOMPLEX DC_J[3][4] ); diff --git a/pygrt/C_extension/include/dynamic/grt.h b/pygrt/C_extension/include/dynamic/grt.h index 2ed060cc..42c51bfc 100755 --- a/pygrt/C_extension/include/dynamic/grt.h +++ b/pygrt/C_extension/include/dynamic/grt.h @@ -49,6 +49,21 @@ void set_num_threads(int num_threads); * @param DSgrn[nr][3] (out)复数数组，90度倾滑的Z、R、T分量频谱结果 * @param SSgrn[nr][3] (out)复数数组，90度走滑的Z、R、T分量频谱结果 * + * @param calc_upar (in)是否计算位移u的空间导数 + * @param EXPgrn_uiz[nr][2] (out)复数数组，爆炸源产生的ui_z(位移u对坐标z的偏导)的Z、R分量频谱结果，下同 + * @param VFgrn_uiz[nr][2] (out)复数数组，垂直力源的Z、R分量频谱结果 + * @param HFgrn_uiz[nr][3] (out)复数数组，水平力源的Z、R、T分量频谱结果 + * @param DDgrn_uiz[nr][2] (out)复数数组，45度倾滑的Z、R分量频谱结果 + * @param DSgrn_uiz[nr][3] (out)复数数组，90度倾滑的Z、R、T分量频谱结果 + * @param SSgrn_uiz[nr][3] (out)复数数组，90度走滑的Z、R、T分量频谱结果 + * @param EXPgrn_uir[nr][2] (out)复数数组，爆炸源产生的ui_r(位移u对坐标r的偏导)的Z、R分量频谱结果，下同 + * @param VFgrn_uir[nr][2] (out)复数数组，垂直力源的Z、R分量频谱结果 + * @param HFgrn_uir[nr][3] (out)复数数组，水平力源的Z、R、T分量频谱结果 + * @param DDgrn_uir[nr][2] (out)复数数组，45度倾滑的Z、R分量频谱结果 + * @param DSgrn_uir[nr][3] (out)复数数组，90度倾滑的Z、R、T分量频谱结果 + * @param SSgrn_uir[nr][3] (out)复数数组，90度走滑的Z、R、T分量频谱结果 + * + * * @param statsstr (in) 积分结果输出目录 * @param nstatsidxs (in) 输出积分结果的特定频点的个数 * @param statsidxs (in) 特定频点的索引值 @@ -68,6 +83,20 @@ void integ_grn_spec_in_C( MYCOMPLEX *DScplx[nr][3], // DSZ, DSR, DST 的实部和虚部 [DS: 90-dip slip] MYCOMPLEX *SScplx[nr][3], // SSZ, SSR, SST 的实部和虚部 [SS: strike slip] + bool calc_upar, + MYCOMPLEX *EXPcplx_uiz[nr][2], // EXZ, EXR 的实部和虚部 + MYCOMPLEX *VFcplx_uiz[nr][2], // VFZ, VFR 的实部和虚部 + MYCOMPLEX *HFcplx_uiz[nr][3], // HFZ, HFR, HFT 的实部和虚部 + MYCOMPLEX *DDcplx_uiz[nr][2], // DDZ, DDR 的实部和虚部 [DD: 45-dip slip] + MYCOMPLEX *DScplx_uiz[nr][3], // DSZ, DSR, DST 的实部和虚部 [DS: 90-dip slip] + MYCOMPLEX *SScplx_uiz[nr][3], // SSZ, SSR, SST 的实部和虚部 [SS: strike slip] + MYCOMPLEX *EXPcplx_uir[nr][2], // EXZ, EXR 的实部和虚部 + MYCOMPLEX *VFcplx_uir[nr][2], // VFZ, VFR 的实部和虚部 + MYCOMPLEX *HFcplx_uir[nr][3], // HFZ, HFR, HFT 的实部和虚部 + MYCOMPLEX *DDcplx_uir[nr][2], // DDZ, DDR 的实部和虚部 [DD: 45-dip slip] + MYCOMPLEX *DScplx_uir[nr][3], // DSZ, DSR, DST 的实部和虚部 [DS: 90-dip slip] + MYCOMPLEX *SScplx_uir[nr][3], // SSZ, SSR, SST 的实部和虚部 [SS: strike slip] + const char *statsstr, // 积分结果输出 MYINT nstatsidxs, // 仅输出特定频点 MYINT *statsidxs @@ -99,6 +128,20 @@ void integ_grn_spec_in_C( * @param DSgrn[nr][3] (out)`GRN` 结构体指针，90度倾滑的Z、R、T分量频谱结果 * @param SSgrn[nr][3] (out)`GRN` 结构体指针，90度走滑的Z、R、T分量频谱结果 * + * @param calc_upar (in)是否计算位移u的空间导数 + * @param EXPgrn_uiz[nr][2] (out)`GRN` 结构体指针，爆炸源产生的ui_z(位移u对坐标z的偏导)的Z、R分量频谱结果，下同 + * @param VFgrn_uiz[nr][2] (out)`GRN` 结构体指针，垂直力源的Z、R分量频谱结果 + * @param HFgrn_uiz[nr][3] (out)`GRN` 结构体指针，水平力源的Z、R、T分量频谱结果 + * @param DDgrn_uiz[nr][2] (out)`GRN` 结构体指针，45度倾滑的Z、R分量频谱结果 + * @param DSgrn_uiz[nr][3] (out)`GRN` 结构体指针，90度倾滑的Z、R、T分量频谱结果 + * @param SSgrn_uiz[nr][3] (out)`GRN` 结构体指针，90度走滑的Z、R、T分量频谱结果 + * @param EXPgrn_uir[nr][2] (out)`GRN` 结构体指针，爆炸源产生的ui_r(位移u对坐标r的偏导)的Z、R分量频谱结果，下同 + * @param VFgrn_uir[nr][2] (out)`GRN` 结构体指针，垂直力源的Z、R分量频谱结果 + * @param HFgrn_uir[nr][3] (out)`GRN` 结构体指针，水平力源的Z、R、T分量频谱结果 + * @param DDgrn_uir[nr][2] (out)`GRN` 结构体指针，45度倾滑的Z、R分量频谱结果 + * @param DSgrn_uir[nr][3] (out)`GRN` 结构体指针，90度倾滑的Z、R、T分量频谱结果 + * @param SSgrn_uir[nr][3] (out)`GRN` 结构体指针，90度走滑的Z、R、T分量频谱结果 + * * @param statsstr (in) 积分结果输出目录 * @param nstatsidxs (in) 输出积分结果的特定频点的个数 * @param statsidxs (in) 特定频点的索引值 @@ -118,6 +161,20 @@ void integ_grn_spec( GRN *DSgrn[nr][3], // DSZ, DSR, DST 的实部和虚部 [DS: 90-dip slip] GRN *SSgrn[nr][3], // SSZ, SSR, SST 的实部和虚部 [SS: strike slip] + bool calc_upar, + GRN *EXPgrn_uiz[nr][2], // EXZ, EXR 的实部和虚部 + GRN *VFgrn_uiz[nr][2], // VFZ, VFR 的实部和虚部 + GRN *HFgrn_uiz[nr][3], // HFZ, HFR, HFT 的实部和虚部 + GRN *DDgrn_uiz[nr][2], // DDZ, DDR 的实部和虚部 [DD: 45-dip slip] + GRN *DSgrn_uiz[nr][3], // DSZ, DSR, DST 的实部和虚部 [DS: 90-dip slip] + GRN *SSgrn_uiz[nr][3], // SSZ, SSR, SST 的实部和虚部 [SS: strike slip] + GRN *EXPgrn_uir[nr][2], // EXZ, EXR 的实部和虚部 + GRN *VFgrn_uir[nr][2], // VFZ, VFR 的实部和虚部 + GRN *HFgrn_uir[nr][3], // HFZ, HFR, HFT 的实部和虚部 + GRN *DDgrn_uir[nr][2], // DDZ, DDR 的实部和虚部 [DD: 45-dip slip] + GRN *DSgrn_uir[nr][3], // DSZ, DSR, DST 的实部和虚部 [DS: 90-dip slip] + GRN *SSgrn_uir[nr][3], // SSZ, SSR, SST 的实部和虚部 [SS: strike slip] + const char *statsstr, // 积分结果输出 MYINT nstatsidxs, // 仅输出特定频点 MYINT *statsidxs diff --git a/pygrt/C_extension/include/dynamic/layer.h b/pygrt/C_extension/include/dynamic/layer.h index 0b1a0482..5fbdaa26 100755 --- a/pygrt/C_extension/include/dynamic/layer.h +++ b/pygrt/C_extension/include/dynamic/layer.h @@ -50,6 +50,13 @@ void calc_R_EV( MYCOMPLEX R_EV[2][2], MYCOMPLEX *R_EVL); +void calc_uiz_R_EV( + MYCOMPLEX xa_rcv, MYCOMPLEX xb_rcv, bool ircvup, + MYREAL k, + const MYCOMPLEX R[2][2], MYCOMPLEX RL, + MYCOMPLEX R_EV[2][2], MYCOMPLEX *R_EVL); + + /** * 计算界面的反射系数RD/RDL/RU/RUL, 透射系数TD/TDL/TU/TUL, 包括时间延迟因子， * 后缀L表示SH波的系数, 其余表示P-SV波的系数, 根据公式(5.4.14)和(5.4.31)计算系数 diff --git a/pygrt/C_extension/include/dynamic/propagate.h b/pygrt/C_extension/include/dynamic/propagate.h index 6d9ec9e1..4002a77b 100755 --- a/pygrt/C_extension/include/dynamic/propagate.h +++ b/pygrt/C_extension/include/dynamic/propagate.h @@ -88,11 +88,18 @@ * @param VF_qwv[3][3] (out)垂直力源核函数 * @param HF_qwv[3][3] (out)水平力源核函数 * @param DC_qwv[3][3] (out)双力偶源核函数 + * @param calc_uiz (in)是否计算ui_z（位移u对坐标z的偏导） + * @param EXP_uiz_qwv[3][3] (out)爆炸源产生的ui_z的核函数，下同 + * @param VF_uiz_qwv[3][3] (out)垂直力源核函数 + * @param HF_uiz_qwv[3][3] (out)水平力源核函数 + * @param DC_uiz_qwv[3][3] (out)双力偶源核函数 * */ void kernel( const MODEL1D *mod1d, MYCOMPLEX omega, MYREAL k, - MYCOMPLEX EXP_qwv[3][3], MYCOMPLEX VF_qwv[3][3], MYCOMPLEX HF_qwv[3][3], MYCOMPLEX DC_qwv[3][3]); + MYCOMPLEX EXP_qwv[3][3], MYCOMPLEX VF_qwv[3][3], MYCOMPLEX HF_qwv[3][3], MYCOMPLEX DC_qwv[3][3], + bool calc_uiz, + MYCOMPLEX EXP_uiz_qwv[3][3], MYCOMPLEX VF_uiz_qwv[3][3], MYCOMPLEX HF_uiz_qwv[3][3], MYCOMPLEX DC_uiz_qwv[3][3]); @@ -108,6 +115,7 @@ void kernel( * @param VF_qwv[3][3] (in)垂直力源核函数 * @param HF_qwv[3][3] (in)水平力源核函数 * @param DC_qwv[3][3] (in)双力偶源核函数 + * @param calc_uir (in)是否计算ui_r（位移u对坐标r的偏导） * @param EXP_J[3][4] (out)爆炸源，该dk区间内的积分值，下同 * @param VF_J[3][4] (out)垂直力源 * @param HF_J[3][4] (out)水平力源 @@ -118,6 +126,7 @@ void int_Pk( MYREAL k, MYREAL r, const MYCOMPLEX EXP_qwv[3][3], const MYCOMPLEX VF_qwv[3][3], const MYCOMPLEX HF_qwv[3][3], const MYCOMPLEX DC_qwv[3][3], + bool calc_uir, MYCOMPLEX EXP_J[3][4], MYCOMPLEX VF_J[3][4], MYCOMPLEX HF_J[3][4], MYCOMPLEX DC_J[3][4] ); diff --git a/pygrt/C_extension/include/dynamic/ptam.h b/pygrt/C_extension/include/dynamic/ptam.h index 44ba22c6..94a8116d 100755 --- a/pygrt/C_extension/include/dynamic/ptam.h +++ b/pygrt/C_extension/include/dynamic/ptam.h @@ -40,6 +40,17 @@ * @param sum_VF_J0[nr][3][4] (out)垂直力源 * @param sum_HF_J0[nr][3][4] (out)水平力源 * @param sum_DC_J0[nr][3][4] (out)双力偶源 + * + * @param calc_upar (in)是否计算位移u的空间导数 + * @param sum_EXP_uiz_J0[nr][3][4] (out)爆炸源产生的ui_z(位移u对坐标z的偏导)，下同 + * @param sum_VF_uiz_J0[nr][3][4] (out)垂直力源 + * @param sum_HF_uiz_J0[nr][3][4] (out)水平力源 + * @param sum_DC_uiz_J0[nr][3][4] (out)双力偶源 + * @param sum_EXP_uir_J0[nr][3][4] (out)爆炸源产生的ui_r(位移u对坐标r的偏导)，下同 + * @param sum_VF_uir_J0[nr][3][4] (out)垂直力源 + * @param sum_HF_uir_J0[nr][3][4] (out)水平力源 + * @param sum_DC_uir_J0[nr][3][4] (out)双力偶源 + * * @param fstats[nr] (out)波数积分过程文件指针 * @param ptam_fstats[nr] (out)峰谷平均法过程文件指针 * @@ -50,6 +61,11 @@ void PTA_method( MYINT nr, MYREAL *rs, MYCOMPLEX sum_EXP_J0[nr][3][4], MYCOMPLEX sum_VF_J0[nr][3][4], MYCOMPLEX sum_HF_J0[nr][3][4], MYCOMPLEX sum_DC_J0[nr][3][4], + bool calc_upar, + MYCOMPLEX sum_EXP_uiz_J0[nr][3][4], MYCOMPLEX sum_VF_uiz_J0[nr][3][4], + MYCOMPLEX sum_HF_uiz_J0[nr][3][4], MYCOMPLEX sum_DC_uiz_J0[nr][3][4], + MYCOMPLEX sum_EXP_uir_J0[nr][3][4], MYCOMPLEX sum_VF_uir_J0[nr][3][4], + MYCOMPLEX sum_HF_uir_J0[nr][3][4], MYCOMPLEX sum_DC_uir_J0[nr][3][4], FILE *(fstats[nr]), FILE *(ptam_fstats[nr])); diff --git a/pygrt/C_extension/src/common/bessel.c b/pygrt/C_extension/src/common/bessel.c index 125e82d2..37195ab9 100755 --- a/pygrt/C_extension/src/common/bessel.c +++ b/pygrt/C_extension/src/common/bessel.c @@ -16,3 +16,11 @@ void bessel012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2){ *bj2 = JN(2, x); } +void besselp012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2){ + MYREAL j0=*bj0; + MYREAL j1=*bj1; + MYREAL j2=*bj2; + *bj0 = -j1; + *bj1 = j0 - RONE/x * j1; + *bj2 = j1 - RTWO/x * j2; +} diff --git a/pygrt/C_extension/src/dynamic/dwm.c b/pygrt/C_extension/src/dynamic/dwm.c index c7a7b083..18819dab 100644 --- a/pygrt/C_extension/src/dynamic/dwm.c +++ b/pygrt/C_extension/src/dynamic/dwm.c @@ -27,6 +27,11 @@ MYREAL discrete_integ( MYINT nr, MYREAL *rs, MYCOMPLEX sum_EXP_J[nr][3][4], MYCOMPLEX sum_VF_J[nr][3][4], MYCOMPLEX sum_HF_J[nr][3][4], MYCOMPLEX sum_DC_J[nr][3][4], + bool calc_upar, + MYCOMPLEX sum_EXP_uiz_J[nr][3][4], MYCOMPLEX sum_VF_uiz_J[nr][3][4], + MYCOMPLEX sum_HF_uiz_J[nr][3][4], MYCOMPLEX sum_DC_uiz_J[nr][3][4], + MYCOMPLEX sum_EXP_uir_J[nr][3][4], MYCOMPLEX sum_VF_uir_J[nr][3][4], + MYCOMPLEX sum_HF_uir_J[nr][3][4], MYCOMPLEX sum_DC_uir_J[nr][3][4], FILE *(fstats[nr])) { MYCOMPLEX EXP_J[3][4], VF_J[3][4], HF_J[3][4], DC_J[3][4]; @@ -40,6 +45,12 @@ MYREAL discrete_integ( MYCOMPLEX (*pVF_qwv)[3] = (sum_VF_J!=NULL)? VF_qwv : NULL; MYCOMPLEX (*pHF_qwv)[3] = (sum_HF_J!=NULL)? HF_qwv : NULL; MYCOMPLEX (*pDC_qwv)[3] = (sum_DC_J!=NULL)? DC_qwv : NULL; + + MYCOMPLEX EXP_uiz_qwv[3][3], VF_uiz_qwv[3][3], HF_uiz_qwv[3][3], DC_uiz_qwv[3][3]; + MYCOMPLEX (*pEXP_uiz_qwv)[3] = (sum_EXP_uiz_J!=NULL)? EXP_uiz_qwv : NULL; + MYCOMPLEX (*pVF_uiz_qwv)[3] = (sum_VF_uiz_J!=NULL)? VF_uiz_qwv : NULL; + MYCOMPLEX (*pHF_uiz_qwv)[3] = (sum_HF_uiz_J!=NULL)? HF_uiz_qwv : NULL; + MYCOMPLEX (*pDC_uiz_qwv)[3] = (sum_DC_uiz_J!=NULL)? DC_uiz_qwv : NULL; MYREAL k = 0.0; MYINT ik = 0; @@ -61,7 +72,8 @@ MYREAL discrete_integ( // printf("w=%15.5e, ik=%d\n", CREAL(omega), ik); // 计算核函数 F(k, w) - kernel(mod1d, omega, k, pEXP_qwv, pVF_qwv, pHF_qwv, pDC_qwv); + kernel(mod1d, omega, k, pEXP_qwv, pVF_qwv, pHF_qwv, pDC_qwv, + calc_upar, pEXP_uiz_qwv, pVF_uiz_qwv, pHF_uiz_qwv, pDC_uiz_qwv); // 震中距rs循环 @@ -77,8 +89,8 @@ MYREAL discrete_integ( // 计算被积函数一项 F(k,w)Jm(kr)k int_Pk(k, rs[ir], - pEXP_qwv, pVF_qwv, pHF_qwv, pDC_qwv, - EXP_J, VF_J, HF_J, DC_J); + pEXP_qwv, pVF_qwv, pHF_qwv, pDC_qwv, false, + EXP_J, VF_J, HF_J, DC_J); // 记录积分结果 if(fstats[ir]!=NULL){ @@ -114,7 +126,42 @@ MYREAL discrete_integ( iendk = iendkrs[ir] = false; } - + + // ---------------- 位移空间导数，EXP_J, VF_J, HF_J, DC_J数组重复利用 -------------------------- + if(calc_upar){ + // ------------------------------- ui_z ----------------------------------- + // 计算被积函数一项 F(k,w)Jm(kr)k + int_Pk(k, rs[ir], + pEXP_uiz_qwv, pVF_uiz_qwv, pHF_uiz_qwv, pDC_uiz_qwv, false, + EXP_J, VF_J, HF_J, DC_J); + + // keps不参与计算位移空间导数的积分，背后逻辑认为u收敛，则uiz也收敛 + for(MYINT m=0; m<3; ++m){ + for(MYINT v=0; v<4; ++v){ + if(sum_EXP_uiz_J!=NULL) sum_EXP_uiz_J[ir][m][v] += EXP_J[m][v]; + if(sum_VF_uiz_J!=NULL) sum_VF_uiz_J[ir][m][v] += VF_J[m][v]; + if(sum_HF_uiz_J!=NULL) sum_HF_uiz_J[ir][m][v] += HF_J[m][v]; + if(sum_DC_uiz_J!=NULL) sum_DC_uiz_J[ir][m][v] += DC_J[m][v]; + } + } + + + // ------------------------------- ui_r ----------------------------------- + // 计算被积函数一项 F(k,w)Jm(kr)k + int_Pk(k, rs[ir], + pEXP_qwv, pVF_qwv, pHF_qwv, pDC_qwv, true, + EXP_J, VF_J, HF_J, DC_J); + + // keps不参与计算位移空间导数的积分，背后逻辑认为u收敛，则uir也收敛 + for(MYINT m=0; m<3; ++m){ + for(MYINT v=0; v<4; ++v){ + if(sum_EXP_uir_J!=NULL) sum_EXP_uir_J[ir][m][v] += EXP_J[m][v]; + if(sum_VF_uir_J!=NULL) sum_VF_uir_J[ir][m][v] += VF_J[m][v]; + if(sum_HF_uir_J!=NULL) sum_HF_uir_J[ir][m][v] += HF_J[m][v]; + if(sum_DC_uir_J!=NULL) sum_DC_uir_J[ir][m][v] += DC_J[m][v]; + } + } + } // END if calc_upar } // END rs loop diff --git a/pygrt/C_extension/src/dynamic/fim.c b/pygrt/C_extension/src/dynamic/fim.c index f591fb21..ceac3d4d 100755 --- a/pygrt/C_extension/src/dynamic/fim.c +++ b/pygrt/C_extension/src/dynamic/fim.c @@ -30,6 +30,11 @@ MYREAL linear_filon_integ( MYINT nr, MYREAL *rs, MYCOMPLEX sum_EXP_J[nr][3][4], MYCOMPLEX sum_VF_J[nr][3][4], MYCOMPLEX sum_HF_J[nr][3][4], MYCOMPLEX sum_DC_J[nr][3][4], + bool calc_upar, + MYCOMPLEX sum_EXP_uiz_J[nr][3][4], MYCOMPLEX sum_VF_uiz_J[nr][3][4], + MYCOMPLEX sum_HF_uiz_J[nr][3][4], MYCOMPLEX sum_DC_uiz_J[nr][3][4], + MYCOMPLEX sum_EXP_uir_J[nr][3][4], MYCOMPLEX sum_VF_uir_J[nr][3][4], + MYCOMPLEX sum_HF_uir_J[nr][3][4], MYCOMPLEX sum_DC_uir_J[nr][3][4], FILE *(fstats[nr])) { for(MYINT ir=0; ir kmax) break; // 计算核函数 F(k, w) - kernel(mod1d, omega, k, pEXP_qwv, pVF_qwv, pHF_qwv, pDC_qwv); + kernel(mod1d, omega, k, pEXP_qwv, pVF_qwv, pHF_qwv, pDC_qwv, + calc_upar, pEXP_uiz_qwv, pVF_uiz_qwv, pHF_uiz_qwv, pDC_uiz_qwv); // 震中距rs循环 iendk = true; @@ -113,7 +121,7 @@ MYREAL linear_filon_integ( // F(k, w)*Jm(kr)k 的近似公式 int_Pk_filon( k, rs[ir], - pEXP_qwv, pVF_qwv, pHF_qwv, pDC_qwv, + pEXP_qwv, pVF_qwv, pHF_qwv, pDC_qwv, false, EXP_J, VF_J, HF_J, DC_J); // 记录积分结果 @@ -151,6 +159,42 @@ MYREAL linear_filon_integ( } + // ---------------- 位移空间导数，EXP_J, VF_J, HF_J, DC_J数组重复利用 -------------------------- + if(calc_upar){ + // ------------------------------- ui_z ----------------------------------- + // 计算被积函数一项 F(k,w)Jm(kr)k + int_Pk_filon(k, rs[ir], + pEXP_uiz_qwv, pVF_uiz_qwv, pHF_uiz_qwv, pDC_uiz_qwv, false, + EXP_J, VF_J, HF_J, DC_J); + + // keps不参与计算位移空间导数的积分，背后逻辑认为u收敛，则uiz也收敛 + for(MYINT m=0; m<3; ++m){ + for(MYINT v=0; v<4; ++v){ + if(sum_EXP_uiz_J!=NULL) sum_EXP_uiz_J[ir][m][v] += EXP_J[m][v]; + if(sum_VF_uiz_J!=NULL) sum_VF_uiz_J[ir][m][v] += VF_J[m][v]; + if(sum_HF_uiz_J!=NULL) sum_HF_uiz_J[ir][m][v] += HF_J[m][v]; + if(sum_DC_uiz_J!=NULL) sum_DC_uiz_J[ir][m][v] += DC_J[m][v]; + } + } + + + // ------------------------------- ui_r ----------------------------------- + // 计算被积函数一项 F(k,w)Jm(kr)k + int_Pk_filon(k, rs[ir], + pEXP_qwv, pVF_qwv, pHF_qwv, pDC_qwv, true, + EXP_J, VF_J, HF_J, DC_J); + + // keps不参与计算位移空间导数的积分，背后逻辑认为u收敛，则uir也收敛 + for(MYINT m=0; m<3; ++m){ + for(MYINT v=0; v<4; ++v){ + if(sum_EXP_uir_J!=NULL) sum_EXP_uir_J[ir][m][v] += EXP_J[m][v]; + if(sum_VF_uir_J!=NULL) sum_VF_uir_J[ir][m][v] += VF_J[m][v]; + if(sum_HF_uir_J!=NULL) sum_HF_uir_J[ir][m][v] += HF_J[m][v]; + if(sum_DC_uir_J!=NULL) sum_DC_uir_J[ir][m][v] += DC_J[m][v]; + } + } + } // END if calc_upar + } // end rs loop @@ -162,39 +206,28 @@ MYREAL linear_filon_integ( // 乘上系数 for(MYINT ir=0; irRe = (MYREAL*)calloc(nf, sizeof(MYREAL)); EXPgrn[ir][i]->Im = (MYREAL*)calloc(nf, sizeof(MYREAL)); } + if(EXPcplx_uiz) { + EXPgrn_uiz[ir][i] = (GRN*)calloc(1, sizeof(GRN)); + EXPgrn_uiz[ir][i]->Re = (MYREAL*)calloc(nf, sizeof(MYREAL)); + EXPgrn_uiz[ir][i]->Im = (MYREAL*)calloc(nf, sizeof(MYREAL)); + } + if(EXPcplx_uir) { + EXPgrn_uir[ir][i] = (GRN*)calloc(1, sizeof(GRN)); + EXPgrn_uir[ir][i]->Re = (MYREAL*)calloc(nf, sizeof(MYREAL)); + EXPgrn_uir[ir][i]->Im = (MYREAL*)calloc(nf, sizeof(MYREAL)); + } + // if(VFcplx) { VFgrn[ir][i] = (GRN*)calloc(1, sizeof(GRN)); VFgrn[ir][i]->Re = (MYREAL*)calloc(nf, sizeof(MYREAL)); VFgrn[ir][i]->Im = (MYREAL*)calloc(nf, sizeof(MYREAL)); } + if(VFcplx_uiz) { + VFgrn_uiz[ir][i] = (GRN*)calloc(1, sizeof(GRN)); + VFgrn_uiz[ir][i]->Re = (MYREAL*)calloc(nf, sizeof(MYREAL)); + VFgrn_uiz[ir][i]->Im = (MYREAL*)calloc(nf, sizeof(MYREAL)); + } + if(VFcplx_uir) { + VFgrn_uir[ir][i] = (GRN*)calloc(1, sizeof(GRN)); + VFgrn_uir[ir][i]->Re = (MYREAL*)calloc(nf, sizeof(MYREAL)); + VFgrn_uir[ir][i]->Im = (MYREAL*)calloc(nf, sizeof(MYREAL)); + } + // if(DDcplx) { DDgrn[ir][i] = (GRN*)calloc(1, sizeof(GRN)); DDgrn[ir][i]->Re = (MYREAL*)calloc(nf, sizeof(MYREAL)); DDgrn[ir][i]->Im = (MYREAL*)calloc(nf, sizeof(MYREAL)); } + if(DDcplx_uiz) { + DDgrn_uiz[ir][i] = (GRN*)calloc(1, sizeof(GRN)); + DDgrn_uiz[ir][i]->Re = (MYREAL*)calloc(nf, sizeof(MYREAL)); + DDgrn_uiz[ir][i]->Im = (MYREAL*)calloc(nf, sizeof(MYREAL)); + } + if(DDcplx_uir) { + DDgrn_uir[ir][i] = (GRN*)calloc(1, sizeof(GRN)); + DDgrn_uir[ir][i]->Re = (MYREAL*)calloc(nf, sizeof(MYREAL)); + DDgrn_uir[ir][i]->Im = (MYREAL*)calloc(nf, sizeof(MYREAL)); + } } + // if(HFcplx) { HFgrn[ir][i] = (GRN*)calloc(1, sizeof(GRN)); HFgrn[ir][i]->Re = (MYREAL*)calloc(nf, sizeof(MYREAL)); HFgrn[ir][i]->Im = (MYREAL*)calloc(nf, sizeof(MYREAL)); } + if(HFcplx_uiz) { + HFgrn_uiz[ir][i] = (GRN*)calloc(1, sizeof(GRN)); + HFgrn_uiz[ir][i]->Re = (MYREAL*)calloc(nf, sizeof(MYREAL)); + HFgrn_uiz[ir][i]->Im = (MYREAL*)calloc(nf, sizeof(MYREAL)); + } + if(HFcplx_uir) { + HFgrn_uir[ir][i] = (GRN*)calloc(1, sizeof(GRN)); + HFgrn_uir[ir][i]->Re = (MYREAL*)calloc(nf, sizeof(MYREAL)); + HFgrn_uir[ir][i]->Im = (MYREAL*)calloc(nf, sizeof(MYREAL)); + } + // if(DScplx) { DSgrn[ir][i] = (GRN*)calloc(1, sizeof(GRN)); DSgrn[ir][i]->Re = (MYREAL*)calloc(nf, sizeof(MYREAL)); DSgrn[ir][i]->Im = (MYREAL*)calloc(nf, sizeof(MYREAL)); } + if(DScplx_uiz) { + DSgrn_uiz[ir][i] = (GRN*)calloc(1, sizeof(GRN)); + DSgrn_uiz[ir][i]->Re = (MYREAL*)calloc(nf, sizeof(MYREAL)); + DSgrn_uiz[ir][i]->Im = (MYREAL*)calloc(nf, sizeof(MYREAL)); + } + if(DScplx_uir) { + DSgrn_uir[ir][i] = (GRN*)calloc(1, sizeof(GRN)); + DSgrn_uir[ir][i]->Re = (MYREAL*)calloc(nf, sizeof(MYREAL)); + DSgrn_uir[ir][i]->Im = (MYREAL*)calloc(nf, sizeof(MYREAL)); + } + // if(SScplx) { SSgrn[ir][i] = (GRN*)calloc(1, sizeof(GRN)); SSgrn[ir][i]->Re = (MYREAL*)calloc(nf, sizeof(MYREAL)); SSgrn[ir][i]->Im = (MYREAL*)calloc(nf, sizeof(MYREAL)); } + if(SScplx_uiz) { + SSgrn_uiz[ir][i] = (GRN*)calloc(1, sizeof(GRN)); + SSgrn_uiz[ir][i]->Re = (MYREAL*)calloc(nf, sizeof(MYREAL)); + SSgrn_uiz[ir][i]->Im = (MYREAL*)calloc(nf, sizeof(MYREAL)); + } + if(SScplx_uir) { + SSgrn_uir[ir][i] = (GRN*)calloc(1, sizeof(GRN)); + SSgrn_uir[ir][i]->Re = (MYREAL*)calloc(nf, sizeof(MYREAL)); + SSgrn_uir[ir][i]->Im = (MYREAL*)calloc(nf, sizeof(MYREAL)); + } } } @@ -112,6 +206,9 @@ void integ_grn_spec_in_C( pymod1d, nf1, nf2, nf, freqs, nr, rs, wI, vmin_ref, keps, ampk, iwk0, k0, Length, print_progressbar, EXPgrn, VFgrn, HFgrn, DDgrn, DSgrn, SSgrn, + calc_upar, + EXPgrn_uiz, VFgrn_uiz, HFgrn_uiz, DDgrn_uiz, DSgrn_uiz, SSgrn_uiz, + EXPgrn_uir, VFgrn_uir, HFgrn_uir, DDgrn_uir, DSgrn_uir, SSgrn_uir, statsstr, nstatsidxs, statsidxs ); //============================================================================== @@ -122,13 +219,31 @@ void integ_grn_spec_in_C( for(int i=0; i<3; ++i){ for(int n=nf1; n<=nf2; ++n){ if(i<2){ + // if(EXPcplx) EXPcplx[ir][i][n] = CMPLX(EXPgrn[ir][i]->Re[n], EXPgrn[ir][i]->Im[n]); + if(EXPcplx_uiz) EXPcplx_uiz[ir][i][n] = CMPLX(EXPgrn_uiz[ir][i]->Re[n], EXPgrn_uiz[ir][i]->Im[n]); + if(EXPcplx_uir) EXPcplx_uir[ir][i][n] = CMPLX(EXPgrn_uir[ir][i]->Re[n], EXPgrn_uir[ir][i]->Im[n]); + // if(VFcplx) VFcplx[ir][i][n] = CMPLX(VFgrn[ir][i]->Re[n], VFgrn[ir][i]->Im[n]); + if(VFcplx_uiz) VFcplx_uiz[ir][i][n] = CMPLX(VFgrn_uiz[ir][i]->Re[n], VFgrn_uiz[ir][i]->Im[n]); + if(VFcplx_uir) VFcplx_uir[ir][i][n] = CMPLX(VFgrn_uir[ir][i]->Re[n], VFgrn_uir[ir][i]->Im[n]); if(DDcplx) DDcplx[ir][i][n] = CMPLX(DDgrn[ir][i]->Re[n], DDgrn[ir][i]->Im[n]); + // + if(DDcplx_uiz) DDcplx_uiz[ir][i][n] = CMPLX(DDgrn_uiz[ir][i]->Re[n], DDgrn_uiz[ir][i]->Im[n]); + if(DDcplx_uir) DDcplx_uir[ir][i][n] = CMPLX(DDgrn_uir[ir][i]->Re[n], DDgrn_uir[ir][i]->Im[n]); } + // if(HFcplx) HFcplx[ir][i][n] = CMPLX(HFgrn[ir][i]->Re[n], HFgrn[ir][i]->Im[n]); + if(HFcplx_uiz) HFcplx_uiz[ir][i][n] = CMPLX(HFgrn_uiz[ir][i]->Re[n], HFgrn_uiz[ir][i]->Im[n]); + if(HFcplx_uir) HFcplx_uir[ir][i][n] = CMPLX(HFgrn_uir[ir][i]->Re[n], HFgrn_uir[ir][i]->Im[n]); + // if(DScplx) DScplx[ir][i][n] = CMPLX(DSgrn[ir][i]->Re[n], DSgrn[ir][i]->Im[n]); + if(DScplx_uiz) DScplx_uiz[ir][i][n] = CMPLX(DSgrn_uiz[ir][i]->Re[n], DSgrn_uiz[ir][i]->Im[n]); + if(DScplx_uir) DScplx_uir[ir][i][n] = CMPLX(DSgrn_uir[ir][i]->Re[n], DSgrn_uir[ir][i]->Im[n]); + // if(SScplx) SScplx[ir][i][n] = CMPLX(SSgrn[ir][i]->Re[n], SSgrn[ir][i]->Im[n]); + if(SScplx_uiz) SScplx_uiz[ir][i][n] = CMPLX(SSgrn_uiz[ir][i]->Re[n], SSgrn_uiz[ir][i]->Im[n]); + if(SScplx_uir) SScplx_uir[ir][i][n] = CMPLX(SSgrn_uir[ir][i]->Re[n], SSgrn_uir[ir][i]->Im[n]); } } } @@ -138,48 +253,226 @@ void integ_grn_spec_in_C( for(int ir=0; irRe); free(EXPgrn[ir][i]->Im); free(EXPgrn[ir][i]); } + if(EXPgrn_uiz) { + free(EXPgrn_uiz[ir][i]->Re); + free(EXPgrn_uiz[ir][i]->Im); + free(EXPgrn_uiz[ir][i]); + } + if(EXPgrn_uir) { + free(EXPgrn_uir[ir][i]->Re); + free(EXPgrn_uir[ir][i]->Im); + free(EXPgrn_uir[ir][i]); + } + // if(VFgrn) { free(VFgrn[ir][i]->Re); free(VFgrn[ir][i]->Im); free(VFgrn[ir][i]); } + if(VFgrn_uiz) { + free(VFgrn_uiz[ir][i]->Re); + free(VFgrn_uiz[ir][i]->Im); + free(VFgrn_uiz[ir][i]); + } + if(VFgrn_uir) { + free(VFgrn_uir[ir][i]->Re); + free(VFgrn_uir[ir][i]->Im); + free(VFgrn_uir[ir][i]); + } + // if(DDgrn) { free(DDgrn[ir][i]->Re); free(DDgrn[ir][i]->Im); free(DDgrn[ir][i]); } + if(DDgrn_uiz) { + free(DDgrn_uiz[ir][i]->Re); + free(DDgrn_uiz[ir][i]->Im); + free(DDgrn_uiz[ir][i]); + } + if(DDgrn_uir) { + free(DDgrn_uir[ir][i]->Re); + free(DDgrn_uir[ir][i]->Im); + free(DDgrn_uir[ir][i]); + } } + // if(HFcplx) { free(HFgrn[ir][i]->Re); free(HFgrn[ir][i]->Im); free(HFgrn[ir][i]); } + if(HFcplx_uiz) { + free(HFgrn_uiz[ir][i]->Re); + free(HFgrn_uiz[ir][i]->Im); + free(HFgrn_uiz[ir][i]); + } + if(HFcplx_uir) { + free(HFgrn_uir[ir][i]->Re); + free(HFgrn_uir[ir][i]->Im); + free(HFgrn_uir[ir][i]); + } + // if(DScplx) { free(DSgrn[ir][i]->Re); free(DSgrn[ir][i]->Im); free(DSgrn[ir][i]); } + if(DScplx_uiz) { + free(DSgrn_uiz[ir][i]->Re); + free(DSgrn_uiz[ir][i]->Im); + free(DSgrn_uiz[ir][i]); + } + if(DScplx_uir) { + free(DSgrn_uir[ir][i]->Re); + free(DSgrn_uir[ir][i]->Im); + free(DSgrn_uir[ir][i]); + } + // if(SScplx) { free(SSgrn[ir][i]->Re); free(SSgrn[ir][i]->Im); free(SSgrn[ir][i]); } + if(SScplx_uiz) { + free(SSgrn_uiz[ir][i]->Re); + free(SSgrn_uiz[ir][i]->Im); + free(SSgrn_uiz[ir][i]); + } + if(SScplx_uir) { + free(SSgrn_uir[ir][i]->Re); + free(SSgrn_uir[ir][i]->Im); + free(SSgrn_uir[ir][i]); + } } } if(EXPgrn) free(EXPgrn); + if(EXPgrn_uiz) free(EXPgrn_uiz); + if(EXPgrn_uir) free(EXPgrn_uir); if(VFgrn) free(VFgrn); + if(VFgrn_uiz) free(VFgrn_uiz); + if(VFgrn_uir) free(VFgrn_uir); if(HFgrn) free(HFgrn); + if(HFgrn_uiz) free(HFgrn_uiz); + if(HFgrn_uir) free(HFgrn_uir); if(DDgrn) free(DDgrn); + if(DDgrn_uiz) free(DDgrn_uiz); + if(DDgrn_uir) free(DDgrn_uir); if(DSgrn) free(DSgrn); + if(DSgrn_uiz) free(DSgrn_uiz); + if(DSgrn_uir) free(DSgrn_uir); if(SSgrn) free(SSgrn); + if(SSgrn_uiz) free(SSgrn_uiz); + if(SSgrn_uir) free(SSgrn_uir); +} + + +/** + * 将计算好的复数形式的积分结果记录到GRN结构体中 + * + * @param iw (in)当前频率索引值 + * @param nr (in)震中距个数 + * @param coef (in)统一系数 + * @param sum_EXP_J[nr][3][4] (in)爆炸源 + * @param sum_VF_J[nr][3][4] (in)垂直力源 + * @param sum_HF_J[nr][3][4] (in)水平力源 + * @param sum_DC_J[nr][3][4] (in)双力偶源 + * @param EXPgrn[nr][2] (out)`GRN` 结构体指针，爆炸源的Z、R分量频谱结果 + * @param VFgrn[nr][2] (out)`GRN` 结构体指针，垂直力源的Z、R分量频谱结果 + * @param HFgrn[nr][3] (out)`GRN` 结构体指针，水平力源的Z、R、T分量频谱结果 + * @param DDgrn[nr][2] (out)`GRN` 结构体指针，45度倾滑的Z、R分量频谱结果 + * @param DSgrn[nr][3] (out)`GRN` 结构体指针，90度倾滑的Z、R、T分量频谱结果 + * @param SSgrn[nr][3] (out)`GRN` 结构体指针，90度走滑的Z、R、T分量频谱结果 + */ +static void recordin_GRN( + MYINT iw, MYINT nr, MYCOMPLEX coef, + MYCOMPLEX sum_EXP_J[nr][3][4], MYCOMPLEX sum_VF_J[nr][3][4], + MYCOMPLEX sum_HF_J[nr][3][4], MYCOMPLEX sum_DC_J[nr][3][4], + GRN *EXPgrn[nr][2], GRN *VFgrn[nr][2], GRN *HFgrn[nr][3], + GRN *DDgrn[nr][2], GRN *DSgrn[nr][3], GRN *SSgrn[nr][3] +) +{ + // 局部变量，将某个频点的格林函数谱临时存放 + MYCOMPLEX (*tmp_EXP)[2] = (MYCOMPLEX(*)[2])calloc(nr, sizeof(*tmp_EXP)); + MYCOMPLEX (*tmp_VF)[2] = (MYCOMPLEX(*)[2])calloc(nr, sizeof(*tmp_VF)); + MYCOMPLEX (*tmp_HF)[3] = (MYCOMPLEX(*)[3])calloc(nr, sizeof(*tmp_HF)); + MYCOMPLEX (*tmp_DD)[2] = (MYCOMPLEX(*)[2])calloc(nr, sizeof(*tmp_DD)); + MYCOMPLEX (*tmp_DS)[3] = (MYCOMPLEX(*)[3])calloc(nr, sizeof(*tmp_DS)); + MYCOMPLEX (*tmp_SS)[3] = (MYCOMPLEX(*)[3])calloc(nr, sizeof(*tmp_SS)); + + for(MYINT ir=0; irRe[iw] = CREAL(mtmp); + EXPgrn[ir][ii]->Im[iw] = CIMAG(mtmp); + } + if(VFgrn!=NULL){ + mtmp = coef*tmp_VF[ir][ii]; // m=0 垂直力源 + VFgrn[ir][ii]->Re[iw] = CREAL(mtmp); + VFgrn[ir][ii]->Im[iw] = CIMAG(mtmp); + } + if(DDgrn!=NULL){ + mtmp = coef*tmp_DD[ir][ii]; // m=0 45-倾滑 + DDgrn[ir][ii]->Re[iw] = CREAL(mtmp); + DDgrn[ir][ii]->Im[iw] = CIMAG(mtmp); + } + } + if(HFgrn!=NULL){ + mtmp = coef*tmp_HF[ir][ii]; // m=1 水平力源 + HFgrn[ir][ii]->Re[iw] = CREAL(mtmp); + HFgrn[ir][ii]->Im[iw] = CIMAG(mtmp); + } + if(DSgrn!=NULL){ + mtmp = coef*tmp_DS[ir][ii]; // m=1 90-倾滑 + DSgrn[ir][ii]->Re[iw] = CREAL(mtmp); + DSgrn[ir][ii]->Im[iw] = CIMAG(mtmp); + } + if(SSgrn!=NULL){ + mtmp = coef*tmp_SS[ir][ii]; // m=2 走滑 + SSgrn[ir][ii]->Re[iw] = CREAL(mtmp); + SSgrn[ir][ii]->Im[iw] = CIMAG(mtmp); + } + + } + } + + free(tmp_EXP); + free(tmp_VF); + free(tmp_HF); + free(tmp_DD); + free(tmp_DS); + free(tmp_SS); } + void integ_grn_spec( PYMODEL1D *pymod1d, MYINT nf1, MYINT nf2, MYINT nf, MYREAL *freqs, MYINT nr, MYREAL *rs, MYREAL wI, @@ -194,6 +487,20 @@ void integ_grn_spec( GRN *DSgrn[nr][3], // DSZ, DSR, DST 的实部和虚部 [DS: 90-dip slip] GRN *SSgrn[nr][3], // SSZ, SSR, SST 的实部和虚部 [SS: strike slip] + bool calc_upar, + GRN *EXPgrn_uiz[nr][2], // EXZ, EXR 的实部和虚部 + GRN *VFgrn_uiz[nr][2], // VFZ, VFR 的实部和虚部 + GRN *HFgrn_uiz[nr][3], // HFZ, HFR, HFT 的实部和虚部 + GRN *DDgrn_uiz[nr][2], // DDZ, DDR 的实部和虚部 [DD: 45-dip slip] + GRN *DSgrn_uiz[nr][3], // DSZ, DSR, DST 的实部和虚部 [DS: 90-dip slip] + GRN *SSgrn_uiz[nr][3], // SSZ, SSR, SST 的实部和虚部 [SS: strike slip] + GRN *EXPgrn_uir[nr][2], // EXZ, EXR 的实部和虚部 + GRN *VFgrn_uir[nr][2], // VFZ, VFR 的实部和虚部 + GRN *HFgrn_uir[nr][3], // HFZ, HFR, HFT 的实部和虚部 + GRN *DDgrn_uir[nr][2], // DDZ, DDR 的实部和虚部 [DD: 45-dip slip] + GRN *DSgrn_uir[nr][3], // DSZ, DSR, DST 的实部和虚部 [DS: 90-dip slip] + GRN *SSgrn_uir[nr][3], // SSZ, SSR, SST 的实部和虚部 [SS: strike slip] + const char *statsstr, // 积分结果输出 MYINT nstatsidxs, // 仅输出特定频点 MYINT *statsidxs @@ -259,14 +566,7 @@ void integ_grn_spec( MYCOMPLEX omega = w - wI*I; // 复数频率 omega = w - i*wI MYCOMPLEX omega2_inv = RONE/omega; // 1/omega^2 omega2_inv = omega2_inv*omega2_inv; - - // 局部变量，将某个频点的格林函数谱临时存放 - MYCOMPLEX (*tmp_EXP)[2] = (MYCOMPLEX(*)[2])calloc(nr, sizeof(*tmp_EXP)); - MYCOMPLEX (*tmp_VF)[2] = (MYCOMPLEX(*)[2])calloc(nr, sizeof(*tmp_VF)); - MYCOMPLEX (*tmp_HF)[3] = (MYCOMPLEX(*)[3])calloc(nr, sizeof(*tmp_HF)); - MYCOMPLEX (*tmp_DD)[2] = (MYCOMPLEX(*)[2])calloc(nr, sizeof(*tmp_DD)); - MYCOMPLEX (*tmp_DS)[3] = (MYCOMPLEX(*)[3])calloc(nr, sizeof(*tmp_DS)); - MYCOMPLEX (*tmp_SS)[3] = (MYCOMPLEX(*)[3])calloc(nr, sizeof(*tmp_SS)); + MYCOMPLEX coef = -dk*fac*omega2_inv; // 最终要乘上的系数 // 局部变量，用于求和 sum F(ki,w)Jm(ki*r)ki // 维度3代表阶数m=0,1,2，维度4代表4种类型的F(k,w)Jm(kr)k的类型，详见int_Pk()函数内的注释 @@ -274,6 +574,17 @@ void integ_grn_spec( MYCOMPLEX (*sum_VF_J)[3][4] = (VFgrn != NULL) ? (MYCOMPLEX(*)[3][4])calloc(nr, sizeof(*sum_VF_J)) : NULL; MYCOMPLEX (*sum_HF_J)[3][4] = (HFgrn != NULL) ? (MYCOMPLEX(*)[3][4])calloc(nr, sizeof(*sum_HF_J)) : NULL; MYCOMPLEX (*sum_DC_J)[3][4] = (DDgrn != NULL || DSgrn != NULL || SSgrn != NULL) ? (MYCOMPLEX(*)[3][4])calloc(nr, sizeof(*sum_DC_J)) : NULL; + + MYCOMPLEX (*sum_EXP_uiz_J)[3][4] = (EXPgrn_uiz != NULL) ? (MYCOMPLEX(*)[3][4])calloc(nr, sizeof(*sum_EXP_uiz_J)) : NULL; + MYCOMPLEX (*sum_VF_uiz_J)[3][4] = (VFgrn_uiz != NULL) ? (MYCOMPLEX(*)[3][4])calloc(nr, sizeof(*sum_VF_uiz_J)) : NULL; + MYCOMPLEX (*sum_HF_uiz_J)[3][4] = (HFgrn_uiz != NULL) ? (MYCOMPLEX(*)[3][4])calloc(nr, sizeof(*sum_HF_uiz_J)) : NULL; + MYCOMPLEX (*sum_DC_uiz_J)[3][4] = (DDgrn_uiz != NULL || DSgrn_uiz != NULL || SSgrn_uiz != NULL) ? (MYCOMPLEX(*)[3][4])calloc(nr, sizeof(*sum_DC_uiz_J)) : NULL; + + MYCOMPLEX (*sum_EXP_uir_J)[3][4] = (EXPgrn_uir != NULL) ? (MYCOMPLEX(*)[3][4])calloc(nr, sizeof(*sum_EXP_uir_J)) : NULL; + MYCOMPLEX (*sum_VF_uir_J)[3][4] = (VFgrn_uir != NULL) ? (MYCOMPLEX(*)[3][4])calloc(nr, sizeof(*sum_VF_uir_J)) : NULL; + MYCOMPLEX (*sum_HF_uir_J)[3][4] = (HFgrn_uir != NULL) ? (MYCOMPLEX(*)[3][4])calloc(nr, sizeof(*sum_HF_uir_J)) : NULL; + MYCOMPLEX (*sum_DC_uir_J)[3][4] = (DDgrn_uir != NULL || DSgrn_uir != NULL || SSgrn_uir != NULL) ? (MYCOMPLEX(*)[3][4])calloc(nr, sizeof(*sum_DC_uir_J)) : NULL; + MODEL1D *local_mod1d = NULL; @@ -293,27 +604,25 @@ void integ_grn_spec( FILE **ptam_fstats = (FILE **)malloc(nr * sizeof(FILE *)); for(MYINT ir=0; ir= 0) || (findElement_MYINT(statsidxs, nstatsidxs, -1) >= 0))){ @@ -350,20 +659,32 @@ void integ_grn_spec( // 常规的波数积分 k = discrete_integ( local_mod1d, dk, kmax, keps, omega, nr, rs, - sum_EXP_J, sum_VF_J, sum_HF_J, sum_DC_J, fstats); + sum_EXP_J, sum_VF_J, sum_HF_J, sum_DC_J, + calc_upar, + sum_EXP_uiz_J, sum_VF_uiz_J, sum_HF_uiz_J, sum_DC_uiz_J, + sum_EXP_uir_J, sum_VF_uir_J, sum_HF_uir_J, sum_DC_uir_J, + fstats); } else { // 基于线性插值的Filon积分 k = linear_filon_integ( local_mod1d, dk, kmax, keps, omega, nr, rs, - sum_EXP_J, sum_VF_J, sum_HF_J, sum_DC_J, fstats); + sum_EXP_J, sum_VF_J, sum_HF_J, sum_DC_J, + calc_upar, + sum_EXP_uiz_J, sum_VF_uiz_J, sum_HF_uiz_J, sum_DC_uiz_J, + sum_EXP_uir_J, sum_VF_uir_J, sum_HF_uir_J, sum_DC_uir_J, + fstats); } // k之后的部分使用峰谷平均法进行显式收敛，建议在浅源地震的时候使用 if(vmin_ref < RZERO){ PTA_method( local_mod1d, k, dk, rmin, rmax, omega, nr, rs, - sum_EXP_J, sum_VF_J, sum_HF_J, sum_DC_J, fstats, ptam_fstats); + sum_EXP_J, sum_VF_J, sum_HF_J, sum_DC_J, + calc_upar, + sum_EXP_uiz_J, sum_VF_uiz_J, sum_HF_uiz_J, sum_DC_uiz_J, + sum_EXP_uir_J, sum_VF_uir_J, sum_HF_uir_J, sum_DC_uir_J, + fstats, ptam_fstats); } // printf("iw=%d, w=%.5e, k=%.5e, dk=%.5e, nk=%d\n", iw, w, k, dk, (int)(k/dk)); @@ -371,53 +692,21 @@ void integ_grn_spec( // 记录到格林函数结构体内 - for(MYINT ir=0; irRe[iw] = CREAL(mtmp); - EXPgrn[ir][ii]->Im[iw] = CIMAG(mtmp); - } - if(VFgrn!=NULL){ - mtmp = mtmp0*tmp_VF[ir][ii]; // m=0 垂直力源 - VFgrn[ir][ii]->Re[iw] = CREAL(mtmp); - VFgrn[ir][ii]->Im[iw] = CIMAG(mtmp); - } - if(DDgrn!=NULL){ - mtmp = mtmp0*tmp_DD[ir][ii]; // m=0 45-倾滑 - DDgrn[ir][ii]->Re[iw] = CREAL(mtmp); - DDgrn[ir][ii]->Im[iw] = CIMAG(mtmp); - } - } - if(HFgrn!=NULL){ - mtmp = mtmp0*tmp_HF[ir][ii]; // m=1 水平力源 - HFgrn[ir][ii]->Re[iw] = CREAL(mtmp); - HFgrn[ir][ii]->Im[iw] = CIMAG(mtmp); - } - if(DSgrn!=NULL){ - mtmp = mtmp0*tmp_DS[ir][ii]; // m=1 90-倾滑 - DSgrn[ir][ii]->Re[iw] = CREAL(mtmp); - DSgrn[ir][ii]->Im[iw] = CIMAG(mtmp); - } - if(SSgrn!=NULL){ - mtmp = mtmp0*tmp_SS[ir][ii]; // m=2 走滑 - SSgrn[ir][ii]->Re[iw] = CREAL(mtmp); - SSgrn[ir][ii]->Im[iw] = CIMAG(mtmp); - } - - } + recordin_GRN( + iw, nr, coef, + sum_EXP_J, sum_VF_J, sum_HF_J, sum_DC_J, + EXPgrn, VFgrn, HFgrn, DDgrn, DSgrn, SSgrn); + if(calc_upar){ + recordin_GRN( + iw, nr, coef, + sum_EXP_uiz_J, sum_VF_uiz_J, sum_HF_uiz_J, sum_DC_uiz_J, + EXPgrn_uiz, VFgrn_uiz, HFgrn_uiz, DDgrn_uiz, DSgrn_uiz, SSgrn_uiz); + recordin_GRN( + iw, nr, coef, + sum_EXP_uir_J, sum_VF_uir_J, sum_HF_uir_J, sum_DC_uir_J, + EXPgrn_uir, VFgrn_uir, HFgrn_uir, DDgrn_uir, DSgrn_uir, SSgrn_uir); } + for(MYINT ir=0; ir,[,...] [-O] [-H/] \n" " [-L] [-V] [-E[/]] \n" " [-K[///]] [-P]\n" -" [-G[///]] [-S,[,...]] [-s]\n" +" [-G[///]] [-S,[,...]] [-e]\n" +" [-s]\n" "\n\n" "Options:\n" "----------------------------------------------------------------\n" @@ -274,7 +278,16 @@ printf("\n" " integration to be output, require 0 <= i <= nf-1,\n" " where nf=nt/2+1. These option is designed to check\n" " the trend of kernel and integrand with wavenumber.\n" -" \n" +"\n" +" -e Compute the spatial derivatives, ui_z and ui_r,\n" +" of displacement u. In filenames, prefix \"r\" means \n" +" ui_r and \"z\" means ui_z. The units of derivatives\n" +" for different sources are: \n" +" + Explosion: 1e-25 /(dyne-cm)\n" +" + Single Force: 1e-20 /(dyne)\n" +" + Double Couple: 1e-25 /(dyne-cm)\n" +" + Moment Tensor: 1e-25 /(dyne-cm)\n" +"\n" " -s Silence all outputs.\n" "\n" " -h Display this help message.\n" @@ -321,7 +334,7 @@ static char* get_basename(char* path) { */ static void getopt_from_command(int argc, char **argv){ int opt; - while ((opt = getopt(argc, argv, ":M:D:N:O:H:L:V:E:K:shR:S:P:G:")) != -1) { + while ((opt = getopt(argc, argv, ":M:D:N:O:H:L:V:E:K:shR:S:P:G:e")) != -1) { switch (opt) { // 模型路径，其中每行分别为 // 厚度(km) Vp(km/s) Vs(km/s) Rho(g/cm^3) Qp Qs @@ -550,6 +563,12 @@ static void getopt_from_command(int argc, char **argv){ free(str_copy); } break; + + // 是否计算位移空间导数 + case 'e': + e_flag = 1; + calc_upar = true; + break; // 帮助 case 'h': @@ -668,8 +687,8 @@ static void print_long_array_in_tabel( * * @param delay 时间延迟 * @param mult 幅值放大系数 - * @param grn 某一道格林函数频谱，GRN结构体指针 - * @param fftw_grn 将GRN结构体的频谱写到FFTW_COMPLEX类型中 + * @param grncplx 复数形式的格林函数频谱 + * @param fftw_grn 将频谱写到FFTW_COMPLEX类型中 * @param out ifft后的时域数据 * @param float_arr 将时域数据写到float类型的数组中 * @param plan FFTW_PLAN @@ -679,7 +698,7 @@ static void print_long_array_in_tabel( static void ifft_one_trace( MYREAL delay, MYREAL mult, MYCOMPLEX *grncplx, _FFTW_COMPLEX *fftw_grn, MYREAL *out, float *float_arr, - _FFTW_PLAN plan, SACHEAD hd, const char *outpath) + _FFTW_PLAN plan, SACHEAD *hd, const char *outpath) { // 赋值复数，包括时移 MYCOMPLEX cfac, ccoef; @@ -708,9 +727,9 @@ static void ifft_one_trace( } // 写入虚频率 - hd.user0 = wI; + hd->user0 = wI; - write_sac(outpath, hd, float_arr); + write_sac(outpath, *hd, float_arr); // FILE *fp = fopen(outpath, "wb"); // fwrite(out, sizeof(float), nt, fp); // fclose(fp); @@ -910,6 +929,36 @@ static void print_outdir_travt(const char *s_output_subdir, const char *s_R, dou } +/** + * 将一条数据反变换回时间域再进行处理，保存到SAC文件 + * + * @param srcname 震源类型 + * @param ch 三分量类型（Z,R,T） + * @param hd SAC头段变量结构体指针 + * @param s_outpath 用于接收保存路径字符串 + * @param s_output_subdir 保存路径所在文件夹 + * @param s_prefix sac文件名以及通道名名称前缀 + * @param sgn 数据待乘符号(-1/1) + * @param grncplx 复数形式的格林函数频谱 + * @param fftw_grn 将频谱写到FFTW_COMPLEX类型中 + * @param out ifft后的时域数据 + * @param float_arr 将时域数据写到float类型的数组中 + * @param plan FFTW_PLAN + * + */ +static void write_one_to_sac( + const char *srcname, const char ch, + SACHEAD *hd, char *s_outpath, const char *s_output_subdir, const char *s_prefix, + const int sgn, MYCOMPLEX *grncplx, fftw_complex *fftw_grn, MYREAL *out, float *float_arr, fftw_plan plan) +{ + char kcmpnm[9]; + snprintf(kcmpnm, sizeof(kcmpnm), "%s%s%c", s_prefix, srcname, ch); + strcpy(hd->kcmpnm, kcmpnm); + sprintf(s_outpath, "%s/%s.sac", s_output_subdir, kcmpnm); + ifft_one_trace(delayT, sgn, grncplx, fftw_grn, out, float_arr, plan, hd, s_outpath); +} + + //==================================================================================== //==================================================================================== @@ -1010,16 +1059,48 @@ int main(int argc, char **argv) { MYCOMPLEX *(*DScplx)[3] = (doDC) ? (MYCOMPLEX*(*)[3])calloc(nr, sizeof(*DScplx)) : NULL; MYCOMPLEX *(*SScplx)[3] = (doDC) ? (MYCOMPLEX*(*)[3])calloc(nr, sizeof(*SScplx)) : NULL; + MYCOMPLEX *(*EXPcplx_uiz)[2] = (calc_upar && doEXP) ? (MYCOMPLEX*(*)[2])calloc(nr, sizeof(*EXPcplx_uiz)) : NULL; + MYCOMPLEX *(*VFcplx_uiz)[2] = (calc_upar && doVF) ? (MYCOMPLEX*(*)[2])calloc(nr, sizeof(*VFcplx_uiz)) : NULL; + MYCOMPLEX *(*HFcplx_uiz)[3] = (calc_upar && doHF) ? (MYCOMPLEX*(*)[3])calloc(nr, sizeof(*HFcplx_uiz)) : NULL; + MYCOMPLEX *(*DDcplx_uiz)[2] = (calc_upar && doDC) ? (MYCOMPLEX*(*)[2])calloc(nr, sizeof(*DDcplx_uiz)) : NULL; + MYCOMPLEX *(*DScplx_uiz)[3] = (calc_upar && doDC) ? (MYCOMPLEX*(*)[3])calloc(nr, sizeof(*DScplx_uiz)) : NULL; + MYCOMPLEX *(*SScplx_uiz)[3] = (calc_upar && doDC) ? (MYCOMPLEX*(*)[3])calloc(nr, sizeof(*SScplx_uiz)) : NULL; + + MYCOMPLEX *(*EXPcplx_uir)[2] = (calc_upar && doEXP) ? (MYCOMPLEX*(*)[2])calloc(nr, sizeof(*EXPcplx_uir)) : NULL; + MYCOMPLEX *(*VFcplx_uir)[2] = (calc_upar && doVF) ? (MYCOMPLEX*(*)[2])calloc(nr, sizeof(*VFcplx_uir)) : NULL; + MYCOMPLEX *(*HFcplx_uir)[3] = (calc_upar && doHF) ? (MYCOMPLEX*(*)[3])calloc(nr, sizeof(*HFcplx_uir)) : NULL; + MYCOMPLEX *(*DDcplx_uir)[2] = (calc_upar && doDC) ? (MYCOMPLEX*(*)[2])calloc(nr, sizeof(*DDcplx_uir)) : NULL; + MYCOMPLEX *(*DScplx_uir)[3] = (calc_upar && doDC) ? (MYCOMPLEX*(*)[3])calloc(nr, sizeof(*DScplx_uir)) : NULL; + MYCOMPLEX *(*SScplx_uir)[3] = (calc_upar && doDC) ? (MYCOMPLEX*(*)[3])calloc(nr, sizeof(*SScplx_uir)) : NULL; + for(int ir=0; ir/////]\n" " [-F//] [-O] \n" " [-D/] [-I] [-J]\n" -" [-P] [-s]\n" +" [-P] [-e] [-s]\n" "\n" "\n\n" "Options:\n" @@ -189,14 +201,17 @@ printf("\n" "\n" " -J Order of differentiation. Default not use\n" "\n" +" -e Compute the spatial derivatives, ui_z and ui_r,\n" +" of displacement u. In filenames, prefix \"r\" means \n" +" ui_r and \"z\" means ui_z. \n" +"\n" " -s Silence all outputs.\n" "\n" " -h Display this help message.\n" "\n\n" "Examples:\n" "----------------------------------------------------------------\n" -" Say you have computed computed Green's functions with following \n" -" command:\n" +" Say you have computed Green's functions with following command:\n" " grt -Mmilrow -N1000/0.01 -D2/0 -Ores -R2,4,6,8,10\n" "\n" " Then you can get synthetic seismograms of Explosion at epicentral\n" @@ -240,7 +255,7 @@ static void check_grn_exist(const char *name){ */ static void getopt_from_command(int argc, char **argv){ int opt; - while ((opt = getopt(argc, argv, ":G:A:S:M:F:T:O:P:D:I:J:hs")) != -1) { + while ((opt = getopt(argc, argv, ":G:A:S:M:F:T:O:P:D:I:J:ehs")) != -1) { switch (opt) { // 格林函数路径 case 'G': @@ -267,6 +282,8 @@ static void getopt_from_command(int argc, char **argv){ fprintf(stderr, "[%s] " BOLD_RED "Error! Azimuth in -A must be in [0, 360].\n" DEFAULT_RESTORE, command); exit(EXIT_FAILURE); } + backazimuth = 180.0 + azimuth; + if(backazimuth >= 360.0) backazimuth -= 360.0; azrad = azimuth * DEG1; saz = sin(azrad); caz = cos(azrad); @@ -287,6 +304,7 @@ static void getopt_from_command(int argc, char **argv){ case 'M': M_flag = 1; computeType = COMPUTE_DC; + sprintf(s_computeType, "%s", "DC"); if(3 != sscanf(optarg, "%lf/%lf/%lf", &strike, &dip, &rake)){ fprintf(stderr, "[%s] " BOLD_RED "Error in -M.\n" DEFAULT_RESTORE, command); exit(EXIT_FAILURE); @@ -309,6 +327,7 @@ static void getopt_from_command(int argc, char **argv){ case 'F': F_flag = 1; computeType = COMPUTE_SF; + sprintf(s_computeType, "%s", "SF"); if(3 != sscanf(optarg, "%lf/%lf/%lf", &fn, &fe, &fz)){ fprintf(stderr, "[%s] " BOLD_RED "Error in -F.\n" DEFAULT_RESTORE, command); exit(EXIT_FAILURE); @@ -319,6 +338,7 @@ static void getopt_from_command(int argc, char **argv){ case 'T': T_flag = 1; computeType = COMPUTE_MT; + sprintf(s_computeType, "%s", "MT"); if(6 != sscanf(optarg, "%lf/%lf/%lf/%lf/%lf/%lf", &Mxx, &Mxy, &Mxz, &Myy, &Myz, &Mzz)){ fprintf(stderr, "[%s] " BOLD_RED "Error in -T.\n" DEFAULT_RESTORE, command); exit(EXIT_FAILURE); @@ -380,6 +400,13 @@ static void getopt_from_command(int argc, char **argv){ } break; + // 是否计算位移空间导数 + case 'e': + e_flag = 1; + calc_upar = true; + calcUTypes = 4; + break; + // 不打印在终端 case 's': s_flag = 1; @@ -436,6 +463,10 @@ static void getopt_from_command(int argc, char **argv){ if( (M_flag==0&&F_flag==0&&T_flag==0) || T_flag == 1){ check_grn_exist("EXR.sac"); check_grn_exist("EXZ.sac"); + if(calc_upar) { + check_grn_exist("zEXR.sac"); check_grn_exist("rEXR.sac"); + check_grn_exist("zEXZ.sac"); check_grn_exist("rEXZ.sac"); + } } if(M_flag == 1){ check_grn_exist("DDR.sac"); @@ -446,6 +477,16 @@ static void getopt_from_command(int argc, char **argv){ check_grn_exist("SSR.sac"); check_grn_exist("SST.sac"); check_grn_exist("SSZ.sac"); + if(calc_upar){ + check_grn_exist("zDDR.sac"); check_grn_exist("rDDR.sac"); + check_grn_exist("zDDZ.sac"); check_grn_exist("rDDZ.sac"); + check_grn_exist("zDSR.sac"); check_grn_exist("rDSR.sac"); + check_grn_exist("zDST.sac"); check_grn_exist("rDST.sac"); + check_grn_exist("zDSZ.sac"); check_grn_exist("rDSZ.sac"); + check_grn_exist("zSSR.sac"); check_grn_exist("rSSR.sac"); + check_grn_exist("zSST.sac"); check_grn_exist("rSST.sac"); + check_grn_exist("zSSZ.sac"); check_grn_exist("rSSZ.sac"); + } } if(F_flag == 1){ check_grn_exist("VFR.sac"); @@ -453,8 +494,15 @@ static void getopt_from_command(int argc, char **argv){ check_grn_exist("HFR.sac"); check_grn_exist("HFT.sac"); check_grn_exist("HFZ.sac"); + if(calc_upar){ + check_grn_exist("zVFR.sac"); check_grn_exist("rVFR.sac"); + check_grn_exist("zVFZ.sac"); check_grn_exist("rVFZ.sac"); + check_grn_exist("zHFR.sac"); check_grn_exist("rHFR.sac"); + check_grn_exist("zHFT.sac"); check_grn_exist("rHFT.sac"); + check_grn_exist("zHFZ.sac"); check_grn_exist("rHFZ.sac"); + } } - + if(O_flag == 1){ // 建立保存目录 @@ -481,14 +529,16 @@ static void getopt_from_command(int argc, char **argv){ * 将某一道合成地震图保存到sac文件 * * @param buffer 输出文件夹字符串(重复使用) + * @param s_prefix2 SAC文件名和通道名前缀 * @param ch 分量名， Z/R/T * @param arr 数据指针 * @param hd SAC头段变量 */ -static void save_to_sac(char *buffer, const char ch, float *arr, SACHEAD hd){ +static void save_to_sac(char *buffer, const char *s_prefix2, const char ch, float *arr, SACHEAD hd){ hd.az = azimuth; - sprintf(hd.kcmpnm, "HH%c", ch); - sprintf(buffer, "%s/%s%c.sac", s_output_dir, s_prefix, ch); + hd.baz = backazimuth; + snprintf(hd.kcmpnm, sizeof(hd.kcmpnm), "%s%s%c", s_prefix2, s_computeType, ch); + sprintf(buffer, "%s/%s%s%c.sac", s_output_dir, s_prefix2, s_prefix, ch); write_sac(buffer, hd, arr); } @@ -508,25 +558,33 @@ static void save_tf_to_sac(char *buffer, float *tfarr, int tfnt, float dt){ /** * 设置每个震源的方向因子 + * + * @param par_theta 方向因子中是否对theta(az)求导 + * @param coef 缩放系数，用于位移空间导数的计算 */ -static void set_source_coef(){ +static void set_source_coef(const bool par_theta, const double coef){ double mult; if(computeType == COMPUTE_SF){ - mult = 1e-15*M0; + mult = 1e-15*M0*coef; } else { - mult = 1e-20*M0; + mult = 1e-20*M0*coef; } if(computeType == COMPUTE_EXP){ - srcCoef[0][0] = srcCoef[1][0] = mult; // Z/R + srcCoef[0][0] = srcCoef[1][0] = (par_theta)? 0.0 : mult; // Z/R srcCoef[2][0] = 0.0; // T } else if(computeType == COMPUTE_SF){ + double A0, A1, A4; + A0 = fz*mult; + A1 = (fn*caz + fe*saz)*mult; + A4 = (- fn*saz + fe*caz)*mult; + // 公式(4.6.20) - srcCoef[0][1] = srcCoef[1][1] = fz*mult; // VF, Z/R - srcCoef[0][2] = srcCoef[1][2] = (fn*caz + fe*saz)*mult; // HF, Z/R + srcCoef[0][1] = srcCoef[1][1] = (par_theta)? 0.0 : A0; // VF, Z/R + srcCoef[0][2] = srcCoef[1][2] = (par_theta)? A4 : A1; // HF, Z/R srcCoef[2][1] = 0.0; // VF, T - srcCoef[2][2] = (- fn*saz + fe*caz)*mult; // HF, T + srcCoef[2][2] = (par_theta)? -A1 : A4; // HF, T } else if(computeType == COMPUTE_DC){ // 公式(4.8.35) @@ -548,12 +606,12 @@ static void set_source_coef(){ A4 = mult * (- cdip2*srak*cthe - cdip*crak*sthe); A5 = mult * (sdip*crak*cthe2 - 0.5*sdip2*srak*sthe2); - srcCoef[0][3] = srcCoef[1][3] = A0; // DD, Z/R - srcCoef[0][4] = srcCoef[1][4] = A1; // DS, Z/R - srcCoef[0][5] = srcCoef[1][5] = A2; // SS, Z/R + srcCoef[0][3] = srcCoef[1][3] = (par_theta)? 0.0 : A0; // DD, Z/R + srcCoef[0][4] = srcCoef[1][4] = (par_theta)? A4 : A1; // DS, Z/R + srcCoef[0][5] = srcCoef[1][5] = (par_theta)? 2.0*A5 : A2; // SS, Z/R srcCoef[2][3] = 0.0; // DD, T - srcCoef[2][4] = A4; // DS, T - srcCoef[2][5] = A5; // DS, T + srcCoef[2][4] = (par_theta)? -A1 : A4; // DS, T + srcCoef[2][5] = (par_theta)? -2.0*A2 : A5; // DS, T } else if(computeType == COMPUTE_MT){ // 公式(4.9.7)但修改了各向同性的量 @@ -573,24 +631,26 @@ static void set_source_coef(){ A4 = mult * (M13*saz - M23*caz); A5 = mult * (-0.5*(M11 - M22)*saz2 + M12*caz2); - srcCoef[0][0] = srcCoef[1][0] = mult*Mexp; // EX, Z/R - srcCoef[0][3] = srcCoef[1][3] = A0; // DD, Z/R - srcCoef[0][4] = srcCoef[1][4] = A1; // DS, Z/R - srcCoef[0][5] = srcCoef[1][5] = A2; // SS, Z/R + srcCoef[0][0] = srcCoef[1][0] = (par_theta)? 0.0 : mult*Mexp; // EX, Z/R + srcCoef[0][3] = srcCoef[1][3] = (par_theta)? 0.0 : A0; // DD, Z/R + srcCoef[0][4] = srcCoef[1][4] = (par_theta)? A4 : A1; // DS, Z/R + srcCoef[0][5] = srcCoef[1][5] = (par_theta)? 2.0*A5 : A2; // SS, Z/R srcCoef[2][0] = 0.0; // EX, T srcCoef[2][3] = 0.0; // DD, T - srcCoef[2][4] = A4; // DS, T - srcCoef[2][5] = A5; // DS, T + srcCoef[2][4] = (par_theta)? -A1 : A4; // DS, T + srcCoef[2][5] = (par_theta)? -2.0*A2 : A5; // DS, T } } +//==================================================================================== +//==================================================================================== +//==================================================================================== int main(int argc, char **argv){ command = argv[0]; getopt_from_command(argc, argv); - set_source_coef(); char *buffer = (char*)malloc(sizeof(char)*(strlen(s_grnpath)+strlen(s_output_dir)+strlen(s_prefix)+100)); float *arr, *arrout=NULL, *convarr=NULL; @@ -602,85 +662,121 @@ int main(int argc, char **argv){ float wI=0.0; // 虚频率 int nt=0; float dt=0.0; + float dist=-12345.0; // 震中距 + double upar_scale=1.0; - for(int c=0; c<3; ++c){ - ch = chs[c]; - for(int k=0; kircvup; MYINT isrc = mod1d->isrc; // 震源所在虚拟层位, isrc>=1 MYINT ircv = mod1d->ircv; // 接收点所在虚拟层位, ircv>=1, ircv != isrc @@ -119,9 +133,13 @@ void kernel( // 自由表面的反射系数 MYCOMPLEX R_tilt[2][2] = INIT_C_ZERO_2x2_MATRIX; // SH波在自由表面的反射系数为1，不必定义变量 - // 接收点处的接收矩阵 + // 接收点处的接收矩阵(转为位移u的(B_m, C_m, P_m)系分量) MYCOMPLEX R_EV[2][2], R_EVL; MYCOMPLEX *const pR_EVL = &R_EVL; + + // 接收点处的接收矩阵(转为位移导数ui_z的(B_m, C_m, P_m)系分量) + MYCOMPLEX uiz_R_EV[2][2], uiz_R_EVL; + MYCOMPLEX *const puiz_R_EVL = &uiz_R_EVL; // 模型参数 @@ -331,7 +349,7 @@ void kernel( source_coef(src_xa, src_xb, src_kaka, src_kbkb, omega, k, pEXP, pVF, pHF, pDC); // 临时中转矩阵 (temperary) - MYCOMPLEX tmpR1[2][2], tmpR2[2][2], tmp2x2[2][2], tmpRL; + MYCOMPLEX tmpR1[2][2], tmpR2[2][2], tmp2x2[2][2], tmpRL, tmp2x2_uiz[2][2], tmpRL_uiz; MYCOMPLEX inv_2x2T[2][2], invT; // 递推RU_FA @@ -360,12 +378,15 @@ void kernel( RU_FB, pRUL_FB, inv_2x2T, &invT); // 公式(5.7.12-14) - cmat2x2_mul(R_EV, inv_2x2T, tmpR1); + // cmat2x2_mul(R_EV, inv_2x2T, tmpR1); cmat2x2_mul(RD_BL, RU_FB, tmpR2); cmat2x2_one_sub(tmpR2); cmat2x2_inv(tmpR2, tmpR2);// (I - xx)^-1 - cmat2x2_mul(tmpR1, tmpR2, tmp2x2); - + cmat2x2_mul(inv_2x2T, tmpR2, tmp2x2); + + if(calc_uiz) cmat2x2_assign(tmp2x2, tmp2x2_uiz); // 为后续计算空间导数备份 + + cmat2x2_mul(R_EV, tmp2x2, tmp2x2); tmpRL = R_EVL * invT / (RONE - RDL_BL * RUL_FB); for(MYINT m=0; m<3; ++m){ if(0==m){ @@ -381,6 +402,27 @@ void kernel( // 剪切位错 if(DC_qwv!=NULL) get_qwv(ircvup, tmp2x2, tmpRL, RD_BL, RDL_BL, DC[m], DC_qwv[m]); } + + + if(calc_uiz){ + calc_uiz_R_EV(rcv_xa, rcv_xb, ircvup, k, RU_FA, RUL_FA, uiz_R_EV, puiz_R_EVL); + cmat2x2_mul(uiz_R_EV, tmp2x2_uiz, tmp2x2_uiz); + tmpRL_uiz = tmpRL / R_EVL * uiz_R_EVL; + for(MYINT m=0; m<3; ++m){ + if(0==m){ + // 爆炸源 + if(EXP_uiz_qwv!=NULL) get_qwv(ircvup, tmp2x2_uiz, tmpRL_uiz, RD_BL, RDL_BL, EXP[m], EXP_uiz_qwv[m]); + // 垂直力源 + if(VF_uiz_qwv!=NULL) get_qwv(ircvup, tmp2x2_uiz, tmpRL_uiz, RD_BL, RDL_BL, VF[m], VF_uiz_qwv[m]); + } + + // 水平力源 + if(1==m && HF_uiz_qwv!=NULL) get_qwv(ircvup, tmp2x2_uiz, tmpRL_uiz, RD_BL, RDL_BL, HF[m], HF_uiz_qwv[m]); + + // 剪切位错 + if(DC_uiz_qwv!=NULL) get_qwv(ircvup, tmp2x2_uiz, tmpRL_uiz, RD_BL, RDL_BL, DC[m], DC_uiz_qwv[m]); + } + } } else { // A震源 B接收 @@ -397,11 +439,15 @@ void kernel( RD_AL, pRDL_AL, inv_2x2T, &invT); // 公式(5.7.26-27) - cmat2x2_mul(R_EV, inv_2x2T, tmpR1); + // cmat2x2_mul(R_EV, inv_2x2T, tmpR1); cmat2x2_mul(RU_FA, RD_AL, tmpR2); cmat2x2_one_sub(tmpR2); cmat2x2_inv(tmpR2, tmpR2);// (I - xx)^-1 - cmat2x2_mul(tmpR1, tmpR2, tmp2x2); + cmat2x2_mul(inv_2x2T, tmpR2, tmp2x2); + + if(calc_uiz) cmat2x2_assign(tmp2x2, tmp2x2_uiz); // 为后续计算空间导数备份 + + cmat2x2_mul(R_EV, tmp2x2, tmp2x2); tmpRL = R_EVL * invT / (RONE - RUL_FA * RDL_AL); for(MYINT m=0; m<3; ++m){ if(0==m){ @@ -418,6 +464,30 @@ void kernel( if(DC_qwv!=NULL) get_qwv(ircvup, tmp2x2, tmpRL, RU_FA, RUL_FA, DC[m], DC_qwv[m]); } + + + if(calc_uiz){ + calc_uiz_R_EV(rcv_xa, rcv_xb, ircvup, k, RD_BL, RDL_BL, uiz_R_EV, puiz_R_EVL); + cmat2x2_mul(uiz_R_EV, tmp2x2_uiz, tmp2x2_uiz); + tmpRL_uiz = tmpRL / R_EVL * uiz_R_EVL; + for(MYINT m=0; m<3; ++m){ + if(0==m){ + // 爆炸源 + if(EXP_uiz_qwv!=NULL) get_qwv(ircvup, tmp2x2_uiz, tmpRL_uiz, RU_FA, RUL_FA, EXP[m], EXP_uiz_qwv[m]); + // 垂直力源 + if(VF_uiz_qwv!=NULL) get_qwv(ircvup, tmp2x2_uiz, tmpRL_uiz, RU_FA, RUL_FA, VF[m], VF_uiz_qwv[m]); + } + + // 水平力源 + if(1==m && HF_uiz_qwv!=NULL) get_qwv(ircvup, tmp2x2_uiz, tmpRL_uiz, RU_FA, RUL_FA, HF[m], HF_uiz_qwv[m]); + + // 剪切位错 + if(DC_uiz_qwv!=NULL) get_qwv(ircvup, tmp2x2_uiz, tmpRL_uiz, RU_FA, RUL_FA, DC[m], DC_uiz_qwv[m]); + + } + } + + } // END if @@ -431,17 +501,38 @@ void int_Pk( const MYCOMPLEX EXP_qwv[3][3], const MYCOMPLEX VF_qwv[3][3], const MYCOMPLEX HF_qwv[3][3], const MYCOMPLEX DC_qwv[3][3], // F(ki,w)Jm(ki*r)ki，维度3代表阶数m=0,1,2，维度4代表4种类型的F(k,w)Jm(kr)k的类型 + bool calc_uir, MYCOMPLEX EXP_J[3][4], MYCOMPLEX VF_J[3][4], - MYCOMPLEX HF_J[3][4], MYCOMPLEX DC_J[3][4] ) + MYCOMPLEX HF_J[3][4], MYCOMPLEX DC_J[3][4]) { MYREAL bj0k, bj1k, bj2k; MYREAL kr = k*r; + MYREAL kr_inv = RONE/kr; + MYREAL kcoef = k; + + MYREAL J1coef, J2coef; bessel012(kr, &bj0k, &bj1k, &bj2k); + if(calc_uir){ + MYREAL j1, j2; + j1 = bj1k; + j2 = bj2k; + besselp012(kr, &bj0k, &bj1k, &bj2k); + kcoef = k*k; + + J1coef = kr_inv * (-kr_inv * j1 + bj1k); + J2coef = kr_inv * (-kr_inv * j2 + bj2k); + } else { + J1coef = bj1k*kr_inv; + J2coef = bj2k*kr_inv; + } + + J1coef *= kcoef; + J2coef *= kcoef; - bj0k = bj0k*k; - bj1k = bj1k*k; - bj2k = bj2k*k; + bj0k *= kcoef; + bj1k *= kcoef; + bj2k *= kcoef; if(EXP_qwv!=NULL){ @@ -460,7 +551,7 @@ void int_Pk( if(HF_qwv!=NULL){ // m=1 水平力源 HF_J[1][0] = HF_qwv[1][0]*bj0k; // q1*J0*k - HF_J[1][1] = - (HF_qwv[1][0] + HF_qwv[1][2])*bj1k/(kr); // - (q1+v1)*J1*k/kr + HF_J[1][1] = - (HF_qwv[1][0] + HF_qwv[1][2])*J1coef; // - (q1+v1)*J1*k/kr HF_J[1][2] = HF_qwv[1][1]*bj1k; // w1*J1*k HF_J[1][3] = - HF_qwv[1][2]*bj0k; // -v1*J0*k } @@ -472,13 +563,13 @@ void int_Pk( // m=1 双力偶源 DC_J[1][0] = DC_qwv[1][0]*bj0k; // q1*J0*k - DC_J[1][1] = - (DC_qwv[1][0] + DC_qwv[1][2])*bj1k/(kr); // - (q1+v1)*J1*k/kr + DC_J[1][1] = - (DC_qwv[1][0] + DC_qwv[1][2])*J1coef; // - (q1+v1)*J1*k/kr DC_J[1][2] = DC_qwv[1][1]*bj1k; // w1*J1*k DC_J[1][3] = - DC_qwv[1][2]*bj0k; // -v1*J0*k // m=2 双力偶源 DC_J[2][0] = DC_qwv[2][0]*bj1k; // q2*J1*k - DC_J[2][1] = - RTWO*(DC_qwv[2][0] + DC_qwv[2][2])*bj2k/(kr); // - (q2+v2)*J2*k/kr + DC_J[2][1] = - RTWO*(DC_qwv[2][0] + DC_qwv[2][2])*J2coef; // - (q2+v2)*J2*k/kr DC_J[2][2] = DC_qwv[2][1]*bj2k; // w2*J2*k DC_J[2][3] = - DC_qwv[2][2]*bj1k; // -v2*J1*k } diff --git a/pygrt/C_extension/src/dynamic/ptam.c b/pygrt/C_extension/src/dynamic/ptam.c index a65f7587..8e8e8436 100755 --- a/pygrt/C_extension/src/dynamic/ptam.c +++ b/pygrt/C_extension/src/dynamic/ptam.c @@ -63,11 +63,131 @@ static void process_peak_or_trough( } +/** + * 在输入被积函数的情况下，对不同震源使用峰谷平均法 + * + * @param ir 震中距索引 + * @param nr 震中距个数 + * @param precoef 积分值系数 + * @param maxNpt 最大峰谷数 + * @param maxnwait 最大等待次数 + * @param k 波数 + * @param dk 波数步长 + * @param EXP_J3 爆炸源对应的被积函数的幅值数组，下同 + * @param VF_J3 垂直力源 + * @param HF_J3 水平力源 + * @param DC_J3 双力偶源 + * @param sum_EXP_J 爆炸源对应的积分值数组，下同 + * @param sum_VF_J 垂直力源 + * @param sum_HF_J 水平力源 + * @param sum_DC_J 双力偶源 + * + * @param kEXPpt 爆炸源对应的积分值峰谷的波数数组，下同 + * @param EXPpt 爆炸源对应的用于存储波峰/波谷点的幅值数组，下同 + * @param iEXPpt 爆炸源对应的用于存储波峰/波谷点的个数数组，下同 + * @param gEXPpt 爆炸源对应的用于存储等待迭次数的数组，下同 + * @param kVFpt 垂直力源 + * @param VFpt 垂直力源 + * @param iVFpt 垂直力源 + * @param gVFpt 垂直力源 + * @param kHFpt 水平力源 + * @param HFpt 水平力源 + * @param iHFpt 水平力源 + * @param gHFpt 水平力源 + * @param kDCpt 双力偶源 + * @param DCpt 双力偶源 + * @param iDCpt 双力偶源 + * @param gDCpt 双力偶源 + * + * + */ +static void ptam_once( + const MYINT ir, const MYINT nr, const MYREAL precoef, + MYINT maxNpt, MYINT maxnwait, MYREAL k, MYREAL dk, + MYCOMPLEX EXP_J3[nr][3][3][4], MYCOMPLEX VF_J3[nr][3][3][4], + MYCOMPLEX HF_J3[nr][3][3][4], MYCOMPLEX DC_J3[nr][3][3][4], + MYCOMPLEX sum_EXP_J[nr][3][4], MYCOMPLEX sum_VF_J[nr][3][4], + MYCOMPLEX sum_HF_J[nr][3][4], MYCOMPLEX sum_DC_J[nr][3][4], + MYREAL kEXPpt[nr][3][4][maxNpt], MYCOMPLEX EXPpt[nr][3][4][maxNpt], MYINT iEXPpt[nr][3][4], MYINT gEXPpt[nr][3][4], + MYREAL kVFpt[nr][3][4][maxNpt], MYCOMPLEX VFpt[nr][3][4][maxNpt], MYINT iVFpt[nr][3][4], MYINT gVFpt[nr][3][4], + MYREAL kHFpt[nr][3][4][maxNpt], MYCOMPLEX HFpt[nr][3][4][maxNpt], MYINT iHFpt[nr][3][4], MYINT gHFpt[nr][3][4], + MYREAL kDCpt[nr][3][4][maxNpt], MYCOMPLEX DCpt[nr][3][4][maxNpt], MYINT iDCpt[nr][3][4], MYINT gDCpt[nr][3][4], + bool *iendk0) +{ + // 赋更新量 + for(MYINT m=0; m<3; ++m){ + for(MYINT v=0; v<4; ++v){ + // EXP_J3, VF_J3, HF_J3, DC_J3转为求和结果 + if(sum_EXP_J!=NULL) { + sum_EXP_J[ir][m][v] += EXP_J3[ir][2][m][v] * precoef; + EXP_J3[ir][2][m][v] = sum_EXP_J[ir][m][v]; + } + if(sum_VF_J!=NULL){ + sum_VF_J[ir][m][v] += VF_J3[ir][2][m][v] * precoef; + VF_J3[ir][2][m][v] = sum_VF_J[ir][m][v]; + } + if(sum_HF_J!=NULL){ + sum_HF_J[ir][m][v] += HF_J3[ir][2][m][v] * precoef; + HF_J3[ir][2][m][v] = sum_HF_J[ir][m][v]; + } + if(sum_DC_J!=NULL){ + sum_DC_J[ir][m][v] += DC_J3[ir][2][m][v] * precoef; + DC_J3[ir][2][m][v] = sum_DC_J[ir][m][v]; + } + + } + } + + // 3点以上，判断波峰波谷 + *iendk0 = true; + for (MYINT m = 0; m < 3; ++m) { + for (MYINT v = 0; v < 4; ++v) { + if (sum_EXP_J != NULL && m == 0 && (v == 0 || v == 2)) { + process_peak_or_trough(ir, m, v, maxNpt, maxnwait, k, dk, EXP_J3, kEXPpt, EXPpt, iEXPpt, gEXPpt, iendk0); + } + if (sum_VF_J != NULL && m == 0 && (v == 0 || v == 2)) { + process_peak_or_trough(ir, m, v, maxNpt, maxnwait, k, dk, VF_J3, kVFpt, VFpt, iVFpt, gVFpt, iendk0); + } + if (sum_HF_J != NULL && m == 1) { + process_peak_or_trough(ir, m, v, maxNpt, maxnwait, k, dk, HF_J3, kHFpt, HFpt, iHFpt, gHFpt, iendk0); + } + if (sum_DC_J != NULL && ((m == 0 && (v == 0 || v == 2)) || m != 0)) { + process_peak_or_trough(ir, m, v, maxNpt, maxnwait, k, dk, DC_J3, kDCpt, DCpt, iDCpt, gDCpt, iendk0); + } + } + } + + + // 左移动点, + for(MYINT m=0; m<3; ++m){ + for(MYINT v=0; v<4; ++v){ + for(MYINT jj=0; jj<2; ++jj){ + if(sum_EXP_J!=NULL) EXP_J3[ir][jj][m][v] = EXP_J3[ir][jj+1][m][v]; + if(sum_VF_J!=NULL) VF_J3[ir][jj][m][v] = VF_J3[ir][jj+1][m][v]; + if(sum_HF_J!=NULL) HF_J3[ir][jj][m][v] = HF_J3[ir][jj+1][m][v]; + if(sum_DC_J!=NULL) DC_J3[ir][jj][m][v] = DC_J3[ir][jj+1][m][v]; + } + + // 点数+1 + if(sum_EXP_J!=NULL) gEXPpt[ir][m][v]++; + if(sum_VF_J!=NULL) gVFpt[ir][m][v]++; + if(sum_HF_J!=NULL) gHFpt[ir][m][v]++; + if(sum_DC_J!=NULL) gDCpt[ir][m][v]++; + } + } +} + + void PTA_method( const MODEL1D *mod1d, MYREAL k0, MYREAL predk, MYREAL rmin, MYREAL rmax, MYCOMPLEX omega, MYINT nr, MYREAL *rs, MYCOMPLEX sum_EXP_J0[nr][3][4], MYCOMPLEX sum_VF_J0[nr][3][4], MYCOMPLEX sum_HF_J0[nr][3][4], MYCOMPLEX sum_DC_J0[nr][3][4], + bool calc_upar, + MYCOMPLEX sum_EXP_uiz_J0[nr][3][4], MYCOMPLEX sum_VF_uiz_J0[nr][3][4], + MYCOMPLEX sum_HF_uiz_J0[nr][3][4], MYCOMPLEX sum_DC_uiz_J0[nr][3][4], + MYCOMPLEX sum_EXP_uir_J0[nr][3][4], MYCOMPLEX sum_VF_uir_J0[nr][3][4], + MYCOMPLEX sum_HF_uir_J0[nr][3][4], MYCOMPLEX sum_DC_uir_J0[nr][3][4], FILE *(fstats[nr]), FILE *(ptam_fstats[nr])) { // 需要兼容对正常收敛而不具有规律波峰波谷的序列 @@ -86,6 +206,13 @@ void PTA_method( MYCOMPLEX (*pHF_qwv)[3] = (sum_HF_J0!=NULL)? HF_qwv : NULL; MYCOMPLEX (*pDC_qwv)[3] = (sum_DC_J0!=NULL)? DC_qwv : NULL; + MYCOMPLEX EXP_uiz_qwv[3][3], VF_uiz_qwv[3][3], HF_uiz_qwv[3][3], DC_uiz_qwv[3][3]; + MYCOMPLEX (*pEXP_uiz_qwv)[3] = (sum_EXP_uiz_J0!=NULL)? EXP_uiz_qwv : NULL; + MYCOMPLEX (*pVF_uiz_qwv)[3] = (sum_VF_uiz_J0!=NULL)? VF_uiz_qwv : NULL; + MYCOMPLEX (*pHF_uiz_qwv)[3] = (sum_HF_uiz_J0!=NULL)? HF_uiz_qwv : NULL; + MYCOMPLEX (*pDC_uiz_qwv)[3] = (sum_DC_uiz_J0!=NULL)? DC_uiz_qwv : NULL; + + static const MYINT maxNpt=PTAM_MAX_PEAK_TROUGH; // 波峰波谷的目标 // 根据波峰波谷的目标也给出一个kmax，+5以防万一 @@ -104,11 +231,31 @@ void PTA_method( MYCOMPLEX (*HF_J3)[3][3][4] = (MYCOMPLEX (*)[3][3][4])calloc(nr, sizeof(*HF_J3)); MYCOMPLEX (*DC_J3)[3][3][4] = (MYCOMPLEX (*)[3][3][4])calloc(nr, sizeof(*DC_J3)); + MYCOMPLEX (*EXP_uiz_J3)[3][3][4] = (MYCOMPLEX (*)[3][3][4])calloc(nr, sizeof(*EXP_uiz_J3)); + MYCOMPLEX (*VF_uiz_J3)[3][3][4] = (MYCOMPLEX (*)[3][3][4])calloc(nr, sizeof(*VF_uiz_J3)); + MYCOMPLEX (*HF_uiz_J3)[3][3][4] = (MYCOMPLEX (*)[3][3][4])calloc(nr, sizeof(*HF_uiz_J3)); + MYCOMPLEX (*DC_uiz_J3)[3][3][4] = (MYCOMPLEX (*)[3][3][4])calloc(nr, sizeof(*DC_uiz_J3)); + + MYCOMPLEX (*EXP_uir_J3)[3][3][4] = (MYCOMPLEX (*)[3][3][4])calloc(nr, sizeof(*EXP_uir_J3)); + MYCOMPLEX (*VF_uir_J3)[3][3][4] = (MYCOMPLEX (*)[3][3][4])calloc(nr, sizeof(*VF_uir_J3)); + MYCOMPLEX (*HF_uir_J3)[3][3][4] = (MYCOMPLEX (*)[3][3][4])calloc(nr, sizeof(*HF_uir_J3)); + MYCOMPLEX (*DC_uir_J3)[3][3][4] = (MYCOMPLEX (*)[3][3][4])calloc(nr, sizeof(*DC_uir_J3)); + // 之前求和的值 - MYCOMPLEX (*sum_EXP_J)[3][4] = (MYCOMPLEX (*)[3][4])calloc(nr, sizeof(*sum_EXP_J)); - MYCOMPLEX (*sum_VF_J)[3][4] = (MYCOMPLEX (*)[3][4])calloc(nr, sizeof(*sum_VF_J)); - MYCOMPLEX (*sum_HF_J)[3][4] = (MYCOMPLEX (*)[3][4])calloc(nr, sizeof(*sum_HF_J)); - MYCOMPLEX (*sum_DC_J)[3][4] = (MYCOMPLEX (*)[3][4])calloc(nr, sizeof(*sum_DC_J)); + MYCOMPLEX (*sum_EXP_J)[3][4] = (sum_EXP_J0!=NULL)? (MYCOMPLEX (*)[3][4])calloc(nr, sizeof(*sum_EXP_J)) : NULL; + MYCOMPLEX (*sum_VF_J)[3][4] = (sum_VF_J0!=NULL)? (MYCOMPLEX (*)[3][4])calloc(nr, sizeof(*sum_VF_J)) : NULL; + MYCOMPLEX (*sum_HF_J)[3][4] = (sum_HF_J0!=NULL)? (MYCOMPLEX (*)[3][4])calloc(nr, sizeof(*sum_HF_J)) : NULL; + MYCOMPLEX (*sum_DC_J)[3][4] = (sum_DC_J0!=NULL)? (MYCOMPLEX (*)[3][4])calloc(nr, sizeof(*sum_DC_J)) : NULL; + + MYCOMPLEX (*sum_EXP_uiz_J)[3][4] = (sum_EXP_uiz_J0!=NULL)? (MYCOMPLEX (*)[3][4])calloc(nr, sizeof(*sum_EXP_uiz_J)) : NULL; + MYCOMPLEX (*sum_VF_uiz_J)[3][4] = (sum_VF_uiz_J0!=NULL)? (MYCOMPLEX (*)[3][4])calloc(nr, sizeof(*sum_VF_uiz_J)) : NULL; + MYCOMPLEX (*sum_HF_uiz_J)[3][4] = (sum_HF_uiz_J0!=NULL)? (MYCOMPLEX (*)[3][4])calloc(nr, sizeof(*sum_HF_uiz_J)) : NULL; + MYCOMPLEX (*sum_DC_uiz_J)[3][4] = (sum_DC_uiz_J0!=NULL)? (MYCOMPLEX (*)[3][4])calloc(nr, sizeof(*sum_DC_uiz_J)) : NULL; + + MYCOMPLEX (*sum_EXP_uir_J)[3][4] = (sum_EXP_uir_J0!=NULL)? (MYCOMPLEX (*)[3][4])calloc(nr, sizeof(*sum_EXP_uir_J)) : NULL; + MYCOMPLEX (*sum_VF_uir_J)[3][4] = (sum_VF_uir_J0!=NULL)? (MYCOMPLEX (*)[3][4])calloc(nr, sizeof(*sum_VF_uir_J)) : NULL; + MYCOMPLEX (*sum_HF_uir_J)[3][4] = (sum_HF_uir_J0!=NULL)? (MYCOMPLEX (*)[3][4])calloc(nr, sizeof(*sum_HF_uir_J)) : NULL; + MYCOMPLEX (*sum_DC_uir_J)[3][4] = (sum_DC_uir_J0!=NULL)? (MYCOMPLEX (*)[3][4])calloc(nr, sizeof(*sum_DC_uir_J)) : NULL; // 存储波峰波谷的位置和值 MYREAL (*kEXPpt)[3][4][maxNpt] = (MYREAL (*)[3][4][maxNpt])calloc(nr, sizeof(*kEXPpt)); @@ -124,11 +271,47 @@ void PTA_method( MYINT (*iHFpt)[3][4] = (MYINT (*)[3][4])calloc(nr, sizeof(*iHFpt)); MYINT (*iDCpt)[3][4] = (MYINT (*)[3][4])calloc(nr, sizeof(*iDCpt)); + MYREAL (*kEXPpt_uiz)[3][4][maxNpt] = (MYREAL (*)[3][4][maxNpt])calloc(nr, sizeof(*kEXPpt_uiz)); + MYREAL (*kVFpt_uiz)[3][4][maxNpt] = (MYREAL (*)[3][4][maxNpt])calloc(nr, sizeof(*kVFpt_uiz)); + MYREAL (*kHFpt_uiz)[3][4][maxNpt] = (MYREAL (*)[3][4][maxNpt])calloc(nr, sizeof(*kHFpt_uiz)); + MYREAL (*kDCpt_uiz)[3][4][maxNpt] = (MYREAL (*)[3][4][maxNpt])calloc(nr, sizeof(*kDCpt_uiz)); + MYCOMPLEX (*EXPpt_uiz)[3][4][maxNpt] = (MYCOMPLEX (*)[3][4][maxNpt])calloc(nr, sizeof(*EXPpt_uiz)); + MYCOMPLEX (*VFpt_uiz)[3][4][maxNpt] = (MYCOMPLEX (*)[3][4][maxNpt])calloc(nr, sizeof(*VFpt_uiz)); + MYCOMPLEX (*HFpt_uiz)[3][4][maxNpt] = (MYCOMPLEX (*)[3][4][maxNpt])calloc(nr, sizeof(*HFpt_uiz)); + MYCOMPLEX (*DCpt_uiz)[3][4][maxNpt] = (MYCOMPLEX (*)[3][4][maxNpt])calloc(nr, sizeof(*DCpt_uiz)); + MYINT (*iEXPpt_uiz)[3][4] = (MYINT (*)[3][4])calloc(nr, sizeof(*iEXPpt_uiz)); + MYINT (*iVFpt_uiz)[3][4] = (MYINT (*)[3][4])calloc(nr, sizeof(*iVFpt_uiz)); + MYINT (*iHFpt_uiz)[3][4] = (MYINT (*)[3][4])calloc(nr, sizeof(*iHFpt_uiz)); + MYINT (*iDCpt_uiz)[3][4] = (MYINT (*)[3][4])calloc(nr, sizeof(*iDCpt_uiz)); + + MYREAL (*kEXPpt_uir)[3][4][maxNpt] = (MYREAL (*)[3][4][maxNpt])calloc(nr, sizeof(*kEXPpt_uir)); + MYREAL (*kVFpt_uir)[3][4][maxNpt] = (MYREAL (*)[3][4][maxNpt])calloc(nr, sizeof(*kVFpt_uir)); + MYREAL (*kHFpt_uir)[3][4][maxNpt] = (MYREAL (*)[3][4][maxNpt])calloc(nr, sizeof(*kHFpt_uir)); + MYREAL (*kDCpt_uir)[3][4][maxNpt] = (MYREAL (*)[3][4][maxNpt])calloc(nr, sizeof(*kDCpt_uir)); + MYCOMPLEX (*EXPpt_uir)[3][4][maxNpt] = (MYCOMPLEX (*)[3][4][maxNpt])calloc(nr, sizeof(*EXPpt_uir)); + MYCOMPLEX (*VFpt_uir)[3][4][maxNpt] = (MYCOMPLEX (*)[3][4][maxNpt])calloc(nr, sizeof(*VFpt_uir)); + MYCOMPLEX (*HFpt_uir)[3][4][maxNpt] = (MYCOMPLEX (*)[3][4][maxNpt])calloc(nr, sizeof(*HFpt_uir)); + MYCOMPLEX (*DCpt_uir)[3][4][maxNpt] = (MYCOMPLEX (*)[3][4][maxNpt])calloc(nr, sizeof(*DCpt_uir)); + MYINT (*iEXPpt_uir)[3][4] = (MYINT (*)[3][4])calloc(nr, sizeof(*iEXPpt_uir)); + MYINT (*iVFpt_uir)[3][4] = (MYINT (*)[3][4])calloc(nr, sizeof(*iVFpt_uir)); + MYINT (*iHFpt_uir)[3][4] = (MYINT (*)[3][4])calloc(nr, sizeof(*iHFpt_uir)); + MYINT (*iDCpt_uir)[3][4] = (MYINT (*)[3][4])calloc(nr, sizeof(*iDCpt_uir)); + // 记录点数，当峰谷找到后，清零 MYINT (*gEXPpt)[3][4] = (MYINT (*)[3][4])calloc(nr, sizeof(*gEXPpt)); MYINT (*gVFpt)[3][4] = (MYINT (*)[3][4])calloc(nr, sizeof(*gVFpt)); MYINT (*gHFpt)[3][4] = (MYINT (*)[3][4])calloc(nr, sizeof(*gHFpt)); MYINT (*gDCpt)[3][4] = (MYINT (*)[3][4])calloc(nr, sizeof(*gDCpt)); + + MYINT (*gEXPpt_uiz)[3][4] = (MYINT (*)[3][4])calloc(nr, sizeof(*gEXPpt_uiz)); + MYINT (*gVFpt_uiz)[3][4] = (MYINT (*)[3][4])calloc(nr, sizeof(*gVFpt_uiz)); + MYINT (*gHFpt_uiz)[3][4] = (MYINT (*)[3][4])calloc(nr, sizeof(*gHFpt_uiz)); + MYINT (*gDCpt_uiz)[3][4] = (MYINT (*)[3][4])calloc(nr, sizeof(*gDCpt_uiz)); + + MYINT (*gEXPpt_uir)[3][4] = (MYINT (*)[3][4])calloc(nr, sizeof(*gEXPpt_uir)); + MYINT (*gVFpt_uir)[3][4] = (MYINT (*)[3][4])calloc(nr, sizeof(*gVFpt_uir)); + MYINT (*gHFpt_uir)[3][4] = (MYINT (*)[3][4])calloc(nr, sizeof(*gHFpt_uir)); + MYINT (*gDCpt_uir)[3][4] = (MYINT (*)[3][4])calloc(nr, sizeof(*gDCpt_uir)); for(MYINT ir=0; ir kmax) break; // 计算核函数 F(k, w) - kernel(mod1d, omega, k, pEXP_qwv, pVF_qwv, pHF_qwv, pDC_qwv); + kernel(mod1d, omega, k, pEXP_qwv, pVF_qwv, pHF_qwv, pDC_qwv, + calc_upar, pEXP_uiz_qwv, pVF_uiz_qwv, pHF_uiz_qwv, pDC_uiz_qwv); iendk=true; for(MYINT ir=0; ir` :ref:`(张海明, 2021) `，默认自动选择 + :param calc_upar: 是否计算位移u的空间导数 :param gf_source: 待计算的震源类型 :param statsfile: 波数k积分（包括Filon积分和峰谷平均法）的过程记录文件，常用于debug或者观察积分过程中 :math:`F(k,\omega)` 和 :math:`F(k,\omega)J_m(kr)k` 的变化 :param statsidxs: 仅输出特定频点的过程记录文件，建议给定频点，否则默认所有频率点的记录文件都输出，很占空间 @@ -289,53 +357,41 @@ def gen_gf_spectra( C_DSgrn = ((c_PGRN*3)*nrs)() if calc_DC else None C_SSgrn = ((c_PGRN*3)*nrs)() if calc_DC else None - - EXPgrn:List[List[PyGreenFunction]] = [] - VFgrn:List[List[PyGreenFunction]] = [] - DDgrn:List[List[PyGreenFunction]] = [] - HFgrn:List[List[PyGreenFunction]] = [] - DSgrn:List[List[PyGreenFunction]] = [] - SSgrn:List[List[PyGreenFunction]] = [] - + # 位移u的空间导数 + C_EXPgrn_uiz = C_VFgrn_uiz = C_HFgrn_uiz = C_DDgrn_uiz = C_DSgrn_uiz = C_SSgrn_uiz = None + C_EXPgrn_uir = C_VFgrn_uir = C_HFgrn_uir = C_DDgrn_uir = C_DSgrn_uir = C_SSgrn_uir = None + if calc_upar: + C_EXPgrn_uiz = ((c_PGRN*2)*nrs)() if calc_EXP else None + C_VFgrn_uiz = ((c_PGRN*2)*nrs)() if calc_VF else None + C_HFgrn_uiz = ((c_PGRN*3)*nrs)() if calc_HF else None + C_DDgrn_uiz = ((c_PGRN*2)*nrs)() if calc_DC else None + C_DSgrn_uiz = ((c_PGRN*3)*nrs)() if calc_DC else None + C_SSgrn_uiz = ((c_PGRN*3)*nrs)() if calc_DC else None + # + C_EXPgrn_uir = ((c_PGRN*2)*nrs)() if calc_EXP else None + C_VFgrn_uir = ((c_PGRN*2)*nrs)() if calc_VF else None + C_HFgrn_uir = ((c_PGRN*3)*nrs)() if calc_HF else None + C_DDgrn_uir = ((c_PGRN*2)*nrs)() if calc_DC else None + C_DSgrn_uir = ((c_PGRN*3)*nrs)() if calc_DC else None + C_SSgrn_uir = ((c_PGRN*3)*nrs)() if calc_DC else None + + + EXPgrn, VFgrn, HFgrn, DDgrn, DSgrn, SSgrn = self._init_grn( + distarr, calc_EXP, calc_VF, calc_HF, calc_DC, + C_EXPgrn, C_VFgrn, C_HFgrn, C_DDgrn, C_DSgrn, C_SSgrn, + nt, dt, freqs, wI) - for ir in range(nrs): - dist = distarr[ir] - EXPgrn.append([]) - VFgrn .append([]) - DDgrn .append([]) - HFgrn .append([]) - DSgrn .append([]) - SSgrn .append([]) - for i, comp in enumerate(['Z', 'R', 'T']): - if i<2: - if calc_EXP: - grn = PyGreenFunction(f'EX{comp}', nt, dt, freqs, wI, dist, depsrc, deprcv) - EXPgrn[ir].append(grn) - C_EXPgrn[ir][i] = pointer(grn.c_grn) - - if calc_VF: - grn = PyGreenFunction(f'VF{comp}', nt, dt, freqs, wI, dist, depsrc, deprcv) - VFgrn[ir].append(grn) - C_VFgrn[ir][i] = pointer(grn.c_grn) - - if calc_DC: - grn = PyGreenFunction(f'DD{comp}', nt, dt, freqs, wI, dist, depsrc, deprcv) - DDgrn[ir].append(grn) - C_DDgrn[ir][i] = pointer(grn.c_grn) - - if calc_HF: - grn = PyGreenFunction(f'HF{comp}', nt, dt, freqs, wI, dist, depsrc, deprcv) - HFgrn[ir].append(grn) - C_HFgrn[ir][i] = pointer(grn.c_grn) - - if calc_DC: - grn = PyGreenFunction(f'DS{comp}', nt, dt, freqs, wI, dist, depsrc, deprcv) - DSgrn[ir].append(grn) - C_DSgrn[ir][i] = pointer(grn.c_grn) - grn = PyGreenFunction(f'SS{comp}', nt, dt, freqs, wI, dist, depsrc, deprcv) - SSgrn[ir].append(grn) - C_SSgrn[ir][i] = pointer(grn.c_grn) - + EXPgrn_uiz, VFgrn_uiz, HFgrn_uiz, DDgrn_uiz, DSgrn_uiz, SSgrn_uiz = ([] for _ in range(6)) + EXPgrn_uir, VFgrn_uir, HFgrn_uir, DDgrn_uir, DSgrn_uir, SSgrn_uir = ([] for _ in range(6)) + if calc_upar: + EXPgrn_uiz, VFgrn_uiz, HFgrn_uiz, DDgrn_uiz, DSgrn_uiz, SSgrn_uiz = self._init_grn( + distarr, calc_EXP, calc_VF, calc_HF, calc_DC, + C_EXPgrn_uiz, C_VFgrn_uiz, C_HFgrn_uiz, C_DDgrn_uiz, C_DSgrn_uiz, C_SSgrn_uiz, + nt, dt, freqs, wI, 'z') + EXPgrn_uir, VFgrn_uir, HFgrn_uir, DDgrn_uir, DSgrn_uir, SSgrn_uir = self._init_grn( + distarr, calc_EXP, calc_VF, calc_HF, calc_DC, + C_EXPgrn_uir, C_VFgrn_uir, C_HFgrn_uir, C_DDgrn_uir, C_DSgrn_uir, C_SSgrn_uir, + nt, dt, freqs, wI, 'r') c_statsfile = None @@ -396,6 +452,9 @@ def gen_gf_spectra( self.c_pymod1d, nf1, nf2, nf, c_freqs, nrs, c_rs, wI, vmin_ref, keps, ampk, iwk0, k0, Length, print_runtime, C_EXPgrn, C_VFgrn, C_HFgrn, C_DDgrn, C_DSgrn, C_SSgrn, + calc_upar, + C_EXPgrn_uiz, C_VFgrn_uiz, C_HFgrn_uiz, C_DDgrn_uiz, C_DSgrn_uiz, C_SSgrn_uiz, + C_EXPgrn_uir, C_VFgrn_uir, C_HFgrn_uir, C_DDgrn_uir, C_DSgrn_uir, C_SSgrn_uir, c_statsfile, nstatsidxs, c_statsidxs ) #================================================================================= @@ -434,6 +493,29 @@ def gen_gf_spectra( stream.append(DSgrn [ir][i].freq2time(delayT, travtP, travtS, sgn )) stream.append(SSgrn [ir][i].freq2time(delayT, travtP, travtS, sgn )) + if calc_upar: + if i<2: + if calc_EXP: + stream.append(EXPgrn_uiz[ir][i].freq2time(delayT, travtP, travtS, sgn )) + stream.append(EXPgrn_uir[ir][i].freq2time(delayT, travtP, travtS, sgn )) + if calc_VF: + stream.append(VFgrn_uiz [ir][i].freq2time(delayT, travtP, travtS, sgn )) + stream.append(VFgrn_uir [ir][i].freq2time(delayT, travtP, travtS, sgn )) + if calc_DC: + stream.append(DDgrn_uiz [ir][i].freq2time(delayT, travtP, travtS, sgn )) + stream.append(DDgrn_uir [ir][i].freq2time(delayT, travtP, travtS, sgn )) + + if calc_HF: + stream.append(HFgrn_uiz [ir][i].freq2time(delayT, travtP, travtS, sgn )) + stream.append(HFgrn_uir [ir][i].freq2time(delayT, travtP, travtS, sgn )) + + if calc_DC: + stream.append(DSgrn_uiz [ir][i].freq2time(delayT, travtP, travtS, sgn )) + stream.append(DSgrn_uir [ir][i].freq2time(delayT, travtP, travtS, sgn )) + stream.append(SSgrn_uiz [ir][i].freq2time(delayT, travtP, travtS, sgn )) + stream.append(SSgrn_uir [ir][i].freq2time(delayT, travtP, travtS, sgn )) + + dataLst.append(stream) diff --git a/pygrt/utils.py b/pygrt/utils.py index ac658f61..bd2ae423 100755 --- a/pygrt/utils.py +++ b/pygrt/utils.py @@ -45,7 +45,194 @@ ] -def gen_syn_from_gf_DC(st:Stream, M0:float, strike:float, dip:float, rake:float, az:float): +def gen_syn_from_gf(st:Stream, calc_upar:bool, compute_type:str, M0:float, az:float, **kwargs): + r""" + 一个发起函数，根据不同震源参数，从格林函数中合成理论地震图 + + :param st: 计算好的时域格林函数, :class:`obspy.Stream` 类型 + :param calc_upar: 是否计算位移u的空间导数 + :param M0: 标量地震矩, 单位dyne*cm + :param compute_type: 计算类型，应为以下之一: + 'COMPUTE_EXP'(爆炸源), 'COMPUTE_SF'(单力源), + 'COMPUTE_DC'(双力偶源), 'COMPUTE_MT'(矩张量源) + """ + chs = ['Z', 'R', 'T'] + sacin_prefixes = ["", "z", "r", ""] # 输入通道名 + sacout_prefixes = ["", "z", "r", "t"] # 输入通道名 + srcName = ["EX", "VF", "HF", "DD", "DS", "SS"] + allchs = [tr.stats.channel for tr in st] + + s_compute_type = compute_type.split("_")[1][:2] + + baz = 180 + az + if baz > 360: + baz -= 360 + + calcUTypes = 4 if calc_upar else 1 + + stall = Stream() + + dist = st[0].stats.sac['dist'] + upar_scale:float = 1.0 + for ityp in range(calcUTypes): + if ityp > 1: + upar_scale = 1e-5 + if ityp == 3: + upar_scale /= dist + + srcCoef = set_source_coef(ityp==3, upar_scale, compute_type, M0, az, **kwargs) + + inpref = sacin_prefixes[ityp] + outpref = sacout_prefixes[ityp] + + for c in range(3): + ch = chs[c] + tr:Trace = st[0].copy() + tr.data[:] = 0.0 + tr.stats.channel = kcmpnm = f'{outpref}{s_compute_type}{ch}' + __check_trace_attr_sac(tr, az=az, baz=baz, kcmpnm=kcmpnm) + for k in range(6): + coef = srcCoef[c, k] + if coef==0.0: + continue + + # 读入数据 + channel = f'{inpref}{srcName[k]}{ch}' + if channel not in allchs: + raise ValueError(f"Failed, channel=\"{channel}\" not exists.") + + tr0 = st.select(channel=channel)[0].copy() + + tr.data += coef*tr0.data + + stall.append(tr) + + return stall + + +def set_source_coef( + par_theta:bool, coef:float, compute_type:str, M0:float, az:float, + fZ=None, fN=None, fE=None, + strike=None, dip=None, rake=None, + MT=None, **kwargs): + r""" + 设置不同震源的方向因子矩阵 + + :param par_theta: 是否求对theta的偏导 + :param coef: 比例系数 + :param compute_type: 计算类型，应为以下之一: + 'COMPUTE_EXP'(爆炸源), 'COMPUTE_SF'(单力源), + 'COMPUTE_DC'(双力偶源), 'COMPUTE_MT'(矩张量源) + :param M0: 地震矩 + :param az: 方位角(度) + + - 其他参数根据计算类型可选: + - 单力源需要: fZ, fN, fE, + - 双力偶源需要: strike, dip, rake + - 矩张量源需要: MT=(Mxx, Mxy, Mxz, Myy, Myz, Mzz) + """ + + # 定义常数: 1度对应的弧度值 + DEG1 = np.pi / 180.0 + + azrad = az*DEG1 + caz = np.cos(azrad) + saz = np.sin(azrad) + + src_coef = np.zeros((3, 6), dtype='f8') + + # 计算乘法因子 + if compute_type == 'COMPUTE_SF': + mult = 1e-15 * M0 * coef + else: + mult = 1e-20 * M0 * coef + + # 根据不同计算类型处理 + if compute_type == 'COMPUTE_EXP': + # 爆炸源情况 + src_coef[0, 0] = src_coef[1, 0] = 0.0 if par_theta else mult # Z/R分量 + src_coef[2, 0] = 0.0 # T分量 + + elif compute_type == 'COMPUTE_SF': + # 单力源情况 + # 计算各向异性系数 + A0 = fZ * mult + A1 = (fN * caz + fE * saz) * mult + A4 = (-fN * saz + fE * caz) * mult + + # 设置震源系数矩阵 (公式4.6.20) + src_coef[0, 1] = src_coef[1, 1] = 0.0 if par_theta else A0 # VF, Z/R + src_coef[0, 2] = src_coef[1, 2] = A4 if par_theta else A1 # HF, Z/R + src_coef[2, 1] = 0.0 # VF, T + src_coef[2, 2] = -A1 if par_theta else A4 # HF, T + + elif compute_type == 'COMPUTE_DC': + # 双力偶源情况 (公式4.8.35) + # 计算各种角度值(转为弧度) + stkrad = strike * DEG1 # 走向角 + diprad = dip * DEG1 # 倾角 + rakrad = rake * DEG1 # 滑动角 + therad = azrad - stkrad # 方位角与走向角差 + + # 计算各种三角函数值 + srak = np.sin(rakrad); crak = np.cos(rakrad) + sdip = np.sin(diprad); cdip = np.cos(diprad) + sdip2 = 2.0 * sdip * cdip; cdip2 = 2.0 * cdip**2 - 1.0 + sthe = np.sin(therad); cthe = np.cos(therad) + sthe2 = 2.0 * sthe * cthe; cthe2 = 2.0 * cthe**2 - 1.0 + + # 计算各向异性系数 + A0 = mult * (0.5 * sdip2 * srak) + A1 = mult * (cdip * crak * cthe - cdip2 * srak * sthe) + A2 = mult * (0.5 * sdip2 * srak * cthe2 + sdip * crak * sthe2) + A4 = mult * (-cdip2 * srak * cthe - cdip * crak * sthe) + A5 = mult * (sdip * crak * cthe2 - 0.5 * sdip2 * srak * sthe2) + + # 设置震源系数矩阵 + src_coef[0, 3] = src_coef[1, 3] = 0.0 if par_theta else A0 # DD, Z/R + src_coef[0, 4] = src_coef[1, 4] = A4 if par_theta else A1 # DS, Z/R + src_coef[0, 5] = src_coef[1, 5] = 2.0 * A5 if par_theta else A2 # SS, Z/R + src_coef[2, 3] = 0.0 # DD, T + src_coef[2, 4] = -A1 if par_theta else A4 # DS, T + src_coef[2, 5] = -2.0 * A2 if par_theta else A5 # DS, T + + elif compute_type == 'COMPUTE_MT': + # 矩张量源情况 (公式4.9.7，修改了各向同性项) + # 初始化矩张量分量 + M11, M12, M13, M22, M23, M33 = MT + + # 计算各向同性部分并减去 + Mexp = (M11 + M22 + M33) / 3.0 + M11 -= Mexp + M22 -= Mexp + M33 -= Mexp + + # 计算方位角的2倍角三角函数 + caz2 = np.cos(2 * azrad) + saz2 = np.sin(2 * azrad) + + # 计算各向异性系数 + A0 = mult * ((2.0 * M33 - M11 - M22) / 6.0) + A1 = mult * (- (M13 * caz + M23 * saz)) + A2 = mult * (0.5 * (M11 - M22) * caz2 + M12 * saz2) + A4 = mult * (M13 * saz - M23 * caz) + A5 = mult * (-0.5 * (M11 - M22) * saz2 + M12 * caz2) + + # 设置震源系数矩阵 + src_coef[0, 0] = src_coef[1, 0] = 0.0 if par_theta else mult * Mexp # EX, Z/R + src_coef[0, 3] = src_coef[1, 3] = 0.0 if par_theta else A0 # DD, Z/R + src_coef[0, 4] = src_coef[1, 4] = A4 if par_theta else A1 # DS, Z/R + src_coef[0, 5] = src_coef[1, 5] = 2.0 * A5 if par_theta else A2 # SS, Z/R + src_coef[2, 0] = 0.0 # EX, T + src_coef[2, 3] = 0.0 # DD, T + src_coef[2, 4] = -A1 if par_theta else A4 # DS, T + src_coef[2, 5] = -2.0 * A2 if par_theta else A5 # DS, T + + + return src_coef + + +def gen_syn_from_gf_DC(st:Stream, M0:float, strike:float, dip:float, rake:float, az:float, calc_upar:bool=False): ''' 双力偶源，角度单位均为度 @@ -55,66 +242,16 @@ def gen_syn_from_gf_DC(st:Stream, M0:float, strike:float, dip:float, rake:float, :param dip: 倾角，以水平面往下为正，0<=dip<=90 :param rake: 滑动角，在断层面相对于走向方向逆时针为正，-180<=rake<=180 :param az: 台站方位角，以北顺时针为正，0<=az<=360 + :param calc_upar: 是否计算位移u的空间导数 :return: - **stream** - 三分量ZRT地震图, :class:`obspy.Stream` 类型 ''' - # 检查通道 - all_sets = {'DDZ', 'DDR', 'DSZ', 'DSR', 'DST', 'SSZ', 'SSR', 'SST'} - chs = [tr.stats.channel for tr in st] - if not all_sets.issubset(chs): - raise ValueError("Failed, No green functions of Double Couple source.") - - - az_deg = az - baz_deg = 180 + az - if baz_deg > 360: - baz_deg -= 360 - - strike = math.radians(strike) - dip = math.radians(dip) - rake = math.radians(rake) - az = math.radians(az) - theta = az - strike - srak = math.sin(rake); crak = math.cos(rake) - sdip = math.sin(dip); cdip = math.cos(dip) - sdip2 = 2*sdip*cdip; cdip2 = 2*cdip*cdip - 1 - sthe = math.sin(theta); cthe = math.cos(theta) - sthe2 = 2*sthe*cthe; cthe2 = 2*cthe*cthe - 1 - - M0 *= 1e-20 - - # 公式(4.8.35) - A0 = M0 * (0.5*sdip2*srak) - A1 = M0 * (cdip*crak*cthe - cdip2*srak*sthe) - A2 = M0 * (0.5*sdip2*srak*cthe2 + sdip*crak*sthe2) - A4 = M0 * (- cdip2*srak*cthe - cdip*crak*sthe) - A5 = M0 * (sdip*crak*cthe2 - 0.5*sdip2*srak*sthe2) - - trZ = st.select(channel='DDZ')[0].copy() - trZ.stats.channel = 'DCZ' - __check_trace_attr_sac(trZ, az=az_deg, baz=baz_deg) - trZ.data = A0*st.select(channel='DDZ')[0].data + \ - A1*st.select(channel='DSZ')[0].data + \ - A2*st.select(channel='SSZ')[0].data - - trR = st.select(channel='DDR')[0].copy() - trR.stats.channel = 'DCR' - __check_trace_attr_sac(trR, az=az_deg, baz=baz_deg) - trR.data = A0*st.select(channel='DDR')[0].data + \ - A1*st.select(channel='DSR')[0].data + \ - A2*st.select(channel='SSR')[0].data - - trT = st.select(channel='DST')[0].copy() - trT.stats.channel = 'DCT' - __check_trace_attr_sac(trT, az=az_deg, baz=baz_deg) - trT.data = A4*st.select(channel='DST')[0].data + \ - A5*st.select(channel='SST')[0].data - - return Stream([trR, trT, trZ]) + return gen_syn_from_gf(st, calc_upar, "COMPUTE_DC", M0, az, strike=strike, dip=dip, rake=rake) -def gen_syn_from_gf_SF(st:Stream, S:float, fN:float, fE:float, fZ:float, az:float): + +def gen_syn_from_gf_SF(st:Stream, S:float, fN:float, fE:float, fZ:float, az:float, calc_upar:bool=False): ''' 单力源，力分量单位均为dyne @@ -124,86 +261,30 @@ def gen_syn_from_gf_SF(st:Stream, S:float, fN:float, fE:float, fZ:float, az:floa :param fE: 东向力，向东为正 :param fZ: 垂向力，向下为正 :param az: 台站方位角，以北顺时针为正，0<=az<=360 + :param calc_upar: 是否计算位移u的空间导数 :return: - **stream** - 三分量ZRT地震图, :class:`obspy.Stream` 类型 ''' - # 检查通道 - all_sets = {'VFZ', 'VFR', 'HFZ', 'HFR', 'HFT'} - chs = [tr.stats.channel for tr in st] - if not all_sets.issubset(chs): - raise ValueError("Failed, No green functions of Single Force source.") - - az_deg = az - baz_deg = 180 + az - if baz_deg > 360: - baz_deg -= 360 - - az = math.radians(az) - saz = math.sin(az); caz = math.cos(az) - - S *= 1e-15 - fN *= S - fE *= S - fZ *= S - - # 公式(4.6.20) - trZ = st.select(channel='VFZ')[0].copy() - trZ.stats.channel = 'SFZ' - __check_trace_attr_sac(trZ, az=az_deg, baz=baz_deg) - trZ.data = (fN*caz + fE*saz)*st.select(channel='HFZ')[0].data + fZ*st.select(channel='VFZ')[0].data - - trR = st.select(channel='VFR')[0].copy() - trR.stats.channel = 'SFR' - __check_trace_attr_sac(trR, az=az_deg, baz=baz_deg) - trR.data = (fN*caz + fE*saz)*st.select(channel='HFR')[0].data + fZ*st.select(channel='VFR')[0].data - - trT = st.select(channel='HFT')[0].copy() - trT.stats.channel = 'SFT' - __check_trace_attr_sac(trT, az=az_deg, baz=baz_deg) - trT.data = (- fN*saz + fE*caz)*st.select(channel='HFT')[0].data + return gen_syn_from_gf(st, calc_upar, "COMPUTE_SF", S, az, fN=fN, fE=fE, fZ=fZ) - return Stream([trR, trT, trZ]) - -def gen_syn_from_gf_EXP(st:Stream, M0:float, az:float): +def gen_syn_from_gf_EXP(st:Stream, M0:float, az:float, calc_upar:bool=False): ''' 爆炸源 :param st: 计算好的时域格林函数, :class:`obspy.Stream` 类型 :param M0: 标量地震矩, 单位dyne*cm :param az: 台站方位角，以北顺时针为正，0<=az<=360 [不用于计算] + :param calc_upar: 是否计算位移u的空间导数 :return: - **stream** - 三分量ZRT地震图, :class:`obspy.Stream` 类型 ''' - # 检查通道 - all_sets = {'EXZ', 'EXR'} - chs = [tr.stats.channel for tr in st] - if not all_sets.issubset(chs): - raise ValueError("Failed, No green functions of Explosion source.") - - az_deg = az - baz_deg = 180 + az - if baz_deg > 360: - baz_deg -= 360 - - M0 *= 1e-20 - - trZ = st.select(channel='EXZ')[0].copy() - trZ.stats.channel = 'EXZ' - __check_trace_attr_sac(trZ, az=az_deg, baz=baz_deg) - trZ.data = M0 * st.select(channel='EXZ')[0].data - - trR = st.select(channel='EXR')[0].copy() - trR.stats.channel = 'EXR' - __check_trace_attr_sac(trR, az=az_deg, baz=baz_deg) - trR.data = M0 * st.select(channel='EXR')[0].data + return gen_syn_from_gf(st, calc_upar, "COMPUTE_EXP", M0, az) - return Stream([trR, trZ]) - -def gen_syn_from_gf_MT(st:Stream, M0:float, MT:np.ndarray, az:float): +def gen_syn_from_gf_MT(st:Stream, M0:float, MT:np.ndarray, az:float, calc_upar:bool=False): ''' 矩张量源，单位为dyne*cm @@ -211,68 +292,12 @@ def gen_syn_from_gf_MT(st:Stream, M0:float, MT:np.ndarray, az:float): :param M0: 标量地震矩 :param MT: 矩张量 (M11, M12, M13, M22, M23, M33),下标1,2,3分别代表北向，东向，垂直向上 :param az: 台站方位角，以北顺时针为正，0<=az<=360 + :param calc_upar: 是否计算位移u的空间导数 :return: - **stream** - 三分量ZRT地震图, :class:`obspy.Stream` 类型 ''' - # 检查通道 - all_sets = {'EXZ', 'EXR', 'DDZ', 'DDR', 'DSZ', 'DSR', 'DST', 'SSZ', 'SSR', 'SST'} - chs = [tr.stats.channel for tr in st] - if not all_sets.issubset(chs): - raise ValueError("Failed, No green functions of Explosion or Double Couple source.") - - az_deg = az - baz_deg = 180 + az - if baz_deg > 360: - baz_deg -= 360 - - MT = np.array(MT) - - # 不使用*=，否则会原地修改MT - MT = MT*M0 - - # 去除各向同性分量，作为爆炸源 - Mexp = (MT[0] + MT[3] + MT[5])/3 - trR, trZ = gen_syn_from_gf_EXP(st, Mexp, az) # Mexp在函数内部会乘1e-20 - trT = trR.copy() - trT.data[:] = 0.0 - - Mexp *= 1e-20 - MT = MT*1e-20 - M11, M12, M13, M22, M23, M33 = MT - M11 -= Mexp - M22 -= Mexp - M33 -= Mexp - - az = math.radians(az) - saz = math.sin(az); caz = math.cos(az) - saz2 = 2*saz*caz; caz2 = 2*caz*caz - 1 - - # 公式(4.9.7)，但去除了各向同性的量 - A0 = ((2*M33 - M11 - M22)/6 ) - A1 = (- (M13*caz + M23*saz)) - A2 = (0.5*(M11 - M22)*caz2+ M12*saz2) - A4 = (M13*saz - M23*caz) - A5 = (-0.5*(M11 - M22)*saz2 + M12*caz2) - - trZ.stats.channel = 'MTZ' - __check_trace_attr_sac(trZ, az=az_deg, baz=baz_deg) - trZ.data += A0*st.select(channel='DDZ')[0].data + \ - A1*st.select(channel='DSZ')[0].data + \ - A2*st.select(channel='SSZ')[0].data - - trR.stats.channel = 'MTR' - __check_trace_attr_sac(trR, az=az_deg, baz=baz_deg) - trR.data += A0*st.select(channel='DDR')[0].data + \ - A1*st.select(channel='DSR')[0].data + \ - A2*st.select(channel='SSR')[0].data - - trT.stats.channel = 'MTT' - __check_trace_attr_sac(trT, az=az_deg, baz=baz_deg) - trT.data += A4*st.select(channel='DST')[0].data + \ - A5*st.select(channel='SST')[0].data - - return Stream([trR, trT, trZ]) + return gen_syn_from_gf(st, calc_upar, "COMPUTE_MT", M0, az, MT=MT) def __check_trace_attr_sac(tr:Trace, **kwargs):