function f=halfpipe(aviname)
% halfpipe movie

if nargin==0, aviname='halfpipe2.avi'; end;
i=1;
for t=0:0.05:10
    halfpipet(t);
    f(i)=getframe(gca);
    i=i+1;
end;
if nargout==0,
    movie2avi(f, aviname, 'quality', 75);
end;


    function halfpipet(t)
        % plot halfpipe and rolling sphere at time t
        
        if nargin==0, t=6; end;
        r=25;
        rk=3;
        m=0.1;
        ts=0.05;
        [x,y]=meshgrid(-40:1:40,0:2:100);
        z=m*y-real(sqrt(r*r-x.*x));
        h=surf(x,y,z,zeros(size(z)));
        light;
        view(-10,30);
        axis equal; axis off
        set(h,'edgealpha',0.3)
        T=[0:ts:t];
        Y=100-T.*T;
        y1=100-t*t;
        X=10*cos(pi/2*T);
        x1=10*cos(pi/2*t);
        z1=m*y1-sqrt(r*r-x1*x1)+rk;
        hold on;
        plot3(X, Y, m*Y-sqrt(r*r-X.*X),'-k','Linewidth', 2);
        [a,b,c]=sphere;
        k=surf(rk*a+x1,rk*b+y1,rk*c+z1,c.*a.*b);
        set (k, 'edgealpha',0.2);
        if t>0
            for d=ts:ts:t
                dy=d*d-(d-ts)*(d-ts);
                dx=10*(cos(pi/2*d)-cos(pi/2*(d-ts)));
                rotate (k, [dy,dx,0], sqrt(dx*dx+dy*dy)/(pi*rk)*180, [x1,y1,z1]);
            end;
        end;
        hold off;
        set(gcf,'Position',[100 50 460 412]);
        set(gca,'Position',[0 0 1 1]);
        set(gcf,'Color','w');
    end

end