clear
close all

rho0 = 1.2;
c0 = 343.1;
Z0 = rho0*c0;
d = 5e-3;
d12 = 12e-2;


Results = load('porous_absorber_000.mat');
Results = Results.Result;

fs = Results.AcqSetup.Fs;
N = Results.AcqSetup.nSamples;

v = -Results.FrequencyResponseFunction(:,1);
p1 = Results.FrequencyResponseFunction(:,2);
p2 = Results.FrequencyResponseFunction(:,3);
f = (0:N/2)'/N*fs;

pM = (p1+p2)/2;
vM = (p1-p2)./(1i*2*pi*f*rho0*d12);
Zm = pM./vM;
k0 = 2*pi*f/c0;
Z_appr = (Zm+1i*tan(k0*d))./(1+1i*Zm.*tan(k0*d))*Z0;
R_appr = (Z_appr - Z0)./(Z_appr + Z0);
alpha_appr = 1-abs(R_appr).^2;

Z_direct = p2./v;
R_direct = (Z_direct - Z0)./(Z_direct + Z0);
alpha_direct = 1-abs(R_direct).^2;

figure
subplot(1,2,1)
plot(f,abs(R_direct));
grid on
xlabel('f -> [Hz]');
ylabel('R -> []');
xlim([20,5e3])
title('Direct measurement')

subplot(1,2,2)
plot(f,alpha_direct);
grid on
xlabel('f -> [Hz]');
ylabel('\alpha -> []');
xlim([20,5e3])
title('Direct measurement')


figure
subplot(1,2,1)
plot(f,abs(R_appr));
grid on
xlabel('f -> [Hz]');
ylabel('R -> []');
xlim([20,5e3])
title('Approximation')


subplot(1,2,2)
plot(f,alpha_appr);
grid on
xlabel('f -> [Hz]');
ylabel('\alpha -> []');
xlim([20,5e3])
title('Approximation')
