function [V,C,T] = sphem(A,type)
% Assume that A exists and is a list of points
% on the surface of a unit sphere.

pagedim = [612,792];
scale = 100;

clear V C T;
% as in the voronoin help page

[V C] = voronoin([ 0 0 0 ; A ]);

V = V(2:length(V),:);
C = C(2:length(C));

for i = 1:length(C),
  C{i} = C{i}(2:length(C{i}));
  C{i} = C{i} - ones(size(C{i}));
end

T = delaunayn([ 0 0 0 ; A ]);

for i = 1:length(T),
  T(i,:) = sort(T(i,:));
end

T = T(:,2:4);
T = T - ones(size(T));

% we'd like to match triangulations to vertices
% for each triangulation, we take the cells, and find the vertex
% that they each have in common
for i = 1:length(T),
  clear s;
  s = sort([C{T(i,1)} C{T(i,2)} C{T(i,3)}]);
  s = s(find([s 0 0] == [0 0 s]));
  T(i,4) = s(1);
  % we're guaranteed to have a single triplicate value
  
  if (length(s) > 1),
    error('unstable data');
  end
end

T = sortrows(T,4);
T = T(:,1:3);

for i = 1:length(T),
  if (dot(cross(A(T(i,1),:),A(T(i,2),:)),A(T(i,3),:)) < 0)
    T(i,:) = [T(i,1) T(i,3) T(i,2)];
  end
end


for i = 1:length(C),
  clear s;
  s = (i == T(C{i}(1),:));
  s = [s(3) s(1:2)] * T(C{i}(1),:)';
  C{i} = [C{i} ; zeros(size(C{i}))];
  C{i}(2,1) = s;
  for j = 2:length(C{i}),
    for k = j:length(C{i}),
      clear s;
      s = sort([C{i}(2,j-1) i T(C{i}(1,k),:)]);
      if sum([s -1] == [-1 s]) == 2
	if (j ~= k)
          C{i} = [C{i}(:,1:j-1) C{i}(:,k) C{i}(:,j+1:k-1) C{i}(:,j) C{i}(:,k+1:length(C{i}))];
	end
	% we can now dispense with k, because we switched columns j and k
        C{i}(2,j) = ((C{i}(2,j-1) ~= T(C{i}(1,j),:)) .* T(C{i}(1,j),:)) * (i ~= T(C{i}(1,j),:))';
      end
    end
  end
end

clear s;
for i = 1:length(V),
  s = sqrt(dot(V(i,:),V(i,:)));
  V(i,:) = V(i,:)/s;
end

for i = 1:length(C),
  if dot(V(C{i}(1,1),:),V(C{i}(1,2),:)) > dot(V(C{i}(1,2),:),V(C{i}(1,3),:))
    C{i} = [C{i}(:,2:length(C{i})) C{i}(:,1)];
  end
end

fid = fopen('out.ps','w');

psheader(fid);

switch type
  case 'voronoi',
    drawarc(fid); cshow(fid);
    pageheader(fid);
    fprintf(fid,'100 690 moveto\n');
    for i = 1:length(C),
      clear s;
      s = [C{i} C{i}];
      fprintf(fid,'0 -20 rmoveto (%d)  ',i);
      for j = 1:length(C{i})
        fprintf(fid,'%f (%d) ', acos(dot(V(s(1,j),:),V(s(1,j+1),:)))*180/pi, s(2,j));
      end
      fprintf(fid,'%d drawarc\n', length(C{i}));
    end
  case 'voronoi2',
    drawarc2(fid); cshow(fid); cshow2(fid);
    pageheader(fid);
    fprintf(fid,'100 690 moveto\n');
    for i = 1:length(C),
      clear s;
      s = [C{i} C{i}];
      fprintf(fid,'0 -20 rmoveto (%d) [ ',i);
      for j = 1:length(C{i})
        fprintf(fid,'%f ', acos(dot(V(s(1,j),:),V(s(1,j+1),:)))*180/pi);
      end
      fprintf(fid,'] [ ',i);
      for j = 1:length(C{i})
        fprintf(fid,'(%d) ', s(2,j));
      end
      fprintf(fid,'] %d drawarc\n', length(C{i}));
    end
  case 'delaunay',
    drawarcendlabels(fid); cshow(fid);
    pageheader(fid);
    fprintf(fid,'100 690 moveto\n');
    for i = 1:length(T)
      fprintf(fid,'0 -20 rmoveto ');
      s = [ T(i,:) T(i,:) ];
      for j = 1:3
        fprintf(fid,'%f (%d) ', acos(dot(A(s(j)),A(s(j+1))))*180/pi, s(j+1));;
      end
      fprintf(fid,'3 drawarc\n');
    end
  case 'plotdelaunay',
    U = '';
    hold off
    newplot;
    hold on;
    axis([-1.2 1.2 -1.2 1.2 -1.2 1.2]);
    for i = 1:length(T)
      s = sort(T(i,:));
      U = [U ; s(1:2); s(2:3); s(1) s(3)];
    end
    U = sortrows(U);
    for i = 1:2:length(U),
      P = [linspace(A(U(i,1),1), A(U(i,2),1), 50)', linspace(A(U(i,1),2), A(U(i,2),2), 50)', linspace(A(U(i,1),3), A(U(i,2),3), 50)'];
      for j = 1:length(P),
        s = sqrt(dot(P(j,:),P(j,:)));
	P(j,:) = P(j,:)/s;
      end
      plot3(P(:,1),P(:,2),P(:,3));
    end,
  case 'plotvoronoi',
    U = '';
    hold off
    newplot;
    hold on;
    axis([-1.2 1.2 -1.2 1.2 -1.2 1.2]);
    for i = 1:length(C)
      s = [C{i}(1,:) C{i}(1,1)];
      while length(s) > 1,
        U = [U ; sort(s(1:2))];
	s = s(2:length(s));
      end
    end
    U = sortrows(U);
    for i = 1:2:length(U),
      P = [linspace(V(U(i,1),1), V(U(i,2),1), 50)', linspace(V(U(i,1),2), V(U(i,2),2), 50)', linspace(V(U(i,1),3), V(U(i,2),3), 50)'];
      for j = 1:length(P),
        s = sqrt(dot(P(j,:),P(j,:)));
	P(j,:) = P(j,:)/s;
      end
      plot3(P(:,1),P(:,2),P(:,3));
    end
  case 'voronoiellipse',
    cshow(fid); cshow2(fid);
    z = 0.9;
    As = [ A(:,1:2) z*A(:,3) ];
    Vs = [ V(:,1:2) z*V(:,3) ];
    pageheader(fid);
    pagepos = [scale+15, pagedim(2)-scale-40 ];
    for i = 1:length(C),
      e = 0.7;
      s = [C{i} C{i}];  % the usual
      fprintf(fid,'gsave %% edge %d\n%d %d translate\n', i, pagepos(1), pagepos(2));
      for j = 1:length(C{i})+1
	c = cross(Vs(s(1,j),:),Vs(s(1,j+1),:));
	d = cross(c,[-c(2),c(1),0]);
	d(3) = d(3)/z;			% we can't normalise to the ellipse easily
	d = d/sqrt(dot(d,d));		% except by scaling it down to the circle,
	ds = [d(1) d(2) d(3)*z];	% scaling it, then scaling the z coordinate again
	% ds is the long axis of the ellipse in the plane that our two points are on
	% a1 and a2 are angles, t? is a temporary vector
        a1 = acos(dot(V(s(1,j),:),d));
        a2 = acos(dot(V(s(1,j+1),:),d));
        t = cross(V(s(1,j),:),V(s(1,j+1),:));
        ta = cross(d,V(s(1,j),:));
        tb = cross(d,V(s(1,j+1),:));
        if (ta(3)*tb(3) > 0) 
	  if (t(3)*ta(3) < 0)
            d = -d;
            ds = -ds;
            a1 = pi - a1;
            a2 = pi - a2;
            % if the segment doesn't cover an end, then move it so that it's to the top
	  end
        else
          if (a1+a2 > pi)
            d = -d;
            ds = -ds;
            a1 = pi - a1;
            a2 = pi - a2;
            % if the segment does cover an end, then move it so that
            % it's across the right half of the x axis
          end
          a1 = -a1;
        end
        as1 = atan2(sin(a1), sqrt(dot(ds,ds))*cos(a1));
        as2 = atan2(sin(a2), sqrt(dot(ds,ds))*cos(a2));
        fprintf(fid,'gsave\n%f rotate\n', (e-as1)*180/pi);
        e = e + as2 - as1;
        fprintf(fid,'%f 1 scale\n',sqrt(dot(ds,ds)));
        if (j == length(C{i})+1)
          fprintf(fid,'0 0 r2 %f %f arcn\n', ((2*a1+a2)/3-0.03)*180/pi, a1*180/pi);
          fprintf(fid,'0 0 r1 %f %f arc\nclosepath stroke\n', a1*180/pi, ((2*a1+a2)/3+0.03)*180/pi);
        elseif (j == 1)
          fprintf(fid,'0 0 r2 %f %f arcn\n', a2*180/pi, ((2*a1+a2)/3-0.03)*180/pi);
          fprintf(fid,'0 0 r1 %f %f arc\nstroke\n', ((2*a1+a2)/3+0.03)*180/pi, a2*180/pi);
	  % this is an edge label
          fprintf(fid,'%f r3 mul %f r3 mul moveto (%d) cshow\n', cos((a1+a2)/2), sin((a1+a2)/2), s(2,j));
	else
          fprintf(fid,'0 0 r2 %f %f arcn\n', a2*180/pi, a1*180/pi);
          fprintf(fid,'0 0 r1 %f %f arc\nstroke\n', a1*180/pi, a2*180/pi);
	  % this is an edge label
          fprintf(fid,'%f r3 mul %f r3 mul moveto (%d) cshow\n', cos((a1+a2)/2), sin((a1+a2)/2), s(2,j));
        end
	if (a2-a1 > 0.15 & j > 1)
	  % This label is massively redundant.  Don't bother with it, if the edge is really short.
          fprintf(fid,'%f r3 mul %f r3 mul moveto (%d) cshow2\n', cos(a1+0.03), sin(a1+0.03), i);
	end
        fprintf(fid,'grestore\n');
        if (j == length(C{i})+1)
	  % this is the x that goes in the center
          fprintf(fid,'r4 2 div dup dup moveto dup -1 mul dup lineto dup -1 mul 2 copy moveto exch lineto stroke\n');
        end
      end  % for loop incrementing j
      fprintf(fid,'grestore\n\n');
      if (pagepos(2) <= scale-30)
        pagepos(2) = pagedim(2)-scale-40;
        if (pagepos(1) + 2*scale + 30 > pagedim(1) - scale - 15)
          fprintf(fid,'showpage\n\n');
	  pageheader(fid);
          pagepos(1) = scale+15;
        else
          pagepos(1) = pagepos(1) + 2*scale + 30;
        end
      else
        pagepos(2) = pagepos(2) - 50;
      end
    end  % for loop incrementing i
    fprintf(fid,'\nshowpage\n');
    % end of voronoiellipse',
  case 'delaunayellipse',
    cshow(fid);
    z = 1.3;
    As = [ A(:,1:2) z*A(:,3) ];
    Vs = [ V(:,1:2) z*V(:,3) ];
    pageheader(fid);
    pagepos = [scale+40, pagedim(2)-scale-40 ];
    for i = 1:length(T),
      e = 0.4;
      s = [T(i,:) T(i,:)];  % the usual
      fprintf(fid,'gsave %% edge %d\n%d %d translate\n', i, pagepos(1), pagepos(2));
      for j = 1:4
	c = cross(As(s(1,j),:),As(s(1,j+1),:));
	d = cross(c,[-c(2),c(1),0]);
	d(3) = d(3)/z;			% we can't normalise to the ellipse easily
	d = d/sqrt(dot(d,d));		% except by scaling it down to the circle,
	ds = [d(1) d(2) d(3)*z];	% scaling it, then scaling the z coordinate again
	% ds is the long axis of the ellipse in the plane that our two points are on
	% a1 and a2 are angles, t? is a temporary vector
        a1 = acos(dot(A(s(1,j),:),d));
        a2 = acos(dot(A(s(1,j+1),:),d));
        t = cross(A(s(1,j),:),A(s(1,j+1),:));
        ta = cross(d,A(s(1,j),:));
        tb = cross(d,A(s(1,j+1),:));
        if (ta(3)*tb(3) > 0) 
	  if (t(3)*ta(3) < 0)
            d = -d;
            ds = -ds;
            a1 = pi - a1;
            a2 = pi - a2;
            % if the segment doesn't cover an end, then move it so that it's to the top
	  end
        else
          if (a1+a2 > pi)
            d = -d;
            ds = -ds;
            a1 = pi - a1;
            a2 = pi - a2;
            % if the segment does cover an end, then move it so that
            % it's across the right half of the x axis
          end
          a1 = -a1;
        end
        as1 = atan2(sin(a1), sqrt(dot(ds,ds))*cos(a1));
        as2 = atan2(sin(a2), sqrt(dot(ds,ds))*cos(a2));
        fprintf(fid,'gsave\n%f rotate\n', (e-as1)*180/pi);
        e = e + as2 - as1;
        fprintf(fid,'%f 1 scale\n',sqrt(dot(ds,ds)));
        if (j == 4)
          fprintf(fid,'0 0 r2 %f %f arcn\n', ((2*a1+a2)/3-0.03)*180/pi, a1*180/pi);
          fprintf(fid,'0 0 r1 %f %f arc\nclosepath stroke\n', a1*180/pi, ((2*a1+a2)/3+0.03)*180/pi);
	elseif (j == 1)
          fprintf(fid,'0 0 r2 %f %f arcn\n', a2*180/pi, ((2*a1+a2)/3-0.03)*180/pi);
          fprintf(fid,'0 0 r1 %f %f arc\nstroke\n', ((2*a1+a2)/3+0.03)*180/pi, a2*180/pi);
	else
          fprintf(fid,'0 0 r2 %f %f arcn\n', a2*180/pi, a1*180/pi);
          fprintf(fid,'0 0 r1 %f %f arc\nstroke\n', a1*180/pi, a2*180/pi);
        end
        if (j ~= 4)
          fprintf(fid,'%f r3 mul %f r3 mul moveto (%d) cshow\n', cos(a2+0.08), sin(a2+0.08), s(j+1));
          fprintf(fid,'%f r3 mul %f r3 mul moveto (%d) cshow\n', cos(a2-0.08), sin(a2-0.08), s(j+1));
        end
        fprintf(fid,'grestore\n');
        if (j == 4)
          fprintf(fid,'r4 2 div dup dup moveto dup -1 mul dup lineto dup -1 mul 2 copy moveto exch lineto stroke\n');
        end
      end  % for loop incrementing j
      fprintf(fid,'grestore\n\n');
      if (pagepos(2) <= scale-40)
        pagepos(2) = pagedim(2)-scale-40;
        if (pagepos(1) + 2*scale + 30 > pagedim(1) - scale - 15)
          fprintf(fid,'showpage\n\n');
	  pageheader(fid);
          pagepos(1) = scale+40;
        else
          pagepos(1) = pagepos(1) + 2*scale + 30;
        end
      else
        pagepos(2) = pagepos(2) - 50;
      end
    end  % for loop incrementing i
    fprintf(fid,'\nshowpage\n');
    % end of delaunayellipse case
  otherwise,
    disp 'unknown plot type'
end

fclose(fid);

%---------------------------
function [] = psheader(fid)
fprintf(fid,'%%!ps\n\n');
fprintf(fid,'%% Generated by $Id: sphem.m,v 1.11 2002/04/14 02:40:53 cfox Exp $ on %s\n\n', date);

fprintf(fid,'/unit %d def /r1 0.90 unit mul def /r2 unit def /r3 0.95 unit mul def /r4 0.04 unit mul def\n', 100);
fprintf(fid,'/Courier findfont 1 r4 mul scalefont setfont\n\n');

function [] = pageheader(fid)
fprintf(fid,'0.1 setlinewidth\n');
%fprintf(fid,'100 690 moveto\n');

function [] = drawarc(fid)
fprintf(fid,'/drawarc { dup 2 mul 1 add copy 0 exch {exch pop add} repeat dup 2 div\n');
fprintf(fid,'dup 90 exch sub exch 90 add currentpoint 4 copy newpath r1 5 -2 roll\n');
fprintf(fid,'exch arcn 4 copy r2 5 -2 roll arc 3 -1 roll pop 4 -1 roll pop 4 -1 roll\n');
fprintf(fid,'{ 5 copy 5 -1 roll 2 div 4 -1 roll add dup sin r3 mul exch cos r3 mul 4\n');
fprintf(fid,'-1 roll add 3 1 roll add moveto cshow 4 -1 roll pop 4 2 roll add 3 copy\n');
fprintf(fid,'dup sin exch cos 4 copy r1 mul 4 -1 roll add 3 1 roll r1 mul add moveto\n');
fprintf(fid,'r2 mul 4 -1 roll add 3 1 roll r2 mul add lineto 3 1 roll } repeat 3 2\n');
fprintf(fid,'roll 3 add dup cos r3 mul exch sin r3 mul 4 copy pop pop 4 -1 roll 3 -1\n');
fprintf(fid,'roll add 3 1 roll add moveto 3 -1 roll cshow stroke moveto } def\n\n');

function [] = drawarcendlabels(fid)

fprintf(fid,'/drawarc { dup 2 mul 1 add copy 0 exch {exch pop add} repeat dup 2 div\n');
fprintf(fid,'dup 90 exch sub exch 90 add 10 add currentpoint 4 copy newpath r1 5 -2\n');
fprintf(fid,'roll exch arcn 4 copy r2 5 -2 roll arc  3 -1 roll pop 4 -1 roll pop 4\n');
fprintf(fid,'-1 roll { 5 -1 roll 4 -1 roll add 4 copy 3 add dup sin r3 mul exch cos\n');
fprintf(fid,'r3 mul 4 -1 roll add 3 1 roll add moveto cshow 4 copy 3 sub dup sin r3\n');
fprintf(fid,'mul exch cos r3 mul 4 -1 roll add 3 1 roll add moveto cshow 4 -1 roll\n');
fprintf(fid,'pop 3 copy dup sin r1 mul exch cos r1 mul 4 -1 roll add 3 1 roll add\n');
fprintf(fid,'moveto 3 copy dup sin r2 mul exch cos r2 mul 4 -1 roll add 3 1 roll add\n');
fprintf(fid,'lineto 3 1 roll } repeat stroke moveto } def\n\n');

function [] = cshow(fid)

fprintf(fid,'/cshow { gsave /Courier findfont 1 r4 mul scalefont setfont 0.1\n');
fprintf(fid,'setlinewidth currentpoint 2 copy newpath 2 copy exch r4 add exch moveto\n');
fprintf(fid,'r4 0 360 arc 3 copy pop pop stringwidth pop -2 div 3 -1 roll add exch\n');
fprintf(fid,'-0.5 r4 mul add moveto show closepath stroke grestore } bind def\n\n');

function [] = cshow2(fid)

fprintf(fid,'/cshow2 { gsave /Courier-Bold findfont 0.8 r4 mul scalefont setfont\n');
fprintf(fid,'currentpoint 3 copy pop pop stringwidth pop -2 div 3 -1 roll add exch\n');
fprintf(fid,'-0.5 r4 mul add moveto show grestore } bind def\n\n');

function [] = drawarc2(fid)
% unit-label [ x1 x2 ... xn ] [ 1l l2 ... ln ] n 
fprintf(fid,'/drawarc { /n exch def /labels exch def /angles exch def /mylabel exch\n');
fprintf(fid,'def currentpoint /cy exch def /cx exch def 0 angles { add } forall 2\n');
fprintf(fid,'div dup dup 90 add 2 add exch 90 exch sub 2 add cx cy r2 4 -1 roll 5 -1\n');
fprintf(fid,'roll newpath arc dup dup 90 add 2 sub exch 90 exch sub 2 sub cx cy r1 5\n');
fprintf(fid,'-1 roll 5 -1 roll arcn closepath stroke 90 exch sub angles 0 get 3 div\n');
fprintf(fid,'sub 0 1 n 1 sub { dup angles exch get 3 copy 2 div exch pop add dup cos\n');
fprintf(fid,'r3 mul cx add exch sin r3 mul cy add moveto exch labels exch get cshow\n');
fprintf(fid,'add dup dup sin exch cos 2 copy r2 mul cx add exch r2 mul cy add moveto\n');
fprintf(fid,'r1 mul cx add exch r1 mul cy add lineto stroke dup 2 add dup cos r3 mul\n');
fprintf(fid,'cx add exch sin r3 mul cy add moveto mylabel cshow2 } for cx cy moveto } def\n\n');