%% RUNME
% Computes and plots the frequency response of a 3d poroacoustic model
% presented in R. Rumpler, P. Göransson, J.-F. Deü. "A finite element
% approach combining a reduced-order system, Padé approximants, and an
% adaptive frequency windowing for fast multi-frequency solution of
% poro-acoustic problems", Int. J. Numer. Methods Eng. 97: 759--784, 2014.
% For more information, see
% http://modelreduction.org/index.php/Porous_absorber.

% Copyright 2023 Quirin Aumann
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in the
% documentation and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 
% "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
% TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
% PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER
% OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
% EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

clear; close all

load('poroacoustic_n386076m1q1')

recompute = true;
if recompute == true
    % Material parameters
    density_fluid           = 1.21;
    viscosity_fluid         = 1.84e-5;
    standard_pressure_fluid = 101e3;
    heat_capacity_fluid     = 1.4;
    prandtl_number_fluid    = 0.71;
    porosity                = 0.96;
    tortuosity              = 1.7;
    flow_resistivity        = 32e3;
    viscous_length          = 90e-6;
    thermal_length          = 165e-6;

    % Frequency-dependent functions.
    viscous_drag = @(s) flow_resistivity * porosity^2 ...
        * sqrt(1 + 4 * s * tortuosity^2 * viscosity_fluid * density_fluid / ...
        (flow_resistivity^2 * viscous_length^2 * porosity^2));

    apparent_mass_density = porosity * density_fluid * (tortuosity - 1);

    alpha = @(s) 1 + 8 * viscosity_fluid / (s * prandtl_number_fluid ...
        * thermal_length^2 * density_fluid) ...
        * sqrt(1 + s * prandtl_number_fluid * thermal_length^2 ...
        * density_fluid / (16 * viscosity_fluid));

    fun_K_p1 = @(s) -porosity * (-apparent_mass_density - viscous_drag(s) / s) ...
        / (porosity * density_fluid + apparent_mass_density ...
        + viscous_drag(s) / s);
    fun_K_p2 = @(s) porosity^2 / (porosity * density_fluid ...
        + apparent_mass_density + viscous_drag(s) / s);
    fun_M_p1 = @(s) s^2 * (viscous_drag(s) / s - ...
        (-apparent_mass_density - viscous_drag(s) / s)^2 / ...
        (porosity * density_fluid + apparent_mass_density ...
        + viscous_drag(s) / s));
    fun_M_p2 = @(s) s^2 * (-(heat_capacity_fluid-1) / (alpha(s) ...
        * heat_capacity_fluid * standard_pressure_fluid));
    fun_M_p3 = @(s) s^2 * fun_K_p1(s);

    mag = zeros(1,length(f));
    for ii=1:length(f)
        fprintf('Frequency step %d, f=%.2f Hz ... ',ii,f)
        tic
        s = 2i*pi*f(ii);
        res = C * ( (...
            s^2 * M + K + ...
            fun_K_p1(s) * K_p1 + ...
            fun_K_p2(s) * K_p2 + ...
            fun_M_p1(s) * M_p1 + ...
            fun_M_p2(s) * M_p2 + ...
            fun_M_p3(s) * M_p3 ...
            ) \ ((-s)^2 * B));
        mag(ii) = 10*log10(abs(res)/1e-9);
        fprintf('finished in %.2f s\n',toc)
    end
end

figure('name','Transfer function')
plot(f,mag)
xlabel('Frequency [Hz]')
ylabel('Magnitude [dB]')
