// S.Ayrinhac (octobre 2014)
// détermine l'aire d'un polygone quelconque
// exemple ci-dessous : 
// theta=linspace(0,2*%pi,4);x=cos(theta);y=sin(theta);N=10;A=Func_aire_poly(x,y,N);

function A=Func_aire_poly(x,y,N)

//N=5;
epsilon=1e-5;
lengthx=length(x);

deltax=abs(min(x)-max(x))*1.1;
deltay=abs(min(y)-max(y))*1.1;
//Nx=floor(deltax/0.1); 
//Ny=floor(deltay/0.1); 

//x00=linspace(min(x),max(x),20);
x0=linspace(min(x),max(x),N);
stepx=abs(x0(2)-x0(1));
x0=x0(2:$-1);
y0=linspace(min(y),max(y),N);
stepy=abs(y0(2)-y0(1));
y0=y0(2:$-1);

// u : aire élémentaire
u=stepx*stepy;
[X,Y]=meshgrid(x0,y0);
compteur=1;
//plot(X,Y,'.r')

//// ligne par ligne
for j=1:size(Y,1)
    hauteurligne=Y(j,1);
    //// on cherche le(s) segment(s) concerné(s)
    m=1;
    for k=1:lengthx-1
        minx=min(x(k),x(k+1)+epsilon); //// !!! attention si égaux
        maxx=max(x(k),x(k+1)+epsilon);
        miny=min(y(k),y(k+1));
        maxy=max(y(k),y(k+1));
        //plot([min(x) max(x)],[hauteurligne hauteurligne],'-c')
        if hauteurligne>=miny
        if hauteurligne<=maxy
           //plot([x(k) x(k+1)],[y(k) y(k+1)],'-g')
           
           pente=(y(k+1)-y(k))/(epsilon+x(k+1)-x(k));    //// pente
           ordon=y(k)-pente*x(k);        //// ordonnée à l'origine
           coeffa(m)=pente;
           coeffb(m)=ordon;
           //plot(x00,pente.*x00+ordon,'m')
           //pause
           m=m+1;
        end // if hauteurligne>=miny
        end // if hauteurligne<=maxy
        
    end
    
    //// on calcule les X des intersections concernées
    n=1;
    xlimite=1;
    for n=1:m-1
        xlimite(n)=(hauteurligne-coeffb(n))/coeffa(n);
        //plot(xlimite(n),hauteurligne,'+m')
    end
    //// on sélectionne les X qui tombent dans les limites
    Xligne=X(j,:);
    Yligne=Y(j,:);
    xlimite=floor(1e7*xlimite)/1e7; //// tronque la précision
    xlimite=unique(xlimite);  //// trie et supprime les doublons
    t=1;
    I=1;
    //// vérification que t est pair //////////////////////////
    if Func_even(t)==0
        t=t-1;
        if t==0
          t=1;
        end
    end
    //////////////////////////////////////////////////////////////////////////////
    m=min(m,length(xlimite));
    xlimite;
    for t=1:2:m-1
       I=find((Xligne>=xlimite(t))&(Xligne<=xlimite(t+1)));
       //plot(Xligne(I),Yligne(I),'+r')
       if length(I)>0
       for w=1:length(I)
           zx(compteur)=Xligne(I(w));
           zy(compteur)=Yligne(I(w));
           //plot(Xligne(I(w)),Yligne(w),'+r')
           compteur=compteur+1;
       end
       end
    end
    //pause
end

// quand la surface ne continet aucun point !
if compteur<=1
  zx=mean(x);
  zy=mean(y);
end

// Aire : nombre de points dedans fois surface élémentaire
A=length(zx)*u;


endfunction


// Fonction qui renvoie 1 si x est pair
function b=Func_even(x)

  b=zeros(1,length(x));
  I=find(modulo(x,2)==0);
  b(I)=1;

endfunction