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
25% The explicit second index works around a numbl slicing bug that can hand
26% back jj(1:lmax+1) as a row vector depending on what else the script calls.
27j = jj(1:lmax + 1, 1) * ((sin(x) / x) / jj(1));
29% ---- y_l: upward recursion ----
30y = zeros(lmax + 1, 1);
31y(1) = -cos(x) / x;
32if lmax >= 1
33 y(2) = -cos(x) / x^2 - sin(x) / x;
34 for n = 1:lmax - 1
35 y(n + 2) = ((2 * n + 1) / x) * y(n + 1) - y(n);
36 end
37end
39h = j - 1i * y;
41% ---- derivatives: h'_l = h_{l-1} - (l+1)/x h_l, h'_0 = -h_1 ----
42hp = zeros(lmax + 1, 1) + 0i;
43h1 = j(min(lmax + 1, 2)) - 1i * y(min(lmax + 1, 2));
44if lmax == 0
45 % need h_1 explicitly for h'_0
46 y1 = -cos(x) / x^2 - sin(x) / x;
47 j1 = sin(x) / x^2 - cos(x) / x;
48 h1 = j1 - 1i * y1;
49end
50hp(1) = -h1;
51for l = 1:lmax
52 hp(l + 1) = h(l) - ((l + 1) / x) * h(l + 1);
53end
54end