%
% MATLAB for DUMMIES 11-12
% Exercice 6
%
% Solution detaillee
%  Vincent Legat, Damien Bouvy, Denis-Gabriel Caprace, Gilles De Neyer, 
%  Kevin Degraux, Mohammad Iranmanesh, Christophe Lorant, Remy Moulin, 
%  Matthew Philippe, Karim Slaoui, Olivier Thiry, Laurent Van Eesbeeck,
%  Thomas Vanderstraeten, Antoine Weynants
%
% -------------------------------------------------------------------------
%
function [tf,v0,n] = shoot(alpha,M)

%
% -1- Verification des donnees
%     C'est un peu bete, mais toujours utile
%

    alpha = abs(alpha);
    if (M <= 0) error('A negative mass and what else, stupid boy...'); end
    tf = 0;

%
% -2- Detection d'un intervalle qui encadre la solution
%     On utilise la solution sans frottement comme borne inferieure
%     On multiplie cette valeur par MAGIC pour obtenir la valeur superieure
%     On repete cette multiplication aussi longtemps que necessaire.
%     Normalement c'est robuste :-)   Et c'est SIMPLE !
%
%     Choisir MAGIC = 2   permet d'obtenir vite un intervalle de depart
%     Choisir MAGIC = 1.2 necessite plus de temps, mais l'intervalle est
%          plus petit et on ira plus vite apres.
%     Choisir MAGIC = 1.1 prendra beaucoup de temps mais l'intervalle sera
%          tres petit
%
%     Heuristiquement, choisir MAGIC = 1.2 n'est pas une mauvaise idee
%
    magic = 1.2;
    a = sqrt(9.81*alpha);  fa = f(a);
    b = magic*a;           fb = f(b);
    n = 2;
    fprintf('\n [a,b] = [%21.14e %21.14e] m/s at %2i iteration  ',a,b,n);

    nmax = 50;
    while (fa*fb > 0) & n <= nmax
        a = b;        fa = fb;
        b = magic*a;  fb = f(b);
        n = n+1;
        fprintf('\n [a,b] = [%21.14e %21.14e] m/s at %2i iteration  ',a,b,n);
    end

%
% -3- Implementation de la methode de bissection 
%
%     Le test d'erreur doit s'effectuer sur l'ordonnee 
%     et non sur l'abscisse de la fonction :-)
%

    tol = 0.01;
    value = tol + 1;

    while abs(value) >= tol & n <= nmax
        delta = (b-a)/2;
        x = a + delta;
        value  = f(x);
        n = n + 1;
        if (value*f(a) > 0)
            a = x;
        else
            b = x;
        end        
        fprintf('\n x = %21.14e m/s (Estimated error %21.14e at %2i iteration ) ',x,abs(value),n);
    end
    v0 = x;

%
% -4- Description du probleme du tir 
%
% -4.1- Fonction objective a annuler
%       f(x) = différence entre le point d'impact et la cible
%       x = norme initiale de la vitesse (hausse de 45 degrés)
%

    function y  = f(x)
        theta = pi/4;
        Xspan  = [0 100];
        u = x* cos(theta);
        v = x* sin(theta);
        options = odeset('RelTol',1e-3,'AbsTol',1e-3,'Event',@events);
        [X U] = ode45(@g,Xspan,[0.0 u 0.0 v],options);
        plot(U(:,1),U(:,3),'-r'); hold on;
        y = alpha - U(end,1);
        tf = X(end);
    end

%
% -4.2- Integration de la trajectoire de l'obus
%       u = (x, dxdt, y, dydt)
%
%

    function dudx = g(x,u)
        dudx = u;
        frott =  0.1 * sqrt(u(2)*u(2) + u(4)*u(4)) / M;
        dudx(1) = u(2);
        dudx(2) = -frott * u(2);
        dudx(3) = u(4);
        dudx(4) = -9.81 - frott * u(4);
    end

%
% -4.3- Critere d'arret de l'integration par la gestion d'evenement
%

    function [value,isterminal,direction] = events(t,u)
        value = u(3);     % Detecte la hauteur  y = 0
        isterminal = 1;   % Arrete l'integration
        direction = -1;   % On considère uniquement le cas 
                          % ou l'obus tombe du ciel vers la terre
    end

end