function sticks
   % Chaotic Swinging Sticks
   % https://www.the-swinging-sticks.com/en
   % https://blogs.mathworks.com/cleve/2024/02/20/chaotic-swinging-sticks/
   % 02/20/2024
      
   function M = swingmass(~,u)      
      % Mass matrix
      c = cos(u(1)-u(2));
      M = [1 0 0 0; 0 1 0 0; 0 0 5 c; 0 0 c 1];
   end
      
   % ------------------------
      
   function f = swingrhs(~,u)
      % Driving force
      g = 1;
      s = sin(u(1)-u(2));
      f = [u(3); u(4); -g*sin(u(1))-s*u(4)^2; -g*sin(u(2))+s*u(3)^2];
   end
      
   % ------------------------
      
   function theta = swinginit(x,y)      
      r = norm([x,y]);
      if r > 2
         alpha = 0;
      else
         alpha = acos(r/2);
      end
      beta = atan2(y,x) + pi/2;
      theta = [beta+alpha; beta-alpha];
   end
      
   % ------------------------
      
   function status = swingplot(t,u,task)
      if ~isequal(task,'done')
      
         theta = u(1:2);
         x = sin(theta);
         y = -cos(theta);
         switch task
                       
            case 'init'
         
               % Initialize plot
         
               plt = plot( ...
                   [-x(1) 3*x(1) NaN 3*x(1)-x(2) 3*x(1)+2*x(2)], ...
                   [-y(1) 3*y(1) NaN 3*y(1)-y(2) 3*y(1)+2*y(2)], ...
                   '-','markersize',14,'linewidth',3);
               plt(2) = line(plt.XData(2),plt.YData(2), ...
                   'marker','.','markersize',24);
               line(0,0,'marker','.','markersize',24)
               axis(6*[-1 1 -1 1])
               axis square
               set(gca,'xtick','','ytick','')
               titl = title(sprintf('t = %8.1f',t));
               set(gcf,'userdata',{plt,titl})
               shg
               gif_frame('sticks_gif.gif',inf')

            case ''
             
               % Update plot
       
               fud = get(gcf,'userdata');
               [plt,titl] = deal(fud{:});

               set(plt(1), ...
                   'xdata',[-2*x(1) 3*x(1) NaN 3*x(1)-x(2) 3*x(1)+2*x(2)], ...
                   'ydata',[-2*y(1) 3*y(1) NaN 3*y(1)-y(2) 3*y(1)+2*y(2)])
               set(plt(2),...
                   'xdata',plt(1).XData(2), ...
                   'ydata',plt(1).YData(2))

               % Display time in title
         
               set(titl,'string',sprintf('t = %8.1f',t));
               gif_frame
               % writeVideo(VidObj,getframe(gcf))
               pause(0.025)
         end
            
         status = false;
      end
   end

   % ------------------------
   
   % Initial position
   
   x0 = 0.1;
   y0 = pi/2;
   theta0 = swinginit(x0,y0);

   % Initialize figure and buttons.
   
   clf reset
   shg
   set(gcf,'name','swinging sticks','menu','none','numbertitle','off') 

   % Start the ode solver
  
   % VidObj = VideoWriter('sticks.mp4','mpeg-4');
   % open(VidObj)

   delt = 0.1;
   u0 = [theta0; 0; 0];
   tspan = [0 10];
   opts = odeset('mass',@swingmass,'outputfcn',@swingplot,'maxstep',delt);
   ode23(@swingrhs,tspan,u0,opts);

   gif_frame('wrap')
   % close(VidObj)
end