diff options
Diffstat (limited to 'metapost/context/base/mpiv/mp-grap.mpiv')
| -rw-r--r-- | metapost/context/base/mpiv/mp-grap.mpiv | 1706 |
1 files changed, 1706 insertions, 0 deletions
diff --git a/metapost/context/base/mpiv/mp-grap.mpiv b/metapost/context/base/mpiv/mp-grap.mpiv new file mode 100644 index 000000000..4fd8ee5bd --- /dev/null +++ b/metapost/context/base/mpiv/mp-grap.mpiv @@ -0,0 +1,1706 @@ +%D \module +%D [ file=mp-grap.mpiv, +%D version=2012.10.16, % 2008.09.08 and earlier, +%D title=\CONTEXT\ \METAPOST\ graphics, +%D subtitle=graph packagesupport, +%D author=Hans Hagen \& Alan Braslau, +%D date=\currentdate, +%D copyright={PRAGMA ADE \& \CONTEXT\ Development Team}] +%C +%C This module is part of the \CONTEXT\ macro||package and is +%C therefore copyrighted by \PRAGMA. See licen-en.pdf for +%C details. + +if known context_grap : endinput ; fi ; + +boolean context_grap ; context_grap := true ; + +% Below is a modified graph.mp + +show numbersystem, numberprecision ; + +%if epsilon/4 = 0 : +if numbersystem <> "double" : + errmessage "The graph macros require the double precision number system." ; + endinput ; +fi + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% $Id : graph.mp,v 1.2 2004/09/19 21 :47 :10 karl Exp $ +% Public domain. + +% Macros for drawing graphs + +% begingraph(width,height) begin a new graph +% setcoords(xtype,ytype) sets up a new coordinate system (log,-linear..) +% setrange(lo,hi) set coord ranges (numeric and string args OK) +% gdraw <file or path> [with...] draw a line in current coord system +% gfill <file or path> [with...] fill a region using current coord system +% gdrawarrow .., gdrawdblarrow.. like gdraw, but with 1 or 2 arrowheads +% augment<path name>(loc) append given coordinates to a polygonal path +% glabel<suffix>(pic,loc) place label pic near graph coords or time loc +% gdotlabel<suffix>(pic,loc) same with dot +% OUT loc value for labels relative to whole graph +% gdata(file,s,text) read coords from file ; evaluate t w/ tokens s[] +% auto.<x or y> default x or y tick locations (for interation) +% tick.<bot|top|..>(fmt,u) draw centered tick from given side at u w/ format +% itick.<bot|top|..>(fmt,u) draw inward tick from given side at u w/ format +% otick.<bot|top|..>(fmt,u) draw outward tick at coord u ; label format fmt +% grid.<bot|top|..>(fmt,u) draw grid line at u with given side labeled +% autogrid([itick|.. bot|..],..) iterate over auto.x, auto.y, drawing tick/grids +% frame.[bot|top..] draw frame (or one side of the frame) +% graph_frame_needed := false ; after begingraph, not to draw a frame at all +% graph_background := color ; fill color for frame, if defined +% endgraph end of graph--the result is a picture + +% option `plot <picture>' draws picture at each path knot, turns off pen +% graph_template.<tickcmd> template paths for tick marks and grid lines +% graph_margin_fraction.low, +% graph_margin_fraction.high fractions determining margins when no setrange +% graph_log_marks[], graph_lin_marks, graph_exp_marks loop text strings used by auto.<x or y> +% graph_minimum_number_of_marks, graph_log_minimum numeric parameters used by auto.<x or y> +% Autoform is the format string used by autogrid +% Autoform_X, Autoform_Y if defined, are used instead + +% Other than the above-documented user interface, all externally visible names +% are of the form X_.<suffix>, Y_.<suffix>, or Z_.<suffix>, or they start +% with `graph_' + +% Used to depend on : + +% input string.mp + +% Private version of a few marith macros, fixed for double math... + +newinternal Mzero ; Mzero := -16384; % Anything at least this small is treated as zero +newinternal mlogten ; mlogten := mlog(10) ; +newinternal largestmantissa ; largestmantissa := 2**52 ; % internal double warningcheck +newinternal singleinfinity ; singleinfinity := 2**128 ; +newinternal doubleinfinity ; doubleinfinity := 2**1024 ; +%Mzero := -largestmantissa ; % Note that we get arithmetic overflows if we set to -doubleinfinity + +% Safely convert a number to mlog form, trapping zero. + +vardef graph_mlog primary x = + if unknown x: whatever + elseif x=0: Mzero + else: mlog(abs x) fi +enddef ; + +vardef graph_exp primary x = + if unknown x: whatever + elseif x<=Mzero: 0 + else: mexp(x) fi +enddef ; + +% and add the following for utility/completeness +% (replacing the definitions in mp-tool.mpiv). + +vardef logten primary x = + if unknown x: whatever + elseif x=0: Mzero + else: mlog(abs x)/mlog(10) fi +enddef ; + +vardef ln primary x = + if unknown x: whatever + elseif x=0: Mzero + else: mlog(abs x)/256 fi +enddef ; + +vardef exp primary x = + if unknown x: whatever + elseif x<= Mzero: 0 + else: (mexp 256)**x fi +enddef ; + +vardef powten primary x = + if unknown x: whatever + elseif x<= Mzero: 0 + else: 10**x fi +enddef ; + +% Convert x from mlog form into a pair whose xpart gives a mantissa and whose +% ypart gives a power of ten. + +vardef graph_Meform(expr x) = + if x<=Mzero : origin + else : + save e, m ; e=floor(x/mlogten)-3; m := mexp(x-e*mlogten) ; + if abs m<1000 : m := m*10 ; e := e-1 ; elseif abs m>=10000 : m := m/10 ; e := e+1 ; fi + (m, e) + fi +enddef ; + +% Modified from above. + +vardef graph_Feform(expr x) = + interim warningcheck :=0 ; + if x=0 : origin + else : + save e, m ; e=floor(if x<0 : -mlog(-x) else : mlog(x) fi/mlogten)-3; m := x/(10**e) ; + if abs m<1000 : m := m*10 ; e := e-1 ; elseif abs m>=10000 : m := m/10 ; e := e+1 ; fi + (m, e) + fi +enddef ; + +vardef graph_error(expr x,s) = + interim showstopping :=0 ; + show x ; errmessage s ; +enddef ; + +%%%%%%%%%%%%%%%%%%%%%%%% Data structures, begingraph %%%%%%%%%%%%%%%%%%%%%%%% + +vardef Z_@# = (X_@#,Y_@#) enddef ; % used in place of plain.mp's z convention + +def graph_suffix(suffix $) = % convert from x or y to X_ or Y_ + if str$="x" : X_ else : Y_ fi +enddef ; + +% New : + +save graph_background ; color graph_background ; % if defined, fill the frame. +save graph_close_file ; boolean graph_close_file ; graph_close_file = false ; + +def begingraph(expr w, h) = + begingroup + save X_, Y_ ; + X_.graph_coordinate_type = + Y_.graph_coordinate_type = linear ; % coordinate system for each axis + Z_.graph_dimensions = (w,h) ; % dimensions of graph not counting axes etc. + %also, Z_.low, Z_.high user-specified coordinate ranges in units used in graph_current_graph + + save graph_finished_graph ; + picture graph_finished_graph ; % the finished part of the graph + graph_finished_graph = nullpicture ; + save graph_current_graph ; + picture graph_current_graph ; % what has been drawn in current coords + graph_current_graph = nullpicture ; + save graph_current_bb ; + picture graph_current_bb ; % picture whose bbox is graph_current_graph's w/ linewidths 0 + graph_current_bb = nullpicture ; + save graph_last_drawn ; + picture graph_last_drawn ; % result of last gdraw or gfill + graph_last_drawn = nullpicture ; + save graph_last_path ; + path graph_last_path ; % last gdraw or gfill path in data coordinates. + save graph_plot_picture ; + picture graph_plot_picture ; % a picture from the `plot' option known when plot allowed + save graph_foreground ; + color graph_foreground ; % drawing color, if set. + save graph_label ; + picture graph_label[] ; % labels to place around the whole graph when it is done + save graph_autogrid_needed ; + boolean graph_autogrid_needed ; % whether autogrid is needed + graph_autogrid_needed = true ; + save graph_frame_needed ; + boolean graph_frame_needed ; % whether frame needs to be drawn + graph_frame_needed = true ; + save graph_number_of_arrowheads ; % number of arrowheads for next gdraw + graph_number_of_arrowheads = 0 ; + + if known graph_background : % new feature! + fill origin--(w,0)--(w,h)--(0,h)--cycle withcolor graph_background ; + fi +enddef ; + +% Additional variables not explained above : +% graph_modified_lower, graph_modified_higher pairs giving bounds used in auto<x or y> +% graph_exponent, graph_comma variables and macros used in auto<x or y> +% graph_modified_bias +% an offset to graph_modified_lower and graph_modified_higher to ease computing exponents +% Some additional variables function as constants. Most can be modified by the +% user to alter the behavior of these macros. +% Not very modifiable : log, linear, +% graph_frame_pair_a, graph_frame_pair_b, graph_margin_pair +% Modifiable : graph_template.suffix, +% graph_log_marks[], graph_lin_marks, graph_exp_marks, +% graph_minimum_number_of_marks, +% graph_log_minimum, Autoform + + +newinternal log, linear ; % coordinate system codes +log :=1 ; linear :=2; + +% note that mp-tool.mpiv defines log as log10. + +%%%%%%%%%%%%%%%%%%%%%% Coordinates : setcoords, setrange %%%%%%%%%%%%%%%%%%%%%% + +% Graph-related user input is `user graph coordinates' as specified by arguments +% to setcoords. +% `Internal graph coordinates' are used for graph_current_graph, graph_current_bb, Z_.low, Z_.high. +% Their meaning depends on the appropriate component of Z_.graph_coordinate_type : +% log means internal graph coords = mlog(user graph coords) +% -log means internal graph coords = -mlog(user graph coords) +% linear means internal graph coords = (user graph coords) +% -linear means internal graph coords = -(user graph coords) + + +vardef graph_set_default_bounds = % Set default Z_.low, Z_.high + forsuffixes $=low,high : + (if known X_$ : whatever else : X_$ fi, if known Y_$ : whatever else : Y_$ fi) + = graph_margin_fraction$[llcorner graph_current_bb,urcorner graph_current_bb] + + graph_margin_pair$ ; + endfor +enddef ; + +pair graph_margin_pair.low, graph_margin_pair.high ; +graph_margin_pair.high = -graph_margin_pair.low = (.00002,.00002) ; + +% Set $, $$, $$$ so that shifting by $ then transforming by $$ and then $$$ maps +% the essential bounding box of graph_current_graph into (0,0)..Z_.graph_dimensions. +% The `essential bounding box' is either what Z_.low and Z_.high imply +% or the result of ignoring pen widths in graph_current_graph. + +vardef graph_remap(suffix $,$$,$$$) = + save p_ ; + graph_set_default_bounds ; + pair p_, $ ; $=-Z_.low; + p_ = (max(X_.high-X_.low,.9), max(Y_.high-Y_.low,.9)) ; + transform $$, $$$ ; + forsuffixes #=$$,$$$ : xpart#=ypart#=xypart#=yxpart#=0 ; endfor + (Z_.high+$) transformed $$ = p_ ; + p_ transformed $$$ = Z_.graph_dimensions ; +enddef ; + +graph_margin_fraction.low=-.07 ; % bbox fraction for default range start +graph_margin_fraction.high=1.07 ; % bbox fraction for default range stop + +def graph_with_pen_and_color(expr q) = + withpen penpart q withcolor + if colormodel q=1 : + false + elseif colormodel q=3 : + (greypart q) + elseif colormodel q=5 : + (redpart q, greenpart q, bluepart q) + elseif colormodel q=7 : + (cyanpart q, magentapart q, yellowpart q, blackpart q) + fi +enddef ; + +% Add picture component q to picture @# and change part p to tp, +% where p is something from q that needs coordinate transformation. +% The type of p is pair or path. +% Pair o is the value of p that makes tp (0,0). This implements the trick +% whereby using 1 instead of 0 for the width or height or the setbounds path +% for a label picture suppresses shifting in x or y. + +%vardef graph_picture_conversion@#(expr q, o)(text tp) = +% save p ; +% if stroked q : +% path p ; p=pathpart q; +% addto @# doublepath tp graph_with_pen_and_color(q) dashed dashpart q ; +% elseif filled q : +% path p ; p=pathpart q; +% addto @# contour tp graph_with_pen_and_color(q) ; +% else : +% interim truecorners :=0 ; +% pair p ; p=llcorner q; +% if urcorner q<>p : p := p + graph_coordinate_multiplication(o-p,urcorner q-p) ; fi +% addto @# also q shifted ((tp)-llcorner q) ; +% fi +%enddef ; + +% This new version makes gdraw clip the result to the window defined with setrange + +vardef graph_picture_conversion@#(expr q, o)(text tp) = + save p ; + save do_clip, tp_clipped ; boolean do_clip ; do_clip := true ; + picture tp_clipped ; tp_clipped := nullpicture; + if stroked q : + path p ; p=pathpart q; + addto tp_clipped doublepath tp graph_with_pen_and_color(q) dashed dashpart q ; + %draw bbox tp_clipped withcolor red ; + elseif filled q : + path p ; p=pathpart q; + addto tp_clipped contour tp graph_with_pen_and_color(q) ; + %draw bbox tp_clipped withcolor green ; + else : + if (urcorner q<>llcorner q) : do_clip := false ; fi % Do not clip the axis labels; + interim truecorners := 0 ; + pair p ; p=llcorner q; + if urcorner q<>p : p := p + graph_coordinate_multiplication(o-p,urcorner q-p) ; fi + addto tp_clipped also q shifted ((tp)-llcorner q) ; + %draw bbox tp_clipped withcolor if do_clip : cyan else : blue fi ; + fi + if do_clip : + clip tp_clipped to origin--(xpart Z_.graph_dimensions,0)--Z_.graph_dimensions-- + (0,ypart Z_.graph_dimensions)--cycle ; + fi + addto @# also tp_clipped ; +enddef ; + +def graph_coordinate_multiplication(expr a,b) = (xpart a*xpart b, ypart a*ypart b) enddef ; + +vardef graph_clear_bounds@# = numeric @#.low, @#.high ; enddef; + +% Finalize anything drawn in the present coordinate system and set up a new +% system as requested + +vardef setcoords(expr tx, ty) = + interim warningcheck :=0 ; + if length graph_current_graph>0 : + save s, S, T ; + graph_remap(s, S, T) ; + for q within graph_current_graph : + graph_picture_conversion.graph_finished_graph(q,-s,p shifted s transformed S transformed T) ; + endfor + graph_current_graph := graph_current_bb := nullpicture ; + fi + graph_clear_bounds.X_ ; graph_clear_bounds.Y_; + X_.graph_coordinate_type := tx ; Y_.graph_coordinate_type := ty; +enddef ; + +% Set Z_.low and Z_.high to correspond to given range of user graph +% coordinates. The text argument should be a sequence of pairs and/or strings +% with 4 components in all. + +vardef setrange(text t) = + interim warningcheck :=0 ; + save r_ ; r_=0; + string r_[]s ; + for x_= + for p_=t : if pair p_ : xpart p_, ypart fi p_, endfor : + r_[incr r_] if string x_ : s fi = x_ ; + if r_>2 : + graph_set_bounds if r_=3 : X_ else : Y_ fi (r_[r_-2] if unknown r_[r_-2] : s fi, x_) ; + fi + exitif r_=4 ; + endfor +enddef ; + +% @# is X_ or Y_ ; l and h are numeric or string + +vardef graph_set_bounds@#(expr l, h) = + graph_clear_bounds@# ; + if @#graph_coordinate_type>0 : + @#low = if unknown l : + whatever + else : + if abs @#graph_coordinate_type=log : graph_mlog fi if string l : scantokens fi l + fi ; + @#high = if unknown h : + whatever + else : + if abs @#graph_coordinate_type=log : graph_mlog fi if string h : scantokens fi h + fi ; + else : + -@#high = if unknown l : + whatever + else : + if abs @#graph_coordinate_type=log : graph_mlog fi if string l : scantokens fi l + fi ; + -@#low = if unknown h : + whatever + else : + if abs @#graph_coordinate_type=log : graph_mlog fi if string h : scantokens fi h + fi ; + fi +enddef ; + +%%%%%%%%%%%%%%%%%%%%%%%%% Converting path coordinates %%%%%%%%%%%%%%%%%%%%%%%%% + +% Find the result of scanning path p and using macros tx and ty to adjust the +% x and y parts of each coordinate pair. Boolean parameter c tells whether to +% force the result to be polygonal. + +vardef graph_scan_path(expr p, c)(suffix tx, ty) = + if (str tx="") and (str ty="") : p + else : + save r_ ; path r_; + r_ := graph_pair_adjust(point 0 of p, tx, ty) + if path p : + for t=1 upto length p : + if c : -- + else : ..controls graph_pair_adjust(postcontrol(t-1) of p, tx, ty) + and graph_pair_adjust(precontrol t of p, tx, ty) .. + fi + graph_pair_adjust(point t of p, tx, ty) + endfor + if cycle p : &cycle fi + fi ; + if pair p : point 0 of fi r_ + fi +enddef ; + +vardef graph_pair_adjust(expr p)(suffix tx, ty) = (tx xpart p, ty ypart p) enddef ; + +% Convert path p from user graph coords to internal graph coords. + +vardef graph_convert_user_path_to_internal primary p = + interim warningcheck :=0 ; + if known p : + graph_scan_path(p, + (abs X_.graph_coordinate_type<>linear) or (abs Y_.graph_coordinate_type<>linear), + if abs X_.graph_coordinate_type=log : graph_mlog fi, + if abs Y_.graph_coordinate_type=log : graph_mlog fi) + transformed (identity + if X_.graph_coordinate_type<0 : xscaled -1 fi + if Y_.graph_coordinate_type<0 : yscaled -1 fi) + fi +enddef ; + +% Convert label location t_ from user graph coords to internal graph coords. +% The label location should be a pair, or two numbers/strings. If t_ is empty +% or a single item of non-pair type, just return t_. Unknown coordinates +% produce unknown components in the result. + +vardef graph_label_convert_user_to_internal(text t_) = + save n_ ; n_=0; + interim warningcheck :=0 ; + if 0 for x_=t_ : +1 if pair x_ : +1 fi endfor <= 1 : + t_ + else : + n_0 = n_1 = 0 ; + point 0 of graph_convert_user_path_to_internal ( + for x_= + for y_=t_ : if pair y_ : xpart y_, ypart fi y_, endfor + 0, 0 : + if known x_ : if string x_ : scantokens fi x_ + else : hide(n_[n_] :=whatever) 0 + fi + exitif incr n_=2 ; + ,endfor) + (n_0,n_1) + fi +enddef ; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Reading data files %%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% Read a line from file f, extract whitespace-separated tokens ignoring any +% initial "%", and return true if at least one token is found. The tokens +% are stored in @#1, @#2, .. with "" in the last @#[] entry. + +% String manipulation routines for MetaPost +% It is harmless to input this file more than once. + +vardef isdigit primary d = + ("0"<=d)and(d<="9") +enddef ; + +% Number of initial characters of string s where `c <character>' is true + +vardef graph_cspan(expr s)(text c) = + 0 + for i=1 upto length s: + exitunless c substring (i-1,i) of s; + + 1 + endfor +enddef ; + +% String s is composed of items separated by white space. Lop off the first +% item and the surrounding white space and return just the item. + +vardef graph_loptok suffix s = + save t, k; + k = graph_cspan(s," ">=); + if k > 0 : + s := substring(k,infinity) of s ; + fi + k := graph_cspan(s," "<); + string t; + t = substring (0,k) of s; + s := substring (k,infinity) of s; + s := substring (graph_cspan(s," ">=),infinity) of s; + t +enddef ; + +vardef graph_read_line@#(expr f) = + save n_, s_ ; string s_; + s_ = readfrom f ; + string @#[] ; + if s_<>EOF : + @#0 := s_ ; + @#1 := graph_loptok s_ ; + n_ = if @#1="%" : 0 else : 1 fi ; + forever : + @#[incr n_] := graph_loptok s_ ; + exitif @#[n_]="" ; + endfor + @#1<>"" + else : false + fi +enddef ; + +% Execute c for each line of data read from file f, and stop at the first +% line with no data. Commands c can use line number i and tokens $1, $2, ... +% and j is the number of fields. + +def gdata(expr f)(suffix $)(text c) = + %boolean flag ; % not used? + for i=1 upto largestmantissa : + exitunless graph_read_line$(f) ; + c + endfor + if graph_close_file : + closefrom f ; + fi +enddef ; + +% Read a path from file f. The path is terminated by blank line or EOF. + +vardef graph_readpath(expr f) = + interim warningcheck :=0 ; + save s ; + gdata(f, s, if i>1 :--fi + if s2="" : ( i, scantokens s1) + else : (scantokens s1, scantokens s2) fi + ) +enddef ; + +% Append coordinates t to polygonal path @#. The coordinates can be numerics, +% strings, or a single pair. + +vardef augment@#(text t) = + interim warningcheck := 0 ; + if not path begingroup @# endgroup : + graph_error(begingroup @# endgroup, "Cannot augment--not a path") ; + else : + def graph_comma= hide(def graph_comma=,enddef) enddef ; + if known @# : @# :=@#-- else : @#= fi + (for p=t : + graph_comma if string p : scantokens fi p + endfor) ; + fi +enddef ; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Drawing and filling %%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% Unknown pair components are set to 0 because glabel and gdotlabel understand +% unknown coordinates as `0 in absolute units'. + +vardef graph_unknown_pair_bbox(expr p) = + interim warningcheck:=0 ; + if known p : addto graph_current_bb doublepath p ; + else : + save x,y ; + z = llcorner graph_current_bb ; + if unknown xpart p : xpart p= else : x := fi 0 ; + if unknown ypart p : ypart p= else : y := fi 0 ; + addto graph_current_bb doublepath (p+z) ; + fi + graph_current_bb := image(fill llcorner graph_current_bb..urcorner graph_current_bb--cycle) ; +enddef ; + +% Initiate a gdraw or gfill command. This must be done before scanning the +% argument, because that could invoke the `if known graph_plot_picture' test in a following +% plot option . + +def graph_addto = + def graph_errorbar_text = enddef ; + color graph_foreground ; + path graph_last_path ; + graph_last_drawn := graph_plot_picture := nullpicture ; addto graph_last_drawn +enddef; + +% Handle the part of a gdraw command that uses path or data file p. + +def graph_draw expr p = + if string p : hide(graph_last_path := graph_readpath(p) ;) + graph_convert_user_path_to_internal graph_last_path + elseif path p or pair p : + hide(graph_last_path := p ;) + graph_convert_user_path_to_internal p + else : graph_error(p,"gdraw argument should be a data file or a path") + origin + fi + withpen currentpen graph_withlist _op_ +enddef ; + +% Handle the part of a gdraw command that uses path or data file p. + +def graph_fill expr p = + if string p : hide(graph_last_path := graph_readpath(p) --cycle ;) + graph_convert_user_path_to_internal graph_last_path + elseif cycle p : hide(graph_last_path := p ;) + graph_convert_user_path_to_internal p + else : graph_error(p,"gfill argument should be a data file or a cyclic path") + origin..cycle + fi graph_withlist _op_ +enddef ; + +def gdraw = graph_addto doublepath graph_draw enddef ; +def gfill = graph_addto contour graph_fill enddef ; + +% This is used in graph_draw and graph_fill to allow postprocessing graph_last_drawn + +def graph_withlist text t_ = t_ ; graph_post_draw; enddef; + +def witherrorbars(text t) text options = + hide( + def graph_errorbar_text = t enddef ; + save pic ; picture pic ; pic := image(draw origin _op_ options ;) ; + if color colorpart pic : graph_foreground := colorpart pic ; fi + ) + options +enddef ; + +% new feature: graph_errorbars + +picture graph_errorbar_picture ; graph_errorbar_picture := image(draw (left--right) scaled .5 ;) ; +%picture graph_xbar_picture ; graph_xbar_picture := image(draw (down--up) scaled .5 ;) ; +%picture graph_ybar_picture ; graph_ybar_picture := image(draw (left--right) scaled .5 ;) ; + +vardef graph_errorbars(text t) = + if known graph_last_path : + save n, p, q ; path p ; pair q ; + save pic ; picture pic[] ; pic0 := nullpicture ; + pic1 := if known graph_xbar_picture : graph_xbar_picture + elseif known graph_errorbar_picture : graph_errorbar_picture rotated 90 + else : nullpicture fi ; + pic2 := if known graph_ybar_picture : graph_ybar_picture + elseif known graph_errorbar_picture : graph_errorbar_picture + else : nullpicture fi ; + if length pic1>0 : + pic1 := pic1 scaled graph_shapesize ; + setbounds pic1 to origin..cycle ; + fi + if length pic2>0 : + pic2 := pic2 scaled graph_shapesize ; + setbounds pic2 to origin..cycle ; + fi + for i=0 upto length graph_last_path : + clearxy ; z = point i of graph_last_path ; + n := 1 ; + for $=t : + if known $ : + q := if path $ : if length $>i : point i of $ else : origin fi + elseif pair $ : $ elseif numeric $ : ($,$) else : origin fi ; + if q<>origin : + p := graph_convert_user_path_to_internal (( + if n=1 : + (-xpart q,0)--(ypart q,0) + else : + (0,-xpart q)--(0,ypart q) + fi ) shifted z) ; + addto pic0 doublepath p ; + if length pic[n]>0 : + if ypart q<>0 : + addto pic0 also pic[n] shifted point 1 of p ; + fi + if xpart q<>0 : + addto pic0 also pic[n] rotated 180 shifted point 0 of p ; + fi + fi + fi + fi + exitif incr n>3 ; + endfor + endfor + if length pic0>0 : + save bg, fg ; color bg, fg ; + bg := if known graph_background : graph_background else : background fi ; + fg := if known graph_foreground : graph_foreground else : black fi ; + addto graph_current_graph also pic0 withpen currentpen scaled 2 _op_ withcolor bg ; + addto graph_current_graph also pic0 withpen currentpen scaled .5 _op_ withcolor fg ; + fi + fi +enddef ; + +% Set graph_plot_picture so the postprocessing step will plot picture p at each path knot. +% Also select nullpen to suppress stroking. + +def plot expr p = + if known graph_plot_picture : + withpen nullpen + hide (graph_plot_picture := image( + if bounded p : for q within p : graph_addto_currentpicture q endfor % Save memory + else : graph_addto_currentpicture p + fi graph_setbounds origin..cycle)) + fi +enddef ; + +% This hides a semicolon that could prematurely end graph_withlist's text argument + +def graph_addto_currentpicture primary p = addto currentpicture also p ; enddef; +def graph_setbounds = setbounds currentpicture to enddef ; + +def gdrawarrow = graph_number_of_arrowheads := 1 ; gdraw enddef; +def gdrawdblarrow = graph_number_of_arrowheads := 2 ; gdraw enddef; + +% Post-process the filled or stroked picture graph_last_drawn as follows : (1) update +% the bounding box information ; (2) transfer it to graph_current_graph unless the pen has +% been set to nullpen to disable stroking ; (3) plot graph_plot_picture at each knot. + +vardef graph_post_draw = + save p ; path p ; p = pathpart graph_last_drawn ; + graph_unknown_pair_bbox(p) ; + if filled graph_last_drawn or not graph_is_null(penpart graph_last_drawn) : + addto graph_current_graph also graph_last_drawn ; + fi + graph_errorbars(graph_errorbar_text) ; + if length graph_plot_picture>0 : + for i=0 upto length p if cycle p : -1 fi : + addto graph_current_graph also graph_plot_picture shifted point i of p ; + endfor + picture graph_plot_picture ; + fi + if graph_number_of_arrowheads>0 : + graph_draw_arrowhead(p, graph_with_pen_and_color(graph_last_drawn)) ; + if graph_number_of_arrowheads>1 : + graph_draw_arrowhead(reverse p, graph_with_pen_and_color(graph_last_drawn)) ; + fi + graph_number_of_arrowheads := 0 ; + fi +enddef ; + +vardef graph_is_null(expr p) = (urcorner p=origin) and (llcorner p=origin) enddef ; + +vardef graph_draw_arrowhead(expr p)(text w) = % Draw arrowhead for path p, with list w + %save r ; r := angle(precontrol infinity of p shifted -point infinity of p) ; + addto graph_current_graph also + image(fill arrowhead (graph_arrowhead_extent(precontrol infinity of p,point infinity of p)) w ; + draw arrowhead (graph_arrowhead_extent(precontrol infinity of p,point infinity of p)) w + undashed ; +%if (r mod 90 <> 0) : % orientation can be wrong due to remapping +% draw textext("\tfxx " & decimal r) shifted point infinity of p withcolor blue ; +%fi + graph_setbounds point infinity of p..cycle ; + ) ; % rotatedabout(point infinity of p,-r) ; +enddef ; + +vardef graph_arrowhead_extent(expr p, q) = + if p<>q : (q - 100pt*unitvector(q-p)) -- fi + q +enddef ; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Drawing labels %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% Argument c is a drawing command that needs an additional argument p that gives +% a location in internal graph coords. Draw in graph_current_graph enclosed in a setbounds +% path. Unknown components of p cause the setbounds path to have width or height 1 instead of 0. +% Then graph_unknown_pair_bbox sets these components to 0 and graph_picture_conversion +% suppresses subsequent repositioning. + +def graph_draw_label(expr p)(suffix $)(text c) = + save sdim_ ; pair sdim_; + sdim_ := (if unknown xpart p : 1+ fi 0, if unknown ypart p : 1+ fi 0) ; + graph_unknown_pair_bbox(p) ; + addto graph_current_graph also + image(c(p) ; graph_setbounds p--p+sdim_--cycle) _op_ +enddef ; + +% Stash the result drawing command c in the graph_label table using with list w and +% an index based on angle mfun_laboff$. + +vardef graph_stash_label(suffix $)(text c) text w = + graph_label[1.5+angle mfun_laboff$ /90] = image(c(origin) w) ; +enddef ; + +def graph_label_location primary p = + if pair p : graph_draw_label(p) + elseif numeric p : graph_draw_label(point p of pathpart graph_last_drawn) + else : graph_stash_label + fi +enddef ; + +% Place label p at user graph coords t using with list w. (t is a time, a pair +% or 2 numerics or strings). + +vardef glabel@#(expr p)(text t) text w = + graph_label_location graph_label_convert_user_to_internal(t) (@#,label@#(p)) w ; enddef; + +% Place label p at user graph coords t using with list w and draw a dot there. +% (t is a time, a pair, or 2 numerics or strings). + +vardef gdotlabel@#(expr p)(text t) text w = + graph_label_location graph_label_convert_user_to_internal(t) (@#,dotlabel@#(p)) w ; enddef; + +def OUT = enddef ; % location text for outside labels + +%%%%%%%%%%%%%%%%%%%%%%%%%% Grid lines, ticks, etc. %%%%%%%%%%%%%%%%%%%%%%%%%% + +% Grid lines and tick marks are transformed versions of the templates below. +% In the template paths, (0,0) is on the edge of the frame and inward is to +% the right. + +path graph_template.tick, graph_template.itick, graph_template.otick, graph_template.grid ; +graph_template.tick = (-3.5bp,0)--(3.5bp,0) ; +graph_template.itick = origin--(7bp,0) ; +graph_template.otick = (-7bp,0)--origin ; +graph_template.grid = origin--(1,0) ; + +vardef tick@#(expr f,u) text w = graph_tick_label(@#,@,false,f,u,w) ; enddef; + +vardef itick@#(expr f,u) text w = graph_tick_label(@#,@,false,f,u,w) ; enddef; + +vardef otick@#(expr f,u) text w = graph_tick_label(@#,@,false,f,u,w) ; enddef; + +vardef grid@#(expr f,u) text w = graph_tick_label(@#,@,true,f,u,w) ; enddef; + + +% Produce a tick or grid mark for label suffix $, graph_template suffix $$, +% coordinate value u, and with list w. Boolean c tells whether graph_template$$ +% needs scaling by X_.graph_dimensions or Y_.graph_dimensions, +% and f gives a format string or a label picture. + +def graph_tick_label(suffix $,$$)(expr c, f, u)(text w) = + graph_draw_label(graph_label_convert_user_to_internal(graph_generate_label_position($,u)),, + draw graph_gridline_picture$($$,c,f,u,w) shifted) +enddef ; + +% Generate label positioning arguments appropriate for label suffix $ and +% coordinate u. + +def graph_generate_label_position(suffix $)(expr u) = + if pair u : u elseif xpart mfun_laboff.$=0 : u,whatever else : whatever,u fi +enddef ; + +% Generate a picture of a grid line labeled with coordinate value u, picture +% or format string f, and with list w. Suffix @# is bot, top, lft, or rt, +% suffix $ identifies entries in the graph_template table, and boolean c tells +% whether to scale graph_template$. + +vardef graph_gridline_picture@#(suffix $)(expr c, f, u)(text w) = + if unknown u : graph_error(u,"Label coordinate should be known") ; nullpicture + else : + save p ; path p; + interim warningcheck :=0 ; + graph_autogrid_needed :=false ; + p = graph_template$ zscaled -mfun_laboff@# + if c : graph_xyscale fi + shifted (((.5 + mfun_laboff@# dotprod (.5,.5)) * mfun_laboff@#) graph_xyscale) ; + image(draw p w ; + label@#(if string f : format(f,u) else : f fi, point 0 of p)) + fi +enddef ; + +def graph_xyscale = xscaled X_.graph_dimensions yscaled Y_.graph_dimensions enddef ; + +% Draw the frame or the part corresponding to label suffix @# using with list w. + +vardef frame@# text w = + graph_frame_needed :=false ; + picture p_ ; + p_ = image(draw + if str@#<>"" : subpath round(angle mfun_laboff@#*graph_frame_pair_a+graph_frame_pair_b) of fi + unitsquare graph_xyscale w) ; + graph_draw_label((whatever,whatever),,draw p_ shifted) ; +enddef ; + +pair graph_frame_pair_a ; graph_frame_pair_a=(1,1)/90; % unitsquare subpath is linear in label angle +pair graph_frame_pair_b ; graph_frame_pair_b=(.75,2.25); + +%%%%%%%%%%%%%%%%%%%%%%%%%% Automatic grid selection %%%%%%%%%%%%%%%%%%%%%%%%%% + +string graph_log_marks[] ; % marking options per decade for logarithmic scales +string graph_lin_marks ; % mark spacing options per decade for linear scales +string graph_exp_marks ; % exponent spacing options for logarithmic scales +newinternal graph_minimum_number_of_marks, graph_log_minimum ; +graph_minimum_number_of_marks := 4 ; % minimum number marks generated by auto.x or auto.y +graph_log_minimum := mlog 3 ; % revert to uniform marks when largest/smallest < this + +def Gfor(text t) = for i=t endfor enddef ; % to shorten the mark templates below + +graph_log_marks[1]="1,2,5" ; +graph_log_marks[2]="1,1.5,2,3,4,5,7" ; +graph_log_marks[3]="1Gfor(6upto10 :,i/5)Gfor(5upto10 :,i/2)Gfor(6upto9 :,i)" ; +graph_log_marks[4]="1Gfor(11upto20 :,i/10)Gfor(11upto25 :,i/5)Gfor(11upto19 :,i/2)" ; +graph_log_marks[5]="1Gfor(21upto40 :,i/20)Gfor(21upto50 :,i/10)Gfor(26upto49 :,i/5)" ; +graph_lin_marks="10,5,2" ; % start with 10 and go down; a final `,1' is appended +graph_exp_marks="20,10,5,2,1" ; + +Ten_to0 = 1 ; +Ten_to1 = 10 ; +Ten_to2 = 100 ; +Ten_to3 = 1000 ; +Ten_to4 = 10000 ; + +% Determine the X_ or Y_ bounds on the range to be covered by automatic grid +% marks. Suffix @# is X_ or Y_. The result is log or linear to specify the +% type of grid spacing to use. Bounds are returned in variables local to +% begingraph..endgraph : pairs graph_modified_lower and graph_modified_higher +% are upper and lower bounds in +% `modified exponential form'. In modified exponential form, (x,y) means +% (x/1000)*10^y, where 1000<=abs x<10000. + +vardef graph_bounds@# = + interim warningcheck :=0 ; + save l, h ; + graph_set_default_bounds ; + if @#graph_coordinate_type>0 : (l,h) else : -(h,l) fi = (@#low, @#high) ; + if abs @#graph_coordinate_type=log : + graph_modified_lower := graph_Meform(l)+graph_modified_bias ; + graph_modified_higher := graph_Meform(h)+graph_modified_bias ; + if h-l >= graph_log_minimum : log else : linear fi + else : + graph_modified_lower := graph_Feform(l)+graph_modified_bias ; + graph_modified_higher := graph_Feform(h)+graph_modified_bias ; + linear + fi +enddef ; + +pair graph_modified_bias ; graph_modified_bias=(0,3); +pair graph_modified_lower, graph_modified_higher ; + +% Scan graph_log_marks[k] and evaluate tokens t for each m where l<=m<=h. + +def graph_scan_marks(expr k, l, h)(text t) = + for m=scantokens graph_log_marks[k] : + exitif m>h ; + if m>=l : t fi + endfor +enddef ; + +% Scan graph_log_marks[k] and evaluate tokens t for each m and e where m*10^e belongs +% between l and h (inclusive), where both l and h are in modified exponent form. + +def graph_scan_mark(expr k, l, h)(text t) = + for e=ypart l upto ypart h : + graph_scan_marks(k, if e>ypart l : 1 else : xpart l/1000 fi, + if e<ypart h : 10 else : xpart h/1000 fi, t) + endfor +enddef ; + +% Select a k for which graph_scan_mark(k,...) gives enough marks. + +vardef graph_select_mark = + save k ; + k = 0 ; + forever : + exitif unknown graph_log_marks[k+1] ; + exitif 0 graph_scan_mark(incr k, graph_modified_lower, graph_modified_higher, +1) + >= graph_minimum_number_of_marks ; + endfor + k +enddef ; + +% Try to select an exponent spacing from graph_exp_marks. If successful, set @# and +% return true + +vardef graph_select_exponent_mark@# = + numeric @# ; + for e=scantokens graph_exp_marks : + @# = e ; + exitif floor(ypart graph_modified_higher/e) - + floor(graph_modified_exponent_ypart(graph_modified_lower)/e) + >= graph_minimum_number_of_marks ; + numeric @# ; + endfor + known @# +enddef ; + +vardef graph_modified_exponent_ypart(expr p) = ypart p if xpart p=1000 : -1 fi enddef ; + +% Compute the mark spacing d between xpart graph_modified_lower and xpart graph_modified_higher. + +vardef graph_tick_mark_spacing = + interim warningcheck :=0 ; + save m, n, d ; + m = graph_minimum_number_of_marks ; + n = 1 for i=1 upto + (mlog(xpart graph_modified_higher-xpart graph_modified_lower) - mlog m)/mlogten : + *10 endfor ; + if n<=1000 : + for x=scantokens graph_lin_marks : + d = n*x ; + exitif 0 graph_generate_numbers(d,+1)>=m ; + numeric d ; + endfor + fi + if known d : d else : n fi +enddef ; + +def graph_generate_numbers(expr d)(text t) = + for m = d*ceiling(xpart graph_modified_lower/d) step d until xpart graph_modified_higher : + t + endfor +enddef ; + +% Evaluate tokens t for exponents e in multiples of d in the range determined +% by graph_modified_lower and graph_modified_higher. + +def graph_generate_exponents(expr d)(text t) = + for e = d*floor(graph_modified_exponent_ypart(graph_modified_lower)/d+1) + step d until d*floor(ypart graph_modified_higher/d) : t + endfor +enddef ; + +% Adjust graph_modified_lower and graph_modified_higher so their exponent parts match +% and they are in true exponent form ((x,y) means x*10^y). Return the new exponent. + +vardef graph_match_exponents = + interim warningcheck := 0 ; + save e ; + e+3 = if graph_modified_lower=graph_modified_bias : ypart graph_modified_higher + elseif graph_modified_higher=graph_modified_bias : ypart graph_modified_lower + else : max(ypart graph_modified_lower, ypart graph_modified_higher) fi ; + forsuffixes $=graph_modified_lower, graph_modified_higher : + $ := (xpart $ for i=ypart $ upto e+2 : /(10) endfor, e) ; + endfor + e +enddef ; + +% Assume e is an integer and either m=0 or 1<=abs(m)<10000. Find m*(10^e) +% and represent the result as a string if its absolute value would be at least +% 4096 or less than .1. It is OK to return 0 as a string or a numeric. + +vardef graph_factor_and_exponent_to_string(expr m, e) = + if (e>3)or(e<-4) : + decimal m & "e" & decimal e + elseif e>=0 : + if abs m<infinity/Ten_to[e] : + m*Ten_to[e] + else : decimal m & "e" & decimal e + fi + else : + save x ; x=m/Ten_to[-e]; + if abs x>=.1 : x else : decimal m & "e" & decimal e fi + fi +enddef ; + +def auto suffix $ = + hide(def graph_comma= hide(def graph_comma=,enddef) enddef) + if graph_bounds.graph_suffix($)=log : + if graph_select_exponent_mark.graph_exponent : + graph_generate_exponents(graph_exponent, + graph_comma graph_factor_and_exponent_to_string(1,e)) + else : + graph_scan_mark(graph_select_mark, graph_modified_lower, graph_modified_higher, + graph_comma graph_factor_and_exponent_to_string(m,e)) + fi + else : + hide(graph_exponent :=graph_match_exponents) + graph_generate_numbers(graph_tick_mark_spacing, + graph_comma graph_factor_and_exponent_to_string(m,graph_exponent)) + fi +enddef ; + +string Autoform ; Autoform = "%g"; + +%vardef autogrid(suffix tx, ty) text w = +% graph_autogrid_needed :=false ; +% if str tx<>"" : for x=auto.x : tx(Autoform,x) w ; endfor fi +% if str ty<>"" : for y=auto.y : ty(Autoform,y) w ; endfor fi +%enddef ; + +% We redefine autogrid, adding the possibility of differing X and Y +% formats. + +% string Autoform_X ; Autoform_X := "@.0e" ; +% string Autoform_Y ; Autoform_Y := "@.0e" ; + +vardef autogrid(suffix tx, ty) text w = + graph_autogrid_needed := false ; + if str tx <> "" : + for x=auto.x : + tx ( + if string Autoform_X : + if Autoform_X <> "" : + Autoform_X + else : + Autoform + fi + else : + Autoform + fi, + x + ) w ; + endfor + fi + if str ty <> "" : + for y=auto.y : + ty ( + if string Autoform_Y : + if Autoform_Y <> "" : + Autoform_Y + else : + Autoform + fi + else : + Autoform + fi, + y + ) w ; + endfor + fi +enddef ; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% endgraph %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +def endgraph = + if graph_autogrid_needed : autogrid(otick.bot, otick.lft) ; fi + if graph_frame_needed : frame ; fi + setcoords(linear,linear) ; + interim truecorners :=1 ; + for b=bbox graph_finished_graph : + setbounds graph_finished_graph to b ; + for i=0 step .5 until 3.5 : + if known graph_label[i] : + addto graph_finished_graph also graph_label[i] shifted point i of b ; + fi + endfor + endfor + graph_finished_graph + endgroup +enddef ; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% We format in luatex (using \mathematics{}) ... +% we could pass via variables and save escaping as that is inefficient + +if unknown context_mlib : + + vardef escaped_format(expr s) = + "" for n=0 upto length(s) : & + if ASCII substring (n,n+1) of s = 37 : + "@" + else : + substring (n,n+1) of s + fi + endfor + enddef ; + + vardef strfmt(expr f, x) = % maybe use mfun_ namespace + "\MPgraphformat{" & escaped_format(f) & "}{" & mfun_tagged_string(x) & "}" + enddef ; + + vardef varfmt(expr f, x) = % maybe use mfun_ namespace + "\MPformatted{" & escaped_format(f) & "}{" & mfun_tagged_string(x) & "}" + enddef ; + + vardef format (expr f, x) = textext(strfmt(f,x)) enddef ; + vardef formatted(expr f, x) = textext(varfmt(f,x)) enddef ; + +fi ; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% A couple of extensions : + +% Define a function plotsymbol() returning a picture : 10 different shapes, +% unfilled outline, interior filled with different shades of the background. +% This allows overlapping points on a plot to be more distinguishable. + +vardef graph_shapesize = (.33BodyFontSize) enddef ; + +path graph_shape[] ; % (internal) symbol path + +graph_shape[0] := (0,0) ; % point +graph_shape[1] := fullcircle ; % circle +graph_shape[2] := (up -- down) scaled .5 ; % vertical bar + +for i = 3 upto 9 : % polygons + graph_shape[i] := + for j = 0 upto i-1 : + (up scaled .5) rotated (360j/i) -- + endfor cycle ; +endfor + +graph_shape[12] := graph_shape[2] rotated +90 ; % horizontal line +graph_shape[22] := graph_shape[2] rotated +45 ; % backslash +graph_shape[32] := graph_shape[2] rotated -45 ; % slash +graph_shape[13] := graph_shape[3] rotated 180 ; % down triangle +graph_shape[23] := graph_shape[3] rotated -90 ; % right triangle +graph_shape[33] := graph_shape[3] rotated +90 ; % left triangle +graph_shape[14] := graph_shape[4] rotated +45 ; % square +graph_shape[15] := graph_shape[5] rotated 180 ; % down pentagon +graph_shape[16] := graph_shape[6] rotated +90 ; % turned hexagon +graph_shape[17] := graph_shape[7] rotated 180 ; +graph_shape[18] := graph_shape[8] rotated +22.5 ; + +numeric l ; + +for j = 5 upto 9 : + l := length(graph_shape[j]) ; + pair p[] ; + for i = 0 upto l : + p[i] = whatever [point i of graph_shape[j], + point (i+2 mod l) of graph_shape[j]] ; + p[i] = whatever [point (i+1 mod l) of graph_shape[j], + point (i+l-1 mod l) of graph_shape[j]] ; + endfor + graph_shape[20+j] := for i = 0 upto l : point i of graph_shape[j]--p[i]--endfor cycle ; +endfor + +path s ; s := graph_shape[4] ; +path q ; q := s scaled .25 ; +numeric l ; l := length(s) ; + +pair p[] ; + +graph_shape[24] := for i = 0 upto l-1 : + hide( + p[i] = whatever [point i of s, point (i+1 mod l) of s] ; + p[i] = whatever [point i of q, point (i-1+l mod l) of q] ; + p[i+l] = whatever [point i of s, point (i+1 mod l) of s] ; + p[i+l] = whatever [point i+1 of q, point (i+2 mod l) of q] ; + ) + point i of q -- p[i] -- p[i+l] -- +endfor cycle ; + +graph_shape[34] := graph_shape[24] rotated 45 ; + +% usage : gdraw p plot plotsymbol( 1,1) ; % a filled circle +% usage : gdraw p plot plotsymbol(14,0) ; % a square +% usage : gdraw p plot plotsymbol( 4,.5) ; % a 50% filled diamond + +def stars(expr f) = plotsymbol(25,f) enddef ; % a 5-point star +def points(expr f) = plotsymbol( 0,f) enddef ; +def circles(expr f) = plotsymbol( 1,f) enddef ; +def crosses(expr f) = plotsymbol(34,f) enddef ; +def squares(expr f) = plotsymbol(14,f) enddef ; +def diamonds(expr f) = plotsymbol( 4,f) enddef ; % a turned square +def uptriangles(expr f) = plotsymbol( 3,f) enddef ; +def downtriangles(expr f) = plotsymbol(13,f) enddef ; +def lefttriangles(expr f) = plotsymbol(33,f) enddef ; +def righttriangles(expr f) = plotsymbol(23,f) enddef ; + +% f (fill) is color, numeric or boolean, otherwise background. +def plotsymbol(expr n, f) text t = + if known graph_shape[n] : + image( + save bg, fg ; color bg, fg ; + bg := if known graph_background : graph_background else : background fi ; + save pic ; picture pic ; pic := image(draw origin _op_ t ;) ; + if color colorpart pic : graph_foreground := colorpart pic ; fi + fg := if known graph_foreground : graph_foreground else : black fi ; + save p ; path p ; p = graph_shape[n] scaled graph_shapesize ; + draw p withcolor bg withpen currentpen scaled 2 ; % halo + currentpen := currentpen scaled .5 ; + if cycle p : + fill p withcolor + if known f : + if color f : + f + elseif numeric f : + f[bg,fg] + elseif boolean f and f : + fg + else + bg + fi + else : + bg + fi ; + fi + draw p _op_ t ; + ) + else : + nullpicture + fi + t +enddef ; + +% standard resistance color code: rainbow sequence (from /usr/share/X11/rgb.txt) +color resistance_color[] ; string resistance_name[] ; +resistance_color0 = (0,0,0) ; resistance_name0 = "black" ; +resistance_color1 = (165/255,42/255,42/255) ; resistance_name1 = "brown" ; +resistance_color2 = (1,0,0) ; resistance_name2 = "red" ; +resistance_color3 = (1,165/255,0) ; resistance_name3 = "orange" ; +resistance_color4 = (1,1,0) ; resistance_name4 = "yellow" ; +resistance_color5 = (0,1,0) ; resistance_name5 = "green" ; +resistance_color6 = (0,0,1) ; resistance_name6 = "blue" ; +resistance_color7 = (148/255,0,211/255) ; resistance_name7 = "darkviolet" ; +resistance_color8 = (190/255,190/255,190/255) ; resistance_name8 = "gray" ; +resistance_color9 = (1,1,1) ; resistance_name9 = "white" ; + +%def rainbow(expr f) = +% ((abs(5f) mod 5) + 2 - floor((abs(5f) mod 5) + 2)) +% [resistance_color[ floor((abs(5f) mod 5) + 2)], +% resistance_color[ceiling((abs(5f) mod 5) + 2)]] +%enddef ; +def rainbow(expr f) = + hide(numeric n_ ; n_ = (abs(5f) mod 5) + 2 ;) + (n_-floor(n_))[resistance_color[floor n_],resistance_color[ceiling n_]] +enddef ; + +% The following extensions are not specific to graph and could be moved to metafun... + +% sort a path. Efficient en memory use, not so efficient in sorting long paths... + +vardef sortpath (suffix $) (text t) = % t can be "xpart", "ypart", "length", "angle", ... + if path $ : + if length $ > 0 : + save n, k ; n := length $ ; + for i=0 upto n : + k := i ; + for j=i+1 upto n : + if t (point j of $) < t (point k of $) : + k := j ; + fi + endfor + if k>i : + $ := if i>0 : subpath (0,i-1) of $ -- fi + point k of $ -- + subpath (i,k-1) of $ + if k<n : -- subpath (k+1,n) of $ fi + ; + fi + endfor + fi + fi +enddef ; + +% convert a polygon path to a smooth path (useful, e.g. as a guide to the eye) + +def smoothpath (suffix $) = + if path $ : + (for i=0 upto length $ : + if i>0 : .. fi + (point i of $) + endfor ) + fi +enddef ; + +% return a path of a function func(x) with abscissa running from f to t over n intervals + +def makefunctionpath (expr f, t, n) (text func) = + (for x=f step ((t-f)/(abs n)) until t : + if x<>f : -- fi + (x, func) + endfor ) +enddef ; + +% shift a path, point by point +% +% example : +% +% p1 := addtopath(p0,(.1normaldeviate,.1normaldeviate)) ; + +vardef addtopath (suffix p) (text t) = + if path p : + (for i=0 upto length p : + if i>0 : -- fi + hide(clearxy ; z = point i of p ;) z shifted t + endfor) + fi +enddef ; + +% return a new path of a function func(z) using the same abscissa as an existing path + +vardef functionpath (suffix p) (text func) = + (for i=0 upto length p : + if i>0 : .. fi + (hide(x := xpart(point i of p))x,func) %(hide(clearxy ; z = point i of p)x,func) + endfor ) +enddef ; + +% least-squares "fit" to a polynomial +% +% example : +% +% path p[] ; +% numeric a[] ; a0 := 1 ; a1 := .1 ; a2 := .01 ; a3 := .001 ; a4 := 0.0001 ; +% p0 := makefunctionpath(0,5,10,polynomial_function(a,4,x)) ; +% p1 := addtopath(p0,(0,.001normaldeviate)) ; +% gdraw p0 ; +% gdraw p1 plot plotsymbol(1,.5) ; +% +% numeric b[] ; +% polynomial_fit(p1, b, 4, 1) ; +% gdraw functionpath(p1,polynomial_function(b,4,x)) ; +% +% numeric c[] ; +% linear_fit(p1, c, 1) ; +% gdraw functionpath(p1,linear_function(c,x)) dashed evenly ; + +% a polynomial function : +% +% y = a0 + a1 * x + a2 * x^2 + ... + a[n] * x^n + +vardef polynomial_function (suffix $) (expr n, x) = + (for j=0 upto n : + $[j]*(x**j) endfor) % no ; +enddef ; + +% find the determinant of a (n+1)*(n+1) matrix ; indices run from 0 to n + +vardef det (suffix $) (expr n) = + hide( + numeric determinant ; determinant := 1 ; + save jj ; numeric jj ; + for k=0 upto n : + if $[k][k]=0 : + jj := -1 ; + for j=0 upto n : + if $[k][j]<>0 : + jj := j ; + exitif true ; + fi + endfor + if jj<0 : + determinant := 0 ; + exitif true ; + fi + for j=k upto n : % interchange the columns + temp := $[j][jj] ; + $[j][jj] := $[j][k] ; + $[j][k] := temp ; + endfor + determinant = -determinant ; + fi + exitif determinant=0 ; + determinant := determinant * $[k][k] ; + if k<n : % subtract row k from lower rows to get a diagonal matrix + for j=k+1 upto n : + for i=k+1 upto n : + $[j][i] := $[j][i]-$[j][k]*$[k][i]/$[k][k] ; + endfor + endfor + fi + endfor ; + ) + determinant % no ; +enddef ; + +numeric fit_chi_squared ; + +% least-squares fit of a polynomial $ of order n to a path p (unweighted) +% +% reference : P. R. Bevington, "Data Reduction and Error Analysis for the Physical +% Sciences", McGraw-Hill, New York 1969. + +vardef polynomial_fit (suffix p, $) (expr n) (text t) = + if not path p : + graph_error(p, "Cannot fit--not a path") ; + elseif length p < n : + graph_error(p, "Cannot fit--not enough points") ; + else : + fit_chi_squared := 0 ; + % calculate sums of the data + save sumx, sumy ; numeric sumx[], sumy[] ; + save w ; numeric w ; + for i=0 upto 2n : + sumx[i] := 0 ; + endfor + for i=0 upto n : + sumy[i] := 0 ; + endfor + for i=0 upto length p : + clearxy ; z = point i of p ; + w := 1 ; % weight + if known t : + if numeric t : + w := 1 if t<>0 : /(abs t) fi ; + elseif pair t : + if t<>origin : + w := 1/(abs t) ; + fi + elseif path t : + if length t>= i: + if point i of t<>origin : + w := 1/(abs point i of t) ; + fi + else : + w := 0 ; + fi ; + fi + fi + x1 := w ; + for j=0 upto 2n : + sumx[j] := sumx[j] + x1 ; + x1 := x1 * x ; + endfor + y1 := y * w ; + for j=0 upto n : + sumy[j] := sumy[j] + y1 ; + y1 := y1 * x ; + endfor + fit_chi_squared := fit_chi_squared + y*y*w ; + endfor + % construct matrices and calculate the polynomial coefficients + save m ; numeric m[][] ; + for j=0 upto n : + for k=0 upto n : + m[j][k] := sumx[j+k] ; + endfor + endfor + save delta ; numeric delta ; + delta := det(m,n) ; % this destroys the matrix m[][], which is OK + if delta = 0 : + fit_chi_squared := 0 ; + for j=0 upto n : + $[j] := 0 ; + endfor + else : + for i=0 upto n : + for j=0 upto n : + for k=0 upto n : + m[j][k] := sumx[j+k] ; + endfor + m[j][i] := sumy[j] ; + endfor + $[i] := det(m,n) / delta ; % matrix m[][] gets destroyed... + endfor + for j=0 upto n : + fit_chi_squared := fit_chi_squared - 2sumy[j]*$[j] ; + for k=0 upto n : + fit_chi_squared := fit_chi_squared + $[j]*$[k]*sumx[j+k] ; + endfor + endfor + % normalize by the number of degrees of freedom + fit_chi_squared := fit_chi_squared / (length(p) - n) ; % length(p)+1-(n+1) + fi + fi +enddef ; + +% y = a0 + a1 * x +% +% of course a line is just a polynomial of order 1 + +vardef linear_function (suffix $) (expr x) = polynomial_function($,1,x) enddef ; +vardef linear_fit (suffix p, $) (text t) = polynomial_fit(p, $, 1, t) ; enddef ; + +% and a constant is polynomial of order 0 + +vardef constant_function (suffix $) (expr x) = polynomial_function($,0,x) enddef ; +vardef constant_fit (suffix p, $) (text t) = polynomial_fit(p, $, 0, t) ; enddef ; + +% y = a1 * exp(a0*x) +% +% exp and ln defined in metafun + +vardef exponential_function (suffix $) (expr x) = $1*exp($0*x) enddef ; + +% since we take a log, this only works for positive ordinates + +vardef exponential_fit (suffix p, $) (text t) = + save a ; numeric a[] ; + save q ; path q[] ; % fit to the log of the ordinate + for i=0 upto length p : + clearxy ; z = point i of p ; + if y>0 : + augment.q0(x,ln(y)) ; + augment.q1( + if known t : + if numeric t : (0,ln(t)) + elseif pair t : (xpart t,ln(ypart t)) + elseif path t : + if length t>=i : + hide(z1 = point i of t;) + (x1,ln(y1)) + else : + origin + fi + fi + else : + (0,1) + fi ) ; + fi + endfor + linear_fit(q0,a,q1) ; + save e ; e := exp(sqrt(fit_chi_squared)) ; + fit_chi_squared := e * e ; + $0 := a1 ; + $1 := exp(a0) ; +enddef ; + +% y = a1 * x**a0 + +vardef power_law_function (suffix $) (expr x) = $1*(x**$0) enddef ; + +% since we take logs, this only works for positive abscissae and ordinates + +vardef power_law_fit (suffix p, $) (text t) = + save a ; numeric a[] ; + save q ; path q[] ; % fit to the logs of the abscissae and ordinates + for i=0 upto length p : + clearxy ; z = point i of p ; + if (x>0) and (y>0) : + augment.q0(ln(x),ln(y)) ; + augment.q1( + if known t : + if numeric t : (0,ln(t)) + elseif pair t : (ln(xpart t),ln(ypart t)) + elseif path t : + if length t>=i : + hide(z1 = point i of t) + (ln(x1),ln(y1)) + else : + origin + fi + fi + else : + (0,1) + fi ) ; + fi + endfor + linear_fit(q0,a,q1) ; + save e ; e := exp(sqrt(fit_chi_squared)) ; + fit_chi_squared := e * e ; + $0 := a1 ; + $1 := exp(a0) ; +enddef ; + +% gaussian : y = a2 * exp(-ln(2)*((x-a0)/a1)^2) +% +% a1 is the hwhm ; sigma := a1/sqrt(2ln(2)) or a1/1.17741 + +newinternal lntwo ; lntwo := ln(2) ; % brrr, why not inline it + +vardef gaussian_function (suffix $) (expr x) = + if $1 = 0 : + if x = $0 : $2 else : 0 fi + else : + $2 * exp(-lntwo*(((x-$0)/$1)**2)) + fi + if known $3 : + + $3 + fi +enddef ; + +% since we take a log, this only works for positive ordinates + +vardef gaussian_fit (suffix p, $) (text t) = + save a ; numeric a[] ; + save q ; path q[] ; % fit to the log of the ordinate + for i=0 upto length p : + clearxy ; z = point i of p ; + if y>0 : + augment.q0(x,ln(y)) ; + augment.q1( + if known t : + if numeric t : (0,ln(t)) + elseif pair t : (xpart t,ln(ypart t)) + elseif path t : + if length t>=i : + hide(z1 = point i of t) + (x1,ln(y1)) + else : + origin + fi + fi + else : + (0,1) + fi ) ; + fi + endfor + polynomial_fit(q0,a,2,q1) ; + save e ; e := exp(sqrt(fit_chi_squared)) ; + fit_chi_squared := e * e ; + $1 := sqrt(-lntwo/a2) ; + $0 := -.5a1/a2 ; + $2 := exp(a0-.25*a1*a1/a2) ; + $3 := 0 ; % polynomial_fit will NOT work with a non-zero background! +enddef ; + +% lorentzian: y = a2 / (1 + ((x - a0)/a1)^2) + +vardef lorentzian_function (suffix $) (expr x) = + if $1 = 0 : + if x = $0 : $2 else : 0 fi + else : + $2 / (1 + ((x - $0)/$1)**2) + fi + if known $3 : + + $3 + fi +enddef ; + +vardef lorentzian_fit (suffix p, $) (text t) = + save a ; numeric a[] ; + save q ; path q ; % fit to the inverse of the ordinate + for i=0 upto length p : + if ypart(point i of p)<>0 : + augment.q(xpart(point i of p), 1/ypart(point i of p)) ; + fi + endfor + polynomial_fit(q,a,2,if t <> 0 : 1/(t) else : 0 fi) ; + fit_chi_squared := 1/fit_chi_squared ; + $0 := -.5a1/a2 ; + $2 := 1/(a0-.25a1*a1/a2) ; + $1 := sqrt((a0-.25a1*a1/a2)/a2) ; + $3 := 0 ; % polynomial_fit will NOT work with a non-zero background! +enddef ; |