function [h]=ver3h(xc,yc,zc,r,ntet)
% [h]=VER3(XC,YC,ZC,R,NTET)
%
% La fonction représente un tube ou cylindre courbé le long
% de la courbe donnée par [Xc Yc Zc] de rayon R.
% La courbe doit au moins comporter 2 points. On tient compte
% de sa périodicité si la courbe est fermée.
% R peut être un scalaire ou un vecteur (de la taille de Xc).
% NTET est la résolution, le nombre de points pour représenter
% un cercle (par défaut, 12).
%
% Une représentation par SURF est effectuée s'il n'y a pas de
% paramètre en sortie. Sinon les coordonnées du ver sont retournées.

% Origine: sans nom sur le serveur de MathWorks.
% Adaptation française: J.-B. Roux, Collège de Saussure, Genève, 1996

% Validation des données
if nargin<4
	help ver3h
	error('Il manque des arguments');
end
if nargin<5
	ntet=12;	% valeur par défaut
end
if ntet<2
	error('Au moins deux points sont nécessaires')
end
n=length(xc);	% nombre de points de la courbe
if length(yc)~=n | length(zc)~=n
	error('Les dimensions en X, Y et Z doivent être identiques')
end

% Initialisation
rayon=r;
if length(r)<n	% vecteur de rayons égaux
	rayon=r*ones(size(xc));
end

tet=(0:2*pi/ntet:2*pi)';	% angles pour un cercle
ct=cos(tet);
st=sin(tet);

X=zeros(n,ntet+1);			% réservation de la place
Y=X;
Z=X;

% Construction des (n-1) cercles normaux à chaque segment
for i=1:n-1
	v=[xc(i+1)-xc(i) yc(i+1)-yc(i) zc(i+1)-zc(i)];
	V=norm(v);		% vecteur tangent à la courbe
	salpha=v(3)/V;	
	calpha=sqrt(v(1)*v(1)+v(2)*v(2))/V;
	sbeta=0;
	cbeta=1;
	if calpha~=0
		sbeta=v(2)/V/calpha;
		cbeta=v(1)/V/calpha;
	end
	x=ct*salpha*cbeta+st*sbeta;
	y=ct*salpha*sbeta-st*cbeta;
	z=-ct*calpha;
	X(i,:)=x'; 
	Y(i,:)=y'; 
	Z(i,:)=z';
end

% Déplacement et mise à l'échelle des cercles
X(n,:)=X(n-1,:)*rayon(n)+xc(n);		% n-ième cercle
Y(n,:)=Y(n-1,:)*rayon(n)+yc(n);
Z(n,:)=Z(n-1,:)*rayon(n)+zc(n);
X1=X(1,:); Y1=Y(1,:); Z1=Z(1,:);
for i=2:n-1
	% rotation du cercle pour faire coïncider avec le préc.
	dx=X1(1)-X(i,:);
	dy=Y1(1)-Y(i,:);
	dz=Z1(1)-Z(i,:);
	[a k]=min(dx.*dx+dy.*dy+dz.*dz);
	kk=[k:ntet 1:k];	% permutation des indices
	X2=X(i,kk); Y2=Y(i,kk); Z2=Z(i,kk);
	X(i,:)=(X1+X2)/2*rayon(i)+xc(i);
	Y(i,:)=(Y1+Y2)/2*rayon(i)+yc(i);
	Z(i,:)=(Z1+Z2)/2*rayon(i)+zc(i);
	X1=X2; Y1=Y2; Z1=Z2;
end
X(n,:)=X(n,kk);		% rotation du n-ième cercle
Y(n,:)=Y(n,kk);
Z(n,:)=Z(n,kk);
X(1,:)=X(1,:)*rayon(1)+xc(1);	% 1er cercle
Y(1,:)=Y(1,:)*rayon(1)+yc(1);
Z(1,:)=Z(1,:)*rayon(1)+zc(1);

% En cas de périodicité on ajuste le 1er et le dernier cercle
d=(xc(1)-xc(n))^2+(yc(1)-yc(n))^2+(zc(1)-zc(n))^2;
if d<1e-16
	X(n,:)=(X(1,:)+X(n,:))/2;
	Y(n,:)=(Y(1,:)+Y(n,:))/2;
	Z(n,:)=(Z(1,:)+Z(n,:))/2;
	X(1,:)=X(n,:);
	Y(1,:)=Y(n,:);
	Z(1,:)=Z(n,:);
end

% Représentation du ver
if nargout==0
	surfl(X,Y,Z)
else
	h=surfl(X,Y,Z);
end

