function fordsphere(nb)
% FORDSPHERE(NB) reprŽsentent les sph�res d'apr�s l'idŽe de L.R. Ford.
% Le principe consiste ˆ prendre tous les fractions complexes
% entre -nb et nb: p/q o� p= j + i k et q = n+ i m.
% Une sph�re est reprŽsentŽe sur le plan complexe ˆ la position
% p/q et le rayon est 1/(2qq') o' q' est le complexe conjuguŽ de q.

% J.-B. Roux, fŽvrier 2003

c=zeros((2*nb+1)^4,3);  % matrice des centres des sph�res
cc=c;               % matrice des indices de couleur
t=1;                % compteur du nombre de sph�res
i_coul=1;           % indice de la couleur
% p = j + i k
% q = n + i m
for j=-nb:nb
    for k=-nb:nb
        for n=-nb:nb
            for m=-nb:nb
                if (n^2 + m^2)~=0   % Žviter les divisions par 0
                    cx=(j*n + k*m)/(n^2 + m^2);
                    cy=(k*n - j*m)/(n^2 + m^2);
                    % Žviter de tomber sur les m�mes sph�res
                    % = rŽduction de la fraction p/q (1/2 = 2/4)
                    li=find(c(:,1)==1|c(:,2)==2);   
                    if isempty(li)
                        if sqrt(cx^2+cy^2)<1    % p/q est < 1
                            cz=1/(2*(n^2+m^2));
                            c(t,:)=[cx cy cz];
                            lc=find(c(1:t-1,3)==cz);
                            if isempty(lc)
                                cc(t)=i_coul;
                                i_coul=i_coul+1;
                            else
                                cc(t)=cc(lc(1));
                            end
                            t=t+1;
                        end
                    end
                end
            end
        end
    end
end
nbs=t-1;                % nombre total de sph�res
c=c(1:nbs,:);

% reprŽsentation graphique
h=plot3(0,0,0);
hold on
delete(h)
[x,y,z]=sphere;

nbc=i_coul-1;           % nombre total de couleurs
sprintf('%d sph�res, %d couleurs \n',nbs,nbc)
cc=cc(1:nbs);           % on purge les zŽros
C=ones(size(x));        % matrice de couleurs pour une sph�re
f=linspace(0,1,nbc);    % couleurs des sph�res
for s=1:nbs
    r=c(s,3);
    if r<0.5
            surf(r*x+c(s,1),r*y+c(s,2),r*z+c(s,3),f(cc(s))*C)
    end
end

axis equal
axis off
hold off
colormap(prism(nbc))    % palette de couleur
light                   % Žclairage
lighting phong          % diffusion
material shiny          % matŽriel de la surface
shading interp          % ombrage interpolŽ
view(0,90)              % on observe de dessus
figure(gcf)
