disp(' ')
disp(' arbit_force.m   ver 1.4   April 27, 2012')
disp(' by Tom Irvine   Email: tomirvine@aol.com')
disp(' ')
disp(' This program calculates the response 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;
%
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)')
else
    disp(' The input time history must have two columns: time(sec) & force(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
%
disp(' ')
fn(1)=input(' Enter the natural frequency (Hz)  ');
%
disp(' ')
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
%
%  Initialize coefficients
%
disp(' ')
disp(' Initialize coefficients')
%for j=1:1
    omegan=2.*pi*fn(1);
    omegad=omegan*sqrt(1.-(damp^2));
    domegan=damp*omegan;
%
    out5 = sprintf(' omegan=%g   omegad=%g ',omegan,omegad);
    disp(out5);    
%    
    cosd=cos(omegad*dt);
    sind=sin(omegad*dt);
%
    out5 = sprintf(' cosd=%g   sind=%g ',cosd,sind);
    disp(out5);    
%
    domegadt=domegan*dt;
%
    eee1=exp(-domegadt);
    eee2=exp(-2.*domegadt);
%
    ecosd=eee1*cosd;
    esind=eee1*sind; 
%
    a1= 2.*ecosd;
    a2=-eee2;
%
    out5 = sprintf(' a1=%g  a2=%g ',a1(1),a2(1));
    disp(out5);     
%
    omeganT=omegan*dt;
%
    phi=2*(damp)^2-1;
    DD1=(omegan/omegad)*phi;
    DD2=2*DD1;
%    
    df1=2*damp*(ecosd-1) +DD1*esind +omeganT;
    df2=-2*omeganT*ecosd +2*damp*(1-eee2) -DD2*esind;
    df3=(2*damp+omeganT)*eee2 +(DD1*esind-2*damp*ecosd);
%     
    VV1=-(damp*omegan/omegad);
%    
    vf1=(-ecosd+VV1*esind)+1;
    vf2=eee2-2*VV1*esind-1;
    vf3=ecosd+VV1*esind-eee2;
%    
%end
%
MD=(mass*omegan^3*dt);
df1=df1/MD;
df2=df2/MD;
df3=df3/MD;
%
VD=(mass*omegan^2*dt);
vf1=vf1/VD;
vf2=vf2/VD;
vf3=vf3/VD;
%
af1=esind/(mass*omegad*dt);
af2=-2*af1;
af3=af1;
%
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(1);
    limit = round( tmax/dt );
    n=limit;
end
%   
ff=zeros(1,n);
for i=1:length(force)
        ff(i)=force(i);        
end     
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Numerical Engine
%
disp(' ')
disp(' Calculating response...');
%
%  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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
t=fix_size(t);
force=fix_size(force);
%
fig_num=1;
figure(fig_num);
fig_num=fig_num+1;
plot(t,force);
if(iu==1)
  ylabel('Force(lbf)');   
else
  ylabel('Force(N)');   
end
xlabel('Time (sec)');
title('Applied Force');
grid on;
zoom on;
%
nn=length(ff);
tmin=t(1);
tmax=(t(1)+(nn-1)*dt);
%
tt=linspace(tmin,tmax,nn);
%
Q=1./(2.*damp);
%
tt=fix_size(tt);
d_resp=fix_size(d_resp);
v_resp=fix_size(v_resp);
a_resp=fix_size(a_resp);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ')
disp(' Displacement Response ')
mu=mean(d_resp);
sd=std(d_resp);
d_rms=sqrt(sd^2+mu^2);
%
out10=sprintf('\n maximum = %10.3g %s',max(d_resp),disp_unit);
out11=sprintf(' minimum = %10.3g %s',min(d_resp),disp_unit);
out12=sprintf(' overall = %10.3g %s RMS',d_rms,disp_unit);
%
disp(out10)
disp(out11)
disp(out12)
%
figure(fig_num);
fig_num=fig_num+1;
disp(' ');
plot(tt,d_resp);
axis([tmin,tmax,min(ylim),max(ylim)]);      
out1=sprintf('Disp (%s)',disp_unit);
ylabel(out1);
xlabel('Time (sec)');
out5 = sprintf(' SDOF Displacement Response   fn=%g Hz   Q=%g ',fn(1),Q);
title(out5);
grid on;
zoom on;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ')
disp(' Velocity Response ')
mu=mean(v_resp);
sd=std(v_resp);
v_rms=sqrt(sd^2+mu^2);
out10=sprintf('\n maximum = %10.3g %s',max(v_resp),vel_unit);
out11=sprintf(' minimum = %10.3g %s',min(v_resp),vel_unit);
out12=sprintf(' overall = %10.3g %s RMS',v_rms,vel_unit);
disp(out10)
disp(out11)
disp(out12)
%
figure(fig_num);
fig_num=fig_num+1;
plot(tt,v_resp);
axis([tmin,tmax,min(ylim),max(ylim)]);  
out1=sprintf('Vel (%s)',vel_unit);
ylabel(out1);
xlabel('Time (sec)');
out5 = sprintf(' SDOF Velocity Response   fn=%g Hz   Q=%g ',fn(1),Q);
title(out5);
grid on;
zoom on;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ')
disp(' Acceleration Response ')
mu=mean(a_resp);
sd=std(a_resp);
a_rms=sqrt(sd^2+mu^2);
out10=sprintf('\n maximum = %10.3g %s',max(a_resp),acc_unit);
out11=sprintf(' minimum = %10.3g %s',min(a_resp),acc_unit);
out12=sprintf(' overall = %10.3g %s RMS',a_rms,acc_unit);
disp(out10)
disp(out11)
disp(out12)
%    
figure(fig_num);
fig_num=fig_num+1;
plot(tt,a_resp);
axis([tmin,tmax,min(ylim),max(ylim)]);    
out1=sprintf('Accel (%s)',acc_unit);
ylabel(out1);
xlabel('Time (sec)');
out5 = sprintf(' SDOF Acceleration Response   fn=%g Hz   Q=%g ',fn(1),Q);
title(out5);
grid on;
zoom on;
%
clear acceleration;
clear velocity;
clear displacement;
%
acceleration=[tt a_resp];
velocity=[tt v_resp];
displacement=[tt d_resp];
%
end