Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
82 changes: 28 additions & 54 deletions pygrt/C_extension/include/grt/dynamic/layer.h
Original file line number Diff line number Diff line change
Expand Up @@ -144,86 +144,60 @@ void grt_delay_GRT_matrix(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *
* 计算该层的连接 P-SV 应力位移矢量与垂直波函数的D矩阵(或其逆矩阵),
* 见公式(5.2.19-20)
*
* @param[in] xa P波归一化垂直波数 \f$ \sqrt{1 - (k_a/k)^2} \f$
* @param[in] xb S波归一化垂直波数 \f$ \sqrt{1 - (k_b/k)^2} \f$
* @param[in] kbkb S波水平波数的平方 \f$ k_b^2=(\frac{\omega}{V_b})^2 \f$
* @param[in] mu 剪切模量
* @param[in] omega 角频率
* @param[in] rho 密度
* @param[in] k 波数
* @param[out] D D矩阵(或其逆矩阵)
* @param[in] mod1d 模型结构体指针
* @param[in] iy 层位索引
* @param[in] inverse 是否生成逆矩阵
* @param[in] liquid_invtype 对于液体层,矩阵会有很多零,至少第二列、第四列和第四行均为零;
* 剩余部分根据所选类型进行讨论:
* [1] 其余6项保留, \f$ 2\mu\Omega \f$ 退化为 \f$ - \rho \omega^2 \f$ ;
* [2] 在 [1] 基础上第一行也置零,这用于满足液体层的边界条件;
* 对应逆矩阵使用伪逆。
* @param[out] D D矩阵(或其逆矩阵)
*
*/
void grt_get_layer_D(
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, int liquid_invtype);
void grt_get_layer_D(const GRT_MODEL1D *mod1d, const size_t iy, bool inverse, int liquid_invtype, cplx_t D[4][4]);

/** 子矩阵 D11,函数参数见 get_layer_D 函数 */
void grt_get_layer_D11(
cplx_t xa, cplx_t xb, real_t k, cplx_t D[2][2]);
void grt_get_layer_D11(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]);

/** 子矩阵 D12,函数参数见 get_layer_D 函数 */
void grt_get_layer_D12(
cplx_t xa, cplx_t xb, real_t k, cplx_t D[2][2]);
void grt_get_layer_D12(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]);

/** 子矩阵 D21,函数参数见 get_layer_D 函数 */
void grt_get_layer_D21(
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]);
void grt_get_layer_D21(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]);

/** 子矩阵 D22,函数参数见 get_layer_D 函数 */
void grt_get_layer_D22(
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]);
void grt_get_layer_D22(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]);

/** 子矩阵 D11_uiz,后缀uiz表示连接位移对z的偏导和垂直波函数,函数参数见 get_layer_D 函数 */
void grt_get_layer_D11_uiz(
cplx_t xa, cplx_t xb, real_t k, cplx_t D[2][2]);
void grt_get_layer_D11_uiz(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]);

/** 子矩阵 D12_uiz,函数参数见 get_layer_D 函数 */
void grt_get_layer_D12_uiz(
cplx_t xa, cplx_t xb, real_t k, cplx_t D[2][2]);
void grt_get_layer_D12_uiz(const GRT_MODEL1D *mod1d, const size_t iy, cplx_t D[2][2]);


/**
* 计算该层的连接 SH 应力位移矢量与垂直波函数的 T 矩阵(或其逆矩阵),
* 见公式(5.2.21-22)
*
* @param[in] xb S波归一化垂直波数 \f$ \sqrt{1 - (k_b/k)^2} \f$
* @param[in] mu 剪切模量
* @param[in] omega 角频率
* @param[in] k 波数
* @param[out] T T矩阵(或其逆矩阵)
* @param[in] mod1d 模型结构体指针
* @param[in] iy 层位索引
* @param[in] inverse 是否生成逆矩阵
* @param[out] T T矩阵(或其逆矩阵)
*
*/
void grt_get_layer_T(
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(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(cplx_t xb1, real_t thk, real_t k, cplx_t E[2][2], bool inverse);



/**
* 【未维护,未使用,仅用于内部代码测试】
* 和 calc_RT_PSV(SH) 函数解决相同问题,但没有使用显式推导的公式,而是直接做矩阵运算,
* 函数接口也类似
*/
void grt_RT_matrix_from_4x4(
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, int *stats);
void grt_get_layer_T(const GRT_MODEL1D *mod1d, const size_t iy, bool inverse, cplx_t T[2][2]);


// /**
// * 【未维护,未使用,仅用于内部代码测试】
// * 和 calc_RT_PSV(SH) 函数解决相同问题,但没有使用显式推导的公式,而是直接做矩阵运算,
// * 函数接口也类似
// */
// void grt_RT_matrix_from_4x4(
// 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, int *stats);
Loading