%   
disp(' ');
disp(' SDOF_base_animation.m   ver 1.1   November 6, 2006 ');
disp(' by Tom Irvine  Email: tomirvine@aol.com ');
disp(' ');
disp(' This program generates an animation of a single-degree-of-freedom ');
disp(' system subjected to sinusoidal base excitation. ');
%
clear all 
%
    disp(' ');
    fn=input(' Enter natural frequency (Hz) ');
    disp(' ');
    f =input(' Enter excitation frequency (Hz) ');	
    disp(' ');
    damp =input(' Enter the damping ratio (typically 0.05) ');	   
%
    if(damp < 0.01)
        disp(' ');
        disp(' Damping changed to mininimum value: 0.01 ');
        damp=0.01;
    end
%
    rho=f/fn;
%
    aa=1.+(2.*damp*rho)^2;
    bb=(1.-rho^2)^2+(2.*damp*rho)^2;
    trans=sqrt(aa/bb);
%    
    aa=2.*damp*rho^3;
    bb=1.+rho^2*(4.*damp^2-1);
    phase=atan(aa/bb);
%
    A=1;
% 
    Tn=1/fn;
    dt1=Tn/32;
    dt2=(1/f)/8;
    dt=dt1;
    if(dt2<dt1)
        dt=dt2;
    end
%
    Tmax=Tn;
    if(Tmax<(1/f))
        Tmax=(1/f);
    end    
%
nframes=4*round(Tmax/dt);
k=1;
while(1)
%
    t=(k-1)*dt;
    th(k)=t;
    yb(k)=A*sin(2*pi*t*f);
%
    if( k>1000 | k>nframes)
        break;
    end
    k=k+1;
end
kend=k;
maxyb=max(yb);
for(k=1:kend)
    t=(k-1)*dt;
    yy(k)=A*trans*sin(2*pi*t*f-phase);
end
ymax=max(yy);
offset=ymax*1.5;
if(offset<4)
    offset=4;
end
y=yy+offset;
yup=max(y)*1.3;
%
bbb=0.5;
xbase=[-0.5,-0.5,0.5,0.5]; % plotting coordinates of base
%
xmass=xbase; % plotting coordinates of mass
%
iflag=1;
while(iflag==1)
%
    ccc=(yup-min(y))/5;
    bbb=ccc/8;
%
    if(bbb<0.5)
        bbb=0.5;
        ccc=bbb*8;
        yup=5*ccc+min(y);
    end
    if(yup< (0.9*max(y)+ccc))
        yup=yup*1.1
    else
        iflag=0;
    end
end
offset=(yup-ccc)/2;
y=yy+offset;
%
ybase=[-bbb,bbb,bbb,-bbb];
ymass=[0,ccc,ccc,0];
%
for k=1:kend
%
    distance=y(k)-(yb(k)+bbb);
%
    x1(k)=0.;
    x2(k)=0.;
    x3(k)=-0.25;
    x4(k)=+0.25;
    x5(k)=-0.25;
    x6(k)=+0.25;
    x7(k)=-0.25;
    x8(k)=+0.25;
    x9(k)=0.;
    x10(k)=0.;        
    dyf=distance/8;
    dyh=dyf/2;
%
    y1(k)=yb(k);
    y2(k)=y1(k)+dyf;
    y3(k)=y2(k)+dyh;
    y4(k)=y3(k)+dyf;
    y5(k)=y4(k)+dyf;
    y6(k)=y5(k)+dyf;
    y7(k)=y6(k)+dyf;
    y8(k)=y7(k)+dyf;
    y9(k)=y8(k)+dyh;
    y10(k)=y9(k)+dyf; 
end
%
y1=y1+bbb;
y2=y2+bbb;
y3=y3+bbb;
y4=y4+bbb;
y5=y5+bbb;
y6=y6+bbb;
y7=y7+bbb;
y8=y8+bbb;
y9=y9+bbb;
y10=y10+bbb;
%
L0=3;
Hf=figure('Units','normalized','Position',[.1,.1,.8,.4]);
%
Ha_f2a2=subplot(1,4,1);
%
Hp_patch_base=patch(xbase,ybase,[0 0.2 0.2]);
Hp_patch_mass=patch(xmass,ymass,[0 0.5 0.5]);
%
k=1;
hspring=line([x1(k) x2(k) x3(k) x4(k) x5(k) x6(k) x7(k) x8(k) x9(k) x10(k)],[y1(k) y2(k) y3(k) y4(k) y5(k) y6(k) y7(k) y8(k) y9(k) y10(k)],'color','k','linewidth',2);
%
axis([-2,2,-2,yup]),grid on 
%
Ha_f2a1=subplot(1,4,2:4);
%
Hls_f2plot1=plot(th(1),yb(1),th(1),y(1));axis([min(th),th(end),-2,yup]),grid on,...
   xlabel('t (sec)')  
%
clear figure;
refresh;
figure(1);
clear F;
set(gcf,'DoubleBuffer','on');
out1=sprintf(' Natural Freq =%7.4g Hz   Excitation Freq =%7.4g Hz   Damp Ratio =%7.4g',fn,f,damp);
for k=1:kend    
%
     plot(th(1:k),yb(1:k),th(1:k),(y(1:k)+ccc/2));axis([min(th),th(end),-2,yup]),grid on,...
   xlabel('time (sec)'),title(out1)  
%
    set(hspring,'XData',[x1(k) x2(k) x3(k) x4(k) x5(k) x6(k) x7(k) x8(k) x9(k) x10(k)],'YData',[y1(k) y2(k) y3(k) y4(k) y5(k) y6(k) y7(k) y8(k) y9(k) y10(k)]);
    set(Hp_patch_base,'YData',ybase+yb(k)); 
    set(Hp_patch_mass,'YData',ymass+y(k));    
    if(k==1)
        pause(3);  
    else
        pause(0.1);        
    end
%    F(k)=getframe;
    drawnow;
end    
pause(10);  