function[x,y,z] = ellipsoide0(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] = ellipsoide(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])
  theta = linspace(0,2.*%pi,n+1)
  phi = linspace(-%pi./2,%pi./2,m+1)
  
//  Evalution des facettes et tracé de l'ellipse :
  [X,Y,Z] = eval3dp(ellipsoide0,theta,phi)
  clf()
  plot3d1(X,Y,Z)

// Matrice des aretes A :  
  A = zeros(4.*n.*m,6)
  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)]
    end
    A(4.*(k-1)+4,:) = [M(4,k),M(8,k),M(12,k),M(1,k),M(5,k),M(9,k)]
  end
// Suppression des doublons dans A :
  for k=[1:2.*n.*m]
    i = vectorfind(A,[A(k,4),A(k,5),A(k,6),A(k,1),A(k,2),A(k,3)],'r')
    A(i,:) = []
//    B = A, B(k,:) = []
//    i = vectorfind(B,[A(k,1),A(k,2),A(k,3),A(k,4),A(k,5),A(k,6)],'r')
//    A(i,:) = []
  end
endfunction


