% air/concrete/air, intensity transmission coefficient

z1 = 415; % air
z2 = 8e6; % concrete
z3 = 415; % air

c1 = 343;
c2 = 3100;
c3 = 343;

f = 20:10:20e3;
k2 = 2 * pi * f / c2;

L = 100e-2; % 1m thick concrete wall
% L = 50e-2; % other thicknesses
% L = 20e-2;
% L = 10e-2;
% L = 5e-2;
% L = 1e-2;
ti = 4 ./ ( 4 * cos(k2 * L).^2 + (z2/z1 + z1/z2)^2 * sin(k2 * L).^2);
plot(f, ti)
xlabel('frequency (Hz)')
ylabel('intensity transmission coefficient')

% 1m thick concrete wall (assuming lossless)
% the first frequency that is transmitted is c2 / 2 / L
f1 = c2 / 2 /L; % f1 = 1550;
% can also try other frequencies
% f1 = 3 * c2 / 2 /L;
% f1 = 2 * c2 / 2 /L;

%f1 = 1549.9999;
%f1 = 1549;

k1 = 2 * pi * f1 / c1;  
k2 = 2 * pi * f1 / c2;  % freq stays the same, sound speed (also wavelength) changes 
k3 = 2 * pi * f1 / c3;

vector = [ 1;  z2; 0; 0];
matrix = [-1 1 1 0; ...
        z2 z1 -z1 0; ...
        0 exp(-j * k2 * L) exp(j * k2 * L) -exp(-j * k3 * L); ...
        0 -z3 * exp(-j * k2 * L) z3 * exp(j * k2 * L) z2 * exp(-j * k3 * L)];
coeff = matrix \ vector;

% vector = [ 1;  1; 0; 0];
% matrix = [-1 1 1 0; ...
%         1 z1/z2 -z1/z2 0; ...
%         0 exp(-j * k2 * L) exp(j * k2 * L) -exp(-j * k3 * L); ...
%         0 -exp(-j * k2 * L)/z2 exp(j * k2 * L)/z2 exp(-j * k3 * L)/z3];
% coeff = matrix \ vector; % solve:  matrix * coeff = vector

rbar = ((1 - z1/z3) * cos(k2 * L) + j * (z2 / z3 - z1/z2) * sin(k2 * L))/...
        ((1 + z1/z3) * cos(k2 * L) + j * (z2 / z3 + z1/z2) * sin(k2 * L));

bbar = 2 / ((1 + z1/z2) * (z2 + z3) / (z3 - z2) * exp(2 * j * k2 * L) + 1 - z1/z2)
abar = (2 * z2 - bbar * (z2 - z1)) / (z2 + z1)
abar = bbar * (z2 + z3) * exp(2 * j * k2 * L) / (z3 - z2)
abar = (z1 + z2 + rbar * (z1 - z2))/2/z1
bbar = abar * exp(-2 * j * k2 * L) * (z3 - z2) / (z3 + z2)
% same thing 
tbar = bbar * exp(j * k3 * L) * exp(j * k2 * L) * (1 + (z2 + z3) / (z3 - z2))
tbar = abar * exp(-j * k2 * L) * 2 * z3 / (z3 + z2) * exp(j * k3 * L)

coeffbar = [rbar; abar; bbar; tbar]
matrix * coeffbar

x1 = -L:(L/1000):0;
x2 = 0:(L/1000):L;
x3 = L:(L/1000):(2 * L);

p1 = (exp(-j * k1 * x1) + rbar * exp(j * k1 * x1));
p2 = (abar * exp(-j * k2 * x2) + bbar * exp(j * k2 * x2));
p3 = (tbar * exp(-j * k3 * x3));
figure(2)
plot(x1, real(p1), 'r', x2, real(p2), 'g', x3, real(p3), 'b')
ylabel('instantaneous pressure')
xlabel('normalized distance (units of L)')
grid
fix_axis

figure(3)
% semilogy(x1, abs(p1), 'r', x2, abs(p2), 'g', x3, abs(p3), 'b')
plot(x1, abs(p1), 'r', x2, abs(p2), 'g', x3, abs(p3), 'b')

ylabel('peak pressure')
xlabel('normalized distance (units of L)')
grid
fix_axis

% amplitude of standing waves in medium 1
max(abs(p1))

% amplitude of traveling waves in medium 3
max(abs(p3))