disp(' ')
disp(' srs.m   ver 3.5   October 4, 2011')
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;
clear x_std;
clear rd_pos;
clear rd_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 ');
disp('   3=Excel file');
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]);
        fclose(fid);
        THM=THM';
end
if(file_choice==2)
        THM = input(' Enter the matrix name:  ');
end
if(file_choice==3)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        THM = xlsread(xfile);
%         
end
%
t=double(THM(:,1));
y=double(THM(:,2));
%
ifn=1;
clear figure(ifn);
figure(ifn);
ifn=ifn+1;
plot(t,y);
grid on;
ylabel('Accel(G)');
xlabel('Time(sec)');
%
tmx=max(t);
tmi=min(t);
nnn = size(y);
n=max(nnn);
%
out1 = sprintf('\n  %d samples ',n);
disp(out1)
%
dt=(tmx-tmi)/(n-1);
sr=1./dt;
%
disp(' ')
disp(' Time Step ');
dtmin=min(diff(t));
dtmax=max(diff(t));
%
out4 = sprintf(' dtmin  = %8.4g sec  ',dtmin);
out5 = sprintf(' dt     = %8.4g sec  ',dt);
out6 = sprintf(' dtmax  = %8.4g sec  ',dtmax);
disp(out4)
disp(out5)
disp(out6)
%
disp(' ')
disp(' Sample Rate ');
out4 = sprintf(' srmin  = %8.4g samples/sec  ',1/dtmax);
out5 = sprintf(' sr     = %8.4g samples/sec  ',1/dt);
out6 = sprintf(' srmax  = %8.4g samples/sec  \n',1/dtmin);
disp(out4)
disp(out5)
disp(out6)
%
ncontinue=1;
if(((dtmax-dtmin)/dt)>0.01)
    disp(' ')
    disp(' Warning:  time step is not constant.  Continue calculation? 1=yes 2=no ')
    ncontinue=input(' ');
end
if(ncontinue==1)
%
    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(' ');
%
disp(' ');
disp(' Include residual? ');
disp('  1=yes  2=no ')
ire=input(' ');
%
if(ire==1)
%
%   Add trailing zeros for residual response
%
    tmax=(tmx-tmi) + 1./fn(1);
else
    tmax=(tmx-tmi);
end
    
    limit = round( tmax/dt );
    n=limit;
    yy=zeros(1,limit);
    for i=1:max(nnn)
        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);
        x_std(j)= std(resp);
%
%%%%%% relative displacement %%%
%
        rd_a1(j)=2.*exp(-domegadt)*cosd; 
        rd_a2(j)=-exp(-2.*domegadt);
        rd_b1(j)=0.;
        rd_b2(j)=-(dt/omegad)*exp(-domegadt)*sind;
        rd_b3(j)=0;        
%
        rd_forward=[ rd_b1(j),  rd_b2(j),  rd_b3(j) ];    
        rd_back   =[     1, -rd_a1(j), -rd_a2(j) ];  
%        
        rd_resp=filter(rd_forward,rd_back,yy);
%        
        rd_pos(j)= max(rd_resp);
        rd_neg(j)= min(rd_resp);      
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
        jnum=j; 
        if  fn(j) > sr/8.
            break
        end
        fn(j+1)=fn(1)*(2. ^ (j*(1./12.)));    
    end
%
    if iunit==1
          rd_pos=rd_pos*386;
          rd_neg=rd_neg*386;
    end
%
    maximaxSRS=max(x_pos,abs(x_neg)); 
%
    fmax=0;
    zmax=0;
    for(i=1:max(size(fn)))
        if(x_pos(i)>zmax)
            zmax=x_pos(i);
            fmax=fn(i);
        end
        if(abs(x_neg(i))>zmax)
            zmax=abs(x_neg(i));
            fmax=fn(i);
        end 
    end
%
    if(iunit==1)
        out5 = sprintf('\n Absolute Peak is %10.5g G at %10.5g Hz ',zmax,fmax);
    else
        out5 = sprintf('\n Absolute Peak is %10.5g m/sec^2 at %10.5g Hz ',zmax,fmax);
    end
    disp(out5)
%
%  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 positive & negative SRS data as');
	    writepfname = fullfile(writepname, writefname);
	    writedata = [fn' x_pos' (abs(x_neg))' ];
	    fid = fopen(writepfname,'w');
	    fprintf(fid,'  %g \t %g \t %g\n',writedata');
	    fclose(fid);
%%
        [writefname, writepname] = uiputfile('*','Save maximax SRS data as');
	    writepfname = fullfile(writepname, writefname);
	    writedata = [fn' maximaxSRS' ];
	    fid = fopen(writepfname,'w');
	    fprintf(fid,'  %g \t %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  
%
    disp(' select plot type:  1=positive & negative   2=maximax ');
    ipt=input(' ');
%
    clear figure(ifn);
    figure(ifn);
    ifn=ifn+1;
%
    if(ipt==1)
        plot(fn,x_pos,fn,abs(x_neg),'-.');
        legend ('positive','negative',2); 
    else
        plot(fn,maximaxSRS);
    end
%
    clear SRS_max;
    SRS_pn=[fn',x_pos',abs(x_neg)'];
    SRS_max=[fn',maximaxSRS'];
%
    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',':','GridLineStyle',':','XScale','log','YScale','log');
%
    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.001
        fmin=0.001;
    end
    if  fn(1) >= 0.01
        fmin=0.01;
    end
    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(' Matlab matrices: ')
    disp('      SRS_pn  - Acceleration SRS positive & negative ')
    disp('      SRS_max - Acceleration SRS maximax             ') 
%
    disp(' ')
    disp(' Plot pseudo velocity? ');
    vchoice=input(' 1=yes   2=no  ' );
    if(vchoice==1)
        disp(' select plot type:  1=positive & negative   2=maximax ');
        ipt=input(' ');     
%
%   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    
%
        clear figure(ifn);
        figure(ifn);
        ifn=ifn+1;  
        clear xpn;
        xpn=max(x_pos,abs(x_neg)); 
        if(ipt==1)
            plot(fn,x_pos,fn,abs(x_neg),'-.');
            legend ('positive','negative',2);   
        else     
            plot(fn,xpn);    
        end
        pseudo_velocity_pn=[fn',x_pos',abs(x_neg)'];
        pseudo_velocity_max=[fn',xpn'];  
%
        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',':','GridLineStyle',':','XScale','log','YScale','log');
%
        srs_max=max(xpn);
        srs_min=min(xpn);
        ymax= 10^(round(log10(srs_max)+0.8));
        ymin= 10^(round(log10(srs_min)-0.6));
%
        axis([fmin,fmax,ymin,ymax]);
 %
        disp(' Matlab matrices: ')
        disp('      pseudo_velocity_pn  - pseudo_velocity SRS positive & negative ')
        disp('      pseudo_velocity_max - pseudo_velocity SRS maximax             ')        
    end
%
    disp(' ');
    disp(' Plot relative displacement? ');
    rdchoice=input(' 1=yes   2=no  ' );
%
    if(rdchoice==1)
        disp(' select plot type:  1=positive & negative   2=maximax ');
        ipt=input(' ');
        clear figure(ifn);
        figure(ifn);
        ifn=ifn+1;
        clear rd_pn;
        rd_pn=max(rd_pos,abs(rd_neg));        
        if(ipt==1)
           plot(fn,rd_pos,fn,abs(rd_neg),'-.');
           legend ('positive','negative',2);        
        else
           plot(fn,rd_pn);
        end       
        if iunit==1
            ylabel('Rel Disp (in)');
        else
            ylabel('Rel Disp (m)');   
        end        
        xlabel('Natural Frequency (Hz)');
        Q=1./(2.*damp);
        out5 = sprintf(' Relative Displacement Shock Response Spectrum Q=%g ',Q);
        title(out5);
        grid;
        set(gca,'MinorGridLineStyle',':','GridLineStyle',':','XScale','log','YScale','log');
%
        rel_disp_pn=[fn',rd_pos',abs(rd_neg)'];
        rel_disp_max=[fn',rd_pn'];  
%
        disp(' Matlab matrices: ')
        disp('      rel_disp_pn  - relative displacement SRS positive & negative ')
        disp('      rel_disp_max - relative displacement SRS maximax             ')
%
    end    
%
    disp(' ')
    disp(' Plot std dev response spectrum? ');
    sdchoice=input(' 1=yes   2=no  ' );
    if(sdchoice==1)
    clear figure(ifn);
    figure(ifn);
        ifn=ifn+1;
        plot(fn,x_std);
        ylabel('Std Dev (G) ');
        xlabel('Natural Frequency (Hz)');
        Q=1./(2.*damp);
        out5 = sprintf(' Std Dev Shock Response Spectrum Q=%g ',Q);
        title(out5);
        grid;
        set(gca,'MinorGridLineStyle',':','GridLineStyle',':','XScale','log','YScale','log'); 
    end     
%
end
%
clear srs_max;
clear srs_min;