Shock Response Spectrum

by

 

08 Apr 2005 (Updated )

This program calculates the SRS of a time history.

srs.m
disp(' ')
disp(' srs.m   ver 2.0   July 3, 2006')
disp(' by Tom Irvine   Email: tomirvine@aol.com')
disp(' ')
disp(' This program calculates the shock response spectrum')
disp(' of an acceleration time history, which is pre-loaded into Matlab.')
disp(' The time history must have two columns: time(sec) & acceleration')
disp(' ')
%
clear t;
clear y;
clear yy;
clear n;
clear fn;
clear a1;
clear a2
clear b1;
clear b2;
clear jnum;
clear THM;
clear resp;
clear x_pos;
clear x_neg;
%
iunit=input(' Enter acceleration unit:   1= G   2= m/sec^2  ');	
%
disp(' ')
disp(' Select file input method ');
disp('   1=external ASCII file ');
disp('   2=file preloaded into Matlab ');
file_choice = input('');
%
if(file_choice==1)
    [filename, pathname] = uigetfile('*.*');
    filename = fullfile(pathname, filename);
%      
    fid = fopen(filename,'r');
    THM = fscanf(fid,'%g %g',[2 inf]);
    THM=THM';
else
    THM = input(' Enter the matrix name:  ');
end
%
t=double(THM(:,1));
y=double(THM(:,2));
%
tmx=max(t);
tmi=min(t);
n = length(y);
%
out1 = sprintf('\n  %d samples \n',n);
disp(out1)
%
dt=(tmx-tmi)/(n-1);
sr=1./dt;
%
out1 = sprintf(' SR  = %g samples/sec    dt = %g sec \n',sr,dt);
disp(out1)
%
fn(1)=input(' Enter the starting frequency (Hz)  ');
if fn(1)>sr/30.
    fn(1)=sr/30.;
end
%
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(' ')
disp(' Select algorithm: ')
disp(' 1=Kelly-Richman  2=Smallwood ');
ialgorithm=input(' ');
%
tmax=(tmx-tmi) + 1./fn(1);
limit = round( tmax/dt );
n=limit;
yy=zeros(1,limit);
for i=1:length(y)
        yy(i)=y(i);
end    
%
disp(' ')
disp(' Calculating response..... ')
%
%  SRS engine
%
for j=1:1000
%
    omega=2.*pi*fn(j);
    omegad=omega*sqrt(1.-(damp^2));
    cosd=cos(omegad*dt);
    sind=sin(omegad*dt);
    domegadt=damp*omega*dt;
%
    if(ialgorithm==1)
        a1(j)=2.*exp(-domegadt)*cosd;
        a2(j)=-exp(-2.*domegadt);
        b1(j)=2.*domegadt;
        b2(j)=omega*dt*exp(-domegadt);
        b2(j)=b2(j)*( (omega/omegad)*(1.-2.*(damp^2))*sind -2.*damp*cosd );
        b3(j)=0;
%
    else
	    E=exp(-damp*omega*dt);
		K=omegad*dt;
		C=E*cos(K);
		S=E*sin(K);
		Sp=S/K;
%
	    a1(j)=2*C;
		a2(j)=-E^2;
		b1(j)=1.-Sp;
		b2(j)=2.*(Sp-C);
		b3(j)=E^2-Sp;
    end
    forward=[ b1(j),  b2(j),  b3(j) ];    
    back   =[     1, -a1(j), -a2(j) ];    
%    
    resp=filter(forward,back,yy);
%
    x_pos(j)= max(resp);
    x_neg(j)= min(resp);
%   
    jnum=j; 
    if  fn(j) > sr/8.
        break
    end
    fn(j+1)=fn(1)*(2. ^ (j*(1./12.)));    
end
%
%  Output options
%
disp(' ')
disp(' Select output option ');
choice=input(' 1=plot only   2=plot & output text file  ' );
disp(' ')
%
if choice == 2 
%%
      [writefname, writepname] = uiputfile('*','Save SRS data as');
	  writepfname = fullfile(writepname, writefname);
	  writedata = [fn' x_pos' (abs(x_neg))' ];
	  fid = fopen(writepfname,'w');
	  fprintf(fid,'  %g  %g  %g\n',writedata');
	  fclose(fid);
%%
%   disp(' Enter output filename ');
%   SRS_filename = input(' ','s');
%
%   fid = fopen(SRS_filename,'w');
%   for j=1:jnum
%        fprintf(fid,'%7.2f %10.3f %10.3f \n',fn(j),x_pos(j),abs(x_neg(j)));
%   end
%   fclose(fid);
end
%
%  Plot SRS
%
disp(' ')
disp(' Plotting output..... ')
%
%  Find limits for plot
%
srs_max = max(x_pos);
if max( abs(x_neg) ) > srs_max
    srs_max = max( abs(x_neg ));
end
srs_min = min(x_pos);
if min( abs(x_neg) ) < srs_min
    srs_min = min( abs(x_neg ));
end  
%
figure(1);
plot(fn,x_pos,fn,abs(x_neg),'-.');
%
if iunit==1
    ylabel('Peak Accel (G)');
else
    ylabel('Peak Accel (m/sec^2)');
end
xlabel('Natural Frequency (Hz)');
Q=1./(2.*damp);
out5 = sprintf(' Acceleration Shock Response Spectrum Q=%g ',Q);
title(out5);
grid;
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
legend ('positive','negative',2);
%
ymax= 10^(round(log10(srs_max)+0.8));
ymin= 10^(round(log10(srs_min)-0.6));
%
fmax=max(fn);
fmin=fmax/10.;
%
fmax= 10^(round(log10(fmax)+0.5));
%
if  fn(1) >= 0.1
    fmin=0.1;
end
if  fn(1) >= 1
    fmin=1;
end
if  fn(1) >= 10
    fmin=10;
end
if  fn(1) >= 100
    fmin=100;
end
axis([fmin,fmax,ymin,ymax]);
%
disp(' ')
disp(' Plot pseudo velocity? ');
vchoice=input(' 1=yes   2=no  ' );
if(vchoice==1)
figure(2);
%
%   Convert to pseudo velocity
%
for j=1:jnum  
    if iunit==1   
       x_pos(j)=386.*x_pos(j)/(2.*pi*fn(j));
       x_neg(j)=386.*x_neg(j)/(2.*pi*fn(j));   
    else
       x_pos(j)=x_pos(j)/(2.*pi*fn(j));
       x_neg(j)=x_neg(j)/(2.*pi*fn(j));   
    end
end    
%
srs_max = max(x_pos);
if max( abs(x_neg) ) > srs_max
    srs_max = max( abs(x_neg ));
end
srs_min = min(x_pos);
if min( abs(x_neg) ) < srs_min
    srs_min = min( abs(x_neg ));
end  
%
plot(fn,x_pos,fn,abs(x_neg),'-.');
%
if iunit==1
    ylabel('Velocity (in/sec)');
else
    ylabel('Velocity (m/sec)');   
end
xlabel('Natural Frequency (Hz)');
Q=1./(2.*damp);
out5 = sprintf(' Pseudo Velocity Shock Response Spectrum Q=%g ',Q);
title(out5);
grid;
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
legend ('positive','negative',2);
%
ymax= 10^(round(log10(srs_max)+0.8));
ymin= 10^(round(log10(srs_min)-0.6));
%
axis([fmin,fmax,ymin,ymax]);
end

Contact us