function out = bem_scatter_sphere(varargin)
%BEM_SCATTER_SPHERE  Sample problem 3: plane-wave scattering by a rigid sphere.
%
%   A plane wave  p_inc = p0 * exp(-j*k*z)  is incident on a rigid sphere.
%   The sphere is acoustically hard, so the TOTAL normal velocity is zero;
%   the scattered field is driven by Neumann data  q_sc = -d(p_inc)/dn.
%   The scattered surface pressure is solved with the direct BEM and the
%   field is compared to the exact Mie (partial-wave) series.
%
%   OUT = BEM_SCATTER_SPHERE('Name', Value, ...)
%
%   Options (defaults in brackets):
%       'a'       : sphere radius (m)             [1.0]
%       'ka'      : nondimensional frequency k*a  [1.0]
%       'p0'      : incident wave amplitude (Pa)  [1.0]
%       'rho'     : fluid density (kg/m^3)        [1.21]
%       'c0'      : sound speed (m/s)             [343]
%       'subdiv'  : icosphere subdivisions        [3]
%       'rfield'  : field radius for sweep (xa)   [5.0] (in units of a)
%       'nterms'  : Mie series terms              [50]
%       'plot'    : true|false make figures       [true]
%
%   OUT struct fields:
%       .N, .ka, .k
%       .p_sc_surf, .q_sc_surf       scattered surface solution
%       .theta, .p_bem_dir, .p_ex_dir  scattered directivity at rfield
%       .p_inc_dir, .p_tot_dir       incident and total pressure at rfield
%       .err_field_pct               field error vs Mie series (%)
%       .ts_bem, .ts_ex              backscatter target strength check
%
%   Example:
%       out = bem_scatter_sphere('ka',1.0);
%       out = bem_scatter_sphere('ka',3.0,'subdiv',4);   % finer mesh

    o = parse_opts(struct('a',1.0,'ka',1.0,'p0',1.0,'rho',1.21,'c0',343, ...
        'subdiv',3,'rfield',5.0,'nterms',50,'plot',true), varargin);

    a = o.a;  k = o.ka / a;
    [~, ~, cen, nrm, ar] = icosphere(o.subdiv, a);
    N = size(cen, 1);

    % ---- incident plane wave and its normal derivative on the surface ----
    pinc     = o.p0 * exp(-1j * k * cen(:,3));
    dpinc_dn = (-1j * k * nrm(:,3)) .* pinc;

    % rigid sphere: scattered Neumann data cancels the incident normal deriv.
    q_sc = -dpinc_dn;
    p_sc_surf = solve_with_q(cen, nrm, ar, k, q_sc);

    % ---- scattered directivity at r = rfield*a ----
    nth = 73;
    theta = linspace(0, pi, nth).';
    rf = o.rfield * a;
    p_bem_dir = zeros(nth,1);
    p_ex_dir  = zeros(nth,1);
    for m = 1:nth
        xf = [rf*sin(theta(m)), 0, rf*cos(theta(m))];
        p_bem_dir(m) = bem_field_point(xf, cen, nrm, ar, k, p_sc_surf, q_sc);
        p_ex_dir(m)  = exact_scatter(a, o.p0, k, rf, theta(m), o.nterms);
    end
    err_field = 100 * norm(p_bem_dir - p_ex_dir) / norm(p_ex_dir);

    % backscatter (theta = pi) target strength check
    ts_bem = abs(p_bem_dir(end));
    ts_ex  = abs(p_ex_dir(end));

    % ---- total pressure = incident + scattered ----
    p_inc_dir = o.p0 * exp(-1j * k * rf * cos(theta));
    p_tot_dir = p_inc_dir + p_bem_dir;

    out = struct('N',N,'ka',o.ka,'k',k, ...
        'p_sc_surf',p_sc_surf,'q_sc_surf',q_sc, ...
        'theta',theta,'p_bem_dir',p_bem_dir,'p_ex_dir',p_ex_dir, ...
        'p_inc_dir',p_inc_dir,'p_tot_dir',p_tot_dir, ...
        'err_field_pct',err_field,'ts_bem',ts_bem,'ts_ex',ts_ex);

    fprintf('[scatter] ka=%.2f  N=%d  field err=%.2f%%  backscatter |p|: BEM=%.4f exact=%.4f\n', ...
        o.ka, N, err_field, ts_bem, ts_ex);

    if o.plot

        % --- Figure 1: Scattered pressure ---
        figure('Color','w','Name',sprintf('Scattered pressure  ka=%.2f',o.ka),'Position',[100 100 950 420]);

        subplot(1,2,1);
        plot(theta*180/pi, abs(p_bem_dir),'b-','LineWidth',2); hold on;
        plot(theta*180/pi, abs(p_ex_dir), 'r--','LineWidth',1.5);
        xlabel('\theta (deg,  0=forward,  180=backscatter)');
        ylabel('|p_{sc}| at r = rfield');
        legend('BEM','Mie series','Location','best');
        grid on;
        title(sprintf('Scattered directivity  (err %.2f%%)', err_field));

        subplot(1,2,2);
        polarplot(theta,  abs(p_bem_dir),'b-','LineWidth',2); hold on;
        polarplot(-theta, abs(p_bem_dir),'b-','LineWidth',2);
        title(sprintf('Scattered |p| (polar)  ka=%.2f', o.ka));

% --- Figure 2: Total pressure ---
        figure('Color','w','Name',sprintf('Total pressure  ka=%.2f',o.ka), ...
               'Position',[100 100 950 420]);

        subplot(1,2,1);
        plot(theta*180/pi, abs(p_tot_dir),'b-','LineWidth',2); hold on;
        plot(theta*180/pi, abs(p_inc_dir),'k:','LineWidth',1.5);
        xlabel('\theta (deg,  0=forward,  180=backscatter)');
        ylabel('|p_{total}| at r = rfield');
        legend('Total (incident + scattered)','Incident only','Location','south');
        grid on;
        title(sprintf('Total pressure directivity  ka=%.2f', o.ka));
        [lim_y, ticks_y] = neat_axes(max(abs(p_tot_dir)));
        ylim([0 lim_y]);
        yticks(ticks_y(ticks_y >= 0));



        subplot(1,2,2);
        h1 = polarplot(theta,  abs(p_tot_dir),'b-','LineWidth',2); hold on;
        polarplot(-theta, abs(p_tot_dir),'b-','LineWidth',2,'HandleVisibility','off');
        h2 = polarplot(theta,  abs(p_inc_dir),'k:','LineWidth',1.5);
        polarplot(-theta, abs(p_inc_dir),'k:','LineWidth',1.5,'HandleVisibility','off');
        legend([h1 h2],{'Total','Incident'},'Location','best');
        title(sprintf('Total |p| (polar)  ka=%.2f', o.ka));


    end
end

% =====================================================================
function p_surf = solve_with_q(cen, nrm, ar, k, q)
% Assemble and solve ( 1/2 I - H ) p = -G q  for a given Neumann data q.
    N = size(cen,1);
    G = zeros(N,N); H = zeros(N,N);
    for i = 1:N
        r  = cen - cen(i,:);
        rn = sqrt(sum(r.^2,2)); rn(i) = 1.0;
        g  = exp(-1j*k*rn)./(4*pi*rn);
        rdotn = sum(r.*nrm,2);
        dgdn  = g.*(rdotn./rn).*(-1j*k - 1.0./rn);
        g(i) = 0; dgdn(i) = 0;
        a_eq = sqrt(ar(i)/pi);
        G(i,:) = (g.*ar).';  G(i,i) = a_eq/2;
        H(i,:) = (dgdn.*ar).'; H(i,i) = 0;
    end
    p_surf = (0.5*eye(N) - H) \ (-(G*q(:)));
end

% =====================================================================
function p = exact_scatter(a,p0,k,r,theta,Nt)
    ka = k*a; ps = 0; x = cos(theta);
    for n = 0:Nt-1
        An = -(2*n+1)*(1j)^(-n)*(sph_djn(n,ka)/sph_dh2(n,ka));
        if n == 0
            Pn = 1;
        else
            Pn = legendre(n, x);
            Pn = Pn(1);
        end
        ps = ps + An*sph_h2(n,k*r)*Pn;
    end
    p = p0*ps;
end

% =====================================================================
function [V,F,cen,nrm,ar] = icosphere(subdiv, radius)
    t = (1+sqrt(5))/2;
    V = [-1 t 0;1 t 0;-1 -t 0;1 -t 0;0 -1 t;0 1 t; ...
          0 -1 -t;0 1 -t;t 0 -1;t 0 1;-t 0 -1;-t 0 1];
    V = V ./ sqrt(sum(V.^2,2));
    F = [1 12 6;1 6 2;1 2 8;1 8 11;1 11 12;2 6 10;6 12 5; ...
         12 11 3;11 8 7;8 2 9;4 10 5;4 5 3;4 3 7;4 7 9;4 9 10; ...
         5 10 6;3 5 12;7 3 11;9 7 8;10 9 2];
    for s = 1:subdiv
        mid = containers.Map('KeyType','char','ValueType','double');
        NF = zeros(size(F,1)*4,3); cnt = 0;
        for f = 1:size(F,1)
            tri = F(f,:);
            ab = get_mid(tri(1),tri(2));
            bc = get_mid(tri(2),tri(3));
            ca = get_mid(tri(3),tri(1));
            NF(cnt+1,:) = [tri(1) ab ca];
            NF(cnt+2,:) = [tri(2) bc ab];
            NF(cnt+3,:) = [tri(3) ca bc];
            NF(cnt+4,:) = [ab bc ca];
            cnt = cnt + 4;
        end
        F = NF;
    end
    V = V * radius;
    cen = (V(F(:,1),:) + V(F(:,2),:) + V(F(:,3),:)) / 3;
    nv  = cross(V(F(:,2),:)-V(F(:,1),:), V(F(:,3),:)-V(F(:,1),:));
    ar  = 0.5 * sqrt(sum(nv.^2,2));
    nrm = nv ./ sqrt(sum(nv.^2,2));
    flip = sum(cen.*nrm,2) < 0; nrm(flip,:) = -nrm(flip,:);
    function idx = get_mid(p1,p2)
        key = sprintf('%d-%d', min(p1,p2), max(p1,p2));
        if isKey(mid,key), idx = mid(key); return; end
        m = (V(p1,:)+V(p2,:))/2; m = m / norm(m);
        V(end+1,:) = m; idx = size(V,1); mid(key) = idx;
    end
end

% =====================================================================
function h  = sph_h2(n,x),  h  = sph_jn(n,x) - 1j*sph_yn(n,x);        end
function dh = sph_dh2(n,x), dh = sph_djn(n,x) - 1j*sph_dyn(n,x);      end
function j  = sph_jn(n,x),  j  = sqrt(pi./(2*x)).*besselj(n+0.5,x);   end
function y  = sph_yn(n,x),  y  = sqrt(pi./(2*x)).*bessely(n+0.5,x);   end
function dj = sph_djn(n,x), dj = sph_jn(n-1,x) - (n+1)./x.*sph_jn(n,x); end
function dy = sph_dyn(n,x), dy = sph_yn(n-1,x) - (n+1)./x.*sph_yn(n,x); end

function o = parse_opts(o, args)
    for i = 1:2:numel(args)
        key = args{i};
        if ~isfield(o, key), error('Unknown option: %s', key); end
        o.(key) = args{i+1};
    end
end
function [lim, ticks] = neat_axes(v)
% Clean upper axis limit and tick positions. Adapts ytick logic by Tom Irvine.
lim   = 0;
ticks = 0;
for i = 1:30
    X = 10^(i-20);
    if v <= 0.13*X
        lim = 0.15*X; ticks = lim*[0  1/3  2/3  1]; return
    end
    if v <= 0.17*X
        lim = 0.2*X;  ticks = lim*[0  0.5  1];       return
    end
    if v <= 0.26*X
        lim = 0.3*X;  ticks = lim*[0  1/3  2/3  1];  return
    end
    if v <= 0.375*X
        lim = 0.4*X;  ticks = lim*[0  0.5  1];       return
    end
    if v <= 0.48*X
        lim = 0.5*X;  ticks = lim*[0  0.2  0.4  0.6  0.8  1]; return
    end
    if v <= 0.96*X
        lim = 1.0*X;  ticks = lim*[0  0.5  1];       return
    end
end
end