path: root/Master/texmf-dist/metafont/mfpic/grafbase.mf

%%
%% This is file `grafbase.mf',
%% generated with the docstrip utility.
%%
%% The original source files were:
%%
%% grafbase.dtx  (with options: `MF')
%% 
%% -------------------------------------------------------------------
%% 
%% Copyright 2002--2012, Daniel H. Luecking
%% 
%% Mfpic may be distributed and/or modified under the conditions of the
%% LaTeX Project Public License, either version 1.3c 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.3c or later is part of all distributions of LaTeX version
%% 2008/12/01 or later.
%% 
%% Mfpic has maintenance status "author-maintained". The Current Maintainer
%% is Daniel H. Luecking. There are several Base Interpreters associated
%% with mfpic: plain TeX, LaTeX, plain Metafont and plain MetaPost.
%% 
if (known grafbaseversion) or (known grafbase):
  message "Grafbase (" & jobname & "): You have loaded grafbase more "
    & "than once! Please make sure that it is loaded only once.";
  endinput;
fi
boolean grafbase; grafbase := true;

string fileversion, filedate;
fileversion := "1.10";  filedate := "2012/12/03";

message " Loading grafbase macros, version " & fileversion & ", " &
  filedate & ".";
message " ";

def GBmsg expr s =  message "Grafbase (" & jobname & "): " & s;
enddef;
def GBwarn expr s =  GBmsg "Warning, " & s;  enddef;
def GBerrmsg (expr s) expr t =  errhelp t;
  errmessage "Grafbase (" & jobname & "): " & s;  errhelp "";
enddef;

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:
    message "";
    GBwarn "Version mismatch: "
         & "mfpic and grafbase versions do not match.";
    message "";
  fi
enddef;

checkversions (110);

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

if not boolean debug:  boolean debug;  debug := false;  fi

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

vardef mftitle expr t =
  if string t:  t;  message t;  fi
enddef;

boolean METAFONT, METAPOST;
METAPOST := known color Carl Philipp Emanuel Bach;
if METAPOST:  METAFONT := false; else: METAFONT := true;  fi

if METAPOST:
  GBerrmsg ("wrong compiler.")
    "This file is for Metafont. For Metapost, use grafbase.mp.";
fi

if unknown mode:
  GBerrmsg ("Metafont mode is unknown.")
    "Set mode to a known mode, perhaps ljfour. "
  & "If you proceed, localfont will be tried. "
  & "If that is unknown, a generic mode will be tried.";
  if known localfont:  mode := localfont;
  else:
    if unknown GBresolution:  GBresolution := 600  fi;
    mode_def GBgeneric =
      mode_param (pixels_per_inch, GBresolution);
      mode_param (blacker, 0);
      mode_param (fillin, 0);
      mode_param (o_correction, 1);
      mode_common_setup_;
    enddef;
    mode := GBgeneric;
  fi
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 warningcheck := 0;
interim turningcheck := 0;

numeric unitlen, xscale, yscale, xneg, xpos, yneg, ypos;

unitlen := 1 bp#;
xscale := 7.2;
yscale := 7.2;
xneg := 0;  xpos := 10;
yneg := 0;  ypos := 10;

newinternal deg, pi, radian;
deg := 1;  pi := 3.14159;
radian := 57.29578;
numeric degree;  degree := deg;

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 ( setboolean (ClipOn) false;  t)
enddef;

boolean showbbox;  showbbox := false;

let color = numeric;  color black, white;
let rgbcolor = numeric;
let cmykcolor = numeric;
black := 0;  white := 1;
def withcolor text t =  enddef;
def _wc_ =  withcolor  enddef;

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

vardef snapto expr t =
  if numeric t:
    if unknown t:  0
    elseif t < 0:  0
    elseif t > 1:  1
    else:          t
    fi
  else:
    GBerrmsg ("Improper expression type.")
      "The argument to `snapto' must be a numeric.";
  fi

enddef;

vardef rgbgray  (expr g) =  (snapto g) * white  enddef;
vardef cmykgray (expr g) =  cmyk(0,0,0,1 - snapto g)  enddef;

vardef knowncolor expr clr =  (known clr) and (iscolor clr)  enddef;

vardef grayscalegray (expr g) =  snapto g  enddef;
vardef gray (expr g) =  grayscalegray (g)  enddef;
vardef cmyk (expr c, m, y, k) =  rgb (1-c-k, 1-m-k, 1-y-k) enddef;
vardef makegray  primary clr =
  if knowncolor clr:  clr  else:  black  fi
enddef;
def makergb  =  makegray  enddef;
def makecmyk =  makegray  enddef;
vardef iscolor expr clr =  (color clr)  enddef;

vardef forceclr (expr c) =
  if unknown c :
    if numeric c:  grayscaleblack
    elseif rgbcolor c:   rgbblack
    elseif cmykcolor c: cmykblack
    else: black
    fi
  elseif numeric c:  gray (c)
  elseif iscolor c:  c
  else:  black
  fi
enddef;
vardef named (suffix c) =  forceclr (c)  enddef;

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

vardef rgbtogray (expr r, g, b) =
  togray(snapto r, snapto g, snapto b)
enddef;
vardef cmyktogray (expr c, m, y, k) =
  rgbtogray (1-c-k,1-m-k,1-y-k)
enddef;

vardef cmyktorgb (expr c,m,y,k) =
  rgb(1-c-k,1-m-k,1-y-k)
enddef;
vardef rgbtocmyk (expr r,g,b) =
  cmyk(1-r,1-g,1-b,0)
enddef;

vardef rgb (expr r, g, b) =
  togray (snapto r, snapto g, snapto b)
enddef;
vardef RGB (expr R, G, B) =
  rgb (R/255, G/255, B/255)
enddef;

def list (suffix v) (text lst) =
  v := 0;  for _itm = lst:  v[incr v] := _itm;  endfor
  if v = 0:
    GBerrmsg ("No list to process!")
      "An attempt was made to produce an array from a "
    & "list of expressions having no valid entries.";
  fi
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 knownnumericarray suffix arr =
  setboolean (_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]) and (numeric arr[_idx]);
    endfor
  fi
  _kna
enddef;

def copyarray (suffix src, dest) =
  for _idx = 1 upto src:  dest[_idx] := src[_idx];  endfor
  dest := src;
enddef;

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

def textpairs =  gsetarray (pair)  enddef;

def setuniquepairs (suffix p) (text t) =
  save p;  pair p[];
  setpairs (_up) (t);
  if _up > 0:
    p := 1;  p1 := _up1;
    for _i = 2 upto _up:
      if _up[_i] <> p[p]:  p[incr p] := _up[_i];  fi
    endfor
  else:
    p := 0;
  fi
enddef;

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

vardef floorpair (expr p) =  (floor (xpart p), floor (ypart p))
enddef;
vardef ceilingpair (expr p) = (ceiling (xpart p), ceiling (ypart p))
enddef;

def hroundpair (expr p) =  (hround (xpart p), hround (ypart p))
enddef;
vardef goodpair (expr p) =  hroundpair(p.t_)  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) =  setpair (_mp) p1;
  for _idx = 2 upto p - 1:  _mp := pairmin (_mp, p[_idx]);  endfor
  pairmin (_mp, p[p])
enddef;
vardef maxpair (suffix p) =  setpair (_mp) p1;
  for _idx = 2 upto p - 1:  _mp := pairmax (_mp, p[_idx]);  endfor
  pairmax (_mp, p[p])
enddef;

primarydef Z xprod W =  (xpart Z * ypart W - xpart W * ypart Z)
enddef;

def force_initial (expr p) (suffix f) =
  hide( setnumeric (_n) length f;
  f := p
    if _n = 0:
      {0,0}
    else:
      ..controls post0 (f) and pre 1 (f).. subpath (1,_n) of f
    fi;)
enddef;

def force_terminal (expr p) (suffix f) =
  hide(setpath (_f) reverse f;
  force_initial (p) (_f);
  f := reverse _f;)
enddef;

def force_equal_ends (suffix f, g) =
  hide(save _p;  pair _p;
  _p := .5[pnt[length f] (f), pnt0(g)];
  force_terminal (_p) (f);  force_initial (_p) (g);)
enddef;

def replace_ends_of_cycle (expr p) (suffix f) =
hide(
if cycle f:
  save _n;  _n := length f;
  f := p
    if _n = 0:  &cycle
    else:  .. controls post0 (f) and pre 1 (f) ..
      if _n = 1:  cycle
      else:  subpath (1, _n - 1) of f ..
        controls post[_n - 1](f) and pre[_n](f) .. cycle
      fi
    fi;
fi)
enddef;

pair thetimes;
numeric _Xtime, _Ytime;
tertiarydef a intersects b =
  begingroup
    thetimes := a intersectiontimes b;
    _Xtime := xpart thetimes;
    _Ytime := ypart thetimes;
    (_Xtime > -1)
  endgroup
enddef;

tertiarydef a misses b =  ((a intersectiontimes b) < origin)  enddef;

vardef makepicture (expr s) =
  if picture s:    s
  elseif path s:   picpath (s)
  else:            nullpicture
  fi
enddef;

vardef onepointpath (expr cyclic, q) =
  q if cyclic:  &cycle  else:  {0,0}  fi
enddef;

vardef fallbackpath (expr cyclic, p) (text t) =
  onepointpath (cyclic, p)
enddef;

def even =  not odd  enddef;

primarydef a divides b =
  ((b mod a) = 0)
enddef;

vardef image (text t) =
  newpicture (currentpicture);
  t;
  currentpicture
enddef;

def beginimage =
  begingroup
    newpicture (currentpicture);
enddef;
def endimage =
  ;  currentpicture
  endgroup
enddef;

def makeimage (suffix name) (expr refpt) =
  setpair (_image_reference_point) zconv (refpt);
  setpicture (name) beginimage
enddef;
def concludeimage  =
  endimage shifted
    -goodpair (_image_reference_point)
enddef;

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

def setnumeric (suffix name) =  save name;  name :=  enddef;
def setboolean   =  setvariable (boolean)  enddef;
def setpair      =  setvariable (pair)  enddef;
def setpath      =  setvariable (path)  enddef;
def setpicture   =  setvariable (picture)  enddef;
def setstring    =  setvariable (string)  enddef;
def settransform =  setvariable (transform)  enddef;
def setpen       =  setvariable (pen)  enddef;
def settension (suffix tn) expr tens =
  setnumeric (tn) if tens > 0:  tens  else:  default_tension  fi;
enddef;
def fixtension (suffix tn) =  if tn < .75:  tn := .75;  fi  enddef;

def newpicture (suffix pic) =  setpicture (pic) nullpicture;  enddef;
def convertpath (suffix g) expr f =  setpath (g) zconv (f);  enddef;

def setarray (text kind) (suffix name) =
  save name;  kind name[];  list (name)
enddef;
def setpairs =  setarray (pair)  enddef;
def gsetarray (text kind) (suffix name) =
  numeric name;  kind name[];  list (name)
enddef;


def setbbox (suffix ll, ur) =
  save ll, ur;  pair ll, ur;  getbbox (ll, ur)
enddef;

def setsplit (suffix s) expr ss =
  setnumeric (s) emax (1, ceiling ss);
enddef;

def setrgbcolor  =  setcolor  enddef;
def setcmykcolor =  setcolor  enddef;
def setcolor = setvariable (color) enddef;
def gsetcolor = gsetvariable (color) enddef;

setcolor(rgbblack) rgb(0,0,0);
setcolor(red)      rgb(1,0,0);
setcolor(green)    rgb(0,1,0);
setcolor(blue)     rgb(0,0,1);
setcolor(rgbwhite) rgb(1,1,1);
setcolor(cmykwhite) cmyk(0,0,0,0);
setcolor(cyan)      cmyk(1,0,0,0); % Maybe these should
setcolor(magenta)   cmyk(0,1,0,0); % be rbg for backward
setcolor(yellow)    cmyk(0,0,1,0); % compatibility?
setcolor(cmykblack) cmyk(0,0,0,1);
setcolor(grayscaleblack) gray(0);
setcolor(grayscalewhite) gray(1);

def setoutputtemplate text garbage = enddef;

vardef GBromannumeral (expr X) =
  save Y, _tmp, U; string U;

  Y.m := X div 1000;        % thousands digit
    _tmp := X - 1000Y.m;    % hundreds digits and lower
  Y.c := _tmp div 100;      % hundreds
    _tmp := _tmp - 100Y.c;  % tens and units
  Y.x := _tmp div 10;       % tens
  Y.i := _tmp - 10Y.x;      % units

  strrepeat("m", Y.m) &
  GBromandigit("c", "d", "m", Y.c) &
  GBromandigit("x", "l", "c", Y.x) &
  GBromandigit("i", "v", "x", Y.i)
enddef;

vardef GBromandigit (expr bot, mid, top, n) =
  if n > 9 :    top & strrepeat(bot, n-10)   % shouldn't happen
  elseif n > 8 :
    bot & top                                % "ix"
  elseif n > 4 :  mid & strrepeat (bot, n-5) % "v"--"viii"
  elseif n > 3 :
    bot & mid                                % "iv"
  else:         strrepeat (bot, n)           % ""--"iii" for 0--3
  fi
enddef;

vardef strrepeat (expr st, rep) =
  "" for i = 1 upto rep:  & st  endfor
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 is";
      show ztr;
      >> "vtr is";
      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;

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;

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

def beginmfpic (expr ch) =
  begingroup
    gcode := ch;
    save 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
  DoClip (active_plane);
  if clipall:  clipto (active_plane) rect (origin, (w_, h_));  fi
  if showbbox:  noclip ( safedraw rect (origin, (w_, h_)) );  fi
  scantokens extra_endchar;
  if proofing > 0:  makebox (proofrule);  fi
  chardx := w_;     % desired width of character in pixels
  shipit;
  if displaying > 0:  makebox (screenrule);  showit;  fi
  endgroup
enddef;

pair label_adjust;
label_adjust := origin;
numeric label_sep, labelpath_sep ;
label_sep :=  0;  labelpath_sep := 0;
def verbatimtex text t =  enddef;

vardef newgblabel (expr hf, vf, BL, r) (text s) (text pts) =
enddef;

vardef gblabel (expr a, b, c, d, r) (text s) (text t) =
  newgblabel (b, d, (c = 0) and (d = 0), r) (s) (t);
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;
  pathdims (xy, lbl) (ll, ur);
  readjustdims (ll, ur) (labelpath_sep)
  _r := if numeric rad:  rad
    elseif not boolean rad:  0
    elseif rad:  emin (xpart(ur-ll), ypart (ur-ll))/sqrt(2)
    else:  0
    fi;
  if _r = 0:
    f :=  rect (ll, ur);
  else:
    save p, q;
    pair p[];  path q;
    p1 := ur - _r*dir(45);
    p3 := ll + _r*dir(45);
    p2 := (xpart p3, ypart p1);
    p4 := (xpart p1, ypart p3);
    q := if _r < 0:  reverse  fi  quartercircle scaled 2_r;
    f :=
        (q shifted p1)--(q rotated 90 shifted p2)
             --(q rotated 180 shifted p3)
             --(q rotated -90 shifted p4)--cycle;

  fi
  readjustdims (ll, ur) (label_sep - labelpath_sep);
  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;
    pathdims (xy, lbl) (ll, ur);
    readjustdims (ll, ur) (labelpath_sep)
    cc := .5[ll, ur];
    (ww, hh) = ur - cc;
    if (ww = 0) or (hh = 0):
      f = (ll--ur);
    else:
      save aa, bb;
      aa := ww ++ if aspect:  ww  else:  hh  fi  *mult;
      bb := hh ++ if aspect:  hh  else:  ww  fi  /mult;
      f := ellipse (cc, aa, bb, 0);
    fi
    readjustdims (ll, ur) (label_sep - labelpath_sep);
    invvconv (thegblabel (ref_shift(a, b, c, ll, ur), r, f))
      shifted loc
  fi
enddef;

def pathdims (expr xy, lbl) (suffix ll, ur) =
  if pair lbl:
    ll := xy;  ur := lbl;
  else:
    ll := ur := origin;
  fi
enddef;

def readjustdims (suffix ll, ur) (expr s) =
  ll := ll - s*(1,1);
  ur := ur + s*(1,1);
enddef;

newinternal reallysmall;  reallysmall := 3epsilon;
newinternal nottoosmall;  nottoosmall := eps/2 + 2epsilon;
def signof (expr X) =  if X < 0:  -  fi  enddef;
def TruncateWarn expr s =
  GBwarn s & " is too large or undefined, so it will be truncated.";
enddef;

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

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

vardef acos primary X =
  if abs X > 1:
    TruncateWarn "Argument of arccosine";
    angle (signof(X) 1, 0)
  else:
    angle (X, 1 +-+ X)
  fi
enddef;
vardef asin primary X =
  if abs X > 1:
    TruncateWarn "Argument of arcsine";
    angle (0, signof(X) 1)
  else:
    angle (1 +-+ X, X)
  fi
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 degrees (expr t) =  t*radian  enddef;
vardef radians (expr t) =  t/radian  enddef;

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

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

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

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

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

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

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

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

vardef acosh primary y =
  if y < 1:
    TruncateWarn "acosh";
    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:
    TruncateWarn "atanh";
    signof (y) infinity
  fi
enddef;

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

primarydef Z zmul W = Z zscaled W enddef;
primarydef Z zdiv W =
  Z zmul ( unitvector (conj W) / (abs W) )
enddef;

vardef Moebius (expr A) primary Z =
  save _D; pair _D;
  _D := (1, 0) + (Z zscaled (conj A));
  (Z + A)/(abs _D) rotated (- angle _D)
enddef;
vardef pshdist (expr Z,W) =  abs(Moebius(-W)(Z))  enddef;
vardef pshdist_hp (expr Z,W) =  abs(Z-W)/abs(Z-conj(W))  enddef;
vardef kelvin (expr Z) =
  save tmp_; tmp_ = abs(Z);
  if tmp_ = 0:
    (infinity, infinity)
  elseif tmp_ < reallysmall:
    infinity*unitvector Z
  else:
    (1/tmp_)*unitvector Z
  fi
enddef;

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

primarydef x**y =
  if y=2:  x*x
  elseif (x = floor x) and (abs y = floor y):
    1 for n=1 upto y:  *x endfor
  else:  takepower y of x
  fi
enddef;
let ^ = **;

transform T_stack[];
numeric 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) )  enddef;
def ecoords =  hide ( T_pop (ztr);  vtr := vectorpart ztr )  enddef;

vardef vectorpart primary T =  T shifted -(origin transformed T)
enddef;

def apply_t (text Transformer) =
  ztr := identity Transformer transformed ztr;
  vtr := vectorpart ztr;
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 transformedpath (text Transformer) expr f =  f Transformer
enddef;

def rotatedpath (expr p, th) =
  transformedpath (
    transformed vtr
    rotatedaround (p transformed vtr, th)
    transformed (inverse vtr)
  )
enddef;
def reflectedpath (expr p, q) =
  transformedpath (
    transformed vtr
    reflectedabout (p transformed vtr, q transformed vtr)
    transformed (inverse vtr)
  )
enddef;

def scaledpath (expr p, s) =
  transformedpath (shifted -p scaled s shifted p)
enddef;
def xscaledpath (expr a, s) =
  transformedpath (shifted (-a, 0) xscaled s shifted (a, 0))
enddef;
def yscaledpath (expr b, s) =
  transformedpath (shifted (0, -b) yscaled s shifted (0, b))
enddef;

def slantedpath =  xslantedpath  enddef;
def xslantedpath (expr b, s) =
  transformedpath (shifted (0, -b) slanted s shifted (0, b))
enddef;
def yslantedpath (expr a, s) =
  transformedpath (shifted (-a, 0) yslant s shifted (0, a))
enddef;

def shiftedpath (expr v) =  transformedpath (shifted v)  enddef;

def xyswappedpath =  transformedpath (xyswap)  enddef;

vardef partialpath (expr a, b) expr f =
  save flag, flo, fhi, lo, hi, n;
  boolean flag;  flag = true;
  convertpath (g) f;
  n := length f;

  flo := snapto emin(a,b);
  if flo = 0:
    lo := 0;
  elseif flo < 1:
    setuplengtharray (cum, tot, idx) g;
    flag := false;
    lo := gettime (cum, idx) (flo*tot);
  else:
    lo := n;
  fi

  fhi := snapto emax (a,b);
  if flo = fhi:
    hi := lo;
  elseif fhi < 1:
    if flag:  setuplengtharray (cum, tot, idx) g;  fi
    hi := gettime (cum, idx) (fhi*tot);
  else:
    hi := n;
  fi

  if a > b:  reverse  fi  subpath (lo, hi) of f
enddef;

vardef gsubpath (expr a, b) expr f =  subpath (a, b) of f  enddef;

def setuplengtharray (suffix cum, tot, idx) =
  save cum, tot, idx;  idx := 0;  tot := makelengtharray (cum)
enddef;

vardef pathtime@# (suffix p) =
  if @# <= 0:  0
  elseif @# >= 1:  length p
  else:
    setuplengtharray (cum, tot, idx) p;
    gettime (cum, idx) (@#*tot)
  fi
enddef;

vardef pathpoint (expr frac) (suffix p) =
  convertpath (_pp) p;  pnt[pathtime[frac] (_pp)] (p)
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 setpicture (t) u;  andto (t, v);  t endgroup
enddef;

def orto (suffix u) (expr v) =
  mono (u);  addto u also v;  cull u keeping (1, infinity);
enddef;
primarydef u picor v =
  begingroup setpicture (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 setpicture (t) u;  xorto (t, v);  t endgroup
enddef;

def subto (suffix u) (expr v) =
  mono (u);  addto u also -v;  cull u keeping (1, infinity);
enddef;
primarydef u picsub v =
  begingroup setpicture (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;

def coloraddon (expr clr) (suffix v) =
  if clr < white:
    _orto (active_plane, v);
  else:
    _subto (active_plane, v);
  fi;
enddef;

def _orto (suffix u, v) =
  mono (u);  mono (v);
  addto u also v;
cull u keeping (1, 2);
enddef;

def _subto (suffix u, v) =
  mono (u);  mono (v);  addto u also -v;  cull u keeping (1, 1);
enddef;

vardef interior expr c =
  newpicture (v);
  addto v contour (c.t_);
  cull v dropping (0,0);
  v
enddef;

vardef interiors suffix cc =
  newpicture (_ints);
  for _idx = 1 upto cc:
    addto _ints also interior cc[_idx]);
  endfor
  mono (_ints);
  _ints
enddef;

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

vardef Clipped (suffix vt) expr c =
  setpicture (_Cl) vt;  clipto (_Cl) c;  _Cl
enddef;
def clip =  Clipped  enddef;

vardef picneg (suffix vt) expr c =
  setpicture (_pn) interior c;
  _subto (_pn, vt);
  _pn
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 =
  newpicture (v);
  if penwd >= minpenwd:
    shpath (v, drawpen) (d);
    mono (v);
  fi
  v
enddef;

def picdot (suffix v) (expr w, p) =
  addto v also
    (w shifted goodpair (p));
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;
  setnumeric (_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_ );
  newpicture (_shp);
  addto _shp contour rect (origin, (_dotwd, _dotht));
  picdot (_shp, _shp, dims);
  dims := 2dims;  mono (_shp);
  if grparam >= .5:  _shp
  else:  (interior (rect (origin, dims))) picsub _shp
  fi
enddef;

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

vardef fillwith (suffix v) (expr pic, dims, ll, ur) =
  newpicture (b);
  save fwdims, _ll, _ur;  pair fwdims, _ll, _ur;
  fwdims := goodpair (dims);
  _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);
enddef;

def thatchf (suffix v) (expr CT, sp, a, b) =
  begingroup
    setnumeric (_sp) signof (ypart b - ypart 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);
  endgroup
enddef;

def axialgradientf (suffix clr, v) (expr theta, sp, a, b) =
  begingroup
    save _hh, _sp, _nn, _y;
    _hh := ypart b - ypart a;
    _sp := signof (_hh) abs(sp);
    _nn := emax (1, round (_hh/_sp));
    _sp := _hh/_nn + signof (_hh) epsilon;
    _nn := _nn-1;
    setpath (_p) rect ((xpart a, 0),(xpart b, _sp));
    _y := ypart a;
    for _i = 0 upto _nn:
      if (clr(_i/_nn)) < white :
        addto v also shaded (clr(_i/_nn)) ( _p shifted (0,_y))
          rotated theta;
      fi
      _y := _y + _sp;
    endfor
    mono (v);
  endgroup
enddef;

def areagradientf (suffix clr, v) (expr sp, tp, a, b) =
begingroup
  save _ww, _hh, _sp, _tp, _nn, _mm, _x, _y;
  _ww := xpart b - xpart a;
  _hh := ypart b - ypart a;
  _sp := signof (_ww) abs(sp);
  _tp := signof (_hh) abs(tp);
  _nn := emax (1, round (_ww/_sp));
  _mm := emax (1, round (_hh/_tp));
  _sp := _ww/_nn + signof (_ww) epsilon;
  _tp := _hh/_mm + signof (_hh) epsilon;
  _mm := _mm-1;  _nn := _nn-1;
  setpath (_p) rect (origin,(_sp,_tp));
  _x := xpart a; y_a := ypart a;
  for _i = 0 upto _nn:
    _y := y_a;
    for _j = 0 upto _mm:
      if (clr(_i/_nn,_j/_mm)) < white:
        addto v also shaded (clr(_i/_nn,_j/_mm)) (_p shifted (_x,_y));
      fi
      _y := _y + _tp;
    endfor
    _x := _x + _sp;
  endfor
    mono (v);
endgroup
enddef;

path unitcircle;
unitcircle := fullcircle scaled 2;
def radialgradientf (suffix clr, v) (expr sp, ctr, rad) =
  begingroup
    save _sp, _r, _nn;
    _nn := emax (1, round (rad/sp));
    _sp := rad/_nn + epsilon;
    _nn := _nn - 1;
    _r := _sp;
    % fill the small center circle first
    if (clr(0)) < white :
      addto v also shaded (clr(0)) (unitcircle scaled _r shifted ctr);
    fi
    for _i = 1 upto _nn:
      if (clr(_i/_nn)) < white :
        addto v also shaded (clr(_i/_nn))
          (unitcircle scaled (_r + _sp) -- reverse unitcircle scaled _r
            --cycle) shifted ctr;
      fi
      _r := _r + _sp;
    endfor
    mono (v);
  endgroup
enddef;

def tile (suffix atile) (expr unit, width, height, clipit) =
  picture atile.pic;  atile.pic := nullpicture;
  pair atile.dims;
  atile.dims := round ((width, height)*unit);
  begingroup
    save active_plane;
    def active_plane =  atile.pic  enddef;
    save ztr, vtr;  transform ztr, vtr;
    ztr := identity scaled unit;  vtr := ztr;
    save xneg, xpos, yneg, ypos;
    xneg := 0;  xpos := width;  yneg := 0;  ypos := height;
    save ClipOn;  boolean ClipOn;
    if clipit:
      ClipOn := true;
      setarray (path) (ClipPath) (rect(origin, atile.dims));
    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.dims) and (pair  atile.dims )
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;

numeric bbox_split;  bbox_split := 4;
def getbbox (suffix ll, ur) expr g =
  setsplit (_s) bbox_split;
  ur := ll := pnt 0 (g);
  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 (ll, ur) subpath ((_j-1)/_s, _j/_s) of g;
  endfor
  if showbbox:  noclip ( safedraw rect (ll, ur) );  fi
enddef;

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

def getradius (suffix rad) expr g =
  setsplit (_s) bbox_split;
  rad := abs (pnt0 (g));
  for _j = 1 upto length g:
    rad := emax(rad, abs(pnt[_j] (g)));
  endfor
  for _j = 1 upto _s*(length g):
    ctrlsradius (rad) subpath ((_j-1)/_s, _j/_s) of g;
  endfor
enddef;

def ctrlsradius (suffix rad) expr p =
  rad := emax( emax (rad, abs(post0 (p))), abs(pre1 (p) ))
enddef;

def safedraw =  colorsafedraw (drawcolor)  enddef;
def colorsafedraw (expr clr) expr d =
  begingroup
    setpicture (v) picpath d;
    DoClip (v);  coloraddon (clr, v);
  endgroup
enddef;

def NoCycle (expr s) expr p =
  GBwarn  s & " cannot be applied to an open path."
            & " The path will be drawn instead.";
  safedraw p;
enddef;

vardef isgray (expr X) =  (X > black) and (X < white)  enddef;

def safefill =  colorsafefill (fillcolor)  enddef;
vardef colorsafefill (expr clr) expr c =
  if cycle c:
    setpicture (v) interior c;
    DoClip (v);
    if isgray (clr):
      _subto (active_plane) (v);
      v := nullpicture;
      v := shaded (clr) c;
    fi
    coloraddon (clr, v);
  else:  NoCycle("fill") c;
  fi
enddef;
def safeunfill expr c =
  if cycle c:  noclip (colorsafefill (background) c);
  else:  NoCycle("unfill") c;
  fi
enddef;

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

def store (suffix fs) expr f =
hide (
  if (not path f) and (not pair f):
    GBerrmsg ("Improper expression type.")
      "The 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 zigzag =  colorzigzag (drawcolor)  enddef;
def colorzigzag (expr clr) =  colorwiggle (false, clr, 0)  enddef;
def sinewave =  colorsinewave (drawcolor)  enddef;
def colorsinewave =  colorwiggle (true)  enddef;

vardef colorwiggle (expr smth, clr, tens, blen, elen, len, wid) expr f
=
  convertpath (g) f;
  setuplengtharray (cumlen, totlen, ct) g;
  save B;
  if cycle f:
    B := 0;
  else:
    B := abs(blen)/_rescale_factor;
    totlen := totlen - B - abs(elen)/_rescale_factor;
  fi
  setnumeric (n) 2*round (totlen/len*_rescale_factor);
  if n < 2:
    colorsafedraw (clr) g;
  else:
    save T, U, X, Y, Z, p;
    pair U, X, Y, Z;  path p;
    T := if cycle f:  0  else:  gettime (cumlen, ct) (B) fi;
    Z := pnt[T] (g);
    p :=if not cycle f:  (subpath (0,T) of g)
          if smth:  {curl 0} ..tension tens..  else:  --  fi
        fi
        for i = 1 upto n:
          hide(
            T := gettime (cumlen, ct) (B+(i/n)*totlen);
            X := Z;  Z := pnt[T] (g);
            Y := .5[X,Z];  U := sgn (Z-X);
          )
          (Y + (U zscaled (0, if even i:  -  fi  wid)))
            if smth:  {U}..tension tens..  else:  --  fi
        endfor
        if cycle f:  cycle
        else:  if smth:  {curl 0}  fi  (subpath (T, length g) of g)
        fi;
    newpicture (v);
    if smth:
      save n, k;
      n := length p;  k = n div 50;
      for i = 0 step 50 until 50*(k-1):
        shpath (v, drawpen) (subpath (i,i+50) of p);
      endfor
      shpath (v, drawpen) (subpath (50k,n) of p);
    else:  shpath (v, drawpen) (p);
    fi
    DoClip(v);  coloraddon (clr, v);
  fi
  f
enddef;

def corkscrew =  colorcorkscrew (drawcolor)  enddef;
vardef colorcorkscrew (expr clr, tens, blen, elen, len, wid) expr f =
  convertpath (g) f;
  setuplengtharray (cumlen, totlen, ct) g;
  save B;
  if cycle f:
    B := 0;
  else:
    B := abs(blen)/_rescale_factor;
    totlen := totlen - B - abs(elen)/_rescale_factor;
  fi
  setnumeric (n) round (totlen/len*_rescale_factor);
  if n < 2:
    colorsafedraw (clr) g;
  else:
    save T, U, X, Y, Z, p;
    pair U, X, Y, Z;  path p;
    T := if cycle f:  0  else:  gettime (cumlen, ct) (B) fi;
    Z := pnt[T] (g);
    p :=if (not cycle f) and (B > 0):  (subpath (0,T) of g)--  fi
        for i = 1 upto n:
          hide(
            T := gettime (cumlen, ct) (B+(i/n)*totlen);
            X := Z;  Z := pnt[T] (g);
            Y := .5[X,Z];  U := sgn (Z-X);
          )
          (X + (U zscaled (0,-wid))){ U}..tension tens..
          (Y + (U zscaled (0, wid))){-U}..tension tens..
        endfor
        if cycle f:  cycle
        else:
          {U}(Z + (U zscaled (0,-wid)))
          if elen <> 0:  --(subpath(T, length g) of g)  fi
        fi;
    newpicture (v);
    save n, k;
    n := length p;  k = n div 50;
    for i = 0 step 50 until 50*(k-1):
      shpath (v, drawpen) (subpath (i,i+50) of p);
    endfor
    shpath (v, drawpen) (subpath (50k,n) of p);
    DoClip(v);  coloraddon (clr, v);
  fi
  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;

vardef shade (expr sp) expr f =
  convertpath (g) f;
  setnumeric (gr) 1 - (.88*abs(shadewd)/sp)**2;
  if not cycle g:  NoCycle("shade") g;
  elseif gr <= 0:  safefill g;
  else:
    setbbox (ll, ur) g;
    ll := floorpair (ll);
    % setpair (dv) ceiling (sp/(sqrt 2))*(1,1);
    % test hex spacing:
    setpair (dv) ( ceiling(.5sp), ceiling(.5sp*sqrt 3) );
    setpicture (sh) setdot (shadedotpath, abs(shadewd));
    newpicture (v);
    fillwith (v) (sh, 2dv, ll, ur);
    newpicture (w);
    addto w also v shifted goodpair (dv);
    DoClip (v);  DoClip (w);
    clipto (v) (g);  clipto (w) (g);
    _orto (active_plane, v);
    v := nullpicture;
    _orto (active_plane, w);
  fi
  f
enddef;

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

vardef polkadot (expr sp) expr f =
  convertpath (g) f;
  if not cycle g:  NoCycle("polkadot") g;
  elseif sp <= emax (2*polkadotwd/3, mindotspace):
    safefill g;
  else:
    setbbox (ll, ur) g;
    save dx, dy, dshift;  pair dshift;
    dx := sp/2;  dy := dx*sqrt 3;
    dshift := (xpart(ur - ll) mod dx, ypart (ur - ll) mod dy)/2;
    save p, dims;  pair p, dims;
    p := ll + dshift;
    dims :=  2(dx, dy);
    setpicture (thepolkadot) setdot (polkadotpath, polkadotwd);
    newpicture (v);
    fillwith (v) (thepolkadot, dims, p, ur);
    fillwith (v) (thepolkadot, dims, p + (dx, dy), ur);
    DoClip (v);  clipto (v) g;
    if isgray (fillcolor):
      _subto (active_plane) (v);
      v := nullpicture;
      thepolkadot :=
        shaded (fillcolor) polkadotpath scaled ceiling (polkadotwd);
      fillwith (v) (thepolkadot, dims, p, ur);
      fillwith (v) (thepolkadot, dims, p + (dx, dy), ur);
      DoClip (v);  clipto (v) g;
    fi
    coloraddon (fillcolor, v);
  fi
  f
enddef;

def thatch =  colorthatch (hatchcolor)  enddef;
vardef colorthatch (expr clr) (expr sp, theta) expr f =
  convertpath (g) f;
  if not cycle g:  NoCycle("hatch") g;
  elseif sp <= abs(hatchwd):  colorsafefill (clr) g;
  else:
    newpicture (v);
    setbbox (ll, ur) g rotated -theta;
    thatchf (v, identity rotated theta, sp, ll, ur);
    DoClip (v);  clipto (v) (g);
    coloraddon (clr, v);
  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;
def colorxhatch (expr clr, sp) =
  colorthatch (clr) (sp, 45) colorthatch (clr) (sp, -45)
enddef;

vardef axialgradient (suffix clr) (expr sp, theta) expr f =
  convertpath (g) f;
  if not cycle g:  NoCycle("axialgradient") g;
  else:
    newpicture (_grd);
    setbbox (ll, ur) g rotated -theta;
    axialgradientf (clr, _grd) (theta, sp, ll, ur);
    DoClip (_grd);  clipto (_grd) (g);
    safeunfill g;
    _orto (active_plane, _grd);
  fi
  f
enddef;

vardef areagradient (suffix clr) (expr sp, tp) expr f =
  convertpath (g) f;
  if not cycle g:  NoCycle("areagradient") g;
  else:
    newpicture (_agr);
    setbbox (ll, ur) g;
    areagradientf (clr, _agr) (sp, tp, ll, ur);
    DoClip (_agr);  clipto (_agr) (g);
    safeunfill g;
    _orto (active_plane, _agr);
  fi
  f
enddef;

vardef radialgradient (suffix clr) (expr sp, ctr) expr f =
  convertpath (g) f;
  if not cycle g:  NoCycle("radialgradient") g;
  else:
    setpair (_ctr) zconv (ctr);
    newpicture (_agr);
    save _rad;
    getradius (_rad) g shifted - _ctr;
    radialgradientf (clr, _agr) (sp, _ctr, _rad);
    DoClip (_agr);  clipto (_agr) (g);
    safeunfill g;
    _orto (active_plane, _agr);
  fi
  f
enddef;

vardef NoTile (suffix atile) expr g =
  GBwarn str atile & " is not a valid tile for tess()."
                   & " The path will be drawn instead.";
  safedraw g;
enddef;
vardef tess (suffix atile) expr c =
  convertpath (_g) c;
  if not cycle _g:  NoCycle("tess") _g;
  elseif not is_tile (atile):  NoTile (atile) _g;
  else:
    setbbox (_ll, _ur) _g;
    newpicture (_ts);
    fillwith (_ts) (atile.pic, atile.dims, _ll, _ur);
    DoClip (_ts);  clipto (_ts) _g;
    _orto (active_plane, _ts);
  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 _rescale_factor:  _rescale_factor := 0.1in;  fi

numeric last_dot_size;  last_dot_size := 0;
vardef gendashed (suffix pat) expr f =
  convertpath (_g) f;
  save _dpat;
  if not mkdasharrays (pat) (_dpat):
    GBwarn "Dash pattern " & str pat
      & " undefined. Path will be drawn instead.";
    safedraw _g;
  elseif _dpat.rep < 2:
    safedraw _g;
  else:
    save _dl;
    forsuffixes _s = start, rep, finish:
      _dl._s := 0;
      for i = 1 upto _dpat._s:
        _dpat._s[i] := _dpat._s[i]/_rescale_factor;
        _dl._s := _dl._s + _dpat._s[i];
      endfor
    endfor

    if _dl.rep = 0:
      GBwarn "Dash pattern " & str pat & " has length 0. "
        & "Path will be drawn instead.";
      safedraw _g;
    else:
      setuplengtharray (_cumlen, _totlen, _ct) _g;
      save _n, _sf, _no_dots;
      boolean _no_dots;  _no_dots := true;
      _sf := scale_adjust (_n, _dl) (_totlen);
      if _n < 0:  safedraw _g;
      else:
        forsuffixes _s = start, rep, finish:
          for _i = 1 upto _dpat._s:
            if (_dpat._s[_i] = 0) and _no_dots:  _no_dots := false;
            else:  _dpat._s[_i] := _dpat._s[_i]*_sf;
            fi
          endfor
          _dl._s := _dl._s*_sf;
        endfor
        if _no_dots:
        else:
          if unknown plot_pic:
            save plot_pic;  path plot_pic;
            plot_pic := dotpath;
          fi;
          last_dot_size :=
            if known plot_pic.size:  plot_pic.size  else:  penwd fi;
          setpicture (dashingdot) makesymbol (plot_pic, last_dot_size);
        fi
        save _t, _d, _v;
        picture _v;  _v := nullpicture;
        _d0 := 0;  _t0 := 0;
        dashit (_dpat.start) (_v);
        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 (_dpat.rep) (_v);
            endfor
            DoClip (_v);
            coloraddon (drawcolor, _v);
            _v := nullpicture;
          endfor
        fi

        _d0 := _totlen - _dl.finish;
        _t0 := gettime (_cumlen, _ct) (_d0);
        dashit (_dpat.finish) (_v);
        DoClip (_v);
        coloraddon (drawcolor, _v);
      fi
    fi
  fi
  f
enddef;

vardef makelengtharray (suffix clen) suffix p =
  setsplit (_s) segment_split;
  numeric clen[];
  clen := _s * length p;  clen0 := 0;
  for _i = 1 upto clen:
    clen[_i] := clen[_i-1] + abs (pnt[_i/_s] (p) - pnt[(_i-1)/_s] (p))
      / _rescale_factor;
  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) =
  setnumeric (_gtl) emax (arr[ct], emin (arr[arr], lngth));
  setsplit (_s) segment_split;
  forever:  exitif inrange (arr[ct], arr[ct+1]) (_gtl);
    next ct;
  endfor
  if arr[ct] = arr[ct+1]:  ct
  else:  ( ct + (_gtl - arr[ct]) / (arr[ct+1] - arr[ct]) )
  fi  /_s
enddef;

def next suffix X =  X :=  X + 1;  enddef;

def dashit (suffix pos) (suffix pic) =
  for _k = 1 upto pos:
    if odd _k:
      if pos[_k] = 0:
        _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:
      _d0 := _d1 + pos[_k];
      _t0 := gettime (_cumlen, _ct) (_d0);
    fi
  endfor
enddef;

def dashpat (suffix pat) (text t) =
  list (pat) (t);
  if (pat = 0) or (odd (pat) and (pat > 1)):
    pat[incr pat] := 0;
  fi
enddef;

vardef mkdasharrays (suffix src, dest) =
  save _bad;  boolean _bad;  _bad := false;
  forsuffixes _s = start, rep, finish:
    numeric dest._s, dest._s[];
    boolean _bad._s;
    if knownnumericarray src._s:
      copyarray (src._s) (dest._s);
      _bad._s := false;
    else:
      _bad := _bad._s := true;
    fi
  endfor
  % _bad = one of the three arrays not copied.
  if _bad:
    if knownnumericarray src:  _bad := false;
      if _bad.rep:      % make dest.rep = src
        copyarray (src) (dest.rep);
      fi
      if _bad.start:    % shrink first dash to get dest.start
        copyarray (src) (dest.start);
        dest.start1 := dash_start*src1;
      fi
      if _bad.finish:   % use partial first dash for dest.finish
        dest.finish := 1;
        dest.finish1 := dash_finish*src1;
      fi
    fi
  fi
  not _bad
enddef;

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

def DASHED =  Dashed  enddef;
 def dashed =  Dashed  enddef;

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

path dotpath;  dotpath := fullcircle;
def dotted =  doplot (dotpath)  enddef;

vardef plotnodes (expr symbol, size) expr f =
  if size > 0:
    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);
  fi
  f
enddef;

def showcontrols =  colorshowcontrols (pointcolor)  enddef;
vardef colorshowcontrols (expr clr, syma, symb, size) expr f =
  save shpre, shpost;
  pair shpre[], shpost[];
  shpre := 0;  shpost := 0;
  for a = 0 upto (length f) if cycle f:  - 1 fi:
    shpre [incr shpre]  := pre [a] (f);
    shpost[incr shpost] := post[a] (f);
    colorsafedraw (clr)
      (zconv (shpre[shpre]--pnt[a](f)--shpost[shpost]));
  endfor
  if size > 0:
    if not numeric syma:  dosymbols (clr, syma, size) (shpre) ;  fi
    if not numeric symb:  dosymbols (clr, symb, size) (shpost);  fi
  fi
  f
enddef;

def doubledraw =  colordoubledraw (drawcolor)  enddef;
vardef colordoubledraw (expr clr, sep) expr f =
  convertpath (g) f;
  colorsafedraw (clr) (parapath ( sep/2) g);
  colorsafedraw (clr) (parapath (-sep/2) g);
  f
enddef;

vardef makesymbol (expr spath, sc) =
  if picture spath :
    setpicture (v) spath;  mono (v);  v
  elseif path spath:  setdot (spath, sc)
  else:
    GBwarn "Undefined symbol for plotting, "
         & "dotpath will be used instead.";
    setdot (dotpath, sc)
  fi
enddef;

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;
      setboolean (clearsymbols) true;
      plotsymbol (Circle, ptwd) (t);
    endgroup
  fi
enddef;

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

def clearopenpath expr f =
  if clearable (f):  safeunfill f & cycle;  fi
enddef;

def plotsymbol =  colorplotsymbol (pointcolor)  enddef;
def colorplotsymbol (expr clr, spath, sc) (text t) =
  if sc > 0:
  begingroup
    setpairs (_cpls) (t);
    if _cpls > 0:  dosymbols (clr, spath, sc) (_cpls);  fi
  endgroup
  fi
enddef;

def dosymbols (expr clr, spath, sc) (suffix arr) =
  if clearsymbols and  clearable (spath):
    addsymbols (background, makesymbol (spath&cycle, sc)) (arr);
  fi
  addsymbols (clr, makesymbol (spath, sc)) (arr);
enddef;

def addsymbols (expr clr, symb) (suffix arr) =
  newpicture (_pls);
  for _idx = 1 upto arr:
    picdot (_pls, symb, zconv (arr[_idx]));
  endfor
  DoClip (_pls);  coloraddon (clr, _pls);
enddef;

def putimage (suffix pic) (text t) =
  newpicture (_pti);
  for _itm = t:
    addto _pti also
      (pic shifted goodpair (zconv (_itm)));
    DoClip (_pti);  addto active_plane also _pti;
    _pti := nullpicture;
  endfor
 mono active_plane
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 xlow  =  xneg + laxis  enddef;
vardef xhigh =  xpos - raxis  enddef;
vardef ylow  =  yneg + baxis  enddef;
vardef yhigh =  ypos - taxis  enddef;

vardef axisline.x =  (xlow, 0)--(xhigh, 0)  enddef;
vardef axisline.y =  (0, ylow)--(0, yhigh)  enddef;
vardef axisline.l =  axisline.y shifted (xlow, 0)  enddef;
vardef axisline.b =  axisline.x shifted (0, ylow)  enddef;
vardef axisline.r =  axisline.y shifted (xhigh, 0)  enddef;
vardef axisline.t =  axisline.x shifted (0, yhigh)  enddef;

vardef axis@# (expr len) =  headpath (len, 0, 0) axisline@#  enddef;

vardef borderrect =
  rect((xlow,ylow),(xhigh,yhigh))
enddef;

vardef between (expr A, B, X) = (A <  X) and (X <  B)  enddef;
vardef inrange (expr A, B, X) = (A <= X) and (X <= B)  enddef;

vardef inbounds (expr Z) =
  inrange (xlow, xhigh) (xpart Z) and inrange (ylow, yhigh) (ypart Z)
enddef;

tertiarydef X isbetween P =  between (xpart P, ypart P, X)  enddef;
tertiarydef X isinrange P =  inrange (xpart P, ypart P, X)  enddef;
tertiarydef P contains X  =  between (xpart P, ypart P, X)  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;  path _tic;
 _ticang := if tp < 0:  inang  else:  90 fi;
 _tp := abs(tp) - 1;
 _U := unitvector (vconv (pdir)) rotated _ticang;
 _tic := (-_U--(0,0)) shifted (_tp*_U) scaled len;
 for _a = t:
   safedraw (_tic shifted zconv (loc + _a*pdir));
 endfor
enddef;

def xmarks =  axismarks ( 90, xtick, origin, right)  enddef;
def ymarks =  axismarks (-90, ytick, origin, up)    enddef;
def lmarks =  axismarks (-90, ltick, (xlow, 0), up)    enddef;
def bmarks =  axismarks ( 90, btick, (0, ylow), right)  enddef;
def rmarks =  axismarks ( 90, rtick, (xhigh, 0), up)    enddef;
def tmarks =  axismarks (-90, ttick, (0, yhigh), right)  enddef;

path griddotpath;  griddotpath := fullcircle;
def grid =  vargrid (0.5bp)  enddef;
vardef vargrid (expr dsize, xsp, ysp) =
  save gdot, gridpic;  picture gdot, gridpic;
  gdot := setdot (griddotpath, dsize);
  gridpic := nullpicture;
  for n = ceiling ((xlow)/xsp) upto floor ((xhigh)/xsp):
  for m = ceiling ((ylow)/ysp) upto floor ((yhigh)/ysp):
    picdot (gridpic, gdot, zconv ((n*xsp, m*ysp)));
  endfor
  endfor
  coloraddon (pointcolor, gridpic);
enddef;
def vgrid =  vargrid  enddef;

def hgridlines (expr ysp) =
  for n = ceiling ((ylow)/ysp) upto floor ((yhigh)/ysp):
    safedraw zconv ((xlow, n*ysp)--(xhigh, n*ysp));
  endfor
enddef;
def vgridlines (expr xsp) =
  for n = ceiling ((xlow)/xsp) upto floor ((xhigh)/xsp):
    safedraw zconv ((n*xsp, ylow)--(n*xsp, yhigh));
  endfor
enddef;
def gridlines (expr xsp, ysp) =
  vgridlines (xsp);  hgridlines (ysp);
enddef;

def vectorfield (expr len, xsp, ysp) (text fcn) (text cond) =
  save _vf, _is_OK;
  vardef _vf (expr x,y) =  ((0,0)--(fcn)) shifted (x,y)  enddef;
  vardef _is_OK (expr x,y) =  cond  enddef;
  mkvectorfield (len, xsp, ysp) (_vf, _is_OK);
enddef;

vardef mkvectorfield (expr len, xsp, ysp) (suffix vf, isOK) =
  for n = ceiling ((xlow)/xsp) upto floor ((xhigh)/xsp):
  for m = ceiling ((ylow)/ysp) upto floor ((yhigh)/ysp):
    if isOK (n*xsp,m*ysp):  arrowdraw (len) (vf(n*xsp,m*ysp));  fi
  endfor
  endfor
enddef;

def plrvectorfield (expr len, rsp, tsp) (text fcn) (text cond) =
  save _vf, _is_OK, _A, _B, _C, _D;
  _A := xlow;  _B := xhigh;
  _C := ylow;  _D := yhigh;
  vardef _vf (expr r,t) =  ((0,0)--(fcn)) shifted (r*dir t)  enddef;
  vardef _is_OK (expr r,t) =
    save _X, _Y;  _X := r*cosd t;  _Y := r*sind t;
    (cond) and between (_A, _B) (_X) and between (_C, _D) (_Y)
  enddef;
  mkplrvectorfield (len, rsp, tsp) (_vf, _is_OK);
enddef;

vardef mkplrvectorfield (expr len, rsp, tsp) (suffix vf, isOK) =
  save rmin, rmax, tmin, tmax;
  getpolarbounds;
  if rmin = 0:
    if isOK (0,tmin):  arrowdraw (len) (vf (0,tmin));  fi
    rmin := rsp;
  fi
  for n = ceiling (rmin/rsp) upto floor (rmax/rsp):
  for m = ceiling (tmin/tsp) upto floor (tmax/tsp):
    if isOK (n*rsp,m*tsp):  arrowdraw (len) (vf (n*rsp,m*tsp));  fi
  endfor
  endfor
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 (origin, 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 plrpatch (expr rstart, rstop, rstep, tstart, tstop, tstep) =
begingroup
  newpicture (v);
  patcharcs (v) (rstart, rstop, rstep, tstart, tstop);
  coloraddon (drawcolor, v);
  v := nullpicture;
  patchrays (v) (tstart, tstop, tstep, rstart, rstop);
  coloraddon (drawcolor, v);
endgroup
enddef;

def gridarcs (expr rstep) =
  beginpolargrid;
    if rmin = 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 polargrid (expr rstep, tstep) =
  gridarcs (rstep);  gridrays (tstep);
enddef;

def polargridpoints (expr dsize, rstep, tstep) =
  beginpolargrid;
    setpicture (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 beginpolargrid =
begingroup;
  save rmax, rmin, tmax, tmin;
  getpolarbounds;
  newpicture (gridpic);
enddef;

def getpolarbounds =
    save p, r, t;
    pair p[];
    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 between (xneg, xpos) (0) and between (yneg, ypos) (0):
      tmin := 0;  tmax := 360;
    elseif (p0 = origin):  tmin :=    0;  tmax :=  90;
    elseif (p1 = origin):  tmin :=  -90;  tmax :=   0;
    elseif (p2 = origin):  tmin := -180;  tmax := -90;
    elseif (p3 = origin):  tmin :=   90;  tmax := 180;
    else:
      tmax := tmin := t0 := angle p0;
      for j = 1 upto 3:
        t := t0 + anglefromto (p0, p[j]);
        if tmax < t:  tmax := t;  fi
        if tmin > t:  tmin := t;  fi
      endfor
      if between (xneg, xpos) (0):
        rmin := emin (abs(yneg), abs(ypos));
      elseif between (yneg, ypos) (0):
        rmin := emin (abs(xneg), abs(xpos));
      else:
        rmin := min (r0, r1, r2, r3);
      fi
    fi
enddef;


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

vardef polarpatch (expr rstart, rstop, rstep, tstart, tstop, tstep) =
  plrpatch (rstart, rstop, rstep, tstart, tstop, tstep);
  safedraw zconv ( arcplr (origin, 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/Bob*_uncle) = 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;
  B := pnt[n + 1] (t);
  C := pnt[n + 2] (t);
  zz = whatever[B,C];
  zz = pnt[n](t) + 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 anglefromto (expr u, v) =
  if (u = origin) or (v = origin):  0
  else:  angle (v rotated (-angle u))
  fi
enddef;

vardef cornerangle (expr A, B, C) =
  if (A = B) or (A = C) :
    if (B = C) :  60
    else:         90
    fi
  else:  anglefromto (B - A, C - A)
  fi
enddef;

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

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

vardef polyline (expr cyclic) (text t) =
  setpairs (_pl) (t);
  if _pl=0:  NoPoints ("polyline", _pl);  fi
  mkpoly (cyclic, _pl)
enddef;

def NoPoints (expr s) (suffix pts) =
  GBwarn  s & " attempted with empty list.";  pts[incr pts] := origin;
enddef;

vardef turtle (text t) =
  setnumeric (_tu) 0;
  setpair (_tmp) origin;
  pair _tu[];
  for _a = t:
    _tmp := _tmp + _a;
    _tu[incr _tu] := _tmp;
  endfor
  if _tu = 0:  NoPoints("turtle", _tu);  fi
  mkpoly (false, _tu)
enddef;

vardef brownianpath (expr start, num, sc) =
  setnumeric (_brp) 1;
  setpair (_tmp) start;
  pair _brp[];  _brp1 := _tmp;
  for _idx := 1 upto num:
    _tmp := _tmp + sc/(sqrt 2)*(normaldeviate,normaldeviate);
    _brp[incr _brp] := _tmp;
  endfor
  mkpoly (false, _brp)
enddef;
vardef randomwalk (expr  start, num, dst) =
  setnumeric (_rdw) 1;
  setpair (_tmp) start;
  pair _rdw[];  _rdw1 := _tmp;
  for _idx := 1 upto num:
    _tmp := _tmp + dst*dir(uniformdeviate(360));
    _rdw[incr _rdw] := _tmp;
  endfor
  mkpoly (false, _rdw)
enddef;
vardef browniangraph (expr num, scst) =
  setnumeric (_brg) 1;
  pair _tmp, _brg[];  _tmp := _brg1 := (0,0);
  for _idx := 1 upto num:
    _tmp := _tmp + scst*(1,normaldeviate);
    _brg[incr _brg] := _tmp;
  endfor
  mkpoly (false, _brg)
enddef;

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

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

vardef mkconvex (expr tens, cyclic) (suffix pts) =
  save _B, _d, _tmp;  pair _d[];
  settension (_tn) tens;  fixtension (_tn);
  if pts < 4:  mktenser (_tn, cyclic) (pts)
  else:
    for _j = 2 upto pts - 1:
      _B[_j] := sqrt(abs((pts[_j]-pts[_j-1])xprod(pts[_j+1]-pts[_j])));
    endfor
    if cyclic:
     _B1     := sqrt(abs((pts1    -  pts[pts])xprod(pts2 - pts1)));
     _B[pts] := sqrt(abs((pts[pts]-pts[pts-1])xprod(pts1 - pts[pts])));
    else:
      _B1     := _B2;
      _B[pts] := _B[pts-1];
    fi
    for _j = 2 upto pts - 1:
      _tmp := _B[_j-1] + _B[_j+1];
      _d[_j] :=
        if _tmp = 0:  origin % signal to use curl1
        else:
          ( _B[_j+1]*(pts[_j] - pts[_j-1]) +
            _B[_j-1]*(pts[_j+1] - pts[_j]) )/_tmp
        fi;
    endfor
    if cyclic:
      _tmp := _B[pts] + _B2;
      _d1 :=
        if _tmp = 0:  origin
        else:
          (_B2*(pts1 - pts[pts]) + _B[pts]*(pts2 - pts1))/_tmp
        fi;
      _tmp := _B[pts-1] + _B1;
      _d[pts] :=
        if _tmp = 0:  origin
        else:
          ( _B1*(pts[pts] - pts[pts-1]) +
            _B[pts-1]*(pts1 - pts[pts]) )/_tmp
        fi;
    else:
      _d1 := origin;  _d[pts] := origin;
    fi
    pts1
    for _j = 1 upto pts-1:
      {if _d[_j] = origin:  curl1  else:  _d[_j] fi}
      ..tension atleast _tn..pts[_j+1]
    endfor
      {if _d[pts] = origin:  curl1  else:  _d[pts] fi}
    if cyclic:  ..tension atleast _tn..cycle  fi
  fi
enddef;

numeric default_tension;  default_tension := 1;
def curve =  tcurve (default_tension)  enddef;
vardef tcurve  (expr tens, cyclic) (text t) =
  setpairs (_tc) (t);
  if _tc=0: NoPoints("curve", _tc); fi
  mksmooth (tens, cyclic, _tc)
enddef;

def ccurve =  tccurve (default_tension)  enddef;
vardef tccurve (expr tens, cyclic) (text t) =
  setuniquepairs (_tcc) (t);
  if _tcc=0: NoPoints("ccurve", _tcc); fi
  mkconvex (tens, cyclic, _tcc)
enddef;

vardef mkbezier (expr tens, cyclic) (suffix pts) =
  settension (_tn) tens;  fixtension (_tn);
  pts1
    for _i = 2 upto pts:  ..tension _tn..pts[_i] endfor
    if cyclic:  ..tension _tn..cycle  else:  {0,0}  fi
enddef;

def bezier =  tbezier (default_tension)  enddef;
vardef tbezier (expr tens, cyclic) (text t) =
  setpairs (_tbs) (t);
  if _tbs=0:  NoPoints ("bezier", _tbs);  fi
  mkbezier (tens, cyclic) (_tbs)
enddef;

vardef mkqbezier (expr cyclic) (suffix pts) =
  pts1
  if pts=1: {0,0}
  else:
    for _i = 2 step 2 until pts - 1:
      ..controls 1/3[pts[_i], pts[_i-1]]
             and 1/3[pts[_i], pts[_i+1]].. pts[_i+1]
    endfor
    if cyclic:
      ..controls 1/3[ pts[pts], pts[pts - 1] ]
        and 1/3[ pts[pts], pts1 ]..cycle
    fi
  fi
enddef;

vardef qbezier (expr cyclic) (text t) =
  setpairs (_qbz) (t);
  if _qbz=0:  NoPoints ("qbezier", _qbz);
  else:
    if (cyclic and odd _qbz) or (not cyclic and even _qbz):
      _qbz[incr _qbz] := _qbz[_qbz-1];
    fi
    mkqbezier (cyclic) (_qbz)
  fi
enddef;

vardef mkcbezier (expr cyclic) (suffix pts) =
  pts1
  if pts=1:  {0,0}
  else:
    for _i = 1 step 3 until pts - 3:
      ..controls pts[_i+1] and pts[_i+2] .. pts[_i+3]
    endfor
    if cyclic:
      ..controls  pts[pts - 1] and pts[pts]..cycle
    fi
  fi
enddef;

vardef cbezier (expr cyclic) (text t) =
  setpairs (_cbz) (t);
  if _cbz=0:  NoPoints ("qbezier", _cbz);
  else:
    % Need 0 mod 3 for cyclic, otherwise 1 mod 3
    setnumeric (_mdt) _cbz mod 3;
    if cyclic:
      if _mdt <> 0:  _cbz[incr _cbz] := _cbz[_cbz-1];  fi
      if _mdt = 1 :  _cbz[incr _cbz] := _cbz1;  fi
    else: % need 1 more, duplicate next to last
      if _mdt = 0:
        _cbz := _cbz + 1;
        _cbz[_cbz]   := _cbz[_cbz-1];
        _cbz[_cbz-1] := _cbz[_cbz-2];
      fi
      if _mdt = 2: % need 2 more, duplicate last 2.
        _cbz := _cbz + 2; % add 2 slots
        _cbz[_cbz]   := _cbz[_cbz-2]; % fill them
        _cbz[_cbz-1] := _cbz[_cbz-2]; % with last node
        _cbz[_cbz-2] := _cbz[_cbz-3]; % orig last slot = orig previous.
      fi
    fi
    mkcbezier (cyclic) (_cbz)
  fi
enddef;

vardef fcncontrol (expr ftens, X, Y, Z) =
  Y if (xpart(Z-Y) <> 0) and (xpart(Y-X) <> 0):
      + xpart(Z-Y)/3/xpart(Z-X)*(Z - X)/ftens  fi
enddef;

vardef mkfcnpath (expr ftens) (suffix q) =
  settension (_tn) ftens;
  if _tn <= 0:  _tn := 1;  fi
  for _i = 1 upto q - 1:
    q[_i]..controls fcncontrol (_tn) (q[_i-1], q[_i],  q[_i+1])
           and fcncontrol (_tn) (q[_i+2], q[_i+1], q[_i])..
  endfor
  q[q]{0,0}
enddef;

def fcncurve =  functioncurve (default_tension)  enddef;
def tfcncurve =  functioncurve  enddef;
vardef functioncurve (expr ftens) (text t) =
  settension (_ftens) ftens;  if _ftens < 1/3:  _ftens := 1/3;  fi
  setuniquepairs (_fc) (t);
  if _fc=0: NoPoints ("functioncurve", _fc); fi
  if _fc > 1:  _fc0 :=  _fc1;  _fc[_fc+1] := _fc[_fc];  fi
  mkfcnpath (_ftens) (_fc)
enddef;

def openqbs  =  qspline (false)  enddef;
def closedqbs =  qspline (true)  enddef;

vardef mkqbs (suffix b) =
  0.5[ b1, b2]
  if b<3:  {0,0}
  else:
    for _i = 2 upto b-1:
      ..controls 1/6[ b[_i], b[_i-1] ] and 1/6[ b[_i], b[_i+1] ]..
      0.5[ b[_i], b[_i+1] ]
    endfor
  fi
enddef;

vardef qspline (expr cyclic) (text t) =
  setpairs (_qs) (t);
  if _qs=0:  NoPoints ("qspline", _qs);  fi
  if _qs=1:  _qs[incr _qs] := _qs1;  fi
  if cyclic:
    _qs[incr _qs] := _qs1;  _qs[incr _qs] := _qs2;
  fi
  mkqbs (_qs) if cyclic:  & cycle  fi
enddef;

vardef mkcbs (suffix b) =
  (b[1]+4b[2]+b[3])/6
  if b < 4:  {0,0}
  else:
    for _i = 3 upto b-1:
      ..controls 1/3[ b[_i-1], b[_i] ] and 1/3[ b[_i], b[_i-1] ]
      .. (b[_i-1] + 4b[_i] + b[_i+1])/6
    endfor
  fi
enddef;

def mkopencbs =  mkcbs  enddef;
vardef mkclosedcbs (suffix b) =
  mkcbs (b) & opencbs (b[b-2],b[b-1],b[b], b1, b2, b3) & cycle
enddef;

def opencbs  =  cspline (false)  enddef;
def closedcbs =  cspline (true)  enddef;

vardef cspline (expr cyclic) (text t) =
  setpairs (_cs) (t);
  if _cs=0:  NoPoints ("cspline", _cs);  fi
  for _idx = _cs upto 2:  _cs[incr _cs] := _cs[_idx];  endfor
  if cyclic:
    for _idx = 1 upto 3:  _cs[incr _cs] := _cs[_idx];  endfor
  fi
  mkcbs (_cs) if cyclic:  & cycle  fi
enddef;

def init_spline_eqns (suffix pts) =
  save _spl_pre, _spl_post;
  pair _spl_pre[], _spl_post[];
  for j= 2 upto pts - 1:
   _spl_post[j] + _spl_pre[j] = 2pts[j];
   _spl_pre[j+1]+2_spl_pre[j] = 2_spl_post[j]+_spl_post[j-1];
  endfor
enddef;

def closed_spline_eqns (suffix pts) =
  _spl_post1 + _spl_pre1 = 2pts1;
  _spl_post[pts] + _spl_pre[pts] = 2pts[pts];
  _spl_pre2 + 2_spl_pre1 = 2_spl_post1 + _spl_post[pts];
  _spl_pre1+2_spl_pre[pts] = 2_spl_post[pts]+_spl_post[pts-1];
enddef;

def relaxed_spline_eqns (suffix pts) =
  _spl_pre2 + pts1 = 2_spl_post1;
  pts[pts] + _spl_post[pts-1] = 2_spl_pre[pts];
enddef;

vardef mksplinepath (expr closed) (suffix pts) =
  pts1..controls _spl_post1 and
  for j = 2 upto pts if not closed:  -1 fi:
    _spl_pre[j]..pts[j]..controls _spl_post[j] and
  endfor
    if closed:  _spl_pre1..cycle  else:  _spl_pre[pts]..pts[pts]  fi
enddef;

def mkspline (expr closed) (suffix pts) =
  init_spline_eqns (pts);
  if closed:  closed_spline_eqns (pts);
  else:     relaxed_spline_eqns (pts);
  fi
  mksplinepath (closed) (pts)
enddef;

vardef dospline (expr closed) (text the_list) =
  setpairs (_sp) (the_list);
  if _sp=0:  NoPoints ("dospline", _sp);  fi
  if _sp=1:  _sp[incr _sp] := _sp1;  fi
  mkspline (closed) (_sp)
enddef;

def init_fcnspl_eqns (suffix pts) =
  save _dx, _sl;  numeric _dx[], _sl[];
  _dx1 := xpart (pts2 - pts1);
  for j = 2 upto pts - 1:
    _dx[j] := xpart (pts[j+1] - pts[j]);
    _sl[j + 1]*_dx[j] + _sl[j-1]*_dx[j-1] + 2_sl[j]*(_dx[j] + _dx[j-1])
            = 3*ypart(pts[j+1] - pts[j-1]);
  endfor
enddef;

def periodic_fcnspl_eqns (suffix pts) =
  _sl1 = _sl[pts];
  _sl2*_dx1 + 2_sl1*_dx1 + 2_sl[pts]*_dx[pts-1] + _sl[pts-1]*_dx[pts-1]
          = 3 * ypart (pts[2] - pts[pts-1]);
enddef;

def relaxed_fcnspl_eqns (suffix pts) =
  _sl2*_dx1 + 2_sl1*_dx1 = 3 * ypart(pts2 - pts1);
  _sl[pts-1]*_dx[pts-1] + 2_sl[pts]*_dx[pts-1]
        = 3 * ypart(pts[pts] - pts[pts-1]);
enddef;

vardef mkfcnsplpath (suffix pts) =
  pts1..controls (pts1 + (1, _sl1)/3*_dx1) and
    for j = 2 upto pts - 1:
      (pts[j] - (1, _sl[j])/3*_dx[j-1]) ..pts[j]..
        controls (pts[j] + (1,_sl[j])/3*_dx[j]) and
    endfor
    (pts[pts] - (1,_sl[pts])*_dx[pts-1]/3)..pts[pts]
enddef;

vardef mkfcnspline (expr periodic) (suffix pts) =
  init_fcnspl_eqns (pts);
  if periodic:  periodic_fcnspl_eqns (pts);
  else:         relaxed_fcnspl_eqns (pts);
  fi
  mkfcnsplpath (pts)
enddef;

vardef fcnspline (expr periodic) (text the_list) =
  setpairs (_fs) (the_list);
  if _fs<2:
    if _fs=0:  NoPoints ("fcnspline", _fs);  fi
    onepointpath (false, _fs1)
  else:
    mkfcnspline (periodic) (_fs)
  fi
enddef;

vardef mkarc (expr center, begpt, endpt, sweep) =
  if (sweep = 0):  begpt--endpt
  else:
    setnumeric (n) ceiling (abs(sweep)/45);
    setpair (d) (begpt - center) rotated (signof (sweep) 90);
    begpt{d}
    for j = 1 upto n-1:
      ..(begpt rotatedabout (center, j/n*sweep)){d rotated (j/n*sweep)}
    endfor ..endpt{d rotated sweep}
  fi
enddef;

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

vardef arcpps (expr begpt, endpt, sweep) =
  if (begpt = endpt) or (sweep = 0):  begpt--endpt
  else:
    setpair (cd) unitvector (endpt-begpt);
    if abs(sweep) <= 45:
      begpt{cd rotated (-sweep/2)}..endpt{cd rotated (sweep/2)}
    elseif abs(sweep) <= 90:
      save m;  pair m;
      m = begpt + whatever*( cd rotated (-sweep/4));
      m = 0.5[begpt, endpt] + whatever*(cd rotated 90);
      begpt{cd rotated (-sweep/2)}..m{cd}..endpt{cd rotated (sweep/2)}
    else:
      setnumeric (ang) 90 - ((sweep/2) mod 180);
      if abs(ang) = 90:
        GBwarn "undefined arc. A line segment will be used instead.";
        begpt--endpt
      else:
        save c;  pair c;
        c = begpt + whatever*(cd rotated ang);
        c = if abs(ang) < 30:
              (0.5)[begpt, endpt] + whatever*(cd rotated 90)
            else:
              endpt + whatever*(-cd rotated -ang)
            fi;
        mkarc (c, begpt, endpt, sweep)
      fi
    fi
  fi
enddef;

vardef arcpp (expr small, begpt, endpt, rad) =
  save full, diam, chord, ang;
  full := signof (rad) 360;
  diam  := 2rad;
  chord := abs(endpt-begpt);
  if chord < abs(diam):
    ang := if not small:  full -  fi  2*asin (chord/diam);
  else:  ang := signof (rad) 180;
  fi
  arcpps (begpt, endpt, ang)
enddef;
def arcppr (expr begpt, endpt, rad, small) =
  arcpp (small, begpt, endpt, rad)
enddef;

vardef arcplr (expr center, frtheta, totheta, rad) =
  if rad = 0:  center--center
  else:
    mkarc (center, center + rad*dir frtheta,
                   center + rad*dir totheta, totheta - frtheta)
  fi
enddef;

vardef arcalt (expr center, radius, frtheta, totheta) =
  arcplr (center, frtheta, totheta, 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) =
  arcpps (one, two, 2*cornerangle (three, one, two))
    & arcpps (two, three, 2*cornerangle (one, two, three))
    & arcpps (three, one, 2*cornerangle (two, three, one))
    & cycle
enddef;

vardef circlepps (expr one, two, sweep) =
  save ang, full;
  full := signof (sweep) 360;
  ang := sweep mod full;
  arcpps (one, two, ang) & arcpps (two, one, full - ang) & cycle
enddef;

vardef circlepp (expr small, one, two, rad) =
  arcpp (small, one, two, rad) & arcpp (not small, two, one, rad)
    & cycle
enddef;

def circleppr (expr one, two, rad, small) =
  circleppr (one, two, rad, small)
enddef;

vardef quarterellipse(expr A,B,C) =
  save T_;
  transform T_;
  (1,0) transformed T_ = A;
  (1,1) transformed T_ = B;
  (0,1) transformed T_ = C;
  quartercircle scaled 2 transformed T_
enddef;

vardef halfellipse (expr A,B,C) =
  save P_; pair P_;
  P_ = (C - A)/2;
  quarterellipse (A, B - P_, B) & quarterellipse (B, B + P_, C)
enddef;

vardef fullellipse (expr C, A, B) =
  save P_; pair P_;
  P_ := 2[A,C];
  halfellipse (A,B,P_) & halfellipse (P_,2[B,C],A) & 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);
  save a, b, c, D, E, F;
  D := abs (B-A) = a + b;
  E := abs (C-B) = b + c;
  F := abs (A-C) = a + c;
  circleppp ((a/D)[A,B], (b/E)[B,C], (c/F)[C,A])
enddef;

vardef excircle 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 a, b, c, D, E, F;
  D := abs (B-A) = a - b;
  E := abs (C-B) = b + c;
  F := abs (C-A) = a - c;
  circleppp ((a/D)[A,B], (b/E)[B,C], (c/F)[A,C])
enddef;

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

vardef pshcircle (expr disk, ctr, rad) =
  if disk:
    if rad >= 1 :
      if rad > 1:
        GBerrmsg ("Impossible radius of pseudohyperbolic circle.")
          "The radius of a pseudohyperbolic circle can be at most 1.";
      fi
      circle ((0,0),1)
    elseif abs(ctr) >= 1 :
      if abs(ctr) > 1:
        GBerrmsg ("Impossible center of pseudohyperbolic circle.")
          "The center of this pseudohyperbolic circle must be in "
          & "the unit disk.";
      fi
      onepointpath (true,ctr)
    else:
      save _r, _dnm;
      _r := abs(ctr);
      _dnm := 1 - _r*_r*rad*rad;
      circle ((1 - rad*rad)/_dnm*ctr, rad*(1 - _r*_r)/_dnm)
    fi
  else:
    if rad >= 1 :
      GBerrmsg ("Impossible radius of pseudohyperbolic circle.")
        "The radius of a pseudohyperbolic circle must be less than 1.";
      onepointpath (true,ctr)
    elseif ypart ctr <= 0:
      if ypart ctr < 0:
        GBerrmsg ("Impossible center of pseudohyperbolic circle.")
          "The center of this pseudohyperbolic circle must be in "
          & "the upper half-plane.";
      fi
      onepointpath (true,ctr)
    else:
      save _y, _dnm;
      _y := ypart ctr;
      _dnm := 1 - rad*rad;
      circle ((xpart ctr, (1 + rad*rad)/_dnm * _y), 2rad/_dnm*_y)
    fi
  fi
enddef;

vardef UHPgeodesic (expr A, B) =
  if xpart A = xpart B:
    A--B
  else:
    save ang_, C_; pair C_;
    if abs(ypart A) < abs(ypart B):
      C_ := conj B;
    else:
      C_ := conj A;
    fi
    if ypart C_ = 0: % both on x-axis
      ang_ := anglefromto(up, B - A);
    else:
      ang_ := anglefromto(A - C_, B - C_);
    fi
    arcpps(A, B, 2ang_)
  fi
enddef;

vardef UDgeodesic (expr A, B) =
  save a_, b_;
  a_ := abs(A); b_ = abs(B);
  if (a_ = 0) or (b_ = 0):
    A--B
  elseif angle A = angle B:
    A--B
  else: % note: A, B and B-A are all nonzero from this point
    save ang_;
    if a_ = 1:
      ang_ := anglefromto (if b_>1: A else: -A fi, B-A)
    elseif b_ = 1:
      ang_ := anglefromto (A-B, if a_>1: B else: -B fi)
    else:
      save C_; pair C_;
      % reflecting A
      if a_ < eps:
        C_ := unitvector A;
        ang_1 := anglefromto(a_*A - C_, a_*B - C_);
      else:
        C_ := (1/a_)*unitvector A;
        ang_1 := anglefromto(A - C_, B - C_);
      fi
      % reflecting B
      if b_ < eps:
        C_ := unitvector B;
        ang_2 := anglefromto(b_*A - C_, b_*B - C_);
      else:
        C_ := (1/b_)*unitvector B;
        ang_2 := anglefromto(A - C_, B - C_);
      fi
      ang_ := if abs(ang_1) < abs(ang_2): ang_1 else: ang_2  fi;
    fi
    arcpps(A, B, 2ang_)
  fi
enddef;

vardef barycenter expr t =
  save m;  m := length t if not cycle t:  + 1 fi;
  pnt0(t)/m  for k = 1 upto m - 1:  + pnt[k](t)/m  endfor
enddef;

vardef sector (expr center, rad, frtheta, totheta) =
  center -- arcalt (center, rad, frtheta, totheta) -- cycle
enddef;

vardef mkbrace (expr S, C, E) =
  save R_, U_, V_, Z_;
  pair U_, V_, Z_[];
  U_ := unitvector (E-S);
  V_ := U_ rotated 90;

  R_ := 0.5*(C-S) dotprod V_;
  if R_ = 0:
    S--C
  else:
    if R_ < 0 :  V_ := -V_; R_ := -R_;  fi
    V_ := R_*V_; U_ := R_*U_;
    Z_1 := S + V_ + U_;
    Z_2 := C - V_ - U_;
    Z_3 := C - V_ + U_;
    Z_4 := E + V_ - U_;
    S{V_}..{U_}Z_1--Z_2{U_}..{V_}C{-V_}..{U_}Z_3--Z_4{U_}..{-V_}E
  fi
enddef;

vardef mkfcn (expr sm, 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 (sm, tens, false, _p)
enddef;

def tfcn (expr sm) =  mkfcn (sm, default_tension)  enddef;

def parafcn (expr sm) =  tparafcn (sm, 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 sm) (expr xmin, xmax, st) (text _fx) =
  save _fp;  vardef _fp (expr _x) =  (_x, _fx(_x))  enddef;
  mkfcn (sm, default_tension) (xmin, xmax, st) (_fp)
enddef;

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

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

def belowfcn (expr sm) =  tbelowfcn (sm, default_tension)  enddef;
vardef tbelowfcn (expr sm, tn, xlo, xhi, st)(text _fx) =
  (xlo,0)--(xhi,0)--
    (reverse tfunction (sm, tn, xlo, xhi, st)(_fx))--cycle
enddef;

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

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

def btwnplrfcn (expr sm) =  tbtwnplrfcn (sm, default_tension)  enddef;
vardef tbtwnplrfcn (expr sm, tn, tlo, thi, st)(text _ft)(text _gt)=
  tplrfcn (sm, tn, tlo, thi, st) (_ft) --
  ( reverse tplrfcn (sm, tn, tlo, thi, st) (_gt) ) -- cycle
enddef;

def plrregion (expr sm) =  tplrregion (sm, default_tension)  enddef;
vardef tplrregion (expr sm, tn, tlo, thi, st) (text _ft) =
  (0,0)--tplrfcn (sm, tn, tlo, thi, st ) (_ft)--cycle
enddef;

numeric tolerancefactor;
tolerancefactor := .02;

vardef mklevelset (expr sm, tens, X, Y, t, a, b, c, d) =
  save _inside_;
  vardef _inside_ (expr U, V) =
    inside_levelset(U, V) and between(a, b)(U) and between(c, d)(V)
  enddef;
  if not _inside_ (X, Y):
    GBwarn "Invalid seed point for levelset.";
    pairmax((a,c), pairmin((X,Y), (b,d)))&cycle
  else:
    save ls, W, A, B, prev, curr, seed;
    pair ls[], prev, curr, seed;
    seed := (X,Y);
    ls := 0;  W := 0;

    save _first_, _next_, get_next;
    vardef _first_ (expr U) =  _inside_ (U, Y)  enddef;
    vardef _next_ (expr ang) =
      _inside_ (X_curr + t * cosd ang, Y_curr + t * sind ang)
    enddef;
    def get_next (expr angA, angB) =
      X_curr := xpart curr;  Y_curr := ypart curr;
      ls[incr ls] := curr + t * dir (solve _next_ (angA, angB));
      prev := curr;  curr := ls[ls];
      W := W + anglefromto (prev - seed, curr - seed);
    enddef;

    interim tolerance := t*tolerancefactor;
    ls[incr ls] := (solve _first_ (X, b), Y);
    curr := ls[ls];
    interim tolerance := radian*tolerancefactor;
    get_next (180, 0);
    for n = 3 upto max_points:
      A := angle (curr - prev);
      get_next (A + 120, A - 120);
      exitif ((abs(W) > 180) or (ls > 10))
        and (abs(ls[ls] - ls1) < 1.2t);
    endfor
    mkpath (sm, tens, true) (ls)
  fi
enddef;

numeric max_points;
max_points := 2000;

def levelset (expr s) =  tlevelset (s, default_tension)  enddef;
vardef tlevelset (expr smth, tens, seed, seg) (text cond) =
  save inside_levelset, _t;
  vardef inside_levelset (expr x, y) =  cond  enddef;
  _t := if seg <= 0:  emax (xpos-xneg, ypos-yneg)/max_points * 20
        else:  seg fi;
  mklevelset (smth, tens, xpart seed, ypart seed, _t)
    (xneg, xpos, yneg, ypos)
enddef;

def RKIV (expr sm) =  tRKIV (sm, default_tension)
enddef;
vardef tRKIV (expr sm, tens, zstart, ds, N) (text _RHS_) =
  save _trj, _ztr, _dz, _ztmp, _ctm;
  pair _trj[],          % The trajectory
       _ztr,            % current point
       _dz[],           % array[4] of displacements
       _ztmp;           % current point for calculating velocity
  _trj := N+1;          % ultimate size of _trj array
  _trj1 := _ztr := zstart;
  save _tt,             % current time
       _dt,             % current time step
       _th;             % current time plus half a step
  _tt := 0;
  for _idx := 2 upto _trj:
    _dt   := ds/emax(1,abs(_RHS_(_tt,_ztr)));
    _th   := _tt + .5_dt;
    _dz1  := _dt*_RHS_(_tt, _ztr);  % displacement for current point
    _ztmp := _ztr + .5_dz1;         % 1st midpoint
  % use _th instead of twice calculating (_tt + .5_dt)
    _dz2  := _dt*_RHS_(_th, _ztmp); % displacement for 1st midpoint
    _ztmp := _ztr + .5_dz2;         % 2nd midpoint
    _dz3  := _dt*_RHS_(_th, _ztmp); % displacement for 2nd midpoint
    _ztmp := _ztr + _dz3;           % temporary end point
  % get time for next loop now since we need it in the next line:
    _tt   := _tt + _dt;
    _dz4  := _dt*_RHS_(_tt, _ztmp); % displacement for end point
  % get next point
    _ztr  := _ztr + (_dz1 + 2_dz2 + 2_dz3 + _dz4)/6;
    _trj[_idx] := _ztr;
  endfor
  mkpath (sm, tens, false, _trj)
enddef;

def xyRKIV (expr sm) = txyRKIV (sm, default_tension)
enddef;
vardef txyRKIV (expr sm, tens, zstart, ds, N) (text _RHS_) =
  save _fgxy, __fgxy;
  vardef __fgxy (expr t, x, y) =  _RHS_  enddef;
  vardef _fgxy (expr t, Z) = __fgxy(t, xpart Z, ypart Z) enddef;
  tRKIV (sm, tens, zstart, ds, N) (_fgxy)
enddef;

def odeRKIV (expr sm) = todeRKIV (sm, default_tension)
enddef;
vardef todeRKIV (expr sm, tens, xstart, ystart, ds, N)
(text _fxy) =
  txyRKIV (sm, tens, (xstart, ystart), ds, N) ((1, _fxy))
enddef;

vardef lclosed expr f =
  f
  if not cycle f:
    if pnt0(f) = pnt[infinity](f):  &  else:  --  fi  cycle
  fi
enddef;

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:  pnt0(f)..tension t..pnt1(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:
    if pnt0(f) = pnt[infinity](f):  &  else:  ..tension t..  fi  cycle
  fi
enddef;

def uclosed =  bclosed  enddef;
def uclosedt =  bclosedt  enddef;

def cbcontrols (suffix b, t) =
  b1 := 2[t3, t2];
  b2 := 2[t2, t1];
  b3 := 2[b1, b2];
  b4 := 2[b2, b3];
enddef;

vardef cbclosed expr f =
  save n;  n := length f;
  if cycle f:  f
  elseif n = 0:  f&cycle
  else:
    save p, q, t;  pair p[], q[], t[];
    t1 := pnt0(f);  t2 := post0(f);  t3 := pre1(f);
    cbcontrols (p, t);       % defines p1 to p4
    t1 := pnt[n](f);  t2 := pre[n](f);  t3 := post[n-1](f);
    cbcontrols (q, t);       % 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;

vardef arccomplement expr p =
  if cycle p:  onepointpath (false, pnt0(p))
  else:
    setnumeric (nn) length p;
    setpairs (pp) (pnt0(p), pnt[.5nn](p), pnt[nn](p));
    arcpps (pp3,pp1,2*cornerangle(pp2,pp3,pp1))
  fi
enddef;

path cuttings;
vardef cutoffbefore (expr b) expr f =
  save t, n;  n := length f;
  if n > 0:
    for k = 1 upto n:
      exitif (subpath (0,k) of f) intersects b;
    endfor
    if _Xtime < 0:
      cuttings := pnt0 (f){0,0};
      f
    else:
      cuttings := subpath (0,_Xtime) of f;
      subpath (_Xtime, n) of f
    fi
  else:  f
  fi
enddef;

vardef cutoffafter (expr b) expr f =
  setpath (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);
  setpair (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;  exitif v <> origin;
    v := pre [ n ] (p) - w;  exitif v <> origin;
    v := pnt [ n ] (p) - w;  exitif v <> origin;
  endfor
  sgn v
enddef;

vardef trivial expr p =  (__dir (p) = origin)  enddef;

newinternal hdwdr, hdten;
boolean hfilled;

def headshape (expr wr, tens, fil) =
  interim hdwdr := wr;
  interim hdten := if tens>0:  tens  else:  default_tension fi;
  if hdten < .75:  hdten := .75;  fi
  setboolean (hfilled) fil;
  mkheadpaths;
enddef;
def mkheadpaths =
  save Arrowhead, Leftharpoon, Rightharpoon;
  path Arrowhead, Leftharpoon, Rightharpoon,
       Arrowhead.clear, Leftharpoon.clear, Rightharpoon.clear;
  Rightharpoon := (0,0){down}..tension hdten..(.5hdwdr,-1);
  Rightharpoon.clear := Rightharpoon--(.5hdwdr,0)--cycle;
  Leftharpoon  := (reverse Rightharpoon) xscaled -1;
  Leftharpoon.clear := (reverse Rightharpoon.clear) xscaled -1;
  Arrowhead := Leftharpoon & Rightharpoon;
  Arrowhead.clear := Leftharpoon.clear & Rightharpoon.clear & cycle;
  if hfilled:
    Arrowhead := Arrowhead--cycle;
    Rightharpoon := Rightharpoon--(0,-1)--cycle;
    Leftharpoon := Leftharpoon--(0,-1)--cycle;
  fi
enddef;
headshape (1,1,false);

def head =  ahead (headcolor)  enddef;

vardef ahead (expr clr, front, back, hwr, tens, filled) =
  settension (_tn) tens;  fixtension (_tn);
  if front <> back:
    setpair (side) (hwr/2) * ((front-back) rotated 90);
    setpath (f) (back + side)..tension _tn..
       {front-back}front{back-front}..tension _tn..(back - side);
    if clearhead:
      safeunfill (back - side)--(front-side)--(front+side)--
        (back+side) & f & cycle;
      colorsafedraw (background) (back - side)--(front-side)--
        (front+side)--(back+side) & f & cycle;
    fi
    if filled:
      f := f--cycle;
      colorsafefill (clr) f;
    fi
    colorsafedraw (clr) f;
  fi
enddef;

def headpath =  Gheadpath (false) (Arrowhead)  enddef;
def headpathx =  Gheadpath (true) (Arrowhead)  enddef;

def colorheadpath =  colorGheadpath (false) (Arrowhead)  enddef;
def colorheadpathx =  colorGheadpath (true) (Arrowhead)  enddef;

def Gheadpath (expr trim) (suffix ah) =
  colorGheadpath (trim) (ah) (headcolor)
enddef;
vardef colorGheadpath
(expr trim) (suffix ah) (expr clr, sc, rot, pos) expr f =
  if (sc <> 0) and (known ah) and (path ah):
    convertpath (_g) f;
    setpair (_P) predirection[length _g] (_g);
    if _P <> origin:
      _P := _P rotated rot;
      setnumeric (_ang) anglefromto (up, _P);
      _P := pnt[length _g] (_g) - pos * _P;
      setpair (_tip) if known ah.tip:  ah.tip  else:  origin fi;

      if trim:
        if known ah.clear:
          safeunfill (ah.clear shifted - _tip)
            scaled sc
            rotated _ang
            shifted _P;
        fi
        setnumeric (_ys) max(bp, penwd, last_dot_size);
        safeunfill cut_path
          xscaled ceiling sc yscaled ceiling _ys
          rotated _ang shifted _P;
      fi
      if cycle ah:  colorsafefill  else:  colorsafedraw  fi  (clr)
        (ah shifted -_tip)
          scaled sc
          rotated _ang
          shifted _P;
    fi
  fi
  f
enddef;

path cut_path;
cut_path := (.5,0)--(.5,.71)--(-.5,.71)--(-.5,0)--cycle;

def tailpath (suffix sh) =  colortailpath (sh) (headcolor)  enddef;
vardef colortailpath (suffix sh) (expr clr, sc, rot, pos) expr f =
  if (sc <> 0) and (known sh) and (path sh):
    convertpath (_g) f;
    setpair(_P) postdirection0 (_g);
    if _P <> origin:
      _P := _P rotated rot;
      if cycle sh:  colorsafefill  else:  colorsafedraw  fi  (clr)
        (sh if known sh.tip:  shifted -sh.tip fi)
          scaled sc
          rotated anglefromto (up, _P)
          shifted (pnt0 (_g) + pos * _P);
    fi
  fi
  f
enddef;

def midpath (suffix sh) =  colormidpath (sh) (headcolor)  enddef;
vardef colormidpath (suffix sh) (expr clr, sc, rot, pos) expr f =
  if (sc <> 0) and (known sh) and (path sh):
    convertpath (_g) f;
    setnumeric (_t) pathtime[pos] (_g);
    setpair (_P) postdirection[_t] (_g);
    if _P <> origin:
      _P := _P rotated rot;
      if cycle sh:  colorsafefill  else:  colorsafedraw  fi  (clr)
        sh scaled sc
           rotated anglefromto (up, _P)
           shifted (pnt[_t] (_g));
    fi
  fi
  f
enddef;

vardef signeddeviate primary X =
  (uniformdeviate 1)[-X,X]
enddef;
vardef scaledeviate (expr W, A) =
  2 ** (signeddeviate W) * dir A
enddef;
vardef polardeviate primary R =
  (uniformdeviate abs(R)) * dir uniformdeviate 360
enddef;
vardef xydeviate primary Z =
  (signeddeviate (xpart Z), signeddeviate (ypart Z))
enddef;

vardef randompair (expr maxshift) =
  if numeric maxshift:  polardeviate (maxshift)
  elseif pair maxshift:  xydeviate (maxshift)
  else:  (0,0)
  fi
enddef;

vardef detrivialized expr f =
  save g;  path p, g[];  g := 0;
  for k = 1 upto length f:
    p := subpath (k-1,k) of f;
    if not trivial p:  g[incr g] := p;  fi
  endfor
  if g = 0:  onepointpath (cycle f, pnt0(f))
  else:  g1 for k = 2 upto g:  &g[k] endfor if cycle f:  &cycle  fi
  fi
enddef;

vardef randompath (expr maxshift, weirdness) expr f =
  save g, n;  path g;
  g := detrivialized f;
  n := length g;
  if n = 0:
    f shifted randompair (maxshift)
  else:
    save X, U, V;
    pair X[], U[], V[];
    if cycle g:  n := n - 1;  fi
    for k = 0 upto n:
      X[k] := pnt[k](g);
      U[k] := X[k] - pre[k](g);
      V[k] := post[k](g) - X[k];
    endfor
    save A, B;
    for k := 0 upto n:
      X[k] := X[k] shifted randompair (maxshift);
      A := anglefromto (U[k],V[k]);
      B := signeddeviate (30weirdness);
      U[k] := X[k] - (U[k] zscaled scaledeviate (weirdness,B));
      B := B - A + A * (2 ** signeddeviate weirdness);
      V[k] := X[k] + (V[k] zscaled scaledeviate (weirdness,B));
    endfor
    X0 for k = 1 upto n:
        .. controls V[k-1] and U[k] .. X[k]
       endfor
    if cycle g:
        .. controls V[n] and U0 .. cycle
    fi
  fi
enddef;

vardef randomlines (expr maxshift) expr f =
  save g, n;  path g;
  g := detrivialized f;
  n := length g;
  if n = 0:
    f shifted randompair (maxshift)
  else:
    if cycle g:  n := n - 1;  fi
    (pnt0(g) shifted randompair (maxshift))
      for k = 1 upto n:
        -- (pnt[k](g) shifted randompair (maxshift))
      endfor
      if cycle g:
        -- cycle
      fi
  fi
enddef;

vardef interpolatedpath (expr t, P) expr Q =
  if not path Q:
    GBerrmsg ("Improper argument to interpolatedpath.")
      "The last argument to interpolatedpath must be a path.";
    if pair P: onepointpath(false, P)
    else:
      if path P:
        P
      else:
        onepointpath (false, origin)
      fi
    fi
  elseif pair P:
    interpolated_pair_path (t, cycle Q, P, Q)
  elseif not path P:
    GBerrmsg ("Improper argument to interpolatedpath.")
      "The second argument to interpolatedpath must be a pair "
      & "or a path.";
    Q
  else:
    if t=0:     Q
    elseif t=1: P
    else:
      save P_, Q_; path P_, Q_;
      P_ := detrivialized P;
      Q_ := detrivialized Q;
      if length P_ = 0:
        interpolated_pair_path (t, cycle Q, pnt0(P_), Q)
      elseif length Q_ = 0:
        interpolated_pair_path (t, cycle Q, pnt0(Q_), P)
      else:
        save G, H, n, m, k, r;
        path G[], H[];
        G := H := 0;
        n := length P_; m := length Q_;
        k := gcd(n, m);
        r := m/k;
        for I=0 upto n-1:
          for J=0 upto r-1:
            G[incr G] := subpath (I+J/r, I+(J+1)/r) of P_;
          endfor
        endfor
        r := n/k;
        for I=0 upto m-1:
          for J=0 upto r-1:
            H[incr H] := subpath (I+J/r, I+(J+1)/r) of Q_;
          endfor
        endfor
        for N = 1 upto G-1:
          force_equal_ends(G[N], G[N+1]);
          force_equal_ends(H[N], H[N+1]);
        endfor
        interpolated_segment (t, G1, H1)
          for N = 2 upto G:  & interpolated_segment (t, G[N], H[N])
          endfor if (pnt0(G1)=pnt1(G[G])) and (cycle Q):  & cycle  fi
      fi
    fi
  fi
enddef;

vardef interpolated_pair_path (expr t, cyclic, P, Q) =
  save N; N := length Q;
  if N=0: onepointpath (cyclic, (t)[pnt0(Q),P])
  else:
    (t)[pnt0(Q),P]..controls (t)[post0(Q),P] and
      for n=1 upto N - 1:
        (t)[pre[n](Q),P]..(t)[pnt[n](Q),P]..controls (t)[post[n](Q),P]
        and
      endfor
      (t)[pre[N](Q),P].. if cyclic:  cycle  else:  (t)[pnt[N](Q),P]  fi
  fi
enddef;

vardef interpolated_segment (expr t, S, T) =
  (t)[ pnt0(S), pnt0(T)]..controls
    (t)[ post0(S), post0(T)] and (t)[ pre1(S), pre1(T)]..
      (t)[ pnt1(S), pnt1(T)]
enddef;

vardef parasegment (expr d, segs, f) =
  if d = 0:  f
  else:
    save u, v, t;  pair u[], v[];
    for n = 0 upto segs:
      t := n/segs;
      u[n] := postdirection [t] (f);
      v[n] := pnt[t] (f) + (u[n] zscaled (0,d));
    endfor
    v0{u0}
    for n = 1 upto segs:  ...v[n]{u[n]} endfor
  fi
enddef;

vardef parapath (expr d) expr f =
  if d = 0:
    f
  else:
    save a, g, h, p, q, s, t, u, v, w;
    path g[], h, p[], q[];
    numeric a, s, t;
    pair u, v, w, w[];
    s := emax(3, emin(segment_split, ceiling(max_points/5/length f)));
    p := 0;
    for i = 1 upto length f:
      h := subpath (i-1, i) of f;
      if not trivial h:
        q[incr p] := h;
        p[p] := parasegment (d, s, h);
      fi
    endfor
    if p = 0:
      f
    else:
      a := if d>0: - fi 180;
      h := p1;
      for i = 1 upto p-1:
        u := predirection 1 (q[i]);
        v := postdirection 0 (q[i+1]);
        w1 := pnt 1  (q[i])  - (u zscaled (0,d));
        w2 := pnt 0 (q[i+1]) - (v zscaled (0,d));
        w3 := pnt [infinity] (h);
        w4 := pnt 0 (p[i+1]);
        g0 := arcpps(w3, w1, a);
        g1 := h & g0;
        g2 := arcpps(w2, w4, a) & p[i+1];
        if (p[i] & g0) intersects reverse g2:
          s := length g2 - _Ytime;
          t := length h - length p[i] + _Xtime;
          g1  := subpath (0, t) of g1;
          g2 := subpath (s, length g2) of g2;
          force_equal_ends (g1, g2);
          h := g1 & g2;
        else:
          h := h .. p[i+1];
        fi
      endfor

      if cycle f:
        u := predirection 1 (q[p]);
        v := postdirection 0 (q[1]);
        w1 := pnt 1  (q[p])  - (u zscaled (0,d));
        w2 := pnt 0 (q[1]) - (v zscaled (0,d));
        w3 := pnt [infinity] (h);
        w4 := pnt 0 (p[1]);
        g3 := arcpps(w3, w1, a);
        g0 := arcpps(w2, w4, a);
        g1 := g0 & h & g3;
        g2 := g0 & p[1];
        if (p[p] & g3) intersects reverse g2:
          s := length g2 - _Ytime;
          t := length g0 + length h - length p[p] + _Xtime;
          g1 := subpath (s, t) of g1;
          force_equal_ends (g1, g1);
          h := g1 & cycle;
        else:
          h := h..cycle;
        fi
      fi
      h
    fi
  fi
enddef;

vardef turnangle@# (expr f) =
  anglefromto(predirection@# (f), postdirection@#(f))
enddef;

def setdatadashes (text lst) =
  save __type;  __type := 0;
  forsuffixes _itm = lst:
    if knownnumericarray _itm :
      copyarray (_itm) (__type[__type]);
      next __type;
    else:  GBwarn "Improper dash pattern in setdatadashes.";
    fi
  endfor
  if __type > 1:
    save dashtype;  dashtype := __type;
    for _j = 0 upto dashtype - 1:
      copyarray (__type[_j]) (dashtype[_j]);
    endfor
  else:
    SetdataWarn "dashes";
  fi
enddef;
def getdashpat expr n =  dashtype[n mod dashtype]  enddef;

def SetdataWarn expr s =
  GBwarn "command setdata"& s &"() failed. Previous values retained.";
enddef;

numeric Solid, Simpledash, Simpledot, Dotdash, Dotdashdot, Dotdashdash;
dashpat (Solid) (0);
dashpat (Simple_dash) (3bp, 4bp);
dashpat (Simple_dot) (0, 4bp);
dashpat (Dot_dash) (0, 4bp, 3bp, 4bp);
dashpat (Dot_dash_dot) (0, 4bp, 3bp, 4bp, 0, 4bp);
dashpat (Dot_dash_dash) (0, 4bp, 3bp, 4bp, 3bp, 4bp);

numeric dashtype, dashtype[], dashtype[][];
def defaultdashes =
  setdatadashes (Solid, Simple_dash, Simple_dot,
    Dot_dash, Dot_dash_dot, Dot_dash_dash);
enddef;
defaultdashes;

def setdatasymbols (text lst) =
  save __type;  path __type[];
  __type := 0;
  for _itm = lst:
    if (known _itm) and (path _itm):
      __type[__type] := _itm;
      next __type;
    else:
      GBwarn "Improper symbol in setdatasymbols().";
    fi
  endfor
  if __type > 1:
    save pointtype;  pointtype := __type;
    path pointtype[];
    for _j = 0 upto pointtype - 1:
      pointtype[_j] := __type[_j];
    endfor
  else:
    SetdataWarn "symbols";
  fi
enddef;
def getsymbol expr n :=  pointtype[n mod pointtype]  enddef;

def DeclareGBSymbols (text S) =
  forsuffixes _itm = S:
    path _itm;
    path _itm.clear;
    pair _itm.tip;
  endfor
enddef;
DeclareGBSymbols(
  Triangle, Square, Circle, Diamond, Star, Plus, Cross,
  Asterisk, Crossbar, Leftbar, Rightbar, Righthook,
  Lefthook, SolidTriangle, SolidSquare, SolidCircle,
  SolidDiamond, SolidStar
);

vardef undo_cycle expr f =  subpath (0, length f) of f  enddef;

SolidTriangle := (up--(dir 210)--(dir -30)--cycle) scaled .78;
Triangle := undo_cycle SolidTriangle;

Triangle.clear := SolidTriangle.clear :=
  ((dir -30)--(cosd 30,1)--(cosd 210,1)--(dir 210)--up--cycle)
  scaled .78;

SolidSquare := (up--(-1,1)--(-1,-1)--(1,-1)--(1,1)--cycle) scaled .443;
Square := undo_cycle SolidSquare;

SolidCircle := fullcircle rotated 90;
Circle := undo_cycle SolidCircle;
Circle.clear := SolidCircle.clear :=
  halfcircle--(-.5,.5)--(.5,.5)--cycle;

SolidDiamond := (up--left--down--right--cycle)
                scaled .522 yscaled 1.44;
Diamond := undo_cycle SolidDiamond;
Diamond.clear := SolidDiamond.clear :=
  (right--(1,1)--(-1,1)--left--up--cycle) scaled .522 yscaled 1.44;

Plus := ((0,0)--up--down--(0,0)--left--right) scaled .65;
Plus.clear := (right--(1,1)--(-1,1)--(left)--cycle) scaled .65;

Cross := ((0,0)--(dir 45)--(dir -135)--(0,0)--(dir -45)--(dir 135))
         scaled .65;
Cross.clear :=
  ((0,0)--(dir -45)--dir(45)--(dir 135)--(dir -135)--cycle) scaled .65;

Asterisk := ((0,0)--up--down--(0,0)--(dir 30)--(dir -150)
    --(0,0)--(dir -30)--(dir 150)) scaled .6;
Asterisk.clear := ((0,0)--(dir -30)--(cosd 30,1)--(cosd 150,1)
                  --(dir -150)--cycle) scaled .6;

Crossbar := ((0,0)--left--right) scaled .65;
Crossbar.clear := rect (right,(-1,.5)) scaled .65;

Leftbar := ((0,0)--left);
Rightbar := ((0,0)--right);
Leftbar.clear := rect((0,0),(-1,.5));
Rightbar.clear := rect((0,0),(1,.5));

Righthook := arcpps((0,0),(1,0),180);
Lefthook  := Righthook xscaled -1;
Righthook.clear := Righthook--cycle;
Lefthook.clear := Lefthook--cycle;

vardef mkstar (expr n, m) (suffix A) =
  save ang;  ang := 360/n;
  A1 := up;  A3 := up rotated ang;
  A2 = (whatever)[A1, A1 rotated ( ang*m)];
  A2 = (whatever)[A3, A3 rotated (-ang*m)];
  for i = 4 upto 2n:
    A[i] := A[i-2] rotated ang;
  endfor
  A := 2n;
  mkpoly (true, A)
enddef;

save _A;  pair _A[];
SolidStar := mkstar (5, 2, _A) scaled .84;
Star := undo_cycle SolidStar;
Star.clear := polyline (true)
    (_A9, _A10, _A1, _A2, _A3, (xpart _A3, 1), (xpart _A9, 1))
    scaled .84;
SolidStar.clear := Star.clear;

forsuffixes S =
 Triangle, Square, Circle, Diamond, Star, Plus, Cross,
 Asterisk, Crossbar, Leftbar, Rightbar, Righthook,
 Lefthook, SolidTriangle, SolidSquare, SolidCircle,
 SolidDiamond, SolidStar :
  S.tip := point 0 of S;
endfor

vardef gcd (expr n, m) =
  save a, b, r;
  a := emax (abs(m), abs(n));
  b := emin (abs(m), abs(n));
  if b > 0:
    forever:
      r := a mod b;
      exitif r < 1;
      a := b;  b := r;
    endfor
    b
  else:
    a
  fi
enddef;

vardef lcm (expr n, m) =
  n/gcd(n, m)*m
enddef;

numeric pointtype;  path pointtype[];
def defaultsymbols =
  setdatasymbols( Circle, Cross, SolidDiamond, Square, Plus,
    Triangle, SolidCircle, Star, SolidTriangle);
enddef;
defaultsymbols;

def computepie (suffix dat) (expr sign, ang, cent, rad) (text data) =
begingroup
  save _tot, _max, _toobig;
  _max := 0;  dat := 0;
  for _val = data:
    dat[incr dat] := _val;
    _max := emax (_max, _val);
  endfor
  if dat=0:  GBwarn  "piechart attempted with empty list.";
    _toobig := 1;
  else:
    _toobig := infinity/dat;
  fi
  if _max > _toobig:
    for _idx = 1 upto dat:
      dat[_idx] := dat[_idx]/_toobig;
    endfor
  fi
  for _idx = 2 upto dat:
    dat[_idx] := dat[_idx - 1] + dat[_idx];
  endfor
  _tot := dat[dat];
  for _idx = dat downto 2:
    dat[_idx] := ang + sign*dat[_idx-1]/_tot*360;
  endfor
  dat1 := ang; dat[dat + 1] := ang + 360sign;
endgroup
enddef;

def piechart (expr sign, ang, cent, rad) (text data) =
  save _dat;
  computepie (_dat) (sign, ang, cent, rad) (data);
  mkpiewedges (_dat, cent, rad);
enddef;

def mkpiewedges (suffix dat) (expr cent, rad) =
  numeric piewedge, piedirection, pieangle, pieangle[];
  pair piecenter, piedirection[];
  path piewedge[];
  piecenter := cent;
  piedirection := pieangle := piewedge := dat;
  for _idx = 1 upto dat:
    pieangle[_idx] := dat[_idx];
    piewedge[_idx] := sector (piecenter, rad, dat[_idx], dat[_idx+1]);
    piedirection[_idx] := dir(0.5[ dat[_idx], dat[_idx+1] ]);
  endfor
enddef;

def namedpiechart (suffix nm) (expr sign, ang, cent, rad) (text data) =
  save _dat;
  computepie (_dat) (sign, ang, cent, rad) (data);
  setnumeric (nm) _dat;
  pair nm.center, nm.direction[];
  path nm.wedge[];
  nm.center := cent;
  for _idx = 1 upto _dat:
    nm.wedge[_idx] := sector (cent, rad, _dat[_idx], _dat[_idx+1]);
    nm.direction[_idx] := dir(0.5[ _dat[_idx], _dat[_idx+1] ]);
  endfor
enddef;

def barchart (expr firstbar, 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 _nn = 1 upto chartbar:
    barstart[_nn]  := firstbar + sep*(_nn-1);
    barlength[_nn] := barend[_nn];
    chartbar[_nn]  :=  rect ((barbegin[_nn], 0), ( barend[_nn], barwd))
          shifted (0, barstart[_nn]) if vert:  xyswap  fi;
  endfor
enddef;

def namedbarchart (suffix nm) (expr first, sep, r, vert) (text data) =
  save nm;
begingroup
  save _bb, _ee, _ww;
  path nm.bar[];
  nm := 0;  _ww := r*sep;
  for _itm = data:
    _ee := if pair _itm:  ypart _itm  else:  _itm  fi;
    _bb := if pair _itm:  xpart _itm  else:  0  fi;
    nm.bar[incr nm] :=  rect ((_bb, 0), ( _ee, _ww) )
          shifted (0, first + sep*(nm-1)) if vert:  xyswap fi;
  endfor
endgroup
enddef;

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

def keepit =
  addto totalpicture also currentpicture;
  mono (totalpicture);
  currentpicture := nullpicture;
  totalnull := totalnull or 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;
      addto _v_ also totalpicture;
      do _v_
    endgroup
  fi
enddef;

boolean noship;  noship := false;
def shipit =  if noship:  else:  mergeit (shipout)  fi  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'.