%% %% 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 %% %% 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'.