function fords(n,style,coul)
% FORDS(N,STYLE,COUL) reprŽsente les sph�res de Ford selon la sŽrie de Farey 
% d'ordre N. STYLE vaut TOUS pour tous les cercles, SANSUN pour ™ter les 1
% DEMI pour afficher la moitiŽ de la sŽrie. COUL est une matrice de couleur
% pour le remplissage des cercles.

% initialisation
switch nargin
    case 0
        help ford
    case 1
        style='tous';
        coul=prism(256);
    case 2
        coul=prism(256);
end

% dŽfinition des cercles
[x,y,z]=sphere;
h=plot3(0,0,0);
hold on
delete(h)

% tableau des fractions non rŽductibles dont le dŽnominateur ne dŽpasse pas N 
% arrangŽes dans l'ordre de grandeur entre 0 et 1
fs=[];
for i=1:n
    for j=i+1:n
        if gcd(i,j)==1  % test de rŽduction
            fs=[fs; i j];
        end
    end
end

nbf=length(fs);
[i,j]=sort(fs(:,1)./fs(:,2));   % tri des fractions dans l'ordre de grandeur

% traitement des styles 
switch upper(style)
    case 'TOUS'
        fs=[0 1; fs(j,:);1 1]';
    case 'SANSUN'
        fs=fs(j,:)';
    case 'DEMI'
        fs=fs(j(1:round(nbf/2)),:)';
    otherwise
        help ford
end

% reprŽsentation des cercles positionnŽs en X ˆ num/dŽnominateur et dont le
% rayon est 1/2*dŽn^2.
C=ones(size(x));
nbs=size(fs,2);
f=linspace(0,1,nbs);
for k=1:size(fs,2)
    i=fs(1,k);
    j=fs(2,k);
    cx=i/j;
    r=1/(2*j^2);
    cy=r;
    surf(r*x+cx,r*y+cy,r*z,f(k)*C)
    surf(r*x+cx,r*y-cy,r*z,f(k)*C)
end

axis equal
axis off
hold off
shading interp
light
lighting phong
material shiny
view(45,45)
colormap(coul)
figure(gcf)
