function Euler_analyt(w_i,I,d_time,t_f,cmp)
%
% Analytical solution to the Euler's equatons of motion
% Using the equations in Wittenburg (p.48:50)
%
% Input:
% w_i     - 3x1 vector with the initial angular velocities
% I       - 3x1 vector with the body principal moments of inertia (I(1) > I(2) > I(3))
% d_time  - step for evaluating the analytical solution
% t_f     - simulation time
% cmp     - if 'Y' the analytical solution is comapared with a numerical one.
%           in this case d_time is used as integration step
%
% Example:
%
% Euler_analyt([0.3;0.4;-0.2],[6;4;2],0.01,5,'Y')
%
% For information see the report:
% D. Dimitrov, "Derivation of the analytical solution for torque free motion of a rigid body".
%
% 2005/08/05
% Dimitar Dimitrov
%

% Check weather  I(1) > I(2) > I(3)
if I(1) <= I(2) || I(2) <= I(3) || I(1) <= I(3)
    disp('You should input I(1) > I(2) > I(3)')
    return
end

I1 = I(1);
I2 = I(2);
I3 = I(3);

% Calculate T and h
T = 1/2*I1*w_i(1)^2 + 1/2*I2*w_i(2)^2 + 1/2*I3*w_i(3)^2;
h = sqrt(I1^2*w_i(1)^2 + I2^2*w_i(2)^2 + I3^2*w_i(3)^2);

% ======================================================
%% Constants definition
% ======================================================
D = h^2/(2*T);

a = sqrt( 2*T*(D - I3)/( I2*(I2-I3) ) );
b = sqrt( 2*T*(I1 - D)/( I2*(I1-I2) ) );

% determining the signs of w_i(1) and w_i(3)
s_1 = sign(w_i(1));
s_3 = sign(w_i(3));

% If the value of w_i(1) or w_i(3) is zero its sign is set to "plus"
if s_1 == 0
    s_1 = 1;
end
if s_3 == 0
    s_3 = 1;
end

% ======================================================
%% Solution of Euler's equations
% ======================================================
% counter
i = 1;

if D > I2
    disp('Case1 (D > I2)')

    k = b/a;
    m = k^2; % the function "ellipj" uses "m" instead of "k" (m = k^2)

    w_1m = sqrt( ( 2*T*(D-I3) )/( I1*(I1-I3) ) );
    w_2m = b;
    w_3m = sqrt( ( 2*T*(I1-D) )/( I3*(I1-I3) ) );

    P = sqrt( ( 2*T*(I1-I2)*(D-I3) ) / ( I1*I2*I3 )  );

    % Upper limit (of integration)
    ul = -s_1*s_3*w_i(2)/(b);

    % Computes Elliptic integrals of the first kind
    t_0_tmp = mfun('EllipticF',ul,k);
    t_0 = -t_0_tmp/P;

    % It can be computed numerically (of course)
    % t_0_tmp = quad(@(x) 1./(sqrt(1-x.^2).*sqrt(1-k^2*x.^2)), 0, ul)

    for time = 0:d_time:t_f

        [SN,CN,DN]  =  ellipj(P*(time - t_0),m);

        w1(i) = s_1*w_1m*DN;
        w2(i) = -s_1*s_3*w_2m*SN;
        w3(i) = s_3*w_3m*CN;

        i = i+1;

    end

end

% Case 3 is included in Case 2 (by including D = I2)
if D <= I2
    if D == I2
        disp('Case3 (D = I2)')
    else
        disp('Case2 (D < I2)')
    end

    k = a/b;
    m = k^2;

    w_1m = sqrt( ( 2*T*(D-I3) )/( I1*(I1-I3) ) );
    w_2m = a;
    w_3m = sqrt( ( 2*T*(I1-D) )/( I3*(I1-I3) ) );

    P = sqrt( ( 2*T*(I1-D)*(I2-I3) ) / ( I1*I2*I3 )  );

    ul = -s_1*s_3*w_i(2)/(a);

    t_0_tmp = mfun('EllipticF',ul,k);
    t_0 = -t_0_tmp/P;

    for time = 0:d_time:t_f

        [SN,CN,DN]  =  ellipj(P*(time - t_0),m);

        w1(i) = s_1*w_1m*CN;
        w2(i) = -s_1*s_3*w_2m*SN;
        w3(i) = s_3*w_3m*DN;

        i = i+1;
    end

end

if strcmp(cmp,'Y')

    % Computes the profile of "w" numerically

    [T,Y] = ode45(@rot_motion_eq,[0:d_time:t_f],w_i);
    w = Y';

    hold on;
    % Plots the numerical solution
    figure(1);plot(w(1,:),'k')
    figure(1);plot(w(2,:),'k')
    figure(1);plot(w(3,:),'k')

end

hold on;
% Plots the analytical solution
figure(1);plot(w1,'r')
figure(1);plot(w2,'b')
figure(1);plot(w3,'g')

% Nested function implementing the Euler equations of motion
    function wd = rot_motion_eq(t,w)
        wd = zeros(3,1);
        wd(1) = ((I2 - I3)/I1)*w(2)*w(3);
        wd(2) = ((I3 - I1)/I2)*w(1)*w(3);
        wd(3) = ((I1 - I2)/I3)*w(1)*w(2);
    end

end