-
BJL(L, X, JL)
Branch order:
low l : explicit formulas/small-x series
very small x : zero for l >= 7
deep pre-peak : exponentially small ascending expansion
far x : large-x oscillatory expansion
transition band : recurrence only for moderate l where needed;
otherwise Airy shoulders, peak polynomial,
or side asymptotics
-
integer intent(in) :: L
-
real(dl) intent(in) :: X
-
real(dl) intent(out) :: JL
BJL_deriv(l, x, jl, djl)
d(j_l)/dx via the standard recurrence dj_l/dx = j_{l-1}(x) - (l + 1)*j_l(x)/x,
given the already-computed jl = j_l(x); handles x = 0 and l = 0 as special cases.
integer intent(in) :: l
real(dl) intent(in) :: x
real(dl) intent(in) :: jl
real(dl) intent(out) :: djl
BJL_postpeak_debye(nu, ax, jl)
Oscillatory Debye asymptotic above the turning point, valid above the
Airy band (u=ax/nu - 1 > 0.42, or eta above the gates in BJL).
Identical to the standard two-correction form but with
sqrt(cotb*cosb)/nu = 1/sqrt(sx*ax), nu/cotb = sx,
exp(-expterm) expanded to second order (expterm < ~1e-3 over the
accepted domain, so the expansion error is < ~1e-9), and
nu*acos(cosb) = nu*(pi/2 - asin(cosb)) evaluated via a polynomial fit
of asin(c)/c in c^2 on [0, 0.4975] (max |asin error| 3.7e-11) when in
range, i.e. for ax > ~1.42*nu; acos is only called close to the
turning point.
real(dl) intent(in) :: nu
real(dl) intent(in) :: ax
real(dl) intent(out) :: jl
BJL_recurrence(L, X, JL)
Stable recurrence evaluation of j_l(x), used by BJL in the moderate-l
transition band where the asymptotic expansions are not yet accurate.
ax > l : direct upward three-term recurrence from j_0, j_1 (stable here).
ax <= l: Miller's algorithm - downward recurrence from an order well
above l, rescaling to avoid overflow, then normalized against
the accurately-known j_0 or j_1 (upward recurrence would be
unstable in this regime).
integer intent(in) :: L
real(dl) intent(in) :: X
real(dl) intent(out) :: JL
BJL_uniform_airy_fast(l, x, nu23, jl)
Two-term corrected Olver uniform Airy approximation:
j_l(x) ~ = pref * [ Ai(tau)
+ eps * (P1(tau) Ai(tau) + Q1(tau) Ai'(tau))
+ eps^2* (P2(tau) Ai(tau) + Q2(tau) Ai'(tau)) ]
nu23 must be (l + 0.5)**(2/3) (passed in to avoid a second power;
the caller already has nu**(1/3)).
Intended domain: the BJL shoulder bands, eta = (x-nu)/nu^(1/3) in
[-2.4,-0.65) U (0.65,3.85] (capped so u = x/nu - 1 is in [-0.26, 0.42])
with l >= 26, i.e. tau = nu^(2/3) zeta in roughly [-4.85, 3.3].
The Olver mapping is evaluated through a single polynomial
P(u) = zeta/u (analytic through the turning point),
with ratio = 4 zeta/(1-z^2) = -4 P(u)/(2 + u), fitted on
u in [-0.26, 0.42] with max |delta zeta| ~ 2e-10. This reproduces the
exact log/acos mapping to ~9e-10 peak-normalized over the gate domain
and is ~1.7x faster. Calls outside the fitted u range are inaccurate;
widen the fit if the eta/u gates in BJL change.
The P1..Q2 corrections are a single weighted least-squares fit of
exact j_l residuals over the gate domain for l >= 26 (so they are
effective coefficients absorbing mapping and higher-order truncation
effects, not the analytic Olver functions); the fit is accurate to
~3.5e-6 peak-normalized at l = 26, improving towards high l, on top
of the ~1e-6 airy_fast floor (see below).
integer intent(in) :: l
real(dl) intent(in) :: x
real(dl) intent(in) :: nu23
real(dl) intent(out) :: jl
ExtendBessels(CP, old_num_xx, requested_xmax)
Extend bessel_horner to the new x range without invalidating the old
Horner coefficients.
The final [150, xmax] range is enlarged so that its endpoint is an integer
multiple of the current final interval beyond 150. Rebuilding BessRanges
with this snapped endpoint gives the same final-range spacing as before,
so the old grid remains a prefix of the new grid.
We recompute an overlap window and overwrite the latter part of the old
intervals in that window. This avoids stitching a newly splined interval
onto an old interval whose right endpoint was previously a spline boundary.
CAMBParams :: CP
integer intent(in) :: old_num_xx
real(dl) intent(in) :: requested_xmax
GenerateBessels(lSamp, CP, requested_xmax)
Build the bessel_horner spline table from scratch: set up the x sampling
(BessRanges, denser at low x and coarsening above x = 25), then for each
sampled l evaluate bjl on the (small-x-truncated) grid and spline it into
Horner-form coefficients used for fast evaluation elsewhere (cmbmain.f90).
lSamples :: lSamp
CAMBParams :: CP
real(dl) intent(in) :: requested_xmax
InitSpherBessels(lSamp, CP, max_bessels_l_index, max_bessels_etak)
Ensure the module-level bessel_horner spline table covers lSamp,
max_bessels_l_index and the current accuracy settings, and reaches
x = max_bessels_etak. Reuses the existing table if it already matches
(extending it in x if only xmax has grown), otherwise regenerates it
from scratch.
lSamples :: lSamp
CAMBParams :: CP
integer intent(in) :: max_bessels_l_index
real(dl) intent(in) :: max_bessels_etak
spline_horner_coeffs(y0, y1, d0, d1, h2over6, c1, c2, c3, c4)
Pack one cubic-spline interval [x_i, x_i+1] (endpoint values y0,y1 and
second derivatives d0,d1; h2over6 = (x_i + 1 - x_i)**2/6) into the
Horner-form coefficients stored in bessel_horner, evaluated elsewhere
(e.g. cmbmain.f90) as jl = c1 + w*(c2 + w*(c3 + w*c4)) with
w = (x_i + 1 - x)/(x_i + 1 - x_i) in [0,1] (so w = 0, jl=c1 = y1 at x_i + 1).
real(dl) intent(in) :: y0
real(dl) intent(in) :: y1
real(dl) intent(in) :: d0
real(dl) intent(in) :: d1
real(dl) intent(in) :: h2over6
real(dl) intent(out) :: c1
real(dl) intent(out) :: c2
real(dl) intent(out) :: c3
real(dl) intent(out) :: c4