clear all;
t=[0:0.1:20];
n=length(t);
a=2.;
x1=linspace(-a,a,2000);
th=linspace(0,2*pi,100);
r=.1;
h=.75;
w=1.5;
rm=.02;
hr=.85;
crcx=w+r*cos(th);
crcy=h+r+r*sin(th);
mouthx=w+rm*cos(th);
mouthy=hr+rm*sin(th);
for i=1:n
   h1=(1+t(i:i))^(-.2)*(1-.9/(1+t(i:i))^.4*x1.^2).^(1/3);
   zer1=find(abs(x1)>(sqrt(10)/3*(1+t(i:i))^.2));
   h1(zer1)=0*h1(zer1);
   jmp=.5*(t(i)>3); 
   jmp=jmp*(mod(i,2*(1+fix(i^2/6/n)))==0);
   plot([1.5,1.5],[.5+jmp,h+jmp])
   hold on;
   if(t(i) <=3) 
      plot([1.45,1.48],[.87,.88])
      plot([1.52,1.55],[.88,.87])
      plot([1.45,1.5,1.55],[.83,.8,.83])
      plot(mouthx,mouthy+jmp)
   else
      plot([1.45,1.48],[.88+jmp,.88+jmp])
      plot([1.52,1.55],[.88+jmp,.88+jmp])
      plot([1.45,1.5,1.55],[.77+jmp,.8+jmp,.77+jmp],'r')
      plot(mouthx,mouthy+jmp,'r')
   end   
   plot([1.35,1.4,1.5,1.6,1.65],[0+jmp,0+jmp,.5+jmp,0+jmp,0+jmp])
   plot([1.35,1.4,1.5,1.6,1.65],[0+jmp,0+jmp,.5+jmp,0+jmp,0+jmp])
   plot([1.4,1.5,1.6],[.5+jmp,.7+jmp,.5+jmp])
   plot(crcx,crcy+jmp)
   fill(x1,h1,'r');
   plot([-a,a],[0,0],'g');
   axis([-a a -1 2] );
   hold off;
   pause(0.001);
end;