disp(' ')
disp(' sdof_ran.m   ver 2.3   November 9, 2013')
disp(' by Tom Irvine   Email: tom@vibrationdata.com')
disp(' ')
disp(' This program calculates the response of a single-degree-of-freedom system')
disp(' to an acceleration power spectral density')
disp(' ')
disp(' The power spectral density may vary arbitrarily with frequency. ')
disp(' The input file must have two columns: frequency(Hz) & accel(G^2/Hz)')
disp(' ')
%
close all;
%
clear f;
clear a;
clear a_vrs;
clear trans;
clear s;
clear damp;
clear THF;
%
fig_num=1;
%
disp(' Select file input method       ');
disp('   1=external ASCII file        ');
disp('   2=file preloaded into Matlab ');
disp('   3=library function           ');
file_choice = input('');
%
if(file_choice==1)
    [filename, pathname] = uigetfile('*.*');
    filename = fullfile(pathname, filename);
    fid = fopen(filename,'r');
    THF = fscanf(fid,'%g %g',[2 inf]);
    THF=THF';
end
if(file_choice==2)
    THF = input(' Enter the matrix name:  ');
end
if(file_choice==1 || file_choice==2)
    size(THF);
    f=THF(:,1);
    a=THF(:,2);
end
%
if(file_choice==3)
%
   [f,a]=PSD_library();
%
end
%
MAX = 5000;
tpi=2.*pi;
%
disp(' ');
natural_frequency=input(' Enter the natural frequency (Hz)  ');
%
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
%
[s,grms]=calculate_PSD_slopes(f,a);
out5 = sprintf('\n Input PSD Overall Level = %10.3f ',grms);
disp(out5);
%
df=1.;
%
[fi,ai]=interpolate_PSD(f,a,s,df);
%
[a_vrs,rd_vrs,trans,opsd]=sdof_ran_engine(fi,ai,damp,df,natural_frequency);
%
Q=1/(2*damp);
%
disp(' ')
disp(' Enter duration (sec) ')
dur=input(' ');
%
c=sqrt(2*log(natural_frequency*dur));
ms=(c+0.5772/c);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
clear length;
n=length(fi);
%
m2=0;
m4=0;
for i=1:n
    m2=m2+opsd(i)*fi(i)^2;
    m4=m4+opsd(i)*fi(i)^4;
end
%
m2=(m2*df);
m4=(m4*df);
%
EP=sqrt(m4/m2);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ')
disp(' SDOF Acceleration Response ')
disp(' ')
out1=sprintf('     = %8.3g GRMS',a_vrs);
out2=sprintf('     = %8.3g G 3-sigma',3.*a_vrs);
disp(out1);
disp(out2);
%
if(natural_frequency>=min(f) && natural_frequency<=max(f))
    out3=sprintf('     = %8.3g G %4.3g-sigma (maximum expected)',ms*a_vrs,ms);
    disp(out3);
end
%
out1=sprintf('\n Number of expected acceleration peaks = %8.4g \n',EP*dur);
disp(out1);
%
omegan=tpi*natural_frequency;
%
disp(' ')
disp(' SDOF Pseudo Velocity Response ')
disp(' ')
out1=sprintf('     = %8.3g inch/sec RMS',omegan*rd_vrs);
out2=sprintf('     = %8.3g inch/sec 3-sigma',omegan*3.*rd_vrs);
disp(out1);
disp(out2);
%
if(natural_frequency>=min(f) && natural_frequency<=max(f))
    out3=sprintf('     = %8.3g inch/sec %4.3g-sigma (maximum expected)',omegan*ms*rd_vrs,ms);
    disp(out3);
end
disp(' ')
out1=sprintf('     = %8.3g mm/sec RMS',omegan*rd_vrs*25.4);
out2=sprintf('     = %8.3g mm/sec 3-sigma',omegan*3.*rd_vrs*25.4);
disp(out1);
disp(out2);
%
disp(' ')
disp(' SDOF Relative Displacement Response ')
disp(' ')
out1=sprintf('     = %8.3g inch RMS',rd_vrs);
out2=sprintf('     = %8.3g inch 3-sigma',3.*rd_vrs);
disp(out1);
disp(out2);
%
if(natural_frequency>=min(f) && natural_frequency<=max(f))
    out3=sprintf('     = %8.3g inch %4.3g-sigma (maximum expected)',ms*rd_vrs,ms);
    disp(out3);
end
disp(' ')
out1=sprintf('     = %8.3g mm RMS',rd_vrs*25.4);
out2=sprintf('     = %8.3g mm 3-sigma',3.*rd_vrs*25.4);
disp(out1);
disp(out2);
%
if(natural_frequency>=min(f) && natural_frequency<=max(f))
    out3=sprintf('     = %8.3g mm %4.3g-sigma (maximum expected)',ms*rd_vrs*25.4,ms);
    disp(out3);
end
%
fi=fix_size(fi);
opsd=fix_size(opsd);
a=fix_size(ai);
trans=fix_size(trans);
%
ppp=[fi opsd];
qqq=[fi a];
%
clear response_psd;
clear power_trans;
%
response_psd=ppp;
power_trans=[fi trans];
%
leg_a=sprintf('Response %5.3g GRMS',a_vrs);
leg_b=sprintf('Input %5.3g GRMS',grms);
%              
x_label=sprintf(' Frequency (Hz)');
y_label=sprintf(' Accel (G^2/Hz)'); 
t_string = sprintf(' Power Spectral Density   fn=%g Hz  Q=%g  ',natural_frequency,Q);   
%
[fig_num]=plot_PSD_two_sdof_ran(fig_num,x_label,y_label,t_string,ppp,qqq,leg_a,leg_b); 
%
ppp=[fi trans];  
x_label=sprintf(' Frequency (Hz)');
y_label=sprintf(' Trans (G^2/G^2)'); 
t_string = sprintf(' Power Transmissibility   fn=%g Hz  Q=%g',natural_frequency,Q);   
[fig_num]=plot_PSD(fig_num,x_label,y_label,t_string,ppp);
%
disp(' ');
disp(' Output files: ');
disp('   response_psd ');
disp('   power_trans ');
%
disp(' ');
disp(' Calculate Dirlik rainflow fatigue cycles?  1=yes 2=no ');
idc=input(' ');
if(idc==1)
    [range,amp,fig_num]=Dirlik_function(response_psd,dur,fig_num);
end