disp(' ');
disp(' single_partition.m   ver 1.3  November 19, 2007 ');
disp(' by Tom Irvine   Email: tomirvine@aol.com ');
disp(' ');
disp(' This script calculates the transmission loss through a single partition.');
disp(' Reference: Beranek & Ver, Noise & Vibration Control Engineering, page 297.');
%
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)
    plateau_height=29;
    plateau_width=11;    
end
if(imat==2)
    E=10.0e+06;      % elastic modulus (psi)
    rho=0.00015024;  % mass/volume (lbf sec^2/in^4)
    plateau_height=29;
    plateau_width=11;
%
    disp(' ');
    disp(' Suggest plateau height between 19 dB & 56 dB, where 19 dB is conservative.' );
    disp(' ');
    disp(' Enter plateau height (dB) ');
    plateau_height=input(' ');  
    disp(' ')
    disp(' Suggest plateau width frequency ratio between 4 & 11.' );
    disp(' ');
    disp(' Enter plateau width ');
    plateau_width=input(' ');   
%
end
if(imat >2)
    disp(' ');
    disp(' Enter Elastic Modulus (lbf/in^2)' );
    E=input(' ');
    disp(' ');
    disp(' Enter mass/volume (lbm/in^3)' );
    rho=input(' ');
    rho=rho/386;
    disp(' ');
    disp(' Suggest plateau height between 19 dB & 56 dB, where 19 dB is conservative.' );
    disp(' ');
    disp(' Enter plateau height (dB) ');
    plateau_height=input(' ');  
    disp(' ')
    disp(' Suggest plateau width frequency ratio between 4 & 11.' );
    disp(' ');
    disp(' Enter plateau width ');
    plateau_width=input(' ');   
end
v=0.33;         % Poisson ratio
%
%
disp(' ');
disp(' Enter thickness (in)' );
h=input(' ');
%
air_rho=1.2103e-007;   % lbf sec^2/in^4  
air_c=13500;           % in/sec
ci_air=air_rho*air_c;  % lbf sec/in^3
%
        CL=sqrt(E/rho);
critical_f=(air_c^2/(2*pi*h))*sqrt(12*rho*(1-v^2)/E);
%
disp(' ');
out1=sprintf(' critical_f = %9.4g Hz',critical_f);
disp(out1);
%
surface_mass=386*rho*h*1550/2.205;
smcf=surface_mass*critical_f;
%
disp(' ');
out1=sprintf(' surface_mass*critical_f = %9.4g Hz kg/m^2 ',smcf);
disp(out1);
out1=sprintf('                         = %9.4g Hz lbm/ft^2 ',smcf*2.205/10.764);
disp(out1);
%
x=log10(smcf);
m=20;   % slope = 20 dB/decade
x1=3;   % sample coordinate at (3, 18) or (1000 Hz kg/m^2, 18 dB)
y1=18;
b=y1-m*x1;
y=m*x+b;
Rzero=y;
%
%Rfield=Rzero-5
x=(plateau_height+5-b)/m;
flow=10^x/surface_mass;
%
disp(' ');
out1=sprintf(' start of plateau frequency (Hz) = %9.4g Hz  ',flow);
disp(out1);
%
fup=flow*plateau_width;
disp(' ');
out1=sprintf('   end of plateau frequency (Hz) = %9.4g Hz  ',fup);
disp(out1);
%
	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;
%
clear ff;
clear TL;
%
k=1;
for(i=7:imax)
    ff(k)=fc(i);
    TL(k)=0.;
    if(fc(i)<flow)
        x=fc(i)*surface_mass;
        x=log10(x);
        y=m*x+b;
        if(y>=5)
            TL(k)=y-5;
            ff(k)=fc(i);
        end
    end
    if(fc(i)>=flow & fc(i)<=fup)
        TL(k)=plateau_height;
        ff(k)=fc(i);  
    end
    if(fc(i)>fup)
        TL(k)=plateau_height+10*log(fc(i)/fup)/log(2);
        ff(k)=fc(i);
    end
    k=k+1;
end
kmax=k-1;
%
clear limit;
limit=0.;
ilimit=0;
disp(' ')
disp(' Limit attenuation? 1=yes 2=no ');
ilimit=input(' ');
%
if(ilimit==1)
    disp(' ')
    disp(' Enter limit (dB) ');
    limit=input(' ');
    for(i=7:kmax)
        if(TL(i)>limit)
            TL(i)=limit;
        end
    end
end
%
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)');
title('Transmission Loss')
%
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
%
disp(' ')
disp(' Matlab filename: transmission_loss ');
%