function[x,y,z] = ellipsoide(theta,phi)
// Paramétrisation cartésienne de l'ellipsoide avec :
// a : demi grand axe
// b : demi axe moyen
// c : demi petit axe
  x = a.*cos(theta).*cos(phi)
  y = b.*sin(theta).*cos(phi)
  z = c.*sin(phi)
endfunction

function[A] = MaillageEllipsoide(a,b,c,n,m)
// Evaluation de l'ellipsoide en un nombre fini de points
// n : discrétisation de theta ([0,2pi])
// m : discrétisation de phi ([-pi,pi])
//1)  Evalution des facettes et tracé de l'ellipse :
  theta = linspace(0,2.*%pi,n+1)
  phi = linspace(-%pi./2,%pi./2,m+1)
  [X,Y,Z] = eval3dp(ellipsoide,theta,phi)
  clf()
  plot3d1(X,Y,Z)

//2) Matrice des aretes A :  
  A = zeros(4.*n.*m,7)
  M = [X;Y;Z]
  for k=[1:n.*m], 
    for i=[1:3]
    A(4.*(k-1)+i,:) = [M(i,k),M(i+4,k),M(i+8,k),M(i+1,k),M(i+5,k),M(i+9,k),0]
    end
    A(4.*(k-1)+4,:) = [M(4,k),M(8,k),M(12,k),M(1,k),M(5,k),M(9,k),0]
  end
  A = clean(A) // Arrondit les coordonnées des points car calculs approchés  

//3) Suppression des doublons dans A et calcul des longueurs des aretes :
// Pour les facettes à 3 côtés aux pôles (aretes reliant un point à lui même)
  i = [vectorfind(A,[0,0,c,0,0,c,0],'r'),vectorfind(A,[0,0,-c,0,0,-c,0],'r')]
  for j=sort(i)
    A(j,:) = []
  end
// Une arete appartient à 2 facettes :
  for k=[1:n.*(2.*m-1)]
    i = vectorfind(A,[A(k,4),A(k,5),A(k,6),A(k,1),A(k,2),A(k,3),0],'r') 
    A(i(1),:) =  []
    A(k,7) = sqrt((A(k,4)-A(k,1))^2+(A(k,5)-A(k,2))^2+(A(k,6)-A(k,3))^2)
end
endfunction

function[] = TraceArete(v,n)
  t=linspace(0,1,n)
  param3d(v(1)+t.*(v(4)-v(1)),v(2)+t.*(v(5)-v(2)),v(3)+t.*(v(6)-v(3)))
endfunction

function[] = TraceEllipsoide(A,n)
  m = size(A)
  for k=1:m(1)
    TraceArete(A(k,:),n)
  end
endfunction

//function[] = MatriceVersGraphe(A)
