function [fn, f_classical] = pipe_acoustics_fea(c, rho, L, nElem, bc)
% PIPE_ACOUSTICS_FEA  1-D acoustic FEA for a gas-filled pipe.
%
% Inputs:
%   c      - Speed of sound (m/s)
%   rho    - Gas mass density (kg/m^3)
%   L      - Pipe length (m)
%   nElem  - Number of equal-length elements
%   bc     - Boundary condition flag
%              1 = Closed-Closed
%              2 = Closed-Open
%              3 = Open-Open
%
% Outputs:
%   fn          - Natural frequencies from FEA (Hz), column vector
%   f_classical - Classical natural frequencies (Hz), column vector
%
% Mode shape plot: first 1-3 modes superimposed on a single figure.
%
% Classical benchmarks:
%   Closed-Closed or Open-Open :  fn = n*c/(2*L),  n = 1,2,3,...
%   Closed-Open               :  fn = (2n-1)*c/(4*L), n = 1,2,3,...
%
% Reference:
%   L. Segerlind, Applied Finite Element Analysis, Wiley, 1984
%   T. Irvine, Acoustic Pressure Oscillation in a Pipe via FEA,
%   Vibrationdata, 2008

% -------------------------------------------------------------------------
% Element length and node count
% -------------------------------------------------------------------------
h      = L / nElem;          % uniform element length (m)
nNodes = nElem + 1;          % total number of nodes

% -------------------------------------------------------------------------
% Elemental matrices
% -------------------------------------------------------------------------
%   [M]_e = (rho*h/6) * [2 1; 1 2]
%   [K]_e = (rho*c^2/h) * [1 -1; -1 1]

Me = (rho * h / 6)   * [2  1;  1  2];
Ke = (rho * c^2 / h) * [1 -1; -1  1];

% -------------------------------------------------------------------------
% Assemble global matrices
% -------------------------------------------------------------------------
K_glob = zeros(nNodes, nNodes);
M_glob = zeros(nNodes, nNodes);

for e = 1:nElem
    dof = [e, e+1];                    % global DOF indices for element e
    K_glob(dof, dof) = K_glob(dof, dof) + Ke;
    M_glob(dof, dof) = M_glob(dof, dof) + Me;
end

% -------------------------------------------------------------------------
% Apply boundary conditions
%   Closed end  -> Neumann (dP/dx = 0): natural BC, no matrix change needed
%   Open end    -> Dirichlet (P = 0):   remove corresponding row/column
% -------------------------------------------------------------------------
switch bc
    case 1  % Closed-Closed: no DOFs removed; filter zero-frequency mode later
        freeDOF = 1:nNodes;
        filterZero = true;

    case 2  % Closed-Open: P = 0 at x = L (last node)
        freeDOF  = 1:(nNodes - 1);
        filterZero = false;

    case 3  % Open-Open: P = 0 at x = 0 (node 1) and x = L (last node)
        freeDOF  = 2:(nNodes - 1);
        filterZero = false;

    otherwise
        error('bc must be 1 (closed-closed), 2 (closed-open), or 3 (open-open).');
end

K_red = K_glob(freeDOF, freeDOF);
M_red = M_glob(freeDOF, freeDOF);

% -------------------------------------------------------------------------
% Solve generalised eigenvalue problem  [K]{P} = omega^2 [M]{P}
% -------------------------------------------------------------------------
[V_red, D] = eig(K_red, M_red);
omega2 = diag(D);

% Sort ascending
[omega2, idx] = sort(omega2);
V_red = V_red(:, idx);

% Remove near-zero (rigid-body) mode for closed-closed case
if filterZero
    keep   = omega2 > 1e-6 * max(omega2);
    omega2 = omega2(keep);
    V_red  = V_red(:, keep);
end

% Remove any negative eigenvalues (numerical noise)
valid  = omega2 > 0;
omega2 = omega2(valid);
V_red  = V_red(:, valid);

fn = sqrt(omega2) / (2 * pi);   % convert rad/s -> Hz

% -------------------------------------------------------------------------
% Classical natural frequencies array
% -------------------------------------------------------------------------
nModes = length(fn);
f_classical = zeros(nModes, 1);

for n = 1:nModes
    switch bc
        case {1, 3}
            f_classical(n) = n * c / (2 * L);
        case 2
            f_classical(n) = (2*n - 1) * c / (4 * L);
    end
end

% -------------------------------------------------------------------------
% Print results and classical benchmarks
% -------------------------------------------------------------------------
bcNames = {'Closed-Closed', 'Closed-Open', 'Open-Open'};
fprintf('\n=== 1-D Pipe Acoustics FEA ===\n');
fprintf('  Boundary condition : %s\n', bcNames{bc});
fprintf('  Speed of sound     : %.2f m/s\n',  c);
fprintf('  Mass density       : %.4f kg/m^3\n', rho);
fprintf('  Pipe length        : %.4f m\n',  L);
fprintf('  Number of elements : %d\n',      nElem);
fprintf('  Element length     : %.6f m\n',  h);
fprintf('\n  Mode   FEA (Hz)   Classical (Hz)   Error (%%)\n');
fprintf('  ----   --------   --------------   ---------\n');

nPrint = min(nModes, 6);
for n = 1:nPrint
    err = 100 * abs(fn(n) - f_classical(n)) / f_classical(n);
    fprintf('  %4d   %8.4f   %14.4f   %9.4f\n', n, fn(n), f_classical(n), err);
end
fprintf('\n');

% -------------------------------------------------------------------------
% Mode shape plot  (up to 3 modes)
% -------------------------------------------------------------------------
nPlot = min(3, nModes);

% Reconstruct full pressure vector at all nodes (applying zero BCs)
x_nodes = linspace(0, L, nNodes);   % node positions (m)

colors  = {'#003f7f', '#00868b', '#cc3300'};   % navy, teal, red
lstyles = {'-',       '--',      '-.'};

figure('Color', 'w', 'Name', '1-D Pipe Acoustic Mode Shapes', ...
       'Position', [100 100 760 460]);

hold on;
legendEntries = cell(1, nPlot);

for m = 1:nPlot

    % Build full pressure mode shape vector
    phi_full = zeros(nNodes, 1);
    phi_red  = V_red(:, m);

    % Normalise so the maximum absolute value = 1
    phi_red  = phi_red / max(abs(phi_red));
    phi_full(freeDOF) = phi_red;

    plot(x_nodes, phi_full, ...
         'Color', colors{m}, ...
         'LineStyle', lstyles{m}, ...
         'LineWidth', 1.5);

    legendEntries{m} = sprintf('Mode %d  (%.4g Hz)', m, fn(m));
end

% Reference zero line
yline(0, 'Color', [0.5 0.5 0.5], 'LineWidth', 0.8, 'LineStyle', ':');

ylim([-1.2 1.2]);

hold off;
grid on;
box on;

xlabel('Position along pipe,  x  (m)', 'FontSize', 12, 'FontWeight', 'bold');
ylabel('Normalised acoustic pressure,  P', 'FontSize', 12, 'FontWeight', 'bold');
title(sprintf('1-D Pipe Acoustic Mode Shapes — %s', bcNames{bc}), ...
    'FontSize', 13, 'FontWeight', 'bold');
subtitle(sprintf('c = %.1f m/s   |   \\rho = %.4f kg/m^3   |   L = %.3f m   |   N_{elem} = %d', ...
    c, rho, L, nElem), ...
    'FontSize', 11);
legend(legendEntries, 'Location', 'best', 'FontSize', 10);

end