%
% Written by:
% -- 
% John L. Weatherwax                2005-08-14
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

clear functions; 
clear; 
close all; 
clc; 

randn('seed',0); rand('seed',0);

t_final = 1.0;


% Part (a) using the square root algorithm:
%
[t_a4,y_a4] = ode45( @sect_5_4_ode_fn_part_a, [0, t_final], [eps;eps;eps] );
P_a = [ y_a4(:,1).^2 + y_a4(:,2).^2, y_a4(:,2) .* y_a4(:,3), y_a4(:,3).^2 ]; % convert this output into components of the covariance matrix P(t) 

% Part (a) using the method from Problem 5.1:
%
[t_a1,y_a1] = ode45( @sect_5_1_ode_fn_part_a, [0, t_final], [0;0;0] );

% Compare the components of P(t) computed two different ways:
figure; subplot( 1, 3, 1 ); 
plot( t_a4, P_a(:,1), '-k' ); hold on;
plot( t_a1, y_a1(:,1), '-r' ); title('P(1,1)');

subplot( 1, 3, 2 ); 
plot( t_a4, P_a(:,2), '-k' ); hold on;
plot( t_a1, y_a1(:,2), '-r' ); title('P(1,2)');

subplot( 1, 3, 3 ); 
plot( t_a4, P_a(:,3), '-k' ); hold on;
plot( t_a1, y_a1(:,3), '-r' ); title('P(2,2)');

pause;
close all; 


% Part (b) using the Chandrasekhar type algorithm: 
%
[t_b4,y_b4] = ode45( @sect_5_4_ode_fn_part_b, [0, t_final], [eps;eps;eps] );
P_b = [ y_b4(:,1).^2 + y_b4(:,2).^2, y_b4(:,2) .* y_b4(:,3), y_b4(:,3).^2 ]; % convert this output into components of the covariance matrix P(t) 

% Part (b) using the method from Problem 5.1:
%
[t_b1,y_b1] = ode45( @sect_5_1_ode_fn_part_b, [0, t_final], [0;0;0] );

% Compare the components of P(t) computed two different ways:
figure; subplot( 1, 3, 1 ); 
plot( t_b4, P_b(:,1), '-k' ); hold on;
plot( t_b1, y_b1(:,1), '-r' ); title('P(1,1)');

subplot( 1, 3, 2 ); 
plot( t_b4, P_b(:,2), '-k' ); hold on;
plot( t_b1, y_b1(:,2), '-r' ); title('P(1,2)');

subplot( 1, 3, 3 ); 
plot( t_b4, P_b(:,3), '-k' ); hold on;
plot( t_b1, y_b1(:,3), '-r' ); title('P(2,2)');

pause;
close all; 


% Part (c) using the Chandrasekhar type algorithm: 
%
[t_c4,y_c4] = ode45( @sect_5_4_ode_fn_part_c, [0, t_final], [eps;eps;eps] );
P_c = [ y_c4(:,1).^2 + y_c4(:,2).^2, y_c4(:,2) .* y_c4(:,3), y_c4(:,3).^2 ]; % convert this output into components of the covariance matrix P(t) 

% Part (c) using the method from Problem 5.1:
%
[t_c1,y_c1] = ode45( @sect_5_1_ode_fn_part_c, [0, t_final], [0;0;0] );

% Compare the components of P(t) computed two different ways:
figure; subplot( 1, 3, 1 ); 
plot( t_c4, P_c(:,1), '-k' ); hold on;
plot( t_c1, y_c1(:,1), '-r' ); title('P(1,1)');

subplot( 1, 3, 2 ); 
plot( t_c4, P_c(:,2), '-k' ); hold on;
plot( t_c1, y_c1(:,2), '-r' ); title('P(1,2)');

subplot( 1, 3, 3 ); 
plot( t_c4, P_c(:,3), '-k' ); hold on;
plot( t_c1, y_c1(:,3), '-r' ); title('P(2,2)');

fn = [ '../../WriteUp/Graphics/Chapter4/sect_5_prob_4_part_c_plots.eps' ];
%saveas( gcf, fn, 'epsc' );