disp(' ')
disp(' sdof_sine_force.m ')
disp(' ver 1.1  February 7, 2012 ')
disp(' by Tom Irvine  Email: tomirvine@aol.com ')
disp(' ')
disp(' This program calculates the vibration of a  ')
disp(' single-degree-of-freedom system subjected to ')
disp(' a sine force. ')
disp(' ')
%
close all;
%
clear omegan;
clear omegad;
clear fn;
clear damp;
clear N;
clear t;
clear d;
clear v;
clear acc;
%
disp(' select units')
iu=input(' 1=English  2=metric ');
%
disp(' ');
if(iu==1)
   disp(' Enter the amplitude (lbf) ');   
else
   disp(' Enter the amplitude (N) ');   
end
F=input(' ');
%
disp(' ');
disp(' Enter the frequency (Hz) ');
freq=input(' ');
beta=2*pi*freq;
beta2=beta^2;
%
if(iu==1)
    disp(' ')
    v0 = input(' input initial velocity (in/sec) ');
    disp(' ')
    d0 = input(' input initial displacement (in) ');
    disp(' ')
    mass=input(' input mass (lbm) ');
    mass=mass/386;
else
    disp(' ')
    v0 = input(' input initial velocity (m/sec) ');
    disp(' ')
    d0 = input(' input initial displacement (m) '); 
    disp(' ')
    mass=input(' input mass (kg) ');    
end
%
disp(' ')
fn=input(' Enter the natural frequency (Hz)  ');
%
omegan=2.*pi*fn;
omegan2=omegan^2;
stiff=mass*omegan2;
M=mass;
K=stiff;
%
if(iu==1)
    out1=sprintf('         stiffness =%8.4g lbf/in \n',stiff);
    out2=sprintf(' static deflection =%8.4g in \n',1/stiff);
else
    out1=sprintf('         stiffness =%8.4g  N/m \n',stiff);
    out2=sprintf(' static deflection =%8.4g m \n',1/stiff);   
end
disp(out1);
disp(out2);
%
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
%
C=2*damp*omegan*mass;
%
omegad=omegan*sqrt(1.-damp^2);
domegan=damp*omegan;
dterm=(1-2*damp^2);
%
dF=(beta2-omegan2)^2+(2*beta*domegan)^2;
FdF=(F/mass)/dF;
%
%
U1=2*beta*domegan;
V1=(beta2-omegan2);
%
U2=(2*domegan*omegad);
V2=(beta2-omegan2*dterm);
%
W=(v0+domegan*d0)/omegad;
%
period=1/fn;
dt=min([period (1/freq)])/64;
sr=1/dt;
%
disp(' ');
disp(' Enter the duration (sec) ');
dur=input(' ');
%
N=round(dur/dt);
%
t=linspace(0,dur,N);
t=t';
%
d=zeros(N,1);
v=zeros(N,1);
a=zeros(N,1);
Fin=zeros(N,1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
for(i=1:N)
%	
    tt=t(i);
%
        Fin(i)=F*sin(beta*tt);
        eee=exp(-domegan*tt); 
        cosbt=cos(beta*tt);  
        sinbt=sin(beta*tt);  
        cosd=cos(omegad*tt);
        sind=sin(omegad*tt);
%
        deee=-domegan*eee;
        dcosbt=-beta*sinbt;  
        dsinbt= beta*cosbt;  
        dcosd=-omegad*sind;
        dsind= omegad*cosd;
%
        d0 =eee*(d0*cosd + W*sind);
        d1 =FdF*(U1*cosbt +V1*sinbt);  
        d2 =FdF*(beta*eee/omegad)*( U2*cosd+ V2*sind);    
%
        d(i)=d0-d1+d2;
%
        v1=deee*(d0*cosd + W*sind) + eee*(d0*dcosd + W*dsind);
        v2=FdF*(U1*dcosbt + V1*dsinbt);
        v3a=FdF*(beta*deee/omegad)*( U2*cosd+ V2*sind);
        v3b=FdF*(beta*eee/omegad)*( U2*dcosd+ V2*dsind); 
%
        v(i)=v1-v2+v3a+v3b;
%
end
%
a=(Fin-C*v-K*d)/M;
%
disp(' ');
figure(1);
plot(t,d);
grid on;
title('Displacement');
if(iu==1)
        ylabel('Displacement (in)');
else
        ylabel('Displacement (m)');
end
xlabel('Time (sec)');
%
figure(2);
plot(t,v);
grid on;   
title('Velocity');
if(iu==1)
   ylabel('Velocity (in/sec)');
else
   ylabel('Velocity (m/sec)');
end
xlabel('Time (sec)');
%
figure(3);
if(iu==1)
   a=a/386.;
end     
plot(t,a);
grid on;  
title('Acceleration');
if(iu==1)
   ylabel('Accel (G)');
else
   ylabel('Accel (m/sec^2)');
end
xlabel('Time (sec)');
%
clear FA;
clear FV;
clear FD;
FA=[ t Fin a];
FV=[ t Fin v];
FD=[ t Fin d];