function [Lappr,Lexa]=longueur_approchee_courbe(N,a,b,x,y,xp,yp,gra,chometh)

% Calcul de la longueur exacte et approchée d'une courbe (par ligne polygonale)
%
% *********************************************************
% % [Lappr,Lexa]=longueur_approchee_courbe(N,a,b,x,y,xp,yp,gra,chometh) :  
% * Variables d'entrées :
%   *N : nombres d'intervalles de discrétisation
%   *a,b : bornes de l'intervalle où vit t
%   *x,y : fonction x(t) et y(t)
% * Variables d'entrées optionnelles :
%   *x',y' : fonction x'(t) et y'(t)
%   * gra : entier dans{0,1}, égal à 1 par défaut pour faire le graphe
%   * chometh : entier dans{0,1}, égal à 1 pour la méthode vectorielle plus
%   rapide
% * Variables de sorties optionnelles :
%   * Lappr : longueur approchée
%   * Lexa : longueur exacte
%
% ************ Fonctions auxiliaires utilisées ************
%
%  Aucune
%
% *********************************************************
%
% 2025 by          Jérôme BASTIEN
%                  Université Claude Bernard Lyon I, Polytech, Laboratoire LIBM, Villeurbanne
%                  E-Mail : jerome.bastien@univ-lyon1.fr
%
%
% exemple
%   [Lappr,Lexa]=longueur_approchee_courbe(100,0,1,@(t) t,@(t) t.^2/2,[],[],1)

chs=0;
if nargin<=5||isempty(xp)
    syms t;
    chs=1;
    xp=inline(vectorize((diff(x(t)))));
    yp=inline(vectorize((diff(y(t)))));
end
if nargin<=7||isempty(gra)
    gra=1;
end
if nargin<=8||isempty(chometh)
    chometh=1;
end
if nargout>=2
    if ~chs
        syms t;
    end
    Lexa=int(sqrt((xp(t))^2+(yp(t))^2),t,0,1);
end
h=(b-a)/N;
Nmax=100;
if gra
    if chometh==0
        xg=zeros(1,N+1);
        yg=xg;
        if N<=Nmax
            xgfin=zeros(1,Nmax+1);
            ygfin=xg;
        end
    end
end
if gra
    if chometh==0
        hfin=(b-a)/Nmax;
        for i=0:Nmax
            xgfin(i+1)=x(a+i*hfin);
            ygfin(i+1)=y(a+i*hfin);
        end
    else
        tfin=linspace(a,b,Nmax+1);
        xgfin=x(tfin);
        ygfin=y(tfin);
    end
end
if chometh==0
    s=0;
    auxix=x(a);
    auxiy=y(a);
    for i=0:N-1
        if gra
            xg(i+1)=auxix;
            yg(i+1)=auxiy;
        end
        auxixn=x(a+(i+1)*h);
        auxiyn=y(a+(i+1)*h);
        if nargout>=1
            s=s+hypot(auxixn-auxix,auxiyn-auxiy);
        end
        auxix=auxixn;
        auxiy=auxiyn;
    end
    if gra
        xg(N+1)=auxix;
        yg(N+1)=auxiy;
    end
else
    if gra||nargout>=1
        tt=linspace(a,b,N+1);
        xs=x(tt);
        ys=y(tt);
    end
    if gra
        if N<=Nmax
            tfin=linspace(a,b,Nmax+1);
            xgfin=x(tfin);
            ygfin=y(tfin);
        end
    end
    if nargout>=1
        s=sum(hypot(xs(1:end-1)-xs(2:end),ys(1:end-1)-ys(2:end)));
    end
end
if nargout>=1
    Lappr=s;
end
if gra
    hold on
    if chometh==0
        if N<=Nmax
            for i=1:N
                plot([xg(i),xg(i+1)],[yg(i),yg(i+1)],xg(i),yg(i),'.');
            end
            for i=1:Nmax
                plot([xgfin(i),xgfin(i+1)],[ygfin(i),ygfin(i+1)],'r');
            end
        else
            for i=1:N
                plot([xg(i),xg(i+1)],[yg(i),yg(i+1)]);
            end
        end
    else
        if N<=Nmax
            plot(xgfin,ygfin,xs,ys);
            plot(xs,ys,'.');
            legend('courbe','ligne brisée','location','best');
        else
            plot(xs,ys);
        end
    end
    axis equal;
    hold off;
end