%!
%%BoundingBox: 0 0 500 500

(ps3d.inc) run

/P {
2 dict begin
/B exch def
/A exch def
% lat long
[A cos B cos mul A sin B cos mul B sin]
end
} def

/page-begin {
gsave
72 dup scale
1 72 div setlinewidth
4 5 translate
1 setlinecap
1 setlinejoin
} def

/page-end {
grestore
showpage
} def

% latitudes in rows 
/sphere-vertex { 1 dict begin
/N exch def
/dA 179.8 N div def
/dB 360 N div 2 div def
/A -89.9 def
% A = latitude
N 1 add { 
  [
  /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
end } def

/N 16 def
/S [ N sphere-vertex ] def

/sphere [
0 1 N 1 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 ] dup normal-function ] 
    [
      [ S i get j get 
      S i 1 add get j 1 add get 
      S i 1 add get j get ] dup normal-function ] 
  } for
} for
] def

% ------------------------------------------------------------

[0 0 6 1] set-eye 
/light [-1 2 1 0 ] normalized def
/Sh [ 0 0.2 0.8 0.9 ] def

/A 0 def
/dA 10 def
/E get-eye def

 {

page-begin
gsave3d

[0 0 1] 90 rotate3d
[0 1 0] 90 rotate3d
[0 0 1] A rotate3d
/Tinv ctm3d 1 get def
/T ctm3d 0 get def

/O E Tinv transform3d def
/L light Tinv transform3d def

0 1 sphere length 1 sub {
  /i exch def
  /f sphere i get def
  % f = [ point-array visibility ]
  /v f 1 get def
  v O dot-product 0 ge {
     /p f 0 get def

/ds [
	0 [ p 0 get aload pop 1] T transform3d render p 0 get L dot-product Sh exch shade
	0 [ p 1 get aload pop 1] T transform3d render p 1 get L dot-product Sh exch shade
	0 [ p 2 get aload pop 1] T transform3d render p 2 get L dot-product Sh exch shade
] def
gsave
<<
	/ShadingType 4
	/ColorSpace [ /DeviceGray ]
	/DataSource ds
>>
shfill
grestore
 
  } if
} for

grestore3d
page-end

/A A dA add def

} loop