% renders surfaces in 3d

/pi 3.141593 def

% needs matrix3d.inc

% a surface is an array of faces
% a face = [ polygon + normal ]
% a polygon = array of vertices

% affine transformation A + a surface s => transformed surface

/surface-transform {
4 dict begin
/s exch def
/A exch def
[
  s {
  /f exch def
  /p f 0 get def
  /n f 1 get def
  [
    [
      p {
        A exch affine-transform-3d
      } forall
    ]
    A 0 get n transform-3d 
  ]
  } forall
]
end
} def

% face

/is-visible-in-projection {
% tests normal has a z-component non-negative
1 get 2 get 0 ge {
  true
} {
  false
} ifelse
} def

% face

/is-visible-in-perspective {
1 dict begin
/f exch def
  f 0 get 0 get 
  f 1 get
  dot-product-3d 
% tests whether v0*normal is non-positive
% normal pointing towards you
0 le {
  true
} {
  false
} ifelse
end
} def

% [x y z] => [x y]

/projection {
1 dict begin
/v exch def
[v 0 get v 1 get]
end
} def

% [x y z] => [-x/z -y/z]

/perspective {
2 dict begin
/v exch def
/z v 2 get neg def
[v 0 get z div v 1 get z div]
end
} def

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

(matrix3d.inc) run

/light-source [-1 1 0.5] normalize-3d def

% normal

/minshade 0.05 def
/maxshade 0.95 def
/shadefactor maxshade minshade sub 4 div def

/shading {
2 dict begin
/n exch def
/d n light-source dot-product-3d 
def
maxshade minshade add 0.5 mul
maxshade minshade sub 0.5 mul 
d mul add
end
} def

% pars /f [si ti] [tf sf] [ns nt] => surface parametrized by f 
% and range [si sf] x [ti tf] - get ns x nt rectangles
% f: pars [s t] => f(s, t)

/mksurface {
16 dict begin
aload pop
/nt exch def
/ns exch def
% [ns nt] ==
aload pop
/tf exch def
/sf exch def
aload pop
/ti exch def
/si exch def
/f exch cvx def
/pars exch def
/sinc sf si sub ns div def
/tinc tf ti sub nt div def
% first make an array of (ns + 1) x (nt + 1) points
[
  /t ti def
  nt 1 add {
  [
    /s si def
    ns 1 add {
	  % (s, t = ) print [s t] ==
      % make the f(s, t) vertex
      pars [s t] f 
      /s s sinc add def
    } repeat
    /t t sinc add def
	]
  } repeat
] /vertex exch def

[
0 1 nt 1 sub {
  /i exch def
  /v0 vertex i get def
  /v1 vertex i 1 add get def
  0 1 ns 1 sub {
	/j exch def
    /J j 1 add def
	% make up face[i, j]
	[
	% list of vertices
	[v0 j get 
	 v0 J get
	 v1 J get
     v1 j get
	 v0 j get ]

	% now the normal vector
	% v1[j] v1[j+1]
	% v0[j] v0[j+1]
	/n 
	  v0 J get
	  v0 j get
	  vector-sub-3d 
	  v1 j get
	  v0 j get
	  vector-sub-3d
	  cross-product	
	def
	/r n vector-length-3d def
	r 0 eq {
	  /n 
	    v1 j get
	    v1 J get
	    vector-sub-3d 
	    v0 J get
	    v1 J get
	    vector-sub-3d
	    cross-product	
	  def
	  /r n vector-length-3d def
	  r 0 eq {
		/n [ 0 0 1] def
		/r 1 def
	  } if
	} if
	[n 0 get r div n 1 get r div n 2 get r div]
	]
  } for

} for
]
end 
} def

% - end of draw3d.inc -------------------------------

false {

72 dup scale
4.25 5.5 translate
0.01 setlinewidth

% [pars] [s t] 

/sphere {
4 dict begin
% t=latitude, s = longitude
aload pop
/t exch 180 mul pi div def
/s exch 180 mul pi div def
aload pop
/R exch def
[
	s cos t cos mul R mul
	s sin t cos mul R mul
	t sin R mul
]
end
} def

[1] /sphere [0 pi -0.5 mul][2 pi mul pi 0.5 mul] [ 4 4 ] mksurface 
/s exch def

s {
/f exch def

} forall 

} if