function tumbling_box
   % Angular momentum of tumbling box.
   % Cleve's Corner Blog, 8/10/2015. 
   % Gil Strang, Differential Equations and Linear Algebra, pp. 176-178.
   % Alar Toomre, MIT.
   % Click anywhere to provide initial condiitons.
   % You can also use the zoom and rotate3d controls in the toolbar.

   function ydot = momentum(t,y,y0)
      ydot = [ y(2)*y(3); -2*y(1)*y(3); y(1)*y(2)];
   end

   function [val,isterm,dir] = gstop(t,y,y0)
      % Event function for periodicity.
      % See C. Moler, Numerical Computing with MATLAB, Section 7.2.
      d = y - y0;
      v = momentum(t,y);
      val = d'*v;
      isterm = 1;
      dir = 1;
   end

   function my_btn_down(h)
      % Interpret mouse click as point on unit sphere.
      % Ref: Mike Garrity, MathWorks, email.
      ax = ancestor(h,'axes');
      cp = get(ax,'CurrentPoint');
      s = cp(2,:) - cp(1,:);
      s = s / norm(s);
      o = cp(1,:);
      r = 1;
      c = [0 0 0];
      q = dot(s,(o-c))^2 - dot(o-c,o-c) + r^2;
      if q >= 0 
         d1 = -dot(s,(o-c)) + sqrt(q);
         d2 = -dot(s,(o-c)) - sqrt(q);      
         pt1 = o + d1*s;
         pt2 = o + d2*s;
         v = get(ax,'CameraTarget') - get(ax,'CameraPosition');
         if dot(v,s) < 0
            y0 = pt1';
         else
            y0 = pt2';
         end
         [t,y] = ode45(@momentum,[0,100],y0,opt,y0);
         line(y(:,1),y(:,2),y(:,3),'color','black')
      end
   end

   opt = odeset('RelTol',1.e-6,'events',@gstop);
   clf
   [X,Y,Z] = sphere(36);
   h = surface(X,Y,Z);
   view(3)
   zoom(2)
   shading interp
   axis equal
   axis vis3d
   axis off
   set(gca,'xdir','rev','ydir','rev')
   % Three critical points
   dkblue = [0 0 .5];
   line(0,0,1.01,'marker','.','color',dkblue,'markersize',20)
   line(0,1.01,0,'marker','.','color',dkblue,'markersize',20)
   line(1.01,0,0,'marker','.','color',dkblue,'markersize',20)
   h.ButtonDownFcn = @(h,~) my_btn_down(h);
end
