Skip to content

FIX: properly handle zero epicentral distance without clamping r #281

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Dengda98 merged 1 commit into from
Aug 5, 2026
Merged

FIX: properly handle zero epicentral distance without clamping r #281

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
4 changes: 2 additions & 2 deletions pygrt/C_extension/include/grt/common/const.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,8 +59,8 @@ typedef double complex cplx_t;
#define GCC_ALWAYS_INLINE __attribute__((always_inline)) ///< gcc编译器不改动内联函数

#define GRT_SWAP(type, a, b) { type temp = a; a = b; b = temp; } ///< 交换两个变量的值
#define GRT_MIN_DISTANCE 1e-5 ///< 最小震中距，用于限制
#define GRT_IS_SMALLE_DISTANCE(r) ((r) <= GRT_MIN_DISTANCE) ///< 判断是否是过小的震中距
#define GRT_ZERO_DISTANCE 1e-8 ///< 判定为零震中距的阈值 (km)
#define GRT_IS_ZERO(r) ((r) <= GRT_ZERO_DISTANCE) ///< 判断震中距是否为零（或数值上视为零）

#define GRT_STRING_FMT "%18s" ///< 字符串输出格式
#define GRT_REAL_FMT "%18.8e" ///< 浮点数输出格式
Expand Down
4 changes: 2 additions & 2 deletions pygrt/C_extension/include/grt/common/coord.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ void grt_rot_zxy2zrt_symtensor2odr(real_t theta, real_t A[6]);
*
* @param[in] theta r轴相对x轴的旋转弧度
* @param[in,out] u 柱坐标下的位移矢量
* @param[in,out] upar 柱坐标下的位移空间偏导
* @param[in] r r轴坐标
* @param[in,out] upar 柱坐标下的位移空间偏导（第三行已是 (1/r)∂_θ 有限部分）
* @param[in] r r 坐标 (cm)；r=0 时联络项 u/r 改用 ∂_r u
*/
void grt_rot_zrt2zxy_upar(const real_t theta, real_t u[3], real_t upar[3][3], const real_t r);
3 changes: 3 additions & 0 deletions pygrt/C_extension/include/grt/dynamic/dyn_postprocess.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,9 @@
* 由动态位移偏导合成应变张量。
* 数组布局：u[分量][采样点]、upar[偏导方向][分量][采样点]、
* res[第二分量][第一分量][采样点]。
*
* ZRT 联络项：r≠0 用 u/r；r=0 改用 ∂_r u（upar[1][*]），
* 与 syn 中轴点处 (1/r)∂_θ 有限部分配套。
*/
void grt_compute_strain(
size_t npts, float dist, float *const u[GRT_CHANNEL_NUM],
Expand Down
3 changes: 3 additions & 0 deletions pygrt/C_extension/include/grt/static/static_postprocess.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,9 @@
*
* 数组布局：u[分量][点]、upar[偏导方向][分量][点]、
* res[第二分量][第一分量][点]。仅写入 res 的上三角分量。
*
* ZRT 联络项：r≠0 用 u/r；r=0 改用 ∂_r u（upar[1][*]），
* 与 syn 中轴点处 (1/r)∂_θ 有限部分配套。
*/
void grt_static_compute_stress(
size_t nx, size_t ny, const real_t *xs, const real_t *ys,
Expand Down
19 changes: 15 additions & 4 deletions pygrt/C_extension/src/common/coord.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -66,6 +66,17 @@ void grt_rot_zrt2zxy_upar(const real_t theta, real_t u[3], real_t upar[3][3], co
real_t cct = ct*ct;
real_t sct = st*ct;

// 变换含联络项 u_r/r、u_θ/r（r 单位 cm）。
// r=0: u/r 联络项改用 ∂_r u_r、∂_r u_θ（s11,s12），与 syn 中 (1/r)∂_θ 有限部分配套。
real_t u1_over_r, u2_over_r;
if(GRT_IS_ZERO(r * 1e-5)){ // cm → km 后再判零
u1_over_r = s11;
u2_over_r = s12;
} else {
u1_over_r = u1/r;
u2_over_r = u2/r;
}

// uz ux uy
// ∂z
// ∂x
Expand All @@ -82,17 +93,17 @@ void grt_rot_zrt2zxy_upar(const real_t theta, real_t u[3], real_t upar[3][3], co
// ∂ uz / ∂ x
upar[1][0] = s10*ct - s20*st;
// ∂ ux / ∂ x
upar[1][1] = s11*cct + s22*sst - (s12+s21)*sct + u1*sst/r + u2*sct/r;
upar[1][1] = s11*cct + s22*sst - (s12+s21)*sct + u1_over_r*sst + u2_over_r*sct;
// ∂ uy / ∂ x
upar[1][2] = s12*cct - s21*sst + (s11-s22)*sct - u1*sct/r + u2*sst/r;
upar[1][2] = s12*cct - s21*sst + (s11-s22)*sct - u1_over_r*sct + u2_over_r*sst;


// ∂ uz / ∂ y
upar[2][0] = s10*st + s20*ct;
// ∂ ux / ∂ y
upar[2][1] = s21*cct - s12*sst + (s11-s22)*sct - u1*sct/r - u2*cct/r;
upar[2][1] = s21*cct - s12*sst + (s11-s22)*sct - u1_over_r*sct - u2_over_r*cct;
// ∂ uy / ∂ y
upar[2][2] = s22*cct + s11*sst + (s12+s21)*sct + u1*cct/r - u2*sct/r;
upar[2][2] = s22*cct + s11*sst + (s12+s21)*sct + u1_over_r*cct - u2_over_r*sct;


// 转矢量
Expand Down
14 changes: 12 additions & 2 deletions pygrt/C_extension/src/dynamic/grt_greenfn.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -295,7 +295,8 @@ printf("\n"
" <length> will be determined automatically\n"
" in program with the criterion (Bouchon, 1980).\n"
" + manually set one POSITIVE <length>, e.g. -L20\n"
" For FIM or SAFIM:\n"
" For FIM or SAFIM (large epicentral distance only;\n"
" zero epicentral distance is not allowed):\n"
" + +l<Flength> defines the dk of the FIM.\n"
" + +a<Ftol> defines the tolerance of the SAFIM.\n"
" you can't set both.\n"
Expand Down Expand Up @@ -845,6 +846,15 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){
fclose(fp);
GRT_SAFE_FREE_PTR(dummy);

// FIM/SAFIM 面向远震中距，公式含 1/r、1/√r，不能用于 r=0（含 kcut 分段）
if(Ctrl->L.FIM.active || Ctrl->L.SAFIM.active){
for(size_t ir=0; ir<Ctrl->R.nr; ++ir){
if(GRT_IS_ZERO(Ctrl->R.rs[ir])){
GRTBadOptionError(L, "FIM/SAFIM cannot be used with zero epicentral distance.");
}
}
}

}


Expand Down Expand Up @@ -884,7 +894,7 @@ int greenfn_main(int argc, char **argv) {
Ctrl->N.winT = Ctrl->N.nt*Ctrl->N.dt;

// 最大震中距
real_t rmax = Ctrl->R.rs[grt_findMax_real_t(Ctrl->R.rs, Ctrl->R.nr)];
real_t rmax = Ctrl->R.rs[grt_findMax_real_t(Ctrl->R.rs, Ctrl->R.nr)];

// 时窗最大截止时刻
real_t tmax = 0.0;
Expand Down
4 changes: 2 additions & 2 deletions pygrt/C_extension/src/dynamic/grt_rotation.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,8 +59,8 @@ void grt_compute_rotation(
const char *chs = rot2ZNE ? GRT_ZNE_CODES : GRT_ZRT_CODES;

for(size_t i=0; i<npts; ++i){
// ZRT 联络项 u_θ/r，1e-5: km→cm
float ut_over_r = u[2][i] / dist * 1e-5f;
// 联络项 u_θ/r（1e-5: km→cm）：r≠0 用 u_θ/r；r=0 改用 ∂_r u_θ
float ut_over_r = GRT_IS_ZERO(dist) ? upar[1][2][i] : (u[2][i] / dist * 1e-5f);

for(int c=0; c<GRT_CHANNEL_NUM; ++c){
for(int c2=c+1; c2<GRT_CHANNEL_NUM; ++c2){
Expand Down
6 changes: 3 additions & 3 deletions pygrt/C_extension/src/dynamic/grt_strain.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,9 +57,9 @@ void grt_compute_strain(
const char *chs = rot2ZNE ? GRT_ZNE_CODES : GRT_ZRT_CODES;

for(size_t i=0; i<npts; ++i){
// ZRT 联络项 u/r，1e-5: km→cm
float ur_over_r = u[1][i] / dist * 1e-5f;
float ut_over_r = u[2][i] / dist * 1e-5f;
// 联络项（1e-5: km→cm）：r≠0 用 u/r；r=0 改用 ∂_r u，与 syn 中 (1/r)∂_θ 有限部分配套
float ur_over_r = GRT_IS_ZERO(dist) ? upar[1][1][i] : (u[1][i] / dist * 1e-5f);
float ut_over_r = GRT_IS_ZERO(dist) ? upar[1][2][i] : (u[2][i] / dist * 1e-5f);

for(int c=0; c<GRT_CHANNEL_NUM; ++c){
for(int c2=c; c2<GRT_CHANNEL_NUM; ++c2){
Expand Down
7 changes: 4 additions & 3 deletions pygrt/C_extension/src/dynamic/grt_stress.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -82,12 +82,13 @@ void grt_compute_stress(
for(size_t i=0; i<nf; ++i) lam_ukk[i] += fwd->W_f[i];
}

// ZRT 联络项 u/r（1e-5: km→cm）；时域先算好再 FFT，避免频域再分支缩放
// 联络项（1e-5: km→cm）：r≠0 用 u/r；r=0 改用 ∂_r u，与 syn 中 (1/r)∂_θ 有限部分配套
// 时域先算好 ur_over_r / ut_over_r，再 FFT，避免频域再分支缩放
float *ur_over_r = (float *)malloc(sizeof(float)*npts);
float *ut_over_r = (float *)malloc(sizeof(float)*npts);
for(size_t i=0; i<npts; ++i){
ur_over_r[i] = u[1][i] / dist * 1e-5f;
ut_over_r[i] = u[2][i] / dist * 1e-5f;
ur_over_r[i] = GRT_IS_ZERO(dist) ? upar[1][1][i] : (u[1][i] / dist * 1e-5f);
ut_over_r[i] = GRT_IS_ZERO(dist) ? upar[1][2][i] : (u[2][i] / dist * 1e-5f);
}

if(!rot2ZNE){
Expand Down
36 changes: 30 additions & 6 deletions pygrt/C_extension/src/dynamic/grt_syn.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -140,6 +140,8 @@ printf("\n"
" -G<grn_path> Green's Functions output directory of module `greenfn`.\n"
"\n"
" -A<azimuth> Azimuth in degree, from source to station.\n"
" Ignored (forced to 0°) when Green's Functions\n"
" have zero epicentral distance.\n"
"\n"
" -S[u]<scale> Scale factor to all kinds of source. \n"
" + For Explosion, Shear and Moment Tensor,\n"
Expand Down Expand Up @@ -535,17 +537,27 @@ static void syn_accum_from_gf(
* syn_upar[偏导方向][分量][采样点]。gf_uiz/gf_uir 在 calc_upar=false
* 时可传 NULL；单个分量指针为 NULL 时跳过该道。
*
* @param[in] azrad 方位角（弧度）
* r=0 时强制 *azrad=0（e_r→N、e_θ→E）并告警；
* 并用 ∂_r 格林函数合成 (1/r)∂_θ 的有限部分（见函数内注释）。
*
* @param[in,out] azrad 方位角（弧度）；r=0 时写回 0
*/
void grt_syn_from_gf(
size_t npts, float dist,
const pfloatChnlGrid gf, const pfloatChnlGrid gf_uiz, const pfloatChnlGrid gf_uir,
GRT_SYN_TYPE computeType, real_t M0, real_t VpVs_ratio, real_t azrad,
GRT_SYN_TYPE computeType, real_t M0, real_t VpVs_ratio, real_t *azrad,
const real_t mchn[GRT_MECHANISM_NUM],
bool rot2ZNE, bool calc_upar,
float *const syn[GRT_CHANNEL_NUM], float *const syn_upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM])
{
const real_t az = azrad;
// r=0：方位角无定义，约定 e_r→N、e_θ→E，故强制 az=0
if(GRT_IS_ZERO(dist)){
GRTRaiseWarning(
"Zero epicentral distance: azimuth is ignored "
"(forced to 0°; e_r→N, e_θ→E).");
*azrad = 0.0;
}
const real_t az = *azrad;

int calcUTypes = calc_upar ? 4 : 1;
realChnlGrid srcRadi = {0};
Expand All @@ -563,14 +575,16 @@ void grt_syn_from_gf(
for(int ityp = 0; ityp < calcUTypes; ++ityp){
real_t upar_scale = 1.0;
// 求位移空间导数时，需调整比例系数（1e-5: km→cm）
// ZRT 协变导数拆两步：此处算 (1/r)∂_θ u，后处理再补 ±u/r。
// ZRT 协变导数拆两步：此处合成 (1/r)∂_θ u，后处理再补 ±u/r。
// r=0 时协变组合仍有限，两项分别换成有限极限（见下方 gf 与后处理）。
if(ityp > 0){
switch (GRT_ZRT_CODES[ityp-1]){
case 'Z': case 'R':
upar_scale = 1e-5;
break;
// (1/r)∂_θ：r≠0 时 scale∝1/r；r=0 时改用 ∂_r GF，scale 仅留 km→cm
case 'T':
upar_scale = 1e-5 / dist;
upar_scale = GRT_IS_ZERO(dist) ? 1e-5 : (1e-5 / dist);
break;
default:
break;
Expand All @@ -582,6 +596,11 @@ void grt_syn_from_gf(
up = gf_uiz;
} else if(ityp == 2){
up = gf_uir;
} else if(ityp == 3 && GRT_IS_ZERO(dist)){
// r=0: 用 ∂_r GF（par_θ 辐射因子）合成 (1/r)∂_θ 的有限部分；
// 后处理中 u/r 联络项改用 ∂_r u，二者合并得有限直角/柱坐标导数。
// 对 u_z：m≥1 ⇒ u_z(0)=0，lim u_z/r=∂_r u_z；m=0 无 ∂_θ；无 ±u_z/r 联络项。
up = gf_uir;
}

memset(srcRadi, 0, sizeof(srcRadi));
Expand Down Expand Up @@ -752,9 +771,14 @@ int syn_main(int argc, char **argv)
grt_syn_from_gf(
(size_t)npts, Ctrl->dist,
gf, calc_upar ? gf_uiz : NULL, calc_upar ? gf_uir : NULL,
Ctrl->computeType, Ctrl->S.M0, Ctrl->VpVs_ratio, Ctrl->A.azrad, Ctrl->mchn,
Ctrl->computeType, Ctrl->S.M0, Ctrl->VpVs_ratio, &Ctrl->A.azrad, Ctrl->mchn,
false, calc_upar, syn, syn_upar);

// C 可能因 r=0 强制 azrad=0，同步方位角头段
Ctrl->A.azimuth = Ctrl->A.azrad / DEG1;
Ctrl->A.backazimuth = Ctrl->A.azimuth + 180.0;
if(Ctrl->A.backazimuth >= 360.0) Ctrl->A.backazimuth -= 360.0;

// 时间函数 / 积分 / 微分（在旋转前，与历史行为一致）
SACTRACE *tfsac = NULL;
if(Ctrl->D.active){
Expand Down
3 changes: 2 additions & 1 deletion pygrt/C_extension/src/integral/dcm.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,8 @@ void grt_dcm_correction(size_t nr, real_t *rs, real_t dk, real_t kcut, K_INTEG *
{
for(size_t ir = 0; ir < nr; ++ir){
real_t r = rs[ir];
if(GRT_IS_SMALLE_DISTANCE(r)) continue;
// 修正系数含 1/r；r=0 跳过（近场极限已在 k_integ 的 Bessel 极限中处理）
if(GRT_IS_ZERO(r)) continue;

real_t c = 1.0 / r;

Expand Down
11 changes: 6 additions & 5 deletions pygrt/C_extension/src/integral/dwm.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,13 +70,14 @@ real_t grt_discrete_integ(
for(size_t ir = 0; ir < nr; ++ir){
if(iendkrs[ir]) continue; // 该震中距下的波数k积分已收敛

// 跳过奇异点
if(depsrc == deprcv && GRT_IS_SMALLE_DISTANCE(rs[ir])) continue;
// 震源-接收点重合且 r=0：格林函数奇异
if(depsrc == deprcv && GRT_IS_ZERO(rs[ir])) continue;

memset(K->SUM, 0, sizeof(cplxIntegGrid));

// 计算被积函数一项 F(k,w)Jm(kr)k
grt_int_Pk(k, rs[ir], (K->applyDCM && GRT_IS_SMALLE_DISTANCE(rs[ir]))? K->QWV_raw : K->QWV, false, K->SUM);
// r=0 时 DCM 解析修正被跳过，故改用未扣除近场的 QWV_raw；近场由 k_integ 极限处理
grt_int_Pk(k, rs[ir], (K->applyDCM && GRT_IS_ZERO(rs[ir]))? K->QWV_raw : K->QWV, false, K->SUM);

iendk0 = true;

Expand All @@ -102,7 +103,7 @@ real_t grt_discrete_integ(
if(K->calc_upar){
// ------------------------------- ui_z -----------------------------------
// 计算被积函数一项 F(k,w)Jm(kr)k
grt_int_Pk(k, rs[ir], (K->applyDCM && GRT_IS_SMALLE_DISTANCE(rs[ir]))? K->QWVz_raw : K->QWVz, false, K->SUM);
grt_int_Pk(k, rs[ir], (K->applyDCM && GRT_IS_ZERO(rs[ir]))? K->QWVz_raw : K->QWVz, false, K->SUM);

// keps不参与计算位移空间导数的积分，背后逻辑认为u收敛，则uiz也收敛
GRT_LOOP_IntegGrid(im, v){
Expand All @@ -111,7 +112,7 @@ real_t grt_discrete_integ(

// ------------------------------- ui_r -----------------------------------
// 计算被积函数一项 F(k,w)Jm(kr)k
grt_int_Pk(k, rs[ir], (K->applyDCM && GRT_IS_SMALLE_DISTANCE(rs[ir]))? K->QWV_raw : K->QWV, true, K->SUM);
grt_int_Pk(k, rs[ir], (K->applyDCM && GRT_IS_ZERO(rs[ir]))? K->QWV_raw : K->QWV, true, K->SUM);

// keps不参与计算位移空间导数的积分，背后逻辑认为u收敛，则uiz也收敛
GRT_LOOP_IntegGrid(im, v){
Expand Down
13 changes: 6 additions & 7 deletions pygrt/C_extension/src/integral/fim.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,9 +29,6 @@ real_t grt_linear_filon_integ(
{
if(k0 + dk0 >= kmax) return k0;

real_t depsrc = mstat->mod1d->depsrc;
real_t deprcv = mstat->mod1d->deprcv;

// 从0开始，存储第二部分Filon积分的结果
K_INTEG *K2 = grt_init_K_INTEG(K->calc_upar, nr);

Expand Down Expand Up @@ -72,8 +69,8 @@ real_t grt_linear_filon_integ(
for(size_t ir = 0; ir < nr; ++ir){
if(iendkrs[ir]) continue; // 该震中距下的波数k积分已收敛

// 跳过奇异点
if(depsrc == deprcv && GRT_IS_SMALLE_DISTANCE(rs[ir])) continue;
// 防御：r=0 应在参数入口已拒绝（FIM 不适于零震中距）
if(GRT_IS_ZERO(rs[ir])) continue;

memset(K2->SUM, 0, sizeof(cplxIntegGrid));

Expand Down Expand Up @@ -180,8 +177,8 @@ real_t grt_linear_filon_integ(
for(size_t ir = 0; ir < nr; ++ir){
real_t r = rs[ir];

// 跳过奇异点
if(depsrc == deprcv && GRT_IS_SMALLE_DISTANCE(rs[ir])) continue;
// 防御：r=0 应在参数入口已拒绝（FIM 不适于零震中距）
if(GRT_IS_ZERO(rs[ir])) continue;

// Gc
grt_int_Pk_filon(k0N, r, true, K->QWV, false, SUM_Gc[iik]);
Expand Down Expand Up @@ -232,6 +229,8 @@ real_t grt_linear_filon_integ(
// 乘上总系数 sqrt(2.0/(PI*r)) / dk0, 除dks0是在该函数外还会再乘dk0, 并将结果加到原数组中
for(size_t ir = 0; ir < nr; ++ir){
real_t r = rs[ir];
// 防御：r=0 应在参数入口已拒绝（总系数含 1/√r）
if(GRT_IS_ZERO(r)) continue;
real_t tmp = sqrt(2.0/(PI*r)) / dk0;

GRT_LOOP_IntegGrid(im, v){
Expand Down
50 changes: 30 additions & 20 deletions pygrt/C_extension/src/integral/k_integ.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -86,34 +86,44 @@ void grt_int_Pk(real_t k, real_t r, const cplxChnlGrid QWV, bool calc_uir, cplxI
{
real_t bjmk[GRT_MORDER_MAX+1] = {0};
real_t kr = k*r;
real_t kr_inv = 1.0/kr;
real_t kcoef = k;

real_t Jmcoef[GRT_MORDER_MAX+1] = {0};

grt_bessel012(kr, &bjmk[0], &bjmk[1], &bjmk[2]);
if(calc_uir){
real_t bjmk0[GRT_MORDER_MAX+1] = {0};
for(int i=0; i<=GRT_MORDER_MAX; ++i) bjmk0[i] = bjmk[i];

if(GRT_IS_SMALLE_DISTANCE(r)){
bjmk[0] = - bjmk0[1];
bjmk[1] = bjmk0[0] - 0.5;
bjmk[2] = bjmk0[1] - 0.25*kr;
Jmcoef[1] = - 0.125 * kr;
// r=0 ⇒ kr≡0（任意 k）。近场项含 J_m(kr)/r、∂_r(J_m/r)，直接取 x=kr→0 极限，
// 避免 1/kr。仅当 GRT_IS_ZERO(r) 时成立，不可对一般小 r、大 k 套用。
if(GRT_IS_ZERO(r)){
if(calc_uir){
// bjmk ← J_m'(0): J0'=0, J1'=1/2, J2'=0
bjmk[0] = 0.0;
bjmk[1] = 0.5;
bjmk[2] = 0.0;
// Jmcoef ← lim (J_m' - J_m/x)/x ；×k² 后即 lim ∂_r(J_m/r)
// m=1 → 0, m=2 → 1/8
Jmcoef[1] = 0.0;
Jmcoef[2] = 0.125;
kcoef = k*k;
} else {
grt_besselp012(kr, &bjmk[0], &bjmk[1], &bjmk[2]);
for(int i=1; i<=GRT_MORDER_MAX; ++i) Jmcoef[i] = kr_inv * (-kr_inv * bjmk0[i] + bjmk[i]);
}

kcoef = k*k;
}
else {
if(GRT_IS_SMALLE_DISTANCE(r)){
// bjmk ← J_m(0): J0=1, J1=J2=0
bjmk[0] = 1.0;
bjmk[1] = 0.0;
bjmk[2] = 0.0;
// Jmcoef ← lim J_m(x)/x ；×k 后即 lim J_m/r
// m=1 → 1/2, m=2 → 0
Jmcoef[1] = 0.5;
Jmcoef[2] = 0.125*kr;
Jmcoef[2] = 0.0;
}
} else {
grt_bessel012(kr, &bjmk[0], &bjmk[1], &bjmk[2]);
if(calc_uir){
real_t bjmk0[GRT_MORDER_MAX+1] = {0};
for(int i=0; i<=GRT_MORDER_MAX; ++i) bjmk0[i] = bjmk[i];
real_t kr_inv = 1.0/kr;
grt_besselp012(kr, &bjmk[0], &bjmk[1], &bjmk[2]);
for(int i=1; i<=GRT_MORDER_MAX; ++i) Jmcoef[i] = kr_inv * (-kr_inv * bjmk0[i] + bjmk[i]);
kcoef = k*k;
} else {
real_t kr_inv = 1.0/kr;
for(int i=1; i<=GRT_MORDER_MAX; ++i) Jmcoef[i] = bjmk[i]*kr_inv;
}
}
Expand Down
Loading
Loading