%
disp(' ');
disp(' cant_beam_arbit_force_RK4.m   ver 1.2   Febuary 6, 2012 ');
disp(' ');
disp(' Transverse vibration of a fixed-free beam, ');
disp(' subjected to an arbitrary force at free end.     ');
disp(' ');
disp(' A modal transient solution is performed using the ');
disp(' Runge-Kutta Fourth order method. ');
disp(' ');
%
disp(' Warning: the RK4 method may gave unstable results ');
disp(' for systems that are too stiff. ');
disp('  ');
%
close all;
clear all;
%
clear nx;
clear nv;
clear ny;
clear na;
%
clear beta;
clear bL;
clear omegan;
clear D;
clear H;
clear HM;
clear freq;
clear damp;
clear cc;
clear kk;
clear df;
clear t;
clear f;
clear tt;
clear ff;
clear accel;
clear force;
clear ncontinue;
%
tpi=2*pi;
%
disp(' The input time history must have two columns: time(sec) & force(lbf)')
disp(' ')
%
[t,force,dt,sr,tmx,tmi,n,ncontinue]=enter_time_history();
%
if(ncontinue==1)
%
iu=1;
%
nmodes=4;
%
[E,I,rho,c,L,A]=geometry_materials();
%
damp=input(' Enter the uniform modal damping ');
%
QS=1/sqrt(rho*A*L);
%
beta(1)=1.87510;
beta(2)=4.69409;
beta(3)=5*pi/2;
beta(4)=7*pi/2;
%
  beta=beta/L;
    bL=beta*L;
%
omegan=beta.*beta;
%
omegan=omegan*sqrt(E*I/(rho*A));
%
max_fn=max(omegan)/tpi;
sra=100*max_fn;
%
out1=sprintf(' sr=%8.4g  sra=%10.5g  ',sr,sra);
disp(out1);
%
disp(' ');
if(sra>sr)
%
    disp(' interpolate ');
%
    dt=1/sra;
%
    nn=round((tmx-tmi)/dt);
    tt=linspace(tmi,tmx,n);
%
% interpolate force
%
    ff = interp1(t,force,tt);
%
%
    clear t
    clear force
%
    t=tt;
    force=ff;
%
end
%
%
clear a;
D=zeros(nmodes,1);
%
for i=1:nmodes
    a=bL(i);
    D(i)=( cos(a)+cosh(a))/(sin(a)+sinh(a) );
end
%
disp(' ');
disp(' Natural Frequency (Hz) ');
for k=1:nmodes
    out1=sprintf(' %8.4g ',omegan(k)/(2*pi));
    disp(out1);
end
%
%%%%%%%%%%%%%%%%%%%%%%%%
  Y = inline('(cosh(arg)-cos(arg))-Di*(sinh(arg)-sin(arg))','Di','arg')
DDY = inline('(cosh(arg)+cos(arg))-Di*(sinh(arg)+sin(arg))','Di','arg')
%%%%%%%%%%%%%%%%%%%%%%%%
%
for j=1:nmodes
%
    arg=beta(j)*L;
    YL(j)=QS*Y(D(j),arg);
%    
end    
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
ndof=nmodes;
%
NT=n;
%
nx=zeros(NT,ndof);
ny=zeros(NT,ndof);
%
X1=zeros(ndof,1);
X2=zeros(ndof,1);
X3=zeros(ndof,1);
X4=zeros(ndof,1);
%
Y1=zeros(ndof,1);
Y2=zeros(ndof,1);
Y3=zeros(ndof,1);
Y4=zeros(ndof,1);
%
cc=zeros(ndof,1);
kk=zeros(ndof,1);
%
xx=zeros(ndof,1);
yy=zeros(ndof,1);
%
nFA=zeros(ndof,1);
nFB=zeros(ndof,1);
%
FF1=zeros(ndof,1);
FF2=zeros(ndof,1);
FF3=zeros(ndof,1);
FF4=zeros(ndof,1);
%
for(j=1:ndof)
    cc(j)=2*damp*omegan(j);
    kk(j)=(omegan(j))^2;
end    
%
progressbar % Create figure and set starting time 
%
sz=size(YL);
if(sz(2)>sz(1))
    YL=YL';
end
%
nFA=YL*force(i);
%
for(i=1:(NT-1))
% 
        progressbar(i/NT) % Update figure
%	    
        h=t(i+1)-t(i);
%
	    hd2=h/2;
        hd6=h/6.;
%
        for(j=1:ndof)
%
            nFB(j)=YL(j)*force(i+1);
%
            FF1(j)= nFA(j);
            FF2(j)= (nFB(j) + nFA(j))/2.;
            FF3(j)= FF2(j);
            FF4(j)= nFB(j);
%
            xx=nx(i,j)';
            yy=ny(i,j)';
%
            X1= xx;
            Y1= yy;
            F1= (FF1(j)-cc(j)*Y1-kk(j)*X1);
%
            X2= X1+Y1*hd2;
            Y2= Y1+F1*hd2;
            F2= (FF2(j)-cc(j)*Y2-kk(j)*X2);
%
            X3= X1+Y2*hd2;
            Y3= Y1+F2*hd2;
            F3= (FF3(j)-cc(j)*Y3-kk(j)*X3);
%
            X4= X1+Y3*h;
            Y4= Y1+F3*h;
            F4= (FF4(j)-cc(j)*Y4-kk(j)*X4);
%
            nx(i+1,j)= (X1+(h/6.)*( Y1 + 2*Y2 + 2*Y3 + Y4 ))';    % displacement
            ny(i+1,j)= (Y1+(h/6.)*( F1 + 2*F2 + 2*F3 + F4 ))';    % velocity
            na(i+1,j)= FF4(j) -cc(j)*ny(i+1,j) -kk(j)*nx(i+1,j); % acceleration  
        end     
%
        nFA=nFB;
%
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
progressbar(1);
%
clear x;
clear y;
clear a;
%
xL=zeros(NT,1);
vL=zeros(NT,1);
aL=zeros(NT,1);
%
progressbar
%
for(i=1:NT)
    progressbar(i/NT) % Update figure    
    xL(i)=(YL'*(nx(i,:))');
    vL(i)=(YL'*(ny(i,:))');
    aL(i)=(YL'*(na(i,:))');    
end
progressbar(1);
%
aL=aL/386;
%
%
sz=size(t);
if(sz(2)>sz(1))
    t=t';
end
sz=size(force);
if(sz(2)>sz(1))
    force=force';
end
%
ijk=1;
while(ijk==1)
%
    disp(' ');
    tmin=input(' Enter plot start time(sec) ');
    tmax=input(' Enter plot end time(sec) ');
%
    figure(1);
    plot(t,xL);
    axis([tmin,tmax,min(ylim),max(ylim)]);
    grid on;
    title('Displacement at x=L');
    ylabel('Disp (in)');
    xlabel('Time (sec)');
%
    figure(2);
    plot(t,vL);
    axis([tmin,tmax,min(ylim),max(ylim)]);    
    grid on;
    title('Velocity at x=L');
    ylabel('Vel (in/sec)');
    xlabel('Time (sec)');
%
    figure(3);
    plot(t,aL);
    axis([tmin,tmax,min(ylim),max(ylim)]);   
    grid on;
    title('Acceleration at x=L');
    ylabel('Accel (G)');
    xlabel('Time (sec)');
%
    disp(' ');
    disp(' Replot with different limits?');
    disp('  1=yes  2=no ');
    ijk=input(' ');
end
%
clear FD;
clear FV;
clear FA;
%
FD=[t force xL];
FV=[t force vL];
FA=[t force aL];
%
mx=max(abs(xL));
mv=max(abs(vL));
ma=max(abs(aL));
%
disp(' ');
disp(' Maximum Response Levels ');
out1=sprintf(' Displacement = %8.4g in',mx);
out2=sprintf('     Velocity = %8.4g in/sec',mv);
out3=sprintf(' Acceleration = %8.4g G',ma);
%
disp(out1);
disp(out2);
disp(out3);
%
end