disp(' ');
disp(' single_partition_plane_wave.m   ver 1.3  September 23, 2009 ');
disp(' by Tom Irvine   Email: tomirvine@aol.com ');
disp(' ');
disp(' This script calculates the transmission loss through a single partition.');
disp(' Air is the medium on either side.');
disp(' Reference: W. Seto, Acoustics, equation (12)');
%
disp(' ');
disp(' Select Material:  1=aluminum  2=graphite epoxy  3=other' );
imat=input(' ');
%
if(imat==1)
    E=9.9e+06;      % elastic modulus (psi)
    rho=0.0002539;  % mass/volume (lbf sec^2/in^4)
    CL=sqrt(E/rho);   
    R2=rho*CL;    
end
if(imat==2)
    E=10.0e+06;      % elastic modulus (psi)
    rho=0.00015024;  % mass/volume (lbf sec^2/in^4)
    CL=sqrt(E/rho);   
    R2=rho*CL;       
end
if(imat >2)
%
        disp(' ');
        disp(' Select mass unit: 1= lbm/in^3  2= lbm/ft^3 ');
        imu=input(' ');
%        
        if(imu==1)
            disp(' Enter mass/volume (lbm/in^3)' );
            rho=input(' ');
        else
            disp(' Enter mass/volume (lbm/ft^3)' );            
            rho=input(' ');
            rho=rho/12^3;
        end  
%        
        rho=rho/386;  
        disp(' ');
        out1=sprintf(' Select method: \n 1= Elastic Modulus  2= speed of sound  3= absorption');
        disp(out1);
        imd=input(' ');
        disp(' ');      
%
        if(imd==1)
            disp(' ');
            disp(' Enter Elastic Modulus (lbf/in^2)' );
            E=input(' ');
            disp(' ');
            CL=sqrt(E/rho);   
            R2=rho*CL;    
        end
        if(imd==2)
            disp(' Enter speed of sound in material (ft/sec): ');
            CL=input(' ');
            CL=CL*12;
            R2=R2*415/22615;
        end  
        if(imd==3)
           disp(' Enter absorption coefficient: ');
           alpha=input(' ');
           aa=1;
           bb=2-4/alpha;
           cc=1;
           R2=(-bb+sqrt(bb^2-4*aa*cc))/(2*aa);  % Crocker page 699, equation (8)
           R2=R2*415/22615;
           disp(' ');
           disp(' Enter speed of sound in material (ft/sec): ');
           CL=input(' ');
           CL=CL*12;  
        end
%
end
%
disp(' ');
disp(' Enter thickness (in)' );
h=input(' ');
%
R1=415/22615;  % Rayls converted to slugs/(in^2 sec)
%
disp(' ');
out1=sprintf('\n CL = %9.5g ft/sec ',CL/12);
disp(out1);
%
out1=sprintf('\n R1 = %9.5g Rayls  \t R2 = %9.5g Rayls ',R1*22615,R2*22615);
disp(out1);
%
out1=sprintf('\n R2/R1 = %9.5g ',R2/R1);
disp(out1);
%
disp(' ');
%
	fc(1)=2.5;
	fc(2)=3.15;
	fc(3)=4.;;
	fc(4)=5.;
	fc(5)=6.3;
	fc(6)=8.;
	fc(7)=10.;
	fc(8)=12.5;
	fc(9)=16.;
	fc(10)=20.;
	fc(11)=25.;
	fc(12)=31.5;
	fc(13)=40.;
	fc(14)=50.;
	fc(15)=63.;
	fc(16)=80.;
	fc(17)=100.;
	fc(18)=125.;
	fc(19)=160.;
	fc(20)=200.;
	fc(21)=250.;
	fc(22)=315.;
	fc(23)=400.;
	fc(24)=500.;
	fc(25)=630.;
	fc(26)=800.;
	fc(27)=1000.;
	fc(28)=1250.;
	fc(29)=1600.;
	fc(30)=2000.;
	fc(31)=2500.;
	fc(32)=3150.;
	fc(33)=4000.;
	fc(34)=5000.;
	fc(35)=6300.;
	fc(36)=8000.;
	fc(37)=10000.;
	fc(38)=12500.;
	fc(39)=16000.;
	fc(40)=20000.;
    imax=40;
%
disp(' ');
disp(' Enter frequency spacing ');
disp(' 1=one-third octave   2= constant one Hz');
ifs=input(' ');
%
clear ff;
clear TL;
%
k=1;
L=h;
%
if(ifs==1)
%
    for(i=7:imax)
        k2=2*pi*fc(i)/CL;
        ff(k)=fc(i);
        den=4*((cos(k2*L))^2) + ((R2/R1)^2)*((sin(k2*L))^2);
        tau=4/den;
        TL(k)=10*log10(1/tau);
        k=k+1;
    end
%
else
%
    for(i=1:20000)
        k2=2*pi*i/CL;
        ff(k)=i;
        den=4*((cos(k2*L))^2) + ((R2/R1)^2)*((sin(k2*L))^2);
        tau=4/den;
        TL(k)=10*log10(1/tau);
        k=k+1;
    end
%
end
%    
kmax=k-1;
%
clear limit;
%
clear transmission_loss;
transmission_loss=[ff' TL'];
%
figure(1);
plot(ff,TL);
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','lin','XminorTick','off','YminorTick','off');
grid('on');
xlabel('Frequency (Hz)');
ylabel('Transmission Loss (dB)');
out1=sprintf('Transmission Loss   Thickness=%8.4g inch',h);
title(out1);
%
if(ifs==1)
%
disp(' ');
disp(' Perform force response analysis?  1=yes  2=no ');
nop=input(' ');
if(nop==1)
disp('  Forced Response Analysis ');
disp(' ');
disp(' Enter one-third octave SPL(dB) input file.')
disp(' Select file input method ');
disp('   1=external ASCII file ');
disp('   2=file preloaded into Matlab ');
file_choice = input('');
%
clear SPL;
%
if(file_choice==1)
    [filename, pathname] = uigetfile('*.*');
    filename = fullfile(pathname, filename);
%      
    fid = fopen(filename,'r');
    SPL = fscanf(fid,'%g %g',[2 inf]);
    SPL=SPL';
else
    SPL = input(' Enter the matrix name:  ');
end
%
clear ft;
clear response;
clear rms;
clear pressure_rms;
%
k=1;
n=size(SPL);
%
reference=2.9e-09;
rms=0;
%
mm=size(TL');
%
disp(' ');
disp(' Freq(Hz)  SPL(dB)  TL(dB)  Response(dB) ');
%
for(i=1:mm(1))
    for(j=1:n(1))
        if(SPL(j,1)==ff(i))
            ft(k)=ff(i);
            response(k)=SPL(j,2)-TL(i);
%            
            out1=sprintf('  %5.1f   %5.1f   %5.1f   %5.1f ',ft(k),SPL(j,2),TL(i),response(k));
            disp(out1);
%            
            pressure_rms=reference*(10.^(response(k)/20.) );
            rms=rms+pressure_rms^2;
            k=k+1;
        end
    end
end    
%
rms=sqrt(rms);
dbrms=20*log10(rms/reference);
%
clear response_SPL;
response_SPL=[ft' response'];
%
figure(2);
plot(ft,response);
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','lin','XminorTick','off','YminorTick','off');
grid('on');
xlabel('Center Frequency (Hz)');
ylabel('1/3 Octave SPL (dB) ');
out1=sprintf('SOUND PRESSURE LEVEL   OASPL = %8.3g dB  (Ref 20 micro Pa)',dbrms);
title(out1);
%
disp(' ');
disp(' Write response SPL to ASCII file? ')
pchoice=input(' 1=yes 2=no ');
if(pchoice == 1) 
%
      [writefname, writepname] = uiputfile('*','Save PSD data as');
	  writepfname = fullfile(writepname, writefname);
	  writedata = [ ft',  response' ];
	  fid = fopen(writepfname,'w');
	  fprintf(fid,'  %g \t %g \n',writedata');
	  fclose(fid);
end  
end
%
if(nop==1)
    disp(' ')
    disp(' Matlab filename: response_SPL ');
end
end
%
disp(' ')
disp(' Matlab filename: transmission_loss ');
%