%%% Date: Mon, 12 Oct 92 22:32:32 EST
%%% Message-Id: <9210121232.AA18641@ee.latrobe.edu.au>
%%% To: ajs@msdrl.com
%%% Subject: graphbase.mf 0.2 gt mod 1 Mon 12 Oct 1992.
%%% Cc: ecsgrt@ee.latrobe.edu.au
%%% Status: OR
%%% 
%%% 12:10 GMT Mon 12 Oct 1992 - Geoffrey Tobin.
%%% 
%%% Dear Anthony,
%%% 
%%% Latest mod of "graphbase.mf".
%%% 
%%% This is still a long way from what it ought to be for Fig 2.1, but
%%% maybe you can use some of the new code.  Much of it is comments!
%%% 
%%% I've modified onedot based on Knuth's "drawdot", written a hatchpath
%%% macro to complement shadepath, attempted a gray level formula to match
%%% human light sensitivity, some macros to draw as well as shade or
%%% hatch, and an arrowpath macro to draw an arrow for arbitrary path and
%%% to draw that path.
%%% 
%%% THe 3-point arc code is included, with arcthreearrow to draw an arrow
%%% as well.
%%% 
%%% Some code for drawing Interpolated Splines is included.  The macros to
%%% test are ispline, isplinearrow, and isplineshade.
%%% 
%%% Backward arrows are on my agenda.
%%% 
%%% All the Best!
%%% Geoffrey Tobin
%%% 
%%%
%%%  File: graphbase.mf
%%%

mode_setup;
message "mfpic version 0.2 graphbase - gt mod 1 - 10:16 GMT Mon 12 Oct 1992";

% set up local environment

def mfpicenv =
begingroup

% miscellaneous utilities

% gt - op_pair operates with "op"
% on both parts of a pair p.

save op_pair;

vardef op_pair (text op) (expr p) =
  (op (xpart p), op (ypart p))
enddef;

save floorpair, ceilingpair;

def floorpair = op_pair (floor) enddef;
def ceilingpair = op_pair (ceiling) enddef;

% gt - Should there be more error-checking,
% eg of types, in these utility routines?
% That would slow them down.

% gt - textpairs converts the text t into the
% array of n pairs, pts, that it contains.

save textpairs;

def textpairs (text t) (suffix pairs_, n_) =
 n_ := 0;
 for q=t:
  pairs_[incr n_] := q;
 endfor;
enddef;

% gt - Watch out!  Need to ensure that "p_", etc.,
% don't clash with any name in the passed text "t".
% That's a nasty error to trace!
%
% A name conflict between local variables and variables
% in a text parameter is especially likely in low-level
% utility macros, such as minpair, maxpair and corner.
%
% Unfortunately, we can't *ensure* it won't happen.
% So I appended the underscore to reduce the
% probability of that happening.
%
% Evidently that's why Knuth uses "u_" in "max" and
% "min" in "plain.mf".

% gt - corner may be used for finding
% a corner of the bounding box of the
% set of points listed in u and t.
% Other uses may be imaginable. (?)

save corner;

vardef corner (text xop) (text yop)
              (expr u)(text t) =
  save p_;
  pair p_;
  p_ := u;
  for q=t:
    p_ := (xop (xpart p_, xpart q),
	   yop (ypart p_, ypart q));
  endfor;
  p_
enddef;

% gt - bottom right, bottom right,
% top left, top right corners.

save blpair, brpair, tlpair, trpair;

def blpair = corner (min) (min) enddef;
def brpair = corner (max) (min) enddef;
def tlpair = corner (min) (max) enddef;
def trpair = corner (max) (max) enddef;

def minpair = blpair enddef;
def maxpair = trpair enddef;

% setup
% gt - sets the graphics coordinates.

save bounds,
  xneg,xpos,yneg,ypos;

def bounds(expr a,b,c,d) = 
 xneg:=a; 
 xpos:=b; 
 yneg:=c; 
 ypos:=d; 
enddef;

% conversion

save xconv;

def xconv(expr xvalue) = 
 ((xvalue-xneg)/(xpos-xneg))*w 
enddef;

save unxconv;

def unxconv(expr pvalue) = 
 ((pvalue/w)*(xpos-xneg)+xneg) 
enddef;

save yconv;

def yconv(expr yvalue) = 
 ((yvalue-yneg)/(ypos-yneg))*h 
enddef;

save ztr;

transform ztr;

save setztr;

def setztr =
 ztr:=identity
 shifted -(xneg,yneg) 
 xscaled (w/(xpos-xneg))
 yscaled (h/(ypos-yneg));
enddef;

% pen width
% in pixel coordinates

save penwd;

newinternal penwd;

% gt - initial value of penwd.

interim penwd := 0.5pt;

% arrowheads
% in pixel coordinates

% hdwdr = arrowhead's ratio of width to length,
% hdten = tension used in drawing its barbs.

save hdwdr, hdten;

newinternal hdwdr, hdten;

% gt - initial values of hdwdr, hdten.

interim hdwdr := 1;
interim hdten := 1;

% draw an arrowhead.

save head, p,side;

def head(expr front, back, width, t) =
 pair p[], side;
 side := (width/2) * 
   ((front-back) rotated 90);
 p1 := back + side;
 p2 := back - side;
 draw front{back-front}..tension t..p1;
 draw front{back-front}..tension t..p2;
enddef;

% draw an arrowhead of length hlen
% for a path f.

save headpath, p;

def headpath(expr f,hlen) =
 pair p[];
 p2:=point infinity of f; 
 p1:=direction infinity of f;
 if p1<>(0,0):
  head(p2,p2-(hlen*unitvector(p1)),
    hdwdr,hdten);
 fi;
enddef;

% shading and hatching routines
% in pixel coordinates

% gt - modified onedot based on
% plain metafont's "drawdot".
% Used in stipple shading.
%
% Note:
% currentpen_path, def_pen_path_,
% t_, and penspeck are defined
% in plain metafont ("plain.mf"
% or "plain.base").

save onedot;

def onedot(expr p)(suffix v) =
  if unknown currentpen_path:
    def_pen_path_
  fi;
  addto v
    contour currentpen_path
    shifted p.t_
    withpen penspeck
enddef;

% gt - draw path f in picture v.
% ("onepath" is the old "onedot",
% but f is intended to be a general path.)
% Used, eg, in hatching and in drawpaths.

save onepath;

def onepath (expr f) (suffix v) =
  addto v doublepath f
    withpen currentpen; 
enddef;

% gt - Paths must be continuous - I think
% - but using suffix, we can pass arrays
% of paths.
%
% My eventual goal is to do as much as
% feasible, and memory-affordable, in
% graphics coordinates, so we can rotate
% and otherwise transform sets of paths
% before drawing.

% gt - draw the n paths f[] in picture v.

save drawpaths;

def drawpaths (expr n) (suffix f, v) =
  for i=1 upto n:
    onepath (f[i], v);
  endfor;
enddef;

% clip picture v to interior of path f.

save clip;

vardef clip(expr f)(suffix v) =
 save vt;
 picture vt;
 vt:=v;
 cull vt keeping (1,infinity);
 addto vt contour f;
 cull vt keeping (2,infinity);
 vt
enddef;

% gt - find bounding box of path f.

save boundingbox, p;

def boundingbox (expr f) (suffix ll, ur) =
  ur := ll := point 0 of f; 
  pair p[];
  for i=0 upto length f:
    p0 := point i of f; 
    p1 := precontrol i of f; 
    p2 := postcontrol i of f;
    ll := minpair (ll, p0, p1, p2);
    ur := maxpair (ur, p0, p1, p2);
  endfor; 
enddef;

% gt - shading.

% gt - I'm not so happy with dot densities
% over a uniform range.
% Here's code to approximate what may
% be the human eye's light sensitivity.
%
% Mind you, this sort of stuff is done much
% faster in C.

save exp;

vardef exp (expr x) =
  mexp (256 * x)
enddef;

% graya scales the spacing sp;
% grayb scales the graylevel g.

save graya, grayb;

newinternal graya, grayb;

% initial values of gray parameters.

interim graya := 0.5 pt;
interim grayb := 3/20;

% setgraypars sets gray parameters.
% experimentation is recommended.

save setgraypars;

def setgraypars (expr a, b) =
  graya := a;
  grayb := b;
enddef;

% gt - grayspace gives the dot spacing
% for graylevel g.
%
% Not sure how this model performs.

save grayspace;

vardef grayspace (expr g) =
  if g <= 1:  % white
    infinity
  elseif g >= 21:  % black
    0
  else:  % gray
    graya / (1 - exp (-g * grayb))
  fi
enddef;

% gt - stipple upright box with lower left
% at ll, upper right at ur, in picture v;
% 2sp is dot spacing (rows offset by sp).
%
% NB: "stipple" means "shade with dots",
% if I understand my English dictionary.
%
% Thomas Leathrum devised the trick whereby
% the dots are arranged on a regular grid
% of mesh size sp by sp with the pixel
% origin as one crosspoint.  This ensures
% that objects shaded with the same stipple
% density may be cleanly overlaid.

save shadebox, sll, mn, m, n, twosp, p;

def shadebox (expr sp, ll, ur) (suffix v) =
  pair sll;
  sll:=sp*(ceilingpair(ll/sp));
  pair mn;
  mn:=floorpair((ur-sll)/sp);
  m:=xpart mn;
  n:=ypart mn;
  twosp:=2sp;
  v:=nullpicture;
  pair p[];
  p2:=sll;
  for i=0 upto m: 
    p3:=p2 if odd i: +(0,sp) fi;
    for j=0 upto n:
      if (not odd (i+j)):
        onedot (p3, v);
        p3:=p3+(0,twosp);
      fi;
    endfor;
    p2:=p2+(sp,0);
  endfor; 
enddef;

% stipple interior of closed path f;
% if spacing not positive, fill.

save shadepath, ll, ur, v;

def shadepath (expr sp, f) =
 if not cycle f: ;
 elseif sp<=0:
   fill f; 
 elseif sp < infinity: 
   pair ll, ur;
   boundingbox (f, ll, ur);
   picture v;
   shadebox (sp, ll, ur, v);
   addto currentpicture 
     also clip(f,v);
 fi;
enddef;

% gt - hatch an upright box in picture v,
% with line separation sep x sep.
%
% Notice the similarity to shadebox.

save hatchbox, llx, lly, urx, ury, sll,
     mn, m, n, f;

def hatchbox (expr sep, ll, ur) (suffix v) =
  llx := xpart ll;
  lly := ypart ll;
  urx := xpart ur;
  ury := ypart ur;
  pair sll;
  sll := sep * ceilingpair (ll/sep);
  pair mn;
  mn := floorpair ((ur-sll)/sep);
  m := xpart mn;
  n := ypart mn;
  v := nullpicture;
  path f;
  f := (xpart sll, lly)--(xpart sll, ury);
  for i=0 upto m:
    onepath (f, v);
    f := f translated (sep, 0);
  endfor;
  f := (llx, ypart sll)--(urx, ypart sll);
  for j=0 upto n:
    onepath (f, v);
    f := f translated (0, sep);
  endfor;
enddef;

save hatchpath, ll, ur, v;

def hatchpath (expr sep, f) =
 if not cycle f: ;
 elseif sep<=0:
   fill f; 
 elseif sep < infinity: 
   pair ll, ur;
   boundingbox (f, ll, ur);
   picture v;
   hatchbox (sep, ll, ur, v);
   addto currentpicture 
     also clip (f, v);
 fi;
enddef;

% gt - shading & hatching macros
% with a syntax like draw, fill,
% unfill and erase.
% sp, sep are in pixel coords,
% f in graphics coordinates;
% f is transformed transparently.

save shade;

def shade (expr sp) expr f =
  shadepath (sp, f transformed ztr);
enddef;

save hatch;

def hatch (expr sep) expr f =
  hatchpath (sep, f transformed ztr);
enddef;

% gt - common combinations.

save drawshade;

def drawshade (expr sp) expr f =
  draw f transformed ztr;
  shade (sp) f;
enddef;

save drawhatch;

def drawhatch (expr sep) expr f =
  draw f transformed ztr;
  hatch (sep) f;
enddef;

% * rest of macros start in graphing 
% coordinates but convert to pixel 
% to draw
% * variables ending in "_px" 
% converted to pixel
% * exceptions are the TeX dimensions
% here called:
% ptwd, hlen, dlen, slen, len, sp, sep
% all of which are in pixel coordinates
% * macros beginning with "mk" operate
% entirely in graphing coordinates

% general path construction

save mkpath;

vardef mkpath(expr smooth, cyclic, n)
  (suffix pts) =
 if smooth:
  if cyclic:
   pts[1]{pts[2]-pts[n]}
  else:
   pts[1]
  fi
  for i=2 upto n-1:
   ..pts[i]{pts[i+1]-pts[i-1]}
  endfor
  if cyclic:
   ..pts[n]{pts[1]-pts[n-1]}..cycle
  else:
   ..pts[n]
  fi
 else:
  for i=1 upto n-1:
    pts[i] --
  endfor
  pts[n]
  if cyclic:
   -- cycle
  fi
 fi
enddef;

% points, lines, and arrows

save pointd, p;

def pointd(expr a,ptwd) = 
 pair p_px;
 p_px:=a transformed ztr;
 fill fullcircle scaled ptwd shifted p_px; 
enddef;

save line;

def line(expr a,b) = 
 draw (a..b) transformed ztr; 
enddef;

% gt - arrowpath draws path f
% with an arrowhead;
% hlen is in pixel coordinates;
% f is in graphics coords;
% f is transformed transparently.
% Compare shade, hatch, etc.,
% and contrast shadepath.

save arrowpath, f_px;

def arrowpath (expr hlen) expr f =
  path f_px;
  f_px := f transformed ztr;
  draw f_px;
  headpath (f_px, hlen);
enddef;

% gt - arrow now uses arrowpath.

save arrow;

def arrow(expr tl,hd,hlen) =
 arrowpath (hlen) tl..hd ;
enddef;

% gt - "px" was too frequent
% in dottedline, and made the code
% hard to read, so I've deleted it.
% Only a and b are in graphics coords.

save dottedline,
  p, v, l, delta, n;

def dottedline (expr a, b, dlen, slen) =
  pair p[];
  p1 := a transformed ztr;
  p3 := b transformed ztr;
  l := length (p3-p1); 
  if (l > 2dlen) and 
    (dlen >= 0) and (slen >= 0): 
  else: 
    pair v;
    v := unitvector (p3-p1);
    n := floor ((l+slen-dlen) / (dlen+slen));
    delta := (l-dlen) / n - (dlen+slen);
    for i=1 upto n:
      p2 := p1 + dlen * v; 
      draw p1..p2; 
      p1 := p2 + (slen+delta) * v;
    endfor; 
  fi;
  draw p1..p3;
enddef;

save dottedarrow;

def dottedarrow(expr tl,hd,dlen,
  slen,hlen) =
 dottedline(tl,hd,dlen,slen); 
 headpath((tl..hd) transformed ztr,hlen);
enddef;

% axes and axis marks

save axes;

def axes(expr hlen) =
 arrow((0,yneg),(0,ypos),hlen); 
 arrow((xneg,0),(xpos,0),hlen);
enddef;

save xmarks;

def xmarks(expr len)(text t) =
 for a=t: 
  draw (xconv(a),yconv(0)-(len/2))..
    (xconv(a),yconv(0)+(len/2)); 
 endfor; 
enddef;

save ymarks;

def ymarks(expr len)(text t) =
 for a=t: 
  draw (xconv(0)-(len/2),yconv(a))..
    (xconv(0)+(len/2),yconv(a)); 
 endfor; 
enddef;

% upright rectangles

save mkrect;

vardef mkrect(expr ll,ur) =
 ll--(xpart ll,ypart ur)--
   ur--(xpart ur,ypart ll)--cycle
enddef;

save rect;

def rect(expr ll,ur) =
 draw mkrect(ll,ur) transformed ztr;
enddef;

save dottedrect;

def dottedrect(expr ll,ur,dlen,slen) =
 dottedline(ll,(xpart ll,ypart ur),
   dlen,slen);
 dottedline((xpart ll,ypart ur),ur,
   dlen,slen);
 dottedline(ur,(xpart ur,ypart ll),
   dlen,slen);
 dottedline((xpart ur,ypart ll),ll,
   dlen,slen);
enddef;

save block;

def block(expr ll,ur) =
 fill mkrect(ll,ur) transformed ztr;
enddef;

% gt - rectshade now uses shade.

save rectshade;

def rectshade(expr sp,ll,ur) =
  shade (sp) mkrect (ll, ur);
enddef;

% circles and ellipses

save mkellipse;

vardef mkellipse(expr center,radx,rady,
  angle) =
 save t;
 transform t; 
 t := identity
   xscaled (2 * radx)
   yscaled (2 * rady)
   rotated angle 
   shifted center;
 fullcircle transformed t
enddef;

save ellipse;

def ellipse(expr center,radx,rady,
  angle) =
 draw 
   mkellipse(center,radx,rady,angle)
   transformed ztr;
enddef;

save circle;

def circle(expr center,rad) =
 ellipse(center,rad,rad,0);
enddef;

% gt - ellshade now uses shade.

save ellshade;

def ellshade (expr sp, center, 
  radx, rady, angle) =
 shade (sp)
   mkellipse (center, radx, rady, angle);
enddef;

save circshade;

def circshade(expr sp, center,rad) =
 ellshade(sp,center,rad,rad,0);
enddef;

% circular arcs

% gt - mkarc now calculates
% n using ceiling, not floor;
% and saves theta, not i.

save mkarc;

vardef mkarc(expr center,from,sweep) =
  if sweep=0:
    from
  else:
   begingroup
    save n, theta, p;
    n := 1 + ceiling (abs (sweep) / 45);
    if n<3: n:=3; fi;
    theta:=sweep/(n-1);
    pair p[];
    p1:=from; 
    for i=2 upto n:
     p[i]:=p[i-1] 
       rotatedabout (center,theta);
    endfor;
    mkpath(true,false,n,p)
   endgroup
  fi
enddef;

% gt - note that when sweep is a multiple
% of 360 degrees, disp is logically
% infinite, not zero; then the center is
% at infinity.  In practice, arccenter
% ought not to be called in that case.

save arccenter;

vardef arccenter(expr from,to,sweep) =
 save midpt, disp;
 pair midpt;
 midpt:=(0.5)[from,to];
 disp:=
   if ((sweep mod 360)=0):
    0
   else:
    cosd(sweep/2)/sind(sweep/2)
   fi;
 midpt+(disp*((to-from) rotated 90)/2)
enddef;

% gt - mkarcto makes an arc given two points
% on the arc and the sweep angle.
% If sweep is a multiple of 360 degrees,
% then the arc is a straight line;
% if sweep is also nonzero, then that
% line should be infinite, but I use
% from--to instead.

save mkarcto;

vardef mkarcto (expr from, to, sweep) =
  if from = to:
    from
  elseif (sweep mod 360) = 0:
    from--to
  else
   begingroup
    save center;
    pair center;
    center:=arccenter (from, to, sweep);
    mkarc (center, from, sweep)
   endgroup
  fi
enddef;

% gt - arc now uses mkarcto.

save arc;

def arc(expr from,to,sweep) =
 draw mkarcto (from, to, sweep)
    transformed ztr;
enddef;

% gt - arcarrow now uses mkarcto
%      and arrowpath.

save arcarrow;

def arcarrow(expr hlen,from,to,sweep) =
  arrowpath (hlen) mkarcto (from, to, sweep);
enddef;

% gt - mkchordto makes a cyclic path from
% the arc from "from" to "to" with a sweep
% angle of "sweep", and its chord from
% "to" to "from".

save mkchordto;

vardef mkchordto (expr from, to, sweep) =
  mkarcto (from, to, sweep) -- cycle
enddef;

% gt - arcshade now uses mkchordto
%      and shade.

save arcshade;

def arcshade(expr sp,from,to,sweep) =
  shade (sp) mkchordto (from, to, sweep);
enddef;

% gt - three-point arcs.

save mkarcthree;

vardef mkarcthree (expr first, mid, last) =
  save p, sweep, n, theta, center;
  pair p[];
  p1 := first;
  sweep := 2 (angle (last-mid) - angle (mid-first));
  if abs (sweep) <= 90:
    n := 3;
    p2 := mid;
    p3 := last;
  else:
    n := 1 + ceiling (abs (sweep) / 45);
    theta := sweep / (n-1);
    pair center;
    center := arcthreecenter (first, mid, last);
    for i=2 upto n:
      p[i] := p[i-1] rotatedabout (center, theta);
    endfor;
  fi;
  mkpath (true, false, n, p)
enddef;

save arcthreecenter;

vardef arcthreecenter (expr first, mid, last) =
  save c, m, d;
  pair c, m[], d[];
  d1 := (mid - first) rotated 90;
  d2 := (last - mid) rotated 90;
  m1 := 0.5 [first, mid];
  m2 := 0.5 [mid, last];
  c = m1 + whatever * d1 = m2 + whatever * d2;
  c
enddef;

save arcthree;

def arcthree (expr first, mid, last) =
  draw mkarcthree (first, mid, last) transformed ztr;
enddef;

save arcthreearrow;

def arcthreearrow (expr hlen, first, mid, last) =
  arrowpath (hlen) mkarcthree (first, mid, last);
enddef;

% modified polar coordinates

% gt - mklinedir makes a path from point "a"
% to a point displaced "len" in direction "theta"
% from "a".

save mklinedir;

vardef mklinedir (expr a, theta, len) =
  a -- (a + len * (dir theta))
enddef;

% gt - linedir now uses mklinedir.

save linedir;

def linedir(expr a,theta,len) =
  draw mklinedir (a, theta, len)
	 transformed ztr;
enddef;

% gt - arrowdir now uses mklinedir
%      and arrowpath.

save arrowdir;

def arrowdir(expr hlen,a,theta,len) =
 arrowpath (hlen)
     mklinedir (a, theta, len);
enddef;

% gt - mkarcth makes an arc path with
% given center, radius "rad", initial
% angle "frtheta", and final angle
% "totheta".

save mkarcth;

vardef mkarcth (expr center,
    frtheta, totheta, rad) =
  save from;
  pair from;
  from := center + rad * (dir frtheta);
  mkarc (center, from, totheta-frtheta)
enddef;

% gt - arcth now uses mkarcth.

save arcth;

def arcth(expr center,
  frtheta,totheta,rad) =
  draw mkarcth (center, frtheta,
		totheta, rad)
    transformed ztr;
enddef;

% gt - arcth now uses mkarcth
%      and arrowpath.

save arctharrow;

def arctharrow(expr hlen,center, 
  frtheta,totheta,rad) =
 arrowpath (hlen)
     mkarcth (center, frtheta,
              totheta, rad);
enddef;

% gt - mkwedge makes a wedge-shaped path
% with apex at "center", radius "rad",
% initial angle "frtheta", and final angle
% "totheta".

save mkwedge;

vardef mkwedge (expr center, frtheta, totheta, rad) =
  center -- mkarcth (from, frtheta, totheta, rad)
         -- cycle
enddef;

% gt - wedge draws a sector of a circle.

save wedge;

def wedge (expr center, frtheta, totheta, rad) =
  draw mkwedge (center, frtheta, totheta, rad)
    transformed ztr;
enddef;

% gt - wedgeshade now uses mkwedge and shade.

save wedgeshade;

def wedgeshade (expr sp, center, 
  frtheta, totheta, rad) =
  shade (sp) mkwedge (center, frtheta, totheta, rad);
enddef;

% gt - drawshadewedge draws and shades a wedge.

save drawshadewedge;

def drawshadewedge (expr sp, center, 
  frtheta, totheta, rad) =
 draw mkwedge (center, frtheta, totheta, rad)
   transformed ztr;
 shade (sp) mkwedge (center, frtheta, totheta, rad);
enddef;

% curves

% gt - watch out for that "text containing a local
% variable's name" conflict!  I dearly wish that
% weren't a danger.
% Perhaps it's not so likely at the level of "mkcurve",
% as the "mk" macros are often fed numeric constants.

save mkcurve;

vardef mkcurve(expr smooth,cyclic)
  (text t) =
 save n_, p_;
 pair p_[];
 textpairs (t) (p_, n_);
 mkpath(smooth,cyclic,n_,p_)
enddef;

save curve;

def curve(expr smooth,cyclic)
  (text t) =
 draw mkcurve(smooth,cyclic,t)
   transformed ztr; 
enddef;

% gt - curvedarrow now uses arrowpath.

save curvedarrow;

def curvedarrow(expr smooth,hlen)
  (text t) =
  arrowpath (hlen)
      mkcurve (smooth, false, t);
enddef;

% shading of cyclic curves

% gt - cycleshade now uses shade.

save cycleshade;

def cycleshade(expr sp,smooth)(text t) =
  shade (sp) mkcurve (smooth,true,t);
enddef;

% gt - interpolated splines with controls.

% gt - mkipath uses the interpolation points,
% p[], and the left and right control points,
% l[] and r[].
% Observe that for cyclic I-splines, l[n] is
% used, not l1, though they are equal; this
% simplifies the algorithm.

save mkipath;

vardef mkipath (expr closed, n)
  (suffix p, l, r) =
  for i=1 upto n-1:
    p[i]..controls r[i] and l[i+1]..
  endfor
  if closed:
    cycle
  else:
    p[n]
  fi
enddef;

% gt - mkisplineA uses the I-spline data,
% in the order that Fig 2.1 gives them,
% points line pl and control line cl,
% stores them in p[], l[] and r[],
% then calls mkipath.
%
% pl should have the form:
%   (x1,y1) ... (xn,yn)
% and cl the form:
%   (lx1,ly1) (rx1,ry1) ... (lxn,lyn) (rxn,ryn)
% which reflect how Fig outputs its data.
%
% Don't feed it the "9999 9999", please!
%
% Perhaps the input should be massaged by a
% preprocessor program (e.g. in C), to separate
% the initially interleaved left and right control
% points, before being given to graphbase?
% That would simplify mkisplineA, and run faster.

save mkisplineA;

vardef mkisplineA (expr closed)
  (text pl) (text cl) =
  save p, l, r, n, i, isleft;
  pair p[], l[], r[];
  boolean isleft;
  textpairs (pl) (p, n);
  i := 1;
  isleft := true;
  for b=cl:
    if isleft:
      l[i] := b;
      isleft := false;
    else:
      r[i] := b;
      i := i+1;
      isleft := true;
    fi;
  endfor;
  mkipath (closed, n, p, l, r)
enddef;

% gt - mkisplineB uses the points line,
% and the separated left and right controls.
%
% See how much simpler this is than
% mkisplineA.

save mkisplineB;

vardef mkisplineB (expr closed)
  (text pl) (text lc) (text rc) =
  save p, l, r, n, i;
  pair p[], l[], r[];
  textpairs (pl) (p, n);
  textpairs (lc) (l, i);
  textpairs (rc) (r, i);
  mkipath (closed, n, p, l, r)
enddef;

% gt - the usual variations.
%
% These use mkisplineA.  I'd prefer
% mkisplineB.

% draw an interpolated spline,
% with points line pl and interleaved
% control line cl.

save ispline;

def ispline (expr closed)
  (text pl) (text cl) =
  draw mkisplineA (closed) (pl) (cl)
    transformed ztr;
enddef;

save isplinearrow;

def isplinearrow (expr hlen, closed)
  (text pl) (text cl) =
  arrowpath (hlen)
      mkisplineA (closed) (pl) (cl);
enddef;

% gt - isplineshade assumes that the
% I-spline is closed.

save isplineshade;

def isplineshade (expr sp)
  (text pl) (text cl) =
  shade (sp)
      mkisplineA (true) (pl) (cl);
enddef;

% functions

% gt - better be on the safe side with
% the function text, so use "_" on local
% variables in "mkfcn".

save mkfcn;

vardef mkfcn(expr smooth,bmin,bmax,bst)
  (suffix bv)(text fcnpr) =
 save p_, i_;
 pair p_[];
 i_ := 0;
 for bv=bmin step bst 
   until bmax+(bst/2):
  p_[incr i_] := fcnpr; 
 endfor;
 mkpath (smooth, false, i_ , p_)
enddef;

save function;

def function(expr smooth,xmin,xmax,st)
  (text fx) =
 draw mkfcn (smooth, xmin, xmax, st,
   x, (x,fx))
   transformed ztr; 
enddef;

save parafcn;

def parafcn(expr smooth,tmin,tmax,st)
  (text ft) =
 draw mkfcn (smooth, tmin, tmax, st,
   t, ft)
   transformed ztr; 
enddef;

% gt - mksfn constructs a path from
% two functions and the verticals
% at either side.
%
% mksfn is used by shadefcn.

save mksfn;

vardef mksfn (expr smooth, xmin, xmax, st)
    (text fcni) (text fcnii) =
  mkfcn(smooth,xmin,xmax,st,x,(x,fcni))
  --
  reverse
  mkfcn(smooth,xmin,xmax,st,x,(x,fcnii))
  -- cycle
enddef;

% gt - description:
% shadefcn shades between two functions over
% the range xmin to xmax, stepping by st,
% with dot spacing sp.
% it does not draw the functions.

% gt - shadefcn now uses mksfn.
% I don't see the connection between the dot
% spacing sp and the function step size st.

save shadefcn, st;

def shadefcn(expr sp, xmin, xmax)
    (text fcni)(text fcnii) =
  st := unxconv (sp);
  shade (sp)
    mksfn (false, xmin, xmax, st) (fcni) (fcnii);
enddef;

% gt - drawshadefcn draws both functions fcni
% and fcnii, and shades between them.

save drawshadefcn;

def drawshadefcn (expr sp, smooth, xmin, xmax, st)
    (text fcni) (text fcnii) =
  function (smooth, xmin, xmax, st) (fcni);
  function (smooth, xmin, xmax, st) (fcnii);
  shadefcn (sp, xmin, xmax) (fcni) (fcnii);
enddef;

enddef;  % mfpicenv

def endmfpicenv =
 endgroup;
enddef;