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[d] = Distance(u,v)
  // u et v vecteurs lignes de R3
  d = sqrt((u(1)-v(1))^2+(u(2)-v(2))^2+(u(3)-v(3))^2)
endfunction

function[A] = MaillageEllipsoide(a,b,c,n,m)
  // Evaluation de l'ellipsoide en un nombre fini de points
  // n : discrétisation de theta ([-pi,pi])
  // m : discrétisation de phi ([-pi/2,pi/2])
  // 1)  Evalution des facettes et tracé de l'ellipse :
  theta = linspace(-%pi,%pi,n+1)
  phi = linspace(-%pi./2,%pi./2,m+1)
  [X,Y,Z] = eval3dp(Ellipsoide,theta,phi)
  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) = Distance([A(k,1),A(k,2),A(k,3)],[A(k,4),A(k,5),A(k,6)])
  end
endfunction

function[S] = NumerotationSommets(A,n,m)
  S = [A(:,1),A(:,2),A(:,3);A(:,4),A(:,5),A(:,6)]
  for k=1:(n.*(m-1)+2)
    i = vectorfind(S,[S(k,1),S(k,2),S(k,3)],'r') 
    p = size(i,2)
    for j=p:-1:2
      S(i(j),:) = []
    end
  end
endfunction

function[D] = Aretes2MatriceDistance(A,S,n,m)
D = %inf*ones(n.*(m-1)+2,n.*(m-1)+2)
for i=1:(n.*(m-1)+2), D(i,i) = 0, end
  for k=1:n.*(2.*m-1)
    i = vectorfind(S,[A(k,1),A(k,2),A(k,3)],'r')
    j = vectorfind(S,[A(k,4),A(k,5),A(k,6)],'r')
    D(i(1),j(1)) = A(k,7)
    D(j(1),i(1)) = A(k,7)
  end
endfunction

function[G] = MatriceDistance2Graphe(D,nom)
  n = size(D,1)
  H = zeros(1,n^2)
  T = zeros(1,n^2)
  L = zeros(1,n^2)
  k = 0
  for i=2:n
    for j=1:(i-1)
      if D(i,j)<>%inf
        k = k+1
        H(k) = i
        T(k) = j
        L(k) = D(i,j)
      end
    end
  end
  G = make_graph(nom,0,n,T(1:k),H(1:k))
  G.edge_length = L(1:k)
endfunction  

function[D,P] = AlgoFloyd(D)
  // Initialisation :
  n = size(D,1)
  P = zeros(n,n)
  for i=1:n
    for j=1:n
      P(i,j)=i
    end
  end
  for k=1:n
    for i=1:n
      for j=1:n
        if D(i,k)+D(k,j)<D(i,j)
              D(i,j) = D(i,k)+D(k,j)
              P(i,j) = P(k,j)
        end
      end
    end
  end
endfunction

function[D,P] = AlgoDijkstra(M,x)
  n = size(M,1)
  D = M(x,:)
  P = x*ones(1,n)
  L = list()
  for k=1:n, L($+1)=k, end
  L(x) = null()             
  while size(L)<>0
    d = %inf
    for k=1:size(L)
      if D(L(k))<d
        d = D(L(k))
        i = k
      end
    end
    if d==%inf, break, end
    for z=L
      if D(z)>d+M(L(i),z)
        D(z) = d+M(L(i),z)
        P(z) = L(i)
      end
    end
    L(i) = null() 
  end
endfunction

function[phi] = DegreSexa2Radians(theta)
  // theta = [degres, minutes, secondes, signe]
  phi = theta(4).*(%pi./180).*(theta(1)+theta(2)./60+theta(3)./3600)
endfunction

function[i] = SommetPlusProche(lat,lon,S)
  theta = DegreSexa2Radians(lon)
  phi = DegreSexa2Radians(lat)
  [x,y,z] = Ellipsoide(theta,phi)
  s = Distance(S(1,:),[x,y,z])
  i = 1
  for k=2:size(S,1)
    d = Distance(S(k,:),[x,y,z])
    if d<s
      i = k
      s = d
    end
  end
endfunction

function[C] = Geodesique(P,j)
// P : vecteur ligne
  C = list(j)
  while P(j)<>C(1)
    C(0) = P(j)
    j = P(j)
  end
endfunction

function[] = TraceArete(v,p,couleur,epaisseur)
  t=linspace(0,1,p)
  param3d1(v(1)+t.*(v(4)-v(1)),v(2)+t.*(v(5)-v(2)),list(v(3)+t.*(v(6)-v(3)),couleur.*ones(1,p)))
  e = gce()
  e.thickness = epaisseur
endfunction

function[] = TraceChemin(C,p,couleur,epaisseur)
  for i=1:(length(C)-1)
    TraceArete([S(C(i),:),S(C(i+1),:)],p,couleur,epaisseur)
  end
endfunction

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

function[] = TraceEllipsoide2(S,M,n,m,p)
  for i=2:n.*(m-1)+2
    for j=1:i-1
      if M(i,j)>0 
        TraceArete([S(i,:),S(j,:)],p,0,1) 
      end
    end
  end
endfunction
