function [Bout,Vout] = bucky
%BUCKY	Connectivity graph of the Buckminister Fuller geodesic dome.
%	B = BUCKY is the 60-by-60 sparse adjacency matrix of the
%	    connectivity graph of the geodesic dome, the soccer ball,
%	    and the carbon-60 molecule.
%	[B,V] = BUCKY also returns xyz coordinates of the vertices.
%
%	BUCKY, with no output arguments, demonstrates itself.

%	C. B. Moler, 2-14-91, 10-8-91, 8-11-92.
%	Copyright (c) 1984-93 by The MathWorks, Inc.

% These all-important parameters were generated by numeric
% search in another program.  I wish I had formulas for them.

alfa = [1.220391131161 0.814053008553 0.540911992917 0.172493462472];
beta = [0.310919101124 0.164146934496];

% Northern hemisphere
B = sparse(30,30);
phi = zeros(30,1);
theta = zeros(30,1);
j = 0;
k = 5;
for t = 0:4
   j = j+1;
   phi(j) = alfa(1);
   theta(j) = t*2*pi/5;
   B(j,k+1) = 1;
   B(j,rem(j,5)+1) = 1;

   k = k+1;
   phi(k) = alfa(2);
   theta(k) = t*2*pi/5;
   B(k,k+1) = 1;
   B(k,k+4) = 1;

   k = k+1;
   phi(k) = alfa(3);
   theta(k) = (t-beta(1))*2*pi/5;
   B(k,k+1) = 1;

   k = k+1;
   phi(k) = alfa(4);
   theta(k) = (t-beta(2))*2*pi/5;
   B(k,k+1) = 1;

   k = k+1;
   phi(k) = alfa(4);
   theta(k) = (t+beta(2))*2*pi/5;
   B(k,k+1) = 1;

   k = k+1;
   phi(k) = alfa(3);
   theta(k) = (t+beta(1))*2*pi/5;
   if k+2 < 30, B(k,k+2) = 1;
   else B(k,7) = 1; end
end
B = B + B';

% Reflect into southern hemisphere
k = 60:-1:31;
phi(k) = -phi;
theta(k) = pi+theta;
B(k,k) = B;

% Connect the two hemispheres
k = find(sum(B)==2);
j = k(1:10);
k = k([15:-1:11 20:-1:16]);
c = sparse(j,k,ones(size(j)),60,60);
B = B + c + c';

% Generate the 3-D coordinates
x = cos(phi).*cos(theta);
y = cos(phi).*sin(theta);
z = sin(phi);
V = [x y z];

if nargout == 2
   Bout = B;
   Vout = V;
elseif nargout == 1
   Bout = B;
else
   help bucky
   disp(' ')
   disp('Press any key to continue after pauses.')
   disp('pause'), pause

   disp(' ')
   disp('The entire graph, with no numbering specified.')
   disp(' ')
   H = sparse(60,60);
   k = 31:60;
   H(k,k) = B(k,k);
   gplot(H,V,'m-')
   axis('equal');
   hold('on')
   gplot(B-H,V,'c-');
   hold('off')
   disp('pause'), pause

   disp(' ')
   disp('Number the nodes in one hemisphere, pentagon by pentagon.')
   disp(' ')
   gplot(B(1:30,1:30),[x(1:30) y(1:30)])
   for j = 1:30
      text(x(j),y(j),int2str(j));
   end
   axis('equal');
   disp('pause'), pause

   disp(' ')
   disp('''Spy'' plot of the resulting sparse matrix.')
   disp(' ')
   spy(B(1:30,1:30))
   disp('pause'), pause

   disp(' ')
   disp('Reflect numbering of one hemisphere into the other.')
   disp(' ')
   % Map northern hemisphere onto 0 <= r <= 1 
   % and southern hemisphere onto 1 <= r <= s.
   s = 1.65;
   r = x.*x + y.*y;
   k = find(z < 0);
   x(k) = x(k)./r(k);
   y(k) = y(k)./r(k);
   r = sqrt(x.*x + y.*y);
   r = min(s,r)./r;
   x = r.*x;
   y = r.*y;
   gplot(B,[x y])
   for j = 1:60
      text(x(j),y(j),int2str(j));
   end
   axis('equal');
   axis([-s s -s s]);
   disp('pause'), pause

   disp(' ')
   disp('''Spy'' plot of the final sparse matrix.')
   disp(' ')
   spy(B)

   disp('End')
end
