HypersphericalBesselAiry Module Source file:hyperspherical_bessels_airy


One-point second-order Airy/Olver approximation for hyperspherical Bessel functions.
It solves the reduced equation u'' = lambda^2 q(chi) u + r(chi) u with
u = S_K(chi) phi_l^nu(chi), normalized to the exact origin amplitude.
The local second-order correction is calibrated to about 1e-4 peak-relative
accuracy in its gate against phi_recurs.
Assume L > 0 and other variables already checked for physical limits

Dependencies

MathUtils Defined in MathUtils
constants Defined in constants
HypersphericalBesselUtils Defined in hyperspherical_bessels_airy

Subroutines

airy_zeta_q(K, beta, chi, turn_chi, zeta, q, sin_k)

Liouville variable zeta (positive below the turning point, negative above) from zeta*(dzeta/dchi)^2 = q, together with q = 1/S_K^2 - beta^2. Within 1e-7 of the turning point both are replaced by their exact linear limits, zeta = -a*(chi-chi_t) and q = -a^3*(chi-chi_t) with a = airy_turn_a. Optionally returns S_K(chi) so callers need not recompute it.
- integer intent(in) :: K
- real(dp) intent(in) :: beta
- real(dp) intent(in) :: chi
- real(dp) intent(in) :: turn_chi
- real(dp) intent(out) :: zeta
- real(dp) intent(out) :: q
- real(dp) intent(out), optional :: sin_k
compute_airy_second_norm_fast(l, K, nu, log_norm)

log of the single scalar that fixes the overall Airy normalization, by matching the chi -> 0 limit of amp*Ai(lambda^(2/3) zeta) to the exact origin behaviour u ~ u0 chi^(l+1). Returned as a log because both factors overflow badly at large l.
- integer intent(in) :: l
- integer intent(in) :: K
- real(dp) intent(in) :: nu
- real(dp) intent(out) :: log_norm
eval_second_scaled_poly(c, zscale, zeta, rescaled, b0, a1)

psi(v) = sum_i c(i) * (v/zscale)**i. Then B0(zeta) = sum_i c(i)*(zeta/zscale)**i/(2*i + 1) and A1 follows from A1' = 0.5*(psi*B0 - B0''), A1(0)=0. This keeps the evaluation scaled, avoiding large powers of 1/zeta. rescaled is false when zscale is just zeta, in which case every power of zeta/zscale is one and the power ladder is skipped; the caller knows this exactly, so it is passed in rather than tested for.
- real(dp) intent(in) :: c(0:AIRY_SECOND_FAST_DEG)
- real(dp) intent(in) :: zscale
- real(dp) intent(in) :: zeta
- logical intent(in) :: rescaled
- real(dp) intent(out) :: b0
- real(dp) intent(out) :: a1
interp_power_from_nodes(x, y, c)

Convert interpolation data (x_i,y_i), i = 0..n, to power coefficients c such that p(x) = sum_i c(i)*x**i. Uses Newton divided differences; for n<=4 this is cheaper and better conditioned than a dense solve.
- real(dp) intent(in) :: x(0:AIRY_SECOND_FAST_DEG)
- real(dp) intent(in) :: y(0:AIRY_SECOND_FAST_DEG)
- real(dp) intent(out) :: c(0:AIRY_SECOND_FAST_DEG)
second_coeffs_onepoint_fast(K, beta, turn_chi, chi, zeta, b0, a1)

Build B0(zeta) and A1(zeta) for this one point only. This version avoids both expensive zeta->chi inversion and the small dense Vandermonde solve. It samples psi directly along a same-side chi segment from the turning point, interpolates in scaled zeta using Newton divided differences, and then evaluates the Olver B0/A1 functionals analytically. Very small |zeta| uses a slightly enlarged same-side fit segment to avoid cancellation in airy_psi_from_s near the turn.
- integer intent(in) :: K
- real(dp) intent(in) :: beta
- real(dp) intent(in) :: turn_chi
- real(dp) intent(in) :: chi
- real(dp) intent(in) :: zeta
- real(dp) intent(out) :: b0
- real(dp) intent(out) :: a1

Functions

real(dp)
airy_u_normalized(l, K, nu, achi, ok, log_norm_in)

Fast second-order Airy/Olver one-point approximation to the reduced u = S_K phi. Only the local psi data for the requested point is built, sampling directly in chi between the turning point and achi. achi must already have been folded into the fundamental domain by normalize_chi; the caller applies the resulting parity sign itself. ok is false where the calibrated gate does not accept the point.
- integer intent(in) :: l
- integer intent(in) :: K
- real(dp) intent(in) :: nu
- real(dp) intent(in) :: achi
- logical intent(out), optional :: ok
- real(dp) intent(in), optional :: log_norm_in