mfpic (4dec12)

git-svn-id: svn://tug.org/texlive/trunk@28444 c570f23f-e606-0410-a88d-b1316a301751

1 files changed, 421 insertions, 62 deletions

diff --git a/Master/texmf-dist/metafont/mfpic/grafbase.mf b/Master/texmf-dist/metafont/mfpic/grafbase.mf
index 87492e00849..25f9a9cf768 100644
--- a/Master/texmf-dist/metafont/mfpic/grafbase.mf
+++ b/Master/texmf-dist/metafont/mfpic/grafbase.mf
@@ -8,7 +8,7 @@
 %% 
 %% -------------------------------------------------------------------
 %% 
-%% Copyright 2002--2011, Daniel H. Luecking
+%% 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
@@ -29,13 +29,14 @@ fi
 boolean grafbase; grafbase := true;
 
 string fileversion, filedate;
-fileversion := "1.06";  filedate := "2011/02/25";
+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 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 "";
@@ -57,7 +58,7 @@ def checkversions (expr g)=
   fi
 enddef;
 
-checkversions (106);
+checkversions (110);
 
 if unknown base_name       :  input plain;
 elseif not string base_name:  input plain;
@@ -81,12 +82,13 @@ vardef mftitle expr t =
   if string t:  t;  message t;  fi
 enddef;
 
-boolean METAPOST;
-METAPOST := known color Geamparalele din Babadag;
+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.";
+    "This file is for Metafont. For Metapost, use grafbase.mp.";
 fi
 
 if unknown mode:
@@ -306,12 +308,13 @@ 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))
+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;
+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;
@@ -333,7 +336,8 @@ vardef maxpair (suffix p) =  setpair (_mp) p1;
   pairmax (_mp, p[p])
 enddef;
 
-primarydef Z xprod W =  (xpart Z * ypart W - xpart W * ypart Z)  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;
@@ -422,7 +426,7 @@ def endimage =
 enddef;
 
 def makeimage (suffix name) (expr refpt) =
-  setpair (_image_reference_point) zconv(refpt);
+  setpair (_image_reference_point) zconv (refpt);
   setpicture (name) beginimage
 enddef;
 def concludeimage  =
@@ -679,7 +683,8 @@ vardef textrectx (expr a, b, c, rot, xy, lbl, rad, loc) =
 
   fi
   readjustdims (ll, ur) (label_sep - labelpath_sep);
-  invvconv (thegblabel (ref_shift(a, b, c, ll, ur), rot, f)) shifted loc
+  invvconv (thegblabel (ref_shift(a, b, c, ll, ur), rot, f))
+    shifted loc
 enddef;
 
 def textovalx =  xellipse (true)  enddef;
@@ -704,7 +709,8 @@ vardef xellipse (expr aspect, a, b, c, r, xy, lbl, mult, loc) =
       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
+    invvconv (thegblabel (ref_shift(a, b, c, ll, ur), r, f))
+      shifted loc
   fi
 enddef;
 
@@ -776,16 +782,16 @@ 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 =  (acos X)/radian  enddef;
-vardef invsin primary X =  (asin X)/radian  enddef;
-vardef invtan primary X =  (atan X)/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;
-def log =  ln  enddef;
+vardef log primary X =  ln (X)  enddef;
 vardef logbase (expr B) primary X =  (mlog X)/(mlog B)  enddef;
-def logtwo =  logbase( 2)  enddef;
-def logten =  logbase(10)  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));
@@ -858,15 +864,35 @@ 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 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;
@@ -894,7 +920,8 @@ 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;
+vardef vectorpart primary T =  T shifted -(origin transformed T)
+enddef;
 
 def apply_t (text Transformer) =
   ztr := identity Transformer transformed ztr;
@@ -1208,6 +1235,80 @@ def thatchf (suffix v) (expr CT, sp, a, b) =
   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;
@@ -1241,6 +1342,7 @@ 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);
@@ -1248,17 +1350,31 @@ def getbbox (suffix ll, ur) expr g =
     ll := pairmin (ll, pnt[_j] (g));  ur := pairmax (ur, pnt[_j] (g));
   endfor
   for _j = 1 upto _s*(length g):
-    ctrlsbbox (subpath ((_j-1)/_s, _j/_s) of g) (ll, ur);
+    ctrlsbbox (ll, ur) subpath ((_j-1)/_s, _j/_s) of g;
   endfor
   if showbbox:  noclip ( safedraw rect (ll, ur) );  fi
 enddef;
 
-numeric bbox_split;  bbox_split := 2;
-def ctrlsbbox (expr p) (suffix ll, ur) =
+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
@@ -1323,7 +1439,8 @@ 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 =
+vardef colorwiggle (expr smth, clr, tens, blen, elen, len, wid) expr f
+=
   convertpath (g) f;
   setuplengtharray (cumlen, totlen, ct) g;
   save B;
@@ -1514,6 +1631,50 @@ 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.";
@@ -1624,8 +1785,8 @@ vardef makelengtharray (suffix clen) suffix p =
   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;
+    clen[_i] := clen[_i-1] + abs (pnt[_i/_s] (p) - pnt[(_i-1)/_s] (p))
+      / _rescale_factor;
   endfor
   clen[clen]
 enddef;
@@ -2201,6 +2362,25 @@ vardef brownianpath (expr start, num, sc) =
   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)
@@ -2237,8 +2417,8 @@ vardef mkconvex (expr tens, cyclic) (suffix pts) =
       _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])));
+     _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];
@@ -2282,12 +2462,16 @@ enddef;
 numeric default_tension;  default_tension := 1;
 def curve =  tcurve (default_tension)  enddef;
 vardef tcurve  (expr tens, cyclic) (text t) =
-  setpairs (_tc) (t);  mksmooth (tens, cyclic, _tc)
+  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);  mkconvex (tens, cyclic, _tcc)
+  setuniquepairs (_tcc) (t);
+  if _tcc=0: NoPoints("ccurve", _tcc); fi
+  mkconvex (tens, cyclic, _tcc)
 enddef;
 
 vardef mkbezier (expr tens, cyclic) (suffix pts) =
@@ -2309,8 +2493,8 @@ vardef mkqbezier (expr cyclic) (suffix pts) =
   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]
+      ..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] ]
@@ -2389,6 +2573,7 @@ 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;
@@ -2560,7 +2745,7 @@ enddef;
 def arccps =  arc  enddef;
 
 vardef arcpps (expr begpt, endpt, sweep) =
-  if begpt = endpt:  begpt--endpt
+  if (begpt = endpt) or (sweep = 0):  begpt--endpt
   else:
     setpair (cd) unitvector (endpt-begpt);
     if abs(sweep) <= 45:
@@ -2649,20 +2834,42 @@ vardef circlepps (expr one, two, sweep) =
 enddef;
 
 vardef circlepp (expr small, one, two, rad) =
- arcpp (small, one, two, rad) & arcpp (not small, two, one, rad) & cycle
+  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);
+    a2 = pnt [ n/4] (p);
     a4 = pnt [3n/4] (p);
   else:
     a2 := a3;
@@ -2716,43 +2923,93 @@ vardef pshcircle (expr disk, ctr, rad) =
     elseif abs(ctr) >= 1 :
       if abs(ctr) > 1:
         GBerrmsg ("Impossible center of pseudohyperbolic circle.")
-          "The center of a pseudohyperbolic circle must be in "
+          "The center of this pseudohyperbolic circle must be in "
           & "the unit disk.";
       fi
       onepointpath (true,ctr)
     else:
-      % compute Euclidean center and radius (and a denominator used twice
-      % in calculations).
       save _r, _dnm;
       _r := abs(ctr);
       _dnm := 1 - _r*_r*rad*rad;
-      circle ( (1 - rad*rad)/_dnm*ctr, rad*(1 - _r*_r)/_dnm)
+      circle ((1 - rad*rad)/_dnm*ctr, rad*(1 - _r*_r)/_dnm)
     fi
   else:
     if rad >= 1 :
-      GBerrmsg ("Impossible pseudohyperbolic circle.")
+      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 pseudohyperbolic circle.")
-          "The center of a pseudohyperbolic circle must be in "
+        GBerrmsg ("Impossible center of pseudohyperbolic circle.")
+          "The center of this pseudohyperbolic circle must be in "
           & "the upper half-plane.";
       fi
       onepointpath (true,ctr)
     else:
-      % compute Euclidean center and radius (and a denominator used twice
-      % in calculations).
-      % Euclidean center at  xpart ctr + (1 + R^2)/(1 - R^2)*ypart ctr
-      % Euclidean radius  4R/(1 - R^2)*ypart ctr
       save _y, _dnm;
       _y := ypart ctr;
       _dnm := 1 - rad*rad;
-      circle ( (xpart ctr, (1 + rad*rad)/_dnm * _y), 2rad/_dnm*_y)
+      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
@@ -2762,6 +3019,26 @@ 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;
@@ -2834,7 +3111,7 @@ 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)
+    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.";
@@ -2902,16 +3179,16 @@ vardef tRKIV (expr sm, tens, zstart, ds, N) (text _RHS_) =
   for _idx := 2 upto _trj:
     _dt   := ds/emax(1,abs(_RHS_(_tt,_ztr)));
     _th   := _tt + .5_dt;
-    _dz1  := _dt*_RHS_(_tt, _ztr);  % displacement based on current point
+    _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 based on 1st midpoint
+    _dz2  := _dt*_RHS_(_th, _ztmp); % displacement for 1st midpoint
     _ztmp := _ztr + .5_dz2;         % 2nd midpoint
-    _dz3  := _dt*_RHS_(_th, _ztmp); % displacement based on 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 right away in next line:
+  % get time for next loop now since we need it in the next line:
     _tt   := _tt + _dt;
-    _dz4  := _dt*_RHS_(_tt, _ztmp); % displacement based on end point
+    _dz4  := _dt*_RHS_(_tt, _ztmp); % displacement for end point
   % get next point
     _ztr  := _ztr + (_dz1 + 2_dz2 + 2_dz3 + _dz4)/6;
     _trj[_idx] := _ztr;
@@ -3165,7 +3442,8 @@ vardef colorGheadpath
   f
 enddef;
 
-path cut_path;  cut_path := (.5,0)--(.5,.71)--(-.5,.71)--(-.5,0)--cycle;
+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 =
@@ -3283,6 +3561,85 @@ vardef randomlines (expr maxshift) expr f =
   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:
@@ -3305,7 +3662,7 @@ vardef parapath (expr d) expr f =
     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)));
+    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;
@@ -3475,8 +3832,8 @@ 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;
+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;
@@ -3512,7 +3869,8 @@ 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;
+    (_A9, _A10, _A1, _A2, _A3, (xpart _A3, 1), (xpart _A9, 1))
+    scaled .84;
 SolidStar.clear := Star.clear;
 
 forsuffixes S =
@@ -3620,7 +3978,8 @@ def barchart (expr firstbar, sep, r, vert)(text data) =
   path chartbar[];
   chartbar := 0;  barwd := r*sep;
   for _itm = data:
-    barend[incr chartbar] := if pair _itm:  ypart _itm  else:  _itm  fi;
+    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;