1
2
3
4
5
6
7
8
9
10
11
12
13
14if known context_grap : endinput ; fi ;
15
16boolean context_grap ; context_grap := true ;
17
18
19
20message ("using number system " & numbersystem & " with precision " & decimal numberprecision) ;
21
22if numbersystem <> "double" :
23 errmessage "The graph macros require the double precision number system." ;
24 endinput ;
25fi
26
27if known contextlmtxmode :
28 def _op_ = base_draw_options enddef ;
29fi ;
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79newinternal Mzero ; Mzero := -16384;
80newinternal mlogten ; mlogten := mlog(10) ;
81newinternal largestmantissa ; largestmantissa := 252 ;
82newinternal singleinfinity ; singleinfinity := 2128 ;
83newinternal doubleinfinity ; doubleinfinity := 21024 ;
84
85
86
87
88vardef graph_mlog primary x =
89 if unknown x: whatever
90 elseif x=0: Mzero
91 else: mlog(abs x) fi
92enddef ;
93
94vardef graph_exp primary x =
95 if unknown x: whatever
96 elseif x<=Mzero: 0
97 else: mexp(x) fi
98enddef ;
99
100
101
102
103vardef logten primary x =
104 if unknown x: whatever
105 elseif x=0: Mzero
106 else: mlog(abs x)mlog(10) fi
107enddef ;
108
109if known contextlmtxmode :
110
111else :
112
113 vardef ln primary x =
114 if unknown x: whatever
115 elseif x=0: Mzero
116 else: mlog(abs x)256 fi
117 enddef ;
118
119 vardef exp primary x =
120 if unknown x: whatever
121 elseif x<= Mzero: 0
122 else: (mexp 256)x fi
123 enddef ;
124
125fi
126
127vardef powten primary x =
128 if unknown x: whatever
129 elseif x<= Mzero: 0
130 else: 10x fi
131enddef ;
132
133
134
135
136vardef graph_Meform(expr x) =
137 if x<=Mzero : origin
138 else :
139 save e, m ; e=floor(xmlogten)-3; m := mexp(xemlogten) ;
140 if abs m<1000 : m := m10 ; e := e-1 ; elseif abs m>=10000 : m := m10 ; e := e+1 ; fi
141 (m, e)
142 fi
143enddef ;
144
145
146
147vardef graph_Feform(expr x) =
148 interim warningcheck :=0 ;
149 if x=0 : origin
150 else :
151 save e, m ; e=floor(if x<0 : mlog(x) else : mlog(x) fimlogten)-3; m := x(10e) ;
152 if abs m<1000 : m := m10 ; e := e-1 ; elseif abs m>=10000 : m := m10 ; e := e+1 ; fi
153 (m, e)
154 fi
155enddef ;
156
157vardef graph_error(expr x,s) =
158 interim showstopping :=0 ;
159 show x ; errmessage s ;
160enddef ;
161
162
163
164vardef Z_@# = (X_@#,Y_@#) enddef ;
165
166def graph_suffix(suffix $) =
167 if str$="x" : X_ else : Y_ fi
168enddef ;
169
170
171
172save graph_background ; color graph_background ;
173
174def begingraph(expr w, h) =
175 begingroup
176 save X_, Y_ ;
177 X_.graph_coordinate_type =
178 Y_.graph_coordinate_type = linear ;
179 Z_.graph_dimensions = (w,h) ;
180
181
182 save graph_finished_graph ;
183 picture graph_finished_graph ;
184 graph_finished_graph = nullpicture ;
185 save graph_current_graph ;
186 picture graph_current_graph ;
187 graph_current_graph = nullpicture ;
188 save graph_current_bb ;
189 picture graph_current_bb ;
190 graph_current_bb = nullpicture ;
191 save graph_last_drawn ;
192 picture graph_last_drawn ;
193 graph_last_drawn = nullpicture ;
194 save graph_last_path ;
195 path graph_last_path ;
196 save graph_plot_picture ;
197 picture graph_plot_picture ;
198 save graph_foreground ;
199 color graph_foreground ;
200 save graph_label ;
201 picture graph_label[] ;
202 save graph_autogrid_needed ;
203 boolean graph_autogrid_needed ;
204 graph_autogrid_needed = true ;
205 save graph_frame_needed ;
206 boolean graph_frame_needed ;
207 graph_frame_needed = true ;
208 save graph_number_of_arrowheads ;
209 graph_number_of_arrowheads = 0 ;
210
211 if known graph_background :
212 addto graph_finished_graph contour
213 origin--(w,0)--(w,h)--(0,h)cycle withcolor graph_background ;
214 fi
215enddef ;
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232newinternal logarithmic, linear ;
233logarithmic :=1 ; linear :=2;
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249vardef graph_set_default_bounds =
250 forsuffixes $=low,high :
251 (if known X_$ : whatever else : X_$ fi, if known Y_$ : whatever else : Y_$ fi)
252 = graph_margin_fraction$[llcorner graph_current_bb,urcorner graph_current_bb]
253 graph_margin_pair$ ;
254 endfor
255enddef ;
256
257pair graph_margin_pair.low, graph_margin_pair.high ;
258graph_margin_pair.high = graph_margin_pair.low = (.00002,.00002) ;
259
260
261
262
263
264
265vardef graph_remap(suffix $,$$,$$$) =
266 save p_ ;
267 graph_set_default_bounds ;
268 pair p_, $ ; $=Z_.low;
269 p_ = (max(X_.highX_.low,.9), max(Y_.highY_.low,.9)) ;
270 transform $$, $$$ ;
271 forsuffixes #=$$,$$$ : xpart#=ypart#=xypart#=yxpart#=0 ; endfor
272 (Z_.high$) transformed $$ = p_ ;
273 p_ transformed $$$ = Z_.graph_dimensions ;
274enddef ;
275
276graph_margin_fraction.low=-.07 ;
277graph_margin_fraction.high=1.07 ;
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295def graph_with_pen_and_color(expr q) =
296 withproperties q
297enddef ;
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325vardef graph_picture_conversion@#(expr q, o)(text tp) =
326 save p ;
327 save do_clip, tp_clipped ; boolean do_clip ; do_clip := true ;
328 picture tp_clipped ; tp_clipped := nullpicture;
329 if stroked q :
330 path p ; p=pathpart q;
331 addto tp_clipped doublepath tp graph_with_pen_and_color(q) dashed dashpart q ;
332
333 elseif filled q :
334 path p ; p=pathpart q;
335 addto tp_clipped contour tp graph_with_pen_and_color(q) ;
336
337 else :
338 if (urcorner q<>llcorner q) : do_clip := false ; fi
339 interim truecorners := 0 ;
340 pair p ; p=llcorner q;
341 if urcorner q<>p : p := p graph_coordinate_multiplication(op,urcorner qp) ; fi
342 addto tp_clipped also q shifted ((tp)llcorner q) ;
343
344 fi
345 if do_clip :
346 clip tp_clipped to origin--(xpart Z_.graph_dimensions,0)Z_.graph_dimensions--
347 (0,ypart Z_.graph_dimensions)cycle ;
348 fi
349 addto @# also tp_clipped ;
350enddef ;
351
352def graph_coordinate_multiplication(expr a,b) = (xpart axpart b, ypart aypart b) enddef ;
353
354vardef graph_clear_bounds@# = numeric @#.low, @#.high ; enddef;
355
356
357
358
359vardef setcoords(expr tx, ty) =
360 interim warningcheck :=0 ;
361 if length graph_current_graph>0 :
362 save s, S, T ;
363 graph_remap(s, S, T) ;
364 for q within graph_current_graph :
365 graph_picture_conversion.graph_finished_graph(q,s,p shifted s transformed S transformed T) ;
366 endfor
367 graph_current_graph := graph_current_bb := nullpicture ;
368 fi
369 graph_clear_bounds.X_ ; graph_clear_bounds.Y_;
370 X_.graph_coordinate_type := tx ; Y_.graph_coordinate_type := ty;
371enddef ;
372
373
374
375
376
377vardef setrange(text t) =
378 interim warningcheck :=0 ;
379 save r_ ; r_=0;
380 string r_[]s ;
381 for x_=
382 for p_=t : if pair p_ : xpart p_, ypart fi p_, endfor :
383 r_[incr r_] if string x_ : s fi = x_ ;
384 if r_>2 :
385 graph_set_bounds if r_=3 : X_ else : Y_ fi (r_[r_-2] if unknown r_[r_-2] : s fi, x_) ;
386 fi
387 exitif r_=4 ;
388 endfor
389enddef ;
390
391
392
393vardef graph_set_bounds@#(expr l, h) =
394 graph_clear_bounds@# ;
395 if @#graph_coordinate_type>0 :
396 @#low = if unknown l :
397 whatever
398 else :
399 if abs @#graph_coordinate_type=logarithmic : graph_mlog fi if string l : scantokens fi l
400 fi ;
401 @#high = if unknown h :
402 whatever
403 else :
404 if abs @#graph_coordinate_type=logarithmic : graph_mlog fi if string h : scantokens fi h
405 fi ;
406 else :
407 @#high = if unknown l :
408 whatever
409 else :
410 if abs @#graph_coordinate_type=logarithmic : graph_mlog fi if string l : scantokens fi l
411 fi ;
412 @#low = if unknown h :
413 whatever
414 else :
415 if abs @#graph_coordinate_type=logarithmic : graph_mlog fi if string h : scantokens fi h
416 fi ;
417 fi
418enddef ;
419
420
421
422
423
424
425
426vardef graph_scan_path(expr p, c)(suffix tx, ty) =
427 if (str tx="") and (str ty="") : p
428 else :
429 save r_ ; path r_;
430 r_ := graph_pair_adjust(point 0 of p, tx, ty)
431 if path p :
432 for t=1 upto length p :
433 if c : --
434 else : ..controls graph_pair_adjust(postcontrol(t-1) of p, tx, ty)
435 and graph_pair_adjust(precontrol t of p, tx, ty) ..
436 fi
437 graph_pair_adjust(point t of p, tx, ty)
438 endfor
439 if cycle p : cycle fi
440 fi ;
441 if pair p : point 0 of fi r_
442 fi
443enddef ;
444
445vardef graph_pair_adjust(expr p)(suffix tx, ty) = (tx xpart p, ty ypart p) enddef ;
446
447
448
449vardef graph_convert_user_path_to_internal primary p =
450 interim warningcheck :=0 ;
451 if known p :
452 graph_scan_path(p,
453 (abs X_.graph_coordinate_type<>linear) or (abs Y_.graph_coordinate_type<>linear),
454 if abs X_.graph_coordinate_type=logarithmic : graph_mlog fi,
455 if abs Y_.graph_coordinate_type=logarithmic : graph_mlog fi)
456 transformed (identity
457 if X_.graph_coordinate_type<0 : xscaled -1 fi
458 if Y_.graph_coordinate_type<0 : yscaled -1 fi)
459 fi
460enddef ;
461
462
463
464
465
466
467vardef graph_label_convert_user_to_internal(text t_) =
468 save n_ ; n_=0;
469 interim warningcheck :=0 ;
470 if 0 for x_=t_ : +1 if pair x_ : +1 fi endfor <= 1 :
471 t_
472 else :
473 n_0 = n_1 = 0 ;
474 point 0 of graph_convert_user_path_to_internal (
475 for x_=
476 for y_=t_ : if pair y_ : xpart y_, ypart fi y_, endfor
477 0, 0 :
478 if known x_ : if string x_ : scantokens fi x_
479 else : hide(n_[n_] :=whatever) 0
480 fi
481 exitif incr n_=2 ;
482 ,endfor) (n_0,n_1)
483 fi
484enddef ;
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501vardef graph_cspan(expr s)(text c) =
502 0
503 for i=1 upto length s:
504 exitunless c substring (i-1,i) of s;
505 1
506 endfor
507enddef ;
508
509
510
511
512vardef graph_loptok suffix s =
513 save t, k;
514 k = graph_cspan(s," ">=);
515 if k > 0 :
516 s := substring(k,infinity) of s ;
517 fi
518 k := graph_cspan(s," "<);
519 string t;
520 t = substring (0,k) of s;
521 s := substring (k,infinity) of s;
522 s := substring (graph_cspan(s," ">=),infinity) of s;
523 t
524enddef ;
525
526vardef graph_read_line@#(expr f) =
527 save n_, s_ ; string s_;
528 s_ = readfrom f ;
529 string @#[] ;
530 if s_<>EOF :
531 @#0 := s_ ;
532 @#1 := graph_loptok s_ ;
533 n_ = if @#1="%" : 0 else : 1 fi ;
534 forever :
535 @#[incr n_] := graph_loptok s_ ;
536 exitif @#[n_]="" ;
537 endfor
538 @#1<>""
539 else : false
540 fi
541enddef ;
542
543
544
545
546
547
548
549
550
551
552
553
554def gdata(expr f)(suffix $)(text c) =
555 for i=1 upto largestmantissa :
556 exitunless graph_read_line$(f) ;
557 c
558 endfor
559enddef ;
560
561
562
563vardef graph_readpath(expr f) =
564 interim warningcheck :=0 ;
565 save s ;
566 gdata(f, s, if i>1 :fi
567 if s2="" : ( i, scantokens s1)
568 else : (scantokens s1, scantokens s2) fi
569 )
570enddef ;
571
572
573
574
575vardef augment@#(text t) =
576 interim warningcheck := 0 ;
577 if not path begingroup @# endgroup :
578 graph_error(begingroup @# endgroup, "Cannot augment--not a path") ;
579 else :
580 def graph_comma= hide(def graph_comma=,enddef) enddef ;
581 if known @# : @# :=@#-- else : @#= fi
582 (for p=t :
583 graph_comma if string p : scantokens fi p
584 endfor) ;
585 fi
586enddef ;
587
588
589
590
591
592
593vardef graph_unknown_pair_bbox(expr p) =
594 interim warningcheck:=0 ;
595 if known p : addto graph_current_bb doublepath p ;
596 else :
597 save x,y ;
598 z = llcorner graph_current_bb ;
599 if unknown xpart p : xpart p= else : x := fi 0 ;
600 if unknown ypart p : ypart p= else : y := fi 0 ;
601 addto graph_current_bb doublepath (pz) ;
602 fi
603 graph_current_bb := image(fill llcorner graph_current_bb..urcorner graph_current_bbcycle) ;
604enddef ;
605
606
607
608
609
610def graph_addto =
611 def graph_errorbar_text = enddef ;
612 color graph_foreground ;
613 path graph_last_path ;
614 graph_last_drawn := graph_plot_picture := nullpicture ; addto graph_last_drawn
615enddef;
616
617
618
619def graph_draw expr p =
620 if string p : hide(graph_last_path := graph_readpath(p) ;)
621 graph_convert_user_path_to_internal graph_last_path
622 elseif path p or pair p :
623 hide(graph_last_path := p ;)
624 graph_convert_user_path_to_internal p
625 else : graph_error(p,"gdraw argument should be a data file or a path")
626 origin
627 fi
628 withpen currentpen graph_withlist _op_
629enddef ;
630
631
632
633def graph_fill expr p =
634 if string p : hide(graph_last_path := graph_readpath(p) cycle ;)
635 graph_convert_user_path_to_internal graph_last_path
636 elseif cycle p : hide(graph_last_path := p ;)
637 graph_convert_user_path_to_internal p
638 else : graph_error(p,"gfill argument should be a data file or a cyclic path")
639 origin..cycle
640 fi graph_withlist _op_
641enddef ;
642
643def gdraw = graph_addto doublepath graph_draw enddef ;
644def gfill = graph_addto contour graph_fill enddef ;
645
646
647
648def graph_withlist text t_ = t_ ; graph_post_draw; enddef;
649
650def witherrorbars(text t) text options =
651 hide(
652 def graph_errorbar_text = t enddef ;
653 save pic ; picture pic ; pic := image(draw origin _op_ options ;) ;
654 if color colorpart pic : graph_foreground := colorpart pic ; fi
655 )
656 options
657enddef ;
658
659
660
661picture graph_errorbar_picture ; graph_errorbar_picture := image(draw (leftright) scaled .5 ;) ;
662
663
664
665vardef graph_errorbars(text t) =
666 if known graph_last_path :
667 save n, p, q ; path p ; pair q ;
668 save pic ; picture pic[] ; pic0 := nullpicture ;
669 pic1 := if known graph_xbar_picture : graph_xbar_picture
670 elseif known graph_errorbar_picture : graph_errorbar_picture rotated 90
671 else : nullpicture fi ;
672 pic2 := if known graph_ybar_picture : graph_ybar_picture
673 elseif known graph_errorbar_picture : graph_errorbar_picture
674 else : nullpicture fi ;
675 if length pic1>0 :
676 pic1 := pic1 scaled graph_shapesize ;
677 setbounds pic1 to origin..cycle ;
678 fi
679 if length pic2>0 :
680 pic2 := pic2 scaled graph_shapesize ;
681 setbounds pic2 to origin..cycle ;
682 fi
683 for i=0 upto length graph_last_path :
684 clearxy ; z = point i of graph_last_path ;
685 n := 1 ;
686 for $=t :
687 if known $ :
688 q := if path $ : if length $>i : point i of $ else : origin fi
689 elseif pair $ : $ elseif numeric $ : ($,$) else : origin fi ;
690 if q<>origin :
691 p := graph_convert_user_path_to_internal ((
692 if n=1 :
693 (xpart q,0)--(ypart q,0)
694 else :
695 (0,xpart q)--(0,ypart q)
696 fi ) shifted z) ;
697 addto pic0 doublepath p ;
698 if length pic[n]>0 :
699 if ypart q<>0 :
700 addto pic0 also pic[n] shifted point 1 of p ;
701 fi
702 if xpart q<>0 :
703 addto pic0 also pic[n] rotated 180 shifted point 0 of p ;
704 fi
705 fi
706 fi
707 fi
708 exitif incr n>3 ;
709 endfor
710 endfor
711 if length pic0>0 :
712 save bg, fg ; color bg, fg ;
713 bg := if known graph_background : graph_background else : background fi ;
714 fg := if known graph_foreground : graph_foreground else : black fi ;
715 addto graph_current_graph also pic0 withpen currentpen scaled 2 _op_ withcolor bg ;
716 addto graph_current_graph also pic0 withpen currentpen scaled .5 _op_ withcolor fg ;
717 fi
718 fi
719enddef ;
720
721
722
723
724def plot expr p =
725 if known graph_plot_picture :
726 withpen nullpen
727 hide (graph_plot_picture := image(
728 if bounded p : for q within p : graph_addto_currentpicture q endfor
729 else : graph_addto_currentpicture p
730 fi graph_setbounds origin..cycle))
731 fi
732enddef ;
733
734
735
736def graph_addto_currentpicture primary p = addto currentpicture also p ; enddef;
737def graph_setbounds = setbounds currentpicture to enddef ;
738
739def gdrawarrow = graph_number_of_arrowheads := 1 ; gdraw enddef;
740def gdrawdblarrow = graph_number_of_arrowheads := 2 ; gdraw enddef;
741
742
743
744
745
746vardef graph_post_draw =
747 save p ; path p ; p = pathpart graph_last_drawn ;
748 graph_unknown_pair_bbox(p) ;
749 if filled graph_last_drawn or not graph_is_null(penpart graph_last_drawn) :
750 addto graph_current_graph also graph_last_drawn ;
751 fi
752 graph_errorbars(graph_errorbar_text) ;
753 if length graph_plot_picture>0 :
754 for i=0 upto length p if cycle p : -1 fi :
755 addto graph_current_graph also graph_plot_picture shifted point i of p ;
756 endfor
757 picture graph_plot_picture ;
758 fi
759 if graph_number_of_arrowheads>0 :
760 graph_draw_arrowhead(p, graph_with_pen_and_color(graph_last_drawn)) ;
761 if graph_number_of_arrowheads>1 :
762 graph_draw_arrowhead(reverse p, graph_with_pen_and_color(graph_last_drawn)) ;
763 fi
764 graph_number_of_arrowheads := 0 ;
765 fi
766enddef ;
767
768vardef graph_is_null(expr p) = (urcorner p=origin) and (llcorner p=origin) enddef ;
769
770vardef graph_draw_arrowhead(expr p)(text w) =
771
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
777
778
779 graph_setbounds point infinity of p..cycle ;
780 ) ;
781enddef ;
782
783vardef graph_arrowhead_extent(expr p, q) =
784 if p<>q : (q 100ptunitvector(qp)) -- fi
785 q
786enddef ;
787
788
789
790
791
792
793
794
795
796def graph_draw_label(expr p)(suffix $)(text c) =
797 save sdim_ ; pair sdim_;
798 sdim_ := (if unknown xpart p : 1 fi 0, if unknown ypart p : 1 fi 0) ;
799 graph_unknown_pair_bbox(p) ;
800 addto graph_current_graph also
801 image(c(p) ; graph_setbounds ppsdim_cycle) _op_
802enddef ;
803
804
805
806
807vardef graph_stash_label(suffix $)(text c) text w =
808 graph_label[1.5angle mfun_laboff$ 90] = image(c(origin) w) ;
809enddef ;
810
811def graph_label_location primary p =
812 if pair p : graph_draw_label(p)
813 elseif numeric p : graph_draw_label(point p of pathpart graph_last_drawn)
814 else : graph_stash_label
815 fi
816enddef ;
817
818
819
820
821vardef glabel@#(expr p)(text t) text w =
822 graph_label_location graph_label_convert_user_to_internal(t) (@#,label@#(p)) w ; enddef;
823
824
825
826
827vardef gdotlabel@#(expr p)(text t) text w =
828 graph_label_location graph_label_convert_user_to_internal(t) (@#,dotlabel@#(p)) w ; enddef;
829
830def OUT = enddef ;
831
832
833
834
835
836
837
838path graph_template.tick, graph_template.itick, graph_template.otick, graph_template.grid ;
839graph_template.tick = (-3.5bp,0)--(3.5bp,0) ;
840graph_template.itick = origin--(7bp,0) ;
841graph_template.otick = (-7bp,0)origin ;
842graph_template.grid = origin--(1,0) ;
843
844vardef tick@#(expr f,u) text w = graph_tick_label(@#,@,false,f,u,w) ; enddef;
845
846vardef itick@#(expr f,u) text w = graph_tick_label(@#,@,false,f,u,w) ; enddef;
847
848vardef otick@#(expr f,u) text w = graph_tick_label(@#,@,false,f,u,w) ; enddef;
849
850vardef grid@#(expr f,u) text w = graph_tick_label(@#,@,true,f,u,w) ; enddef;
851
852
853
854
855
856
857
858def graph_tick_label(suffix $,$$)(expr c, f, u)(text w) =
859 graph_draw_label(graph_label_convert_user_to_internal(graph_generate_label_position($,u)),,
860 draw graph_gridline_picture$($$,c,f,u,w) shifted)
861enddef ;
862
863
864
865
866def graph_generate_label_position(suffix $)(expr u) =
867 if pair u : u elseif xpart mfun_laboff.$=0 : u,whatever else : whatever,u fi
868enddef ;
869
870
871
872
873
874
875vardef graph_gridline_picture@#(suffix $)(expr c, f, u)(text w) =
876 if unknown u : graph_error(u,"Label coordinate should be known") ; nullpicture
877 else :
878 save p ; path p;
879 interim warningcheck :=0 ;
880 graph_autogrid_needed :=false ;
881 p = graph_template$ zscaled mfun_laboff@#
882 if c : graph_xyscale fi
883 shifted (((.5 mfun_laboff@# dotprod (.5,.5)) mfun_laboff@#) graph_xyscale) ;
884 image(draw p w ;
885 label@#(if string f : format(f,u) else : f fi, point 0 of p) w)
886 fi
887enddef ;
888
889def graph_xyscale = xscaled X_.graph_dimensions yscaled Y_.graph_dimensions enddef ;
890
891
892
893vardef frame@# text w =
894 graph_frame_needed :=false ;
895 picture p_ ;
896 p_ = image(draw
897 if str@#<>"" : subpath round(angle mfun_laboff@#graph_frame_pair_agraph_frame_pair_b) of fi
898 unitsquare graph_xyscale w) ;
899 graph_draw_label((whatever,whatever),,draw p_ shifted) ;
900enddef ;
901
902pair graph_frame_pair_a ; graph_frame_pair_a=(1,1)90;
903pair graph_frame_pair_b ; graph_frame_pair_b=(.75,2.25);
904
905
906
907string graph_log_marks[] ;
908string graph_lin_marks ;
909string graph_exp_marks ;
910newinternal graph_minimum_number_of_marks, graph_log_minimum ;
911graph_minimum_number_of_marks := 4 ;
912graph_log_minimum := mlog 3 ;
913
914def Gfor(text t) = for i=t endfor enddef ;
915
916graph_log_marks[1]="1,2,5" ;
917graph_log_marks[2]="1,1.5,2,3,4,5,7" ;
918graph_log_marks[3]="1Gfor(6upto10 :,i/5)Gfor(5upto10 :,i/2)Gfor(6upto9 :,i)" ;
919graph_log_marks[4]="1Gfor(11upto20 :,i/10)Gfor(11upto25 :,i/5)Gfor(11upto19 :,i/2)" ;
920graph_log_marks[5]="1Gfor(21upto40 :,i/20)Gfor(21upto50 :,i/10)Gfor(26upto49 :,i/5)" ;
921graph_lin_marks="10,5,2" ;
922graph_exp_marks="20,10,5,2,1" ;
923
924Ten_to0 = 1 ;
925Ten_to1 = 10 ;
926Ten_to2 = 100 ;
927Ten_to3 = 1000 ;
928Ten_to4 = 10000 ;
929
930
931
932
933
934
935
936
937
938vardef graph_bounds@# =
939 interim warningcheck :=0 ;
940 save l, h ;
941 graph_set_default_bounds ;
942 if @#graph_coordinate_type>0 : (l,h) else : (h,l) fi = (@#low, @#high) ;
943 if abs @#graph_coordinate_type=logarithmic :
944 graph_modified_lower := graph_Meform(l)graph_modified_bias ;
945 graph_modified_higher := graph_Meform(h)graph_modified_bias ;
946 if hl >= graph_log_minimum : logarithmic else : linear fi
947 else :
948 graph_modified_lower := graph_Feform(l)graph_modified_bias ;
949 graph_modified_higher := graph_Feform(h)graph_modified_bias ;
950 linear
951 fi
952enddef ;
953
954pair graph_modified_bias ; graph_modified_bias=(0,3);
955pair graph_modified_lower, graph_modified_higher ;
956
957
958
959def graph_scan_marks(expr k, l, h)(text t) =
960 for m=scantokens graph_log_marks[k] :
961 exitif m>h ;
962 if m>=l : t fi
963 endfor
964enddef ;
965
966
967
968
969def graph_scan_mark(expr k, l, h)(text t) =
970 for e=ypart l upto ypart h :
971 graph_scan_marks(k, if e>ypart l : 1 else : xpart l1000 fi,
972 if e<ypart h : 10 else : xpart h1000 fi, t)
973 endfor
974enddef ;
975
976
977
978vardef graph_select_mark =
979 save k ;
980 k = 0 ;
981 forever :
982 exitif unknown graph_log_marks[k+1] ;
983 exitif 0 graph_scan_mark(incr k, graph_modified_lower, graph_modified_higher, +1)
984 >= graph_minimum_number_of_marks ;
985 endfor
986 k
987enddef ;
988
989
990
991
992vardef graph_select_exponent_mark@# =
993 numeric @# ;
994 for e=scantokens graph_exp_marks :
995 @# = e ;
996 exitif floor(ypart graph_modified_highere)
997 floor(graph_modified_exponent_ypart(graph_modified_lower)e)
998 >= graph_minimum_number_of_marks ;
999 numeric @# ;
1000 endfor
1001 known @#
1002enddef ;
1003
1004vardef graph_modified_exponent_ypart(expr p) = ypart p if xpart p=1000 : -1 fi enddef ;
1005
1006
1007
1008vardef graph_tick_mark_spacing =
1009 interim warningcheck :=0 ;
1010 save m, n, d ;
1011 m = graph_minimum_number_of_marks ;
1012 n = 1 for i=1 upto
1013 (mlog(xpart graph_modified_higherxpart graph_modified_lower) mlog m)mlogten :
1014 10 endfor ;
1015 if n<=1000 :
1016 for x=scantokens graph_lin_marks :
1017 d = nx ;
1018 exitif 0 graph_generate_numbers(d,+1)>=m ;
1019 numeric d ;
1020 endfor
1021 fi
1022 if known d : d else : n fi
1023enddef ;
1024
1025def graph_generate_numbers(expr d)(text t) =
1026 for m = dceiling(xpart graph_modified_lowerd) step d until xpart graph_modified_higher :
1027 t
1028 endfor
1029enddef ;
1030
1031
1032
1033
1034def graph_generate_exponents(expr d)(text t) =
1035 for e = dfloor(graph_modified_exponent_ypart(graph_modified_lower)d+1)
1036 step d until dfloor(ypart graph_modified_higherd) : t
1037 endfor
1038enddef ;
1039
1040
1041
1042
1043vardef graph_match_exponents =
1044 interim warningcheck := 0 ;
1045 save e ;
1046 e+3 = if graph_modified_lower=graph_modified_bias : ypart graph_modified_higher
1047 elseif graph_modified_higher=graph_modified_bias : ypart graph_modified_lower
1048 else : max(ypart graph_modified_lower, ypart graph_modified_higher) fi ;
1049 forsuffixes $=graph_modified_lower, graph_modified_higher :
1050 $ := (xpart $ for i=ypart $ upto e+2 : (10) endfor, e) ;
1051 endfor
1052 e
1053enddef ;
1054
1055
1056
1057
1058
1059vardef graph_factor_and_exponent_to_string(expr m, e) =
1060 if (e>3)or(e<-4) :
1061 decimal m & "e" & decimal e
1062 elseif e>=0 :
1063 if abs m<infinityTen_to[e] :
1064 mTen_to[e]
1065 else : decimal m & "e" & decimal e
1066 fi
1067 else :
1068 save x ; x=mTen_to[e];
1069 if abs x>=.1 : x else : decimal m & "e" & decimal e fi
1070 fi
1071enddef ;
1072
1073def auto suffix $ =
1074 hide(def graph_comma= hide(def graph_comma=,enddef) enddef)
1075 if graph_bounds.graph_suffix($)=logarithmic :
1076 if graph_select_exponent_mark.graph_exponent :
1077 graph_generate_exponents(graph_exponent,
1078 graph_comma graph_factor_and_exponent_to_string(1,e))
1079 else :
1080 graph_scan_mark(graph_select_mark, graph_modified_lower, graph_modified_higher,
1081 graph_comma graph_factor_and_exponent_to_string(m,e))
1082 fi
1083 else :
1084 hide(graph_exponent :=graph_match_exponents)
1085 graph_generate_numbers(graph_tick_mark_spacing,
1086 graph_comma graph_factor_and_exponent_to_string(m,graph_exponent))
1087 fi
1088enddef ;
1089
1090string Autoform ; Autoform = "%g";
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104vardef autogrid(suffix tx, ty) text w =
1105 graph_autogrid_needed := false ;
1106 if str tx <> "" :
1107 for x=auto.x :
1108 tx (
1109 if string Autoform_X :
1110 if Autoform_X <> "" :
1111 Autoform_X
1112 else :
1113 Autoform
1114 fi
1115 else :
1116 Autoform
1117 fi,
1118 x
1119 ) w ;
1120 endfor
1121 fi
1122 if str ty <> "" :
1123 for y=auto.y :
1124 ty (
1125 if string Autoform_Y :
1126 if Autoform_Y <> "" :
1127 Autoform_Y
1128 else :
1129 Autoform
1130 fi
1131 else :
1132 Autoform
1133 fi,
1134 y
1135 ) w ;
1136 endfor
1137 fi
1138enddef ;
1139
1140
1141
1142def endgraph =
1143 if graph_autogrid_needed : autogrid(otick.bot, otick.lft) ; fi
1144 if graph_frame_needed : frame ; fi
1145 setcoords(linear,linear) ;
1146 interim truecorners :=1 ;
1147 for b=bbox graph_finished_graph :
1148 setbounds graph_finished_graph to b ;
1149 for i=0 step .5 until 3.5 :
1150 if known graph_label[i] :
1151 addto graph_finished_graph also graph_label[i] shifted point i of b ;
1152 fi
1153 endfor
1154 endfor
1155 graph_finished_graph
1156 endgroup
1157enddef ;
1158
1159
1160
1161
1162
1163
1164if known contextlmtxmode :
1165
1166elseif unknown context_mlib :
1167
1168 vardef escaped_format(expr s) =
1169 "" for n=0 upto length(s) : &
1170 if ASCII substring (n,n+1) of s = 37 :
1171 "@"
1172 else :
1173 substring (n,n+1) of s
1174 fi
1175 endfor
1176 enddef ;
1177
1178 vardef strfmt(expr f, x) =
1179 "\MPgraphformat{" & escaped_format(f) & "}{" & mfun_tagged_string(x) & "}"
1180 enddef ;
1181
1182 vardef varfmt(expr f, x) =
1183 "\MPformatted{" & escaped_format(f) & "}{" & mfun_tagged_string(x) & "}"
1184 enddef ;
1185
1186 vardef format (expr f, x) = textext(strfmt(f,x)) enddef ;
1187 vardef formatted(expr f, x) = textext(varfmt(f,x)) enddef ;
1188
1189fi ;
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199vardef graph_shapesize = (.33BodyFontSize) enddef ;
1200
1201path graph_shape[] ;
1202
1203graph_shape[0] := (0,0) ;
1204graph_shape[1] := fullcircle ;
1205graph_shape[2] := (up -- down) scaled .5 ;
1206
1207for i = 3 upto 9 :
1208 graph_shape[i] :=
1209 for j = 0 upto i-1 :
1210 (up scaled .5) rotated (360ji) --
1211 endfor cycle ;
1212endfor
1213
1214graph_shape[12] := graph_shape[2] rotated +90 ;
1215graph_shape[22] := graph_shape[2] rotated +45 ;
1216graph_shape[32] := graph_shape[2] rotated -45 ;
1217graph_shape[13] := graph_shape[3] rotated 180 ;
1218graph_shape[23] := graph_shape[3] rotated -90 ;
1219graph_shape[33] := graph_shape[3] rotated +90 ;
1220graph_shape[14] := graph_shape[4] rotated +45 ;
1221graph_shape[15] := graph_shape[5] rotated 180 ;
1222graph_shape[16] := graph_shape[6] rotated +90 ;
1223graph_shape[17] := graph_shape[7] rotated 180 ;
1224graph_shape[18] := graph_shape[8] rotated +22.5 ;
1225
1226numeric l ;
1227
1228for j = 5 upto 9 :
1229 l := length(graph_shape[j]) ;
1230 pair p[] ;
1231 for i = 0 upto l :
1232 p[i] = whatever [point i of graph_shape[j],
1233 point (i+2 mod l) of graph_shape[j]] ;
1234 p[i] = whatever [point (i+1 mod l) of graph_shape[j],
1235 point (il-1 mod l) of graph_shape[j]] ;
1236 endfor
1237 graph_shape[20j] := for i = 0 upto l : point i of graph_shape[j]p[i]endfor cycle ;
1238endfor
1239
1240path s ; s := graph_shape[4] ;
1241path q ; q := s scaled .25 ;
1242numeric l ; l := length(s) ;
1243
1244pair p[] ;
1245
1246graph_shape[24] := for i = 0 upto l-1 :
1247 hide(
1248 p[i] = whatever [point i of s, point (i+1 mod l) of s] ;
1249 p[i] = whatever [point i of q, point (i-1l mod l) of q] ;
1250 p[il] = whatever [point i of s, point (i+1 mod l) of s] ;
1251 p[il] = whatever [point i+1 of q, point (i+2 mod l) of q] ;
1252 )
1253 point i of q -- p[i] -- p[il] --
1254endfor cycle ;
1255
1256graph_shape[34] := graph_shape[24] rotated 45 ;
1257
1258
1259
1260
1261
1262def stars(expr f) = plotsymbol(25,f) enddef ;
1263def points(expr f) = plotsymbol( 0,f) enddef ;
1264def circles(expr f) = plotsymbol( 1,f) enddef ;
1265def crosses(expr f) = plotsymbol(34,f) enddef ;
1266def squares(expr f) = plotsymbol(14,f) enddef ;
1267def diamonds(expr f) = plotsymbol( 4,f) enddef ;
1268def uptriangles(expr f) = plotsymbol( 3,f) enddef ;
1269def downtriangles(expr f) = plotsymbol(13,f) enddef ;
1270def lefttriangles(expr f) = plotsymbol(33,f) enddef ;
1271def righttriangles(expr f) = plotsymbol(23,f) enddef ;
1272
1273
1274def plotsymbol(expr n, f) text t =
1275 if known graph_shape[n] :
1276 image(
1277 save bg, fg ; color bg, fg ;
1278 bg := if known graph_background : graph_background else : background fi ;
1279 save pic ; picture pic ; pic := image(draw origin _op_ t ;) ;
1280 if color colorpart pic : graph_foreground := colorpart pic ; fi
1281 fg := if known graph_foreground : graph_foreground else : black fi ;
1282 save p ; path p ; p = graph_shape[n] scaled graph_shapesize ;
1283 draw p withcolor bg withpen currentpen scaled 2 ;
1284 if cycle p :
1285 fill p withcolor
1286 if known f :
1287 if color f :
1288 f
1289 elseif numeric f :
1290 f[bg,fg]
1291 elseif boolean f and f :
1292 fg
1293 else
1294 bg
1295 fi
1296 else :
1297 bg
1298 fi ;
1299 fi
1300 draw p _op_ t ;
1301 )
1302 else :
1303 nullpicture
1304 fi
1305 t
1306enddef ;
1307
1308
1309color resistance_color[] ; string resistance_name[] ;
1310resistance_color0 = (0,0,0) ; resistance_name0 = "black" ;
1311resistance_color1 = (165255,42255,42255) ; resistance_name1 = "brown" ;
1312resistance_color2 = (1,0,0) ; resistance_name2 = "red" ;
1313resistance_color3 = (1,165255,0) ; resistance_name3 = "orange" ;
1314resistance_color4 = (1,1,0) ; resistance_name4 = "yellow" ;
1315resistance_color5 = (0,1,0) ; resistance_name5 = "green" ;
1316resistance_color6 = (0,0,1) ; resistance_name6 = "blue" ;
1317resistance_color7 = (148255,0,211255) ; resistance_name7 = "darkviolet" ;
1318resistance_color8 = (190255,190255,190255) ; resistance_name8 = "gray" ;
1319resistance_color9 = (1,1,1) ; resistance_name9 = "white" ;
1320
1321
1322
1323
1324
1325
1326def rainbow(expr f) =
1327 hide(numeric n_ ; n_ = (abs(5f) mod 5) 2 ;)
1328 (n_floor(n_))[resistance_color[floor n_],resistance_color[ceiling n_]]
1329enddef ;
1330
1331
1332
1333
1334
1335vardef sortpath (suffix $) (text t) =
1336 if path $ :
1337 if length $ > 0 :
1338 save n, k ; n := length $ ;
1339 for i=0 upto n :
1340 k := i ;
1341 for j=i+1 upto n :
1342 if t (point j of $) < t (point k of $) :
1343 k := j ;
1344 fi
1345 endfor
1346 if k>i :
1347 $ := if i>0 : subpath (0,i-1) of $ -- fi
1348 point k of $ --
1349 subpath (i,k-1) of $
1350 if k<n : -- subpath (k+1,n) of $ fi
1351 ;
1352 fi
1353 endfor
1354 fi
1355 fi
1356enddef ;
1357
1358
1359
1360def smoothpath (suffix $) =
1361 if path $ :
1362 (for i=0 upto length $ :
1363 if i>0 : .. fi
1364 (point i of $)
1365 endfor )
1366 fi
1367enddef ;
1368
1369
1370
1371def makefunctionpath (expr f, t, n) (text func) =
1372 (for x=f step ((tf)(abs n)) until t :
1373 if x<>f : -- fi
1374 (x, func)
1375 endfor )
1376enddef ;
1377
1378
1379
1380
1381
1382
1383
1384vardef addtopath (suffix p) (text t) =
1385 if path p :
1386 (for i=0 upto length p :
1387 if i>0 : -- fi
1388 hide(clearxy ; z = point i of p ;) z shifted t
1389 endfor)
1390 fi
1391enddef ;
1392
1393
1394
1395vardef functionpath (suffix p) (text func) =
1396 (for i=0 upto length p :
1397 if i>0 : .. fi
1398 (hide(x := xpart(point i of p))x,func)
1399 endfor )
1400enddef ;
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425vardef polynomial_function (suffix $) (expr n, x) =
1426 (for j=0 upto n : $[j](xj) endfor)
1427enddef ;
1428
1429
1430
1431vardef det (suffix $) (expr n) =
1432 hide(
1433 numeric determinant ; determinant := 1 ;
1434 save jj ; numeric jj ;
1435 for k=0 upto n :
1436 if $[k][k]=0 :
1437 jj := -1 ;
1438 for j=0 upto n :
1439 if $[k][j]<>0 :
1440 jj := j ;
1441 exitif true ;
1442 fi
1443 endfor
1444 if jj<0 :
1445 determinant := 0 ;
1446 exitif true ;
1447 fi
1448 for j=k upto n :
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 :
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
1467enddef ;
1468
1469numeric fit_chi_squared ;
1470
1471
1472
1473
1474
1475
1476vardef polynomial_fit (suffix p, $) (expr n) (text t) =
1477 if not path p :
1478 graph_error(p, "Cannot fit--not a path") ;
1479 elseif length p < n :
1480 graph_error(p, "Cannot fit--not enough points") ;
1481 else :
1482 fit_chi_squared := 0 ;
1483
1484 save sumx, sumy ; numeric sumx[], sumy[] ;
1485 save w ; numeric w ;
1486 for i=0 upto 2n :
1487 sumx[i] := 0 ;
1488 endfor
1489 for i=0 upto n :
1490 sumy[i] := 0 ;
1491 endfor
1492 for i=0 upto length p :
1493 clearxy ; z = point i of p ;
1494 w := 1 ;
1495 if known t :
1496 if numeric t :
1497 w := 1 if t<>0 : (abs t) fi ;
1498 elseif pair t :
1499 if t<>origin :
1500 w := 1(abs t) ;
1501 fi
1502 elseif path t :
1503 if length t>= i:
1504 if point i of t<>origin :
1505 w := 1(abs point i of t) ;
1506 fi
1507 else :
1508 w := 0 ;
1509 fi ;
1510 fi
1511 fi
1512 x1 := w ;
1513 for j=0 upto 2n :
1514 sumx[j] := sumx[j] x1 ;
1515 x1 := x1 x ;
1516 endfor
1517 y1 := y w ;
1518 for j=0 upto n :
1519 sumy[j] := sumy[j] y1 ;
1520 y1 := y1 x ;
1521 endfor
1522 fit_chi_squared := fit_chi_squared yyw ;
1523 endfor
1524
1525 save m ; numeric m[][] ;
1526 for j=0 upto n :
1527 for k=0 upto n :
1528 m[j][k] := sumx[jk] ;
1529 endfor
1530 endfor
1531 save delta ; numeric delta ;
1532 delta := det(m,n) ;
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[jk] ;
1543 endfor
1544 m[j][i] := sumy[j] ;
1545 endfor
1546 $[i] := det(m,n) delta ;
1547 endfor
1548 for j=0 upto n :
1549 fit_chi_squared := fit_chi_squared 2sumy[j]$[j] ;
1550 for k=0 upto n :
1551 fit_chi_squared := fit_chi_squared $[j]$[k]sumx[jk] ;
1552 endfor
1553 endfor
1554
1555 fit_chi_squared := fit_chi_squared (length(p) n) ;
1556 fi
1557 fi
1558enddef ;
1559
1560
1561
1562
1563
1564vardef linear_function (suffix $) (expr x) = polynomial_function($,1,x) enddef ;
1565vardef linear_fit (suffix p, $) (text t) = polynomial_fit(p, $, 1, t) ; enddef ;
1566
1567
1568
1569vardef constant_function (suffix $) (expr x) = polynomial_function($,0,x) enddef ;
1570vardef constant_fit (suffix p, $) (text t) = polynomial_fit(p, $, 0, t) ; enddef ;
1571
1572
1573
1574
1575
1576vardef exponential_function (suffix $) (expr x) = $1exp($0x) enddef ;
1577
1578
1579
1580vardef exponential_fit (suffix p, $) (text t) =
1581 save a ; numeric a[] ;
1582 save q ; path q[] ;
1583 for i=0 upto length p :
1584 clearxy ; z = point i of p ;
1585 if y>0 :
1586 augment.q0(x,ln(y)) ;
1587 augment.q1(
1588 if known t :
1589 if numeric t : (0,ln(t))
1590 elseif pair t : (xpart t,ln(ypart t))
1591 elseif path t :
1592 if length t>=i :
1593 hide(z1 = point i of t;)
1594 (x1,ln(y1))
1595 else :
1596 origin
1597 fi
1598 fi
1599 else :
1600 (0,1)
1601 fi ) ;
1602 fi
1603 endfor
1604 linear_fit(q0,a,q1) ;
1605 save e ; e := exp(sqrt(fit_chi_squared)) ;
1606 fit_chi_squared := e e ;
1607 $0 := a1 ;
1608 $1 := exp(a0) ;
1609enddef ;
1610
1611
1612
1613vardef power_law_function (suffix $) (expr x) = $1(x$0) enddef ;
1614
1615
1616
1617vardef power_law_fit (suffix p, $) (text t) =
1618 save a ; numeric a[] ;
1619 save q ; path q[] ;
1620 for i=0 upto length p :
1621 clearxy ; z = point i of p ;
1622 if (x>0) and (y>0) :
1623 augment.q0(ln(x),ln(y)) ;
1624 augment.q1(
1625 if known t :
1626 if numeric t : (0,ln(t))
1627 elseif pair t : (ln(xpart t),ln(ypart t))
1628 elseif path t :
1629 if length t>=i :
1630 hide(z1 = point i of t)
1631 (ln(x1),ln(y1))
1632 else :
1633 origin
1634 fi
1635 fi
1636 else :
1637 (0,1)
1638 fi ) ;
1639 fi
1640 endfor
1641 linear_fit(q0,a,q1) ;
1642 save e ; e := exp(sqrt(fit_chi_squared)) ;
1643 fit_chi_squared := e e ;
1644 $0 := a1 ;
1645 $1 := exp(a0) ;
1646enddef ;
1647
1648
1649
1650
1651
1652newinternal lntwo ; lntwo := ln(2) ;
1653
1654vardef gaussian_function (suffix $) (expr x) =
1655 if $1 = 0 :
1656 if x = $0 : $2 else : 0 fi
1657 else :
1658 $2 exp(lntwo(((x$0)$1)2))
1659 fi
1660 if known $3 :
1661 $3
1662 fi
1663enddef ;
1664
1665
1666
1667vardef gaussian_fit (suffix p, $) (text t) =
1668 save a ; numeric a[] ;
1669 save q ; path q[] ;
1670 for i=0 upto length p :
1671 clearxy ; z = point i of p ;
1672 if y>0 :
1673 augment.q0(x,ln(y)) ;
1674 augment.q1(
1675 if known t :
1676 if numeric t : (0,ln(t))
1677 elseif pair t : (xpart t,ln(ypart t))
1678 elseif path t :
1679 if length t>=i :
1680 hide(z1 = point i of t)
1681 (x1,ln(y1))
1682 else :
1683 origin
1684 fi
1685 fi
1686 else :
1687 (0,1)
1688 fi ) ;
1689 fi
1690 endfor
1691 polynomial_fit(q0,a,2,q1) ;
1692 save e ; e := exp(sqrt(fit_chi_squared)) ;
1693 fit_chi_squared := e e ;
1694 $1 := sqrt(lntwoa2) ;
1695 $0 := -.5a1a2 ;
1696 $2 := exp(a0-.25a1a1a2) ;
1697 $3 := 0 ;
1698enddef ;
1699
1700
1701
1702vardef lorentzian_function (suffix $) (expr x) =
1703 if $1 = 0 :
1704 if x = $0 : $2 else : 0 fi
1705 else :
1706 $2 (1 ((x $0)$1)2)
1707 fi
1708 if known $3 :
1709 $3
1710 fi
1711enddef ;
1712
1713vardef lorentzian_fit (suffix p, $) (text t) =
1714 save a ; numeric a[] ;
1715 save q ; path q ;
1716 for i=0 upto length p :
1717 if ypart(point i of p)<>0 :
1718 augment.q(xpart(point i of p), 1ypart(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 := 1fit_chi_squared ;
1723 $0 := -.5a1a2 ;
1724 $2 := 1(a0-.25a1a1a2) ;
1725 $1 := sqrt((a0-.25a1a1a2)a2) ;
1726 $3 := 0 ;
1727enddef ;
1728 |