function [h, hp] = maraca_hankel(lmax, x)
% MARACA_HANKEL Spherical Hankel functions of the second kind and derivatives.
% h(l+1) = h_l^(2)(x) = j_l(x) - i y_l(x), hp(l+1) = d/dx h_l^(2)(x), for
% l = 0..lmax, real x > 0. With the e^{+i omega t} convention used throughout,
% h^(2)(kr) is the outgoing wave.
%
% j_l uses downward (Miller) recursion, stable for x below l; y_l uses upward
% recursion, stable everywhere. Verified against the closed forms for l = 0, 1.

if x <= 0
  error('maraca_hankel: x must be positive');
end

% ---- j_l: downward recursion, normalized by j_0 = sin(x)/x ----
Lstart = lmax + 15 + ceil(x);
jj = zeros(Lstart + 2, 1);            % jj(n+1) holds j_n up to a common scale
jj(Lstart + 2) = 0;
jj(Lstart + 1) = 1e-30;
for n = Lstart:-1:1
  jj(n) = ((2 * n + 1) / x) * jj(n + 1) - jj(n + 2);
  if abs(jj(n)) > 1e250                % renormalize to avoid overflow
    jj = jj * 1e-250;
  end
end
% The explicit second index pins the column orientation. numbl <= 0.4.19
% guessed a row here when the JIT could not prove jj was a column (its
% length is only known at runtime), silently transposing the result; fixed
% upstream, but kept until a release carrying the fix is out.
j = jj(1:lmax + 1, 1) * ((sin(x) / x) / jj(1));

% ---- y_l: upward recursion ----
y = zeros(lmax + 1, 1);
y(1) = -cos(x) / x;
if lmax >= 1
  y(2) = -cos(x) / x^2 - sin(x) / x;
  for n = 1:lmax - 1
    y(n + 2) = ((2 * n + 1) / x) * y(n + 1) - y(n);
  end
end

h = j - 1i * y;

% ---- derivatives: h'_l = h_{l-1} - (l+1)/x h_l, h'_0 = -h_1 ----
hp = zeros(lmax + 1, 1) + 0i;
h1 = j(min(lmax + 1, 2)) - 1i * y(min(lmax + 1, 2));
if lmax == 0
  % need h_1 explicitly for h'_0
  y1 = -cos(x) / x^2 - sin(x) / x;
  j1 = sin(x) / x^2 - cos(x) / x;
  h1 = j1 - 1i * y1;
end
hp(1) = -h1;
for l = 1:lmax
  hp(l + 1) = h(l) - ((l + 1) / x) * h(l + 1);
end
end