clc; clear all; close all hidden

%Rohr mit zwei Stufen

l_2     = 0.2;              %m Länge des mittleren Abschnittes
l_3     = 0.2;              %m
c       = 343;              %m/s Schallgeschw.
S_1     = 0.15^2*pi/4;      %m^2 Querschnitt erstes Rohr
S_2     = 0.11^2*pi/4;      %m^2
S_3     = 0.07^2*pi/4;      %m^2
S_4     = 0.03^2*pi/4;      %m^2
roh     = 1.183892;         %kg/m^3
freq    = linspace(80,2000,2000);%1/s Frequenz
% freq   = 80:100:2000;     %Hz
% freq = 1000:1000:2000
omega   = freq*2*pi;        %Kreisfrequenz


for j = 1:length(omega)
    k(j)  = omega(j)./c;
    z_4_0      = roh * c;
    z_3_l_3    = S_3/S_4 * z_4_0;
    z_3_0(j)   = (z_3_l_3/z_4_0*cos(k(j)*l_3) + i*sin(k(j)*l_3)) / (i*z_3_l_3/z_4_0.*sin(k(j)*l_3) + sin(k(j)*l_3))*z_4_0;
    z_2_l_2(j)    = S_2/S_3 * z_3_0(j);
    z_2_0(j)   = (z_2_l_2(j)/z_4_0.*cos(k(j)*l_2) + i*sin(k(j)*l_2)) ./ (i*z_2_l_2(j)/z_4_0.*sin(k(j)*l_2) + sin(k(j)*l_2))*z_4_0;
    r_E(j) = (S_1/S_2*z_2_0(j)/z_4_0 - 1) / (S_1/S_2*z_2_0(j)/z_4_0 + 1);
    tau(j) = 1 - (abs(r_E(j))).^2;
    R(j) = -10*log10(tau(j));
end


plot(freq,R)
title('Matlab')
xlabel('Frequenz [Hz]')
ylabel('Transmission Loss')