disp(' ')
disp(' srs_force.m   ver 1.0   December 20, 2013')
disp(' by Tom Irvine   Email: tom@vibrationdata.com')
disp(' ')
disp(' This program calculates the SRS of a single-degree-of-freedom')
disp(' system to an arbitrary force input time history. ')
%
close all;
%
clear omega;
clear omegad;
clear dt;
clear sr;
clear damp;
clear a1;
clear a2;
clear b1;
clear b2;
clear c1;
clear c2;
clear sum;
clear t;
clear tt;
clear y;
clear yy;
clear THM;
clear x_pos;
clear x_neg;
clear d_pos;
clear d_neg;
clear a_resp;
clear v_resp;
clear x_resp;
%
clear forward;
clear back;
clear d_forward;
clear d_back;
clear v_forward;
clear v_back;
clear a_forward;
clear a_back;
%
clear length;
clear mass;
clear df;
clear force;
clear ff;
%
close all;
%
fig_num=1;
%
disp(' ');
disp(' Enter units:  1=English 2=metric ');
iu=input(' ');
%
disp(' ')
if(iu==1)
    disp(' The input time history must have two columns: time(sec) & force(lbf)')
    force_unit='lbf';
else
    disp(' The input time history must have two columns: time(sec) & force(N)')
    force_unit='N';
end
disp(' ')
%
[t,force,dt,sr,tmx,tmi,n,ncontinue]=enter_time_history();
%
if(ncontinue==1)
%
disp(' ')
if(iu==1)
    mass=input(' Enter the mass (lbm)  ');
    mass=mass/386;
    disp_unit='inch';
    vel_unit='in/sec';
    acc_unit='G';
else
    mass=input(' Enter the mass (kg)  '); 
    disp_unit='m';
    vel_unit='m/sec';
    acc_unit='G';   
end
%
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
%
j=1;
while(1)
    if (fn(j) > sr/8.)
        break
    end
    fn(j+1)=fn(1)*(2. ^ (j*(1./12.))); 
    j=j+1;
end
%
%  Initialize coefficients
%
disp(' ')
disp(' Initialize coefficients')
%
num=length(fn);
%
disp(' ');
disp(' Include residual? ');
disp('  1=yes  2=no ')
ire=input(' ');
%
if(ire==1)
%
%   Add trailing zeros for residual response
%
    disp(' ')
    disp(' Add trailing zeros for residual response')
    tmax=(tmx-tmi) + 2./fn;
    limit = round( tmax/dt );
    n=limit;
end
%   
ff=zeros(1,n);
for i=1:length(force)
        ff(i)=force(i);        
end     
%
force=ff;
ff=ff/mass;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Numerical Engine
%
%
accel_srs=zeros(num,3);
velox_srs=zeros(num,3);
displ_srs=zeros(num,3);
trans_force_srs=zeros(num,3);
%
disp(' ')
disp(' Calculating response...');
%
for i=1:num
%
omegan=2*pi*fn(i);
[a1,a2,df1,df2,df3,vf1,vf2,vf3,af1,af2,af3]=...
                   ramp_invariant_filter_coefficients(1,omegan,damp,dt);
%
%  displacement
%
    d_forward=[   df1,  df2, df3 ];
    d_back   =[     1, -a1, -a2 ];
    d_resp=filter(d_forward,d_back,ff);
%    
%  velocity
%
    v_forward=[   vf1,  vf2, vf3 ];
    v_back   =[     1, -a1, -a2 ];
    v_resp=filter(v_forward,v_back,ff);
%    
%  acceleration
%   
    a_forward=[   af1,  af2, af3 ];
    a_back   =[     1, -a1, -a2 ]; 
    a_resp=filter(a_forward,a_back,ff);
%
    if(iu==1)
        a_resp=a_resp/386;  
    else
        a_resp=a_resp/9.81;     
    end
%
%  transmitted foree
%
    fnf=fn(i);
    [a1,a2,b1,b2,b3,rd_a1,rd_a2,rd_b1,rd_b2,rd_b3]=...
                                   srs_coefficients(fnf,damp,dt);
    back   =[     1, -a1, -a2 ];
    aa_forward=[ b1,  b2,  b3 ];
%    
    ft_resp=filter(aa_forward,back,force);
%
    accel_srs(i,2)=max(a_resp);
    accel_srs(i,3)=abs(min(a_resp));   
%
    velox_srs(i,2)=max(v_resp);
    velox_srs(i,3)=abs(min(v_resp));   
%
    displ_srs(i,2)=max(d_resp);
    displ_srs(i,3)=abs(min(d_resp));
%
    trans_force_srs(i,2)=max(ft_resp);
    trans_force_srs(i,3)=abs(min(ft_resp));
%
    accel_srs(i,1)=fn(i);
    velox_srs(i,1)=fn(i);
    displ_srs(i,1)=fn(i);
    trans_force_srs(i,1)=fn(i);
%
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
figure(fig_num);
fig_num=fig_num+1;
plot(t,force); 
xlabel('Time (sec)');
%
if(iu==1)
    ylabel('Force (lbf)');
else
    ylabel('Force (N)');
end
title('Applied Force');
%
grid;
zoom on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(fig_num);
fig_num=fig_num+1;
plot(fn,accel_srs(:,2),fn,accel_srs(:,3));
xlabel('Natural Frequency(Hz)');
grid on;
set(gca,'MinorGridLineStyle',':','GridLineStyle',':','XScale','log','YScale','log');
grid on;
ylabel('Peak Accel (G)');
if(iu==1)
    out5 = sprintf(' Acceleration SRS: mass=%g lbm Q=%g ',mass*386,Q);
else
    out5 = sprintf(' Acceleration SRS: mass=%g kg  Q=%g ',mass,Q);    
end
title(out5);
legend('positive','negative',1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(fig_num);
fig_num=fig_num+1;
plot(fn,velox_srs(:,2),fn,velox_srs(:,3));
xlabel('Natural Frequency(Hz)');
grid on;
set(gca,'MinorGridLineStyle',':','GridLineStyle',':','XScale','log','YScale','log');
grid on;
if(iu==1)
    ylabel('Peak Vel (in/sec)');
else
    ylabel('Peak Vel (m/sec)');
end
if(iu==1)
    out5 = sprintf(' Velocity SRS: mass=%g lbm Q=%g ',mass*386,Q);
else
    out5 = sprintf(' Velocity SRS: mass=%g kg  Q=%g ',mass,Q);    
end
title(out5);
legend('positive','negative',1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(fig_num);
fig_num=fig_num+1;
plot(fn,displ_srs(:,2),fn,displ_srs(:,3));
xlabel('Natural Frequency(Hz)');
grid on;
set(gca,'MinorGridLineStyle',':','GridLineStyle',':','XScale','log','YScale','log');
grid on;
if(iu==1)
    ylabel('Peak Disp (in)');
else
    ylabel('Peak Disp (mm)');
end
if(iu==1)
    out5 = sprintf(' Displacement SRS: mass=%g lbm Q=%g ',mass*386,Q);
else
    out5 = sprintf(' Displacement SRS: mass=%g kg  Q=%g ',mass,Q);    
end
title(out5);
legend('positive','negative',1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(fig_num);
fig_num=fig_num+1;
plot(fn,trans_force_srs(:,2),fn,trans_force_srs(:,3));
xlabel('Natural Frequency(Hz)');
grid on;
set(gca,'MinorGridLineStyle',':','GridLineStyle',':','XScale','log','YScale','log');
grid on;
if(iu==1)
    ylabel('Force (lbf)');
else
    ylabel('Force (N)');
end
if(iu==1)
    out5 = sprintf('Transmitted Force SRS: mass=%g lbm Q=%g ',mass*386,Q);
else
    out5 = sprintf('Transmitted Force SRS: mass=%g kg  Q=%g ',mass,Q);    
end
title(out5);
legend('positive','negative',1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end