disp(' ')
disp(' srs_tripartite.m   ver 2.2   June 27, 2013')
disp(' by Tom Irvine   Email: tom@vibrationdata.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(' ')
%
close all;
%
clear length;
clear t;
clear y;
clear yy;
clear n;
clear fn;
clear a1;
clear a2
clear b1;
clear b2;
clear THM;
clear resp;
clear acc_pos;
clear acc_neg;
clear pv_pos;
clear pv_neg;
clear pseudo_velocity_pn;
clear pseudo_velocity_max;
clear RD_SRS_max;
clear PV_SRS_max;
%
%
while(1)
    disp(' Select units ');
    disp('  1=English: accel(G),       vel(in/sec),  disp(in)  ');
    disp('  2=metric : accel(G),        vel(m/sec),  disp(mm)  ');
    disp('  3=metric : accel(m/sec^2),  vel(m/sec),  disp(mm)  ');   
    iunit=input(' ');
    if(iunit==1 || iunit==2 || iunit==3)
        break;
    end
end	
%
[t,y,dt,sr,tmx,tmi,n,ncontinue]=enter_time_history();
%
fig_num=1;
figure(fig_num);
fig_num=fig_num+1;
plot(t,y);
grid on;
if(iunit==1 || iunit==2)
    ylabel('Accel(G)');
else
    ylabel('Accel(m/sec^2)');    
end
xlabel('Time(sec)');
%
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(' Include residual? ');
    disp('  1=yes  2=no ')
    ire=input(' ');
%
    disp(' ');
    disp(' Select method for calculating psuedo velocity ');
    disp('  1=APVSS from absolute acceleration ');
    disp('  2=PVSS  from relative displacement ');
    imethod=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 );
    yy=zeros(1,limit);
    for i=1:length(y)
        yy(i)=y(i);
    end    
%
    disp(' ')
    disp(' Calculating response..... ')
%
%  SRS engine
%
j=1;
while(1)
    if (fn(j) > sr/8.)
        break
    end
    fn(j+1)=fn(1)*(2. ^ (j*(1./12.))); 
    j=j+1;
end
%
[a1,a2,b1,b2,b3,rd_a1,rd_a2,rd_b1,rd_b2,rd_b3]=...
                                   srs_coefficients(fn,damp,dt);
%
    ns=max(size(fn));    
%
    fn=fix_size(fn);
    omegan=2*pi*fn;    
%
    acc_pos=zeros(ns,1);
    acc_neg=zeros(ns,1);  
     pv_pos=zeros(ns,1);
     pv_neg=zeros(ns,1);
     rd_pos=zeros(ns,1);
     rd_neg=zeros(ns,1);     
%
    PV_SRS_max=zeros(ns,2);
    RD_SRS_max=zeros(ns,2);
%
    RD_SRS_max(:,1)=fn;
    PV_SRS_max(:,1)=fn;   
     A_SRS_max(:,1)=fn; 
%
    if(imethod==1)
%        
        for j=1:ns
%
            forward=[ b1(j),  b2(j),  b3(j) ];    
            back   =[     1, -a1(j), -a2(j) ];    
%    
            acc_resp=filter(forward,back,yy);
%
            acc_pos(j)= max(acc_resp);
            acc_neg(j)= min(acc_resp);
%    
        end
%
        acc_pos=fix_size(acc_pos);
        acc_neg=fix_size(acc_neg);
%    
        maximaxSRS=max(acc_pos,abs(acc_neg)); 
        A_SRS_pn=[fn,acc_pos,abs(acc_neg)];
        A_SRS_max=[fn,maximaxSRS]; 
%
        for i=1:ns
            PV_SRS_max(i,2)=A_SRS_max(i,2)/omegan(i);
            RD_SRS_max(i,2)=A_SRS_max(i,2)/omegan(i)^2;
        end
%
    else
%
        for j=1:ns
%
            forward=[ rd_b1(j),  rd_b2(j),  rd_b3(j) ];    
            back   =[     1, -rd_a1(j), -rd_a2(j) ];    
%    
            rd_resp=filter(forward,back,yy);
%
            rd_pos(j)= max(rd_resp);
            rd_neg(j)= min(rd_resp);
%
        end
%
        rd_pos=fix_size(rd_pos);
        rd_neg=fix_size(rd_neg);
%
        if(iunit==1)
          scale=386;
        end
        if(iunit==2)
          scale=9.81*1000;
        end
        if(iunit==3)
           scale=1000;
        end
        rd_pos=rd_pos*scale;
        rd_neg=rd_neg*scale;        
%
        maximaxSRS=max(rd_pos,abs(rd_neg)); 
        RD_SRS_pn=[fn,rd_pos,abs(rd_neg)];
        RD_SRS_max=[fn,maximaxSRS]; 
%
        if(iunit==1)
          vscale=1;
          ascale=1/386;
        end
        if(iunit==2)
          vscale=1/1000;
          ascale=1/(9.81*1000);
        end
        if(iunit==3)
          vscale=1/1000;
          ascale=1/1000;
        end
%
        for i=1:ns
            pv_pos(i)=vscale*rd_pos(i)*omegan(i);
            pv_neg(i)=vscale*rd_neg(i)*omegan(i);
            PV_SRS_max(i,2)=vscale*RD_SRS_max(i,2)*omegan(i);
%            
            A_SRS_max(i,2)=ascale*RD_SRS_max(i,2)*omegan(i)^2;
        end
%
        pseudo_velocity_pn=[fn,pv_pos,pv_neg];
        pseudo_velocity_max=PV_SRS_max;
%
    end
%    
    [zmax,fzmax]=find_max(A_SRS_max);
%      
    [v_zmax,v_fzmax]=find_max(PV_SRS_max);    
    [rd_zmax,rd_fzmax]=find_max(RD_SRS_max);  
%
    disp(' ');
    disp(' Absolute Peak Values ');
    disp(' ');
    if(iunit==1)
        out5 = sprintf('    Acceleration: %9.4g G at %9.4g Hz ',zmax,fzmax);
        out6 = sprintf(' Pseudo Velocity: %9.4g in/sec at %9.4g Hz ',v_zmax*386,v_fzmax);
        out7 = sprintf('   Relative Disp: %9.4g in at %9.4g Hz ',rd_zmax*386,rd_fzmax);        
    end
    if(iunit==2)
        out5 = sprintf('    Acceleration: %9.4g G at %9.4g Hz ',zmax,fzmax);
        out6 = sprintf(' Pseudo Velocity: %9.4g m/sec at %9.4g Hz ',v_zmax*9.81,v_fzmax);  
        out7 = sprintf('   Relative Disp: %9.4g mm at %9.4g Hz ',rd_zmax*9810,rd_fzmax);         
    end
    if(iunit==3)
        out5 = sprintf('    Acceleration: %9.4g m/sec^2 at %9.4g Hz ',zmax,fzmax);
        out6 = sprintf(' Pseudo Velocity: %9.4g m/sec at %9.4g Hz ',v_zmax,v_fzmax);  
        out7 = sprintf('   Relative Disp: %9.4g mm at %9.4g Hz ',rd_zmax*1000,rd_fzmax);         
    end    
%
    disp(out5)
    disp(out6)
    disp(out7)
%    
    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
%
    tpi_fmin=(2*pi*fmin);
    tpi_fmax=(2*pi*fmax);
%
    disp(' '); 
    disp(' Matlab matrices: ')
    disp('      A_SRS_pn  - Acceleration SRS positive & negative ')
    disp('      A_SRS_max - Acceleration SRS maximax             ') 
    disp(' ')
%
    disp('      pseudo_velocity_pn  - pseudo_velocity SRS positive & negative ')
    disp('      pseudo_velocity_max - pseudo_velocity SRS maximax             ') 
%
%  Output options
%
%  Plot SRS
%
        disp(' ');
        disp(' select plot type:  1=positive & negative   2=maximax ');
        ipt=input(' ');   
%
        if(imethod==1) % convert to pseudo velocity
%
            if iunit==1   
                acc_pos=386.*acc_pos;
                acc_neg=386.*acc_neg;   
            end
            if iunit==2
                acc_pos=9.81*acc_pos;
                acc_neg=9.81*acc_neg;   
            end
            for j=1:ns 
                pv_pos(j)=acc_pos(j)/omegan(j);
                pv_neg(j)=acc_neg(j)/omegan(j);  
            end
%
            pseudo_velocity_pn=[fn,pv_pos,abs(pv_neg)];
            xpn=max(pv_pos,abs(pv_neg));
            pseudo_velocity_max=[fn,xpn];
        end
%
        Q=1/(2*damp);
%        
    [~]=plot_srs_tripartite(fn,pv_pos,pv_neg,Q,fig_num,ipt,fmin,fmax,iunit);
%
end