disp(' ')
disp(' sdof_force_accel.m  version 1.1    January 13, 2011 ')
disp(' By Tom Irvine   Email:  tomirvine@aol.com ')
disp(' ')
disp(' This program calculates the response of')
disp(' a single-degree-of-freedom system subjected')
disp(' to combined sinusoidal force and base excitation.')
disp(' ')
disp(' Assume zero initial displacement and zero initial velocity.')
%
clear TT;
clear acc;
clear rd;
clear fn;
clear omegan;
clear m;
clear Q;
clear A;
clear F;
%
tpi=2.*pi;
%
disp(' ');
fn = input(' Enter the natural frequency (Hz) ');
Q = input(' Enter amplification factor Q ');
damp=1./(2.*Q);
disp(' ');
m = input(' Enter the mass(lbm) ');
m=m/386;
%
disp(' ');
dur=input(' Enter duration (sec)  ');
disp(' ');
SR=input(' Enter the sample rate (sample/sec)  ');
%
dt=1/SR;
nt=(dur*SR);
%
clear acc;
clear rd;
%
%   read in acceleration and force
%
disp(' ');
A = input(' Enter base acceleration amplitude (G) ');
A=A*386;
alpha=input(' Enter base acceleration frequency (Hz) ');
alpha=alpha*tpi;
%
disp(' ');
F = input(' Enter force amplitude (lbf) ');
F=F/m;
beta=input(' Enter force frequency (Hz) ');
beta=beta*tpi;
%
disp(' ');
%
omegan=tpi*fn;
omegad=omegan*sqrt(1-damp^2);
alpha2=alpha^2;
beta2=beta^2;
%
omegan2=omegan^2;
domegan=damp*omegan;
damp2=damp^2;
%
dterm=(1-2*damp2);
%    
omratio=omegan/omegad;
dA=(alpha2-omegan2)^2+(2*alpha*domegan)^2;
dF=(beta2-omegan2)^2+(2*beta*domegan)^2;
%
acc=zeros(nt+1,1);
rd=zeros(nt+1,1);
%
for(i=0:nt)
%		
	t=dt*i;
    TT(i+1)=t;
%
    ee=exp(-domegan*t); 
    cosat=cos(alpha*t);
    cosbt=cos(beta*t);  
    sinat=sin(alpha*t);
    sinbt=sin(beta*t);  
    cosd=cos(omegad*t);
    sind=sin(omegad*t);   
%
    rd_alpha1=(A/dA)*((2*alpha*domegan)*cosat+(alpha2-omegan2)*sinat);
    rd_beta1 =(F/dF)*((2*beta*domegan)*cosbt +(beta2-omegan2)*sinbt);
%  
    rd_alpha2=(A/dA)*(alpha*ee/omegad)*((2*domegan*omegad)*cosd+(alpha2-omegan2*dterm)*sind);    
    rd_beta2 =(F/dF)*(beta*ee/omegad)*( (2*domegan*omegad)*cosd+ (beta2-omegan2*dterm)*sind);    
%
    rd(i+1)=rd_alpha1-rd_alpha2-rd_beta1+rd_beta2;
%
    ra_alpha1=(A/dA)*(alpha2)*(-(alpha2-omegan2)*sinat-(2*alpha*domegan)*cosat);
    ra_beta1 =(F/dF)*(beta2)*( -(beta2-omegan2 )*sinbt-(2*beta*domegan )*cosbt);
%    
    ra_alpha2=(A/dA)*(alpha*omegan*ee)*((2*damp*alpha2)*cosd+ omratio*(-omegan2+alpha2*dterm)*sind);    
    ra_beta2 =(F/dF)*(beta*omegan*ee)*( (2*damp*beta2)*cosd + omratio*(-omegan2+beta2*dterm )*sind);    
%
    ra=ra_alpha1+ra_alpha2-ra_beta1-ra_beta2;
%
    acc(i+1)=ra+A*sin(alpha*t);
%
end
%
acc=acc/386;
out1 = sprintf('\n maximum acceleration = %12.4g G',max(acc));
out2 = sprintf(' minimum acceleration = %12.4g G \n',min(acc));
disp(out1);
disp(out2);
%
disp(' Plot the acceleration response time history ?');
choice = input(' 1=yes  2= no  ');
if(choice==1)
            figure(1);   
            plot(TT,acc);
            xlabel('Time (sec)');  
            ylabel('Accel (G)');    
            grid;
            set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','lin','YScale','lin'); 
            out5 = sprintf(' SDOF Response, fn=%g Hz, Q=%g  ',fn,Q);
            title(out5);         
end    
%
out1 = sprintf('\n\n maximum relative disp = %12.4g inch',max(rd));
out2 = sprintf(' minimum relative disp = %12.4g inch ',min(rd));
disp(out1);
disp(out2);
%
disp(' ');  
disp(' Plot the relative displacement time history ?');
choice = input(' 1=yes  2= no  ');
if(choice==1)
            figure(2); 
            plot(TT,rd);
            xlabel('Time (sec)');  
            ylabel('Relative Disp (inch)');    
            grid;
            set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','lin','YScale','lin'); 
            out5 = sprintf(' SDOF Response, fn=%g Hz, Q=%g  ',fn,Q);
            title(out5);         
end            
%
clear acceleration;
clear relative_displacement;
acceleration=[TT',acc];
relative_displacement=[TT',rd];
%
disp(' ');
disp(' Output Files: ');
disp('  acceleration ');
disp('  relative_displacement ');
%