% -------------------------------------------------------------------------
%
% MATLAB for DUMMIES 13-14
% Exercice 2
%
% Solution detaillee
%  Vincent Legat, Leopold Cambier, Victor Colognesi, Corentin Damman,
%  Francois Heremans, Benoit Legat, Gauthier Limpens, Lucas Nyssens,
%  Ignace Ransquin, Gaetan Ransonnet, Alexandre Sanchez Falcon,
%  Adrien Scheuer, Harold Taeter, Joachim Van Verdeghem, Paul-Emile Bernard
%
% -------------------------------------------------------------------------
%
function soluce2()

%
% -1- Test propose dans l'enonce !
%     Attention, le test n'etait qu'un cas particulier des specifications
%     de la fonction demandee :-)
%
%     En particulier, rien dans l'enonce prevoyait que T = 1:n+1 !
%     Il fallait donc que votre fonction soit capable de traiter un
%     vecteur T tout-a-fait general : c'est loin d'etre toujours le cas.
%
%     La fonction de test ne permettait pas non plus de verifier que la
%     longueur du vecteur de sortie etait egale au vecteur d'entree "t"
%     Pas mal de programmes echouent sur ce test aussi.
%
%

close all; clear all; figure;
axis([0 10 0 10]); axis off; hold on;
X = []; Y = [];  n = 0;
disp('Left mouse button picks points.')
disp('Right mouse button picks last point.')
button = 1;
while button == 1
    n = n + 1;
    [X(n),Y(n),button] = ginput(1);
    plot(X(n),Y(n),'r.','Markersize',25);
end

X = [X X(1)];
Y = [Y Y(1)];
dX = [X(2)-X(n) X(3:n+1)-X(1:n-1) X(2)-X(n)]/2;
dY = [Y(2)-Y(n) Y(3:n+1)-Y(1:n-1) Y(2)-Y(n)]/2;


T = 1:n+1;
t = 1:0.001:n+1;
plot(hermite(t,T,X,dX),hermite(t,T,Y,dY),'b-','LineWidth',2);

%
% -2- Le test qu'on a realise pour corriger votre programme
%     Eh oui, souvent cela ne fonctionne pas pour pas mal d'etudiants.
%     Nous, on pensait naivement que le devoir etait facile : beh non :-)
%

n = 10; m = 1000;
T = [0:3*pi/(2*(n+1)):3*pi/2 2*pi];
t = linspace(T(1),T(end), m);
X  =  cos(T); Y  =  sin(T);
dX = -sin(T); dY =  cos(T);
figure;
plot(X,Y,'r.','Markersize',25); hold on;
plot(hermite(t,T,X,dX),hermite(t,T,Y,dY),'b-','LineWidth',2);
axis([-1.5 1.5 -1.5 1.5 ]); axis equal;


end

%
%  Il est utile d'observer qu'il existe une difference de rapidite
%  non negligeable entre "s.*(s+1)" et "s.^2+s" !
%  L'implementation efficace evite donc de faire appel a des puissances.
%  Voila, voila : pourquoi il faut factoriser l'expression :-)
%
%  Pour etre plus ou moins rapide, il faut essayer d'eviter des
%  boucles sur la longueur du vecteur "x"...
%  Eviter toute boucle n'est pas possible a mon avis en raison de la
%  detection de l'intervalle adequat :-)
%

function [uh] = hermite(x,X,U,dU)
%
% -1- Calcul des coefficients de m-1 polynomes sur chaque intervalle
%
m = length(X);
h = diff(X);
DU = diff(U)./h;
A  = (  3*DU - 2*dU(1:m-1) - dU(2:m))./h;
B  = (- 2*DU + dU(1:m-1)   + dU(2:m))./(h.*h);
%
% -2- Detection du polynome a choisir pour chaque
%     valeur x : i tel que  X(i) <= x < X(i+1)
%
i = ones(size(x));
for j = 2:m-1
   i(X(j) <= x) = j;
end
%
%  -3- Calcul de l'interpolant
%
s = x - X(i);
uh = U(i) + s.*(dU(i) + s.*(A(i) + s.*B(i)));
end