%
disp(' ');
disp(' cant_beam_arbit_force.m   ver 1.3   February 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 a ramp');
disp(' invariant digital recursive filtering relationship. ');
disp(' ');
%
close all;
%
clear beta;
clear bL;
clear omegan;
clear D;
clear H;
clear HM;
clear freq;
clear damp;
clear cc;
clear df;
clear t;
clear f;
clear accel;
clear force;
clear ncontinue;
clear YL;
%
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();
NT=n;
%
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=12*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
%
sz=size(t);
if(sz(2)>sz(1))
    t=t';
end
sz=size(force);
if(sz(2)>sz(1))
    force=force';
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')
%%%%%%%%%%%%%%%%%%%%%%%%
%
YL=zeros(nmodes,1);
%
a1=zeros(nmodes,1);
a2=zeros(nmodes,1);
%
df1=zeros(nmodes,1);
df2=zeros(nmodes,1);
df3=zeros(nmodes,1);
%
vf1=zeros(nmodes,1);
vf2=zeros(nmodes,1);
vf3=zeros(nmodes,1);
%
af1=zeros(nmodes,1);
af2=zeros(nmodes,1);
af3=zeros(nmodes,1);
%
mass=1;
%
for j=1:nmodes
%
    arg=beta(j)*L;
    YL(j)=QS*Y(D(j),arg);
%
    omegad=omegan(j)*sqrt(1.-(damp^2));
    domegan=damp*omegan(j);
%    
    cosd=cos(omegad*dt);
    sind=sin(omegad*dt); 
%
    domegadt=domegan*dt;
%
    eee1=exp(-domegadt);
    eee2=exp(-2.*domegadt);
%
    ecosd=eee1*cosd;
    esind=eee1*sind; 
%
    a1(j)=2.*ecosd;
    a2(j)=-eee2;
%
    omeganT=omegan(j)*dt;
    phi=(2*damp)^2-1;
    DD1=(omegan(j)/omegad)*phi;
    DD2=2*DD1;
%    
    df1(j)=2*damp*(ecosd-1) +DD1*esind +omeganT;
    df2(j)=-2*omeganT*ecosd +2*damp*(1-eee2) -DD2*esind;
    df3(j)=(2*damp+omeganT)*eee2 +(DD1*esind-2*damp*ecosd);
%     
    VV1=(damp*omegan(j)/omegad);
%    
    vf1(j)=(-ecosd+VV1*esind)+1;
    vf2(j)=eee2-2*VV1*esind-1;
    vf3(j)=ecosd+VV1*esind-eee2;
%
    MD=(mass*omegan(j)^3*dt);
    df1(j)=df1(j)/MD;
    df2(j)=df2(j)/MD;
    df3(j)=df3(j)/MD;
%
    VD=(mass*omegan(j)^2*dt);
    vf1(j)=vf1(j)/VD;
    vf2(j)=vf2(j)/VD;
    vf3(j)=vf3(j)/VD;
%
    af1(j)=esind/(mass*omegad*dt);
    af2(j)=-2*af1(j);
    af3(j)=af1(j);
%
end    
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
nx=zeros(NT,nmodes);
nv=zeros(NT,nmodes);
na=zeros(NT,nmodes);
%
xL=zeros(NT,1);
vL=zeros(NT,1);
aL=zeros(NT,1);
%
for(j=1:nmodes)
%
    nodal_force=YL(j)*force;
%        
%  displacement
%
    d_forward=[   df1(j),  df2(j), df3(j) ];
    d_back   =[     1, -a1(j), -a2(j) ];
    d_resp=filter(d_forward,d_back,nodal_force);
%    
%  velocity
%
    v_forward=[   vf1(j),  vf2(j), vf3(j) ];
    v_back   =[     1, -a1(j), -a2(j) ];
    v_resp=filter(v_forward,v_back,nodal_force);
%    
%  acceleration
%   
    a_forward=[   af1(j),  af2(j), af3(j) ];
    a_back   =[     1, -a1(j), -a2(j) ]; 
    a_resp=filter(a_forward,a_back,nodal_force);
%
    nx(:,j)=d_resp;  % displacement
    nv(:,j)=v_resp;  % velocity
    na(:,j)=a_resp;  % acceleration  
%
end
%    
xL= (YL'*nx')';
vL= (YL'*nv')';
aL= (YL'*na')';
%    
aL=aL/386;
%
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