mp-grap.mpiv /size: 59 Kb    last modification: 2021-10-28 13:50
1%D \module
2%D   [       file=mp-grap.mpiv,
3%D        version=2012.10.16, % 2008.09.08 and earlier,
4%D          title=\CONTEXT\ \METAPOST\ graphics,
5%D       subtitle=graph packagesupport,
6%D         author=Hans Hagen \& Alan Braslau,
7%D           date=\currentdate,
8%D      copyright={PRAGMA ADE \& \CONTEXT\ Development Team}]
9%C
10%C This module is part of the \CONTEXT\ macro||package and is
11%C therefore copyrighted by \PRAGMA. See licen-en.pdf for
12%C details.
13
14if known context_grap : endinput ; fi ;
15
16boolean context_grap ; context_grap := true ;
17
18% Below is a modified graph.mp
19
20message ("using number system " & numbersystem & " with precision " & decimal numberprecision) ;
21
22if numbersystem <> "double" :
23    errmessage "The graph macros require the double precision number system." ;
24    endinput ;
25fi
26
27if known contextlmtxmode :
28    def _op_ = base_draw_options enddef ;
29fi ;
30
31%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
32
33% $Id : graph.mp,v 1.2 2004/09/19 21 :47 :10 karl Exp $
34% Public domain.
35
36% Macros for drawing graphs
37
38% begingraph(width,height)       begin a new graph
39% setcoords(xtype,ytype)         sets up a new coordinate system (logarithmic,-linear..)
40% setrange(lo,hi)                set coord ranges (numeric and string args OK)
41% gdraw <file or path> [with...] draw a line in current coord system
42% gfill <file or path> [with...] fill a region using current coord system
43% gdrawarrow .., gdrawdblarrow.. like gdraw, but with 1 or 2 arrowheads
44% augment<path name>(loc)        append given coordinates to a polygonal path
45% glabel<suffix>(pic,loc)        place label pic near graph coords or time loc
46% gdotlabel<suffix>(pic,loc)     same with dot
47% OUT                            loc value for labels relative to whole graph
48% gdata(file,s,text)             read coords from file ; evaluate t w/ tokens s[]
49% auto.<x or y>                  default x or y tick locations (for interation)
50% tick.<bot|top|..>(fmt,u)       draw centered tick from given side at u w/ format
51% itick.<bot|top|..>(fmt,u)      draw inward tick from given side at u w/ format
52% otick.<bot|top|..>(fmt,u)      draw outward tick at coord u ; label format fmt
53% grid.<bot|top|..>(fmt,u)       draw grid line at u with given side labeled
54% autogrid([itick|.. bot|..],..) iterate over auto.x, auto.y, drawing tick/grids
55% frame.[bot|top..]              draw frame (or one side of the frame)
56% graph_frame_needed := false ;  after begingraph, not to draw a frame at all
57% graph_background := color ;    fill color for frame, if defined
58% endgraph                       end of graph--the result is a picture
59
60% option `plot <picture>'        draws picture at each path knot, turns off pen
61% graph_template.<tickcmd>            template paths for tick marks and grid lines
62% graph_margin_fraction.low,
63% graph_margin_fraction.high     fractions determining margins when no setrange
64% graph_log_marks[], graph_lin_marks, graph_exp_marks    loop text strings used by auto.<x or y>
65% graph_minimum_number_of_marks, graph_log_minimum                numeric parameters used by auto.<x or y>
66% Autoform                       is the format string used by autogrid
67% Autoform_X, Autoform_Y         if defined, are used instead
68
69% Other than the above-documented user interface, all externally visible names
70% are of the form X_.<suffix>, Y_.<suffix>, or Z_.<suffix>, or they start
71% with `graph_'
72
73% Used to depend on :
74
75% input string.mp
76
77% Private version of a few marith macros, fixed for double math...
78
79newinternal Mzero ;           Mzero := -16384; % Anything at least this small is treated as zero
80newinternal mlogten ;         mlogten         := mlog(10) ;
81newinternal largestmantissa ; largestmantissa := 2**52 ; % internal double warningcheck
82newinternal singleinfinity ;  singleinfinity  := 2**128 ;
83newinternal doubleinfinity ;  doubleinfinity  := 2**1024 ;
84%Mzero := -largestmantissa ; % Note that we get arithmetic overflows if we set to -doubleinfinity
85
86% Safely convert a number to mlog form, trapping zero.
87
88vardef graph_mlog primary x =
89  if unknown x: whatever
90  elseif x=0: Mzero
91  else: mlog(abs x) fi
92enddef ;
93
94vardef graph_exp  primary x =
95  if unknown x: whatever
96  elseif x<=Mzero: 0
97  else: mexp(x) fi
98enddef ;
99
100% and add the following for utility/completeness
101% (replacing the definitions in mp-tool.mpiv).
102
103vardef logten primary x =
104  if unknown x: whatever
105  elseif x=0: Mzero
106  else: mlog(abs x)/mlog(10) fi
107enddef ;
108
109if known contextlmtxmode :
110    % already defined
111else :
112
113    vardef ln primary x =
114      if unknown x: whatever
115      elseif x=0: Mzero
116      else: mlog(abs x)/256 fi
117    enddef ;
118
119    vardef exp primary x =
120      if unknown x: whatever
121      elseif x<= Mzero: 0
122      else: (mexp 256)**x fi
123    enddef ;
124
125fi
126
127vardef powten primary x =
128  if unknown x: whatever
129  elseif x<= Mzero: 0
130  else: 10**x fi
131enddef ;
132
133% Convert x from mlog form into a pair whose xpart gives a mantissa and whose
134% ypart gives a power of ten.
135
136vardef graph_Meform(expr x) =
137  if x<=Mzero : origin
138  else :
139    save e, m ; e=floor(x/mlogten)-3; m := mexp(x-e*mlogten) ;
140    if abs m<1000 : m := m*10 ; e := e-1 ; elseif abs m>=10000 : m := m/10 ; e := e+1 ; fi
141    (m, e)
142  fi
143enddef ;
144
145% Modified from above.
146
147vardef graph_Feform(expr x) =
148  interim warningcheck :=0 ;
149  if x=0 : origin
150  else :
151    save e, m ; e=floor(if x<0 : -mlog(-x) else : mlog(x) fi/mlogten)-3; m := x/(10**e) ;
152    if abs m<1000 : m := m*10 ; e := e-1 ; elseif abs m>=10000 : m := m/10 ; e := e+1 ; fi
153    (m, e)
154  fi
155enddef ;
156
157vardef graph_error(expr x,s) =
158  interim showstopping :=0 ;
159  show x ; errmessage s ;
160enddef ;
161
162%%%%%%%%%%%%%%%%%%%%%%%% Data structures, begingraph %%%%%%%%%%%%%%%%%%%%%%%%
163
164vardef Z_@# = (X_@#,Y_@#) enddef ; % used in place of plain.mp's z convention
165
166def graph_suffix(suffix $) = % convert from x or y to X_ or Y_
167  if str$="x" : X_ else : Y_ fi
168enddef ;
169
170% New :
171
172save graph_background ; color graph_background ; % if defined, fill the frame.
173
174def begingraph(expr w, h) =
175  begingroup
176  save X_, Y_ ;
177  X_.graph_coordinate_type =
178  Y_.graph_coordinate_type = linear ; % coordinate system for each axis
179  Z_.graph_dimensions = (w,h) ;       % dimensions of graph not counting axes etc.
180  %also, Z_.low, Z_.high  user-specified coordinate ranges in units used in graph_current_graph
181
182  save    graph_finished_graph ;
183  picture graph_finished_graph ; % the finished part of the graph
184          graph_finished_graph  = nullpicture ;
185  save    graph_current_graph ;
186  picture graph_current_graph ;  % what has been drawn in current coords
187          graph_current_graph   = nullpicture ;
188  save    graph_current_bb ;
189  picture graph_current_bb ;     % picture whose bbox is graph_current_graph's w/ linewidths 0
190          graph_current_bb      = nullpicture ;
191  save    graph_last_drawn ;
192  picture graph_last_drawn ;     % result of last gdraw or gfill
193          graph_last_drawn      = nullpicture ;
194  save    graph_last_path ;
195  path    graph_last_path ;      % last gdraw or gfill path in data coordinates.
196  save    graph_plot_picture ;
197  picture graph_plot_picture ;   % a picture from the `plot' option known when plot allowed
198  save    graph_foreground ;
199  color   graph_foreground ;     % drawing color, if set.
200  save    graph_label ;
201  picture graph_label[] ;        % labels to place around the whole graph when it is done
202  save    graph_autogrid_needed ;
203  boolean graph_autogrid_needed ; % whether autogrid is needed
204          graph_autogrid_needed  = true ;
205  save    graph_frame_needed ;
206  boolean graph_frame_needed ;    % whether frame needs to be drawn
207          graph_frame_needed     = true ;
208  save    graph_number_of_arrowheads ; % number of arrowheads for next gdraw
209          graph_number_of_arrowheads = 0 ;
210
211  if known graph_background : % new feature!
212    addto graph_finished_graph contour
213      origin--(w,0)--(w,h)--(0,h)--cycle withcolor graph_background ;
214  fi
215enddef ;
216
217% Additional variables not explained above :
218% graph_modified_lower, graph_modified_higher  pairs giving bounds used in auto<x or y>
219% graph_exponent, graph_comma                  variables and macros used in auto<x or y>
220% graph_modified_bias
221%   an offset to graph_modified_lower and graph_modified_higher to ease computing exponents
222% Some additional variables function as constants.  Most can be modified by the
223% user to alter the behavior of these macros.
224% Not very modifiable :  logarithmic, linear,
225%                        graph_frame_pair_a, graph_frame_pair_b, graph_margin_pair
226% Modifiable :           graph_template.suffix,
227%                        graph_log_marks[], graph_lin_marks, graph_exp_marks,
228%                        graph_minimum_number_of_marks,
229%                        graph_log_minimum, Autoform
230
231
232newinternal logarithmic, linear ;        % coordinate system codes
233logarithmic :=1 ; linear :=2;
234
235% note that mp-tool.mpiv defines log as log10.
236
237%%%%%%%%%%%%%%%%%%%%%% Coordinates : setcoords, setrange %%%%%%%%%%%%%%%%%%%%%%
238
239% Graph-related user input is `user graph coordinates' as specified by arguments
240% to setcoords.
241% `Internal graph coordinates' are used for graph_current_graph, graph_current_bb, Z_.low, Z_.high.
242% Their meaning depends on the appropriate component of Z_.graph_coordinate_type :
243%  logarithmic means internal graph coords =  mlog(user graph coords)
244% -logarithmic means internal graph coords = -mlog(user graph coords)
245%  linear      means internal graph coords =      (user graph coords)
246% -linear      means internal graph coords =     -(user graph coords)
247
248
249vardef graph_set_default_bounds =         % Set default Z_.low, Z_.high
250  forsuffixes $=low,high :
251    (if known X_$ : whatever else : X_$ fi, if known Y_$ : whatever else : Y_$ fi)
252        = graph_margin_fraction$[llcorner graph_current_bb,urcorner graph_current_bb] +
253          graph_margin_pair$ ;
254  endfor
255enddef ;
256
257pair graph_margin_pair.low, graph_margin_pair.high ;
258graph_margin_pair.high = -graph_margin_pair.low = (.00002,.00002) ;
259
260% Set $, $$, $$$ so that shifting by $ then transforming by $$ and then $$$ maps
261% the essential bounding box of graph_current_graph into (0,0)..Z_.graph_dimensions.
262% The `essential bounding box' is either what Z_.low and Z_.high imply
263% or the result of ignoring pen widths in graph_current_graph.
264
265vardef graph_remap(suffix $,$$,$$$) =
266  save p_ ;
267  graph_set_default_bounds ;
268  pair p_, $ ; $=-Z_.low;
269  p_ = (max(X_.high-X_.low,.9), max(Y_.high-Y_.low,.9)) ;
270  transform $$, $$$ ;
271  forsuffixes #=$$,$$$ : xpart#=ypart#=xypart#=yxpart#=0 ; endfor
272  (Z_.high+$) transformed $$ = p_ ;
273  p_ transformed $$$ = Z_.graph_dimensions ;
274enddef ;
275
276graph_margin_fraction.low=-.07 ;         % bbox fraction for default range start
277graph_margin_fraction.high=1.07 ;        % bbox fraction for default range stop
278
279%def graph_with_pen_and_color(expr q) =
280%  withpen penpart q
281%  withcolor
282%  if colormodel q=nocolormodel :
283%    false
284%  elseif colormodel q=graycolormodel :
285%    (greypart q)
286%  elseif colormodel q=rgbcolormodel :
287%    (redpart q, greenpart q, bluepart q)
288%  elseif colormodel q=cmykcolormodel :
289%    (cyanpart q, magentapart q, yellowpart q, blackpart q)
290%  fi
291%  withprescript  prescriptpart q
292%  withpostscript postscriptpart q
293%enddef ;
294
295def graph_with_pen_and_color(expr q) =
296  withproperties q
297enddef ;
298
299
300% Add picture component q to picture @# and change part p to tp,
301% where p is something from q that needs coordinate transformation.
302% The type of p is pair or path.
303% Pair o is the value of p that makes tp (0,0).  This implements the trick
304% whereby using 1 instead of 0 for the width or height or the setbounds path
305% for a label picture suppresses shifting in x or y.
306
307%vardef graph_picture_conversion@#(expr q, o)(text tp) =
308%  save p ;
309%  if stroked q :
310%    path p ; p=pathpart q;
311%    addto @# doublepath tp graph_with_pen_and_color(q) dashed dashpart q ;
312%  elseif filled q :
313%    path p ; p=pathpart q;
314%    addto @# contour tp graph_with_pen_and_color(q) ;
315%  else :
316%    interim truecorners :=0 ;
317%    pair p ; p=llcorner q;
318%    if urcorner q<>p : p := p + graph_coordinate_multiplication(o-p,urcorner q-p) ; fi
319%    addto @# also q shifted ((tp)-llcorner q) ;
320%  fi
321%enddef ;
322
323% This new version makes gdraw clip the result to the window defined with setrange
324
325vardef graph_picture_conversion@#(expr q, o)(text tp) =
326  save p ;
327  save do_clip, tp_clipped ; boolean do_clip ; do_clip := true ;
328  picture tp_clipped ; tp_clipped := nullpicture;
329  if stroked q :
330    path p ; p=pathpart q;
331    addto tp_clipped doublepath tp graph_with_pen_and_color(q) dashed dashpart q ;
332    %draw bbox tp_clipped withcolor red ;
333  elseif filled q :
334    path p ; p=pathpart q;
335    addto tp_clipped contour tp graph_with_pen_and_color(q) ;
336    %draw bbox tp_clipped withcolor green ;
337  else :
338    if (urcorner q<>llcorner q) : do_clip := false ; fi % Do not clip the axis labels;
339    interim truecorners := 0 ;
340    pair p ; p=llcorner q;
341    if urcorner q<>p : p := p + graph_coordinate_multiplication(o-p,urcorner q-p) ; fi
342    addto tp_clipped also q shifted ((tp)-llcorner q) ;
343    %draw bbox tp_clipped withcolor if do_clip : cyan else : blue fi ;
344  fi
345  if do_clip :
346    clip tp_clipped to origin--(xpart Z_.graph_dimensions,0)--Z_.graph_dimensions--
347                                 (0,ypart Z_.graph_dimensions)--cycle ;
348  fi
349  addto @# also tp_clipped ;
350enddef ;
351
352def graph_coordinate_multiplication(expr a,b) = (xpart a*xpart b, ypart a*ypart b) enddef ;
353
354vardef graph_clear_bounds@# =  numeric @#.low, @#.high ;  enddef;
355
356% Finalize anything drawn in the present coordinate system and set up a new
357% system as requested
358
359vardef setcoords(expr tx, ty) =
360  interim warningcheck :=0 ;
361  if length graph_current_graph>0 :
362    save s, S, T ;
363    graph_remap(s, S, T) ;
364    for q within graph_current_graph :
365      graph_picture_conversion.graph_finished_graph(q,-s,p shifted s transformed S transformed T) ;
366    endfor
367    graph_current_graph := graph_current_bb := nullpicture ;
368  fi
369  graph_clear_bounds.X_ ; graph_clear_bounds.Y_;
370  X_.graph_coordinate_type := tx ; Y_.graph_coordinate_type := ty;
371enddef ;
372
373% Set Z_.low and Z_.high to correspond to given range of user graph
374% coordinates.  The text argument should be a sequence of pairs and/or strings
375% with 4 components in all.
376
377vardef setrange(text t) =
378  interim warningcheck :=0 ;
379  save r_ ; r_=0;
380  string r_[]s ;
381  for x_=
382      for p_=t : if pair p_ : xpart p_, ypart fi p_, endfor :
383    r_[incr r_] if string x_ : s fi = x_ ;
384    if r_>2 :
385      graph_set_bounds if r_=3 : X_ else : Y_ fi (r_[r_-2] if unknown r_[r_-2] : s fi, x_) ;
386    fi
387    exitif r_=4 ;
388  endfor
389enddef ;
390
391% @# is X_ or Y_ ; l and h are numeric or string
392
393vardef graph_set_bounds@#(expr l, h) =
394  graph_clear_bounds@# ;
395  if @#graph_coordinate_type>0 :
396     @#low  = if unknown l :
397                whatever
398              else :
399                if abs @#graph_coordinate_type=logarithmic : graph_mlog fi if string l : scantokens fi l
400              fi ;
401     @#high = if unknown h :
402                whatever
403              else :
404                if abs @#graph_coordinate_type=logarithmic : graph_mlog fi if string h : scantokens fi h
405              fi ;
406  else :
407    -@#high = if unknown l :
408                whatever
409              else :
410                if abs @#graph_coordinate_type=logarithmic : graph_mlog fi if string l : scantokens fi l
411              fi ;
412    -@#low  = if unknown h :
413                whatever
414              else :
415                if abs @#graph_coordinate_type=logarithmic : graph_mlog fi if string h : scantokens fi h
416              fi ;
417  fi
418enddef ;
419
420%%%%%%%%%%%%%%%%%%%%%%%%% Converting path coordinates %%%%%%%%%%%%%%%%%%%%%%%%%
421
422% Find the result of scanning path p and using macros tx and ty to adjust the
423% x and y parts of each coordinate pair.  Boolean parameter c tells whether to
424% force the result to be polygonal.
425
426vardef graph_scan_path(expr p, c)(suffix tx, ty) =
427  if (str tx="") and (str ty="") :  p
428  else :
429    save r_ ; path r_;
430    r_ := graph_pair_adjust(point 0 of p, tx, ty)
431    if path p :
432      for t=1 upto length p :
433        if c : --
434        else : ..controls graph_pair_adjust(postcontrol(t-1) of p, tx, ty)
435          and graph_pair_adjust(precontrol t of p, tx, ty) ..
436        fi
437        graph_pair_adjust(point t of p, tx, ty)
438      endfor
439      if cycle p : &cycle fi
440    fi ;
441    if pair p : point 0 of fi r_
442  fi
443enddef ;
444
445vardef graph_pair_adjust(expr p)(suffix tx, ty) = (tx xpart p, ty ypart p) enddef ;
446
447% Convert path p from user graph coords to internal graph coords.
448
449vardef graph_convert_user_path_to_internal primary p =
450  interim warningcheck :=0 ;
451  if known p :
452    graph_scan_path(p,
453      (abs X_.graph_coordinate_type<>linear) or (abs Y_.graph_coordinate_type<>linear),
454      if abs X_.graph_coordinate_type=logarithmic : graph_mlog fi,
455      if abs Y_.graph_coordinate_type=logarithmic : graph_mlog fi)
456      transformed  (identity
457        if X_.graph_coordinate_type<0 : xscaled -1 fi
458        if Y_.graph_coordinate_type<0 : yscaled -1 fi)
459  fi
460enddef ;
461
462% Convert label location t_ from user graph coords to internal graph coords.
463% The label location should be a pair, or two numbers/strings.  If t_ is empty
464% or a single item of non-pair type, just return t_.  Unknown coordinates
465% produce unknown components in the result.
466
467vardef graph_label_convert_user_to_internal(text t_) =
468  save n_ ; n_=0;
469  interim warningcheck :=0 ;
470  if 0 for x_=t_ : +1 if pair x_ : +1 fi endfor <= 1 :
471    t_
472  else :
473    n_0 = n_1 = 0 ;
474    point 0 of graph_convert_user_path_to_internal (
475      for x_=
476        for y_=t_ : if pair y_ : xpart y_, ypart fi y_, endfor
477        0, 0 :
478        if known x_ : if string x_ : scantokens fi x_
479        else : hide(n_[n_] :=whatever) 0
480        fi
481        exitif incr n_=2 ;
482      ,endfor) + (n_0,n_1)
483  fi
484enddef ;
485
486%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Reading data files %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
487
488% Read a line from file f, extract whitespace-separated tokens ignoring any
489% initial "%", and return true if at least one token is found.  The tokens
490% are stored in @#1, @#2, .. with "" in the last @#[] entry.
491
492% String manipulation routines for MetaPost
493% It is harmless to input this file more than once.
494
495% vardef isdigit primary d =
496%     ("0"<=d)and(d<="9")
497% enddef ;
498
499% Number of initial characters of string s where `c <character>' is true
500
501vardef graph_cspan(expr s)(text c) =
502    0
503    for i=1 upto length s:
504        exitunless c substring (i-1,i) of s;
505        + 1
506    endfor
507enddef ;
508
509% String s is composed of items separated by white space.  Lop off the first
510% item and the surrounding white space and return just the item.
511
512vardef graph_loptok suffix s =
513    save t, k;
514    k = graph_cspan(s," ">=);
515    if k > 0 :
516        s := substring(k,infinity) of s ;
517    fi
518    k := graph_cspan(s," "<);
519    string t;
520    t = substring (0,k) of s;
521    s := substring (k,infinity) of s;
522    s := substring (graph_cspan(s," ">=),infinity) of s;
523    t
524enddef ;
525
526vardef graph_read_line@#(expr f) =
527  save n_, s_ ; string s_;
528  s_ = readfrom f ;
529  string @#[] ;
530  if s_<>EOF :
531    @#0 := s_ ;
532    @#1 := graph_loptok s_ ;
533    n_ = if @#1="%" : 0 else : 1 fi ;
534    forever :
535      @#[incr n_] := graph_loptok s_ ;
536      exitif @#[n_]="" ;
537    endfor
538    @#1<>""
539  else : false
540  fi
541enddef ;
542
543% Execute c for each line of data read from file f, and stop at the first
544% line with no data.  Commands c can use line number i and tokens $1, $2, ...
545% and j is the number of fields.
546
547% def gdata(expr f)(suffix $)(text c) =
548%   for i=1 upto largestmantissa :
549%     exitunless graph_read_line$(f) ;
550%     c
551%   endfor ;
552% enddef ;
553
554def gdata(expr f)(suffix $)(text c) =
555  for i=1 upto largestmantissa :
556    exitunless graph_read_line$(f) ;
557    c
558  endfor
559enddef ;
560
561% Read a path from file f. The path is terminated by blank line or EOF.
562
563vardef graph_readpath(expr f) =
564  interim warningcheck :=0 ;
565  save s ;
566  gdata(f, s, if i>1 :--fi
567      if s2="" : (            i, scantokens s1)
568      else :     (scantokens s1, scantokens s2) fi
569  )
570enddef ;
571
572% Append coordinates t to polygonal path @#.  The coordinates can be numerics,
573% strings, or a single pair.
574
575vardef augment@#(text t) =
576  interim warningcheck := 0 ;
577  if not path begingroup @# endgroup :
578    graph_error(begingroup @# endgroup, "Cannot augment--not a path") ;
579  else :
580    def graph_comma= hide(def graph_comma=,enddef) enddef ;
581    if known @# :  @# :=@#--  else :  @#=  fi
582    (for p=t :
583       graph_comma if string p : scantokens fi p
584     endfor) ;
585  fi
586enddef ;
587
588%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Drawing and filling %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
589
590% Unknown pair components are set to 0 because glabel and gdotlabel understand
591% unknown coordinates as `0 in absolute units'.
592
593vardef graph_unknown_pair_bbox(expr p) =
594  interim warningcheck:=0 ;
595  if known p : addto graph_current_bb doublepath p ;
596  else :
597    save x,y ;
598    z = llcorner graph_current_bb ;
599    if unknown xpart p : xpart p= else : x := fi 0 ;
600    if unknown ypart p : ypart p= else : y := fi 0 ;
601    addto graph_current_bb doublepath (p+z) ;
602  fi
603  graph_current_bb := image(fill llcorner graph_current_bb..urcorner graph_current_bb--cycle) ;
604enddef ;
605
606% Initiate a gdraw or gfill command.  This must be done before scanning the
607% argument, because that could invoke the `if known graph_plot_picture' test in a following
608% plot option .
609
610def graph_addto =
611  def graph_errorbar_text = enddef ;
612  color graph_foreground ;
613  path graph_last_path ;
614  graph_last_drawn := graph_plot_picture := nullpicture ; addto graph_last_drawn
615enddef;
616
617% Handle the part of a gdraw command that uses path or data file p.
618
619def graph_draw expr p =
620  if string p :    hide(graph_last_path := graph_readpath(p) ;)
621                   graph_convert_user_path_to_internal graph_last_path
622  elseif path p or pair p :
623                   hide(graph_last_path := p ;)
624                   graph_convert_user_path_to_internal p
625  else : graph_error(p,"gdraw argument should be a data file or a path")
626        origin
627  fi
628  withpen currentpen graph_withlist _op_
629enddef ;
630
631% Handle the part of a gdraw command that uses path or data file p.
632
633def graph_fill expr p =
634  if string p :    hide(graph_last_path := graph_readpath(p) --cycle ;)
635                   graph_convert_user_path_to_internal graph_last_path
636  elseif cycle p : hide(graph_last_path := p ;)
637                   graph_convert_user_path_to_internal p
638  else : graph_error(p,"gfill argument should be a data file or a cyclic path")
639        origin..cycle
640  fi graph_withlist _op_
641enddef ;
642
643def gdraw = graph_addto doublepath graph_draw enddef ;
644def gfill = graph_addto contour graph_fill enddef ;
645
646% This is used in graph_draw and graph_fill to allow postprocessing graph_last_drawn
647
648def graph_withlist text t_ = t_ ; graph_post_draw; enddef;
649
650def witherrorbars(text t) text options =
651    hide(
652        def graph_errorbar_text = t enddef ;
653        save pic ; picture pic ; pic := image(draw origin _op_ options ;) ;
654        if color colorpart pic : graph_foreground := colorpart pic ; fi
655    )
656    options
657enddef ;
658
659% new feature: graph_errorbars
660
661picture graph_errorbar_picture ; graph_errorbar_picture := image(draw (left--right) scaled .5 ;) ;
662%picture graph_xbar_picture ;    graph_xbar_picture := image(draw (down--up)    scaled .5 ;) ;
663%picture graph_ybar_picture ;    graph_ybar_picture := image(draw (left--right) scaled .5 ;) ;
664
665vardef graph_errorbars(text t) =
666    if known graph_last_path :
667        save n, p, q ; path p ; pair q ;
668        save pic ; picture pic[] ; pic0 := nullpicture ;
669        pic1 := if     known graph_xbar_picture :     graph_xbar_picture
670                elseif known graph_errorbar_picture : graph_errorbar_picture rotated 90
671                else :                                nullpicture fi ;
672        pic2 := if     known graph_ybar_picture :     graph_ybar_picture
673                elseif known graph_errorbar_picture : graph_errorbar_picture
674                else :                                nullpicture fi ;
675        if length pic1>0 :
676            pic1 := pic1 scaled graph_shapesize ;
677            setbounds pic1 to origin..cycle ;
678        fi
679        if length pic2>0 :
680            pic2 := pic2 scaled graph_shapesize ;
681            setbounds pic2 to origin..cycle ;
682        fi
683        for i=0 upto length graph_last_path :
684            clearxy ; z = point i of graph_last_path ;
685            n := 1 ;
686            for $=t :
687                if known $ :
688                    q := if path $ : if length $>i : point i of $ else : origin fi
689                       elseif pair $ : $ elseif numeric $ : ($,$) else : origin fi ;
690                    if q<>origin :
691                        p := graph_convert_user_path_to_internal ((
692                            if n=1 :
693                                (-xpart q,0)--(ypart q,0)
694                            else :
695                                (0,-xpart q)--(0,ypart q)
696                            fi ) shifted z) ;
697                        addto pic0 doublepath p ;
698                        if length pic[n]>0 :
699                           if ypart q<>0 :
700                               addto pic0 also pic[n] shifted point 1 of p ;
701                           fi
702                           if xpart q<>0 :
703                               addto pic0 also pic[n] rotated 180 shifted point 0 of p ;
704                           fi
705                        fi
706                    fi
707                fi
708                exitif incr n>3 ;
709            endfor
710        endfor
711        if length pic0>0 :
712            save bg, fg ; color bg, fg ;
713            bg := if known graph_background : graph_background else : background fi ;
714            fg := if known graph_foreground : graph_foreground else : black fi ;
715            addto graph_current_graph also pic0 withpen currentpen scaled 2  _op_ withcolor bg ;
716            addto graph_current_graph also pic0 withpen currentpen scaled .5 _op_ withcolor fg ;
717        fi
718    fi
719enddef ;
720
721% Set graph_plot_picture so the postprocessing step will plot picture p at each path knot.
722% Also select nullpen to suppress stroking.
723
724def plot expr p =
725  if known graph_plot_picture :
726    withpen nullpen
727    hide (graph_plot_picture := image(
728        if bounded p : for q within p : graph_addto_currentpicture q endfor    % Save memory
729        else : graph_addto_currentpicture p
730        fi graph_setbounds origin..cycle))
731  fi
732enddef ;
733
734% This hides a semicolon that could prematurely end graph_withlist's text argument
735
736def graph_addto_currentpicture primary p = addto currentpicture also p ; enddef;
737def graph_setbounds = setbounds currentpicture to enddef ;
738
739def gdrawarrow =    graph_number_of_arrowheads := 1 ; gdraw enddef;
740def gdrawdblarrow = graph_number_of_arrowheads := 2 ; gdraw enddef;
741
742% Post-process the filled or stroked picture graph_last_drawn as follows : (1) update
743% the bounding box information ; (2) transfer it to graph_current_graph unless the pen has
744% been set to nullpen to disable stroking ; (3) plot graph_plot_picture at each knot.
745
746vardef graph_post_draw =
747  save p ; path p ; p = pathpart graph_last_drawn ;
748  graph_unknown_pair_bbox(p) ;
749  if filled graph_last_drawn or not graph_is_null(penpart graph_last_drawn) :
750    addto graph_current_graph also graph_last_drawn ;
751  fi
752  graph_errorbars(graph_errorbar_text) ;
753  if length graph_plot_picture>0 :
754    for i=0 upto length p if cycle p : -1 fi :
755      addto graph_current_graph also graph_plot_picture shifted point i of p ;
756    endfor
757    picture graph_plot_picture ;
758  fi
759  if graph_number_of_arrowheads>0 :
760    graph_draw_arrowhead(p, graph_with_pen_and_color(graph_last_drawn)) ;
761    if graph_number_of_arrowheads>1 :
762      graph_draw_arrowhead(reverse p, graph_with_pen_and_color(graph_last_drawn)) ;
763    fi
764    graph_number_of_arrowheads := 0 ;
765  fi
766enddef ;
767
768vardef graph_is_null(expr p) = (urcorner p=origin) and (llcorner p=origin) enddef ;
769
770vardef graph_draw_arrowhead(expr p)(text w) = % Draw arrowhead for path p, with list w
771  %save r ; r := angle(precontrol infinity of p shifted -point infinity of p) ;
772  addto graph_current_graph also
773    image(fill arrowhead (graph_arrowhead_extent(precontrol infinity of p,point infinity of p)) w ;
774          draw arrowhead (graph_arrowhead_extent(precontrol infinity of p,point infinity of p)) w
775            undashed ;
776%if (r mod 90 <> 0) : % orientation can be wrong due to remapping
777%  draw textext("\tfxx " & decimal r) shifted point infinity of p withcolor blue ;
778%fi
779          graph_setbounds point infinity of p..cycle ;
780    ) ; % rotatedabout(point infinity of p,-r) ;
781enddef ;
782
783vardef graph_arrowhead_extent(expr p, q) =
784  if p<>q : (q - 100pt*unitvector(q-p)) -- fi
785  q
786enddef ;
787
788%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Drawing labels %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
789
790% Argument c is a drawing command that needs an additional argument p that gives
791% a location in internal graph coords.  Draw in graph_current_graph enclosed in a setbounds
792% path. Unknown components of p cause the setbounds path to have width or height 1 instead of 0.
793% Then graph_unknown_pair_bbox sets these components to 0 and graph_picture_conversion
794% suppresses subsequent repositioning.
795
796def graph_draw_label(expr p)(suffix $)(text c) =
797  save sdim_ ; pair sdim_;
798  sdim_ := (if unknown xpart p : 1+ fi 0, if unknown ypart p : 1+ fi 0) ;
799  graph_unknown_pair_bbox(p) ;
800  addto graph_current_graph also
801    image(c(p) ; graph_setbounds p--p+sdim_--cycle) _op_
802enddef ;
803
804% Stash the result drawing command c in the graph_label table using with list w and
805% an index based on angle mfun_laboff$.
806
807vardef graph_stash_label(suffix $)(text c) text w =
808  graph_label[1.5+angle mfun_laboff$ /90] = image(c(origin) w) ;
809enddef ;
810
811def graph_label_location primary p =
812  if pair p : graph_draw_label(p)
813  elseif numeric p : graph_draw_label(point p of pathpart graph_last_drawn)
814  else : graph_stash_label
815  fi
816enddef ;
817
818% Place label p at user graph coords t using with list w. (t is a time, a pair
819% or 2 numerics or strings).
820
821vardef glabel@#(expr p)(text t) text w  =
822  graph_label_location graph_label_convert_user_to_internal(t)  (@#,label@#(p)) w ; enddef;
823
824% Place label p at user graph coords t using with list w and draw a dot there.
825% (t is a time, a pair, or 2 numerics or strings).
826
827vardef gdotlabel@#(expr p)(text t) text w =
828  graph_label_location graph_label_convert_user_to_internal(t)  (@#,dotlabel@#(p)) w ; enddef;
829
830def OUT = enddef ;       % location text for outside labels
831
832%%%%%%%%%%%%%%%%%%%%%%%%%% Grid lines, ticks, etc. %%%%%%%%%%%%%%%%%%%%%%%%%%
833
834% Grid lines and tick marks are transformed versions of the templates below.
835% In the template paths, (0,0) is on the edge of the frame and inward is to
836% the right.
837
838path graph_template.tick, graph_template.itick, graph_template.otick, graph_template.grid ;
839graph_template.tick  = (-3.5bp,0)--(3.5bp,0) ;
840graph_template.itick = origin--(7bp,0) ;
841graph_template.otick = (-7bp,0)--origin ;
842graph_template.grid  = origin--(1,0) ;
843
844vardef  tick@#(expr f,u) text w =  graph_tick_label(@#,@,false,f,u,w) ;  enddef;
845
846vardef itick@#(expr f,u) text w =  graph_tick_label(@#,@,false,f,u,w) ;  enddef;
847
848vardef otick@#(expr f,u) text w =  graph_tick_label(@#,@,false,f,u,w) ;  enddef;
849
850vardef  grid@#(expr f,u) text w =  graph_tick_label(@#,@,true,f,u,w) ;  enddef;
851
852
853% Produce a tick or grid mark for label suffix $, graph_template suffix $$,
854% coordinate value u, and with list w.  Boolean c tells whether graph_template$$
855% needs scaling by X_.graph_dimensions or Y_.graph_dimensions,
856% and f gives a format string or a label picture.
857
858def graph_tick_label(suffix $,$$)(expr c, f, u)(text w) =
859  graph_draw_label(graph_label_convert_user_to_internal(graph_generate_label_position($,u)),,
860                   draw graph_gridline_picture$($$,c,f,u,w) shifted)
861enddef ;
862
863% Generate label positioning arguments appropriate for label suffix $ and
864% coordinate u.
865
866def graph_generate_label_position(suffix $)(expr u) =
867  if pair u : u elseif xpart mfun_laboff.$=0 : u,whatever else : whatever,u fi
868enddef ;
869
870% Generate a picture of a grid line labeled with coordinate value u, picture
871% or format string f, and with list w.  Suffix @# is bot, top, lft, or rt,
872% suffix $ identifies entries in the graph_template table, and boolean c tells
873% whether to scale graph_template$.
874
875vardef graph_gridline_picture@#(suffix $)(expr c, f, u)(text w) =
876  if unknown u : graph_error(u,"Label coordinate should be known") ; nullpicture
877  else :
878    save p ; path p;
879    interim warningcheck :=0 ;
880    graph_autogrid_needed :=false ;
881    p = graph_template$ zscaled -mfun_laboff@#
882        if c : graph_xyscale fi
883        shifted (((.5 + mfun_laboff@# dotprod (.5,.5)) * mfun_laboff@#) graph_xyscale) ;
884    image(draw p w ;
885        label@#(if string f : format(f,u) else : f fi, point 0 of p) w)
886  fi
887enddef ;
888
889def graph_xyscale = xscaled X_.graph_dimensions yscaled Y_.graph_dimensions enddef ;
890
891% Draw the frame or the part corresponding to label suffix @# using with list w.
892
893vardef frame@# text w =
894  graph_frame_needed :=false ;
895  picture p_ ;
896  p_ = image(draw
897    if str@#<>"" :  subpath round(angle mfun_laboff@#*graph_frame_pair_a+graph_frame_pair_b) of  fi
898    unitsquare graph_xyscale  w) ;
899  graph_draw_label((whatever,whatever),,draw p_ shifted) ;
900enddef ;
901
902pair graph_frame_pair_a ; graph_frame_pair_a=(1,1)/90; % unitsquare subpath is linear in label angle
903pair graph_frame_pair_b ; graph_frame_pair_b=(.75,2.25);
904
905%%%%%%%%%%%%%%%%%%%%%%%%%% Automatic grid selection %%%%%%%%%%%%%%%%%%%%%%%%%%
906
907string graph_log_marks[] ;           % marking options per decade for logarithmic scales
908string graph_lin_marks ;             % mark spacing options per decade for linear scales
909string graph_exp_marks ;             % exponent spacing options for logarithmic scales
910newinternal graph_minimum_number_of_marks, graph_log_minimum ;
911graph_minimum_number_of_marks := 4 ; % minimum number marks generated by auto.x or auto.y
912graph_log_minimum := mlog 3 ;        % revert to uniform marks when largest/smallest < this
913
914def Gfor(text t) = for i=t endfor enddef ; % to shorten the mark templates below
915
916graph_log_marks[1]="1,2,5" ;
917graph_log_marks[2]="1,1.5,2,3,4,5,7" ;
918graph_log_marks[3]="1Gfor(6upto10 :,i/5)Gfor(5upto10 :,i/2)Gfor(6upto9 :,i)" ;
919graph_log_marks[4]="1Gfor(11upto20 :,i/10)Gfor(11upto25 :,i/5)Gfor(11upto19 :,i/2)" ;
920graph_log_marks[5]="1Gfor(21upto40 :,i/20)Gfor(21upto50 :,i/10)Gfor(26upto49 :,i/5)" ;
921graph_lin_marks="10,5,2" ;           % start with 10 and go down; a final `,1' is appended
922graph_exp_marks="20,10,5,2,1" ;
923
924Ten_to0 =     1 ;
925Ten_to1 =    10 ;
926Ten_to2 =   100 ;
927Ten_to3 =  1000 ;
928Ten_to4 = 10000 ;
929
930% Determine the X_ or Y_ bounds on the range to be covered by automatic grid
931% marks.  Suffix @# is X_ or Y_.  The result is logarithmic or linear to specify the
932% type of grid spacing to use.  Bounds are returned in variables local to
933% begingraph..endgraph : pairs graph_modified_lower and graph_modified_higher
934%  are upper and lower bounds in
935% `modified exponential form'.  In modified exponential form, (x,y) means
936% (x/1000)*10^y, where 1000<=abs x<10000.
937
938vardef graph_bounds@# =
939  interim warningcheck :=0 ;
940  save l, h ;
941  graph_set_default_bounds ;
942  if @#graph_coordinate_type>0 : (l,h) else : -(h,l) fi = (@#low, @#high) ;
943  if abs @#graph_coordinate_type=logarithmic :
944    graph_modified_lower  := graph_Meform(l)+graph_modified_bias ;
945    graph_modified_higher := graph_Meform(h)+graph_modified_bias ;
946    if h-l >= graph_log_minimum : logarithmic else : linear fi
947  else :
948    graph_modified_lower  := graph_Feform(l)+graph_modified_bias ;
949    graph_modified_higher := graph_Feform(h)+graph_modified_bias ;
950    linear
951  fi
952enddef ;
953
954pair graph_modified_bias ; graph_modified_bias=(0,3);
955pair graph_modified_lower, graph_modified_higher ;
956
957% Scan graph_log_marks[k] and evaluate tokens t for each m where l<=m<=h.
958
959def graph_scan_marks(expr k, l, h)(text t) =
960  for m=scantokens graph_log_marks[k] :
961    exitif m>h ;
962    if m>=l : t fi
963  endfor
964enddef ;
965
966% Scan graph_log_marks[k] and evaluate tokens t for each m and e where m*10^e belongs
967% between l and h (inclusive), where both l and h are in modified exponent form.
968
969def graph_scan_mark(expr k, l, h)(text t) =
970  for e=ypart l upto ypart h :
971    graph_scan_marks(k, if e>ypart l : 1 else : xpart l/1000 fi,
972        if e<ypart h : 10 else : xpart h/1000 fi,  t)
973  endfor
974enddef ;
975
976% Select a k for which graph_scan_mark(k,...) gives enough marks.
977
978vardef graph_select_mark =
979  save k ;
980  k = 0 ;
981  forever :
982    exitif unknown graph_log_marks[k+1] ;
983    exitif 0 graph_scan_mark(incr k, graph_modified_lower, graph_modified_higher, +1)
984           >= graph_minimum_number_of_marks ;
985  endfor
986  k
987enddef ;
988
989% Try to select an exponent spacing from graph_exp_marks.  If successful, set @# and
990% return true
991
992vardef graph_select_exponent_mark@# =
993  numeric @# ;
994  for e=scantokens graph_exp_marks :
995    @# = e ;
996    exitif floor(ypart graph_modified_higher/e) -
997           floor(graph_modified_exponent_ypart(graph_modified_lower)/e)
998           >= graph_minimum_number_of_marks ;
999    numeric @# ;
1000  endfor
1001  known @#
1002enddef ;
1003
1004vardef graph_modified_exponent_ypart(expr p) = ypart p  if xpart p=1000 : -1 fi  enddef ;
1005
1006% Compute the mark spacing d between xpart graph_modified_lower and xpart graph_modified_higher.
1007
1008vardef graph_tick_mark_spacing =
1009  interim warningcheck :=0 ;
1010  save m, n, d ;
1011  m = graph_minimum_number_of_marks ;
1012  n = 1 for i=1 upto
1013                (mlog(xpart graph_modified_higher-xpart graph_modified_lower) - mlog m)/mlogten :
1014        *10 endfor ;
1015  if n<=1000 :
1016    for x=scantokens graph_lin_marks :
1017      d = n*x ;
1018      exitif 0 graph_generate_numbers(d,+1)>=m ;
1019      numeric d ;
1020    endfor
1021  fi
1022  if known d : d else : n fi
1023enddef ;
1024
1025def graph_generate_numbers(expr d)(text t) =
1026  for m = d*ceiling(xpart graph_modified_lower/d) step d until xpart graph_modified_higher :
1027    t
1028  endfor
1029enddef ;
1030
1031% Evaluate tokens t for exponents e in multiples of d in the range determined
1032% by graph_modified_lower and graph_modified_higher.
1033
1034def graph_generate_exponents(expr d)(text t) =
1035  for e = d*floor(graph_modified_exponent_ypart(graph_modified_lower)/d+1)
1036      step d until d*floor(ypart graph_modified_higher/d) :  t
1037  endfor
1038enddef ;
1039
1040% Adjust graph_modified_lower and graph_modified_higher so their exponent parts match
1041% and they are in true exponent form ((x,y) means x*10^y).  Return the new exponent.
1042
1043vardef graph_match_exponents =
1044  interim warningcheck := 0 ;
1045  save e ;
1046  e+3 = if graph_modified_lower=graph_modified_bias : ypart graph_modified_higher
1047        elseif graph_modified_higher=graph_modified_bias : ypart graph_modified_lower
1048        else : max(ypart graph_modified_lower, ypart graph_modified_higher) fi ;
1049  forsuffixes $=graph_modified_lower, graph_modified_higher :
1050    $ := (xpart $ for i=ypart $ upto e+2 : /(10) endfor, e) ;
1051  endfor
1052  e
1053enddef ;
1054
1055% Assume e is an integer and either m=0 or 1<=abs(m)<10000.  Find m*(10^e)
1056% and represent the result as a string if its absolute value would be at least
1057% 4096 or less than .1.  It is OK to return 0 as a string or a numeric.
1058
1059vardef graph_factor_and_exponent_to_string(expr m, e) =
1060  if (e>3)or(e<-4) :
1061    decimal m & "e" & decimal e
1062  elseif e>=0 :
1063    if abs m<infinity/Ten_to[e] :
1064          m*Ten_to[e]
1065    else : decimal m & "e" & decimal e
1066    fi
1067  else :
1068    save x ; x=m/Ten_to[-e];
1069    if abs x>=.1 : x else : decimal m & "e" & decimal e fi
1070  fi
1071enddef ;
1072
1073def auto suffix $ =
1074  hide(def graph_comma= hide(def graph_comma=,enddef) enddef)
1075  if graph_bounds.graph_suffix($)=logarithmic :
1076    if graph_select_exponent_mark.graph_exponent :
1077      graph_generate_exponents(graph_exponent,
1078        graph_comma graph_factor_and_exponent_to_string(1,e))
1079    else :
1080      graph_scan_mark(graph_select_mark, graph_modified_lower, graph_modified_higher,
1081        graph_comma graph_factor_and_exponent_to_string(m,e))
1082    fi
1083  else :
1084    hide(graph_exponent :=graph_match_exponents)
1085    graph_generate_numbers(graph_tick_mark_spacing,
1086      graph_comma graph_factor_and_exponent_to_string(m,graph_exponent))
1087  fi
1088enddef ;
1089
1090string Autoform ; Autoform = "%g";
1091
1092%vardef autogrid(suffix tx, ty) text w =
1093%  graph_autogrid_needed :=false ;
1094%  if str tx<>"" : for x=auto.x : tx(Autoform,x) w ; endfor fi
1095%  if str ty<>"" : for y=auto.y : ty(Autoform,y) w ; endfor fi
1096%enddef ;
1097
1098% We redefine autogrid, adding the possibility of differing X and Y
1099% formats.
1100
1101% string Autoform_X ; Autoform_X := "@.0e" ;
1102% string Autoform_Y ; Autoform_Y := "@.0e" ;
1103
1104vardef autogrid(suffix tx, ty) text w =
1105    graph_autogrid_needed := false ;
1106    if str tx <> "" :
1107        for x=auto.x :
1108            tx (
1109                if string Autoform_X :
1110                    if Autoform_X <> "" :
1111                        Autoform_X
1112                    else :
1113                        Autoform
1114                    fi
1115                else :
1116                    Autoform
1117                fi,
1118                x
1119            ) w ;
1120        endfor
1121    fi
1122    if str ty <> "" :
1123        for y=auto.y :
1124            ty (
1125                if string Autoform_Y :
1126                    if Autoform_Y <> "" :
1127                        Autoform_Y
1128                    else :
1129                        Autoform
1130                    fi
1131                else :
1132                    Autoform
1133                fi,
1134                y
1135            ) w ;
1136        endfor
1137    fi
1138enddef ;
1139
1140%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% endgraph %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1141
1142def endgraph =
1143  if graph_autogrid_needed : autogrid(otick.bot, otick.lft) ; fi
1144  if    graph_frame_needed : frame ; fi
1145  setcoords(linear,linear) ;
1146  interim truecorners :=1 ;
1147  for  b=bbox graph_finished_graph :
1148    setbounds graph_finished_graph to b ;
1149    for i=0 step .5 until 3.5 :
1150      if known graph_label[i] :
1151        addto graph_finished_graph also graph_label[i] shifted point i of b ;
1152      fi
1153    endfor
1154  endfor
1155  graph_finished_graph
1156  endgroup
1157enddef ;
1158
1159%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1160
1161% We format in luatex (using \mathematics{}) ...
1162% we could pass via variables and save escaping as that is inefficient
1163
1164if known contextlmtxmode :
1165% already defined
1166elseif unknown context_mlib :
1167
1168    vardef escaped_format(expr s) =
1169        "" for n=0 upto length(s) : &
1170            if ASCII substring (n,n+1) of s = 37 :
1171                "@"
1172            else :
1173                substring (n,n+1) of s
1174            fi
1175        endfor
1176    enddef ;
1177
1178    vardef strfmt(expr f, x) = % maybe use mfun_ namespace
1179        "\MPgraphformat{" & escaped_format(f) & "}{" & mfun_tagged_string(x) & "}"
1180    enddef ;
1181
1182    vardef varfmt(expr f, x) = % maybe use mfun_ namespace
1183        "\MPformatted{" & escaped_format(f) & "}{" & mfun_tagged_string(x) & "}"
1184    enddef ;
1185
1186    vardef format   (expr f, x) = textext(strfmt(f,x)) enddef ;
1187    vardef formatted(expr f, x) = textext(varfmt(f,x)) enddef ;
1188
1189fi ;
1190
1191%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1192
1193% A couple of extensions :
1194
1195% Define a function plotsymbol() returning a picture : 10 different shapes,
1196% unfilled outline, interior filled with different shades of the background.
1197% This allows overlapping points on a plot to be more distinguishable.
1198
1199vardef graph_shapesize = (.33BodyFontSize) enddef ;
1200
1201path graph_shape[] ; % (internal) symbol path
1202
1203graph_shape[0] := (0,0) ;                      % point
1204graph_shape[1] := fullcircle ;                 % circle
1205graph_shape[2] := (up -- down) scaled .5 ;     % vertical bar
1206
1207for i = 3 upto 9 :                          % polygons
1208    graph_shape[i] :=
1209        for j = 0 upto i-1 :
1210            (up scaled .5) rotated (360j/i) --
1211        endfor cycle ;
1212endfor
1213
1214graph_shape[12] := graph_shape[2] rotated +90 ;   % horizontal line
1215graph_shape[22] := graph_shape[2] rotated +45 ;   % backslash
1216graph_shape[32] := graph_shape[2] rotated -45 ;   % slash
1217graph_shape[13] := graph_shape[3] rotated 180 ;   % down  triangle
1218graph_shape[23] := graph_shape[3] rotated -90 ;   % right triangle
1219graph_shape[33] := graph_shape[3] rotated +90 ;   % left  triangle
1220graph_shape[14] := graph_shape[4] rotated +45 ;   % square
1221graph_shape[15] := graph_shape[5] rotated 180 ;   % down pentagon
1222graph_shape[16] := graph_shape[6] rotated +90 ;   % turned hexagon
1223graph_shape[17] := graph_shape[7] rotated 180 ;
1224graph_shape[18] := graph_shape[8] rotated +22.5 ;
1225
1226numeric l ;
1227
1228for j = 5 upto 9 :
1229    l := length(graph_shape[j]) ;
1230    pair p[] ;
1231    for i = 0 upto l :
1232        p[i] = whatever [point i             of graph_shape[j],
1233                         point (i+2   mod l) of graph_shape[j]] ;
1234        p[i] = whatever [point (i+1   mod l) of graph_shape[j],
1235                         point (i+l-1 mod l) of graph_shape[j]] ;
1236    endfor
1237    graph_shape[20+j] := for i = 0 upto l : point i of graph_shape[j]--p[i]--endfor cycle ;
1238endfor
1239
1240path s    ; s := graph_shape[4] ;
1241path q    ; q := s scaled .25 ;
1242numeric l ; l := length(s) ;
1243
1244pair p[] ;
1245
1246graph_shape[24] := for i = 0 upto l-1 :
1247     hide(
1248         p[i]   = whatever [point i   of s, point (i+1   mod l) of s] ;
1249         p[i]   = whatever [point i   of q, point (i-1+l mod l) of q] ;
1250         p[i+l] = whatever [point i   of s, point (i+1   mod l) of s] ;
1251         p[i+l] = whatever [point i+1 of q, point (i+2   mod l) of q] ;
1252     )
1253     point i of q -- p[i] -- p[i+l] --
1254endfor cycle ;
1255
1256graph_shape[34] := graph_shape[24] rotated 45 ;
1257
1258% usage : gdraw p plot plotsymbol( 1,1) ;  % a filled circle
1259% usage : gdraw p plot plotsymbol(14,0) ;  % a square
1260% usage : gdraw p plot plotsymbol( 4,.5) ; % a 50% filled diamond
1261
1262def          stars(expr f) = plotsymbol(25,f) enddef ; % a 5-point star
1263def         points(expr f) = plotsymbol( 0,f) enddef ;
1264def        circles(expr f) = plotsymbol( 1,f) enddef ;
1265def        crosses(expr f) = plotsymbol(34,f) enddef ;
1266def        squares(expr f) = plotsymbol(14,f) enddef ;
1267def       diamonds(expr f) = plotsymbol( 4,f) enddef ; % a turned square
1268def    uptriangles(expr f) = plotsymbol( 3,f) enddef ;
1269def  downtriangles(expr f) = plotsymbol(13,f) enddef ;
1270def  lefttriangles(expr f) = plotsymbol(33,f) enddef ;
1271def righttriangles(expr f) = plotsymbol(23,f) enddef ;
1272
1273% f (fill) is color, numeric or boolean, otherwise background.
1274def plotsymbol(expr n, f) text t =
1275    if known graph_shape[n] :
1276        image(
1277            save bg, fg ; color bg, fg ;
1278            bg := if known graph_background : graph_background else : background fi ;
1279            save pic ; picture pic ; pic := image(draw origin _op_ t ;) ;
1280            if color colorpart pic : graph_foreground := colorpart pic ; fi
1281            fg := if known graph_foreground : graph_foreground else : black fi ;
1282            save p ; path p ; p = graph_shape[n] scaled graph_shapesize ;
1283            draw p withcolor bg withpen currentpen scaled 2 ; % halo
1284            if cycle p :
1285                fill p withcolor
1286                    if known f :
1287                        if color f :
1288                            f
1289                        elseif numeric f :
1290                            f[bg,fg]
1291                        elseif boolean f and f :
1292                            fg
1293                        else
1294                            bg
1295                        fi
1296                    else :
1297                        bg
1298                    fi ;
1299            fi
1300            draw p _op_ t ;
1301        )
1302    else :
1303        nullpicture
1304    fi
1305    t
1306enddef ;
1307
1308% standard resistance color code: rainbow sequence (from /usr/share/X11/rgb.txt)
1309color resistance_color[] ;                      string resistance_name[] ;
1310resistance_color0 = (0,0,0) ;                   resistance_name0 = "black" ;
1311resistance_color1 = (165/255,42/255,42/255) ;   resistance_name1 = "brown" ;
1312resistance_color2 = (1,0,0) ;                   resistance_name2 = "red" ;
1313resistance_color3 = (1,165/255,0) ;             resistance_name3 = "orange" ;
1314resistance_color4 = (1,1,0) ;                   resistance_name4 = "yellow" ;
1315resistance_color5 = (0,1,0) ;                   resistance_name5 = "green" ;
1316resistance_color6 = (0,0,1) ;                   resistance_name6 = "blue" ;
1317resistance_color7 = (148/255,0,211/255) ;       resistance_name7 = "darkviolet" ;
1318resistance_color8 = (190/255,190/255,190/255) ; resistance_name8 = "gray" ;
1319resistance_color9 = (1,1,1) ;                   resistance_name9 = "white" ;
1320
1321%def rainbow(expr f) =
1322%    ((abs(5f) mod 5) + 2 - floor((abs(5f) mod 5) + 2))
1323%       [resistance_color[  floor((abs(5f) mod 5) + 2)],
1324%        resistance_color[ceiling((abs(5f) mod 5) + 2)]]
1325%enddef ;
1326def rainbow(expr f) =
1327    hide(numeric n_ ; n_ = (abs(5f) mod 5) + 2 ;)
1328    (n_-floor(n_))[resistance_color[floor n_],resistance_color[ceiling n_]]
1329enddef ;
1330
1331% The following extensions are not specific to graph and could be moved to metafun...
1332
1333% sort a path. Efficient en memory use, not so efficient in sorting long paths...
1334
1335vardef sortpath (suffix $) (text t) = % t can be "xpart", "ypart", "length", "angle", ...
1336    if path $ :
1337        if length $ > 0 :
1338            save n, k ; n := length $ ;
1339            for i=0 upto n :
1340                k := i ;
1341                for j=i+1 upto n :
1342                    if t (point j of $) < t (point k of $) :
1343                        k := j ;
1344                    fi
1345                endfor
1346                if k>i :
1347                    $ := if i>0 : subpath (0,i-1) of $ -- fi
1348                         point k of $ --
1349                         subpath (i,k-1) of $
1350                         if k<n : -- subpath (k+1,n) of $ fi
1351                    ;
1352                fi
1353            endfor
1354        fi
1355    fi
1356enddef ;
1357
1358% convert a polygon path to a smooth path (useful, e.g. as a guide to the eye)
1359
1360def smoothpath (suffix $) =
1361    if path $ :
1362    (for i=0 upto length $ :
1363        if i>0 : .. fi
1364        (point i of $)
1365     endfor )
1366    fi
1367enddef ;
1368
1369% return a path of a function func(x) with abscissa running from f to t over n intervals
1370
1371def makefunctionpath (expr f, t, n) (text func) =
1372    (for x=f step ((t-f)/(abs n)) until t :
1373        if x<>f : -- fi
1374        (x, func)
1375    endfor )
1376enddef ;
1377
1378% shift a path, point by point
1379%
1380% example :
1381%
1382%   p1 := addtopath(p0,(.1normaldeviate,.1normaldeviate)) ;
1383
1384vardef addtopath (suffix p) (text t) =
1385    if path p :
1386    (for i=0 upto length p :
1387        if i>0 : -- fi
1388        hide(clearxy ; z = point i of p ;) z shifted t
1389    endfor)
1390    fi
1391enddef ;
1392
1393% return a new path of a function func(z) using the same abscissa as an existing path
1394
1395vardef functionpath (suffix p) (text func) =
1396    (for i=0 upto length p :
1397        if i>0 : .. fi
1398        (hide(x := xpart(point i of p))x,func) %(hide(clearxy ; z = point i of p)x,func)
1399     endfor )
1400enddef ;
1401
1402% least-squares "fit" to a polynomial
1403%
1404% example :
1405%
1406%   path p[] ;
1407%   numeric a[] ; a0 := 1 ; a1 := .1 ; a2 := .01 ; a3 := .001 ; a4 := 0.0001 ;
1408%   p0 := makefunctionpath(0,5,10,polynomial_function(a,4,x)) ;
1409%   p1 := addtopath(p0,(0,.001normaldeviate)) ;
1410%   gdraw p0 ;
1411%   gdraw p1 plot plotsymbol(1,.5) ;
1412%
1413%   numeric b[] ;
1414%   polynomial_fit(p1, b, 4, 1) ;
1415%   gdraw functionpath(p1,polynomial_function(b,4,x)) ;
1416%
1417%   numeric c[] ;
1418%   linear_fit(p1, c, 1) ;
1419%   gdraw functionpath(p1,linear_function(c,x)) dashed evenly ;
1420
1421% a polynomial function :
1422%
1423% y = a0 + a1 * x + a2 * x^2 + ... + a[n] * x^n
1424
1425vardef polynomial_function (suffix $) (expr n, x) =
1426    (for j=0 upto n : + $[j]*(x**j) endfor) % no ;
1427enddef ;
1428
1429% find the determinant of a (n+1)*(n+1) matrix ; indices run from 0 to n
1430
1431vardef det (suffix $) (expr n) =
1432    hide(
1433        numeric determinant ; determinant := 1 ;
1434        save jj ; numeric jj ;
1435        for k=0 upto n :
1436            if $[k][k]=0 :
1437                jj := -1 ;
1438                for j=0 upto n :
1439                    if $[k][j]<>0 :
1440                        jj := j ;
1441                        exitif true ;
1442                    fi
1443                endfor
1444                if jj<0 :
1445                    determinant := 0 ;
1446                    exitif true ;
1447                fi
1448                for j=k upto n : % interchange the columns
1449                    temp := $[j][jj] ;
1450                    $[j][jj] := $[j][k] ;
1451                    $[j][k] := temp ;
1452                endfor
1453                determinant = -determinant ;
1454            fi
1455            exitif determinant=0 ;
1456            determinant := determinant * $[k][k] ;
1457            if k<n : % subtract row k from lower rows to get a diagonal matrix
1458                for j=k+1 upto n :
1459                    for i=k+1 upto n :
1460                        $[j][i] := $[j][i]-$[j][k]*$[k][i]/$[k][k] ;
1461                    endfor
1462                endfor
1463            fi
1464        endfor ;
1465    )
1466    determinant % no ;
1467enddef ;
1468
1469numeric fit_chi_squared ;
1470
1471% least-squares fit of a polynomial $ of order n to a path p (unweighted)
1472%
1473% reference : P. R. Bevington, "Data Reduction and Error Analysis for the Physical
1474% Sciences", McGraw-Hill, New York 1969.
1475
1476vardef polynomial_fit (suffix p, $) (expr n) (text t) =
1477    if not path p :
1478        graph_error(p, "Cannot fit--not a path") ;
1479    elseif length p < n :
1480        graph_error(p, "Cannot fit--not enough points") ;
1481    else :
1482        fit_chi_squared := 0 ;
1483        % calculate sums of the data
1484        save sumx, sumy ; numeric sumx[], sumy[] ;
1485        save w ; numeric w ;
1486        for i=0 upto 2n :
1487            sumx[i] := 0 ;
1488        endfor
1489        for i=0 upto n :
1490            sumy[i] := 0 ;
1491        endfor
1492        for i=0 upto length p :
1493            clearxy ; z = point i of p ;
1494            w := 1 ; % weight
1495            if known t :
1496                if  numeric t :
1497                    w := 1 if t<>0 : /(abs t) fi ;
1498                elseif pair t :
1499                    if t<>origin :
1500                        w := 1/(abs t) ;
1501                    fi
1502                elseif path t :
1503                    if length t>= i:
1504                        if point i of t<>origin :
1505                            w := 1/(abs point i of t) ;
1506                        fi
1507                    else :
1508                        w := 0 ;
1509                    fi ;
1510                fi
1511            fi
1512            x1 := w ;
1513            for j=0 upto 2n :
1514                sumx[j] := sumx[j] + x1 ;
1515                x1 := x1 * x ;
1516            endfor
1517            y1 := y * w ;
1518            for j=0 upto n :
1519               sumy[j] := sumy[j] + y1 ;
1520               y1 := y1 * x ;
1521            endfor
1522            fit_chi_squared := fit_chi_squared + y*y*w ;
1523        endfor
1524        % construct matrices and calculate the polynomial coefficients
1525        save m ; numeric m[][] ;
1526        for j=0 upto n :
1527            for k=0 upto n :
1528                m[j][k] := sumx[j+k] ;
1529            endfor
1530        endfor
1531        save delta ; numeric delta ;
1532        delta := det(m,n) ; % this destroys the matrix m[][], which is OK
1533        if delta = 0 :
1534            fit_chi_squared := 0 ;
1535            for j=0 upto n :
1536                $[j] := 0 ;
1537            endfor
1538        else :
1539            for i=0 upto n :
1540                for j=0 upto n :
1541                    for k=0 upto n :
1542                        m[j][k] := sumx[j+k] ;
1543                    endfor
1544                    m[j][i] := sumy[j] ;
1545                endfor
1546                $[i] := det(m,n) / delta ; % matrix m[][] gets destroyed...
1547            endfor
1548            for j=0 upto n :
1549                fit_chi_squared := fit_chi_squared - 2sumy[j]*$[j] ;
1550                for k=0 upto n :
1551                    fit_chi_squared := fit_chi_squared + $[j]*$[k]*sumx[j+k] ;
1552                endfor
1553            endfor
1554            % normalize by the number of degrees of freedom
1555            fit_chi_squared := fit_chi_squared / (length(p) - n) ; % length(p)+1-(n+1)
1556        fi
1557    fi
1558enddef ;
1559
1560% y = a0 + a1 * x
1561%
1562% of course a line is just a polynomial of order 1
1563
1564vardef linear_function (suffix $) (expr x) = polynomial_function($,1,x)   enddef ;
1565vardef linear_fit   (suffix p, $) (text t) = polynomial_fit(p, $, 1, t) ; enddef ;
1566
1567% and a constant is polynomial of order 0
1568
1569vardef constant_function (suffix $) (expr x) = polynomial_function($,0,x)   enddef ;
1570vardef constant_fit   (suffix p, $) (text t) = polynomial_fit(p, $, 0, t) ; enddef ;
1571
1572% y = a1 * exp(a0*x)
1573%
1574% exp and ln defined in metafun
1575
1576vardef exponential_function (suffix $) (expr x) = $1*exp($0*x) enddef ;
1577
1578% since we take a log, this only works for positive ordinates
1579
1580vardef exponential_fit (suffix p, $) (text t) =
1581    save a ; numeric a[] ;
1582    save q ; path q[] ; % fit to the log of the ordinate
1583    for i=0 upto length p :
1584        clearxy ; z = point i of p ;
1585        if y>0 :
1586            augment.q0(x,ln(y)) ;
1587            augment.q1(
1588                if known t :
1589                    if numeric t : (0,ln(t))
1590                    elseif pair t : (xpart t,ln(ypart t))
1591                    elseif path t :
1592                        if length t>=i :
1593                            hide(z1 = point i of t;)
1594                            (x1,ln(y1))
1595                        else :
1596                            origin
1597                        fi
1598                    fi
1599                else :
1600                    (0,1)
1601                fi ) ;
1602        fi
1603    endfor
1604    linear_fit(q0,a,q1) ;
1605    save e ; e := exp(sqrt(fit_chi_squared)) ;
1606    fit_chi_squared := e * e ;
1607    $0 := a1 ;
1608    $1 := exp(a0) ;
1609enddef ;
1610
1611% y = a1 * x**a0
1612
1613vardef power_law_function (suffix $) (expr x) = $1*(x**$0) enddef ;
1614
1615% since we take logs, this only works for positive abscissae and ordinates
1616
1617vardef power_law_fit (suffix p, $) (text t) =
1618    save a ; numeric a[] ;
1619    save q ; path q[] ; % fit to the logs of the abscissae and ordinates
1620    for i=0 upto length p :
1621        clearxy ; z = point i of p ;
1622        if (x>0) and (y>0) :
1623            augment.q0(ln(x),ln(y)) ;
1624            augment.q1(
1625                if known t :
1626                    if numeric t : (0,ln(t))
1627                    elseif pair t : (ln(xpart t),ln(ypart t))
1628                    elseif path t :
1629                        if length t>=i :
1630                            hide(z1 = point i of t)
1631                            (ln(x1),ln(y1))
1632                        else :
1633                            origin
1634                        fi
1635                    fi
1636                else :
1637                    (0,1)
1638                fi ) ;
1639        fi
1640    endfor
1641    linear_fit(q0,a,q1) ;
1642    save e ; e := exp(sqrt(fit_chi_squared)) ;
1643    fit_chi_squared := e * e ;
1644    $0 := a1 ;
1645    $1 := exp(a0) ;
1646enddef ;
1647
1648% gaussian : y = a2 * exp(-ln(2)*((x-a0)/a1)^2)
1649%
1650% a1 is the hwhm ; sigma := a1/sqrt(2ln(2)) or a1/1.17741
1651
1652newinternal lntwo ; lntwo := ln(2) ; % brrr, why not inline it
1653
1654vardef gaussian_function (suffix $) (expr x) =
1655    if $1 = 0 :
1656        if x = $0 : $2 else : 0 fi
1657    else :
1658        $2 * exp(-lntwo*(((x-$0)/$1)**2))
1659    fi
1660    if known $3 :
1661        + $3
1662    fi
1663enddef ;
1664
1665% since we take a log, this only works for positive ordinates
1666
1667vardef gaussian_fit (suffix p, $) (text t) =
1668    save a ; numeric a[] ;
1669    save q ; path q[] ; % fit to the log of the ordinate
1670    for i=0 upto length p :
1671        clearxy ; z = point i of p ;
1672        if y>0 :
1673            augment.q0(x,ln(y)) ;
1674            augment.q1(
1675                if known t :
1676                    if numeric t : (0,ln(t))
1677                    elseif pair t : (xpart t,ln(ypart t))
1678                    elseif path t :
1679                        if length t>=i :
1680                            hide(z1 = point i of t)
1681                            (x1,ln(y1))
1682                        else :
1683                            origin
1684                        fi
1685                    fi
1686                else :
1687                    (0,1)
1688                fi ) ;
1689        fi
1690    endfor
1691    polynomial_fit(q0,a,2,q1) ;
1692    save e ; e := exp(sqrt(fit_chi_squared)) ;
1693    fit_chi_squared := e * e ;
1694    $1 := sqrt(-lntwo/a2) ;
1695    $0 := -.5a1/a2 ;
1696    $2 := exp(a0-.25*a1*a1/a2) ;
1697    $3 := 0 ; % polynomial_fit will NOT work with a non-zero background!
1698enddef ;
1699
1700% lorentzian: y = a2 / (1 + ((x - a0)/a1)^2)
1701
1702vardef lorentzian_function (suffix $) (expr x) =
1703    if $1 = 0 :
1704        if x = $0 : $2 else : 0 fi
1705    else :
1706        $2 / (1 + ((x - $0)/$1)**2)
1707    fi
1708    if known $3 :
1709        + $3
1710    fi
1711enddef ;
1712
1713vardef lorentzian_fit (suffix p, $) (text t) =
1714    save a ; numeric a[] ;
1715    save q ; path q ; % fit to the inverse of the ordinate
1716    for i=0 upto length p :
1717        if ypart(point i of p)<>0 :
1718            augment.q(xpart(point i of p), 1/ypart(point i of p)) ;
1719        fi
1720    endfor
1721    polynomial_fit(q,a,2,if t <> 0 : 1/(t) else : 0 fi) ;
1722    fit_chi_squared := 1/fit_chi_squared ;
1723    $0 := -.5a1/a2 ;
1724    $2 := 1/(a0-.25a1*a1/a2) ;
1725    $1 := sqrt((a0-.25a1*a1/a2)/a2) ;
1726    $3 := 0 ; % polynomial_fit will NOT work with a non-zero background!
1727enddef ;
1728