disp(' ')
disp(' unit_step_force.m ')
disp(' ver 1.1  September 21, 2011 ')
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 unit step force. ')
disp(' ')
%
clear omegan;
clear omegad;
clear fn;
clear damp;
clear N;
clear t;
clear x;
clear v;
%
disp(' select units')
iunits=input(' 1=English  2=metric ');
%
if(iunits==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;
stiff=mass*omegan^2;
%
if(iunits==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
%
omegad=omegan*sqrt(1.-damp^2);
%
%   State Space Solution
%
period=1/fn;
dt=period/64;
sr=1/dt;
dur=15*period;
N=dur/dt;
%
t=linspace(0,dur,N);
t=t';
%
x=zeros(N,1);
v=zeros(N,1);
%
%
domegan=damp*omegan;
t1=-domegan*t;
t2=omegad*t;  
%    
if(damp<1)
%
    clear c;
    c=(v0+domegan*d0)/omegad;
    d=domegan/omegad;  
%
    x=exp(t1).*( d0*cos(t2) + c*sin(t2));
    v=-domegan*x + omegad*exp(t1).*( -d0*sin(t2) + c*cos(t2));
%
    x=x+(1. -exp(t1).*(cos(t2) + d*sin(t2)) )/stiff;
    v=v+domegan*exp(t1).*(cos(t2)  + d*sin(t2))/stiff;
    v=v- omegad*exp(t1).*(-sin(t2) + d*cos(t2))/stiff;
%
end
%
if(damp==1)
    x=exp(t1).*(d0+(v0+omegan.*d0).*t);
    v=-domegan*x + exp(t1).*(v0+omegan.*d0);
%
    alpha=omegan*t;
    x=x+(1-(1+alpha).*exp(-alpha))/stiff;
    v=v+(omegan*(1+alpha).*exp(-alpha))/stiff;
    v=v-omegan*exp(-alpha)/stiff;
%
end
%
if(damp>1)
%    
    t3=omegan*t;
    dr=sqrt(damp^2-1);
    C=1/(2*omegan*dr);
    alpha=-damp+dr;
    beta=-damp-dr;
%    
    A= v0+omegan*d0*(-beta);
    B=-v0+omegan*d0*(alpha);
%    
    x=C*(A*exp(alpha*t3)+B*exp(beta*t3));
    v=C*omegan*(A*alpha*exp(alpha*t3)+B*beta*exp(beta*t3));
%
    s1=omegan*(-damp+dr);
    s2=omegan*(-damp-dr);
%
    H=stiff*(s2-s1);
%
    es1=exp(s1*t);
    es2=exp(s2*t);
%
    des1=s1*exp(s1*t);
    des2=s2*exp(s2*t);
%
    x=x+(1/stiff)-(-s1*es1 +s2*es2)/H  +(2*damp*omegan)*(es1-es2)/H; 
    v=v-(-s1*des1 +s2*des2)/H  +(2*damp*omegan)*(des1-des2)/H;
%
end
%
disp(' ');
disp(' plot displacement?');
k=input(' 1=yes  2=no ');
if(k==1)
figure(1);
    plot(t,x);
    grid on;
    if(iunits==1)
        ylabel('Displacement (in)');
    else
        ylabel('Displacement (m)');
    end
    xlabel('Time (sec)');
end
%
disp(' ');
disp(' plot velocity?');
k=input(' 1=yes  2=no ');
if(k==1)
    figure(2);
    plot(t,v);
    grid on;   
    if(iunits==1)
        ylabel('Velocity (in/sec)');
    else
        ylabel('Velocity (m/sec)');
    end
    xlabel('Time (sec)');
end
%
disp(' ');
disp(' plot acceleration?');
k=input(' 1=yes  2=no ');
if(k==1)
    figure(3);
    clear a;
    a= (1/mass)-(2.*domegan)*v -(omegan^2)*x;
    if(iunits==1)
        a=a/386.;
    end     
    plot(t,a);
    grid on;  
    if(iunits==1)
        ylabel('Accel (G)');
    else
        ylabel('Accel (m/sec^2)');
    end
    xlabel('Time (sec)');
end
%
disp(' ');
disp(' plot phase portraits?');
k=input(' 1=yes  2=no ');
if(k==1)
    figure(4);
    plot(x,v);
    grid on;
    if(iunits==1)
        ylabel('Velocity (in/sec)');
        xlabel('Displacement (in)');       
    else
        ylabel('Velocity (m/sec)');
        xlabel('Displacement (m)');        
    end
%
    figure(5);
    plot(x,a);
    grid on;    
    if(iunits==1)
        ylabel('Accel (G)');
        xlabel('Displacement (in)');       
    else
        ylabel('Accel (m/sec^2)');
        xlabel('Displacement (m)');        
    end
%
end