#!/usr/bin/perl -w
# GPL copyleft, Dror Bar-Natan;
# provided you don't make nasty comments about the code below.

print("makefont running...");

use Getopt::Long;

$fig2mf = "fig2dev -L mf -x 0.6 -X 5.4 -y 0.728 -Y 5.82 -m 2.338";

GetOptions(
  "source=s" => \$source, "s=s" => \$source,
  "fontname=s" => \$fn, "fn=s" => \$fn,
  "f2m_opts_file=s" => \$opts_file, "f2m_opts=s" => \$opts_file
) || die "failed reading command line options";
$mf = $fn.".mf";
$sty = $fn.".sty";
print(
  "source directory is ", $source, ";\n",
  "metafont output to ", $mf, ";\n",
  "sty output to ", $sty, ";\n"
);
open(MF, ">".$mf) || die "can't open ".$mf.": $!";
open(STY, ">".$sty) || die "can't open ".$sty.": $!";

if ($opts_file) {
  open(FORMATTING, $opts_file) || die "can't open ".$opts_file.": $!";
  while(<FORMATTING>) {
    chomp;
    @line = split(/:/);
    $f2m_opts{$line[0]} = $line[1];
  }
  print("fig2mf options read from ".$opts_file.".\n");
}

print MF "message \"This is ".$fn.".mf, made on ".localtime()."\";\n";
while (<DATA>) {print MF};

print STY "% This is `".$sty."'
% GPL copyleft
\\NeedsTeXFormat{LaTeX2e}% LaTeX 2.09 can't be used (nor non-LaTeX)
[1994/12/01]% LaTeX date must be December 1994 or later
\\ProvidesPackage{".$fn."}[2000]
\\DeclareFontFamily{U}{".$fn."}{}
\\DeclareSymbolFont{".$fn."}{U}{".$fn."}{m}{n}
\\DeclareFontShape{U}{".$fn."}{m}{n}{<1-100>".$fn."}{}
";

while(defined($fig = glob($source."/[0-9][0-9][0-9]*.fig"))) {
  $num = $fig; $num =~ s#$source/##; $num = substr($num,0,3);
  $name = $fig; $name =~ s#$source/[0-9]{3}##; $name =~ s/.fig//;
  if (!exists $f2m_opts{$name}) {$f2m_opts{$name} = "";}
  open(FIG, $fig2mf." -C ".$num." ".$f2m_opts{$name}." ".$fig." |");
  while (<FIG>) {
    s/incr code/code/;
    s/input grafbase.mf; //;
    s/end\.//;
    print MF
  };
  close(FIG);
  print STY
    "\\DeclareMathSymbol{\\".$name."}{\\mathord}{"
    .$fn."}{".$num."}\n";
  print $num, " ", $name, " ", $f2m_opts{$name}, "\n";
}
print MF "end.\n";

__END__

%%%  First we include the grafbase.mf macro package, from
%%%  CTAN/graphics/mfpic. Then the font characters in order.

%%%
%%%  File: grafbase.mf
%%%

message "mfpic version 0.2.10.2 alfa, Tue 23 Jan 1996";
message "";

% nonstopmode;

% Don't complain when variables get too large.
interim warningcheck := 0;

% Determine whether Metafont or MetaPost is being used.
% MetaPost has colors.
boolean metapost;
metapost = if known blue: true else: false fi;

if metapost:
 % ensure that Plain MetaPost is loaded.
 if unknown base_name:
  input mfplain;
 elseif not string base_name:
  input mfplain;
 elseif base_name <> "mfplain":
  input mfplain;
 fi
else:  % METAFONT
 % ensure that Plain Metafont is loaded.
 if unknown base_name:
  input plain;
 elseif not string base_name:
  input plain;
 elseif base_name <> "plain":
  input plain;
 fi
fi

% Only for MetaPost:
if metapost:
% Dubious, but one of the most common printer modes.
  mode = cx;
fi

% intercept the mode variable, before mode_setup can set proof mode.
if unknown mode:
  errhelp ("Please use \mode=localfont; or a mode known on your system.");
  errmessage ("Unknown METAFONT mode");
fi

mode_setup;

% Font identifier, and font coding scheme.
% NB:  These become ALL-CAPITALS in the PL (and TFM) files.

font_identifier := "MFpic graphics";
font_coding_scheme := "Arbitrary";

% Design size: this is somewhat arbitrary, but it needs to be large.

interim designsize := 128pt#;

%% Global Variables.

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

% message "DEBUG " & if debug: "ON" else: "OFF" fi;

if debug:
  message "DEBUG:  pixels_per_inch = " & decimal pixels_per_inch;
fi

% when the numeric variable `grafbase' is known,
% then grafbase.mf has been input.
% (idea taken from DEK's cmbase.mf)

if known grafbase:
  errhelp ("Please make sure that grafbase.mf is loaded only once.");
  errmessage ("You have loaded grafbase.mf more than once");
else:
  boolean grafbase;
  grafbase := true;
fi

% METAFONT title string t, and message t.

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

% character size variables, and their default values.

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

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

% degrees for angles.  In METAFONT, one degree is the unit of angle.

deg := 1;

% Drawing Pen.

% pen width for drawing, in _pixel_ coordinates.
newinternal penwd;

% pen for drawing.
pen drawpen;

% Hatching Pen.

% pen width for hatching, in _pixel_ coordinates.
numeric hatchwd;
hatchwd := 1pt;

% pen for hatching.
pen hatchpen;
hatchpen := pencircle scaled hatchwd yscaled aspect_ratio;

% Storing a Path.

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

% Clipping.

% clipping state.
boolean ClipOn;
% clip by default: this helps Metafont and TeX.
ClipOn := true;
ClipOn := false;

% clipping paths.
path ClipPath[];
numeric ClipPath;
ClipPath = 0;  % initially no clipping paths, and no candidates known.

%% Frank Michielsen's macros, slightly modified.
% These are:  list, pathlength, timelist, dashlist, dotlist.
% Global Parameters:  segmentsplit, dashsize, dashgap, dashstart,
% dashfinish, dotgap.

def list (suffix v) (text lst) =
  begingroup
    v := 0;
    for item= lst: 
      incr v;
      v[v]:= item;
    endfor;  
  endgroup
enddef;

if unknown segmentsplit: segmentsplit:= 5; fi

def pathlength (expr p) =
  begingroup save l;
    l:= 0;
    for i= 1 upto length(p):
      for j= 1 upto segmentsplit:
        l:= l + length( (point ((i-1) + j/segmentsplit) of p) - 
                        (point ((i-1) + (j-1)/segmentsplit) of p) );
      endfor
    endfor
    l
  endgroup
enddef;

def timelist (suffix tlst) (suffix llst) (expr p) =
  begingroup save t,sl,tl,ct;
  numeric ct,t,t[],sl[],tl[]; 
    t[0]:= 0;
    for i= 1 upto length(p):
      for j= 1 upto segmentsplit:
        t[ (i-1)*segmentsplit + j ]:= i-1 + j/segmentsplit;
      endfor
    endfor;
    t:= length(p) * segmentsplit;
% sl = segment-length
% tl = total-length
    tl[0]:= 0;  
    for i= 1 upto t: 
      sl[i]:= length( point t[i] of p - point t[i-1] of p);
      tl[i]:= tl[i-1] + sl[i];
    endfor
% ct = current time  
    ct:= 0;
    for i= 1 upto llst:
% convert negative lengths to beginning of path
      if llst[i] <= 0: 
        tlst[i]:= 0; ct:= 0;
% convert too large lengths to end of path
      elseif llst[i] >= tl[t]: 
        tlst[i]:= t; 
% find the startpoint [[tc]] of the sub-segment 
% which contains the specified point and solve by 
% linear interpolation.
      elseif llst[i] < tl[ct]: 
        forever: 
          ct:= ct - 1;
          exitif llst[i] >= tl[ct];
        endfor 
        tlst[i]:= ct + (llst[i] - tl[ct])/sl[ct+1];  
      else: % llst[i] >= tl[ct]
        forever: 
          ct:= ct + 1;
          exitif tl[ct] > llst[i];
        endfor
        ct:= ct - 1;
        tlst[i]:= ct + (llst[i] - tl[ct])/sl[ct+1];
      fi
      tlst[i]:= tlst[i]/segmentsplit;
    endfor
  tlst:= llst;
  endgroup
enddef;

if unknown dashsize:     dashsize  := 3pt; fi
if unknown dashgap:      dashgap := 1dashsize; fi
if unknown dashfinish:   dashfinish := .5; fi
if unknown dashstart:    dashstart := .5; fi

def dashlist (suffix dp) (expr p) =
  begingroup save l,x,n,dsh,spc,scl,llst;
    numeric lt,n,dsh,spc,scl,llst,llst[],tlst,tlst[];

    lt:= pathlength(p);
    n:= ( lt - (dashstart + dashfinish -1)*dashsize )/(dashsize+ dashgap);

% if the path is too short, we do not dash it.
    if n < 1:
       dp[1]:= p; dp:=1; 
% otherwise we must adjust the dashsize and dashgap to make it fit
    else: 
      n:= round(n);
      scl:= lt/( (dashstart + dashfinish + n -1)*dashsize + n*dashgap );
      dsh:= scl*dashsize;
      spc:= scl*dashgap;
% and generate a list which contains the subdivision of the path 
      llst[1]:= 0;
      llst[2]:= dashstart*dsh;
      for i=2 upto n:
        llst[2i-1]:= llst[2i -2] + spc;
        llst[2i]:= llst[2i-1] + dsh;
      endfor;
      llst[2n+1]:= llst[2n] + spc;
      llst[2n+2]:= llst[2n+1] + dashfinish*dsh;
      llst:= 2n+2;
% then convert the list of pathlengths to a list of pathtimes
      timelist(tlst)(llst)(p);
% and keep only the subpaths corresponding to the dashes
      for i=1 upto n+1:
        dp[i]:= subpath(tlst[2*i-1],tlst[2*i]) of p;
      endfor  
      dp:= n+1;
    fi;  
  endgroup
enddef;

if unknown dotgap: dotgap:= 3pt; fi

def dotlist (suffix dp) (expr p) =
  begingroup save lt,sl,n,tlst,llst;
    numeric tlst,tlst[],llst,llst[];
    lt:= pathlength(p);
    n:= lt/dotgap;
    n:= round(n);
    if n = 0:
      llst[1]:= lt/2;
      llst:= 1;
    else:
      sl:= lt/n;
      llst[1]:= 0;
      for i= 2 upto n+1: llst[i]:= llst[i-1] + sl; endfor
      llst:= n+1;
    fi

    timelist(tlst)(llst)(p);
    
    for i=1 upto n+1: dp[i]:= point tlst[i] of p; endfor
    dp:= n+1
 endgroup;
enddef;

% That was Frank Michielsen's contribution.

%% Utility Macros.

% (Reference for map macro: Abelson & Sussman,
% _Structure & Interpretation of Computer Programs_,
% page 250.)

% Generate new text list by applying proc to each item in given
% text list.  The new list must be separated by commas, since I
% wish to use it as the list in a `for' loop; this requirement
% is difficult to reconcile with other needs.

% NOTE WELL:  Text arguments are particularly pernicious,
% because they can inadvertently include a local name
% used in the body of the macro.  Also, METAFONT then gives
% rather uninformative error messages.

def map (text proc) (text list) =
  for a=list:
    , proc (a)
  endfor
enddef;

% apply procedure proc to each member of array p[] with p members.

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

% convert text pairs list t to array p.
% best used inside a vardef.
% note:  t must be a comma-separated list,
% because the `for' needs it.

def textpairs (suffix p) (text t) =
 save p; pair p[]; p:=0;
 for a=t:
   p[incr p] := a;
 endfor
enddef;

% chpair : `change pair'.
% apply macro proc (which maps numeric -> numeric)
% to each part of pair p, return resultant pair.

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

vardef floorpair (expr p) =
  chpair (floor) (p)
enddef;

vardef ceilingpair (expr p) =
  chpair (ceiling) (p)
enddef;

vardef hroundpair (expr p) =
  chpair (hround) (p)
enddef;

% Return the pair comprising the minimum x and minimum y coordinates
% of all pairs in the array p[], where p itself is a numeric count of
% the members in p[].

vardef minpair (suffix p) =
 save P_; pair P_;
 save x_, y_; numeric x_, y_;
 P_ := p1;
 for i=2 upto p:
   x_ := min (xpart P_, xpart p[i]);
   y_ := min (ypart P_, ypart p[i]);
   P_ := (x_, y_);
 endfor
 P_
enddef;

% Similarly, but for the maximum coordinates.

vardef maxpair (suffix p) =
 save P_; pair P_;
 save x_, y_; numeric x_, y_;
 P_ := p1;
 for i=2 upto p:
    x_ := max (xpart P_, xpart p[i]);
    y_ := max (ypart P_, ypart p[i]);
    P_ := (x_, y_);
 endfor
 P_
enddef;


%% Coordinate Conversion.
% (affine transformation)
% graph -> sharp
% xneg, xpos, yneg, ypos are in _graph_ coordinates.
% plain MF's beginchar defines the _sharp_ values charwd and charht
% as well the _pixel_ quantities w and h.
% (Reference:  _The METAFONTbook_, page 275, definition of beginchar.)
% _graph_ and _sharp_ coordinates coincide iff ztr is the identity
% transform.

transform ztr, invztr;

% currenttransform (.t_) takes care of the aspect ratio.
% (Also, I may try coordinate transformations of currenttransform.
%  Consider bcoords ... ecoords.)

def setztr =
 if debug:
  message " DEBUG ";
  message "charwd = " & decimal charwd & "pt#";
  message "charht = " & decimal charht & "pt#";
  message "w = " & decimal w & " pixels";
  message "h = " & decimal h & " pixels";
  message "xneg = " & decimal xneg;
  message "xneg = " & decimal xpos;
  message "yneg = " & decimal yneg;
  message "yneg = " & decimal ypos;
  message "xscale = " & decimal xscale;
  message "yscale = " & decimal yscale;
  message "unitlen = " & decimal unitlen & "pt#";
  message "hppp = " & decimal hppp;
  message " GUBED ";
 fi
 ztr := identity
   shifted (-(xneg,yneg))
   xscaled (xscale*unitlen*hppp)
   yscaled (yscale*unitlen*hppp);
 invztr := inverse ztr;
 if debug:
   message " DEBUG ";
   message "ztr: ";
   show ztr;
   message " GUBED ";
 fi
enddef;

% zconv converts a variety of expressions from graph -> sharp coords.
% The expressions include pairs, paths, and transforms.
% This is an affine transform.

% Clipping greatly complicates zconv, because of the many cases.
% Another disadvantage is that a clipped point distorts paths.

boolean SimpleZconv;
SimpleZconv = true;

if SimpleZconv:

vardef zconv (expr a) =
  a transformed ztr
enddef;

else:  % Complicated zconv.

vardef zconv (expr a) =
  if ClipOn:  % Clipping on the character boundary, if I know how.
    if pair a:
      if (xpart a < xneg) or (xpart a > xpos) or
         (ypart a < yneg) or (ypart a > ypos):
        origin
      else:  % A pair that is within the outer boundary.
        a transformed ztr
      fi
    else:  % Not a pair.  (Paths should also be clipped, but how?)
      a transformed ztr
    fi
  else:  % No clipping.
    a transformed ztr
  fi
enddef;

fi  % SimpleZconv.

% invzconv converts a variety of expressions from sharp -> graph coords.
% The expressions include pairs, paths, and transforms.
% This is an affine transform.

vardef invzconv (expr a) =
 a transformed invztr
enddef;

% vconv converts a vector v from graph -> sharp coords.
% This is a linear (ie, vector) transform.

vardef vconv (expr v) =
 (v transformed ztr) - (origin transformed ztr)
enddef;

% invvconv converts a vector from sharp -> graph coords.
% This is a linear (ie, vector) transform.

vardef invvconv (expr v) =
 (v transformed invztr) - (origin transformed invztr)
enddef;

%% Initial Setup.

% active_plane is the active drawing plane.
% currentpicture is unknown at this stage
% (because it's set in beginchar),
% so use a def, not a picture assignment.

def active_plane = currentpicture enddef;

def initpic =
 setztr;
 def active_plane = currentpicture enddef;
 % Set drawpen to standard shape and size.
 interim penwd := 0.5pt;
 drawpen := pencircle scaled penwd yscaled aspect_ratio;
 % Sets currentpen to drawpen.
 pickup drawpen;
 if ClipOn:
   ClipPath := 1;
   % Specified character boundary.
   ClipPath[1] := rect (origin, (w,h));
 fi
 if debug:
   message " DEBUG ";
   message "Drawing Nominal Bounding Box Around MF Picture";
   safedraw rect (origin, (w,h));
   message " GUBED ";
 fi
enddef;

%% Compatibility with older graphbase.mf (for fig2dev's genmf.c)

def mfpicenv = enddef;
def endmfpicenv = enddef;

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

%% Character Wrapper.

def beginmfpic (expr ch) =
  beginchar (ch, (xpos-xneg)*xscale*unitlen, (ypos-yneg)*yscale*unitlen, 0);
  initpic;
enddef;

def endmfpic =
 if debug:
  message " DEBUG ";
  message "TFM charwd = " & decimal charwd & "pt#";
  message "TFM charht = " & decimal charht & "pt#";
  message " GUBED ";
 fi
 endchar;
enddef;


%% Extra Trigonometric and Hyperbolic Functions.
% tand, cotd, acos, asin;
% exp, ln, cosh, sinh, acosh, asinh.
% NOTE:  acos and asin return angle in degrees.

def tand (expr x) = (sind (x) / cosd (x)) enddef;
def cotd (expr x) = (cosd (x) / sind (x)) enddef;

def acos (expr x) = angle ((x, 1+-+x)) enddef;
def asin (expr y) = angle ((1+-+y, y)) enddef;

def exp primary X = (mexp (256 * X)) enddef;
def ln  primary X = (mlog (X) / 256) enddef;

vardef cosh primary X =
  save temp;
  temp = exp X;
  (temp + 1/temp) / 2;
enddef;

vardef sinh primary X =
  save temp;
  temp = exp X;
  (temp - 1/temp) / 2;
enddef;

def acosh (expr y) = ln (y + (y+-+1)) enddef;
def asinh (expr y) = ln (y + (y++1)) enddef;


%% Coordinate Systems and Transformations.

% Coordinate Nesting.
% Make this local, by using a (global) transform stack;
% and transparent to METAFONT syntax, by using hide.
% (For a primary expression, would use gobble instead of hide.)

transform T_stack[];
T_stack := 0;

def T_push (expr T) =
  T_stack[incr T_stack] := T;
enddef;

def T_pop (suffix $) =
  if T_stack > 0:
    $ := T_stack[T_stack];
    T_stack := T_stack - 1;
  fi
enddef;

def bcoords =
  hide (T_push (ztr))
enddef;

def ecoords =
  hide (T_pop (ztr))
enddef;

% Coordinate Changes.
% Example:  `apply_t (rotated theta);'
%   where Transformer is `rotated theta'.

def apply_t (text Transformer) =
  ztr := identity Transformer transformed ztr;
%  currenttransform := identity Transformer transformed currenttransform;
enddef;

let xslant = slanted;  % (x+sy, y).

def yslant primary s =  % (x, y+sx).
  transformed
   begingroup
     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
      transform T_;
      xpart T_ = ypart T_ = 0;
      xxpart T_ = yypart T_ = xpart p;
      xypart T_ = yxpart T_ = ypart p;
      T_
    endgroup
enddef;

def xyswap =  % swap X and Y axes.
  zslant (0,1)
enddef;

def boost primary X =  % boosts for special relativity.
  zslant (cosh X, sinh X)
enddef;

% Path Rotation.
% Rotate path f about point p by angle th,
% where f and p are in _graph_ coordinates,
% and th is in _sharp_ coordinates.

vardef rotatedpath (expr p,th) expr f =
 f transformed ztr rotatedaround (p transformed ztr, th) transformed invztr
enddef;


%% Bitmaps, Clipping and Rendering.

% Bitmap to Bitmap --- Bitwise Operations.

% NB: `cull' is *unknown* to METAPOST, plain METAPOST, and mfplain.

if metapost:
  % p being a pair.
  def keeping p = enddef;
  % v being a picture.
  def cull expr v = enddef;
fi

% Change a picture u to {0,1} (`monochrome') form.
% Apply this before the other monochrome operations.

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

% bitwise and.

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

% inclusive or.

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

% exclusive or.

def xorto (suffix u) (expr v) =
 addto u also v;
 cull u keeping (1,1);
enddef;
 
primarydef u picxor v =
 begingroup
  picture t;
  t=u+v;
  cull t keeping (1,1);
  t
 endgroup
enddef;
 
% (nonsymmetric) difference.

def subto (suffix u) (expr v) =
 addto u also -v;
 cull u keeping (1,1);
enddef;
 
primarydef u picsub v =
 begingroup
  picture t;
  t=u-v;
  cull t keeping (1,1);
  t
 endgroup
enddef;
 
% Contour to Bitmap --- Clipping and Filling.

% (safely filled) interior of contour c,
% where c is in _pixel_ coordinates;
% adapted from _The METAFONTbook_'s "safefill".

% plain.mf defines `.t_' as `transformed currenttransform';
% currenttransform compensates for the aspect ratio of _pixels_,
% and is set in mode_setup, which is called near the top of
% grafbase.mf; see above.
% (Reference:  _The METAFONTbook_, page 269.)

vardef interior expr c =
  save vs; picture vs; vs=nullpicture;
  interim turningcheck := 0;
  addto vs contour (c.t_) withpen nullpen;
  cull vs dropping (0,0);
  vs
enddef;

% (safely filled) interiors of contours cc[],
% where cc[] are in _pixel_ coordinates;

vardef interiors suffix cc =
 save vss; picture vss; vss=nullpicture;
 for i=1 upto cc:
  orto (vss, interior cc[i]);
 endfor
 vss
enddef;

% derived bitmap operations.

% clip picture vt to interior of cycle c,
% where c is in _pixel_ coordinates.

def clipto (suffix vt) expr c =
 mono (vt);
 andto (vt, interior c);
enddef;

% clip picture vt to interiors of cycles in path array cc,
% where c is in _pixel_ coordinates.

def clipsto (suffix vt, cc) =
 mono (vt);
 andto (vt, interiors cc);
enddef;


% clip copy of picture vt to interior of cycle c,
% where c is in _pixel_ coordinates,
% return result (which is a `subpicture' of vt).

vardef clip (suffix vt) expr c =
  mono (vt);
  vt picand (interior c)
enddef;

% fill region inside c in picture vt,
% where c is in _pixel_ coordinates.

vardef picfill (suffix vt) expr c =
  mono (vt);
  orto (vt, interior c);
enddef;

% unfill region inside c in picture vt,
% where c is in _pixel_ coordinates.

vardef picunfill (suffix vt) expr c =
  mono (vt);
  subto (vt, interior c);
enddef;

% (return) reverse video of vt inside c,
% where c is in _pixel_ coordinates.

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

% Rendering Paths --- Drawing and Filling.

% draw path f in picture v using pen q,
% where f is in _pixel_ coordinates.

def shpath (suffix v) (expr q) (expr f) =
  interim turningcheck := 0;
  addto v doublepath (f.t_)
    withpen q;
enddef;

% draw path d safely, return picture of drawing,
% where d is in _pixel_ coordinates.
% (Courtesy of Uwe Bonnes.)

numeric minpenwd;
minpenwd := 0.01pt;

vardef picpath expr d =
 save vs; picture vs; vs=nullpicture;
 if penwd > minpenwd:
   shpath (vs, drawpen) (d);
   mono (vs);
 fi
 vs
enddef;

% drawing a path d safely in picture vt,
% where d is in _pixel_ coordinates.

def picdraw (suffix vt) expr d =
 addto vt also (picpath d);
enddef;

% drawing a path d safely,
% where d is in _pixel_ coordinates.

def safedraw expr d =
  picdraw (active_plane) d;
  if ClipOn:
    clipsto (active_plane, ClipPath);
  fi
enddef;

% filling a region, the interior of c, safely,
% where c is in _pixel_ coordinates.

def safefill expr c =
  if cycle c:
    picfill (active_plane) c;
    if ClipOn:
      clipsto (active_plane, ClipPath);
    fi
  else:  % so we can see something
    safedraw c;
  fi
enddef;

% erasing a region (the interior of c) safely,
% where c is in _pixel_ coordinates.
% - really works this time -

def safeunfill expr c =
  if cycle c:
    picunfill (active_plane) c;
    if ClipOn:
      clipsto (active_plane, ClipPath);
    fi
  else:  % so we can see something
    safedraw c;
  fi
enddef;

% draw path f and return path f,
% where f is in _graph_ coordinates.

vardef drawn expr f =
 safedraw zconv (f);
 f
enddef;

% fill contour c and return contour c,
% where c is in _graph_ coordinates.

vardef filled expr c =
 safefill zconv (c);
 c
enddef;

% unfill contour c and return contour c,
% where c is in _graph_ coordinates.

vardef unfilled expr c =
 safeunfill zconv (c);
 c
enddef;


%% Shading and Hatching.

% Shading
% - expressed mainly in _pixel_ coordinates.

% A closed path representing one dot of unit size
% (one pixel across is recommended).
path dotpath;

% Initially, set dotpath to a circle.
dotpath := fullcircle;

% Return picture that renders specified apath path to specified scale.
% We use ceiling to ensure visibility.

vardef setdot (expr apath) (expr scale) =
  if cycle apath:
    interior
  else:
    if debug:
      errhelp ("I hope you don't mind");
      errmessage ("setdot with open `boundary', drawing curve instead");
    fi
    picpath
  fi
  (apath scaled ceiling (scale))
enddef;

% Picture containing one dot.
picture onedot;

% Initially, set onedot to picture of initial dotpath scaled 0.5pt.
onedot := setdot (dotpath, 0.5pt);

% Draw picture w at position p in picture v,
% where p is in _pixel_ coordinates.

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

% Draw the onedot picture
% at position p (in _pixel_ coordinates)
% in active drawing plane.

def pixdot (expr p) =
  picdot (active_plane, onedot, p);
enddef;

% Draw the onedot picture
% at position p (in _graph_ coordinates)
% in active drawing plane.

def ourdot (expr p) =
  pixdot (zconv (p));
enddef;


% Whether to show bounding boxes (for debugging purposes).

boolean showbbox;
showbbox := false;


% Calculate tight bounding box points (ll, ur) for path g.
% g is in _pixel_ coordinates.

vardef tightbbox (expr g) (suffix ll, ur) =
  % vertical tangent:
  def vtang (expr x) =
    (x, -infinity)--(x, infinity)
  enddef;
  % horizontal tangent:
  def htang (expr y) =
    (-infinity, y)--(infinity, y)
  enddef;
  % xlimit(x) = true iff vtang(x) does not intersect path g:
  vardef xlimit (expr x) =
    (xpart (g intersectiontimes (vtang(x))) < 0)
  enddef;
  % ylimit(y) = true iff htang(y) does not intersect path g:
  vardef ylimit (expr y) =
    (xpart (g intersectiontimes (htang(y))) < 0)
  enddef;
  save minx, maxx, miny, maxy;
  numeric minx, maxx, miny, maxy;
  minx := (solve xlimit (-infinity, xpart point 1 of g));
  maxx := (solve xlimit ( infinity, xpart point 1 of g));
  miny := (solve ylimit (-infinity, ypart point 1 of g));
  maxy := (solve ylimit ( infinity, ypart point 1 of g));
  ll := (minx, miny);
  ur := (maxx, maxy);
  if showbbox:
    safedraw (rect (ll, ur));
  fi
enddef;

% calculate tight bounding box points (ll, ur) for path array g.
% g is in _pixel_ coordinates.

vardef tbbox (suffix g) (suffix ll, ur) =
  save gll, gur;
  pair gll, gur;
  for i=1 upto g:
    tightbbox (g[i], gll[i], gur[i]);
  endfor
  gll = gur = g;
  ll := minpair (gll);
  ur := maxpair (gur);
  if showbbox:
    safedraw (rect (ll, ur));
  fi
enddef;

% calculate loose bounding box points (ll, ur) for path g.
% g is in _pixel_ coordinates.

vardef bbox (expr g) (suffix ll, ur) =
  ur := ll := point 0 of g;
  save p; pair p[];
  for i=0 upto length g:
    p1 := point i of g;
    p2 := precontrol i of g;
    p3 := postcontrol i of g;
    p := 4;
    p4 := ll;
    ll := minpair (p);
    p4 := ur;
    ur := maxpair (p);
  endfor
  if showbbox:
    safedraw (rect (ll, ur));
  fi
enddef;


% shade interior of path f with dots,
% of pre-specified shape and size (see setdot macro),
% spaced sp apart,
% where f is in _graph_ coordinates,
% and sp is in _pixel_ coordinates.
% Return f.

vardef shade (expr sp) expr f =
 save g; path g; g := zconv (f);
 if not cycle g:  % so we can see something
  safedraw g;
  message "attempt to shade open curve; drawing curve instead";
 elseif sp<=0:
  safefill g;
 else:
  save ll, ur; pair ll, ur;
  bbox (g, ll, ur);
  ll := sp * (ceilingpair (ll/sp));
  save mn; pair mn;
  mn := floorpair ((ur-ll)/sp);
  m  := xpart mn;
  n  := ypart mn;
  twosp := 2*sp;
  save p; pair p[];
  p2 := ll;
  picture v;
  v := nullpicture;
  for i=0 upto m:
   p3 := p2 if odd i: + (0,sp) fi;
   for j=0 upto n:
    if (not odd (i+j)):
     picdot (v, onedot, p3);
     p3 := p3 + (0,twosp);
    fi
   endfor
   p2:=p2+(sp,0);
  endfor
  addto active_plane also
    clip(v) g;
 fi
 f
enddef;


% Given:
% a picture variable v_c,
% an affine transform CT,
% a spacing sp,
% and an upright box (with edges xa, xb, ya, yb)
% specified in the coordinate system that CT determines.
% Then:
% hatch the interior of the box in the picture,
% using lines that in that coordinate system are horizontal
% and spaced sp units apart.
% Notes:
% 1.  sp must be strictly positive.
% 2.  CT, xa, xb, ya, yb, and sp are in _pixel_ coordinates.
% 3.  xa, xb, ya, and yb are the left, right, bottom and top
%     edges of the box.

vardef thatchf (suffix v_c) (expr CT, sp, xa, xb, ya, yb) =
  % Number of hatch lines.
  save m; numeric m;
  m := ceiling (abs (yb - ya) / sp);
  save p; pair p[];  % In Pixel Coords.
  % Hatch lines have same alignment, spacing and orientation.
  % Displacement from one hatch line to the next.
  p1 := sp * up;
  % Y coordinate of first hatch line.
  save yc; numeric yc;
  yc := sp * ceiling (ya/sp);
  % Endpoints of first hatch line.
  p2 := (xa, yc);
  p3 := (xb, yc);
  for i=1 upto m:
   % The next source line does all the drawing.
   shpath (v_c, hatchpen) ((p2--p3) transformed CT);
   p2 := p2+p1;
   p3 := p3+p1;
  endfor
enddef;

% Hatch interior of path f
% with lines at angle theta,
% spaced sp apart.
% where f is in _graph_ coordinates,
% and sp is in _pixel_ coordinates.

vardef thatch (expr sp, theta) expr f =
 % g is in _pixel_ coordinates.
 save g; path g; g := zconv (f);
 if not cycle g:  % so we can see something
  message "attempt to hatch open curve; drawing curve instead";
  safedraw g;
 elseif sp<=0:
  safefill g;
 else:
  save v_c; picture v_c;
  v_c := nullpicture;
  save CT; transform CT;
  CT := identity rotated theta;
  save ll, ur; pair ll, ur;
  bbox (g transformed inverse CT, ll, ur);
  save xa, xb, ya, yb; numeric xa, xb, ya, yb;
  (xa, ya) = ll;
  (xb, yb) = ur;
  % Hatch interior of bounding box.
  thatchf (v_c, CT, sp, xa, xb, ya, yb);
  % Clip to interior of path g.
  addto active_plane also
    clip(v_c) g;
%    v_c;
 fi
 f
enddef;


% Hatch interior of path f with horizontal lines spaced sp apart.
% where f is in _graph_ coordinates,
% and sp is in _pixel_ coordinates.

vardef hhatch (expr sp) expr f =
  thatch (sp, 0 deg) f
enddef;

% Hatch interior of path f with vertical lines spaced sp apart.
% where f is in _graph_ coordinates,
% and sp is in _pixel_ coordinates.

vardef vhatch (expr sp) expr f =
  thatch (sp, 90 deg) f
enddef;

% Left-hatch interior of path f with lines spaced sp apart,
% where f is in _graph_ coordinates,
% and sp is in _pixel_ coordinates.
% Return f.

vardef lhatch (expr sp) expr f =
  thatch (sp, -45 deg) f
enddef;

% Right-hatch interior of path f with lines spaced sp apart,
% where f is in _graph_ coordinates,
% and sp is in _pixel_ coordinates.
% Return f.

vardef rhatch (expr sp) expr f =
  thatch (sp, 45 deg) f
enddef;

% Cross-hatch interior of path f with lines spaced sp apart,
% where f is in _graph_ coordinates,
% and sp is in _pixel_ coordinates.

vardef xhatch (expr sp) expr f =
 lhatch (sp) rhatch (sp) f
enddef;


%% Tiles.

% environment for a tile `atile', drawn with unit length
% `unit', with width `width' units and height `height' units.
% If `clipon' is true, clip the drawing to the tile's
% (width x height) boundary.

def tile (suffix atile) (expr unit, width, height, clipon) =
 picture atile.pic; atile.pic=nullpicture;
 % for reference by `tess' macro:
 numeric atile.wd, atile.ht;
 atile.wd=width;
 atile.ht=height;
 begingroup
  save active_plane;
  def active_plane = atile.pic enddef;
  save wtr; transform wtr;
  wtr = identity scaled unit;
  save ClipOn; boolean ClipOn;
  if clipon:
   ClipOn := true;
   save ClipPath; path ClipPath[];
   ClipPath = 1;
   ClipPath[1] =
     ((0,0)--(width,0)--(width,height)--(0,height)--cycle)
     transformed wtr;
  else:
   ClipOn := false;
  fi
enddef;

def endtile =
 endgroup
enddef;

% test whether atile is really a tile.

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

% tesselation of interior of closed path c,
% using tile `atile'.
% c is in _graph_ units.
% return c.

vardef tess (suffix atile) expr c =
 save g; path g; g := zconv (c);
 if cycle g:
  save ll, ur; pair ll, ur;
  bbox (g, ll, ur);
  save a, b; numeric a, b;
  if not is_tile (atile):
    errhelp ("Please check the faulty definition of the `tile'");
    errmessage ("tess given invalid tile");
  fi
  a = atile.wd;
  b = atile.ht;
  save m[], n[];
  numeric m[], n[];
  m1 = floor   ((xpart ll) / a);
  m2 = ceiling ((xpart ur) / a);
  n1 = floor   ((ypart ll) / b);
  n2 = ceiling ((ypart ur) / b);
  save vs;
  picture vs;
  vs := nullpicture;
  % fill vs with rectangular tesselation large enough
  % to cover interior of g:
  for i=m1 upto m2:
   for j=n1 upto n2:
    addto vs also (atile.pic shifted (i*a, j*b);
   endfor
  endfor
  % trim tesselation to interior of g:
  clipto (vs) g;
  addto active_plane also vs;
 else:  % so we can see something
  safedraw g;
  message "attempt to tesselate open curve; drawing curve instead";
 fi
 c
enddef;

%% Dots and Dashes.

% Draw dashed curve along path f;
% f is in _graph_ coordinates;
% return f.

vardef dashed (expr dlen, dgap) expr f =
 save dashsize; numeric dashsize; dashsize := dlen;
 save dashgap; numeric dashgap; dashgap := dgap;
 save g; path g; g := zconv (f);
 save dp; numeric dp; path dp[];
 dashlist (dp) (g);
 for i=1 upto dp:
   safedraw dp[i];
 endfor
 f
enddef;

% Draw dotted curve along path f;
% f is in _graph_ coordinates;
% return f.

vardef dotted (expr dsize, dgap) expr f =
 save dotgap; numeric dotgap; dotgap := dgap;
 save g; path g; g := zconv (f);
 save dp; numeric dp; pair dp[];
 dotlist (dp) (g);
 onedot := setdot (dotpath, dsize);
 for i=1 upto dp:
   pixdot (dp[i]);
 endfor
 f
enddef;


%% Points.

% ptwd is a numeric scale factor,
% b is a pair;
% return a path representing
% a point
% with diameter ptwd
% located at b.
% All quantities are in _pixel_ units. 

vardef bpoint (expr ptwd, b) =
  fullcircle scaled ptwd shifted b
enddef;

% Draw discs with _pixel_ diameter ptwd
% at the _graph_ locations listed in the text t.
% If filled is true, make the discs black,
% otherwise make them white with a black boundary.

vardef pointd (expr ptwd, filled) (text t) =
 save f; path f;  % in _graph_ coords.
 for a=t:
  f := bpoint (ptwd, zconv (a));
  if filled:
    safefill f;
  else:
    safeunfill f;
    safedraw f;
  fi
 endfor
enddef;

%% Arrows.
% arrowheads are in _pixel_ coordinates.

newinternal hdwdr, hdten;
interim hdwdr := 1;
interim hdten := 1;

boolean hfilled;
hfilled := false;

% wr is in _pixel_ coordinates;
% tens is a pure number;
% fil is _boolean_.

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

% initial arrowhead shape.
headshape (1, 1, false);

% front, back and hdwdr are in _pixel_ coords.

vardef head (expr front, back, hdwdr, tens, filled) =
 save p, side;
 pair p[], side;
 side := (hdwdr/2) *
   ((front-back) rotated 90);
 p1 := back + side;
 p2 := back - side;
 save f; path f;
 f :=
   p1..tension tens..{front-back}front{back-front}..tension tens..p2;
 if filled:
   safefill (f--cycle);
 fi
 safedraw f;
enddef;

% hlen, hrot and hback are in _pixel_ coords.
% f is in _graph_ coords.
% return f.

vardef headpath (expr hlen, hrot, hback) expr f =
 save g; path g; g := zconv (f);
 save P; pair P[];
 P2 := point     (length g) of g;
 P1 := direction (length g) of g;
 if P1 <> (0,0) :
  P3 := ((unitvector (P1)) rotated hrot);
  P4 := P2 - (hback * P3);
  P5 := P4 - (hlen  * P3);
  head (P4, P5, hdwdr, hdten, hfilled);
  if debug:
    message " DEBUG ";
    message "P1 to P5 :";
    show P1, P2, P3, P4, P5;
    message " GUBED ";
  fi
 fi
 f
enddef;

% hlen is in _pixel_ coordinates;
% f is in _graph_ coords.

def arrowdraw (expr hlen) (expr f) =
 safedraw zconv (f);
 store (trash) headpath (hlen, 0, 0pt) f;
enddef;


%% Axes, Axis Tic Marks, and Grid.

% hlen is in _pixel_ coords.
% axes are described by xneg, xpos, yneg and ypos,
% which are in _graph_ coords.

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

% Axis Tic Marks are at right angles to axes,
% and have length len.

% len is in _pixel_ Y coords.
% t is a list of _graph_ X coords.

vardef xmarks (expr len) (text t) =
 % P is position on axis, tic is tic mark.
 save tic, P;
 pair tic, P;
 tic := (len/2) * unitvector (vconv (right)) rotated 90;
 for a=t:
   P := zconv (a * right);
   safedraw (((-tic)..tic) shifted P);
 endfor
enddef;

% len is in _pixel_ X coords.
% t is a list of _graph_ Y coords.

vardef ymarks (expr len) (text t) =
 save tic, P;
 pair tic, P;
 tic := (len/2) * unitvector (vconv (up)) rotated 90;
 for a=t:
   P := zconv (a * up);
   safedraw (((-tic)..tic) shifted P);
 endfor
enddef;

% Grid.

% Draw dots at each grid coordinate.
% xspace, yspace = spacings between grid coordinates, in _graph_ coords.

vardef grid (expr xspace, yspace) =
 for x=xneg step xspace until xpos:
 for y=yneg step yspace until ypos:
   ourdot ((x, y));
 endfor
 endfor
enddef;


%% Upright Rectangles.
% coordinate-independent.

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


%% Path Construction.
% coordinate-independent.

% general path construction

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


%% Polylines, including Polygons.
% coordinate-independent.

% polyline construction
% pts is both an array of pairs or paths,
% and a count of its members (from 1 to pts). 

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

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


%% Smooth Curves.
% coordinate-independent.

% smooth path construction

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

vardef curve (expr cyclic) (text t) =
 textpairs (p,t);
 mksmooth (cyclic, p)
enddef;

% quadratic B-splines.
% p[] == B-spline control points,
% p in number.

vardef openqbs (text t) =
  textpairs (p,t);
  for i=1 upto p-2:
    0.5[p[i],p[i+1]]
    ..controls 1/6[p[i+1],p[i]]
           and 1/6[p[i+1],p[i+2]]..
  endfor
  0.5[p[p-1],p[p]]
enddef;

vardef closedqbs (text t) =
  textpairs (p,t);
  p[p+1]:=p1;
  p[p+2]:=p2;
  for i=1 upto p:
    0.5[p[i],p[i+1]]
    ..controls 1/6[p[i+1],p[i]]
           and 1/6[p[i+1],p[i+2]]..
  endfor
  cycle
enddef;

% cubic B-splines.

vardef mkopencbs (suffix b) =
  for i=1 upto b-3:
    (b[i]+4b[i+1]+b[i+2])/6
    ..controls 1/3[b[i+1],b[i+2]]
           and 2/3[b[i+1],b[i+2]]..
  endfor
  (b[b-2]+4b[b-1]+b[b])/6
enddef;

vardef opencbs (text t) =
  textpairs (b,t);
  mkopencbs (b)
enddef;

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

vardef closedcbs (text t) =
  textpairs (b,t);
  mkclosedcbs (b)
enddef;


%% Path Closure.
% coordinate-independent.
%% QUESTION:  How to close open paths while retaining their original shape?

% path f closed by line segment.

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

% close path f in manner of mksmooth.

vardef sclosed expr f =
 if cycle f:
  f
 else:
  save n; numeric n;
  n := length f;  % number of Bezier segments in f
  if n <= 1:  % single Bezier segment (n=1), or single point (n=0)
   f..cycle
  else:  % (n >= 2)
   (point 0 of f){(point 1 of f)-(point infinity of f)}
   ..(subpath (1,n-1) of f)
   ..(point infinity of f){(point 0 of f)-(point (n-1) of f)}
   ..cycle
  fi
 fi
enddef;

% path f closed by bezier.

vardef bclosed expr f =
 f
 if not cycle f:
   ..cycle
 fi
enddef;

% conversion of Bezier segment key points, z,
% to cubic B-spline control points, b.

def ztob (suffix z, b) =
  save b; pair b[];
  b := 4;
  b1 = 6z1-7z2+2z3;
  b2 =     2z2- z3;
  b3 =    - z2+2z3;
  b4 =     2z2-7z3+6z4;
enddef;

% closure of path f by a cubic B-spline.

vardef cbclosed expr f =
 if cycle f:
  f
 else:
  save p; pair p[];
  p1 := point 0 of f;
  p2 := postcontrol 0 of f;
  p3 := precontrol 1 of f;
  p4 := point 1 of f;
  ztob (p,a);
  n := length f;
  p1 := point (n-1) of f;
  p2 := postcontrol (n-1) of f;
  p3 := precontrol infinity of f;
  p4 := point infinity of f;
  ztob (p,b);
  b1 := b3;
  b2 := b4;
  b3 := a1;
  b4 := a2;
  f..mkopencbs(b)..cycle
 fi
enddef;


%% Circles and Ellipses.
% coordinate-independent.

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

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


%% Circular Arcs.
% coordinate-independent.

% center-point-sweep form of arc

vardef arc (expr center, from, sweep) =
 save f; path f;
 f := from;
 if (center<>from) and (sweep<>0):
   n := ceiling (abs(sweep)/45) + 1;
   if n<3: n := 3; fi
   theta := sweep/(n-1);
   save p,q; pair p,q;
   p := from;
   for i=2 upto n:
     p := p rotatedabout (center, theta);
     q := (p-center) rotated 90;
     if theta<0:  q := -q;  fi
     f := f..p{unitvector q};
   endfor
 fi
 f
enddef;

% (0.5)[from,to] is mid point of chord across arc.
% cotd (sweep/2) is the displacement of ...

def arccenter (expr from, to, sweep) =
 ((0.5)[from,to]
  if (sweep mod 360) <> 0:
    + ((cotd (sweep/2)) * ((to-from) rotated 90) / 2)
  fi)
enddef;

% point-point-sweep form of arc

vardef arcpps (expr from, to, sweep) =
 save center; pair center;
 center := arccenter (from,to,sweep);
 arc (center, from, sweep)
enddef;

% modified polar --

% center, angle, angle, radius

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

% center-point-sweep form

def arccps (expr center, from, theta) =
 (arc (center, from, theta))
enddef;

% point-point-point form

vardef arcppp (expr first, second, third) =
 save sweep; numeric sweep;
 sweep:=2*(angle(third-second)-angle(second-first));
 sweep:=sweep mod 720;
 if sweep > 360: sweep:=sweep-720; fi
 critical:=5;
 if abs(sweep) <= critical:  % center may blow out
  save p; pair p[];
  p:=3;
  p1:=first;
  p2:=second;
  p3:=third;
  mksmooth (false,p)
 else:
  save m, d, center;
  pair m[], d[], center;
  m1:=(0.5)[first,second]; d1:=(second-first) rotated 90;
  m2:=(0.5)[second,third]; d2:=(third-second) rotated 90;
  center = m1+whatever*d1 = m2+whatever*d2;
  arc (center,first,sweep)
 fi
enddef;


%% Polar Coordinates.
% coordinate-independent.

% polar to rectangular conversion:  (r,th) -> (x,y) = r * dir (th).

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


%% Turtle.
% coordinate-independent.

% Argument t is a list of graph coord. pairs.
% First item in list t is starting point,
% rest are successive moves.

vardef turtle (text t) =
 save p; pair p[]; p:=0;
 for a=t:
  p[incr p] := if p<>1:  p[p-1] +  fi  a;
 endfor
 mkpoly (false, p)
enddef;


%% Sectors.
% coordinate-independent.

vardef sector (expr center, rad, frtheta, totheta) =
 save from; pair from;
 if debug:
   message " DEBUG ";
   message "sector's center, rad, frtheta, totheta: ";
   show center, rad, frtheta, totheta;
   message " GUBED ";
 fi
 from:=center+rad*(dir frtheta);
 if debug:
   message " DEBUG ";
   message "sector's from: ";
   show from;
   message " GUBED ";
 fi
 center -- arc (center,from,totheta-frtheta) -- cycle
enddef;


%% Utility Functions.
% coordinate-independent.

% the identity function.
def id (expr x) =
  x
enddef;


%% Plotting of functions.
% coordinate-independent.

% In these macros, if the boolean argument `smooth' then the path
% returned will be a B\'ezier, otherwise it will be a polyline.
% In mkfcn, for example, the parameter values determining the points
% through which the path interpolates begin at bmin, and step by bst
% until bmax is reached or exceeded.

% mkfcn: convert function `pf : (x, y) = pf (t)' to a path.
% pf : numeric -> pair.
% domain interpreted as parameter values,
% codomain as cartesian coordinate pairs.

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

let tfcn = mkfcn;

% parafcn: like tfcn, but the text argument is a literal pair
% expression `(fa(t), fb(t))' in the literal `t'.

vardef parafcn (expr smooth) (expr bmin, bmax, bst) (text pf) =
 save p; pair p[]; p := 0;
 for t = bmin step bst until bmax+(bst/2):
   p[incr p] := pf;
 endfor
 mkpath (smooth, false, p)
enddef;

% xfcn: convert function `fx : y = f (x)' to a path.
% fx : numeric -> numeric.
% domain interpreted as X coordinate,
% codomain as Y coordinate.
% Return the path.

vardef xfcn (expr smooth) (expr xmin, xmax, st) (text fx) =
   def fp (expr x) =
     (x, fx(x))
   enddef;
   mkfcn (smooth) (xmin, xmax, st) (fp)
enddef;

% function: like xfcn, but the text argument is a literal numeric
% expression `fx(x)' in the literal `x'.

vardef function (expr smooth) (expr xmin, xmax, st) (text fx) =
   def fp (expr x) =
     (x, fx)
   enddef;
   mkfcn (smooth) (xmin, xmax, st) (fp)
enddef;

% rfcn: convert function `ft : r = ft (\theta)' to a path.
% ft : numeric -> numeric.
% domain interpreted as angle in degrees,
% codomain as radial distance from origin.

vardef rfcn (expr smooth) (expr tmin, tmax, st) (text ft) =
   def fq (expr t) =
     (ft(t) * (dir t))
   enddef;
   mkfcn (smooth) (tmin, tmax, st) (fq)
enddef;

% plrfcn: like rfcn, but the text argument is a literal numeric
% expression `ft(t)' in the literal `t' which represents the
% angle \theta.

vardef plrfcn (expr smooth) (expr tmin, tmax, st) (text ft) =
   def fq (expr t) =
     (ft * (dir t))
   enddef;
   mkfcn (smooth) (tmin, tmax, st) (fq)
enddef;

%% Overlays - taken from MFbook, p 295.  (Bruce Leban)

picture totalpicture;
boolean totalnull, currentnull;

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

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

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

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

def shipit =
  mergeit (shipout)
enddef;

def showit_ =
  mergeit (show_)
enddef;

def show_ suffix v =
  display v inwindow currentwindow
enddef;

%%%
%%%  end  grafbase.mf
%%%