function [Oout Rout Pout normout]=alphahull(X,r0)
d=size(X,2);
N=size(X,1);

%Delaunay Triang, Radius and Center
Tri=delaunayn(X);
[V C]=voronoin(X);
[R O]=rayonscir(Tri,X);
Nt=size(Tri,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Bout, Pout
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 1/Bout, ceux qui sont "loins
R00=r0;
ii=find(R>=R00);
Oout=O(ii,:);
Rout=R(ii);
Nout=size(Rout,1);

%% 2/ Pout
[K VV]=convhulln(X);
Pout=X(K(:,1),:); %pt de la face
normout=zeros(size(K,1),d);

for i=1:size(K,1)
    MatFace=zeros(d-1,d);
    for j=1:d-1
        MatFace(j,:)=X(K(i,j+1),:)-X(K(i,1),:);
    end
    MatFace;
    eta=null(MatFace);
    neta=norm(eta);
    MatFace*eta;
    XX=X-ones(N,1)*Pout(i,:);
    PV=sum(XX*eta);
    
    if PV>0
        normout(i,:)=-eta'/neta;
        etafin=-eta;
    else
        normout(i,:)=eta'/neta;
        etafin=eta;
    end
end


%% 3/ Bout ceux qui sont de rayon R0
Kt=sort(K,2);

Trit=sort(Tri,2);
%fait les faces et met les rayons
Nt=size(Trit,1);
Ft=[];
Isimp=[];
for j=1:d+1
    Tf=Trit;
    Tf(:,j)=[];
    Ft=[Ft;Tf];
    Isimp=[Isimp;[1:1:Nt]'];%fait les faces et numérote les simplexes.
end

[simp,IA]=unique(Ft,'rows','first');
[simp,IB]=unique(Ft,'rows','last');
[Rf Of]=rayonscir(simp,X);
%simp est la face, Isimp(IA) est le numéro d'un simplexe qui contiens la face Isimp(IB) est le numéro de l'autre simplexe qui contiens la face
angle=[0:2*pi/200:2*pi]';

for i=1:size(simp,1);
    %un simplexe avec un rayon < r0 et l'autre avec un rayon >r0
    R1=R(Isimp(IA(i)));
    R2=R(Isimp(IB(i)));
    if and(min(R1,R2)<r0,max(R1,R2)>=r0)
        if R1<R2
            etaout=O(Isimp(IB(i)),:)-Of(i,:);
            etaout=etaout/norm(etaout);
            newOout=Of(i,:)+sqrt(r0^2-Rf(i)^2)*etaout;
            Oout=[Oout;newOout];
            Rout=[Rout;r0];
        else
            etaout=O(Isimp(IA(i)),:)-Of(i,:);
            etaout=etaout/norm(etaout);
            newOout=Of(i,:)+sqrt(r0^2-Rf(i)^2)*etaout;
            Oout=[Oout;newOout];
            Rout=[Rout;r0];
        end
    end
    
    %un seul simplexe
    if IA(i)==IB(i)
        if R(Isimp(IA(i)))<=r0
            [simpl ii1 iiK]=intersect(simp(i,:),Kt,'rows');
            
            etaout=normout(iiK,:);
            newOout=Of(i,:)+sqrt(r0^2-Rf(i)^2)*etaout;
            Oout=[Oout;newOout];
            Rout=[Rout;r0];
        end
    end
end

