%%
%% This is file `grafbase.mf',
%% generated with the docstrip utility.
%%
%% The original source files were:
%%
%% grafbase.dtx  (with options: `MF')
%% 
%% -------------------------------------------------------------------
%% 
%% Copyright 2002--2003, Daniel H. Luecking
%% 
%% Mfpic consists of the 3 files mfpic.dtx, grafbase.dtx and mfpic.ins
%% and the 5 files they generate: mfpic.tex, mfpic.sty, grafbase.mf,
%% grafbase.mp, and dvipsnam.mp.
%% 
%% Documentation, examples, and ancillary files are separate. See
%% readme.1st for a complete list.
%% 
%% Mfpic may be distributed and/or modified under the conditions of the
%% LaTeX Project Public License, either version 1.2 of this license or (at
%% your option) any later version. The latest version of this license is in
%% 
%%    http://www.latex-project.org/lppl.txt
%% 
%% and version 1.2 or later is part of all distributions of LaTeX version
%% 1999/12/01 or later.
%% 
%% With respect to the proposed draft LPPL-1.3: mfpic has maintenance
%% status "maintained". The Current Maintainer is Daniel H. Luecking. There
%% are several Base Interpreters: TeX, LaTeX, MetaPost and Metafont.
%% 
string fileversion, filedate;
fileversion := "0.7a beta"; filedate := "2004/04/16";
def GBmsg expr s = message "Grafbase: " & s; enddef;
def GBerrmsg (expr s) expr t =
  errhelp t;
  errmessage "Grafbase: "& s;
  errhelp "";
enddef;
if (known grafbaseversion) or (known grafbase) :
  GBmsg "You have loaded grafbase more than once! " &
    "Please make sure that it is loaded only once.";
  endinput;
fi

boolean grafbase; grafbase := true;
boolean MFPIC; MFPIC := false;

def checkversions (expr g)=
  numeric grafbaseversion; grafbaseversion := g;
  if unknown mfpicversion :  % no mfpic, or  < 0.63
    GBmsg "Recent mfpic not detected.";
  elseif g = mfpicversion :
    MFPIC := true;
  else:
    GBerrmsg ("version mismatch")
      "The installation may be broken: mfpic and grafbase " &
      "versions do not match.";
  fi
enddef;

checkversions (70);

if unknown base_name :         input plain;
elseif not string base_name :  input plain;
elseif base_name <> "plain" :  input plain;
fi

boolean METAPOST;

if known color Maurits Cornelis Escher :  METAPOST := true;
else:  METAPOST := false;
fi

if METAPOST :
  GBerrmsg ("wrong compiler")
    "You may have input to Metapost a file designed for Metafont. " &
    "Instead of the file grafbase.mf, Metapost should be using " &
    "grafbase.mp. Make sure Metapost can find it.";
fi

if (unknown debug) or (not boolean debug) :
  boolean debug; debug := false;
fi

def GBdebug =
  begingroup
    save >>; def >> =  message  enddef;
    >> "Grafbase DEBUG: ";
enddef;
def GBenddebug =
    >> "End DEBUG";
  endgroup
enddef;

vardef mftitle expr t =
  t; message t;
enddef;

if unknown mode :
  GBerrmsg ("unknown Metafont mode")
    "Please use \mode:=localfont; or a mode known on your " &
    "system. If you continue, ljfour mode will be tried.";
  mode := lfjour;
fi

mode_setup;
if debug :
  GBdebug;
    >> "pixels_per_inch = " & decimal pixels_per_inch;
  GBenddebug;
fi
font_identifier := "MFpic graphics";
font_coding_scheme := "Arbitrary";
interim designsize := 128pt#;

if unknown aspect_ratio: aspect_ratio := 1; fi
if unknown hppp : hppp := 1 fi;
if unknown currenttransform :
  transform currenttransform;
  currenttransform := identity yscaled aspect_ratio;
fi

interim turningcheck := 0;
interim warningcheck := 0;
numeric unitlen, xscale, yscale, xneg, xpos, yneg, ypos;

unitlen := 1 bp#;
xscale := 7.2;  % (xscale * unitlen) = 1/10 inch
yscale := 7.2;  % (yscale * unitlen) = 1/10 inch
xneg := 0; xpos := 10;
yneg := 0; ypos := 10;

newinternal radian, pi, deg;
deg := 1; pi := 3.14159;
radian := 180/pi;


newinternal penwd; penwd := 0.5pt;
pen drawpen;

def resizedrawpen (expr s) =
  interim penwd := s;
  setvariable (pen) (drawpen) (pencircle scaled penwd);
  save currentpen; pen currentpen; pickup drawpen;
enddef;

numeric hatchwd; hatchwd := 0.5bp;
pen hatchpen; hatchpen := pencircle scaled hatchwd;

boolean clipall; clipall := false;
boolean ClipOn; ClipOn := false;
path ClipPath[]; numeric ClipPath; ClipPath = 0;
boolean truebbox; truebbox := false;

def DoClip (suffix v) =
  if ClipOn and (ClipPath > 0) : clipsto (v, ClipPath); fi
enddef;

def noclip (text t) =
  hide( save ClipOn; boolean ClipOn; ClipOn := false; t)
enddef;

boolean showbbox; showbbox := false;

let color = numeric; color black, white;
black := 0; white := 1;
def _wc_ = killtext enddef;

color currentcolor, fillcolor, drawcolor, hatchcolor,
  headcolor, pointcolor, tlabelcolor, background;
currentcolor := fillcolor := drawcolor := hatchcolor :=
  headcolor := pointcolor := tlabelcolor := black;
background := white;

vardef snapto expr t =
  if unknown t           :  0
  elseif not (numeric t) :  0
  elseif t < 0           :  0
  elseif t > 1           :  1
  else                   :  t
  fi
enddef;

vardef gray (expr g) = (snapto g)*white enddef;

vardef makeclr (expr r, g, b) =
  gray (sqrt((2r*r + 4g*g + b*b)/7))
enddef;

vardef rgb (expr r, g, b) =
  makeclr (snapto r, snapto g, snapto b)
enddef;

vardef cmyk (expr c, m, y, k) =
  rgb (1-c-k, 1-m-k, 1-y-k)
enddef;

vardef RGB (expr R, G, B) =
  rgb (R/255, G/255, B/255)
enddef;

vardef named (suffix c) = forceclr (c) enddef;
vardef forceclr (expr c) =
  if unknown c     :  black
  elseif numeric c :  gray (c)
  elseif color c   :  c
  else             :  black
  fi
enddef;

color red, green, blue, cyan, magenta, yellow;
red     := rgb (1, 0, 0);
green   := rgb (0, 1, 0);
blue    := rgb (0, 0, 1);
cyan    := rgb (0, 1, 1);
magenta := rgb (1, 0, 1);
yellow  := rgb (1, 1, 0);

vardef list (suffix v) (text lst) =
  v := 0;  for _itm = lst: v[incr v] := _itm; endfor
enddef;

def map (text proc) (text lst) =
  hide(_map := 0;)
  for _a = lst :
    if _map = 0 :  hide(_map := 1;)  else:  ,  fi
    proc(_a)
  endfor
enddef;

vardef knownarray suffix arr =
  save _kna; boolean _kna;
  _kna := (known arr) and (numeric arr);
  if _kna :
    _kna := (arr = floor arr) and (arr >= 1);
    for _idx = 1 upto arr :
      exitif not _kna;
      _kna := known arr[_idx];
    endfor
  fi
  _kna
enddef;

def copyarray(suffix from, to) =
  to := 0;
  for _idx = 1 upto from:
    to[incr to] := from[_idx];
  endfor
enddef;

def maparr (text proc) (suffix p) =
  for _idx = 1 upto p: proc (p[_idx]); endfor
enddef;

def textpairs (suffix p) (text t) =
  numeric p; pair p[]; list (p) (t);
enddef;

vardef chpair (text proc) (expr p) =
  (proc (xpart p), proc (ypart p))
enddef;

vardef floorpair (expr p) = chpair (floor) (p) enddef;
vardef ceilingpair (expr p) = chpair (ceiling) (p) enddef;
vardef hroundpair (expr p) = chpair (hround) (p) enddef;

vardef emin (expr a, b) = if a < b :  a  else:  b  fi enddef;
vardef emax (expr a, b) = if a > b :  a  else:  b  fi enddef;

vardef pairmin (expr z, w) =
    ( emin (xpart z,  xpart w), emin (ypart z, ypart w ) )
enddef;
vardef pairmax (expr z, w) =
    ( emax (xpart z,  xpart w), emax (ypart z, ypart w ) )
enddef;

vardef minpair (suffix p) =
  save _mp; pair _mp; _mp := p1;
  for _idx = 2 upto p - 1 :
    _mp := pairmin(_mp, p[_idx]);
  endfor
  pairmin (_mp, p[p])
enddef;
vardef maxpair (suffix p) =
  save _mp; pair _mp; _mp := p1;
  for _idx = 2 upto p - 1: _mp := pairmax(_mp, p[_idx]); endfor
  pairmax (_mp, p[p])
enddef;

transform ztr, vtr;
def setztr =
  if debug :
    GBdebug;
      >> "charwd = " & decimal charwd & "pt#";
      >> "charht = " & decimal charht & "pt#";
      >> "w_ = " & decimal w_ & " pixels";
      >> "h_ = " & decimal h_ & " pixels";
      >> "unitlen = " & decimal unitlen & "pt#";
      >> "hppp = " & decimal hppp;
      >> "xneg = " & decimal xneg;
      >> "xpos = " & decimal xpos;
      >> "yneg = " & decimal yneg;
      >> "ypos = " & decimal ypos;
      >> "xscale = " & decimal xscale;
      >> "yscale = " & decimal yscale;
    GBenddebug;
  fi
  save ztr, vtr;
  transform ztr, vtr;
  vtr := identity xscaled (xscale) yscaled (yscale)
    scaled (unitlen*hppp);
  ztr := identity shifted (-(xneg, yneg)) transformed vtr;
  if debug :
    GBdebug;
      >> "ztr: ";
      show ztr;
      >> "vtr: ";
      show vtr;
    GBenddebug;
  fi
enddef;

vardef zconv (expr a) = a transformed ztr enddef;
vardef invzconv (expr a) = a transformed (inverse ztr) enddef;
vardef vconv (expr v) = v transformed vtr enddef;
vardef invvconv (expr v) = v transformed (inverse vtr) enddef;

def active_plane = currentpicture enddef;

vardef image(text t) =
  save currentpicture; picture currentpicture;
  currentpicture := nullpicture;
  t;
  currentpicture
enddef;

def initpic =
  setztr;
  resizedrawpen (penwd);
  if ClipOn : ClipPath := 1;
    ClipPath1 := rect (origin, (w_, h_));
  fi
  if debug :
    GBdebug;
      >> "Drawing nominal bounding box around picture";
    GBenddebug;
    noclip ( safedraw rect (origin, (w_, h_)) );
  fi
enddef;

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

def setvariable (text kind) (suffix name) (expr value) =
  save name; kind name; name := value;
enddef;

string extra_beginmfpic; extra_beginmfpic := "";
string   extra_endmfpic;   extra_endmfpic := "";

def beginmfpic (expr ch) =
  begingroup
    gcode := ch;
    save w_, h_, d_; numeric w_, h_, d_;
    charwd := (xpos-xneg)*xscale*unitlen;
    charht := (ypos-yneg)*yscale*unitlen;
    chardp := 0;
    charcode := if known ch :  byte ch  else:  0  fi;
    w_ := hround(charwd*hppp);
    h_ := vround(charht*hppp);
    d_ := vround(chardp*hppp);
    charic := 0; clearxy; clearit; clearpen;
    scantokens extra_beginchar;
    initpic;
    scantokens extra_beginmfpic;
enddef;

def endmfpic =
  scantokens extra_endmfpic;
  if debug :
    GBdebug;
      >> "TFM charwd = " & decimal charwd & "pt#";
      >> "TFM charht = " & decimal charht & "pt#";
    GBenddebug;
  fi
  if proofing>0 : makebox(proofrule); fi
  chardx := w_;     % desired width of character in pixels
  DoClip (active_plane);
  if clipall : clipto (active_plane) rect(origin, (w_, h_)); fi
  if showbbox : noclip ( safedraw rect (origin, (w_, h_)) ); fi
  shipit;
  if displaying > 0 : makebox(screenrule); showit; fi
  endgroup
enddef;

pair label_adjust; label_adjust := (0, 0);
numeric label_sep;    label_sep :=  0;

def verbatimtex = killtext enddef;

def newgblabel (expr hf, vf, BL, r) (text s) (text pts) = enddef;
vardef ref_shift (expr hf, vf, BL, ll, ur) =
  - ( (hf)[xpart ll, xpart ur],
      (vf)[if BL: 0 else: (ypart ll) fi, ypart ur] )
enddef;

vardef thegblabel (expr z, r, p) =
  ((p shifted z) rotated r) shifted label_adjust
enddef;

vardef textrect (expr lbl, rad, loc) =
  textrectx (.5, .5, false, 0) (origin, lbl, rad, loc)
enddef;
vardef textoval (expr lbl, mult, loc) =
  xellipse (true, .5, .5, false, 0) (origin, lbl, mult, loc)
enddef;
vardef textellipse (expr lbl, rat, loc) =
  xellipse (false, .5, .5, false, 0) (origin, lbl, rat, loc)
enddef;

boolean roundends; roundends := true;
vardef textrectx (expr a, b, c, rot, xy, lbl, rad, loc) =
  save ll, ur, _r, f, zz;
  pair ll, ur, zz;  path f;
  labeldims (xy, lbl) (ll, ur);
  _r :=
    if boolean rad :
      if rad :  emin (xpart (ur-ll), ypart (ur-ll))/sqrt(2)
      else:     0
      fi
    elseif numeric rad :  rad
    else:                 0
    fi;
  if _r = 0 :
    f :=  rect(ll, ur);
  else:
    save p, q;
    pair p[];  path q;
    p1 := ur - _r*dir(45);    % center of upper right arc
    p3 := ll + _r*dir(45);    % lower left
    p2 := (xpart p3, ypart p1); % upper left
    p4 := (xpart p1, ypart p3); % lower right
    q  := quartercircle scaled 2_r;
    if _r > 0:
      f  := (q shifted p1) -- (q rotated 90 shifted p2)
            -- (q rotated 180 shifted p3)
            -- (q rotated -90 shifted p4)
            -- cycle;
    else:
      f  := (q shifted p1) -- (q rotated -90 shifted p4)
            -- (q rotated 180 shifted p3)
            -- (q rotated 90 shifted p2)
            -- cycle;
    fi
  fi
  invvconv(thegblabel(ref_shift(a, b, c, ll, ur), rot, f)) shifted loc
enddef;

def textovalx = xellipse (true) enddef;
def textellipsex = xellipse (false) enddef;

vardef xellipse (expr aspect, a, b, c, r, xy, lbl, mult, loc) =
  if mult = 0 :
    textrectx (a, b, c, r) (xy, lbl, 0, loc)
  else:
    save ll, ur, cc, ww, hh, f;
    pair ll, ur, cc;  path f;
    labeldims (xy, lbl) (ll, ur);
    cc := .5[ll, ur]; % center
    (ww, hh) = ur - cc;
    if (ww = 0) or (hh = 0) : % make a line:
      f = (ll--ur);
    else:
      save aa, bb, mm;
      mm := if aspect : ww/hh*mult else: mult fi;
      aa := ww ++ hh*mm;
      bb := aa/mm;
      f := ellipse(cc, aa, bb, 0);
    fi
    invvconv(thegblabel (ref_shift (a, b, c, ll, ur), r, f)) shifted loc
  fi
enddef;


def labeldims (expr xy, lbl) (suffix ll, ur) =
  if pair lbl :
    ll := xy; ur := lbl;
  else:
    ll := ur := (0, 0);
  fi
  ll := ll - label_sep*(1, 1);
  ur := ur + label_sep*(1, 1);
enddef;

newinternal nottoosmall; nottoosmall := eps/2 + 2epsilon;
newinternal reallysmall; reallysmall := 3epsilon;
def signof (expr X) = if X < 0 : - fi enddef;
def TruncateWarn expr s =  GBmsg s & " too large; truncating";  enddef;

vardef secd primary X =
  save temp; temp := cosd(X);
  if abs(temp) < reallysmall :
    TruncateWarn "Secant";
    temp := signof (temp) reallysmall;
  fi
  1/temp
enddef;
vardef tand primary X =  sind(X)*secd(X)  enddef;

vardef cscd primary X =
  save temp; temp := sind(X);
  if abs(temp) < reallysmall :
    TruncateWarn "Cosecant";
    temp := signof(temp) reallysmall;
  fi
  1/temp
enddef;
vardef cotd primary X =  cosd(X)*cscd(X)  enddef;

vardef acos primary X =  angle (X, 1 +-+ X)  enddef;
vardef asin primary X =  angle (1 +-+ X, X)  enddef;
vardef atan primary X =  angle (1, X)  enddef;

vardef sin primary X =  sind (X*radian)  enddef;
vardef cos primary X =  cosd (X*radian)  enddef;
vardef tan primary X =  tand (X*radian)  enddef;
vardef cot primary X =  cotd (X*radian)  enddef;
vardef sec primary X =  secd (X*radian)  enddef;
vardef csc primary X =  cscd (X*radian)  enddef;

vardef invcos primary X =  (acos X)/radian  enddef;
vardef invsin primary X =  (asin X)/radian  enddef;
vardef invtan primary X =  (atan X)/radian  enddef;

vardef exp primary X = mexp (256 * X) enddef;
vardef ln  primary X = (mlog X) / 256 enddef;
def log = ln enddef;
vardef logbase (expr B) primary X = (mlog X)/(mlog B) enddef;
def logtwo = logbase( 2) enddef;
def logten = logbase(10) enddef;

vardef Arg primary Z = (angle Z)/radian enddef;
vardef Log primary Z = (ln(abs(Z)), Arg (Z)) enddef;
vardef cis primary T  = dir(radian*T) enddef;
vardef zexp primary Z = (exp (xpart Z))*(cis(ypart Z)) enddef;
vardef sgn primary Z =
  if Z = origin :  origin  else:  unitvector Z  fi
enddef;

vardef cosh primary X =
  save temp; temp := 2 exp (-abs(X));
  if temp < reallysmall :
    TruncateWarn "Cosh";
    temp := reallysmall;
  fi
  1/temp + temp/4
enddef;

vardef sinh primary X =
  save temp; temp := 2 exp (-abs(X));
  if temp < reallysmall :
    TruncateWarn "Sinh";
    temp := reallysmall;
  fi
  signof (X) (1/temp - temp/4)
enddef;

vardef sech primary X =
  save temp; temp := exp(-(abs (X)));
  2temp/(1 + temp*temp)
enddef;

vardef tanh primary X =
  save temp; temp := exp(-2(abs (X)));
  signof (X) (1 - temp)/(1 + temp)
enddef;

vardef csch primary X =
  save temp; temp := exp(-(abs (X)));
  if abs(1 - temp*temp) < reallysmall :
    TruncateWarn "Csch";
    signof (X) 2temp / reallysmall
  else:
    signof (X) 2temp / (1 - temp*temp)
  fi
enddef;

vardef coth primary X =
  save temp; temp := tanh(X);
  if abs(temp) < reallysmall :
    TruncateWarn "Coth";
    temp := signof (temp) reallysmall;
  fi
  1/temp
enddef;

vardef acosh primary y =
  if y < 1 :
    GBerrmsg ("Undefined function: acosh " & decimal y)
      "If you proceed, a value of 0 will be used. " &
      "Expect more errors later.";
    0
  else:
    ln (y + (y+-+1))
  fi
enddef;

vardef asinh primary y = ln (y + (y++1)) enddef;

vardef atanh primary y =
  if abs (y) < 1 :
    (ln(1+y) - ln(1-y))/2
  else:
    GBerrmsg ("Undefined function: atanh " & decimal y)
      "If you proceed, a value of plus or minus infinity " &
      "will be used. Expect more errors later.";
    signof (y) infinity
  fi
enddef;

vardef polar (expr p) = (xpart p) * dir (ypart p) enddef;
def id (expr x) = x enddef;

transform T_stack[]; T_stack := 0;
def T_push (expr T) = T_stack[incr T_stack] := T; enddef;
def T_pop (suffix $) =
  if T_stack > 0 :
    $ := T_stack[T_stack]; T_stack := T_stack - 1;
  fi
enddef;

def bcoords = hide ( T_push (ztr); T_push (vtr) ) enddef;
def ecoords = hide ( T_pop (vtr); T_pop (ztr) ) enddef;

def apply_t (text Transformer) =
  ztr := identity Transformer transformed ztr;
  vtr := ztr shifted - zconv(origin);
enddef;

def xslant = slanted enddef;  % (x+sy, y).
def yslant primary s =  % (x, y+sx).
  transformed
    begingroup
      save _T; transform _T;
      origin transformed _T = origin;
      (1, 0) transformed _T = (1, s);
      (0, 1) transformed _T = (0, 1);
      _T
    endgroup
enddef;

def zslant primary p =  % (xu+yv, xv+yu), where p = (u, v).
  transformed
    begingroup
      save _T; transform _T;
      xpart _T = ypart _T = 0;
      xxpart _T = yypart _T = xpart p;
      xypart _T = yxpart _T = ypart p;
      _T
    endgroup
enddef;

def xyswap =  zslant (0, 1) enddef;
def boost primary X = zslant (cosh X, sinh X) enddef;

vardef rotatedpath (expr p, th) expr f =
  f transformed vtr rotatedaround (p transformed vtr, th)
    transformed (inverse vtr)
enddef;
vardef scaledpath (expr p, s) expr f =
  f shifted -p scaled s shifted p
enddef;
vardef xslantedpath (expr b, s) expr f =
  f shifted (0, -b) slanted s shifted (0, b)
enddef;
def slantedpath = xslantedpath enddef;
vardef yslantedpath (expr a, s) expr f =
  f shifted (-a, 0) yslant s shifted (0, a)
enddef;
vardef xscaledpath (expr a, s) expr f =
  f shifted (-a, 0) xscaled s shifted (a, 0)
enddef;
vardef yscaledpath (expr b, s) expr f =
  f shifted (0, -b) yscaled s shifted (0, b)
enddef;
vardef shiftedpath (expr v) expr f = f shifted v enddef;
vardef reflectedpath (expr p, q) expr f =
  f transformed vtr
    reflectedabout (p transformed vtr, q transformed vtr)
    transformed (inverse vtr)
enddef;
vardef xyswappedpath expr f = f xyswap enddef;
vardef transformedpath (text Transformer) expr f =
    f Transformer
enddef;

vardef partialpath (expr a, b) expr f =
  save p; path p;
  p := zconv (f) scaled (1/unit_of_length);
  save cumlen, totlen, idx, ta, tb;
  totlen := makelengtharray(cumlen) p;
  idx := 0;
  if a <= b:
    ta := gettime (cumlen, idx) (a*totlen);
    tb := gettime (cumlen, idx) (b*totlen);
  else:
    tb := gettime (cumlen, idx) (b*totlen);
    ta := gettime (cumlen, idx) (a*totlen);
  fi
  subpath (ta, tb) of f
enddef;
vardef gsubpath (expr a, b) expr f = subpath (a, b) of f enddef;

def mono (suffix u) = cull u keeping (1, infinity); enddef;

def andto (suffix u) (expr v) =
  mono (u); addto u also v; cull u keeping (2, 2);
enddef;
primarydef u picand v =
  begingroup  save t; picture t;
    t := u; andto (t, v); t
  endgroup
enddef;

def orto (suffix u) (expr v) =
  mono (u); addto u also v; cull u keeping (1, 2);
enddef;
primarydef u picor v =
  begingroup  save t; picture t;
    t := u; orto (t, v); t
  endgroup
enddef;

def xorto (suffix u) (expr v) =
  mono (u); addto u also v; cull u keeping (1, 1);
enddef;
primarydef u picxor v =
  begingroup  save t; picture t;
    t := u; xorto (t, v); t
  endgroup
enddef;

def subto (suffix u) (expr v) =
  mono (u); addto u also -v; cull u keeping (1, 1);
enddef;
primarydef u picsub v =
  begingroup save t; picture t;
    t := u; mono (t); subto (t, v); t
  endgroup
enddef;

def coloraddto (expr clr) (suffix u) (expr v) =
  if clr < white :
    orto (u, v);
  else:
    subto (u) (v);
  fi;
enddef;

vardef interior expr c =
  save v; picture v; v := nullpicture;
  addto v contour (c.t_);
  cull v dropping (0, 0);
  v
enddef;
vardef interiors suffix cc =
  save _ints; picture _ints; _ints := nullpicture;
  for _idx = 1 upto cc:
    addto _ints also interior cc[_idx]);
  endfor
  mono (_ints);
  _ints
enddef;

def clipto (suffix vt) expr c =
  andto (vt, interior c);
enddef;
def clipsto (suffix vt, cc) =
  andto (vt, interiors cc);
enddef;

vardef Clipped (suffix vt) expr c =
  save _Cl; picture _Cl; _Cl := vt; clipto (_Cl) c; _Cl
enddef;
def clipped = Clipped enddef;
def clip = Clipped enddef;

vardef picneg (suffix vt) expr c =
  mono (vt);
  (interior c) picsub vt
enddef;

def shpath (suffix v) (expr q, f) =
  addto v doublepath (f.t_) withpen (q.t_);
enddef;

numeric minpenwd;
minpenwd := 1;      % 1 pixel
vardef picpath expr d =
  save v; picture v; v := nullpicture;
  if penwd >= minpenwd :
    shpath (v, drawpen) (d);
    mono (v);
  fi
  v
enddef;

def picdot (suffix v) (expr w, p) =
  addto v also
    (w shifted hroundpair (p.t_));
enddef;

vardef setdot (expr apath, sc) =
  if cycle apath : interior
  else           : picpath
  fi
    (apath scaled emax(ceiling (sc), minpenwd))
enddef;

numeric shadepicsize; shadepicsize := 0.8bp;
vardef shadepic (suffix dims) (expr grparam) =
  pair dims;
  save _frac; _frac := 2*emin(grparam, 1 - grparam);
  save _hp, _vp, _dotwd, _dotht;
  if aspect_ratio < 1 :
    _vp := emax (2, hround(shadepicsize.o_));
    _hp := hround (_vp._o_);
    _dotwd := hround (_hp*sqrt _frac);
    _dotht := if _dotwd = 0 : 0
              else: hround (_hp*_vp*_frac/_dotwd)
              fi;
  else:
    _hp := emax (2, hround (shadepicsize));
    _vp := hround (_hp.o_);
    _dotht := hround (_vp*sqrt _frac);
    _dotwd := if _dotht = 0 : 0
              else: hround (_hp*_vp*_frac/_dotht)
              fi;
  fi
  dims := ( _hp, _vp._o_ );
  save _shp; picture _shp; _shp := nullpicture;
  addto _shp contour rect ((0,0), (_dotwd, _dotht));
  picdot (_shp, _shp, dims);
  dims := 2dims; mono (_shp);
  if grparam >= .5 :  _shp
  else : (interior (rect ((0,0), dims))) picsub _shp
  fi
enddef;

vardef shaded (expr clr) expr c =
  if cycle c :
    if (clr <= black) or (clr >= white) :
      interior c
    else:
      save vsh, shpic, shdims, ll, ur;
      picture vsh, shpic; pair shdims, ll, ur;
      shpic := shadepic (shdims) (clr);
      gbbox (c, ll, ur);
      vsh := filledwith (shpic) (shdims, ll, ur);
      Clipped (vsh) c
    fi
  else: picpath c % should we? or just make it null?
  fi
enddef;

vardef filledwith (expr pic, dims, ll, ur) =
  save b, v; picture b, v;
  b := v := nullpicture;
  save fwdims, _ll, _ur; pair fwdims, _ll, _ur;
  fwdims := hroundpair (dims.t_);
  _ll   := floorpair  (ll.t_);
  _ur   := ur.t_;
  for s = xpart _ll step xpart fwdims until xpart _ur:
    addto b also pic shifted (s, 0);
  endfor
  for s = ypart _ll step ypart fwdims until ypart _ur:
    addto v also b shifted (0, s);
  endfor
  mono (v);
  v
enddef;

vardef thatchf (suffix v) (expr CT, sp, a, b) =
  save _sp;
  _sp = signof (ypart(b - a)) abs(sp);
  for _y = _sp*( ceiling ((ypart a)/_sp) ) step _sp until ypart b:
    shpath (v, hatchpen)
      ( ( (xpart a, _y)--(xpart b, _y) ) transformed CT );
  endfor
  mono(v);
enddef;

def tile (suffix atile) (expr unit, width, height, clipit) =
  picture atile.pic; atile.pic := nullpicture;
  numeric atile.wd, atile.ht;
  (atile.wd, atile.ht) = round ((width, height)*unit);
  boolean atile.clipon; atile.clipon := clipit;
  begingroup
    save active_plane;
    def active_plane = atile.pic enddef;
    save ztr, vtr; transform ztr, vtr;
    ztr := identity scaled unit; vtr := ztr;
    save ClipOn; boolean ClipOn;
    if clipit :
      ClipOn := true;
      save ClipPath; path ClipPath[];
      ClipPath = 1;
      ClipPath[1] = rect(origin, (atile.wd, atile.ht));
    else:
      ClipOn := false;
    fi
enddef;
def endtile =
    DoClip (active_plane);
  endgroup
enddef;

vardef is_tile (suffix atile) =
  (known atile.pic   ) and (picture atile.pic   ) and
  (known atile.wd    ) and (numeric atile.wd    ) and
  (known atile.ht    ) and (numeric atile.ht    ) and
  (known atile.clipon) and (boolean atile.clipon)
enddef;

vardef pnt@#  (expr p) = point       @# of p enddef;
vardef pre@#  (expr p) = precontrol  @# of p enddef;
vardef post@# (expr p) = postcontrol @# of p enddef;

vardef gbbox (expr g) (suffix ll, ur) =
  save _s; _s := emax(1, ceiling bbox_split);
  ur := ll := pnt 0 (g);
  if (length g) > 0 :
    for _j = 1 upto length g:
      ll := pairmin(ll, pnt[_j] (g)); ur := pairmax(ur, pnt[_j] (g));
    endfor
    for _j = 1 upto _s*(length g):
      ctrlsbbox (subpath ((_j-1)/_s, _j/_s) of g) (ll, ur);
    endfor
  fi
  if showbbox : noclip ( safedraw rect (ll, ur) ); fi
enddef;

numeric bbox_split; bbox_split := 2;
def ctrlsbbox (expr p) (suffix ll, ur) =
  ll := pairmin (ll, post0 (p)); ll := pairmin (ll, pre 1 (p));
  ur := pairmax (ur, post0 (p)); ur := pairmax (ur, pre 1 (p));
enddef;

def safedraw = colorsafedraw (drawcolor) enddef;
vardef colorsafedraw (expr clr) expr d =
  save v; picture v; v := picpath d;
  DoClip (v);
  coloraddto (clr) (active_plane, v);
enddef;

def NoCycleWarn expr s =
  GBmsg  s & " cannot be applied to an open path. "
           & "The path will be drawn instead.";
enddef;

def safefill = colorsafefill (fillcolor) enddef;
vardef colorsafefill (expr clr) expr c =
  if cycle c :
    save v; picture v; v := interior c;
    DoClip (v);
    if (clr > black) and (clr < white) :
      subto (active_plane) (v);
      v := nullpicture;
      v := shaded (clr) c;
    fi
    coloraddto (clr) (active_plane, v);
  else: NoCycleWarn "fill"; safedraw c;
  fi
enddef;
def safeunfill expr c =
  if cycle c : noclip (colorsafefill (background) c);
  else: NoCycleWarn "unfill"; safedraw c;
  fi
enddef;

def safeclip expr c =
  if cycle c : clipto (active_plane) c;
  else: NoCycleWarn "clip"; safedraw c;
  fi
enddef;

def store (suffix fs) expr f =
hide(
  if (not path f) and (not pair f) :
    GBerrmsg ("Second argument to `store' must be a path or pair")
      "";
  fi
  if not path fs : path fs; fi
  fs := f
)
enddef;
vardef stored (suffix fs) expr f = store (fs) f; f enddef;

def drawn = colordrawn (drawcolor) enddef;
vardef colordrawn (expr clr) expr f =
  colorsafedraw (clr) (zconv (f)); f
enddef;

def filled = colorfilled (fillcolor) enddef;
vardef colorfilled (expr clr) expr c =
  colorsafefill (clr) zconv (c); c
enddef;
vardef unfilled expr c = safeunfill zconv (c); c enddef;
vardef Clip expr c = safeclip zconv(c); c enddef;

numeric shadewd; shadewd := 0.5bp;
path shadedotpath;
shadedotpath := fullcircle; % unitsquare;

vardef shade (expr sp) expr f =
  save g; path g; g := zconv (f);
  save gr; numeric gr;
  gr :=  1 - (.88*abs(shadewd)/sp)**2;
  if not cycle g : NoCycleWarn "shade"; safedraw g;
  elseif gr <= 0 : safefill g;
  else:
    save ll, ur; pair ll, ur;
    gbbox (g, ll, ur);
    ll := floorpair (ll);
    save sh, v; picture sh, v;
    save dx; dx := ceiling(sp/(sqrt 2));
    sh := setdot (shadedotpath, abs(shadewd));
    v := filledwith (sh, 2(dx, dx), ll, ur);
    DoClip(v);
    orto (active_plane, Clipped (v) g);
    sh := sh shifted hroundpair ((dx, dx).t_);
    v := nullpicture;
    v := filledwith (sh, 2(dx, dx), ll, ur);
    DoClip(v);
    orto (active_plane, Clipped (v) g);
  fi
  f
enddef;

polkadotwd := 5bp;
mindotspace := 1bp;
path polkadotpath; polkadotpath := fullcircle;

vardef polkadot (expr sp) expr f =
  save g; path g; g := zconv (f);
  if not cycle g : NoCycleWarn "polkadot"; safedraw g;
  elseif sp <= emax (2*polkadotwd/3, mindotspace) :
    safefill g;
  else:
    save ll, ur; pair ll, ur;
    gbbox (g, ll, ur);
    save dx, dy; dx := sp/2; dy := dx*(sqrt 3);
    hshift := ((xpart (ur - ll)) mod dx)/2;
    vshift := ((ypart (ur - ll)) mod dy)/2;
    save p, dims; pair p, dims;
    p := ll + (hshift, vshift);
    dims :=  2(dx, dy);
    save v, thepolkadot; picture v, thepolkadot;
    thepolkadot := setdot (polkadotpath, polkadotwd);
    v := filledwith (thepolkadot, dims, p, ur);
    p := p + (dx, dy);
    orto (v, filledwith (thepolkadot, dims, p, ur));
    DoClip (v); clipto (v) g;
    if (fillcolor > black) and (fillcolor < white): % gray
      subto (active_plane) (v);
      thepolkadot :=
        shaded (fillcolor) polkadotpath scaled ceiling(polkadotwd);
      v := filledwith (thepolkadot, dims, p, ur);
      p := p - (dx, dy);
      orto (v, filledwith (thepolkadot, dims, p, ur));
      DoClip (v); clipto (v) g;
    fi
    coloraddto (fillcolor) (active_plane) (v);
  fi
  f
enddef;

def thatch = colorthatch (hatchcolor) enddef;
vardef colorthatch (expr clr) (expr sp, theta) expr f =
  save g; path g; g := zconv (f);
  if not cycle g : NoCycleWarn "hatch"; safedraw g;
  elseif sp <= abs(hatchwd) : colorsafefill (clr) g;
  else:
    save v; picture v; v := nullpicture;
    save CT; transform CT; CT := identity rotated theta;
    save ll, ur; pair ll, ur;
    gbbox (g transformed inverse CT, ll, ur);
    thatchf (v, CT, sp, ll, ur);
    DoClip(v);
    coloraddto (clr) (active_plane) (Clipped (v) g);
  fi
  f
enddef;

def hhatch (expr sp) = thatch (sp, 0)   enddef;
def vhatch (expr sp) = thatch (sp, 90)  enddef;
def lhatch (expr sp) = thatch (sp, -45) enddef;
def rhatch (expr sp) = thatch (sp, 45)  enddef;

def xhatch = colorxhatch (hatchcolor) enddef;
vardef colorxhatch (expr clr, sp) expr f =
  colorthatch (clr) (sp, 45) colorthatch (clr) (sp, -45) f
enddef;

vardef tess (suffix atile) expr c =
  save _g; path _g; _g := zconv (c);
  if not is_tile (atile) :
    GBerrmsg ("Tile parameter " & str atile & " of tess() is invalid")
      "This tile may be undefined or incorrectly defined. " &
      "If you proceed, tess() will be abandoned and the curve " &
      "merely drawn."; safedraw _g;
  elseif not cycle _g : NoCycleWarn "tess"; safedraw _g;
  else:
    save _ll, _ur; pair _ll, _ur;
    gbbox (_g, _ll, _ur);
    save _ts; picture _ts;
    _ts := filledwith (atile.pic, (atile.wd, atile.ht), _ll, _ur);
    DoClip (_ts);
    orto (active_plane, Clipped (_ts) _g);
  fi
  c
enddef;

if unknown segment_split : segment_split := 8; fi
if unknown dashsize :     dashsize  := 3bp; fi
if unknown dashgap :      dashgap := dashsize + 2penwd; fi
if unknown dash_finish :   dash_finish := .5; fi
if unknown dash_start :    dash_start := .5; fi
if unknown unit_of_length : unit_of_length := 0.1in; fi

vardef gendashed (suffix pat) expr f =
  save _g; path _g; _g := zconv(f);
  if (unknown pat.rep) : % no "pattern"
    GBmsg "Dash pattern " & str pat & " undefined. " &
      "Path will be drawn instead.";
    safedraw _g;
  elseif pat.rep < 2 : % no "spaces"
    safedraw _g;
  else:
    save _dl, _tmppat;
    forsuffixes _s = start, rep, finish :
      _dl._s := 0; _tmppat._s := pat._s;
      for i = 1 upto pat._s :
        _tmppat._s[i] := pat._s[i]/unit_of_length;
        _dl._s := _dl._s + _tmppat._s[i];
      endfor
    endfor
    if _dl.rep = 0 :
      GBmsg "Dash pattern " & str pat & " has length 0. " &
      "Path will be drawn instead.";
      safedraw _g;
    else:
      save _p; path _p;
      _p := _g scaled (1/unit_of_length);
      save _cumlen, _totlen, _n, _sf;
      _totlen := makelengtharray(_cumlen) _p;
      _sf := scale_adjust (_n, _dl)(_totlen);
      if _n < 0 : safedraw _g;
      else:
        forsuffixes _s = start, rep, finish :
          for _i = 1 upto _tmppat._s :
            _tmppat._s[_i] := _tmppat._s[_i]*_sf;
          endfor
          _dl._s := _dl._s*_sf;
        endfor
        save dashingdot; picture dashingdot;
        if known plot_pic :  dashingdot := makesymbol(plot_pic, penwd);
        else:  dashingdot := makesymbol(dotpath,  penwd);
        fi
        save _ct, _t, _d, _v;
        picture _v; _v := nullpicture;
        _ct := 0;
        %   Begin with pat.start
        _d0 := 0; _t0 := 0;
        dashit (_tmppat.start) (_v);
        % then pat.rep
        if _n > 0 :
          save _m; _m := ceiling sqrt(_n);
          for _j = 0 step _m until _n - 1 :
            for _i = 0 upto _m - 1 :
              exitif (_i + _j) > _n - 1;
              _d0 := _dl.start + (_j + _i)*_dl.rep;
              _t0 := gettime(_cumlen, _ct) (_d0);
              dashit (_tmppat.rep) (_v);
            endfor
            % add _m patterns and reset.
            DoClip(_v);
            mono (_v)
            coloraddto (drawcolor) (active_plane, _v);
            _v := nullpicture;
          endfor
        fi
        %   and finally, pat.finish
        _d0 := _totlen - _dl.finish;
        _t0 := gettime(_cumlen, _ct) (_d0);
        dashit (_tmppat.finish) (_v);
        DoClip(_v);
        mono(_v)
        coloraddto (drawcolor) (active_plane, _v);
      fi
    fi
  fi
  f
enddef;

vardef makelengtharray (suffix clen) expr p =
  save _s; _s := emax (1, ceiling segment_split);
  clen := _s*length p; clen[0] := 0;
  for _i = 1 upto clen :
    clen[_i] := clen[_i-1] + abs (pnt[_i/_s] (p) - pnt[(_i-1)/_s] (p));
  endfor
  clen[clen]
enddef;

vardef scale_adjust (suffix n, pl) (expr lngth) =
  n := (lngth - pl.start - pl.finish)/pl.rep;
  n := if n < 0 : -1 else: round(n) fi;
  lngth/(pl.start + emax(n, 0)*pl.rep + pl.finish)
enddef;

vardef gettime (suffix arr, ct) (expr lngth) =
  save _gtl, _s;
  _s := emax(1, ceiling segment_split);
  _gtl := emax (arr[ct], emin (arr[arr], lngth));
  forever:
    exitif ( (arr[ct] <= _gtl) and (_gtl <= arr[ct+1]) );
    ct := ct + 1; % need to exit *before* incrementing
  endfor
  if arr[ct] = arr[ct+1] :  ct
  else: ( ct + (_gtl - arr[ct]) / (arr[ct+1] - arr[ct]) )
  fi /_s
enddef;

def dashit (suffix pos) (suffix pic) =
  for _k = 1 upto pos:
    if odd _k : % draw a dash of length pos[_k]
      if pos[_k] = 0 : % point required
        _d1 := _d0; _t1 := _t0;
        picdot (pic, dashingdot, pnt [_t0] (_g));
      else:
        _d1 := _d0 + pos[_k];
        _t1 := gettime (_cumlen, _ct) (_d1);
        shpath (pic, drawpen) (subpath (_t0, _t1) of _g);
      fi
    else: % find the start of the next dash:
      _d0 := _d1 + pos[_k];
      _t0 := gettime(_cumlen, _ct) (_d0);
    fi
  endfor
enddef;

def dashpat (suffix pat) (text t) =
  pat.rep := 0;
  for _itm = t:
    pat.rep[incr pat.rep] := _itm;
  endfor;
  if odd (pat.rep) and (pat.rep > 1):
    pat.rep[incr pat.rep] := 0;
  fi
  pat.start := 1;
  pat.start[1] := pat.rep[1]*dash_start;
  for _idx = 2 upto pat.rep :
    pat.start[incr pat.start] := pat.rep[_idx];
  endfor
  pat.finish := 1;
  pat.finish[1] := pat.rep[1]*dash_finish;
enddef;

vardef DASHED (expr dlen, dgap) expr f =
  save dashes; dashpat (dashes) (dlen, dgap);
  gendashed (dashes) f
enddef;

def dashed = DASHED enddef;

vardef doplot (expr spath, sc, dgap) expr f =
  save dots; dashpat (dots) (0, dgap);
  save plot_pic; picture plot_pic;
  plot_pic := makesymbol (spath, sc);
  gendashed (dots) f
enddef;

path dotpath; dotpath := fullcircle;
vardef dotted (expr dsize, dgap) expr f =
  doplot (dotpath, dsize, dgap) f
enddef;

vardef plotnodes (expr symbol, size) expr f =
  save _pln; pair _pln[];
  _pln := 0;
  for _a = 0 upto (length f) if cycle f : - 1 fi :
    _pln[incr _pln] := pnt[_a] (f);
  endfor
  dosymbols (drawcolor, symbol, size) (_pln);
  f
enddef;

vardef makesymbol (expr spath, sc) =
  if path spath        : setdot (spath, sc)
  elseif picture spath :
    save v; picture v; v:= spath; mono (v); v
  else:
    GBmsg "Undefined symbol for plotting, using dotpath instead.";
    setdot (dotpath, sc)
  fi
enddef;

path Triangle, Square, Circle, Diamond, Star, Plus, Cross,
    Asterisk, SolidTriangle, SolidSquare, SolidCircle,
    SolidDiamond, SolidStar;

Triangle := (for n = 0 upto 2:
  (up rotated 120n)-- endfor up) scaled .78;
SolidTriangle := Triangle & cycle;

Square := (for n = 0 upto 3:
  dir (90n + 45)-- endfor dir 45) scaled .63;
SolidSquare := Square & cycle;

Circle := halfcircle & halfcircle rotated 180;
SolidCircle := Circle & cycle;

Diamond := (Square rotated 45) xscaled (1/1.2) yscaled 1.2;
SolidDiamond := Diamond & cycle;

Plus := (origin for n = 0 upto 3:
    --(up rotated 90n)--origin endfor) scaled .65;
Cross := Plus rotated 45;
Asterisk := (origin for n = 0 upto 5:
    --(up rotated 60n)--origin endfor) scaled .6;

pair zz;
zz = (whatever)[up, up rotated 144];
zz = (whatever)[up rotated 72, up rotated -72];
Star := (for n = 0 upto 4:
    (up rotated 72n)--(zz rotated 72n)-- endfor up) scaled .84;
SolidStar := Star & cycle;
save zz;

numeric dashtype;
forsuffixes s = start, rep, finish :
  numeric dashtype[].s, dashtype[].s[];
endfor
def defaultdashes =
  dashpat (dashtype0) (0);                           % solid
  dashpat (dashtype1) (3bp, 4bp);                    % dashed
  dashpat (dashtype2) (0, 4bp);                      % dotted
  dashpat (dashtype3) (0, 4bp, 3bp, 4bp);          % dot-dash
  dashpat (dashtype4) (0, 4bp, 3bp, 4bp, 0, 4bp);% dot-dash-dot
  dashpat (dashtype5) (0, 4bp, 3bp, 4bp, 3bp, 4bp);% dot-dash-dash
  dashtype := 6;
enddef;
defaultdashes;

vardef isdashpat suffix pat =
  (knownarray pat.start)  and
  (knownarray pat.finish) and
  (knownarray pat.rep)
enddef;

def setdatadashes (text lst) =
  save dashtype; dashtype := 0;
  forsuffixes _itm = lst :
    if isdashpat _itm :
      forsuffixes _s = start, rep, finish :
        copyarray (_itm._s, dashtype[dashtype]._s);
      endfor
      dashtype := dashtype + 1;
    else: GBmsg "Improper dash pattern in setdatadashes.";
    fi
  endfor
  if dashtype = 0 :
    SetdataWarn "dashes";
    defaultdashes;
  fi
enddef;
def getdashpat expr n = dashtype[n mod dashtype] enddef;

def SetdataWarn expr s =
  GBmsg "Command setdata"& s &"() failed; using defaults."
enddef;

def setdatasymbols (text lst) =
  save pointtype; path pointtype[]; pointtype := 0;
  for _itm = lst :
    if (known _itm) and (path _itm):
      pointtype[pointtype] := _itm;
      pointtype := pointtype + 1;
    else:
      GBmsg "Improper path in setdatasymbols().";
    fi
  endfor
  if pointtype = 0:
    SetdataWarn "symbols";
    defaultsymbols;
  fi
enddef;
def getsymbol expr n := pointtype[n mod pointtype] enddef;

numeric pointtype; path pointtype[];
def defaultsymbols =
  pointtype0 := Circle;
  pointtype1 := Cross;
  pointtype2 := SolidDiamond;
  pointtype3 := Square;
  pointtype4 := Plus;
  pointtype5 := Triangle;
  pointtype6 := SolidCircle;
  pointtype7 := Star;
  pointtype8 := SolidTriangle;
  pointtype := 9;
enddef;
defaultsymbols;

vardef bpoint (expr ptwd, b) =
  fullcircle scaled ptwd shifted b
enddef;
def pointd (expr ptwd, filled) (text t) =
  if filled :
    plotsymbol (SolidCircle, ptwd) (t);
  else :
    begingroup;
      save clearsymbols; boolean clearsymbols; clearsymbols := true;
      plotsymbol (Circle, ptwd) (t);
    endgroup
  fi
enddef;

boolean clearsymbols; clearsymbols := false;
vardef clearable (expr pth) =
  false
  if path pth :
  if (not cycle pth) and (length pth > 0):
  if ( pnt0 (pth) = pnt[length pth] (pth) ) :
    or true
  fi fi fi
enddef;

def plotsymbol = colorplotsymbol (pointcolor) enddef;
vardef colorplotsymbol (expr clr, spath, sc) (text t) =
  save _cpls;
  textpairs (_cpls) (t);
  dosymbols (clr, spath, sc) (_cpls);
enddef;

vardef dosymbols (expr clr, spath, sc) (suffix arr) =
  save one_symbol, _pls; picture one_symbol, _pls;
  if clearsymbols and  clearable (spath):
    addsymbols (background, spath&cycle, sc) (arr);
  fi
  addsymbols (clr, spath, sc) (arr);
enddef;

def addsymbols (expr clr, spath, sc) (suffix arr) =
  one_symbol := makesymbol (spath, sc);
  _pls := nullpicture;
  for _idx = 1 upto arr:
    picdot (_pls, one_symbol, zconv(arr[_idx]));
  endfor
  DoClip (_pls);
  mono (_pls);
  coloraddto (clr) (active_plane) (_pls);
enddef;

vardef lclosed expr f = f  if not cycle f : --cycle fi enddef;

numeric default_tension;
default_tension := 1;

def sclosed  = sclosedt (default_tension) enddef;
vardef sclosedt (expr t) expr f =
  if cycle f :  f
  else:  save n; n := length f;
    if n = 0     : f&cycle
    elseif n = 1 : f..tension t..cycle
    else         :
      (pnt0 (f)) { (pnt1(f)) - (pnt[n] (f)) }..tension t
         ..(subpath (1, n-1) of f)..tension t
         ..(pnt[n](f)) { pnt0(f) - pnt[n-1](f) }
         ..tension t..cycle
    fi
  fi
enddef;

def bclosed = bclosedt (default_tension) enddef;
vardef bclosedt (expr t) expr f =
  f  if not cycle f : ..tension t..cycle fi
enddef;

def uclosed  = bclosed enddef;

def bsplinecontrols (suffix b) expr f =
  b  := 4;
  b1 := 2[pre 1(f), post0(f)];
  b2 := 2[post0(f), pnt 0(f)];
  b3 := 2[b1, b2];
  b4 := 2[b2, b3];
enddef;

vardef cbclosed expr f =
  if cycle f : f
  elseif (length f)=0 : f&cycle
  else:
    save p, q; pair p[], q[];
    bsplinecontrols (p) f;         % defines p1 to p4
    bsplinecontrols (q) reverse f; % defines q1 to q4
    f..controls q2 and q3..opencbs (q1,q4,p4,p1)
     ..controls p3 and p2..cycle
  fi
enddef;

vardef qbclosed expr f =
  if cycle f : f
  else: save n; n := length f;
    if n = 0 : f&cycle
    else:
      save p; pair p[]; p := 4;
      p1 := (3/2)[pnt[n](f), pre[n](f)];
      p2 := 2[p1, pnt[n](f)];
      p4 := (3/2)[pnt 0 (f), post0 (f)];
      p3 := 2[p4, pnt 0 (f)];
      f & mkqbs (p) & cycle
    fi
  fi
enddef;

vardef makesector expr p =
  (pathcenter p)--p--cycle
enddef;

path cuttings;
vardef cutoffbefore (expr b) expr f =
  save w, t, u, n; n:= length f;
  pair w;
  for k = 1 upto n :
    w := (subpath (0,k) of f) intersectiontimes b;
    exitif w > left;
  endfor
  if debug :
    GBdebug;
      >> "Intersectiontimes:";
      show w;
    GBenddebug;
  fi
  t := xpart w;
  if t < 0:
    cuttings := pnt0 (f);
    f
  else:
    cuttings := subpath (0,t) of f;
    subpath (t, n) of f
  fi
enddef;

vardef cutoffafter (expr b) expr f =
  save g; path g;
  g := cutoffbefore (b) reverse f;
  cuttings := reverse cuttings;
  reverse g
enddef;

vardef trimmedpath (expr btrim, etrim) expr f =
  save g, h; path g, h;
  g := invvconv (fullcircle scaled 2btrim) shifted pnt0(f);
  h := invvconv (fullcircle scaled 2etrim) shifted pnt[length f] (f);
  cutoffafter (h) cutoffbefore (g) f
enddef;

vardef predirection@# (expr p) =
  - postdirection[length p - @#] (reverse p)
enddef;

vardef postdirection@# (expr p) =
  save _n; _n := length (p);
  save v; pair v; v := __dir (subpath (@#, @# + _n) of p);
  if v = origin :
    v := - __dir (subpath (@#, @# - _n) of p);
  fi
  v
enddef;

vardef __dir (expr p) =
  save v, w; pair v, w; w := pnt0 (p);
  v := origin;
  for n = 1 upto length (p) :
    v := post[n-1] (p) - w;
    if v = origin :
      v := pre[n] (p) - w;
      if v = origin :
        v := pnt[n] (p) - w;
      fi
    fi
    exitif v <> origin;
  endfor
  v
enddef;

newinternal hdwdr, hdten;
interim hdwdr := 1; interim hdten := 1;
boolean hfilled; hfilled := false;

def headshape (expr wr, tens, fil) =
 interim hdwdr := wr; interim hdten := tens;
 save hfilled; boolean hfilled; hfilled := fil;
enddef;

def head = ahead (headcolor) enddef;

vardef ahead (expr clr, front, back, hwr, tens, filled) =
  if front <> back :
    save side; pair side;
    side := (hwr/2) * ((front-back) rotated 90);
    save f; path f;
    f := (back + side)..tension tens..
         {front-back}front{back-front}..tension tens..
         (back - side)
      if filled : --cycle; colorsafefill (clr) f fi;
    colorsafedraw (clr) f;
  fi
enddef;

def headpath = colorheadpath (headcolor) enddef;
vardef colorheadpath (expr clr, hlen, hrot, hback) expr f =
  if hlen <> 0 :
    save g; path g; g := zconv (f);
    save P; pair P[];
    P2 := pnt[length g] (g);
    P1 := predirection[length g] (g);
    if P1 <> (0, 0) :
      P3 := (unitvector P1) rotated hrot;
      P4 := P2 - (hback * P3);
      P5 := P4 - (hlen  * P3);
      ahead (clr, P4, P5, hdwdr, hdten, hfilled);
    fi
  fi
  f
enddef;

def arrowdraw (expr hlen) (expr f) =
  store (curpath) headpath (hlen, 0, 0) drawn f;
enddef;

def xaxis (expr hlen) = arrowdraw (hlen) ((xneg, 0)--(xpos, 0)); enddef;
def yaxis (expr hlen) = arrowdraw (hlen) ((0, yneg)--(0, ypos)); enddef;
def axes (expr hlen) = xaxis (hlen);  yaxis (hlen); enddef;

laxis := baxis := raxis := taxis := 0;

vardef axisline.x = (xneg + laxis, 0)--(xpos - raxis, 0) enddef;
vardef axisline.y = (0, yneg + baxis)--(0, ypos - taxis) enddef;
vardef axisline.l = axisline.y shifted (xneg + laxis, 0) enddef;
vardef axisline.b = axisline.x shifted (0, yneg + baxis) enddef;
vardef axisline.r = axisline.y shifted (xpos - raxis, 0) enddef;
vardef axisline.t = axisline.x shifted (0, ypos - taxis) enddef;
vardef axis@# (expr len) =
  headpath (len, 0, 0) axisline@#
enddef;

numeric inside, outside, centered, onleft, onright, ontop, onbottom;
inside   := -2;
outside  := -1;
onright  :=  1;
onleft   :=  2;
centered := .5[onright, onleft];
onbottom := onright;
ontop    := onleft;
ltick := rtick := ttick := btick := inside;
xtick := ytick := centered;

vardef axismarks (expr inang, tp, loc, pdir) (expr len) (text t) =
 save _tp, _U, _P, _tic, _ticang;
 pair _U, _P, _tic[];
 _ticang := if tp<0 : inang else: 90 fi;
 _tp := abs(tp) - 1;
 _U := unitvector (vconv (pdir)) rotated _ticang;
 _tic1 := (_tp - 1) * len * _U; % start of mark
 _tic2 := _tp * len * _U;       % end of mark
 for _a = t:
   safedraw ((_tic1--_tic2) shifted zconv (loc + _a*pdir));
 endfor
enddef;

def xmarks =  axismarks ( 90, xtick, (0, 0), right) enddef;
def ymarks =  axismarks (-90, ytick, (0, 0), up)    enddef;
def lmarks =  axismarks (-90, ltick, (xneg + laxis, 0), up)    enddef;
def bmarks =  axismarks ( 90, btick, (0, yneg + baxis), right) enddef;
def rmarks =  axismarks ( 90, rtick, (xpos - raxis, 0), up)    enddef;
def tmarks =  axismarks (-90, ttick, (0, ypos - taxis), right) enddef;

path griddotpath; griddotpath := fullcircle;
def grid = vgrid (0.5bp) enddef;
vardef vgrid (expr dsize, xspace, yspace) =
  save gdot, gridpic; picture gdot, gridpic;
  gdot := setdot (griddotpath, dsize);
  gridpic := nullpicture;
  for n = ceiling(xneg/xspace) upto floor(xpos/xspace):
  for m = ceiling(yneg/yspace) upto floor(ypos/yspace):
      picdot (gridpic, gdot, zconv((n*xspace, m*yspace)));
  endfor
  endfor
  mono (gridpic);
  coloraddto (pointcolor) (active_plane) (gridpic);
enddef;

def hgridlines (expr ysp) =
  for n = ceiling((yneg + baxis)/ysp) upto floor((ypos - taxis)/ysp) :
    safedraw zconv((xneg + laxis, n*ysp)--(xpos - raxis, n*ysp));
  endfor
enddef;
def vgridlines (expr xsp) =
  for n = ceiling((xneg + laxis)/xsp) upto floor((xpos - raxis)/xsp) :
    safedraw zconv((n*xsp, yneg + baxis)--(n*xsp, ypos - taxis));
  endfor
enddef;
def gridlines (expr xsp, ysp) =
  vgridlines (xsp); hgridlines (ysp);
enddef;

vardef plrpatch (expr rstart, rstop, rstep, tstart, tstop, tstep) =
  save v; picture v; v := nullpicture;
  patcharcs (v) (rstart, rstop, rstep, tstart, tstop);
  coloraddto (drawcolor) (active_plane, v);
  v := nullpicture;
  patchrays (v) (tstart, tstop, tstep, rstart, rstop);
  coloraddto (drawcolor) (active_plane, v);
enddef;
def patcharcs (suffix X) (expr rstart, rstop, rstep, tstart, tstop) =
  for rad = (if rstart=0: rstep else: rstart fi) step rstep until rstop:
    orto (X, picpath zconv (arcplr ((0, 0), tstart, tstop, rad)) );
  endfor
enddef;
def patchrays (suffix X) (expr tstart, tstop, tstep, rstart, rstop) =
  for _ang = tstart step tstep until tstop:
    orto (X) (picpath zconv ((rstart*dir _ang)--(rstop*dir _ang)));
  endfor
enddef;

def polargrid (expr rstep, tstep) =
  gridarcs (rstep); gridrays (tstep);
enddef;

def polargridpoints (expr dsize, rstep, tstep) =
  beginpolargrid;
    save gdot; picture gdot; gdot := setdot (griddotpath, dsize);
    if rmin = 0:
      picdot (gridpic, gdot, zconv(origin));
      rmin := rstep;
    fi
    for n = ceiling (rmin/rstep) upto floor (rmax/rstep) :
    for m = ceiling (tmin/tstep) upto floor (tmax/tstep) :
      picdot ( gridpic, gdot, zconv ( polar((n*rstep, m*tstep)) ) );
    endfor
    endfor
  endpolargrid (pointcolor, .5dsize);
enddef;

def gridarcs (expr rstep) =
  beginpolargrid;
    if rmin = 0 :  % add "circle" of radius 0
      picdot (gridpic, setdot(griddotpath, penwd), zconv(origin));
    fi
    rmin := rstep * floor(rmin/rstep + 1);
    rmax := rstep*ceiling(rmax/rstep - 1);
    patcharcs (gridpic) (rmin, rmax, rstep, tmin, tmax);
  endpolargrid (drawcolor, .5penwd);
enddef;

def gridrays (expr tstep) =
  beginpolargrid;
    tmin := tstep*ceiling(tmin/tstep);
    tmax := tstep * floor(tmax/tstep);
    patchrays (gridpic) (tmin, tmax, tstep, rmin, rmax);
  endpolargrid (drawcolor, .5penwd);
enddef;

def beginpolargrid =
  begingroup;
    save p, r, t, rmax, rmin, tmax, tmin;
    pair p[];
    % Four corners:
    p0 := (xneg, yneg); p1 := (xneg, ypos);
    p2 := (xpos, ypos); p3 := (xpos, yneg);
    r0 := abs(p0); rmax := r0;
    for j = 1 upto 3 :
      r[j] := abs(p[j]);
      if rmax < r[j] : rmax := r[j]; fi
    endfor
    rmin := 0;
    if (xneg < 0) and (xpos > 0) and (yneg < 0) and (ypos > 0) :
      tmin := 0; tmax := 360;
    elseif (p0 = (0,0)) : tmin :=    0; tmax :=  90;
    elseif (p1 = (0,0)) : tmin :=  -90; tmax :=   0;
    elseif (p2 = (0,0)) : tmin := -180; tmax := -90;
    elseif (p3 = (0,0)) : tmin :=   90; tmax := 180;
    else :
      tmax := tmin := t0 := angle p0;
      for j = 1 upto 3:
        t := t0 + angle (p[j] rotated -t0);
        if tmax < t : tmax := t; fi
        if tmin > t : tmin := t; fi
      endfor
      if (xneg < 0) and (xpos > 0) :        % (1)
        rmin := emin(abs(yneg), abs(ypos));
      elseif (yneg < 0) and (ypos > 0) :    % (2)
        rmin := emin(abs(xneg), abs(xpos));
      else :                                % (3)
        rmin := min(r0, r1, r2, r3);
      fi
    fi
    save gridpic; picture gridpic; gridpic := nullpicture;
enddef;


def endpolargrid (expr clr, size)=
    clipto (gridpic) rect ( zconv((xneg, yneg)) - size*(1,1),
                            zconv((xpos, ypos)) + size*(1,1) );
    coloraddto (clr) (active_plane) (gridpic);
  endgroup
enddef;

vardef polarpatch (expr rstart, rstop, rstep, tstart, tstop, tstep) =
  plrpatch (rstart, rstop, rstep, tstart, tstop, tstep);
  safedraw zconv ( arcplr ((0, 0), tstart, tstop, rstop) );
  safedraw zconv ( ((rstart, 0)--(rstop, 0)) rotated tstop );
enddef;

vardef rect (expr ll, ur) =
  ll--(xpart ur, ypart ll)--ur--(xpart ll, ypart ur)--cycle
enddef;
vardef triangle (expr A, B, C) = A--B--C--cycle enddef;

vardef regularpolygon (expr n) (suffix Bob) (text eqns) =
  pair Bob[]; Bob := emax(round (abs (n)), 2);
  eqns;
  for _uncle = 1 upto Bob - 1 :
    (Bob1 - Bob0) rotated (360*_uncle/Bob) = Bob[_uncle+1] - Bob0;
  endfor
  mkpoly (true) (Bob)
enddef;

vardef altitudept expr n of t =
  save A, B, C, zz; pair A, B, C, zz;
  A := pnt[n]     (t);
  B := pnt[n + 1] (t); % wraps around a cyclic path
  C := pnt[n + 2] (t);
  zz = whatever[B,C];
  zz = A + whatever*((C-B) rotated 90);
  zz
enddef;

vardef altitude expr n of t =
  (pnt[n](t))--(altitudept n of t)
enddef;

vardef medianpt expr n of t =
  0.5[pnt[n + 1] (t), pnt[n + 2] (t)]
enddef;

vardef median expr n of t =
  (pnt[n](t))--(medianpt n of t)
enddef;

vardef anglebisectorpt expr n of t =
  save A, B, C; pair A, B, C;
  A := pnt[n    ] (t);
  B := pnt[n + 1] (t);
  C := pnt[n + 2] (t);
  save zz; pair zz;
  zz = whatever[B,C];
  zz = A + whatever*((B-A) rotated (.5*cornerangle (A,B,C)));
  zz
enddef;

vardef anglebisector expr n of t =
  (pnt[n](t))--(anglebisectorpt n of t)
enddef;

vardef cornerangle (expr A, B, C) =
  if (A = B) and (B = C)    : 60
  elseif (B = C)            :  0
  elseif (A = B) or (A = C) : 90
  else:  angle ((C - A) rotated (-angle (B - A)))
  fi
enddef;

vardef mkpath (expr smooth, tens, cyclic) (suffix pts) =
  if smooth :  mksmooth (tens, cyclic, pts)
  else      :  mkpoly (cyclic, pts)
  fi
enddef;

vardef mkpoly (expr cyclic) (suffix pts) =
  for _i = 1 upto pts-1:  pts[_i]-- endfor
  pts[pts] if cyclic : -- cycle fi
enddef;

vardef polyline (expr cyclic) (text t) =
  save _pl; textpairs (_pl) (t);  mkpoly (cyclic, _pl)
enddef;

vardef mksmooth (expr tens, cyclic) (suffix pts) =
  pts1
  if pts = 1 :
    if cyclic : &cycle fi
  else:
    if cyclic :
      {pts[2]-pts[pts]}
    fi
    for _i = 2 upto pts-1:
      ..tension tens..pts[_i]{pts[_i+1]-pts[_i-1]}
    endfor
    ..tension tens..pts[pts]
    if cyclic :
      {pts[1]-pts[pts-1]}..tension tens..cycle
    fi
  fi
enddef;

def curve = tcurve (default_tension) enddef;
vardef tcurve  (expr tens, cyclic) (text t) =
  save _tc;  textpairs (_tc) (t);  mksmooth (tens, cyclic, _tc)
enddef;

vardef mkbezier (expr tens, cyclic) (suffix pts) =
  for _i = 1 upto pts-1 : pts[_i]..tension tens.. endfor
    pts[pts] if cyclic : ..tension tens..cycle fi
enddef;

def bezier = tbezier (default_tension) enddef;
vardef tbezier (expr tens, cyclic) (text t) =
  save _tsb;  textpairs (_tsb) (t);  mkbezier (tens, cyclic) (_tsb)
enddef;

vardef mkqbezier (expr cyclic) (suffix pts) =
  save _mqb;  _mqb := pts;
  if (cyclic and odd pts) or not (cyclic or odd pts):
    pts[incr _mqb] := pts[pts];
  fi
  if cyclic : pts[incr _mqb] := pts1; fi
  pts1
  for _i = 2 step 2 until _mqb - 1 :
    ..controls  1/3[pts[_i],pts[_i-1]] and 1/3[pts[_i], pts[_i+1]]
    ..pts[_i+1]
  endfor
  if cyclic : &cycle fi
enddef;

vardef qbezier (expr cyclic) (text t) =
  save _qbz;  textpairs (_qbz) (t);  mkqbezier (cyclic) (_qbz)
enddef;

vardef openqbs (text t) =
  save _oq;  textpairs (_oq) (t);  mkqbs (_oq)
enddef;

vardef closedqbs (text t) =
  save _cq; textpairs (_cq) (t);
  _cq[incr _cq] := _cq1;  _cq[incr _cq] := _cq2;
  mkqbs (_cq) & cycle
enddef;

vardef mkqbs (suffix b) =
  for _i = 1 upto b-2:
    0.5[b[_i], b[_i+1]]
      ..controls 1/6[b[_i+1], b[_i]] and 1/6[b[_i+1], b[_i+2]]..
  endfor
  0.5[b[b-1], b[b]]
enddef;

vardef mkopencbs (suffix b) =
  for _i = 1 upto b-3:
    (b[_i]+4b[_i+1]+b[_i+2])/6
      ..controls 1/3[b[_i+1], b[_i+2]] and 2/3[b[_i+1], b[_i+2]]..
  endfor
  (b[b-2]+4b[b-1]+b[b])/6
enddef;
vardef mkclosedcbs (suffix b) =
  mkopencbs (b) & opencbs (b[b-2],b[b-1],b[b], b1, b2, b3) & cycle
enddef;

vardef opencbs (text t) =
  save _oc;  textpairs (_oc) (t);  mkopencbs (_oc)
enddef;
vardef closedcbs (text t) =
  save _clc; textpairs (_clc) (t);  mkclosedcbs (_clc)
enddef;

vardef fcncontrol (expr ftens, X, Y, Z) =
  save dl, dr, before, after;  pair before, after;
  before := Y - X;  after := Z - Y;
  dl := xpart (before);  dr := xpart (after);
  if (dr = 0) or (dl = 0):
    Y + abs(dr)/ftens * sgn before
  else:
    Y + abs(dr)/ftens * unitvector (before*dr/dl + after*dl/dr)
  fi
enddef;

vardef  mkfcnpath (expr ftens) (suffix q) =
  for _i = 1 upto q - 1:
    q[_i]..controls fcncontrol (ftens) (q[_i-1], q[_i],  q[_i+1])
           and fcncontrol (ftens) (q[_i+2], q[_i+1], q[_i])..
  endfor
  q[q]
enddef;

def fcncurve = functioncurve (emax(1.2default_tension, eps)) enddef;
vardef functioncurve (expr ftens) (text t) =
  save _fc; textpairs (_fc) (t);
  if _fc > 1 :  _fc0 :=  _fc1;  _fc[_fc+1] := _fc[_fc]; fi
  mkfcnpath (ftens)(_fc)
enddef;

vardef turtle (text t) =
  save _tu; pair _tu[]; _tu := 0; _tu0 := (0, 0);
  for _a = t: _tu[incr _tu] := _tu[_tu - 1] + _a; endfor
  mkpoly (false, _tu)
enddef;

vardef mkarc (expr center, from, to, sweep) =
  save n, d; pair d;
  n := ceiling (abs(sweep)/45);
  d := (from - center) rotated (signof (sweep) 90);
  from{d}
  for j = 1 upto n-1 :
    ..(from rotatedabout (center, j/n*sweep)){d rotated (j/n*sweep)}
  endfor ..to{d rotated sweep}
enddef;

vardef arc (expr center, from, sweep) =
  if (center = from) or (sweep = 0) :
    from--from
  else:
    save to; pair to;
    to := from rotatedabout (center, sweep);
    mkarc (center, from, to, sweep)
  fi
enddef;
def arccps = arc enddef;

vardef arccenter (expr from, to, sweep) =
  save ang, c;
  pair c;
  ang := 90 - (sweep mod 360)/2; %  -90 < ang <= 90
  if (abs(ang) = 90) or (from = to) :
    GBmsg "The central point of this arc is undefined. " &
      "Using midpoint of chord instead.";
    0.5[from, to]
  else:
    save cd; pair cd;  cd := to - from;
    c = from + whatever*(cd rotated ang);
    if abs(ang) < 30 :
      c = (0.5)[from, to] + whatever*(cd rotated 90);
    else:
      c = to + whatever*(-cd rotated -ang);
    fi
    c
  fi
enddef;

vardef midarc (expr from, to, sweep) =
  save m, cd; pair m, cd;
  cd := to - from;
  m = from + whatever*( cd rotated (-sweep/4));
  m = 0.5[from, to] + whatever*(cd rotated 90);
  m
enddef;

vardef arcpps (expr from, to, sweep) =
  if ((sweep mod 360) = 0) or (from = to) :
    GBmsg "Undefined arc. A line segment will be used instead.";
    from--to
  elseif abs(sweep) <= 90 :
    save cd; pair cd; cd := to - from;
    if abs(sweep) <= 45 :
      from{cd rotated (-sweep/2)}..to{cd rotated (sweep/2)}
    else:
      from{cd rotated (-sweep/2)}..midarc(from, to, sweep){cd}
        ..to{cd rotated (sweep/2)}
    fi
  else:
    save center; pair center;
    center := arccenter (from, to, sweep);
    mkarc (center, from, to, sweep)
  fi
enddef;

vardef arcplr (expr center, frtheta, totheta, rad) =
  if rad = 0 :
    center -- center
  else:
    save from, to; pair from, to;
    from := center + rad*dir frtheta;
    to   := center + rad*dir totheta;
    if frtheta = totheta :
      from--to
    else:
      mkarc (center, from, to, totheta - frtheta)
    fi
  fi
enddef;
vardef arcalt (expr center, radius, anglefrom, angleto) =
  arcplr (center, anglefrom, angleto, radius)
enddef;

vardef arcppp (expr first, second, third) =
  arcpps (first, second, 2*cornerangle (third, first, second)) &
  arcpps (second, third, 2*cornerangle (first, second, third))
enddef;

vardef ellipse (expr center, radx, rady, angle) =
  fullcircle xscaled (2*radx) yscaled (2*rady) rotated angle
    shifted center
enddef;
vardef circle (expr center, rad) =
  fullcircle scaled (2*rad) shifted center
enddef;

vardef circlecp (expr center, point) =
  mkarc (center, point, point, 360) & cycle
enddef;
vardef circleppp (expr one, two, three) =
  save ang; numeric ang[];
  ang0 := cornerangle(three, one, two);
  ang1 := cornerangle(one, two, three);
  ang2 := cornerangle(two, three, one);
  arcpps (one, two, 2ang0) & arcpps (two, three, 2ang1) &
    arcpps (three, one, 2ang2) & cycle
enddef;
vardef circlepps (expr one, two, sweep) =
  save ang, full; numeric ang[], full;
  full := signof (sweep) 360;
  ang1 := sweep mod (full);
  ang2 := full - ang1;
  arcpps (one, two, ang1) & arcpps (two, one, ang2) & cycle
enddef;

vardef pathcenter expr p =
  save a, cntr, n; pair cntr, a[];
  n := length p;
  a1 = pnt   0   (p);
  a3 = pnt [n/2] (p);
  if cycle p :
    a2 = pnt [n/4] (p);
    a4 = pnt [3n/4] (p);
  else:
    a2 := a3;
    a4 := pnt[n] (p);
  fi
  cntr = .5[a1, a3] + whatever*((a3 - a1) rotated 90);
  cntr = .5[a2, a4] + whatever*((a4 - a2) rotated 90);
  cntr
enddef;

vardef circumcircle expr t =
  circleppp (pnt0 (t), pnt1 (t), pnt2 (t))
enddef;

vardef incircle expr t =
  save A, B, C; pair A, B, C;
  A := pnt0 (t);
  B := pnt1 (t);
  C := pnt2 (t);
  % Find the tangent points on the sides. E.g., a is the common
  % distance from A to the tangent points on the adjacent sides.
  save a, b, c;
  a + b     = abs (B-A);
      b + c = abs (C-B);
  a     + c = abs (A-C);
  circleppp (A + a*unitvector (B-A),
             B + b*unitvector (C-B),
             C + c*unitvector (A-C))
enddef;

vardef excircle expr n of t =
  save A, B, C; pair A, B, C;
  A := pnt[n]     (t);
  B := pnt[n + 1] (t); % wraps around
  C := pnt[n + 2] (t);
  save a, b, c;
  a - b     = abs (B-A);
      b + c = abs (C-B);
  a     - c = abs (C-A);
  circleppp (A + a*unitvector(B-A),
             B + b*unitvector(C-B),
             C + c*unitvector(C-A))
enddef;

vardef ninepointcircle expr t =
  circleppp (medianpt 0 of t, medianpt 1 of t, medianpt 2 of t)
enddef;

vardef circumcenter expr t = pathcenter circumcircle t enddef;
vardef incenter expr t = pathcenter incircle t enddef;
vardef excenter expr n of t = pathcenter excircle n of t enddef;
vardef ninepointcenter expr t = pathcenter ninepointcircle t enddef;

vardef barycenter expr t =
  save n, m; n := length t; m := n + 1;
  save xxx;
  xxx : = pnt0 (t)/m  for k = 1 upto n-1 :  + pnt[k] (t)/m  endfor;
  if cycle t:  xxx*(1 + 1/n)
  else:  xxx + pnt[n] (t)/m
  fi
enddef;
vardef sector (expr center, rad, frtheta, totheta) =
  center -- arcalt (center, rad, frtheta, totheta) -- cycle
enddef;

vardef mkfcn (expr smooth, tens) (expr bmin, bmax, bst) (text pf) =
  save _p; pair _p[]; _p := 0;
  save _dx, _n, _r; numeric _dx, _n, _r;
  if bmax = bmin :  _n := 1;
  else:
    _r := bmax - bmin;
    _dx := max (abs(bst), nottoosmall*abs(_r), epsilon);
    _n := emax (round(abs(_r)/_dx), 1);
  fi
  for _i = 0 upto _n:  _p[incr _p] := pf(bmin + _i/_n*_r);  endfor
  mkpath (smooth, tens, false, _p)
enddef;
def tfcn (expr smooth) = mkfcn (smooth, default_tension) enddef;

def parafcn (expr smooth) = tparafcn (smooth, default_tension) enddef;
vardef tparafcn (expr sm, tn) (expr bmin, bmax, bst) (text pf) =
  save _fp;  vardef _fp (expr t) = pf enddef;
  mkfcn (sm, tn) (bmin, bmax, bst) (_fp)
enddef;

vardef xfcn (expr smooth) (expr xmin, xmax, st) (text _fx) =
  save _fp;  vardef _fp (expr _x) = (_x, _fx(_x)) enddef;
  mkfcn (smooth, default_tension) (xmin, xmax, st) (_fp)
enddef;

def function (expr smooth) = tfunction (smooth, default_tension) enddef;
vardef tfunction (expr smooth, tens) (expr xmin, xmax, st) (text _fx) =
  save _fp;  vardef _fp (expr x) = (x, _fx) enddef;
  mkfcn (smooth, tens) (xmin, xmax, st) (_fp)
enddef;

def btwnfcn (expr sm) = tbtwnfcn (sm, default_tension) enddef;
vardef tbtwnfcn (expr sm, tn)(expr xlo, xhi, st)(text _fx)(text _gx) =
  tfunction (sm, tn) (xlo, xhi, st) (_fx) --
  ( reverse tfunction (sm, tn) (xlo, xhi, st) (_gx) ) -- cycle
enddef;

vardef rfcn (expr smooth) (expr tmin, tmax, st) (text ft) =
  save _fq;  vardef _fq (expr t) = (ft(t)) * (dir t) enddef;
  mkfcn (smooth, default_tension) (tmin, tmax, st) (_fq)
enddef;

def plrfcn (expr smooth) = tplrfcn (smooth, default_tension) enddef;
vardef tplrfcn (expr smooth, tens) (expr tmin, tmax, st) (text ft) =
  save _fq;  vardef _fq (expr t) = (ft) * (dir t) enddef;
  mkfcn (smooth, tens) (tmin, tmax, st) (_fq)
enddef;

vardef piechart (expr sign, ang, cent, rad) (text data) =
  save _sum, _tot;
  numeric piewedge; piewedge := 0;
  numeric pieangle, pieangle[]; pieangle0 := 0;
  for _val = data :
    pieangle[incr piewedge] := pieangle[piewedge - 1] + _val;
  endfor
  _tot := pieangle[piewedge];
  pair piecenter; piecenter := cent;
  path piewedge[];
  numeric piedirection; pair piedirection[];
  pieangle[piewedge + 1] = ang + sign*360;
  for _n = piewedge downto 1 :
    pieangle[_n] := ang + sign*pieangle[_n - 1]/_tot*360;
    piewedge[_n] =
        sector(cent, rad, pieangle[_n], pieangle[_n+1]);
    piedirection[_n] := dir(0.5[ pieangle[_n], pieangle[_n+1] ]);
  endfor
  piedirection := pieangle := piewedge;
enddef;

def barchart (expr start, sep, r, vert)(text data) =
  numeric barbegin, barbegin[],
          barend, barend[],
          barlength, barlength[],
          barstart, barstart[],
          chartbar, barwd;
  path chartbar[];
  chartbar := 0; barwd := r*sep;
  for _itm = data :
    barend[incr chartbar] := if pair _itm: ypart _itm else: _itm fi;
    barbegin[chartbar]    := if pair _itm: xpart _itm else:   0  fi;
  endfor
  barbegin := barend := barlength := barstart := chartbar;
  for _n = 1 upto chartbar :
    barstart[_n]  := start + sep*(_n-1);
    barlength[_n] := barend[_n];
    chartbar[_n]  :=  rect ((barbegin[_n], 0), ( barend[_n], barwd) )
          shifted (0, barstart[_n]) if vert: xyswap fi;
  endfor
enddef;

picture totalpicture;
boolean totalnull, currentnull;
def clearit =
  currentpicture := totalpicture := nullpicture;
  currentnull := totalnull := true;
enddef;

def keepit =
  mono (currentpicture);
  addto totalpicture also currentpicture;
  currentpicture := nullpicture;
  totalnull := currentnull;
  currentnull := true;
enddef;

def addto_currentpicture =
  currentnull := false;
  addto currentpicture
enddef;

def mergeit (text do) =
  if totalnull :
    do currentpicture
  elseif currentnull :
    do totalpicture
  else:
    begingroup
      save _v_; picture _v_;
      _v_ := currentpicture;
      mono (_v_);
      addto _v_ also totalpicture;
      do _v_
    endgroup
  fi
enddef;

def shipit_ =
  mergeit (shipout)
enddef;
def shipit = shipit_ enddef;

def showit_ =
  mergeit (show_)
enddef;
def show_ suffix v =
  display v inwindow currentwindow
enddef;

numeric gcode; gcode := 0;

%  end grafbase.mf
endinput.
%%
%% End of file `grafbase.mf'.