%!

% longitude latitude -> point on unit sphere
/P { 1 dict begin
  /l exch def
  /L exch def
  % L=longitude l=latitude
  [L cos l cos mul 
   L sin l cos mul 
   l sin ]
end } def

% N -> latitudes in N+1 rows
/sphere-vertex { 1 dict begin
/N exch def

/dA 180 N div def
/dB 360 N div 2 div def

% A = latitude
[
  [ 
    [0 0 -1]
  ]
/A -90 dA add def
N 1 sub { 
  [
    /B 0 def
    % B = longitude
    2 N mul 1 add {
      B A P 
      /B B dB add def
    } repeat
  ]
  /A A dA add def
  } repeat
  [
    [0 0 1]
  ]
]
end } def

% N -> array of faces
/sphere { 1 dict begin
/N exch def

/S N sphere-vertex def
% S now is an array of vertices, arranged in latitudes
[
  % the triangles at the south pole
  0 1 2 N mul 1 sub { /j exch def
    [
      [ S 0 get 0 get 
        S 1 get j 1 add get 
	    S 1 get j get 
      ] dup normal-function
    ]
  } for
  % the rectangular regions in the middle
  1 1 N 2 sub { /i exch def
    0 1 2 N mul 1 sub { /j exch def
      [
	    [ 
	      S i get j get 
	      S i get j 1 add get
          S i 1 add get j 1 add get 
          S i 1 add get j get 
	    ] dup normal-function 
      ] 
    } for
  } for
  % the triangles at the north pole
  0 1 2 N mul 1 sub { /j exch def
    [
      [ S N 1 sub get j get 
        S N 1 sub get j 1 add get 
        S N get 0 get
      ] dup normal-function
    ]
  } for
] 
end } def