inertiaJ = diag([10.00; 10.00; 2.00]);
inertiaJinv = inv(inertiaJ);
alpha = 1000; % positive large constant for CSM
 
rotation_quaternion_ODE = @(t,q) rotation_quaternion_ODE_params (t,q, inertiaJ, inertiaJinv, alpha);
%interval = 0:0.0001:20;
interval = [0,20];
sinit = [1;0;0;0; 0;0;0;0];
[time, s] = ode45(rotation_quaternion_ODE, interval, sinit);

q0 = s(:,1); q1 = s(:,2); q2 = s(:,3); q3 = s(:,4);
coshalf = q0; sinhalf = sqrt(q1.^2+q2.^2+q3.^2);
angle = atan2(2*sinhalf.*coshalf, coshalf.^2-sinhalf.^2);
ux = q1./sinhalf; uy = q2./sinhalf; uz = q3./sinhalf;

zeta = zeta_axis(q0, q1, q2, q3);
plot3(zeta(:,1), zeta(:,2), zeta(:,3));
axis([-1,1, -1,1, -1,1]);
pbaspect([1 1 1]);
title('\zeta-axis');
pause;

plot(time,zeta(:,2),'-', time,zeta(:,3),'--');
axis([0,5, -1,1]);
title('\zeta_{y} and \zeta_{z} (time 0s - 5s)');
pause;
axis([5,10, -1,1]);
title('\zeta_{y} and \zeta_{z} (time 5s - 10s)');
pause;
