disp(' ')
disp(' sdof_free_Dahl.m ')
disp(' ver 1.0  April 2, 2012 ')
disp(' by Tom Irvine  Email: tomirvine@aol.com ')
disp(' ')
disp(' This program calculates the free vibration of a  ')
disp(' single-degree-of-freedom system with Dahl damping ');
disp(' subjected to initial conditions. ');
disp(' ')
disp(' The solution method is Newmark-beta, which is a numerical method.')
disp(' ')
%
close all;
%
clear omegan;
clear omegad;
clear fn;
clear damp;
clear N;
clear t;
clear x;
clear w;
clear m;
clear k;
clear M;
clear C;
clear K;
clear U;
clear Ud;
clear Udd;
%
disp(' ');
disp(' Select units: ');
disp('   1=English  2=metric ');
iu=input(' ');
%
if(iu==1)
%
    disp(' ')
    v0 = input(' input initial velocity (in/sec) ');
    disp(' ')
    d0 = input(' input initial displacement (in) ');   
%
    disp(' ')
    m = input(' input mass (lbm) ');
    m = m/386;
%
    disp(' ')
    k = input(' input stiffness (lbf/in) ');
%
else
%
    disp(' ')
    v0 = input(' input initial velocity (meters/sec) ');
    disp(' ')
    d0 = input(' input initial displacement (meters) ');   
%
    disp(' ')
    m = input(' input mass (kg) ');
%
    disp(' ')
    k = input(' input stiffness (N/m) ');
%
end
%
disp(' ');
disp(' Enter coefficient of sliding friction mu ');
mu=input(' ');
%
if(iu==1)
    N=m*386;
else
    N=m*9.81;
end
%
Fc=mu*N;
%
disp(' ');
if(iu==1)
    disp(' Enter coefficient sigma-zero (lbf/in)');
else
    disp(' Enter coefficient sigma-zero (N/m)');   
end
sigma=input(' ');
%
omegan=sqrt(k/m);    
%
%   State Space Solution
%
fn=omegan/(2*pi); 
%
out1=sprintf('\n  fn=%8.4g Hz ',fn);
disp(out1);
%
period=1/fn;
%
dt=period/256;
sr=1/dt;
dur=30*period;
N=dur/dt;
num=N;
%
% out6 = sprintf('\n dt = %f sec    dur = %f sec ',dt,dur);
% disp(out6)
%
t=linspace(0,dur,N);
t=t';
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
[a0,a1,a2,a3,a4,a5,a6,a7]=Newmark_coefficients(dt);
%
M=m;
C=0;
K=k;
%
KH=K+a0*M+a1*C;
%
clear length;
%
U=zeros(num,1);
Ud=zeros(num,1);
Udd=zeros(num,1);
tf=zeros(num,1);
%
U(1)=d0;
Ud(1)=v0;
Udd(1)=0;
%
alpha=sigma/Fc;
%
for(i=2:num);
%
        F= -Fc*(1-exp(-alpha*abs(U(i-1))))*sign(Ud(i-1));
%
        FH=F+M*(a0*U(i-1)+a2*Ud(i-1)+a3*Udd(i-1))+C*(a1*U(i-1)+a4*Ud(i-1)+a5*Udd(i-1));
        Un=FH/KH;
        Uddn=a0*(Un-U(i-1))-a2*Ud(i-1)-a3*Udd(i-1);
        Udn=Ud(i-1)+a6*Udd(i-1)+a7*Uddn;
%        
        U(i)=Un;
        Ud(i)=Udn;
        Udd(i)=Uddn;
%
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
dmax=max(U);
vmax=max(Ud);
%
if(iu==1)
    out6 = sprintf('\n max disp = %9.3g in   max vel = %9.3g in/sec',dmax,vmax);
else
    out6 = sprintf('\n max disp = %9.3g meters   max vel = %9.3g m/sec',dmax,vmax);
end
disp(out6);
%
figure(1);
plot(t,U);
title('Displacement');
xlabel('Time (sec)'); 
if(iu==1)
       ylabel('Disp (inch)');
else
       ylabel('Disp (meters)');
end
grid on;    
%
figure(2);
plot(t,Ud);
title('Velocity');
xlabel('Time (sec)'); 
if(iu==1)
        ylabel('Vel (in/sec)');
else
        ylabel('Vel (meters/sec)');        
end
grid on;   
%