summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/metapost/context/base/mp-grap.mpiv
diff options
context:
space:
mode:
Diffstat (limited to 'Master/texmf-dist/metapost/context/base/mp-grap.mpiv')
-rw-r--r--Master/texmf-dist/metapost/context/base/mp-grap.mpiv602
1 files changed, 406 insertions, 196 deletions
diff --git a/Master/texmf-dist/metapost/context/base/mp-grap.mpiv b/Master/texmf-dist/metapost/context/base/mp-grap.mpiv
index a6a4e21e8c9..6d69c0b1ed0 100644
--- a/Master/texmf-dist/metapost/context/base/mp-grap.mpiv
+++ b/Master/texmf-dist/metapost/context/base/mp-grap.mpiv
@@ -41,6 +41,7 @@ fi
% OUT loc value for labels relative to whole graph
% gdata(file,s,text) read coords from file ; evaluate t w/ tokens s[]
% auto.<x or y> default x or y tick locations (for interation)
+% tick.<bot|top|..>(fmt,u) draw centered tick from given side at u w/ format
% itick.<bot|top|..>(fmt,u) draw inward tick from given side at u w/ format
% otick.<bot|top|..>(fmt,u) draw outward tick at coord u ; label format fmt
% grid.<bot|top|..>(fmt,u) draw grid line at u with given side labeled
@@ -51,11 +52,11 @@ fi
% endgraph end of graph--the result is a picture
% option `plot <picture>' draws picture at each path knot, turns off pen
-% Gtemplate.<tickcmd> template paths for tick marks and grid lines
+% graph_template.<tickcmd> template paths for tick marks and grid lines
% graph_margin_fraction.low,
% graph_margin_fraction.high fractions determining margins when no setrange
-% Glmarks[], Gumarks, Gemarks loop text strings used by auto.<x or y>
-% Gmarks, Gminlog numeric parameters used by auto.<x or y>
+% graph_log_marks[], graph_lin_marks, graph_exp_marks loop text strings used by auto.<x or y>
+% graph_minimum_number_of_marks, graph_log_minimum numeric parameters used by auto.<x or y>
% Autoform is the format string used by autogrid
% Autoform_X, Autoform_Y if defined, are used instead
@@ -64,40 +65,64 @@ fi
% with `graph_'
% Depends on :
+
input string.mp
% Private version of a few marith macros, fixed for double math...
-newinternal mlogten ; mlogten := mlog(10) ;
-newinternal doubleinfinity ; doubleinfinity := 2**1024 ;
-% Note that we get arithmetic overflows if we set to -doubleinfinity below.
+
+newinternal Mzero ; Mzero := -16384; % Anything at least this small is treated as zero
+newinternal mlogten ; mlogten := mlog(10) ;
+newinternal largestmantissa ; largestmantissa := 2**52 ; % internal double warningcheck
+newinternal singleinfinity ; singleinfinity := 2**128 ;
+newinternal doubleinfinity ; doubleinfinity := 2**1024 ;
+Mzero := -largestmantissa ; % Note that we get arithmetic overflows if we set to -doubleinfinity
% Safely convert a number to mlog form, trapping zero.
+
vardef graph_mlog primary x =
- if unknown x: whatever elseif x=0: -.5doubleinfinity else: mlog(abs x) fi
+ if unknown x: whatever
+ elseif x=0: Mzero
+ else: mlog(abs x) fi
enddef ;
+
vardef graph_exp primary x =
- if unknown x: whatever else: mexp(x) fi
+ if unknown x: whatever
+ elseif x<=Mzero: 0
+ else: mexp(x) fi
enddef ;
% and add the following for utility/completeness
% (replacing the definitions in mp-tool.mpiv).
+
vardef logten primary x =
- if unknown x: whatever elseif x=0: -.5doubleinfinity else: mlog(abs x)/mlog(10) fi
+ if unknown x: whatever
+ elseif x=0: Mzero
+ else: mlog(abs x)/mlog(10) fi
enddef ;
+
vardef ln primary x =
- if unknown x: whatever elseif x=0: -.5doubleinfinity else: mlog(abs x)/256 fi
+ if unknown x: whatever
+ elseif x=0: Mzero
+ else: mlog(abs x)/256 fi
enddef ;
+
vardef exp primary x =
- if unknown x: whatever else: (mexp 256)**x fi
+ if unknown x: whatever
+ elseif x<= Mzero: 0
+ else: (mexp 256)**x fi
enddef ;
+
vardef powten primary x =
- if unknown x: whatever else: 10**x fi
+ if unknown x: whatever
+ elseif x<= Mzero: 0
+ else: 10**x fi
enddef ;
% Convert x from mlog form into a pair whose xpart gives a mantissa and whose
% ypart gives a power of ten.
+
vardef graph_Meform(expr x) =
- if x<=-doubleinfinity : origin
+ if x<=Mzero : origin
else :
save e, m ; e=floor(x/mlogten)-3; m := mexp(x-e*mlogten) ;
if abs m<1000 : m := m*10 ; e := e-1 ; elseif abs m>=10000 : m := m/10 ; e := e+1 ; fi
@@ -106,6 +131,7 @@ vardef graph_Meform(expr x) =
enddef ;
% Modified from above.
+
vardef graph_Feform(expr x) =
interim warningcheck :=0 ;
if x=0 : origin
@@ -118,7 +144,7 @@ enddef ;
vardef graph_error(expr x,s) =
interim showstopping :=0 ;
- show x ; errmessage s;
+ show x ; errmessage s ;
enddef ;
%%%%%%%%%%%%%%%%%%%%%%%% Data structures, begingraph %%%%%%%%%%%%%%%%%%%%%%%%
@@ -130,6 +156,7 @@ def graph_suffix(suffix $) = % convert from x or y to X_ or Y_
enddef ;
% New :
+
save graph_background ; color graph_background ; % if defined, fill the frame.
save graph_close_file ; boolean graph_close_file ; graph_close_file = false ;
@@ -153,8 +180,12 @@ def begingraph(expr w, h) =
save graph_last_drawn ;
picture graph_last_drawn ; % result of last gdraw or gfill
graph_last_drawn = nullpicture ;
+ save graph_last_path ;
+ path graph_last_path ; % last gdraw or gfill path in data coordinates.
save graph_plot_picture ;
picture graph_plot_picture ; % a picture from the `plot' option known when plot allowed
+ save graph_foreground ;
+ color graph_foreground ; % drawing color, if set.
save graph_label ;
picture graph_label[] ; % labels to place around the whole graph when it is done
save graph_autogrid_needed ;
@@ -180,17 +211,20 @@ enddef ;
% user to alter the behavior of these macros.
% Not very modifiable : log, linear,
% graph_frame_pair_a, graph_frame_pair_b, graph_margin_pair
-% Modifiable : Gtemplate.suffix, Glmarks[], Gumarks, Gemarks, Gmarks,
-% Gminlog, Autoform
+% Modifiable : graph_template.suffix,
+% graph_log_marks[], graph_lin_marks, graph_exp_marks,
+% graph_minimum_number_of_marks,
+% graph_log_minimum, Autoform
newinternal log, linear ; % coordinate system codes
log :=1 ; linear :=2;
+
% note that mp-tool.mpiv defines log as log10.
%%%%%%%%%%%%%%%%%%%%%% Coordinates : setcoords, setrange %%%%%%%%%%%%%%%%%%%%%%
-% Graph-related usr input is `user graph coordinates' as specified by arguments
+% Graph-related user input is `user graph coordinates' as specified by arguments
% to setcoords.
% `Internal graph coordinates' are used for graph_current_graph, graph_current_bb, Z_.low, Z_.high.
% Their meaning depends on the appropriate component of Z_.graph_coordinate_type :
@@ -207,14 +241,15 @@ vardef graph_set_default_bounds = % Set default Z_.low, Z_.high
graph_margin_pair$ ;
endfor
enddef ;
+
pair graph_margin_pair.low, graph_margin_pair.high ;
graph_margin_pair.high = -graph_margin_pair.low = (.00002,.00002) ;
+% Set $, $$, $$$ so that shifting by $ then transforming by $$ and then $$$ maps
+% the essential bounding box of graph_current_graph into (0,0)..Z_.graph_dimensions.
+% The `essential bounding box' is either what Z_.low and Z_.high imply
+% or the result of ignoring pen widths in graph_current_graph.
-% Set $, $$, $$$ so that shifting by $ then transforming by $$ and then $$$
-% maps the essential bounding box of graph_current_graph into (0,0)..Z_.graph_dimensions. The
-% `essential bounding box' is either what Z_.low and Z_.high imply or the
-% result of ignoring pen widths in graph_current_graph.
vardef graph_remap(suffix $,$$,$$$) =
save p_ ;
graph_set_default_bounds ;
@@ -225,10 +260,10 @@ vardef graph_remap(suffix $,$$,$$$) =
(Z_.high+$) transformed $$ = p_ ;
p_ transformed $$$ = Z_.graph_dimensions ;
enddef ;
+
graph_margin_fraction.low=-.07 ; % bbox fraction for default range start
graph_margin_fraction.high=1.07 ; % bbox fraction for default range stop
-
def graph_with_pen_and_color(expr q) =
withpen penpart q withcolor
if colormodel q=1 :
@@ -242,12 +277,13 @@ def graph_with_pen_and_color(expr q) =
fi
enddef ;
-% Add picture component q to picture @# and change part p to tp, where p is
-% something from q that needs coordinate transformation. The type of p is pair
-% or path.
+% Add picture component q to picture @# and change part p to tp,
+% where p is something from q that needs coordinate transformation.
+% The type of p is pair or path.
% Pair o is the value of p that makes tp (0,0). This implements the trick
-% whereby using 1 instead of 0 for th the width or height or the setbounds path
+% whereby using 1 instead of 0 for the width or height or the setbounds path
% for a label picture suppresses shifting in x or y.
+
%vardef graph_picture_conversion@#(expr q, o)(text tp) =
% save p ;
% if stroked q :
@@ -259,44 +295,47 @@ enddef ;
% else :
% interim truecorners :=0 ;
% pair p ; p=llcorner q;
-% if urcorner q<>p : p :=p+graph_coordinate_multiplication(o-p,urcorner q-p) ; fi
+% if urcorner q<>p : p := p + graph_coordinate_multiplication(o-p,urcorner q-p) ; fi
% addto @# also q shifted ((tp)-llcorner q) ;
% fi
%enddef ;
-% TH : new version from code found at sarovar tracker. This makes
-% grdaw clip the result to the window defined with setrange
+
+% This new version makes gdraw clip the result to the window defined with setrange
+
vardef graph_picture_conversion@#(expr q, o)(text tp) =
- save p, tp_geclipt ;
- picture tp_geclipt ; tp_geclipt :=nullpicture;
+ save p ;
+ save do_clip, tp_clipped ; boolean do_clip ; do_clip := true ;
+ picture tp_clipped ; tp_clipped := nullpicture;
if stroked q :
path p ; p=pathpart q;
- %%% --- SDV added
- addto tp_geclipt doublepath tp graph_with_pen_and_color(q) dashed dashpart q ;
- clip tp_geclipt to origin--(xpart Z_.graph_dimensions,0)--Z_.graph_dimensions--(0, ypart Z_.graph_dimensions)--cycle ;
- addto @# also tp_geclipt ;
- %%%
- %%% --- SDV deleted
- %%addto @# doublepath tp graph_with_pen_and_color(q) dashed dashpart q ;
- %%%
+ addto tp_clipped doublepath tp graph_with_pen_and_color(q) dashed dashpart q ;
+ %draw bbox tp_clipped withcolor red ;
elseif filled q :
path p ; p=pathpart q;
- addto @# contour tp graph_with_pen_and_color(q) ;
+ addto tp_clipped contour tp graph_with_pen_and_color(q) ;
+ %draw bbox tp_clipped withcolor green ;
else :
- interim truecorners :=0 ;
+ if (urcorner q<>llcorner q) : do_clip := false ; fi % Do not clip the axis labels;
+ interim truecorners := 0 ;
pair p ; p=llcorner q;
- if urcorner q<>p : p :=p+graph_coordinate_multiplication(o-p,urcorner q-p) ; fi
- addto @# also q shifted ((tp)-llcorner q) ;
+ if urcorner q<>p : p := p + graph_coordinate_multiplication(o-p,urcorner q-p) ; fi
+ addto tp_clipped also q shifted ((tp)-llcorner q) ;
+ %draw bbox tp_clipped withcolor if do_clip : cyan else : blue fi ;
+ fi
+ if do_clip :
+ clip tp_clipped to origin--(xpart Z_.graph_dimensions,0)--Z_.graph_dimensions--
+ (0,ypart Z_.graph_dimensions)--cycle ;
fi
+ addto @# also tp_clipped ;
enddef ;
def graph_coordinate_multiplication(expr a,b) = (xpart a*xpart b, ypart a*ypart b) enddef ;
-
vardef graph_clear_bounds@# = numeric @#.low, @#.high ; enddef;
-
% Finalize anything drawn in the present coordinate system and set up a new
% system as requested
+
vardef setcoords(expr tx, ty) =
interim warningcheck :=0 ;
if length graph_current_graph>0 :
@@ -311,10 +350,10 @@ vardef setcoords(expr tx, ty) =
X_.graph_coordinate_type := tx ; Y_.graph_coordinate_type := ty;
enddef ;
-
% Set Z_.low and Z_.high to correspond to given range of user graph
% coordinates. The text argument should be a sequence of pairs and/or strings
% with 4 components in all.
+
vardef setrange(text t) =
interim warningcheck :=0 ;
save r_ ; r_=0;
@@ -329,8 +368,8 @@ vardef setrange(text t) =
endfor
enddef ;
-
% @# is X_ or Y_ ; l and h are numeric or string
+
vardef graph_set_bounds@#(expr l, h) =
graph_clear_bounds@# ;
if @#graph_coordinate_type>0 :
@@ -358,15 +397,12 @@ vardef graph_set_bounds@#(expr l, h) =
fi
enddef ;
-
-
-
-
%%%%%%%%%%%%%%%%%%%%%%%%% Converting path coordinates %%%%%%%%%%%%%%%%%%%%%%%%%
% Find the result of scanning path p and using macros tx and ty to adjust the
% x and y parts of each coordinate pair. Boolean parameter c tells whether to
% force the result to be polygonal.
+
vardef graph_scan_path(expr p, c)(suffix tx, ty) =
if (str tx="") and (str ty="") : p
else :
@@ -385,10 +421,11 @@ vardef graph_scan_path(expr p, c)(suffix tx, ty) =
if pair p : point 0 of fi r_
fi
enddef ;
-vardef graph_pair_adjust(expr p)(suffix tx, ty) = (tx xpart p, ty ypart p) enddef ;
+vardef graph_pair_adjust(expr p)(suffix tx, ty) = (tx xpart p, ty ypart p) enddef ;
% Convert path p from user graph coords to internal graph coords.
+
vardef graph_convert_user_path_to_internal primary p =
interim warningcheck :=0 ;
graph_scan_path(p,
@@ -400,11 +437,11 @@ vardef graph_convert_user_path_to_internal primary p =
if Y_.graph_coordinate_type<0 : yscaled -1 fi)
enddef ;
-
% Convert label location t_ from user graph coords to internal graph coords.
% The label location should be a pair, or two numbers/strings. If t_ is empty
% or a single item of non-pair type, just return t_. Unknown coordinates
% produce unknown components in the result.
+
vardef graph_label_convert_user_to_internal(text t_) =
save n_ ; n_=0;
interim warningcheck :=0 ;
@@ -424,13 +461,12 @@ vardef graph_label_convert_user_to_internal(text t_) =
fi
enddef ;
-
-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Reading data files %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Read a line from file f, extract whitespace-separated tokens ignoring any
% initial "%", and return true if at least one token is found. The tokens
% are stored in @#1, @#2, .. with "" in the last @#[] entry.
+
vardef graph_read_line@#(expr f) =
save n_, s_ ; string s_;
s_ = readfrom f ;
@@ -448,9 +484,9 @@ vardef graph_read_line@#(expr f) =
fi
enddef ;
-
% Execute c for each line of data read from file f, and stop at the first
% line with no data. Commands c can use line number i and tokens $1, $2, ...
+
def gdata(expr f)(suffix $)(text c) =
boolean flag ;
for i=1 upto infinity :
@@ -462,8 +498,8 @@ def gdata(expr f)(suffix $)(text c) =
fi
enddef ;
-
% Read a path from file f. The path is terminated by blank line or EOF.
+
vardef graph_readpath(expr f) =
interim warningcheck :=0 ;
save s ;
@@ -473,13 +509,13 @@ vardef graph_readpath(expr f) =
)
enddef ;
-
% Append coordinates t to polygonal path @#. The coordinates can be numerics,
% strings, or a single pair.
+
vardef augment@#(text t) =
interim warningcheck := 0 ;
if not path begingroup @# endgroup :
- Gerr(begingroup @# endgroup, "Cannot augment--not a path") ;
+ graph_error(begingroup @# endgroup, "Cannot augment--not a path") ;
else :
def graph_comma= hide(def graph_comma=,enddef) enddef ;
if known @# : @# :=@#-- else : @#= fi
@@ -489,13 +525,13 @@ vardef augment@#(text t) =
fi
enddef ;
-
-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Drawing and filling %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Unknown pair components are set to 0 because glabel and gdotlabel understand
% unknown coordinates as `0 in absolute units'.
+
vardef graph_unknown_pair_bbox(expr p) =
+ interim warningcheck:=0 ;
if known p : addto graph_current_bb doublepath p ;
else :
save x,y ;
@@ -507,30 +543,38 @@ vardef graph_unknown_pair_bbox(expr p) =
graph_current_bb := image(fill llcorner graph_current_bb..urcorner graph_current_bb--cycle) ;
enddef ;
-
% Initiate a gdraw or gfill command. This must be done before scanning the
% argument, because that could invoke the `if known graph_plot_picture' test in a following
% plot option .
+
def graph_addto =
+ def graph_errorbar_text = enddef ;
+ color graph_foreground ;
+ path graph_last_path ;
graph_last_drawn := graph_plot_picture := nullpicture ; addto graph_last_drawn
enddef;
+% Handle the part of a gdraw command that uses path or data file p.
-% Handle the part of a Gdraw command that uses path or data file p.
def graph_draw expr p =
- if string p : graph_convert_user_path_to_internal graph_readpath(p)
- elseif path p or pair p : graph_convert_user_path_to_internal p
+ if string p : hide(graph_last_path := graph_readpath(p) ;)
+ graph_convert_user_path_to_internal graph_last_path
+ elseif path p or pair p :
+ hide(graph_last_path := p ;)
+ graph_convert_user_path_to_internal p
else : graph_error(p,"gdraw argument should be a data file or a path")
origin
fi
withpen currentpen graph_withlist _op_
enddef ;
+% Handle the part of a gdraw command that uses path or data file p.
-% Handle the part of a Gdraw command that uses path or data file p.
def graph_fill expr p =
- if string p : graph_convert_user_path_to_internal graph_readpath(p) --cycle
- elseif cycle p : graph_convert_user_path_to_internal p
+ if string p : hide(graph_last_path := graph_readpath(p) --cycle ;)
+ graph_convert_user_path_to_internal graph_last_path
+ elseif cycle p : hide(graph_last_path := p ;)
+ graph_convert_user_path_to_internal p
else : graph_error(p,"gfill argument should be a data file or a cyclic path")
origin..cycle
fi graph_withlist _op_
@@ -539,17 +583,88 @@ enddef ;
def gdraw = graph_addto doublepath graph_draw enddef ;
def gfill = graph_addto contour graph_fill enddef ;
-
% This is used in graph_draw and graph_fill to allow postprocessing graph_last_drawn
+
def graph_withlist text t_ = t_ ; graph_post_draw; enddef;
+def witherrorbars(text t) text options =
+ hide(
+ def graph_errorbar_text = t enddef ;
+ save pic ; picture pic ; pic := image(draw origin _op_ options ;) ;
+ if color colorpart pic : graph_foreground := colorpart pic ; fi
+ )
+ options
+enddef ;
+
+% new feature: graph_errorbars
+
+picture graph_errorbar_picture ; graph_errorbar_picture := image(draw (left--right) scaled .5 ;) ;
+%picture graph_xbar_picture ; graph_xbar_picture := image(draw (down--up) scaled .5 ;) ;
+%picture graph_ybar_picture ; graph_ybar_picture := image(draw (left--right) scaled .5 ;) ;
+
+vardef graph_errorbars(text t) =
+ if known graph_last_path :
+ save n, p, q ; path p ; pair q ;
+ save pic ; picture pic[] ; pic0 := nullpicture ;
+ pic1 := if known graph_xbar_picture : graph_xbar_picture
+ elseif known graph_errorbar_picture : graph_errorbar_picture rotated 90
+ else : nullpicture fi ;
+ pic2 := if known graph_ybar_picture : graph_ybar_picture
+ elseif known graph_errorbar_picture : graph_errorbar_picture
+ else : nullpicture fi ;
+ if length pic1>0 :
+ pic1 := pic1 scaled graph_shapesize ;
+ setbounds pic1 to origin..cycle ;
+ fi
+ if length pic2>0 :
+ pic2 := pic2 scaled graph_shapesize ;
+ setbounds pic2 to origin..cycle ;
+ fi
+ for i=0 upto length graph_last_path :
+ clearxy ; z = point i of graph_last_path ;
+ n := 1 ;
+ for $=t :
+ if known $ :
+ q := if path $ : if length $>i : point i of $ else : origin fi
+ elseif pair $ : $ elseif numeric $ : ($,$) else : origin fi ;
+ if q<>origin :
+ p := graph_convert_user_path_to_internal ((
+ if n=1 :
+ (-xpart q,0)--(ypart q,0)
+ else :
+ (0,-xpart q)--(0,ypart q)
+ fi ) shifted z) ;
+ addto pic0 doublepath p ;
+ if length pic[n]>0 :
+ if ypart q<>0 :
+ addto pic0 also pic[n] shifted point 1 of p ;
+ fi
+ if xpart q<>0 :
+ addto pic0 also pic[n] rotated 180 shifted point 0 of p ;
+ fi
+ fi
+ fi
+ fi
+ exitif incr n>3 ;
+ endfor
+ endfor
+ if length pic0>0 :
+ save bg, fg ; color bg, fg ;
+ bg := if known graph_background : graph_background else : background fi ;
+ fg := if known graph_foreground : graph_foreground else : black fi ;
+ addto graph_current_graph also pic0 withpen currentpen scaled 2 _op_ withcolor bg ;
+ addto graph_current_graph also pic0 withpen currentpen scaled .5 _op_ withcolor fg ;
+ fi
+ fi
+enddef ;
% Set graph_plot_picture so the postprocessing step will plot picture p at each path knot.
% Also select nullpen to suppress stroking.
+
def plot expr p =
if known graph_plot_picture :
withpen nullpen
- hide (graph_plot_picture :=image(
+ hide (graph_plot_picture := image(
if bounded p : for q within p : graph_addto_currentpicture q endfor % Save memory
else : graph_addto_currentpicture p
fi graph_setbounds origin..cycle))
@@ -557,24 +672,24 @@ def plot expr p =
enddef ;
% This hides a semicolon that could prematurely end graph_withlist's text argument
+
def graph_addto_currentpicture primary p = addto currentpicture also p ; enddef;
def graph_setbounds = setbounds currentpicture to enddef ;
-
-def gdrawarrow = graph_number_of_arrowheads :=1 ; gdraw enddef;
-def gdrawdblarrow = graph_number_of_arrowheads :=2 ; gdraw enddef;
-
+def gdrawarrow = graph_number_of_arrowheads := 1 ; gdraw enddef;
+def gdrawdblarrow = graph_number_of_arrowheads := 2 ; gdraw enddef;
% Post-process the filled or stroked picture graph_last_drawn as follows : (1) update
% the bounding box information ; (2) transfer it to graph_current_graph unless the pen has
% been set to nullpen to disable stroking ; (3) plot graph_plot_picture at each knot.
+
vardef graph_post_draw =
- save p ;
- path p ; p=pathpart graph_last_drawn;
+ save p ; path p ; p = pathpart graph_last_drawn ;
graph_unknown_pair_bbox(p) ;
if filled graph_last_drawn or not graph_is_null(penpart graph_last_drawn) :
addto graph_current_graph also graph_last_drawn ;
fi
+ graph_errorbars(graph_errorbar_text) ;
if length graph_plot_picture>0 :
for i=0 upto length p if cycle p : -1 fi :
addto graph_current_graph also graph_plot_picture shifted point i of p ;
@@ -586,17 +701,23 @@ vardef graph_post_draw =
if graph_number_of_arrowheads>1 :
graph_draw_arrowhead(reverse p, graph_with_pen_and_color(graph_last_drawn)) ;
fi
- graph_number_of_arrowheads :=0 ;
+ graph_number_of_arrowheads := 0 ;
fi
enddef ;
-vardef graph_is_null(expr p) = (urcorner p=origin) and (llcorner p=origin) enddef ;
+vardef graph_is_null(expr p) = (urcorner p=origin) and (llcorner p=origin) enddef ;
vardef graph_draw_arrowhead(expr p)(text w) = % Draw arrowhead for path p, with list w
+ %save r ; r := angle(precontrol infinity of p shifted -point infinity of p) ;
addto graph_current_graph also
- image(filldraw arrowhead(
- graph_arrowhead_extent(precontrol infinity of p, point infinity of p)) w ;
- graph_setbounds point infinity of p..cycle) ;
+ image(fill arrowhead (graph_arrowhead_extent(precontrol infinity of p,point infinity of p)) w ;
+ draw arrowhead (graph_arrowhead_extent(precontrol infinity of p,point infinity of p)) w
+ undashed ;
+%if (r mod 90 <> 0) : % orientation can be wrong due to remapping
+% draw textext("\tfxx " & decimal r) shifted point infinity of p withcolor blue ;
+%fi
+ graph_setbounds point infinity of p..cycle ;
+ ) ; % rotatedabout(point infinity of p,-r) ;
enddef ;
vardef graph_arrowhead_extent(expr p, q) =
@@ -604,8 +725,6 @@ vardef graph_arrowhead_extent(expr p, q) =
q
enddef ;
-
-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Drawing labels %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Argument c is a drawing command that needs an additional argument p that gives
@@ -613,6 +732,7 @@ enddef ;
% path. Unknown components of p cause the setbounds path to have width or height 1 instead of 0.
% Then graph_unknown_pair_bbox sets these components to 0 and graph_picture_conversion
% suppresses subsequent repositioning.
+
def graph_draw_label(expr p)(suffix $)(text c) =
save sdim_ ; pair sdim_;
sdim_ := (if unknown xpart p : 1+ fi 0, if unknown ypart p : 1+ fi 0) ;
@@ -621,14 +741,13 @@ def graph_draw_label(expr p)(suffix $)(text c) =
image(c(p) ; graph_setbounds p--p+sdim_--cycle) _op_
enddef ;
-
% Stash the result drawing command c in the graph_label table using with list w and
-% an index based on angle laboff$.
+% an index based on angle mfun_laboff$.
+
vardef graph_stash_label(suffix $)(text c) text w =
- graph_label[1.5+angle laboff$ /90] = image(c(origin) w) ;
+ graph_label[1.5+angle mfun_laboff$ /90] = image(c(origin) w) ;
enddef ;
-
def graph_label_location primary p =
if pair p : graph_draw_label(p)
elseif numeric p : graph_draw_label(point p of pathpart graph_last_drawn)
@@ -636,32 +755,33 @@ def graph_label_location primary p =
fi
enddef ;
-
% Place label p at user graph coords t using with list w. (t is a time, a pair
% or 2 numerics or strings).
+
vardef glabel@#(expr p)(text t) text w =
graph_label_location graph_label_convert_user_to_internal(t) (@#,label@#(p)) w ; enddef;
-
% Place label p at user graph coords t using with list w and draw a dot there.
% (t is a time, a pair, or 2 numerics or strings).
+
vardef gdotlabel@#(expr p)(text t) text w =
graph_label_location graph_label_convert_user_to_internal(t) (@#,dotlabel@#(p)) w ; enddef;
-
def OUT = enddef ; % location text for outside labels
-
-
%%%%%%%%%%%%%%%%%%%%%%%%%% Grid lines, ticks, etc. %%%%%%%%%%%%%%%%%%%%%%%%%%
% Grid lines and tick marks are transformed versions of the templates below.
% In the template paths, (0,0) is on the edge of the frame and inward is to
% the right.
-path Gtemplate.itick, Gtemplate.otick, Gtemplate.grid ;
-Gtemplate.itick = origin--(7bp,0) ;
-Gtemplate.otick = (-7bp,0)--origin ;
-Gtemplate.grid = origin--(1,0) ;
+
+path graph_template.tick, graph_template.itick, graph_template.otick, graph_template.grid ;
+graph_template.tick = (-3.5bp,0)--(3.5bp,0) ;
+graph_template.itick = origin--(7bp,0) ;
+graph_template.otick = (-7bp,0)--origin ;
+graph_template.grid = origin--(1,0) ;
+
+vardef tick@#(expr f,u) text w = graph_tick_label(@#,@,false,f,u,w) ; enddef;
vardef itick@#(expr f,u) text w = graph_tick_label(@#,@,false,f,u,w) ; enddef;
@@ -670,75 +790,82 @@ vardef otick@#(expr f,u) text w = graph_tick_label(@#,@,false,f,u,w) ; enddef;
vardef grid@#(expr f,u) text w = graph_tick_label(@#,@,true,f,u,w) ; enddef;
-% Produce a tick or grid mark for label suffix $, Gtemplate suffix $$,
-% coordinate value u, and with list w. Boolean c tells whether Gtemplate$$
+% Produce a tick or grid mark for label suffix $, graph_template suffix $$,
+% coordinate value u, and with list w. Boolean c tells whether graph_template$$
% needs scaling by X_.graph_dimensions or Y_.graph_dimensions,
% and f gives a format string or a label picture.
+
def graph_tick_label(suffix $,$$)(expr c, f, u)(text w) =
- graph_draw_label(graph_label_convert_user_to_internal(graph_generate_label_position($,u)),,draw graph_gridline_picture$($$,c,f,u,w) shifted)
+ graph_draw_label(graph_label_convert_user_to_internal(graph_generate_label_position($,u)),,
+ draw graph_gridline_picture$($$,c,f,u,w) shifted)
enddef ;
-
% Generate label positioning arguments appropriate for label suffix $ and
% coordinate u.
+
def graph_generate_label_position(suffix $)(expr u) =
- if xpart laboff.$=0 : u,whatever else : whatever,u fi
+ if pair u : u elseif xpart mfun_laboff.$=0 : u,whatever else : whatever,u fi
enddef ;
-
% Generate a picture of a grid line labeled with coordinate value u, picture
% or format string f, and with list w. Suffix @# is bot, top, lft, or rt,
-% suffix $ identifies entries in the Gtemplate table, and boolean c tells
-% whether to scale Gtemplate$.
+% suffix $ identifies entries in the graph_template table, and boolean c tells
+% whether to scale graph_template$.
+
vardef graph_gridline_picture@#(suffix $)(expr c, f, u)(text w) =
if unknown u : graph_error(u,"Label coordinate should be known") ; nullpicture
else :
save p ; path p;
interim warningcheck :=0 ;
graph_autogrid_needed :=false ;
- p = Gtemplate$ zscaled -laboff@#
- if c : Gxyscale fi
- shifted (((.5 + laboff@# dotprod (.5,.5)) * laboff@#) Gxyscale) ;
+ p = graph_template$ zscaled -mfun_laboff@#
+ if c : graph_xyscale fi
+ shifted (((.5 + mfun_laboff@# dotprod (.5,.5)) * mfun_laboff@#) graph_xyscale) ;
image(draw p w ;
label@#(if string f : format(f,u) else : f fi, point 0 of p))
fi
enddef ;
-def Gxyscale = xscaled X_.graph_dimensions yscaled Y_.graph_dimensions enddef ;
+def graph_xyscale = xscaled X_.graph_dimensions yscaled Y_.graph_dimensions enddef ;
% Draw the frame or the part corresponding to label suffix @# using with list w.
+
vardef frame@# text w =
graph_frame_needed :=false ;
picture p_ ;
p_ = image(draw
- if str@#<>"" : subpath round(angle laboff@#*graph_frame_pair_a+graph_frame_pair_b) of fi
- unitsquare Gxyscale w) ;
+ if str@#<>"" : subpath round(angle mfun_laboff@#*graph_frame_pair_a+graph_frame_pair_b) of fi
+ unitsquare graph_xyscale w) ;
graph_draw_label((whatever,whatever),,draw p_ shifted) ;
enddef ;
-pair graph_frame_pair_a ; graph_frame_pair_a=(1,1)/90; % unitsquare subpath is linear in label angle
-pair graph_frame_pair_b ; graph_frame_pair_b=(.75,2.25);
-
-
+pair graph_frame_pair_a ; graph_frame_pair_a=(1,1)/90; % unitsquare subpath is linear in label angle
+pair graph_frame_pair_b ; graph_frame_pair_b=(.75,2.25);
%%%%%%%%%%%%%%%%%%%%%%%%%% Automatic grid selection %%%%%%%%%%%%%%%%%%%%%%%%%%
-string Glmarks[] ; % marking options per decade for logarithmic scales
-string Gumarks ; % mark spacing options per decade for linear scales
-string Gemarks ; % exponent spacing options for logarithmic scales
-newinternal Gmarks, Gminlog ;
-Gmarks := 4 ; % minimum number marks generated by auto.x or auto.y
-Gminlog := mlog 3 ; % revert to uniform marks when largest/smallest < this
+string graph_log_marks[] ; % marking options per decade for logarithmic scales
+string graph_lin_marks ; % mark spacing options per decade for linear scales
+string graph_exp_marks ; % exponent spacing options for logarithmic scales
+newinternal graph_minimum_number_of_marks, graph_log_minimum ;
+graph_minimum_number_of_marks := 4 ; % minimum number marks generated by auto.x or auto.y
+graph_log_minimum := mlog 3 ; % revert to uniform marks when largest/smallest < this
+
+def Gfor(text t) = for i=t endfor enddef ; % to shorten the mark templates below
-def Gfor(text t) = for i=t endfor enddef ; % to shorten the mark templates below
-Glmarks[1]="1,2,5" ;
-Glmarks[2]="1,1.5,2,3,4,5,7" ;
-Glmarks[3]="1Gfor(6upto10 :,i/5)Gfor(5upto10 :,i/2)Gfor(6upto9 :,i)" ;
-Glmarks[4]="1Gfor(11upto20 :,i/10)Gfor(11upto25 :,i/5)Gfor(11upto19 :,i/2)" ;
-Glmarks[5]="1Gfor(21upto40 :,i/20)Gfor(21upto50 :,i/10)Gfor(26upto49 :,i/5)" ;
-Gumarks="10,5,2" ; % start with 10 and go down; a final `,1' is appended
-Gemarks="20,10,5,2,1" ;
+graph_log_marks[1]="1,2,5" ;
+graph_log_marks[2]="1,1.5,2,3,4,5,7" ;
+graph_log_marks[3]="1Gfor(6upto10 :,i/5)Gfor(5upto10 :,i/2)Gfor(6upto9 :,i)" ;
+graph_log_marks[4]="1Gfor(11upto20 :,i/10)Gfor(11upto25 :,i/5)Gfor(11upto19 :,i/2)" ;
+graph_log_marks[5]="1Gfor(21upto40 :,i/20)Gfor(21upto50 :,i/10)Gfor(26upto49 :,i/5)" ;
+graph_lin_marks="10,5,2" ; % start with 10 and go down; a final `,1' is appended
+graph_exp_marks="20,10,5,2,1" ;
+Ten_to0 = 1 ;
+Ten_to1 = 10 ;
+Ten_to2 = 100 ;
+Ten_to3 = 1000 ;
+Ten_to4 = 10000 ;
% Determine the X_ or Y_ bounds on the range to be covered by automatic grid
% marks. Suffix @# is X_ or Y_. The result is log or linear to specify the
@@ -747,6 +874,7 @@ Gemarks="20,10,5,2,1" ;
% are upper and lower bounds in
% `modified exponential form'. In modified exponential form, (x,y) means
% (x/1000)*10^y, where 1000<=abs x<10000.
+
vardef graph_bounds@# =
interim warningcheck :=0 ;
save l, h ;
@@ -755,28 +883,29 @@ vardef graph_bounds@# =
if abs @#graph_coordinate_type=log :
graph_modified_lower := graph_Meform(l)+graph_modified_bias ;
graph_modified_higher := graph_Meform(h)+graph_modified_bias ;
- if h-l >= Gminlog : log else : linear fi
+ if h-l >= graph_log_minimum : log else : linear fi
else :
graph_modified_lower := graph_Feform(l)+graph_modified_bias ;
graph_modified_higher := graph_Feform(h)+graph_modified_bias ;
linear
fi
enddef ;
+
pair graph_modified_bias ; graph_modified_bias=(0,3);
pair graph_modified_lower, graph_modified_higher ;
+% Scan graph_log_marks[k] and evaluate tokens t for each m where l<=m<=h.
-% Scan Glmarks[k] and evaluate tokens t for each m where l<=m<=h.
def graph_scan_marks(expr k, l, h)(text t) =
- for m=scantokens Glmarks[k] :
+ for m=scantokens graph_log_marks[k] :
exitif m>h ;
if m>=l : t fi
endfor
enddef ;
-
-% Scan Gmark[k] and evaluate tokens t for each m and e where m*10^e belongs
+% Scan graph_log_marks[k] and evaluate tokens t for each m and e where m*10^e belongs
% between l and h (inclusive), where both l and h are in modified exponent form.
+
def graph_scan_mark(expr k, l, h)(text t) =
for e=ypart l upto ypart h :
graph_scan_marks(k, if e>ypart l : 1 else : xpart l/1000 fi,
@@ -784,27 +913,29 @@ def graph_scan_mark(expr k, l, h)(text t) =
endfor
enddef ;
-
% Select a k for which graph_scan_mark(k,...) gives enough marks.
+
vardef graph_select_mark =
save k ;
k = 0 ;
forever :
- exitif unknown Glmarks[k+1] ;
- exitif 0 graph_scan_mark(incr k, graph_modified_lower, graph_modified_higher, +1) >= Gmarks ;
+ exitif unknown graph_log_marks[k+1] ;
+ exitif 0 graph_scan_mark(incr k, graph_modified_lower, graph_modified_higher, +1)
+ >= graph_minimum_number_of_marks ;
endfor
k
enddef ;
-
-% Try to select an exponent spacing from Gemarks. If successful, set @# and
+% Try to select an exponent spacing from graph_exp_marks. If successful, set @# and
% return true
+
vardef graph_select_exponent_mark@# =
numeric @# ;
- for e=scantokens Gemarks :
+ for e=scantokens graph_exp_marks :
@# = e ;
exitif floor(ypart graph_modified_higher/e) -
- floor(graph_modified_exponent_ypart(graph_modified_lower)/e) >= Gmarks ;
+ floor(graph_modified_exponent_ypart(graph_modified_lower)/e)
+ >= graph_minimum_number_of_marks ;
numeric @# ;
endfor
known @#
@@ -812,17 +943,17 @@ enddef ;
vardef graph_modified_exponent_ypart(expr p) = ypart p if xpart p=1000 : -1 fi enddef ;
-
% Compute the mark spacing d between xpart graph_modified_lower and xpart graph_modified_higher.
+
vardef graph_tick_mark_spacing =
interim warningcheck :=0 ;
save m, n, d ;
- m = Gmarks ;
+ m = graph_minimum_number_of_marks ;
n = 1 for i=1 upto
(mlog(xpart graph_modified_higher-xpart graph_modified_lower) - mlog m)/mlogten :
*10 endfor ;
if n<=1000 :
- for x=scantokens Gumarks :
+ for x=scantokens graph_lin_marks :
d = n*x ;
exitif 0 graph_generate_numbers(d,+1)>=m ;
numeric d ;
@@ -831,25 +962,24 @@ vardef graph_tick_mark_spacing =
if known d : d else : n fi
enddef ;
-
def graph_generate_numbers(expr d)(text t) =
for m = d*ceiling(xpart graph_modified_lower/d) step d until xpart graph_modified_higher :
t
endfor
enddef ;
-
% Evaluate tokens t for exponents e in multiples of d in the range determined
% by graph_modified_lower and graph_modified_higher.
+
def graph_generate_exponents(expr d)(text t) =
for e = d*floor(graph_modified_exponent_ypart(graph_modified_lower)/d+1)
step d until d*floor(ypart graph_modified_higher/d) : t
endfor
enddef ;
-
% Adjust graph_modified_lower and graph_modified_higher so their exponent parts match
% and they are in true exponent form ((x,y) means x*10^y). Return the new exponent.
+
vardef graph_match_exponents =
interim warningcheck := 0 ;
save e ;
@@ -862,10 +992,10 @@ vardef graph_match_exponents =
e
enddef ;
-
% Assume e is an integer and either m=0 or 1<=abs(m)<10000. Find m*(10^e)
% and represent the result as a string if its absolute value would be at least
% 4096 or less than .1. It is OK to return 0 as a string or a numeric.
+
vardef graph_factor_and_exponent_to_string(expr m, e) =
if (e>3)or(e<-4) :
decimal m & "e" & decimal e
@@ -880,7 +1010,6 @@ vardef graph_factor_and_exponent_to_string(expr m, e) =
fi
enddef ;
-
def auto suffix $ =
hide(def graph_comma= hide(def graph_comma=,enddef) enddef)
if graph_bounds.graph_suffix($)=log :
@@ -898,7 +1027,6 @@ def auto suffix $ =
fi
enddef ;
-
string Autoform ; Autoform = "%g";
%vardef autogrid(suffix tx, ty) text w =
@@ -949,12 +1077,11 @@ vardef autogrid(suffix tx, ty) text w =
fi
enddef ;
-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% endgraph %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def endgraph =
if graph_autogrid_needed : autogrid(otick.bot, otick.lft) ; fi
- if graph_frame_needed : frame ; fi
+ if graph_frame_needed : frame ; fi
setcoords(linear,linear) ;
interim truecorners :=1 ;
for b=bbox graph_finished_graph :
@@ -971,27 +1098,35 @@ enddef ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-Ten_to0 = 1 ;
-Ten_to1 = 10 ;
-Ten_to2 = 100 ;
-Ten_to3 = 1000 ;
-Ten_to4 = 10000 ;
+% We format in luatex (using \mathematics{}) ...
+% we could pass via variables and save escaping as that is inefficient
-vardef escaped_format(expr s) =
- "" for n=1 upto length(s) : &
- if ASCII substring (n,n+1) of s = 37 :
- "@"
- else :
- substring (n,n+1) of s
- fi
- endfor
-enddef ;
+if unknown context_mlib :
-vardef strfmt(expr f, x) =
- "\MPgraphformat{" & escaped_format(f) & "}{" & (if string x : x else : decimal x fi) & "}"
-enddef ;
+ vardef escaped_format(expr s) =
+ "" for n=0 upto length(s) : &
+ if ASCII substring (n,n+1) of s = 37 :
+ "@"
+ else :
+ substring (n,n+1) of s
+ fi
+ endfor
+ enddef ;
+
+ vardef strfmt(expr f, x) = % maybe use mfun_ namespace
+ "\MPgraphformat{" & escaped_format(f) & "}{" & mfun_tagged_string(x) & "}"
+ enddef ;
+
+ vardef varfmt(expr f, x) = % maybe use mfun_ namespace
+ "\MPformatted{" & escaped_format(f) & "}{" & mfun_tagged_string(x) & "}"
+ enddef ;
+
+ vardef format (expr f, x) = textext(strfmt(f,x)) enddef ;
+ vardef formatted(expr f, x) = textext(varfmt(f,x)) enddef ;
-vardef format(expr f, x) = textext(strfmt(f, x)) enddef ;
+fi ;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A couple of extensions :
@@ -999,7 +1134,7 @@ vardef format(expr f, x) = textext(strfmt(f, x)) enddef ;
% unfilled outline, interior filled with different shades of the background.
% This allows overlapping points on a plot to be more distinguishable.
-vardef graph_shapesize = .33BodyFontSize enddef ;
+vardef graph_shapesize = (.33BodyFontSize) enddef ;
path graph_shape[] ; % (internal) symbol path
@@ -1073,25 +1208,35 @@ def downtriangles(expr f) = plotsymbol(13,f) enddef ;
def lefttriangles(expr f) = plotsymbol(33,f) enddef ;
def righttriangles(expr f) = plotsymbol(23,f) enddef ;
+% f (fill) is color, numeric or boolean, otherwise background.
def plotsymbol(expr n, f) text t =
if known graph_shape[n] :
image(
save bg, fg ; color bg, fg ;
bg := if known graph_background : graph_background else : background fi ;
save pic ; picture pic ; pic := image(draw origin _op_ t ;) ;
- fg := if color colorpart pic : colorpart pic else : black fi ;
+ if color colorpart pic : graph_foreground := colorpart pic ; fi
+ fg := if known graph_foreground : graph_foreground else : black fi ;
save p ; path p ; p = graph_shape[n] scaled graph_shapesize ;
draw p withcolor bg withpen currentpen scaled 2 ; % halo
- if cycle p : fill p withcolor
- if color f and known f :
- f
- elseif numeric f and known f :
- f[bg,fg]
- else :
- bg
- fi ;
+ currentpen := currentpen scaled .5 ;
+ if cycle p :
+ fill p withcolor
+ if known f :
+ if color f :
+ f
+ elseif numeric f :
+ f[bg,fg]
+ elseif boolean f and f :
+ fg
+ else
+ bg
+ fi
+ else :
+ bg
+ fi ;
fi
- draw p withpen currentpen _op_ t ;
+ draw p _op_ t ;
)
else :
nullpicture
@@ -1164,7 +1309,7 @@ enddef ;
def makefunctionpath (expr f, t, n) (text func) =
(for x=f step ((t-f)/(abs n)) until t :
- if x<>f : .. fi
+ if x<>f : -- fi
(x, func)
endfor )
enddef ;
@@ -1269,9 +1414,9 @@ numeric fit_chi_squared ;
vardef polynomial_fit (suffix p, $) (expr n) (text t) =
if not path p :
- Gerr(p, "Cannot fit--not a path") ;
+ graph_error(p, "Cannot fit--not a path") ;
elseif length p < n :
- Gerr(p, "Cannot fit--not enough points") ;
+ graph_error(p, "Cannot fit--not enough points") ;
else :
fit_chi_squared := 0 ;
% calculate sums of the data
@@ -1285,7 +1430,24 @@ vardef polynomial_fit (suffix p, $) (expr n) (text t) =
endfor
for i=0 upto length p :
clearxy ; z = point i of p ;
- w := if length(t) > 0 : t else : 1 fi ; % weight
+ w := 1 ; % weight
+ if known t :
+ if numeric t :
+ w := 1 if t<>0 : /(abs t) fi ;
+ elseif pair t :
+ if t<>origin :
+ w := 1/(abs t) ;
+ fi
+ elseif path t :
+ if length t>= i:
+ if point i of t<>origin :
+ w := 1/(abs point i of t) ;
+ fi
+ else :
+ w := 0 ;
+ fi ;
+ fi
+ fi
x1 := w ;
for j=0 upto 2n :
sumx[j] := sumx[j] + x1 ;
@@ -1356,13 +1518,29 @@ vardef exponential_function (suffix $) (expr x) = $1*exp($0*x) enddef ;
vardef exponential_fit (suffix p, $) (text t) =
save a ; numeric a[] ;
- save q ; path q ; % fit to the log of the ordinate
+ save q ; path q[] ; % fit to the log of the ordinate
for i=0 upto length p :
- if ypart(point i of p)>0 :
- augment.q(xpart(point i of p),ln(ypart(point i of p))) ;
+ clearxy ; z = point i of p ;
+ if y>0 :
+ augment.q0(x,ln(y)) ;
+ augment.q1(
+ if known t :
+ if numeric t : (0,ln(t))
+ elseif pair t : (xpart t,ln(ypart t))
+ elseif path t :
+ if length t>=i :
+ hide(z1 = point i of t;)
+ (x1,ln(y1))
+ else :
+ origin
+ fi
+ fi
+ else :
+ (0,1)
+ fi ) ;
fi
endfor
- linear_fit(q,a,t) ;
+ linear_fit(q0,a,q1) ;
save e ; e := exp(sqrt(fit_chi_squared)) ;
fit_chi_squared := e * e ;
$0 := a1 ;
@@ -1377,13 +1555,29 @@ vardef power_law_function (suffix $) (expr x) = $1*(x**$0) enddef ;
vardef power_law_fit (suffix p, $) (text t) =
save a ; numeric a[] ;
- save q ; path q ; % fit to the logs of the abscissae and ordinates
+ save q ; path q[] ; % fit to the logs of the abscissae and ordinates
for i=0 upto length p :
- if (xpart(point i of p)>0) and (ypart(point i of p)>0) :
- augment.q(ln(xpart(point i of p)),ln(ypart(point i of p))) ;
+ clearxy ; z = point i of p ;
+ if (x>0) and (y>0) :
+ augment.q0(ln(x),ln(y)) ;
+ augment.q1(
+ if known t :
+ if numeric t : (0,ln(t))
+ elseif pair t : (ln(xpart t),ln(ypart t))
+ elseif path t :
+ if length t>=i :
+ hide(z1 = point i of t)
+ (ln(x1),ln(y1))
+ else :
+ origin
+ fi
+ fi
+ else :
+ (0,1)
+ fi ) ;
fi
endfor
- linear_fit(q,a,t) ;
+ linear_fit(q0,a,q1) ;
save e ; e := exp(sqrt(fit_chi_squared)) ;
fit_chi_squared := e * e ;
$0 := a1 ;
@@ -1411,13 +1605,29 @@ enddef ;
vardef gaussian_fit (suffix p, $) (text t) =
save a ; numeric a[] ;
- save q ; path q ; % fit to the log of the ordinate
+ save q ; path q[] ; % fit to the log of the ordinate
for i=0 upto length p :
- if ypart(point i of p)>0 :
- augment.q(xpart(point i of p), ln(ypart(point i of p))) ;
+ clearxy ; z = point i of p ;
+ if y>0 :
+ augment.q0(x,ln(y)) ;
+ augment.q1(
+ if known t :
+ if numeric t : (0,ln(t))
+ elseif pair t : (xpart t,ln(ypart t))
+ elseif path t :
+ if length t>=i :
+ hide(z1 = point i of t)
+ (x1,ln(y1))
+ else :
+ origin
+ fi
+ fi
+ else :
+ (0,1)
+ fi ) ;
fi
endfor
- polynomial_fit(q,a,2,if t > 0 : ln(t) else : 0 fi) ;
+ polynomial_fit(q0,a,2,q1) ;
save e ; e := exp(sqrt(fit_chi_squared)) ;
fit_chi_squared := e * e ;
$1 := sqrt(-lntwo/a2) ;