%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sample Matlab program for Homework 2.
%
% A Matlab script file that uses the function crossmat.m
% and plotBox.m to produce plots of a spinning frame
%
% Here's how to get started: 
%        1) Copy myprogram2A.m (this file), crossmat.m 
%           and plotBox.m into your working directory.
%        2) Invoke Matlab by typing "matlab" at
%           the unix prompt. 
%        3) Type "myprogram2A" in the Matlab window
%           to run this file.
% Note: the program may take a few seconds to run.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

e1 = [1 0 0]' ;    %define e1-coord vector (column) 
e2 = [0 1 0]' ;    %define e2-coord vector (column) 
e3 = [0 0 1]' ;    %define e3-coord vector (column) 

Id = eye(3,3) ;    %initialize 3x3 identity matrix

figure(1); clf ;   %create figure window #1 and clear its contents

a=1 ; b=1 ; c=1 ;        %set dimensions for plotBox function
origin=[0 0 0]' ;        %set origin for plotBox function
lim = 1.5*max([a,b,c]) ; %set limits for x,y,z axes

plotlimits = [-lim lim -lim lim -lim lim];

plot3(0,0,0) ;                    %initialize 3D plot
plotBox(origin,e1,e2,e3,a,b,c) ;  %plot box defined by e-vecs
axis(plotlimits) ;                %set limits for plot axes
axis square  ;                    %set scales for plot axes 
axis manual  ;                    %freeze scales in plot
pause(1.0)   ;                    %pause to draw and see plot


omega = (1/sqrt(3))*[1 1 1]' ;  %set angular velocity 

T = 2*pi ; N = 50 ; dt = T/N ; 

time  = zeros(N+1,1) ; %define array for time values
ddata = zeros(N+1,1) ; %define array for determinant data

t = 0;    %initial time
Qn = Id ; %initial condition for Q

time(1) = t ;           %store initial t value
ddata(1) = det(Qn) ;  %store initial determinant

for i=1:N
   t = t + dt;                        %increment time
   Qn = Qn*(Id + crossmat(dt*omega)); %compute next matrix w/Euler method

   d1=Qn*e1 ; d2=Qn*e2 ; d3=Qn*e3 ;   %compute d-vecs 
   clf;                               %clear figure 
   plotBox(origin,d1,d2,d3,a,b,c) ;   %plot box defined by d-vecs
   axis(plotlimits) ;                 %set limits for plot axes
   axis square ;                      %set scales for plot axes 
   axis manual ;                      %freeze scales in plot
   pause(0.1) ;                       %briefly pause to draw and see plot

   time(i+1) = t ;                    %store t value
   ddata(i+1)= det(Qn) ;              %store det value
end

figure(2) ;clf; plot(time,ddata) ;    %plot det(Qn) vs t 
xlabel('time: t') ;                   %set label for x axis
ylabel('determinant: det(Qn)') ;      %set label for y axis

Qexact = Id*expm(crossmat(T*omega)) ; %compute exact solution at time T
Qerror = norm(Qn-Qexact, 2)   ;       %compute matrix 2-norm of error

fprintf('When N=%g we get Error=%g \n', N, Qerror) ;  %print N, Qerror to screen

