function animate_rotation(angle,axis)
% Function to animate rotation through angle (radians) about axis (unit
% vector)

time = 20;
axis = axis/sqrt(dot(axis,axis));
omega = axis*(angle/time);

y0 = reshape(eye(3),[9,1]);
options = odeset('RelTol',0.0001);
sol = ode45(@(t,w) odefunc(t,w,omega),[0,time],y0,options);
rotation_movie(sol,[0,time]);

end

function  dwdt = odefunc(t,w,omega)
   R= transpose(reshape(w,[3,3])); %Current rotation matrix
   spin = [0,-omega(3),omega(2);omega(3),0,-omega(1);-omega(2),omega(1),0];
   Rdot = spin*R; % Rate of change of rotation matrix
   dwdt = reshape(transpose(Rdot),[9,1]); 
end

function rotation_movie(sol,time_interval)
%   Function to animate a rotation history calculated by ode45
%   sol:  ode solution variable object from ode45
%   time_interval[tstart,tstop]:  
%
    close all
    n_frames = 100; % # Frames in the movie
    a = [4,2,1];    % Dimensions [x,y,x] of the prism
%    
    initial_vertices = [0,0,0; a(1),0,0; a(1),a(2),0;...
                        0,a(2),0; 0,0,a(3); a(1),0,a(3);... 
                        a(1),a(2),a(3);0,a(2),a(3)];
%  Shift origin to center of prism
    for j=1:3
      initial_vertices(:,j) = initial_vertices(:,j)-a(j)/2;
    end
    face_matrix = [1,2,6,5;2,3,7,6;3,4,8,7;4,1,5,8;1,2,3,4;5,6,7,8];

    ax_length = 1.25*max(a)/2;
    figure1=figure;
    set(figure1,'Color','white');

    for i=1:n_frames
       t = time_interval(1) + (time_interval(2)-time_interval(1))*(i-1)/(n_frames-1);
       if (t>time_interval(2)) t = time_interval(2); end
       y = deval(sol,t);
       R= transpose(reshape(y,[3,3]));
       new_vertices = transpose(R*transpose(initial_vertices));
       clf;
       axes1 = axes('Parent',figure1,'Position',[0.3804 0.3804 0.44 0.46]);
       patch('Vertices',new_vertices,'Faces',face_matrix,'FaceVertexCdata',hsv(6),'FaceColor','flat');
       view(axes1,[142.1 26]);
       axis([-ax_length,ax_length,-ax_length,ax_length,-ax_length,ax_length])
%       axis([0,2*ax_length,0,2*ax_length,-1.5*ax_length,ax_length/2])
       grid(axes1,'on');
       set(axes1,'XAxisLocation','origin','YAxisLocation','origin');
       pause(0.05)
    end
end
