clear

Result = load("absorber_micfroflown_mic_000.mat");
Result = Result.Result;

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

Pmic0 = Result.FrequencyResponseFunction(:,1);
Pmf0 = Result.FrequencyResponseFunction(:,2);
Vmf0 = -(Result.FrequencyResponseFunction(:,3));

Result = load("absorber_micfroflown_mic_001.mat");
Result = Result.Result;

Pmic1 = Result.FrequencyResponseFunction(:,1);
Pmf1 = Result.FrequencyResponseFunction(:,2);
Vmf1 = -Result.FrequencyResponseFunction(:,3);


Result = load("absorber_micfroflown_mic_002.mat");
Result = Result.Result;

Pmic2 = Result.FrequencyResponseFunction(:,1);
Pmf2 = Result.FrequencyResponseFunction(:,2);
Vmf2 = -Result.FrequencyResponseFunction(:,3);
%%
rho0 = 1.2;
c0 = 343.1;
dr = 0.1;
f = (0:N/2)'/N*fs;
w = 2*pi*f;

P1_appr = (Pmic0+Pmic2)./2;
V1_appr = (Pmic2-Pmic0)./(1i*w*rho0*dr) * 2;

figure
subplot(1,2,1)
semilogx(f, 20*log10(abs(V1_appr)));
xlim([20,10e3])
grid on
xlabel('f -> [Hz]')
ylabel('V -> [m/s]')
hold on
semilogx(f, 20*log10(abs(Vmf1)));
legend('Approximation','Microflown');
title('Velocity')

subplot(1,2,2)
semilogx(f, 20*log10(abs(P1_appr)));
xlim([20,10e3])
grid on
xlabel('f -> [Hz]')
ylabel('P -> [Pa]')
hold on
semilogx(f, 20*log10(abs(Pmf1)));
semilogx(f, 20*log10(abs(Pmic1)));
legend('Approximation','Microflown','Direct measurement');
title('Pressure')
%%
Z0 = rho0*c0;
d0 = 1.5e-2;
d1 = 6.5e-2;
d2 = 11.5e-2;
Zm_mf_0 = Pmf0./Vmf0;
Zm_mf_1 = Pmf1./Vmf1;
Zm_mf_2 = Pmf2./Vmf2; 
Zm_twomic = Pmic1./V1_appr;

Zs_mf_0 = (Zm_mf_0 +  1i*Z0*tan(w*d0/c0))  ./  ( Z0 + 1i*Zm_mf_0.*tan(w*d0/c0)  );
Zs_mf_1 = (Zm_mf_1 +  1i*Z0*tan(w*d1/c0))  ./  ( Z0 + 1i*Zm_mf_1.*tan(w*d1/c0)  );
Zs_mf_2 = (Zm_mf_2 +  1i*Z0*tan(w*d2/c0))  ./  ( Z0 + 1i*Zm_mf_2.*tan(w*d2/c0)  );
Zs_mic = ( Zm_twomic +  1i*Z0*tan(w*d1/c0))  ./  ( Z0 + 1i*Zm_twomic.*tan(w*d1/c0)  );

R_mf_0 = (Zs_mf_0 - Z0)./(Zs_mf_0 + Z0);
R_mf_1 = (Zs_mf_1 - Z0)./(Zs_mf_1 + Z0);
R_mf_2 = (Zs_mf_2 - Z0)./(Zs_mf_2 + Z0);
R_mic = (Zs_mic - Z0)./(Zs_mic + Z0);

alpha_mf_0 = 1 - (abs( R_mf_0 )).^2;
alpha_mf_1 = 1 - (abs( R_mf_1  )).^2;
alpha_mf_2 = 1 - (abs( R_mf_2  )).^2;
alpha_mic = 1 - (abs( R_mic )).^2;

figure

subplot(1,2,1)
semilogx(f, abs(R_mf_0));
grid on 
xlabel('f -> [Hz]')
ylabel('R -> []')
hold on
semilogx(f, abs(R_mf_1));
semilogx(f, abs(R_mf_2));
semilogx(f, abs(R_mic));
legend('Microflown_1','Microflown_2','Microflown_3','Microphone')
xlim([100,1e3])
title('Reflection coefficient')



subplot(1,2,2)
semilogx(f, alpha_mf_0);
grid on 
xlabel('f -> [Hz]')
ylabel('\alpha -> []')
hold on
semilogx(f, alpha_mf_1);
semilogx(f, alpha_mf_2);
semilogx(f, alpha_mic);
legend('Microflown_1','Microflown_2','Microflown_3','Microphone')
xlim([100,10e3])
title('Absorption coefficient')
