disp(' ')
disp(' sdof_spl.m  version 1.4, December 31, 2013 ')
disp(' by Tom Irvine ')
disp(' Email:  tom@vibrationdata.com ')
disp(' ')
%
disp(' This program calculates the response of an SDOF system to')
disp(' an applied acoustic pressure field which is uniform and ');
disp(' fully correlated over a surface area. ');
disp(' ')
%
clear db;
clear rms;
clear sum;
clear ms;
clear oadb;
clear THM;
clear fc;
clear spl;
clear length;
clear area;
clear FT_psd;
clear apsd;
clear vpsd;
clear dpsd;
%
close all
%
fig_num=1;
%
[THM,oadb,delta,fig_num]=enter_spl(fig_num);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
[fc,psd,applied_pressure_psd,tstring,fig_num,iu]=...
                                     convert_spl_to_psd(THM,delta,fig_num);
%     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
if(iu==1)
   area=input(' Enter surface area (in^2) ');
else
   area=input(' Enter surface area (m^2) '); 
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
clear THM;
clear F;
clear a;
clear fi;
clear ai;
clear length;
%
F=applied_pressure_psd(:,1);
a=applied_pressure_psd(:,2)*area^2;
%
tpi=2.*pi;
%
nnn=length(F);
%
f1=F(1);
f2=F(nnn);
%
[m0,frms]=cm0(a,F,nnn);
%
figure(fig_num);
fig_num=fig_num+1;
plot(F,a);
grid on;
xlabel('Frequency (Hz)');
if(iu==1)
    ylabel('Force (lbf^2/Hz)');
else
    ylabel('Force (Pa^2/Hz)');    
end    
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log','XminorTick','on','YminorTick','on');
%
if(iu==1)
    out1=sprintf('Applied Force PSD  Overall Level = %10.3f lbf RMS',frms);
    out5 = sprintf('\n Input PSD Overall Level = %10.3f lbf RMS',frms);
else
    out1=sprintf('Applied Force PSD  Overall Level = %10.3f N RMS',frms);
    out5 = sprintf('\n Input PSD Overall Level = %10.3f N RMS',frms);    
end    
%
title(out1);
disp(out5);
%
ocn=1./24.;
%
fi(1)=f1;
%
j=2;
while(1)
    fi(j)=fi(j-1)*2^(ocn);
    if(fi(j)>f2)
        fi(j)=f2;
        break;
    end
    j=j+1;
end
%
n_ref=length(fi);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
[ai]=interp_psd_oct(F,a,fi);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
[fn,mass,stiffness,omegan,omegan2]=enter_fn_mass_stiffness(iu,tpi);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
idamp=input(' Enter damping format:  1= damping ratio   2= Q  ');	
%
disp(' ')
if(idamp==1)
    damp=input(' Enter damping ratio (typically 0.05) ');
else
    Q=input(' Enter the amplification factor (typically Q=10) ');
    damp=1/(2*Q);
end
%
disp('  Also calculate relative damage index?  1=yes 2=no  ');
ird=input(' ');
disp(' ');
%
if(ird==1)
    bex=input(' Enter fatigue exponent  ');
    duration=input(' Enter duration (sec)  ');
    disp(' ');
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
omega=zeros(n_ref,1);
omega2=zeros(n_ref,1);
%
for j=1:n_ref   
    omega(j)=tpi*fi(j);
    omega2(j)=omega(j)^2;
end
%
dpsd=zeros(n_ref,1);
vpsd=zeros(n_ref,1);
apsd=zeros(n_ref,1);
FT_psd=zeros(n_ref,1);
%
%      Absolute Displacement
%
for j=1:n_ref
%        
	num=( (omegan2/stiffness)^2 );
	den=( (omegan2-omega2(j))^2) + ( (2*damp*omega(j)*omegan)^2);
	t= num/den;
    dpsd(j)=( t*ai(j) );
%
end
[m0,drms]=cm0(dpsd,fi,n_ref);
%
%      Velocity
%
for j=1:n_ref
%     
    num=( (omega(j)*omegan2/stiffness)^2 );
    den=( (omegan2-omega2(j))^2) + ( (2*damp*omega(j)*omegan)^2);
    t= num/den;
    vpsd(j)=( t*ai(j) );
%
end
[m0,vrms]=cm0(vpsd,fi,n_ref);
%
%      Acceleration
%
for j=1:n_ref
%
    num=( (omega2(j)*omegan2/stiffness)^2 );
    den=( (omegan2-omega2(j))^2) + ( (2*damp*omega(j)*omegan)^2);
    t= num/den;
    apsd(j)=( t*ai(j) );
%
end
%
if(iu==1)
    apsd=apsd/(386^2);
else
    apsd=apsd/(9.81^2);        
end 
%
[a_m0,grms]=cm0(apsd,fi,n_ref);
%
%      Transmitted Force
%
for j=1:n_ref
%
    rho=fi(j)/fn;
    num= 1.+(2*damp*rho)^2;
    den= (1-rho^2)^2 + (2*damp*rho)^2;
    t= num/den;
    FT_psd(j)=( t*ai(j) );
%
end
[FT_m0,FT_rms]=cm0(FT_psd,fi,n_ref);
%
fi=fix_size(fi);
apsd=fix_size(apsd);
%
figure(fig_num);
fig_num=fig_num+1;
plot(fi,apsd);
grid on;
t_string=sprintf(' Acceleration PSD   %8.4g GRMS Overall   ',grms);
title(t_string);
xlabel('Frequency (Hz)');
ylabel('Accel (G^2/Hz)');
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log','XminorTick','on','YminorTick','on');
%
%
figure(fig_num);
fig_num=fig_num+1;
plot(fi,FT_psd);
grid on;
xlabel('Frequency (Hz)');
if(iu==1)
    ylabel('Force (lbf^2/Hz)');
    t_string=sprintf(' Transmitted Force PSD   %8.4g  lbf RMS Overall   ',FT_rms);
else
    ylabel('Force (N^2/Hz)');
    t_string=sprintf(' Transmitted Force PSD   %8.4g  N RMS Overall   ',FT_rms);
end    
title(t_string);
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log','XminorTick','on','YminorTick','on');
%
%
if(ird==1)
    [a_m1]=cm1(apsd,fi,n_ref);
    [a_m2]=cm2(apsd,fi,n_ref);
    [a_m4]=cm4(apsd,fi,n_ref);
%
    a_EP=sqrt(a_m4/a_m2);
%
    [a_damage]=Dirlik_basic(duration,bex,a_m0,a_m1,a_m2,a_m4,a_EP,grms);
%
    [FT_m1]=cm1(FT_psd,fi,n_ref);
    [FT_m2]=cm2(FT_psd,fi,n_ref);
    [FT_m4]=cm4(FT_psd,fi,n_ref);
%
    FT_EP=sqrt(FT_m4/FT_m2);
%
    [ft_damage]=Dirlik_basic(duration,bex,FT_m0,FT_m1,FT_m2,FT_m4,FT_EP,FT_rms);
%
end    
%
% *************************************************************************
% 
disp(' ');
disp(' Output Overall Levels ');
disp(' ');
%	    
if(iu==1)
    out1=sprintf('  Displacement = %12.4g inches       RMS',drms);
    out2=sprintf('      Velocity = %12.4g inches/sec   RMS',vrms);
    out3=sprintf('  Acceleration = %12.4g G            RMS \n',grms);
    out4=sprintf('  Transmitted Force    = %12.4g lbf RMS',FT_rms);
    out5=sprintf('  Transmitted Pressure = %12.4g psi RMS',FT_rms/area);
else
    out1=sprintf('  Displacement = %12.4g mm     RMS',drms*1000);
    out2=sprintf('      Velocity = %12.4g m/sec  RMS',vrms);
    out3=sprintf('  Acceleration = %12.4g G      RMS \n',grms);
    out4=sprintf('  Transmitted Force    = %12.4g N RMS ',FT_rms);
    out5=sprintf('  Transmitted Pressure = %12.4g Pa RMS',FT_rms/area);
end
%
disp(out1);
disp(out2);
disp(out3);
disp(out4);
disp(out5);
%
%
if(ird==1)
%
    disp(' ');
    disp(' Relative Fatigue Damage Indices (peak-valley)/2 ');
    disp(' with respect to:');
%
    out9=sprintf('\n      Acceleration    = %12.4g '  ,a_damage);
    disp(out9);            
%
    out11=sprintf(' Transmitted Force    = %12.4g '  ,ft_damage);
    disp(out11); 
%
    out13=sprintf(' Transmitted Pressure = %12.4g '  ,ft_damage/(area^bex));
    disp(out13); 
%
end