%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
message ("using number system " & numbersystem & " with precision " & decimal numberprecision) ;
if numbersystem <> "double" :
errmessage "The graph macros require the double precision number system." ;
endinput ;
fi
if known contextlmtxmode :
def _op_ = base_draw_options enddef ;
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 (logarithmic,-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 ;
if known contextlmtxmode :
% already defined
else :
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 ;
fi
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.
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!
addto graph_finished_graph contour
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 : logarithmic, 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 logarithmic, linear ; % coordinate system codes
logarithmic :=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 :
% logarithmic means internal graph coords = mlog(user graph coords)
% -logarithmic 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=nocolormodel :
% false
% elseif colormodel q=graycolormodel :
% (greypart q)
% elseif colormodel q=rgbcolormodel :
% (redpart q, greenpart q, bluepart q)
% elseif colormodel q=cmykcolormodel :
% (cyanpart q, magentapart q, yellowpart q, blackpart q)
% fi
% withprescript prescriptpart q
% withpostscript postscriptpart q
%enddef ;
def graph_with_pen_and_color(expr q) =
withproperties q
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=logarithmic : graph_mlog fi if string l : scantokens fi l
fi ;
@#high = if unknown h :
whatever
else :
if abs @#graph_coordinate_type=logarithmic : graph_mlog fi if string h : scantokens fi h
fi ;
else :
-@#high = if unknown l :
whatever
else :
if abs @#graph_coordinate_type=logarithmic : graph_mlog fi if string l : scantokens fi l
fi ;
-@#low = if unknown h :
whatever
else :
if abs @#graph_coordinate_type=logarithmic : 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=logarithmic : graph_mlog fi,
if abs Y_.graph_coordinate_type=logarithmic : 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) =
% for i=1 upto largestmantissa :
% exitunless graph_read_line$(f) ;
% c
% endfor ;
% enddef ;
def gdata(expr f)(suffix $)(text c) =
for i=1 upto largestmantissa :
exitunless graph_read_line$(f) ;
c
endfor
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) w)
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 logarithmic 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=logarithmic :
graph_modified_lower := graph_Meform(l)+graph_modified_bias ;
graph_modified_higher := graph_Meform(h)+graph_modified_bias ;
if h-l >= graph_log_minimum : logarithmic 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($)=logarithmic :
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 known contextlmtxmode :
% already defined
elseif 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
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 ;