%
disp(' ');
disp(' transmission_five_media.m   ver 1.1  September 23, 2009');
disp(' ');
disp(' by Tom Irvine ');
disp(' ');
disp(' Plane Acoustic Wave Transmission Loss through Five Media ');
disp(' ');
%
clear R;
clear c;
clear cc5;
clear rho;
clear f;
clear alpha_t;
clear k;
clear L;
clear TL;
clear aaa;
%
c=zeros(5,1);
L=zeros(5,1);
R=zeros(5,1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Acoustic Impedance values Pa s/m
%
for(i=1:5)
    out1=sprintf('\n Select Medium %d:  1=air  2=aluminum  3=graphite epoxy  4=other',i);
    disp(out1);
    imat=input(' ');
%
    if(i>=2 & i<=4)
        disp(' ');
        disp(' Enter thickness (inch): ');
        L(i)=input(' ');
    end
%
    if(imat==1)  % air
       R(i)=0.018351;    % slugs/(in^2 sec) 
       c(i)=1128*12;
    end
%
    if(imat==2)  % aluminum
        R(i)=614.65;     % slugs/(in^2 sec)
        c(i)=200000;
    end  
%
    if(imat==3)  % graphite epoxy
        R(i)=465.19;     % slugs/(in^2 sec)
        c(i)=21499*12;      
    end    
%
    if(imat==4)  % other
        disp(' ');
        disp(' Enter speed of sound in material (ft/sec): ');
        c(i)=input(' ');
        c(i)=c(i)*12;
%
        disp(' ');
        disp(' Select method: 1=mass density  2=absorption coefficient ');
        imd=input(' ');
        disp(' ');
%
        if(imd==1)
            disp(' Enter mass density (lbm/in^3): ');
            rho=input(' ');
            rho=rho/386;
            R(i)=rho*c(i);    
        else
           disp(' Enter absorption coefficient: ');
           alpha=input(' ');
           aa=1;
           bb=2-4/alpha;
           cc=1;
           R(i)=(-bb+sqrt(bb^2-4*aa*cc))/(2*aa);  % Crocker page 699, equation (8)
           R(i)=R(i)*415/22615;
        end    
%
    end    
%
end
%
disp(' ');
disp('      Acoustic      Speed ');  
disp('      Impedance     Sound    Thickness');  
disp('      (Rayls)      (ft/sec)   (inch)  ');  
%
for(i=1:5)
    out1=sprintf(' %d \t %8.4g \t %8.4g \t %6.3g',i,R(i)*22615,c(i)/12,L(i));
    disp(out1);
end
disp(' ');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
j=sqrt(-1);
%
for(i=1:5000)
%
    omega=2*pi*i;
%    
    k2=omega/c(2);
    k3=omega/c(3);
    k4=omega/c(4);   
%
    M=zeros(9,9);
%    
    M(1,1)=1;
%    
    M(2,1)=1;
    M(2,2)=1;
    M(2,3)=-1;
    M(2,4)=-1;
%    
    M(3,3)=exp(-j*k2*L(2));
    M(3,4)=exp(+j*k2*L(2));
    M(3,5)=-1;
    M(3,6)=-1;
%    
    M(4,5)=exp(-j*k3*L(3));
    M(4,6)=exp(+j*k3*L(3));
    M(4,7)=-1;
    M(4,8)=-1;
%    
    M(5,7)=exp(-j*k4*L(4));
    M(5,8)=exp(+j*k4*L(4));
    M(5,9)=-1;
%    
    M(6,1)=-R(1);
    M(6,2)= R(1);    
    M(6,3)= R(2);
    M(6,4)=-R(2);
%    
    M(7,3)=-R(2)*exp(-j*k2*L(2));
    M(7,4)= R(2)*exp(+j*k2*L(2));
    M(7,5)= R(3);
    M(7,6)=-R(3);    
%
    M(8,5)=-R(3)*exp(-j*k3*L(3));
    M(8,6)= R(3)*exp(+j*k3*L(3));
    M(8,7)= R(4);
    M(8,8)=-R(4);
%
    M(9,7)=-R(4)*exp(-j*k4*L(4));
    M(9,8)= R(4)*exp(+j*k4*L(4));
    M(9,9)= R(5);     
%    
    B=zeros(9,1);
    B(1,1)=1.;
%
    A=M\B;
%
    TL(i,1)=i;
    f(i)=i;
    TL(i,2)=abs(A(9));
    alpha_t(i) = (R(1)/R(5))*(TL(i,2))^2;   
    aaa(i)=-10.*log10(alpha_t(i)); 
end
%
transmission_loss=[f' aaa'];
%
disp(' ');
disp(' Transmission loss variable is named:  transmission_loss ');
%
plot(f,aaa);
%
axis([1,10000,0,60]);
set(gca,'ytick',[0 5 10 15 20 25 30 35 40 45 50 55 60])
set(gca,'YTickLabel',{'0';'5';'10';'15';'20';'25';'30';'35';'40';'45';'50';'55';'60'})
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','lin');
grid on;
xlabel('Frequency (Hz)');
ylabel('Transmission Loss (dB)');
title('Transmission Loss through Five Media');
%