%---------------------------------------------------------
% D_SOLUTION.M
% Solving the FE system of equations for a given excitation
% Fiala Peter
% 2008 BME, Laboratory of Acoustics
%---------------------------------------------------------
clear all;
close all;
clc;
%% Load necessary parameters
load data/parameters Lx Ly Lz Le
load data/boxmesh box DOF freeDOF
load data/boxmodes phi om
%% Excitation force vector
coord = [Lx, 0.5*Ly, 0.5*Lz];    % load coordinate
node = feutil(['findnode rad < ' num2str([Le/2 coord])], box.Node, box.Elt);
DOFind = fe_c(freeDOF, node+0.01, 'ind');
f = zeros(length(freeDOF),1);
f(DOFind) = 1;
phif = phi.' * f;
%% Solution in the frequency domain
Nom = 1024;
fmax = 1000;
eta = 0.01;
omvec = linspace(0, 2*pi*fmax, Nom+1);
Lambda = diag(om.^2);
I = eye(length(om));
alpha = zeros(length(om), length(omvec));
wb = waitbar(0, 'Computing modal coordinates');
for iom = 1 : length(omvec)
    alpha(:,iom) = (Lambda*(1+i*eta)-omvec(iom).^2*I)\phif;
    waitbar(iom/length(omvec), wb);
end
close(wb);
save data/solution omvec coord f alpha
%% frequency content of the modal coordinates
figure(2);
plot(omvec/2/pi, 20*log10(abs(alpha(1:20,:))));
xlabel('Frequency [Hz]');
ylabel('Modal coordinates [dB]');
%% response shape at a single frequency
figure(1);
freq = 800;
[res, n] = min(abs(omvec-2*pi*freq));
plot_free_model(SDT2Free(box), freeDOF, phi*alpha(:,n), 1);
xlabel('x [m]');
ylabel('y [m]');
zlabel('z [m]');
axis equal;
axis tight;
view(15, 25);
%% frequency content of the response
figure(2);
plot(omvec/2/pi, 20*log10(abs(phi(DOFind,:)*alpha)));
xlabel('Frequency [Hz]');
ylabel('Displacement [dB re 1 m/Hz]');