function out = bem_panel_acoustic_excite(varargin)
%BEM_PANEL_ACOUSTIC_EXCITE  Problem 4: oblique plane wave driving a simply-supported panel.
%
%  Excitation is specified as a 1/3-octave band sound pressure level (SPL)
%  spectrum, referenced to 20 microPa.
%
%  USAGE
%    out = bem_panel_acoustic_excite('theta',45,'eta',0.02,'SPL',spl_table)
%    out = bem_panel_acoustic_excite('theta','random','SPL',spl_table)
%
%  where spl_table is an Nx2 matrix  [fc_Hz, SPL_dB]
%  with fc_Hz being standard 1/3-octave centre frequencies (e.g. 10 12.5 16
%  20 25 ... 800 Hz) and SPL_dB the corresponding band levels re 20 uPa.
%
%  If 'SPL' is omitted the function defaults to a flat 100 dB spectrum
%  (100 dB re 20 uPa in every 1/3-octave band) as a demonstration.
%
%  CONVERSION  1/3-octave band SPL --> continuous pressure PSD
%    p2_band(k) = (20e-6)^2 * 10^(SPL(k)/10)          [Pa^2]
%    bw(k)      = fc(k) * (2^(1/6) - 2^(-1/6))        [Hz]  (exact 1/3-oct BW)
%    Gpp(fc(k)) = p2_band(k) / bw(k)                   [Pa^2/Hz]
%  Gpp is then log-log interpolated onto the internal fine frequency grid.
%
%  PHYSICS
%    Blocked pressure on baffled panel: 2*p0*exp(-j*(kx*x+ky*y)).
%    FRF  |Hvv(f)|^2  [(m/s)^2/Pa^2]  computed by modal decomposition
%    with closed-form spatial integrals for simply-supported modes.
%    Svv(f)  = |Hvv(f)|^2 * Gpp(f)            [(m/s)^2/Hz]
%    Saa(f)  = (2*pi*f/g)^2 * Svv(f)          [G^2/Hz]
%
%  INPUTS (name-value pairs)
%    'theta'      incidence angle from normal, deg  [default 45]
%                 OR the string 'random' for diffuse-field average
%    'eta'        structural loss factor             [default 0.02]
%    'SPL'        Nx2 matrix [fc_Hz, SPL_dB]        [default: flat 100 dB]
%    'fmax'       upper frequency limit, Hz          [default 800]
%    'Nm'         number of modes per side           [default 5]
%    'theta_max'  Paris limit for random incidence   [default 78 deg]
%    'Ntheta'     quadrature angles (random)         [default 36]
%    'plot'       1 = plot (default), 0 = suppress
%
%  Panel & fluid inputs (name-value pairs, all optional with defaults):
%    'a'      panel length, m               [default 0.60]
%    'b'      panel width, m               [default 0.40]
%    'h'      panel thickness, m           [default 0.002]
%    'rho_s'  panel material density, kg/m^3  [default 2700]
%    'E'      Young's modulus, Pa           [default 70e9]
%    'nu'     Poisson's ratio               [default 0.33]
%    'c0'     speed of sound in fluid, m/s  [default 343]
%
%  OUTPUT struct fields
%    f         frequency vector, Hz  (2000x1)
%    Hvv_sq    space-avg velocity FRF |H|^2, (m/s)^2/Pa^2
%    Gpp       pressure PSD on fine grid, Pa^2/Hz
%    Svv_psd   one-sided velocity PSD, (m/s)^2/Hz
%    Saa_psd   one-sided acceleration PSD, G^2/Hz
%    vel_rms   overall velocity RMS, m/s
%    grms      overall acceleration GRMS, G
%    fc        1/3-octave centre frequencies used, Hz
%    spl       SPL values used, dB re 20 uPa
%    f_mn      panel natural frequencies, Hz (Nm x Nm)
%    theta     angle used (degrees or 'random')
%    eta       loss factor

%% ---- default 1/3-octave table (flat 100 dB, 10 Hz – 1000 Hz) -----------
fc_default = [10 12.5 16 20 25 31.5 40 50 63 80 ...
             100 125 160 200 250 315 400 500 630 800 1000]';
spl_default = 100 * ones(size(fc_default));
default_spl = [fc_default, spl_default];

%% ---- parse inputs -------------------------------------------------------
p = inputParser;
addParameter(p,'theta',     45);
addParameter(p,'eta',       0.02);
addParameter(p,'SPL',       default_spl);
addParameter(p,'fmax',      800);
addParameter(p,'Nm',        5);
addParameter(p,'theta_max', 78);
addParameter(p,'Ntheta',    36);
addParameter(p,'plot',      1);
% Panel & fluid parameters
addParameter(p,'a',         0.60);
addParameter(p,'b',         0.40);
addParameter(p,'h',         0.002);
addParameter(p,'rho_s',     2700);
addParameter(p,'E',         70e9);
addParameter(p,'nu',        0.33);
addParameter(p,'c0',        343);
parse(p, varargin{:});

theta_in  = p.Results.theta;
eta       = p.Results.eta;
SPL_in    = p.Results.SPL;
fmax      = p.Results.fmax;
Nm        = p.Results.Nm;
theta_max = p.Results.theta_max;
Ntheta    = p.Results.Ntheta;
do_plot   = p.Results.plot;
% Panel & fluid parameters
a     = p.Results.a;
b     = p.Results.b;
h     = p.Results.h;
rho_s = p.Results.rho_s;
E     = p.Results.E;
nu    = p.Results.nu;
c0    = p.Results.c0;

% validate SPL table
if ~(ismatrix(SPL_in) && size(SPL_in,2)==2)
    error('''SPL'' must be an Nx2 matrix [fc_Hz, SPL_dB].');
end
fc_bands = SPL_in(:,1);
spl_vals = SPL_in(:,2);
fmax = fc_bands(end);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[fig_num] = spl_plot(1, 1, fc_bands, spl_vals);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% incidence mode
if ischar(theta_in) || isstring(theta_in)
    if strcmpi(theta_in,'random')
        random_inc = true;
    else
        error('theta must be a number (degrees) or ''random''.');
    end
else
    random_inc   = false;
    theta_single = theta_in;
end

%% ---- fixed constants ----------------------------------------------------
g_n   = 9.80665;   % m/s^2 per G
p_ref = 20e-6;     % Pa  (SPL reference)

D  = E * h^3 / (12*(1-nu^2));
ms = rho_s * h;

%% ---- natural frequencies ------------------------------------------------
m_vec = (1:Nm)';
n_vec = (1:Nm)';
[MM, NN] = meshgrid(m_vec, n_vec);
omega_mn = pi^2 * sqrt(D/ms) .* ((MM/a).^2 + (NN/b).^2);
f_mn     = omega_mn / (2*pi);


%% ---- 1/96-octave frequency vector --------------------------------------

f0 = 10;                    % starting frequency
oct_step = 1/96;            % octave increment

Nf = floor(96*log2(fmax/f0)) + 1;

f = f0 * 2.^((0:Nf-1)'*oct_step);

omega = 2*pi*f;


omega = 2*pi*f;
Nf    = length(f);

%% ---- 1/3-octave SPL --> continuous Gpp on fine grid --------------------
bw_factor = 2^(1/6) - 2^(-1/6);
p2_band   = (p_ref^2) .* 10.^(spl_vals/10);
bw        = fc_bands * bw_factor;
Gpp_bands = p2_band ./ bw;

% log-log interpolation onto the fine grid (extrapolate flat at ends)
log_fc  = log10(fc_bands);
log_Gpp = log10(Gpp_bands);
Gpp = 10.^interp1(log_fc, log_Gpp, log10(f), 'linear', 'extrap');

%% ---- modal mass ---------------------------------------------------------
M_mn = ms * a * b / 4;

%% ---- angle set ----------------------------------------------------------
if random_inc
    theta_vec = linspace(0, theta_max, Ntheta+1);
    theta_vec = 0.5*(theta_vec(1:end-1) + theta_vec(2:end));
    dtheta    = theta_max / Ntheta * pi/180;
    w_raw     = sin(theta_vec * pi/180) * dtheta;
    w_norm    = w_raw / sum(w_raw);
    label_str = sprintf('Random incidence (\\theta_{max}=%g\xb0)', theta_max);
else
    theta_vec = theta_single;
    w_norm    = 1;
    label_str = sprintf('\\theta = %g\xb0', theta_single);
end

%% ---- main frequency loop: compute |Hvv(f)|^2  ---------------------------
Hvv_sq = zeros(Nf, 1);

for fi = 1:Nf
    k0_f   = omega(fi) / c0;
    Hvv_fi = 0;

    for ti = 1:length(theta_vec)
        theta_r = theta_vec(ti) * pi/180;
        kx      = k0_f * sin(theta_r);

        F = zeros(Nm, Nm);
        for m = 1:Nm
            km      = m*pi/a;
            denom_x = kx^2 - km^2;
            if abs(denom_x) < 1e-6*km^2
                Ix = -1j*a/2 * (-1)^m;
            else
                Ix = km * (1 - (-1)^m * exp(-1j*kx*a)) / denom_x;
            end
            for n = 1:Nm
                if mod(n,2)==0
                    Iy = 0;
                else
                    Iy = 2*b / (n*pi);
                end
                F(m,n) = 2 * Ix * Iy;
            end
        end

        denom_mn = M_mn .* (omega_mn.^2 - omega(fi)^2 + 1j*eta.*omega_mn.^2);
        V_mn     = 1j * omega(fi) * F ./ denom_mn;

        Hvv_fi = Hvv_fi + w_norm(ti) * 0.25 * sum(abs(V_mn(:)).^2);
    end

    Hvv_sq(fi) = Hvv_fi;
end

%% ---- response PSDs  -----------------------------------------------------
Svv_psd = Hvv_sq .* Gpp;
Saa_psd = (omega ./ g_n).^2 .* Svv_psd;

%% ---- RMS ----------------------------------------------------------------
vel_rms = sqrt(trapz(f, Svv_psd));
grms    = sqrt(trapz(f, Saa_psd));

%% ---- output -------------------------------------------------------------
out.f        = f;
out.Hvv_sq   = Hvv_sq;
out.Gpp      = Gpp;
out.Svv_psd  = Svv_psd;
out.Saa_psd  = Saa_psd;
out.vel_rms  = vel_rms;
out.grms     = grms;
out.fc       = fc_bands;
out.spl      = spl_vals;
out.f_mn     = f_mn;
out.theta    = theta_in;
out.eta      = eta;
% Echo panel/fluid parameters in output struct
out.a        = a;
out.b        = b;
out.h        = h;
out.rho_s    = rho_s;
out.E        = E;
out.nu       = nu;
out.c0       = c0;

%% ---- helpers ------------------------------------------------------------
f_flat = sort(f_mn(f_mn(:) <= fmax));

    function psd_logplot(fv, psd, lcolor, ttl, sub, ylab, flo, fhi, fnat)
        psd_max   = max(psd(psd > 0));
        psd_floor = psd_max * 1e-6;
        psd_plot  = max(psd, psd_floor);

        loglog(fv, psd_plot, 'Color', lcolor, 'LineWidth', 1.5);
        hold on;
        % for k = 1:length(fnat)
        %    xline(fnat(k), '--', 'Color', [0.8 0.4 0], 'LineWidth', 0.8);
        % end
        ylim([psd_floor, psd_max * 10]);

        [xtt, xTT, iflag] = xtick_label(flo, fhi);
        if iflag
            set(gca, 'XTick', xtt, 'XTickLabel', xTT);
        end
        xlim([flo, fhi]);
        xlabel('Frequency (Hz)', 'FontSize', 12);
        ylabel(ylab,             'FontSize', 12);
        title(ttl,    'FontSize', 13);
        subtitle(sub, 'FontSize', 11);
        legend('Modal solution', 'Location', 'northeast');
        grid on;
        set(gca, 'FontSize', 11, 'XMinorTick', 'off');
    end

    function apply_log_yticks(ylo, yhi)
        if ylo <= 0 || yhi <= 0 || ylo >= yhi; return; end
        ymin_dec = floor(log10(ylo));
        ymax_dec = ceil(log10(yhi));
        yt = 10.^(ymin_dec:ymax_dec);
        yt = yt(yt >= ylo & yt <= yhi);
        if ~isempty(yt)
            set(gca, 'YTick', yt, 'YTickLabel', ...
                arrayfun(@(y) sprintf('%.0e',y), yt, 'UniformOutput', false));
            ylim([ylo, yhi]);
        end
    end

%% ---- plots --------------------------------------------------------------
if do_plot
    flo = f(1);
    md  = 5;

    p2_total    = sum(p2_band(fc_bands <= fmax));
    spl_overall = 10*log10(p2_total / p_ref^2);

    figure('Name','Velocity PSD','Color','w');
    psd_logplot(f, Svv_psd, [0 0.4 0.8], ...
        sprintf('Sound-Excited Panel \x2014 Velocity PSD\n%s,  \\eta = %.3f', label_str, eta), ...
        sprintf('Overall RMS = %.3g m/s', vel_rms), ...
        'Velocity PSD  ((m/s)^{2}/Hz)', ...
        flo, fmax, f_flat);
    [ymax, ymin] = ymax_ymin_md(max(Svv_psd), min(Svv_psd(Svv_psd>0)), md);
    apply_log_yticks(ymin, ymax);

    figure('Name','Acceleration PSD','Color','w');
    psd_logplot(f, Saa_psd, [0.8 0.1 0.1], ...
        sprintf('Sound-Excited Panel \x2014 Acceleration PSD\n%s,  \\eta = %.3f', label_str, eta), ...
        sprintf('Overall GRMS = %.3g G', grms), ...
        'Acceleration PSD  (G^{2}/Hz)', ...
        flo, fmax, f_flat);
    [ymax, ymin] = ymax_ymin_md(max(Saa_psd), min(Saa_psd(Saa_psd>0)), md);
    apply_log_yticks(ymin, ymax);

    figure('Name','Velocity FRF','Color','w');
    psd_logplot(f, Hvv_sq, [0.1 0.6 0.1], ...
        sprintf('Sound-Excited Panel \x2014 Velocity FRF\n%s,  \\eta = %.3f', label_str, eta), ...
        '|H_{vv}|^{2}  space-averaged', ...
        '|H_{vv}|^{2}  ((m/s)^{2}/Pa^{2})', ...
        flo, fmax, f_flat);
    [ymax, ymin] = ymax_ymin_md(max(Hvv_sq), min(Hvv_sq(Hvv_sq>0)), md);
    apply_log_yticks(ymin, ymax);
end

end