disp(' ');
disp(' RK4_force.m   ver 1.0  October 19, 2011 ');
disp(' by Tom Irvine   Email: tomirvine@aol.com');
%
disp(' ')
disp(' This program calculates the response of a single-degree-of-freedom')
disp(' system to an arbitrary force time history using the ')
disp(' Runge-Kutta fourth order method.')
disp(' ')
disp(' The input time history must have two columns: time(sec) & force(lbf)')
disp(' ')
%
close all;
clear t;
clear y;
clear x;
clear z;
clear a;
clear accel;
clear length;
clear omegan;
%
fig_num=1;
%
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);
a=THM(:,2);
num=length(t);
%
figure(fig_num);
fig_num=fig_num+1;
plot(t,a);
grid on;
xlabel('Time(sec)');
ylabel('Force(lbf)');
title('Applied Force ');
%
mass=input(' Enter the mass (lbm)  ');
mass=mass/386;
%
disp(' ')
fn=input(' Enter the natural frequency (Hz)  ');
%
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
%
omegan=2*pi*fn;
omegad=omegan*sqrt(1-(damp^2));
domegan=damp*omegan;
%
two_domegan=2*domegan;
%
omn2=omegan^2;
%
stiff=mass*omn2;
%
ccc=two_domegan*mass;
%
out1=sprintf('\n damping coefficient = %8.4g lbf/(in/sec) ',ccc);
disp(out1);
out1=sprintf('           stiffness = %8.4g lbf/in ',stiff);
disp(out1);
%
x=zeros(num,1);
y=zeros(num,1);
accel=zeros(num,1);
%
disp(' ');
d0=input(' Enter the initial displacement (in) ');
disp(' ');
v0=input(' Enter the initial velocity (in/sec) ');
%
x(1)=d0;
y(1)=v0;
%
a=a/mass;
%
for(i=1:(num-1))
% 
	    h=t(i+1)-t(i);
%
	    hd2=h/2;
%
	    a1=           a(i);
	    a2= (a(i+1) + a(i))/2.;
	    a3= a2;
	    a4=  a(i+1);
%
        T1= t(i);
        X1= x(i);
        Y1= y(i);
        F1=a1-two_domegan*Y1-omn2*X1;
%
        T2=t(i)+hd2;
        X2= x(i)+Y1*hd2;
        Y2= y(i)+F1*hd2;
        F2= a2-two_domegan*Y2-omn2*X2;
%
        T3=t(i)+hd2;
        X3= x(i)+Y2*hd2;
        Y3= y(i)+F2*hd2;
        F3= a3-two_domegan*Y3-omn2*X3;
%
        T4=t(i)+h;
        X4= x(i)+Y3*h;
        Y4= y(i)+F3*h;
        F4= a4-two_domegan*Y4-omn2*X4;
%
        x(i+1)=x(i)+(h/6.)*( Y1 + 2*Y2 + 2*Y3 + Y4 );   % displacement
        y(i+1)=y(i)+(h/6.)*( F1 + 2*F2 + 2*F3 + F4 );   % velocity
%
        accel(i+1)=-two_domegan*Y4-omn2*X4;
%
end
%
%
clear accel_resp;
clear vel_resp;
clear disp_resp;
clear tf_resp;
%
accel_resp=accel/386;
vel_resp=y;
disp_resp=x;
tf_resp=ccc*vel_resp + stiff*disp_resp;
%
%
disp(' ')
disp(' ')
disp(' Acceleration Response ')
mu=mean(accel_resp);
sd=std(accel_resp);
arms=sqrt(sd^2+mu^2);
peak=max(abs(accel_resp));
%
clear aar;
aar=accel_resp-mu;
aar=aar.^4;
n=num;
kurtosis=sum(aar)/(n*sd^4);
%
out10=sprintf('\n      maximum = %10.2f G',max(accel_resp));
  out11=sprintf('      minimum = %10.2f G',min(accel_resp));
  out12=sprintf('      overall = %10.2f GRMS',arms);
out13=sprintf('\n         mean = %10.2f G',mu);
  out14=sprintf('      std dev = %10.2f G',sd);
  out15=sprintf(' crest factor = %10.2f',peak/sd);
  out16=sprintf('     kurtosis = %10.2f',kurtosis);  
%
disp(out10)
disp(out11)
disp(out12)
disp(out13)
disp(out14)
disp(out15)
disp(out16)
%
disp(' ')
disp(' Velocity Response ')
mu=mean(vel_resp);
sd=std(vel_resp);
vel_rms=sqrt(sd^2+mu^2);
out10=sprintf('\n maximum = %10.3g in/sec',max(vel_resp));
out11=sprintf(' minimum = %10.3g in/sec',min(vel_resp));
out12=sprintf(' overall = %10.3g in/sec RMS',vel_rms);
disp(out10)
disp(out11)
disp(out12)
%
disp(' ')
disp(' Displacement Response ')
mu=mean(disp_resp);
sd=std(disp_resp);
rd_rms=sqrt(sd^2+mu^2);
out10=sprintf('\n maximum = %10.3g inch',max(disp_resp));
out11=sprintf(' minimum = %10.3g inch',min(disp_resp));
out12=sprintf(' overall = %10.3g inch RMS',rd_rms);
disp(out10)
disp(out11)
disp(out12)
%
nn=length(accel_resp);
%
tmin=t(1);
tmax=max(t);
%
tt=linspace(tmin,tmax,nn);
ymin=1.3*min(accel_resp);
ymax=1.3*max(accel_resp);
if(abs(ymin)>abs(ymax))
    ymax=abs(ymin); 
end
if(abs(ymax)>abs(ymin))
    ymin=-abs(ymax); 
end
%    
%  Output 
%
disp(' ')
disp('Output Acceleration Response Text File?')
choice=input(' 1=yes  2=no  ' );
disp(' ')
%
if choice == 1 
%
      [writefname, writepname] = uiputfile('*.dat','Save acceleration data as');
	  writepfname = fullfile(writepname, writefname);
	  writedata = [tt' accel_resp'];
	  fid = fopen(writepfname,'w');
	  fprintf(fid,'  %g     %g\n',writedata');
	  fclose(fid);
%    
%   disp(' Enter output filename ');
% TH_filename = input(' ','s');
%
%  fid = fopen(TH_filename,'w');
%   for i=1:nn
%        fprintf(fid,'%12.5f %10.3f \n',tt(i),accel_resp(i));
%   end
%   fclose(fid);
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp('Plot Acceleration Time History?')
choice=input(' 1=yes  2=no  ' );
disp(' ')
if choice == 1 
    figure(fig_num);
    fig_num=fig_num+1;
    disp(' Plotting output..... ')
    plot(tt,accel_resp);
    axis([tmin,tmax,ymin,ymax]);
    ylabel('Accel (G)');
    xlabel('Time (sec)');
    Q=1./(2.*damp);
    out5 = sprintf(' SDOF Acceleration Response   fn=%g Hz   Q=%g ',fn,Q);
    title(out5);
    grid;
    zoom on;
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp('Plot Velocity Time History?')
choice=input(' 1=yes  2=no  ' );
disp(' ')
if choice == 1 
    figure(fig_num);
    fig_num=fig_num+1;
    disp(' Plotting output..... ')
    plot(tt,vel_resp);
    ymin=1.3*min(vel_resp);
    ymax=1.3*max(vel_resp);
    ymin=ymin;
    ymax=ymax;
    if(abs(ymin)>abs(ymax))
        ymax=abs(ymin); 
    end
    if(abs(ymax)>abs(ymin))
        ymin=-abs(ymax); 
    end  
    axis([tmin,tmax,ymin,ymax]);
    ylabel('Vel (in/sec)');
    xlabel('Time (sec)');
    Q=1./(2.*damp);
    out5 = sprintf(' SDOF Velocity Response   fn=%g Hz   Q=%g ',fn,Q);
    title(out5);
    grid;
    zoom on;
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ')
disp('Plot Displacement Time History?')
choice=input(' 1=yes  2=no  ' );
disp(' ')
if choice == 1
    figure(fig_num);
    fig_num=fig_num+1;
    disp(' Plotting output..... ')
    plot(tt,disp_resp);
    ymin=1.3*min(disp_resp);
    ymax=1.3*max(disp_resp);
    ymin=ymin;
    ymax=ymax;
    if(abs(ymin)>abs(ymax))
        ymax=abs(ymin); 
    end
    if(abs(ymax)>abs(ymin))
        ymin=-abs(ymax); 
    end    
    axis([tmin,tmax,ymin,ymax]);    
    ylabel('Disp (inch)');
    xlabel('Time (sec)');
    Q=1./(2.*damp);
    out5 = sprintf(' SDOF Displacement Response   fn=%g Hz   Q=%g ',fn,Q);
    title(out5);
    grid;
    zoom on;
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ')
disp('Plot Transmitted Force Time History?')
choice=input(' 1=yes  2=no  ' );
disp(' ')
if choice == 1
    figure(fig_num);
    fig_num=fig_num+1;
    disp(' Plotting output..... ')
    plot(tt,tf_resp);
    ymin=1.3*min(tf_resp);
    ymax=1.3*max(tf_resp);
    ymin=ymin;
    ymax=ymax;
    if(abs(ymin)>abs(ymax))
        ymax=abs(ymin); 
    end
    if(abs(ymax)>abs(ymin))
        ymin=-abs(ymax); 
    end    
    axis([tmin,tmax,ymin,ymax]);    
    ylabel('Force (lbf)');
    xlabel('Time (sec)');
    Q=1./(2.*damp);
    out5 = sprintf(' Transmitted Force   fn=%g Hz   Q=%g ',fn,Q);
    title(out5);
    grid;
    zoom on;
end