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 14
if
known
context_grap
:
endinput
;
fi
;
15 16
boolean
context_grap
;
context_grap
:
=
true
;
17 18
% Below is a modified graph.mp
19 20
message
(
"
using number system
"
&
numbersystem
&
"
with precision
"
&
decimal
numberprecision
)
;
21 22
if
numbersystem
<
>
"
double
"
:
23
errmessage
"
The graph macros require the double precision number system.
"
;
24
endinput
;
25
fi
26 27
if
known
contextlmtxmode
:
28
def
_op_
=
base_draw_options
enddef
;
29
fi
;
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 79
newinternal
Mzero
;
Mzero
:
=
-16384
;
% Anything at least this small is treated as zero
80
newinternal
mlogten
;
mlogten
:
=
mlog
(
10
)
;
81
newinternal
largestmantissa
;
largestmantissa
:
=
2
*
*
52
;
% internal double warningcheck
82
newinternal
singleinfinity
;
singleinfinity
:
=
2
*
*
128
;
83
newinternal
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 88
vardef
graph_mlog
primary
x
=
89
if
unknown
x
:
whatever
90
elseif
x
=
0
:
Mzero
91
else
:
mlog
(
abs
x
)
fi
92
enddef
;
93 94
vardef
graph_exp
primary
x
=
95
if
unknown
x
:
whatever
96
elseif
x
<
=
Mzero
:
0
97
else
:
mexp
(
x
)
fi
98
enddef
;
99 100
% and add the following for utility/completeness
101
% (replacing the definitions in mp-tool.mpiv).
102 103
vardef
logten
primary
x
=
104
if
unknown
x
:
whatever
105
elseif
x
=
0
:
Mzero
106
else
:
mlog
(
abs
x
)
/
mlog
(
10
)
fi
107
enddef
;
108 109
if
known
contextlmtxmode
:
110
% already defined
111
else
:
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 125
fi
126 127
vardef
powten
primary
x
=
128
if
unknown
x
:
whatever
129
elseif
x
<
=
Mzero
:
0
130
else
:
10
*
*
x
fi
131
enddef
;
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 136
vardef
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
143
enddef
;
144 145
% Modified from above.
146 147
vardef
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
155
enddef
;
156 157
vardef
graph_error
(
expr
x
,
s
)
=
158
interim
showstopping
:
=
0
;
159
show
x
;
errmessage
s
;
160
enddef
;
161 162
%%%%%%%%%%%%%%%%%%%%%%%% Data structures, begingraph %%%%%%%%%%%%%%%%%%%%%%%%
163 164
vardef
Z_
@#
=
(
X_
@#
,
Y_
@#
)
enddef
;
% used in place of plain.mp's z convention
165 166
def
graph_suffix
(
suffix
$
)
=
% convert from x or y to X_ or Y_
167
if
str
$
=
"
x
"
:
X_
else
:
Y_
fi
168
enddef
;
169 170
% New :
171 172
save
graph_background
;
color
graph_background
;
% if defined, fill the frame.
173 174
def
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
215
enddef
;
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 232
newinternal
logarithmic
,
linear
;
% coordinate system codes
233
logarithmic
:
=
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 249
vardef
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
255
enddef
;
256 257
pair
graph_margin_pair
.
low
,
graph_margin_pair
.
high
;
258
graph_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 265
vardef
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
;
274
enddef
;
275 276
graph_margin_fraction
.
low
=
-.07
;
% bbox fraction for default range start
277
graph_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 295
def
graph_with_pen_and_color
(
expr
q
)
=
296
withproperties
q
297
enddef
;
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 325
vardef
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
;
350
enddef
;
351 352
def
graph_coordinate_multiplication
(
expr
a
,
b
)
=
(
xpart
a
*
xpart
b
,
ypart
a
*
ypart
b
)
enddef
;
353 354
vardef
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 359
vardef
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
;
371
enddef
;
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 377
vardef
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
389
enddef
;
390 391
% @# is X_ or Y_ ; l and h are numeric or string
392 393
vardef
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
418
enddef
;
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 426
vardef
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
443
enddef
;
444 445
vardef
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 449
vardef
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
460
enddef
;
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 467
vardef
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
484
enddef
;
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 501
vardef
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
507
enddef
;
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 512
vardef
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
524
enddef
;
525 526
vardef
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
541
enddef
;
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 554
def
gdata
(
expr
f
)
(
suffix
$
)
(
text
c
)
=
555
for
i
=
1
upto
largestmantissa
:
556
exitunless
graph_read_line
$
(
f
)
;
557
c
558
endfor
559
enddef
;
560 561
% Read a path from file f. The path is terminated by blank line or EOF.
562 563
vardef
graph_readpath
(
expr
f
)
=
564
interim
warningcheck
:
=
0
;
565
save
s
;
566
gdata
(
f
,
s
,
if
i
>
1
:
-
-
fi
567
if
s
2
=
"
"
:
(
i
,
scantokens
s
1
)
568
else
:
(
scantokens
s
1
,
scantokens
s
2
)
fi
569
)
570
enddef
;
571 572
% Append coordinates t to polygonal path @#. The coordinates can be numerics,
573
% strings, or a single pair.
574 575
vardef
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
586
enddef
;
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 593
vardef
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
)
;
604
enddef
;
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 610
def
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
615
enddef
;
616 617
% Handle the part of a gdraw command that uses path or data file p.
618 619
def
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_
629
enddef
;
630 631
% Handle the part of a gdraw command that uses path or data file p.
632 633
def
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_
641
enddef
;
642 643
def
gdraw
=
graph_addto
doublepath
graph_draw
enddef
;
644
def
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 648
def
graph_withlist
text
t_
=
t_
;
graph_post_draw
;
enddef
;
649 650
def
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
657
enddef
;
658 659
% new feature: graph_errorbars
660 661
picture
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 665
vardef
graph_errorbars
(
text
t
)
=
666
if
known
graph_last_path
:
667
save
n
,
p
,
q
;
path
p
;
pair
q
;
668
save
pic
;
picture
pic
[
]
;
pic
0
:
=
nullpicture
;
669
pic
1
:
=
if
known
graph_xbar_picture
:
graph_xbar_picture
670
elseif
known
graph_errorbar_picture
:
graph_errorbar_picture
rotated
90
671
else
:
nullpicture
fi
;
672
pic
2
:
=
if
known
graph_ybar_picture
:
graph_ybar_picture
673
elseif
known
graph_errorbar_picture
:
graph_errorbar_picture
674
else
:
nullpicture
fi
;
675
if
length
pic
1
>
0
:
676
pic
1
:
=
pic
1
scaled
graph_shapesize
;
677
setbounds
pic
1
to
origin
.
.
cycle
;
678
fi
679
if
length
pic
2
>
0
:
680
pic
2
:
=
pic
2
scaled
graph_shapesize
;
681
setbounds
pic
2
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
pic
0
doublepath
p
;
698
if
length
pic
[
n
]
>
0
:
699
if
ypart
q
<
>
0
:
700
addto
pic
0
also
pic
[
n
]
shifted
point
1
of
p
;
701
fi
702
if
xpart
q
<
>
0
:
703
addto
pic
0
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
pic
0
>
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
pic
0
withpen
currentpen
scaled
2
_op_
withcolor
bg
;
716
addto
graph_current_graph
also
pic
0
withpen
currentpen
scaled
.5
_op_
withcolor
fg
;
717
fi
718
fi
719
enddef
;
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 724
def
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
732
enddef
;
733 734
% This hides a semicolon that could prematurely end graph_withlist's text argument
735 736
def
graph_addto_currentpicture
primary
p
=
addto
currentpicture
also
p
;
enddef
;
737
def
graph_setbounds
=
setbounds
currentpicture
to
enddef
;
738 739
def
gdrawarrow
=
graph_number_of_arrowheads
:
=
1
;
gdraw
enddef
;
740
def
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 746
vardef
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
766
enddef
;
767 768
vardef
graph_is_null
(
expr
p
)
=
(
urcorner
p
=
origin
)
and
(
llcorner
p
=
origin
)
enddef
;
769 770
vardef
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) ;
781
enddef
;
782 783
vardef
graph_arrowhead_extent
(
expr
p
,
q
)
=
784
if
p
<
>
q
:
(
q
-
100
pt
*
unitvector
(
q
-
p
)
)
--
fi
785
q
786
enddef
;
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 796
def
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_
802
enddef
;
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 807
vardef
graph_stash_label
(
suffix
$
)
(
text
c
)
text
w
=
808
graph_label
[
1.5
+
angle
mfun_laboff
$
/
90
]
=
image
(
c
(
origin
)
w
)
;
809
enddef
;
810 811
def
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
816
enddef
;
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 821
vardef
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 827
vardef
gdotlabel
@#
(
expr
p
)
(
text
t
)
text
w
=
828
graph_label_location
graph_label_convert_user_to_internal
(
t
)
(
@#
,
dotlabel
@#
(
p
)
)
w
;
enddef
;
829 830
def
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 838
path
graph_template
.
tick
,
graph_template
.
itick
,
graph_template
.
otick
,
graph_template
.
grid
;
839
graph_template
.
tick
=
(
-3.5
bp
,
0
)
--
(
3.5
bp
,
0
)
;
840
graph_template
.
itick
=
origin
--
(
7
bp
,
0
)
;
841
graph_template
.
otick
=
(
-7
bp
,
0
)
-
-
origin
;
842
graph_template
.
grid
=
origin
--
(
1
,
0
)
;
843 844
vardef
tick
@#
(
expr
f
,
u
)
text
w
=
graph_tick_label
(
@#
,
@
,
false
,
f
,
u
,
w
)
;
enddef
;
845 846
vardef
itick
@#
(
expr
f
,
u
)
text
w
=
graph_tick_label
(
@#
,
@
,
false
,
f
,
u
,
w
)
;
enddef
;
847 848
vardef
otick
@#
(
expr
f
,
u
)
text
w
=
graph_tick_label
(
@#
,
@
,
false
,
f
,
u
,
w
)
;
enddef
;
849 850
vardef
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 858
def
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
)
861
enddef
;
862 863
% Generate label positioning arguments appropriate for label suffix $ and
864
% coordinate u.
865 866
def
graph_generate_label_position
(
suffix
$
)
(
expr
u
)
=
867
if
pair
u
:
u
elseif
xpart
mfun_laboff
.
$
=
0
:
u
,
whatever
else
:
whatever
,
u
fi
868
enddef
;
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 875
vardef
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
887
enddef
;
888 889
def
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 893
vardef
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
)
;
900
enddef
;
901 902
pair
graph_frame_pair_a
;
graph_frame_pair_a
=
(
1
,
1
)
/
90
;
% unitsquare subpath is linear in label angle
903
pair
graph_frame_pair_b
;
graph_frame_pair_b
=
(
.75
,
2.25
)
;
904 905
%%%%%%%%%%%%%%%%%%%%%%%%%% Automatic grid selection %%%%%%%%%%%%%%%%%%%%%%%%%%
906 907
string
graph_log_marks
[
]
;
% marking options per decade for logarithmic scales
908
string
graph_lin_marks
;
% mark spacing options per decade for linear scales
909
string
graph_exp_marks
;
% exponent spacing options for logarithmic scales
910
newinternal
graph_minimum_number_of_marks
,
graph_log_minimum
;
911
graph_minimum_number_of_marks
:
=
4
;
% minimum number marks generated by auto.x or auto.y
912
graph_log_minimum
:
=
mlog
3
;
% revert to uniform marks when largest/smallest < this
913 914
def
Gfor
(
text
t
)
=
for
i
=
t
endfor
enddef
;
% to shorten the mark templates below
915 916
graph_log_marks
[
1
]
=
"
1,2,5
"
;
917
graph_log_marks
[
2
]
=
"
1,1.5,2,3,4,5,7
"
;
918
graph_log_marks
[
3
]
=
"
1Gfor(6upto10 :,i/5)Gfor(5upto10 :,i/2)Gfor(6upto9 :,i)
"
;
919
graph_log_marks
[
4
]
=
"
1Gfor(11upto20 :,i/10)Gfor(11upto25 :,i/5)Gfor(11upto19 :,i/2)
"
;
920
graph_log_marks
[
5
]
=
"
1Gfor(21upto40 :,i/20)Gfor(21upto50 :,i/10)Gfor(26upto49 :,i/5)
"
;
921
graph_lin_marks
=
"
10,5,2
"
;
% start with 10 and go down; a final `,1' is appended
922
graph_exp_marks
=
"
20,10,5,2,1
"
;
923 924
Ten_to
0
=
1
;
925
Ten_to
1
=
10
;
926
Ten_to
2
=
100
;
927
Ten_to
3
=
1000
;
928
Ten_to
4
=
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 938
vardef
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
952
enddef
;
953 954
pair
graph_modified_bias
;
graph_modified_bias
=
(
0
,
3
)
;
955
pair
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 959
def
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
964
enddef
;
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 969
def
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
974
enddef
;
975 976
% Select a k for which graph_scan_mark(k,...) gives enough marks.
977 978
vardef
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
987
enddef
;
988 989
% Try to select an exponent spacing from graph_exp_marks. If successful, set @# and
990
% return true
991 992
vardef
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
@#
1002
enddef
;
1003 1004
vardef
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 1008
vardef
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
1023
enddef
;
1024 1025
def
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
1029
enddef
;
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 1034
def
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
1038
enddef
;
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 1043
vardef
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
1053
enddef
;
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 1059
vardef
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
1071
enddef
;
1072 1073
def
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
1088
enddef
;
1089 1090
string
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 1104
vardef
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
1138
enddef
;
1139 1140
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% endgraph %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1141 1142
def
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
1157
enddef
;
1158 1159
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1160 1161
% We format in luatex (using \mathematics{}) ...
1162
% we could pass via variables and save escaping as that is inefficient
1163 1164
if
known
contextlmtxmode
:
1165
% already defined
1166
elseif
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 1189
fi
;
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 1199
vardef
graph_shapesize
=
(
.33
BodyFontSize
)
enddef
;
1200 1201
path
graph_shape
[
]
;
% (internal) symbol path
1202 1203
graph_shape
[
0
]
:
=
(
0
,
0
)
;
% point
1204
graph_shape
[
1
]
:
=
fullcircle
;
% circle
1205
graph_shape
[
2
]
:
=
(
up
--
down
)
scaled
.5
;
% vertical bar
1206 1207
for
i
=
3
upto
9
:
% polygons
1208
graph_shape
[
i
]
:
=
1209
for
j
=
0
upto
i
-1
:
1210
(
up
scaled
.5
)
rotated
(
360
j
/
i
)
--
1211
endfor
cycle
;
1212
endfor
1213 1214
graph_shape
[
12
]
:
=
graph_shape
[
2
]
rotated
+90
;
% horizontal line
1215
graph_shape
[
22
]
:
=
graph_shape
[
2
]
rotated
+45
;
% backslash
1216
graph_shape
[
32
]
:
=
graph_shape
[
2
]
rotated
-45
;
% slash
1217
graph_shape
[
13
]
:
=
graph_shape
[
3
]
rotated
180
;
% down triangle
1218
graph_shape
[
23
]
:
=
graph_shape
[
3
]
rotated
-90
;
% right triangle
1219
graph_shape
[
33
]
:
=
graph_shape
[
3
]
rotated
+90
;
% left triangle
1220
graph_shape
[
14
]
:
=
graph_shape
[
4
]
rotated
+45
;
% square
1221
graph_shape
[
15
]
:
=
graph_shape
[
5
]
rotated
180
;
% down pentagon
1222
graph_shape
[
16
]
:
=
graph_shape
[
6
]
rotated
+90
;
% turned hexagon
1223
graph_shape
[
17
]
:
=
graph_shape
[
7
]
rotated
180
;
1224
graph_shape
[
18
]
:
=
graph_shape
[
8
]
rotated
+22.5
;
1225 1226
numeric
l
;
1227 1228
for
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
;
1238
endfor
1239 1240
path
s
;
s
:
=
graph_shape
[
4
]
;
1241
path
q
;
q
:
=
s
scaled
.25
;
1242
numeric
l
;
l
:
=
length
(
s
)
;
1243 1244
pair
p
[
]
;
1245 1246
graph_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
]
--
1254
endfor
cycle
;
1255 1256
graph_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 1262
def
stars
(
expr
f
)
=
plotsymbol
(
25
,
f
)
enddef
;
% a 5-point star
1263
def
points
(
expr
f
)
=
plotsymbol
(
0
,
f
)
enddef
;
1264
def
circles
(
expr
f
)
=
plotsymbol
(
1
,
f
)
enddef
;
1265
def
crosses
(
expr
f
)
=
plotsymbol
(
34
,
f
)
enddef
;
1266
def
squares
(
expr
f
)
=
plotsymbol
(
14
,
f
)
enddef
;
1267
def
diamonds
(
expr
f
)
=
plotsymbol
(
4
,
f
)
enddef
;
% a turned square
1268
def
uptriangles
(
expr
f
)
=
plotsymbol
(
3
,
f
)
enddef
;
1269
def
downtriangles
(
expr
f
)
=
plotsymbol
(
13
,
f
)
enddef
;
1270
def
lefttriangles
(
expr
f
)
=
plotsymbol
(
33
,
f
)
enddef
;
1271
def
righttriangles
(
expr
f
)
=
plotsymbol
(
23
,
f
)
enddef
;
1272 1273
% f (fill) is color, numeric or boolean, otherwise background.
1274
def
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
1306
enddef
;
1307 1308
% standard resistance color code: rainbow sequence (from /usr/share/X11/rgb.txt)
1309
color
resistance_color
[
]
;
string
resistance_name
[
]
;
1310
resistance_color
0
=
(
0
,
0
,
0
)
;
resistance_name
0
=
"
black
"
;
1311
resistance_color
1
=
(
165
/
255
,
42
/
255
,
42
/
255
)
;
resistance_name
1
=
"
brown
"
;
1312
resistance_color
2
=
(
1
,
0
,
0
)
;
resistance_name
2
=
"
red
"
;
1313
resistance_color
3
=
(
1
,
165
/
255
,
0
)
;
resistance_name
3
=
"
orange
"
;
1314
resistance_color
4
=
(
1
,
1
,
0
)
;
resistance_name
4
=
"
yellow
"
;
1315
resistance_color
5
=
(
0
,
1
,
0
)
;
resistance_name
5
=
"
green
"
;
1316
resistance_color
6
=
(
0
,
0
,
1
)
;
resistance_name
6
=
"
blue
"
;
1317
resistance_color
7
=
(
148
/
255
,
0
,
211
/
255
)
;
resistance_name
7
=
"
darkviolet
"
;
1318
resistance_color
8
=
(
190
/
255
,
190
/
255
,
190
/
255
)
;
resistance_name
8
=
"
gray
"
;
1319
resistance_color
9
=
(
1
,
1
,
1
)
;
resistance_name
9
=
"
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 ;
1326
def
rainbow
(
expr
f
)
=
1327
hide
(
numeric
n_
;
n_
=
(
abs
(
5
f
)
mod
5
)
+
2
;
)
1328
(
n_
-
floor
(
n_
)
)
[
resistance_color
[
floor
n_
]
,
resistance_color
[
ceiling
n_
]
]
1329
enddef
;
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 1335
vardef
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
1356
enddef
;
1357 1358
% convert a polygon path to a smooth path (useful, e.g. as a guide to the eye)
1359 1360
def
smoothpath
(
suffix
$
)
=
1361
if
path
$
:
1362
(
for
i
=
0
upto
length
$
:
1363
if
i
>
0
:
..
fi
1364
(
point
i
of
$
)
1365
endfor
)
1366
fi
1367
enddef
;
1368 1369
% return a path of a function func(x) with abscissa running from f to t over n intervals
1370 1371
def
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
)
1376
enddef
;
1377 1378
% shift a path, point by point
1379
%
1380
% example :
1381
%
1382
% p1 := addtopath(p0,(.1normaldeviate,.1normaldeviate)) ;
1383 1384
vardef
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
1391
enddef
;
1392 1393
% return a new path of a function func(z) using the same abscissa as an existing path
1394 1395
vardef
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
)
1400
enddef
;
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 1425
vardef
polynomial_function
(
suffix
$
)
(
expr
n
,
x
)
=
1426
(
for
j
=
0
upto
n
:
+
$
[
j
]
*
(
x
*
*
j
)
endfor
)
% no ;
1427
enddef
;
1428 1429
% find the determinant of a (n+1)*(n+1) matrix ; indices run from 0 to n
1430 1431
vardef
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 ;
1467
enddef
;
1468 1469
numeric
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 1476
vardef
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
2
n
:
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
x
1
:
=
w
;
1513
for
j
=
0
upto
2
n
:
1514
sumx
[
j
]
:
=
sumx
[
j
]
+
x
1
;
1515
x
1
:
=
x
1
*
x
;
1516
endfor
1517
y
1
:
=
y
*
w
;
1518
for
j
=
0
upto
n
:
1519
sumy
[
j
]
:
=
sumy
[
j
]
+
y
1
;
1520
y
1
:
=
y
1
*
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
-
2
sumy
[
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
1558
enddef
;
1559 1560
% y = a0 + a1 * x
1561
%
1562
% of course a line is just a polynomial of order 1
1563 1564
vardef
linear_function
(
suffix
$
)
(
expr
x
)
=
polynomial_function
(
$
,
1
,
x
)
enddef
;
1565
vardef
linear_fit
(
suffix
p
,
$
)
(
text
t
)
=
polynomial_fit
(
p
,
$
,
1
,
t
)
;
enddef
;
1566 1567
% and a constant is polynomial of order 0
1568 1569
vardef
constant_function
(
suffix
$
)
(
expr
x
)
=
polynomial_function
(
$
,
0
,
x
)
enddef
;
1570
vardef
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 1576
vardef
exponential_function
(
suffix
$
)
(
expr
x
)
=
$
1
*
exp
(
$
0
*
x
)
enddef
;
1577 1578
% since we take a log, this only works for positive ordinates
1579 1580
vardef
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
.
q
0
(
x
,
ln
(
y
)
)
;
1587
augment
.
q
1
(
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
(
z
1
=
point
i
of
t
;
)
1594
(
x
1
,
ln
(
y
1
)
)
1595
else
:
1596
origin
1597
fi
1598
fi
1599
else
:
1600
(
0
,
1
)
1601
fi
)
;
1602
fi
1603
endfor
1604
linear_fit
(
q
0
,
a
,
q
1
)
;
1605
save
e
;
e
:
=
exp
(
sqrt
(
fit_chi_squared
)
)
;
1606
fit_chi_squared
:
=
e
*
e
;
1607
$
0
:
=
a
1
;
1608
$
1
:
=
exp
(
a
0
)
;
1609
enddef
;
1610 1611
% y = a1 * x**a0
1612 1613
vardef
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 1617
vardef
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
.
q
0
(
ln
(
x
)
,
ln
(
y
)
)
;
1624
augment
.
q
1
(
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
(
z
1
=
point
i
of
t
)
1631
(
ln
(
x
1
)
,
ln
(
y
1
)
)
1632
else
:
1633
origin
1634
fi
1635
fi
1636
else
:
1637
(
0
,
1
)
1638
fi
)
;
1639
fi
1640
endfor
1641
linear_fit
(
q
0
,
a
,
q
1
)
;
1642
save
e
;
e
:
=
exp
(
sqrt
(
fit_chi_squared
)
)
;
1643
fit_chi_squared
:
=
e
*
e
;
1644
$
0
:
=
a
1
;
1645
$
1
:
=
exp
(
a
0
)
;
1646
enddef
;
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 1652
newinternal
lntwo
;
lntwo
:
=
ln
(
2
)
;
% brrr, why not inline it
1653 1654
vardef
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
1663
enddef
;
1664 1665
% since we take a log, this only works for positive ordinates
1666 1667
vardef
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
.
q
0
(
x
,
ln
(
y
)
)
;
1674
augment
.
q
1
(
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
(
z
1
=
point
i
of
t
)
1681
(
x
1
,
ln
(
y
1
)
)
1682
else
:
1683
origin
1684
fi
1685
fi
1686
else
:
1687
(
0
,
1
)
1688
fi
)
;
1689
fi
1690
endfor
1691
polynomial_fit
(
q
0
,
a
,
2
,
q
1
)
;
1692
save
e
;
e
:
=
exp
(
sqrt
(
fit_chi_squared
)
)
;
1693
fit_chi_squared
:
=
e
*
e
;
1694
$
1
:
=
sqrt
(
-
lntwo
/
a
2
)
;
1695
$
0
:
=
-.5
a
1
/
a
2
;
1696
$
2
:
=
exp
(
a
0
-.25
*
a
1
*
a
1
/
a
2
)
;
1697
$
3
:
=
0
;
% polynomial_fit will NOT work with a non-zero background!
1698
enddef
;
1699 1700
% lorentzian: y = a2 / (1 + ((x - a0)/a1)^2)
1701 1702
vardef
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
1711
enddef
;
1712 1713
vardef
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
:
=
-.5
a
1
/
a
2
;
1724
$
2
:
=
1
/
(
a
0
-.25
a
1
*
a
1
/
a
2
)
;
1725
$
1
:
=
sqrt
(
(
a
0
-.25
a
1
*
a
1
/
a
2
)
/
a
2
)
;
1726
$
3
:
=
0
;
% polynomial_fit will NOT work with a non-zero background!
1727
enddef
;
1728