% - currentpoint ------------------------

/pathcount {
0 
{pop pop 1 add}
{pop pop 1 add}
{pop pop pop pop pop pop 1 add}
{1 add}
pathforall
} def

/thereisacurrentpoint {
0
pathcount 
ne 
{true}
{false}
ifelse
} def

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

/pwd 0.032 def

/pixel {
2 copy
gsave
% display coordinates:
% 2 copy exch == == (===) ==
newpath 
moveto
pwd 2 div neg pwd 2 div neg rmoveto
pwd 0 rlineto
0 pwd rlineto
pwd neg 0 rlineto
0 pwd neg rlineto
closepath
fill
grestore
} def

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

% stack: [parameters] f t0 t1 N

/mkpath {
16 dict begin
/N exch def
/t1 exch def
/t0 exch def
/f exch cvx def
/pars exch def

/h t1 t0 sub N div def
/h3 h 0.333333 mul def

/currentloc pars t0 
f def
pars t0 f 0 get
aload pop 

thereisacurrentpoint
{
	lineto
}
{
	moveto
}
ifelse

N {				% x y = currentpoint
% currentpoint pixel pop

currentloc 0 get 0 get		% x0 dx0
currentloc 1 get 0 get
h3 mul
add

currentloc 0 get 1 get
currentloc 1 get 1 get
h3 mul 
add

/t0 t0 h add def
/currentloc pars t0 f def
currentloc 0 get 0 get
currentloc 1 get 0 get
h3 mul sub
currentloc 0 get 1 get
currentloc 1 get 1 get
h3 mul sub
currentloc 0 get 0 get
currentloc 0 get 1 get
% pixel
curveto
} repeat
end
} def

% an addition to the 2D drawing: produces a polygon
% with exactly the same call as mkpath - for debugging
% when you suspect your derivative formula is incorrect

/mkpolypath {
10 dict begin
/N exch def
/tf exch def
/ti exch def
/f exch cvx def
/pars exch def

/h tf ti sub N div def
pars ti f 		% [[x  y][x' y']]
0 get aload pop moveto
/t ti h add def
N {
pars t f 
0 get aload pop lineto
/t t h add def
} repeat
end
} def

% - testing ---------------------------------------

false {

/Pi 3.1415926536 def

% in:  [parameters] t 
% out: [[xt yt][dxt dyt]]

/circle {
1 dict begin
/t exch def
aload pop /R exch def
/t t 180 mul Pi div def
[
[ t cos R mul t sin R mul ]
[ t sin neg R mul t cos R mul ]
]
end
} def

72 72 scale
4.25 5.5 translate
0.02 setlinewidth
1 setlinecap
1 setlinejoin

newpath
[2] /circle 0 2 Pi mul 8 drawpath
closepath
stroke

} if