disp(' ') disp(' sdof_ran.m ver 1.5 July 21, 2011') disp(' by Tom Irvine Email: tomirvine@aol.com') disp(' ') disp(' This program calculates the response of a single-degree-of-freeom 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(' ') % clear f; clear a; clear s; clear THF; % 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) disp(' Select input PSD '); disp(' 1=MIL-STD-1540B QTP '); disp(' 2=MIL-STD-1540B ATP '); disp(' 3=NAVMAT P9495 '); disp(' 4=NASA CEQATR AVT '); disp(' 5=NASA-STD-7001 '); ilibrary=input(' '); % if(ilibrary==1 | ilibrary==2 ) % MIL-STD-1540B ATP f(1)=20.; a(1)=0.0053; f(2)=150.; a(2)=0.04; f(3)=600.; a(3)=0.04; f(4)=2000.; a(4)=0.0036; n=4; Lflag=1; end if(ilibrary==1) % MIL-STD-1540B QTP a(1)=a(1)*4.; a(2)=a(2)*4.; a(3)=a(3)*4.; a(4)=a(4)*4.; end if(ilibrary==3) % NAVMAT P9495 f(1)=20.; a(1)=0.01; f(2)=80.; a(2)=0.04; f(3)=350.; a(3)=0.04; f(4)=2000.; a(4)=0.007; end if(ilibrary==4) % NASA CEQATR AVT f(1)=20.; a(1)=0.01; f(2)=80.; a(2)=0.04; f(3)=800.; a(3)=0.04; f(4)=2000.; a(4)=0.01; end if(ilibrary==5) % NASA-STD-7001 f(1)=20.; a(1)=0.01; f(2)=80.; a(2)=0.04; f(3)=500.; a(3)=0.04; f(4)=2000.; a(4)=0.01; end % 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 % n=length(f); % [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(' '); ms=sqrt(2*log(natural_frequency*dur)); % 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 % 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(' '); %%% disp(' Plot the output PSD? ') %%% choice = input(' 1=yes 2=no '); choice=1; % if(choice == 1) % figure(1); plot(fi,opsd,f,a,'-.'); out1=sprintf('Response %5.3g GRMS',a_vrs); out2=sprintf('Input %5.3g GRMS',grms); % disp(out1) % disp(out2) legend(out1,out2,2); xlabel(' Frequency (Hz)'); ylabel(' Accel (G^2/Hz)'); at = sprintf(' Power Spectral Density fn=%g Hz Q=%g ',natural_frequency,Q); title(at); set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log'); grid; % maxa=max(a); mina=min(a); % if(max(opsd)>maxa) maxa=max(opsd); end % maxf=max(f); minf=min(f); % if min(f) >= 0.1 minf=0.1; end if min(f) >= 1 minf=1; end if min(f) >= 10 minf=10; end if min(f) >= 20 minf=20; end if min(f) >= 100 minf=100; end % ttt=1; if(max(f)<=10000 & minf >=10) maxf=10000; set(gca,'xtick',[10 100 1000 10000]); ttt=1; end if(max(f)<=5000 & minf >=10) maxf=5000; set(gca,'xtick',[10 100 1000 5000]); ttt=2; end if(max(f)<=2000 & minf >=10) maxf=2000; set(gca,'xtick',[10 100 1000 2000]); ttt=3; end if(max(f)<=1000 & minf >=10) maxf=1000; set(gca,'xtick',[10 100 1000]); ttt=4; end if(max(f)<=500 & minf >=10) maxf=1000; set(gca,'xtick',[10 100 1000]); ttt=5; end % if(max(f)<=10000 & minf >=20) maxf=10000; set(gca,'xtick',[20 100 1000 10000]); ttt=6; end % if(max(f)<=2000 & minf >=20) maxf=2000; set(gca,'xtick',[20 100 1000 2000]); ttt=7; end % clear aaa; % ymax= 10^(round(log10(maxa)+0.8)); ymin= 10^(round(log10(mina)-0.6)); for(j=20:-1:-20) aaa=(10^j); if(ymaxaaa) break; else end end ymin=aaa; axis([minf,maxf,ymin,ymax]); % end % disp(' '); disp(' Plot the power transmissibility function? ') choice = input(' 1=yes 2=no '); % if(choice == 1) % figure(2); plot(fi,trans); xlabel(' Frequency (Hz)'); ylabel(' Trans (G^2/G^2)'); at = sprintf(' Power Transmissibility fn=%g Hz Q=%g',natural_frequency,Q); title(at); set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log'); grid; % maxa=max(trans); mina=min(trans); % clear aaa; % ymax= 10^(round(log10(maxa)+0.8)); ymin= 10^(round(log10(mina)-0.6)); % for(j=20:-1:-20) aaa=(10^j); if(ymaxaaa) break; else end end ymin=aaa; % if(ttt==1) set(gca,'xtick',[10 100 1000 10000]) end if(ttt==2) set(gca,'xtick',[10 100 1000 5000]); end if(ttt==3) set(gca,'xtick',[10 100 1000 2000]); end if(ttt==4) set(gca,'xtick',[10 100 1000]); end if(ttt==5) set(gca,'xtick',[10 100 1000]); end if(ttt==6) set(gca,'xtick',[20 100 1000 10000]); end if(ttt==7) set(gca,'xtick',[20 100 1000 2000]); end % axis([minf,maxf,ymin,ymax]); % end