aa03146Maraca shake simulator: MATLAB model and interactive web appmagland 1function [h, hp] = maraca_hankel(lmax, x)
2% MARACA_HANKEL Spherical Hankel functions of the second kind and derivatives.
3% h(l+1) = h_l^(2)(x) = j_l(x) - i y_l(x), hp(l+1) = d/dx h_l^(2)(x), for
4% l = 0..lmax, real x > 0. With the e^{+i omega t} convention used throughout,
5% h^(2)(kr) is the outgoing wave.
6%
7% j_l uses downward (Miller) recursion, stable for x below l; y_l uses upward
8% recursion, stable everywhere. Verified against the closed forms for l = 0, 1.
10if x <= 0
11 error('maraca_hankel: x must be positive');
12end
14% ---- j_l: downward recursion, normalized by j_0 = sin(x)/x ----
15Lstart = lmax + 15 + ceil(x);
16jj = zeros(Lstart + 2, 1); % jj(n+1) holds j_n up to a common scale
17jj(Lstart + 2) = 0;
18jj(Lstart + 1) = 1e-30;
19for n = Lstart:-1:1
20 jj(n) = ((2 * n + 1) / x) * jj(n + 1) - jj(n + 2);
21 if abs(jj(n)) > 1e250 % renormalize to avoid overflow
22 jj = jj * 1e-250;
23 end
24end
ee7257cNote the upstream numbl fix beside the slice-orientation workaroundmagland 25% The explicit second index pins the column orientation. numbl <= 0.4.19
26% guessed a row here when the JIT could not prove jj was a column (its
27% length is only known at runtime), silently transposing the result; fixed
28% upstream, but kept until a release carrying the fix is out.
aa03146Maraca shake simulator: MATLAB model and interactive web appmagland 29j = jj(1:lmax + 1, 1) * ((sin(x) / x) / jj(1));
31% ---- y_l: upward recursion ----
32y = zeros(lmax + 1, 1);
33y(1) = -cos(x) / x;
34if lmax >= 1
35 y(2) = -cos(x) / x^2 - sin(x) / x;
36 for n = 1:lmax - 1
37 y(n + 2) = ((2 * n + 1) / x) * y(n + 1) - y(n);
38 end
39end
41h = j - 1i * y;
43% ---- derivatives: h'_l = h_{l-1} - (l+1)/x h_l, h'_0 = -h_1 ----
44hp = zeros(lmax + 1, 1) + 0i;
45h1 = j(min(lmax + 1, 2)) - 1i * y(min(lmax + 1, 2));
46if lmax == 0
47 % need h_1 explicitly for h'_0
48 y1 = -cos(x) / x^2 - sin(x) / x;
49 j1 = sin(x) / x^2 - cos(x) / x;
50 h1 = j1 - 1i * y1;
51end
52hp(1) = -h1;
53for l = 1:lmax
54 hp(l + 1) = h(l) - ((l + 1) / x) * h(l + 1);
55end
56end