function out = bem_radiate_sphere_v5(varargin)
%BEM_RADIATE_SPHERE_V4  Sample problems 1 & 2: radiation from a sphere.
%
%   Fix log:
%     v5: exact_oscillating sign flipped (-j not +j)
%     v4: TWO SIGN CORRECTIONS
%       (a) bem_acoustic_solve: A = 0.5*I - H  (was +H)
%           Standard exterior Neumann BIE: (0.5I - K)*p = G*q
%       (b) bem_field_point: p = sum((g*q - dgdn*p)*ar)  (was dgdn*p - g*q)
%           Standard exterior representation: p = G*q - K*p
%       Both signs were wrong; they partially cancelled for the monopole
%       (pulsating) mode, accidentally giving ~1% error there while giving
%       ~230% error for the dipole (oscillating) mode.
%     v3: vn = U*nrm(:,3) for oscillating (no effect -- not the bug)
%     v2: exact_oscillating /k fix (reverted in v4 -- not needed)
%
%   Example:
%       out = bem_radiate_sphere_v5('mode','pulsating','ka',1.0);
%       out = bem_radiate_sphere_v5('mode','oscillating','ka',1.0);

    o = parse_opts(struct('mode','pulsating','a',1.0,'ka',1.0,'U',1.0, ...
        'rho',1.21,'c0',343,'subdiv',3,'rfield',3.0,'plot',true), varargin);

    a = o.a;  k = o.ka / a;  omega = k * o.c0;
    [V, F, cen, nrm, ar] = icosphere(o.subdiv, a);
    N = size(cen, 1);

    switch lower(o.mode)
        case 'pulsating'
            vn = o.U * ones(N, 1);
        case 'oscillating'
            vn = o.U * nrm(:,3);
        otherwise
            error('mode must be pulsating or oscillating');
    end

    [p_surf, q_surf] = bem_acoustic_solve(cen, nrm, ar, k, vn, o.rho, o.c0);

    th_surf = acos(max(-1,min(1, cen(:,3)/a)));
    switch lower(o.mode)
        case 'pulsating'
            p_ex_surf = exact_pulsating(a, o.U, k, o.rho, o.c0, a) * ones(N,1);
        case 'oscillating'
            p_ex_surf = exact_oscillating(a, o.U, k, o.rho, o.c0, a, th_surf);
    end
    err_surf = 100 * norm(p_surf - p_ex_surf) / norm(p_ex_surf);

    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_surf, q_surf);
        switch lower(o.mode)
            case 'pulsating'
                p_ex_dir(m) = exact_pulsating(a, o.U, k, o.rho, o.c0, rf);
            case 'oscillating'
                p_ex_dir(m) = exact_oscillating(a, o.U, k, o.rho, o.c0, rf, theta(m));
        end
    end
    err_field = 100 * norm(p_bem_dir - p_ex_dir) / norm(p_ex_dir);

    nr = 60;
    rr = linspace(1.05*a, 10*a, nr).';
    p_bem_rad = zeros(nr,1);  p_ex_rad = zeros(nr,1);
    for m = 1:nr
        xf = [0, 0, rr(m)];
        p_bem_rad(m) = bem_field_point(xf, cen, nrm, ar, k, p_surf, q_surf);
        switch lower(o.mode)
            case 'pulsating'
                p_ex_rad(m) = exact_pulsating(a, o.U, k, o.rho, o.c0, rr(m));
            case 'oscillating'
                p_ex_rad(m) = exact_oscillating(a, o.U, k, o.rho, o.c0, rr(m), 0);
        end
    end

    out = struct('N',N,'ka',o.ka,'k',k,'mode',o.mode, ...
        'p_surf',p_surf,'q_surf',q_surf,'err_surf_pct',err_surf, ...
        'theta',theta,'p_bem_dir',p_bem_dir,'p_ex_dir',p_ex_dir, ...
        'err_field_pct',err_field, ...
        'th_surf',th_surf,'p_ex_surf',p_ex_surf, ...
        'rr',rr,'p_bem_rad',p_bem_rad,'p_ex_rad',p_ex_rad);

    fprintf('[%s] ka=%.2f  N=%d  surface err=%.2f%%  field err=%.2f%%\n', ...
        o.mode, o.ka, N, err_surf, err_field);

    if o.plot
        figA = figure('Color','w','Name',sprintf('Radiation: %s ka=%.2f',o.mode,o.ka));
        p0 = figA.Position;
        figA.Position = [p0(1)-0.3*p0(3), p0(2)-0.3*p0(4), 1.2*p0(3), 1.2*p0(4)];
        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)'); ylabel('|p| at r = rfield');
        legend('BEM','Analytic','Location','best'); grid on;
        title(sprintf('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('Directivity (polar, |p|)');

        figB = figure('Color','w','Name',sprintf('Diagnostics: %s ka=%.2f',o.mode,o.ka));
        p0 = figB.Position;
        figB.Position = [p0(1)-0.3*p0(3), p0(2)-0.3*p0(4), 1.4*p0(3), 1.4*p0(4)];
        [ths, is] = sort(th_surf*180/pi);
        subplot(2,2,1);
        plot(ths, abs(p_surf(is)),'b.','MarkerSize',8); hold on;
        plot(ths, abs(p_ex_surf(is)),'r--','LineWidth',1.5);
        xlabel('\theta (deg)'); ylabel('|p_{surf}|');
        legend('BEM','Analytic','Location','best'); grid on;
        title(sprintf('Surface pressure (err %.2f%%)',err_surf));
        subplot(2,2,2);
        loglog(rr/a, abs(p_bem_rad),'b-','LineWidth',2); hold on;
        loglog(rr/a, abs(p_ex_rad),'r--','LineWidth',1.5);
        loglog(rr/a, abs(p_ex_rad(1))*(rr(1)./rr),'k:','LineWidth',1);
        xlabel('r / a'); ylabel('|p| along z'); grid on;
        legend('BEM','Analytic','1/r ref','Location','best');
        title('Radial decay');
        yall = abs([p_bem_rad; p_ex_rad]); yall = yall(yall > 0);
        ylim([10^floor(log10(min(yall))), 10^ceil(log10(max(yall)))]);
        subplot(2,2,3);
        plot(theta*180/pi, unwrap(angle(p_bem_dir))*180/pi,'b-','LineWidth',2); hold on;
        plot(theta*180/pi, unwrap(angle(p_ex_dir))*180/pi,'r--','LineWidth',1.5);
        xlabel('\theta (deg)'); ylabel('phase(p) (deg)'); grid on;
        legend('BEM','Analytic','Location','best');
        title('Field phase at r = rfield');
        subplot(2,2,4);
        plot(ths, angle(p_surf(is))*180/pi,'b.','MarkerSize',8); hold on;
        plot(ths, angle(p_ex_surf(is))*180/pi,'r--','LineWidth',1.5);
        xlabel('\theta (deg)'); ylabel('phase(p_{surf}) (deg)'); grid on;
        legend('BEM','Analytic','Location','best');
        title('Surface phase');

        figure('Color','w','Name',sprintf('Surface map: %s ka=%.2f',o.mode,o.ka));
        patch('Faces',F,'Vertices',V, ...
              'FaceVertexCData',real(p_surf),'FaceColor','flat', ...
              'EdgeColor',[0.5 0.5 0.5],'EdgeAlpha',0.3);
        axis equal vis3d off; view(135,20); camlight; lighting gouraud;
        colormap(jet); cb = colorbar; cb.Label.String = 'Re(p_{surf})';
        title(sprintf('Surface Re(p):  %s, ka=%.2f',o.mode,o.ka));

        % Text just above sphere
        % Sphere radius
        rmax = max(sqrt(sum(V.^2,2)));
        text(0,0,1.35*rmax,...
            sprintf('Real Surface Pressure (Pa) : %s, ka=%.2f',o.mode,o.ka),...
            'HorizontalAlignment','center',...
            'FontSize',11,...
            'FontWeight','bold');


    end
end

% =====================================================================
function [p, q] = bem_acoustic_solve(cen, nrm, ar, k, vn, rho, c0)
% Exterior Neumann BIE: (0.5*I - H)*p = G*q
% Standard form: (c*I - K)*p = V*q  where c=0.5 on smooth surface
    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
    % FIX v4(a): correct sign -- standard exterior BIE is (0.5I - H)*p = G*q
    A = 0.5*eye(N) - H;
    b = 1j * rho * c0 * k * (G * vn);
    p = A \ b;
    q = 1j * k * rho * c0 * vn;
end

function p = bem_field_point(xf, cen, nrm, ar, k, p_surf, q_surf)
% Exterior representation: p(x) = G*q - K*p = integral[g*q - dgdn*p] dS
    r     = xf - cen;
    rn    = sqrt(sum(r.^2,2));
    g     = exp(-1j*k*rn) ./ (4*pi*rn);
    rdotn = sum(r.*nrm,2);
    dgdn  = g .* (rdotn./rn) .* (-1j*k - 1.0./rn);
    % FIX v4(b): correct sign -- standard representation is G*q - K*p
    p = sum((g.*q_surf - dgdn.*p_surf) .* ar);
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 p = exact_pulsating(a,U,k,rho,c0,r)
    jka = 1j*k*a;
    p = rho*c0*U*(jka/(1+jka))*(a./r).*exp(-1j*k*(r-a));
end

function p = exact_oscillating(a,U,k,rho,c0,r,theta)
% n=1 partial wave, outgoing, e^{-jwt} convention
% BC: dp/dr|_{r=a} = j*omega*rho*U*cos(theta)
% => p = j*rho*c0*U*cos(theta) * h1(kr)/h1'(ka)
    p = -1j*rho*c0*U*cos(theta).*(sph_h2(1,k*r) ./ sph_dh2(1,k*a));
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
