disp(' ');
disp(' arbit_newmark_force.m,  ver 1.1, December 27, 2011 ');
disp(' by Tom Irvine,  Email: tomirvine@aol.com ');
disp(' This program calculates the response of an SDOF system ');
disp(' to an arbitrary force via the Newmark-beta method.  ');
%
%
clear t;
clear f;
clear U;
clear Ud;
clear Udd;
clear tf;
%
tpi=2*pi;
%
disp(' ');
disp(' Select Units ');
disp('   1=English  2=metric ');
iunits=input(' ');
%
disp(' ');
disp(' The applied force should have a constant dt. ');
disp(' The base input file must contain two columns: ');
%
if(iunits==1)
    disp(' time(sec)  force(lbf)   ');
else
    disp(' time(sec)  force(N)    ');
end
%
disp(' ');
disp(' Select file input method ');
disp('   1=external ASCII file ');
disp('   2=file preloaded into Matlab ');
disp('   3=Excel file ');
file_choice = input('');
%
if(file_choice==1)
        [filename, pathname] = uigetfile('*.*');
        filename = fullfile(pathname, filename);
        fid = fopen(filename,'r');
        THM = fscanf(fid,'%g %g',[2 inf]);
        THM=THM';
end
if(file_choice==2)
        THM = input(' Enter the matrix name:  ');
end
if(file_choice==3)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        THM = xlsread(xfile);
%         
end
%
t=THM(:,1);
f=THM(:,2);
%
%
tmx=max(t);
tmi=min(t);
n = length(f);
dt=(tmx-tmi)/(n-1);
sr=1./dt;
%
disp(' ')
disp(' Time Step ');
dtmin=min(diff(t));
dtmax=max(diff(t));
%
out4 = sprintf(' dtmin  = %8.4g sec  ',dtmin);
out5 = sprintf(' dt     = %8.4g sec  ',dt);
out6 = sprintf(' dtmax  = %8.4g sec  ',dtmax);
disp(out4)
disp(out5)
disp(out6)
%
disp(' ')
disp(' Sample Rate ');
out4 = sprintf(' srmin  = %8.4g samples/sec  ',1/dtmax);
out5 = sprintf(' sr     = %8.4g samples/sec  ',1/dt);
out6 = sprintf(' srmax  = %8.4g samples/sec  \n',1/dtmin);
disp(out4)
disp(out5)
disp(out6)
%
%**************************************************************************************************
%
disp(' ');
if(iunits==1)
    disp(' Enter the mass (lbm) ');	
else
    disp(' Enter the mass (kg) ');
end
%
mass=input(' ');
%
if(iunits==1)
    mass=mass/386.;
end
disp(' ');
disp(' Enter natural frequency (Hz) ');
fn=input(' ');
%
omegan=tpi*fn;
stiff=(omegan^2)*mass;
%
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);
domegan=damp*omegan;
%
[a0,a1,a2,a3,a4,a5,a6,a7]=Newmark_coefficients(dt)
%
M=mass;
K=stiff;
C=2*damp*omegan*mass;
%
KH=K+a0*M+a1*C;
%
U=0.;
Ud=0.;
Udd=0.;
tf=0.;
%
clear length;
num=n;
%
U=zeros(num,1);
Ud=zeros(num,1);
Udd=zeros(num,1);
tf=zeros(num,1);
%
for(i=2:num);
%
	    F=f(i);
%
        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
%
tf=K*U + C*Ud;
%
%
if(iunits==2)
		out1=sprintf('\n\n maximum disp = %12.3g mm ',max(U));
		out2=sprintf(' minimum disp = %12.3g mm ',min(U));
else
		out1=sprintf('\n maximum disp = %12.3g inch ',max(U));
		out2=sprintf(' minimum disp = %12.3g inch ',min(U));
end
%
if(iunits==2)
		out3=sprintf('\n maximum velocity = %12.3g mm/sec ',max(Ud));
		out4=sprintf(' minimum velocity = %12.3g mm/sec ',min(Ud));
else
		out3=sprintf('\n maximum velocity = %12.3g inch/sec ',max(Ud));
		out4=sprintf(' minimum velocity = %12.3g inch/sec ',min(Ud));
end
%
if(iunits==2)
		out5=sprintf('\n maximum accel = %12.3g m/sec^2 ',max(Udd));
		out6=sprintf(' minimum accel = %12.3g m/sec^2 ',min(Udd));
else
		out5=sprintf('\n maximum accel = %12.3g G ',max(Udd)/386.);
		out6=sprintf(' minimum accel = %12.3g G ',min(Udd)/386.);
end
%
if(iunits==2)
		out7=sprintf('\n maximum transmitted force = %12.3g N ',max(tf));
		out8=sprintf(' minimum transmitted force  = %12.3g N ',min(tf));
else
		out7=sprintf('\n maximum transmitted force  = %12.3g lbf ',max(tf));
		out8=sprintf(' minimum transmitted force  = %12.3g lbf ',min(tf));
end
%   
if(iunits==2)
	out_d='Disp(m)';
    out_v='Vel(m/sec)';
    out_a='Acc(m/sec^2)';
    out_tf='Force(N)';
else
	out_d='Disp(in)';
    out_v='Vel(in/sec)';
    out_a='Acc(in/sec^2)';
    out_tf='Force(lbf)';
end
%
disp(out1);
disp(out2);
disp(out3);
disp(out4);
disp(out5);
disp(out6);
disp(out7);
disp(out8);
%
figure(1);
plot(t,U);
grid on;
title('Displacement');
ylabel(out_d);
xlabel('Time(sec)');
%
figure(2);
plot(t,Ud);
grid on;
title('Velocity');
ylabel(out_v);
xlabel('Time(sec)');
%
figure(3);
plot(t,Udd);
grid on;
title('Acceleration');
ylabel(out_a);
xlabel('Time(sec)');
%
figure(4);
plot(t,tf);
grid on;
title('Transmitted Force');
ylabel(out_tf);
xlabel('Time(sec)');