%!PS
% PS Data Point Power Curve Fitter
% =================================
% by Don Lancaster
% Copyright c 2003 by Don Lancaster & Synergetics, Box 809, Thatcher, AZ, 85552
% (928) 428-4073 Email: don@tinaja.com Website: http://www.tinaja.com
% Consulting services available http://www.tinaja.com/info01.html
% All commercial rights and all electronic media rights ~fully~ reserved.
% Linking usually welcome. Reposting expressly forbidden. Version 1.1
% PostScript-as-language tutorial and utilites fit power series curves to 3, 4,
% 5, or 6 data points. Uses Gauss-Jordan elimination. Extendable to more points.
% Modify in editor or wp. Then send to distiller.
% Additional tutorial info at http://www.tinaja.com/gurgrm01.asp
% =========
% This file requires the previous download of gonzo.psl
% available from https://www.tinaja.com/pssamp1.shtml
% Make sure the following line agrees with your own gonzo.psl location
(C:/Users/don/Desktop/gonzo/gonzo.psl) run % use internal gonzo
% ==========
/guru { gonzo begin
ps.util.1 begin printerror nuisance begin} def
guru % activate gonzo utilities
% =========
[/CropBox [0 0 400 400] % set the wierd size
/PAGES pdfmark
[ {Catalog} <<
/ViewerPreferences << /FitWindow true /CenterWindow true >>
/PageLayout /OneColumn % continous
/Pagemode /UseNone % no thumbs
/View [/XYZ null null 1 ] % force 100%
>>
/PUT pdfmark
%%%%%%%%%%% QUADRATIC THREE-POINT FIT %%%%%%%%%%%%%%
save /snap1 exch store
% Fit a quadratic to three data points
% points are 0,0 xx1,yy1, and 1,1
% seek general solution.
% y = ax2 + bx + c There is no c since y=0 at x=0
% 1 = a + b at x = 1 EQUATION AA
% a(xx1)2 + b(xx1) = (yy1) % evaluating remaining point
% a + b/(xx1) = (yy1)/(xx1)2
% 1 - b + b/(xx1) = (yy1)/(xx1)2
% - b + b/(xx1) = (yy1)/(xx1)2 - 1
% b ((1 - (xx1))/(xx1) = (yy1)/(xx1)2 - 1
% b ((1 - (xx1))/(xx1) = (yy1)/(xx1)2 - (xx1)/(xx1)
% b ((1 - (xx1)) = (yy1)/(xx1) - (xx1)
/bb {y1 x1 div x1 sub 1 x1 sub div} store
/aa {1 bb sub} store
% ========== quadratic demo plot =====
/x1 0.5 store
/y1 0.75 store
% in this example, a y1 > 0.75 overshoots. No overshoot if aa < 1.
100 100 10 setgrid
20 20 showgrid
0 0 mt dot
1 20 mul 1 20 mul mt dot
x1 20 mul y1 20 mul mt dot
line1 0 0 mt 20 20 lineto stroke
line2
0 0 moveto
0 0.01 1 { /xx exch store
xx 20 mul
xx dup mul aa mul
xx bb mul add
20 mul lineto
} for
stroke
showpage
snap1 restore
%%%%%%%%%%% CUBIC FOUR-POINT FIT %%%%%%%%%%%%%%
save /snap2 exch store
% fit a cubic to four data points
% seek general solution.
% ax3 + bx2 + cx + d no d since y=0 at x=0
% 1 = a + b + c at x = 1 EQUATION AA
% a(j)3 + b(j)2 + c(j) = u evaluating first point Eqn U1
% a(k)3 + b(k)2 + b(k) = v evaluating second point Eqn V1
% scale to clean a
% a + b/j + c/j2 = u/j3 eqn u1a
% a + b/k + c/k2 = v/k3 eqn v1a
% substitute AA to get rid of a
% 1 - b - c + b/j + c/j2 = u/j3
% 1 - b - c + b/k + c/k2 = v/k3
% b ( 1/j - 1) + c (1/j2 - 1) = u/j3 - 1 rearranging
% c ( 1/k - 1) + c (1/k2 - 1) = v/k3 - 1
% b ((1-j)/j) + c ((1-j2)/j2) = (u - j3)/j3 denomenatoring
% c ((1-k)/k) + c ((1-k2)/k2) = (v - k3)/k3
% b (1-j) + c (1-j2)/j = (u - j3)/j2 multiplying by j or k
% b (1-k) + c (1-k2)/k = (v - k3)/k2
% b + c (1-j2)/j(1-j) = (u-j3)/j2(1-j) scaling by 1/ 1-j
% b + c (1-k2)/k(1-k) = (v-k3)/k2(1-k) scaling by 1/ 1-k
% m = (1-j2)/j(1-j) % simplify variables
% n = (u-j3)/j2(1-j)
% p = (1-k2)/k(1-k)
% q = (v-k3)/k2(1-k)
% b + cm = n substitute for sanity
% b + cp = q
% cm - cp = n - q % subtracting the equations to cancel b
% c (m-p) = (n-q)
% c = (n-q)/(m-p) < CCCCCCCC
% b + ((n-q)/(m-p))m = n % and resolving for b
% b = n - ((n-q)/(m-p))m < BBBBBBBBB
% a = 1 - n + ((n-q)/(m-p))m - (n-q)/(m-p) < --- AAAAAA
% m = (1-j2)/j(1-j)
% n = (u-j3)/j2(1-j)
% p = (1-k2)/k(1-k)
% q = (v-k3)/k2(1-k)
/mm { 1 jj dup mul sub jj 1 jj sub mul div} store
/pp { 1 kk dup mul sub kk 1 kk sub mul div} store
/nn { uu jj dup dup mul mul sub jj dup mul 1 jj sub mul div} store
/qq { vv kk dup dup mul mul sub kk dup mul 1 kk sub mul div} store
% c = (n-q)/(m-p)
% b = n - ((n-q)/(m-p))m < BBBBBBBBB
% a = 1 - n + ((n-q)/(m-p))m - (n-q)/(m-p) < --- AAAAAA
/cc {nn qq sub mm pp sub div} store
/bb { nn nn qq sub mm pp sub div mm mul sub } store
/aa { 1 bb sub cc sub } store
% ========== cubic demo plot =====
100 100 10 setgrid
20 20 showgrid
/jj 0.3333 store
/kk 0.6667 store
/delta 0.15 store
/uu {jj delta sub} store
/vv {kk delta sub} store
0 0 mt dot
1 20 mul 1 20 mul mt dot
% jj 20 mul uu 20 mul mt dot
% kk 20 mul vv 20 mul mt dot
line1 0 0 mt 20 20 lineto stroke
line2
0 0.05 0.5 {/delta exch store
jj 20 mul uu 20 mul mt dot
kk 20 mul vv 20 mul mt dot
0 0 moveto
0 0.01 1 { /xx exch store
xx 20 mul
xx dup dup mul mul aa mul
xx dup mul bb mul add
xx cc mul add
20 mul lineto
} for
stroke
} for
stroke
showpage
snap2 restore
%%%%%%%%%%% QUARTIC FIVE-POINT FIT %%%%%%%%%%%%%%
save /snap3 exch store
% fit a quartic to five data points
% points are 0,0 x1,y1, x2,y2, x3,y3 and 1,1
% seek general solution.
% y = ax4 + bx3 + cx2 + dx + e no e since y=0 at x=0
% 1 = a + b + c + d at x = 1 EQUATION AA
% a = 1 - b - c - d % egn I
% a(x1)4 + b(x1)3 + c(x1)2 + d(x1) = (y1) % evaluating first point
% a(x2)4 + b(x2)3 + c(x2)2 + d(x2) = (y2) % evaluating second point
% a(x3)4 + b(x3)3 + c(x3)2 + d(x3) = (y3) % evaluating third point
% scale to clean a
% a + b/(x1) + c/(x1)2 + d(x1)3 = (y1)/(x1)4 % divide by (x1)4
% a + b/(x2) + c/(x2)2 + d(x2)3 = (y2)/(x2)4 % divide by (x2)4
% a + b/(x3) + c/(x3)2 + d(x3)3 = (y1)/(x3)4 % divide by (x3)4
% reduce a with eqn I
% 1 - b - c - d + b/(x1) + c/(x1)2 + d(x1)3 = (y1)/(x1)4 % substituting for a
% 1 - b - c - d + b/(x2) + c/(x2)2 + d(x2)3 = (y2)/(x1)4 % substituting for a
% 1 - b - c - d + b/(x3) + c/(x3)2 + d(x3)3 = (y3)/(x1)4 % substituting for a
% rearranging
% b (1-(x1))/(x1) + c (1-(x1)2)/(x1)2 + d (1-(x1)3)/(x1)3 = (y1)/(x1)4 - 1
% b (1-(x2))/(x2) + c (1-(x2)2)/(x2)2 + d (1-(x2)3)/(x2)3 = (y2)/(x2)4 - 1
% b (1-(x3))/(x3) + c (1-(x3)2)/(x3)2 + d (1-(x3)3)/(x3)3 = (y3)/(x3)4 - 1
% simplifying variables
% f1 = (1-(x1))/(x1)
% f2 = (1-(x1)2)/(x1)2
% f3 = (1-(x1)3)/(x1)3
% f4 = (y1)/(x1)4 - 1
/f1 { 1 x1 sub x1 div } store
/f2 { 1 x1 dup mul sub x1 dup mul div } store
/f3 { 1 x1 dup dup mul mul sub x1 dup dup mul mul div } store
/f4 { y1 x1 dup mul dup mul div 1 sub } store
% g1 = (1-(x2))/(x2)
% g2 = (1-(x2)2)/(x2)2
% g3 = (1-(x2)3)/(x2)3
% g4 = (y2)/(x2)4 - 1
/g1 { 1 x2 sub x2 div } store
/g2 { 1 x2 dup mul sub x2 dup mul div } store
/g3 { 1 x2 dup dup mul mul sub x2 dup dup mul mul div } store
/g4 { y2 x2 dup mul dup mul div 1 sub } store
% h1 = (1-(x3))/(x3)
% h2 = (1-(x3)2)/(x3)2
% h3 = (1-(x3)3)/(x3)3
% h4 = (y3)/(x3)4 - 1
/h1 { 1 x3 sub x3 div } store
/h2 { 1 x3 dup mul sub x3 dup mul div } store
/h3 { 1 x3 dup dup mul mul sub x3 dup dup mul mul div } store
/h4 { y3 x3 dup mul dup mul div 1 sub } store
% b(f1) + c(f2) + d(f3) = (f4)
% b(g1) + c(g2) + d(g3) = (g4)
% b(h1) + c(h2) + d(h3) = (h4)
% unify b
% b + c(f2/f1) + d(f3/f1) = (f4)/(f1)
% b + c(g2/g1) + d(g3/g1) = (g4)/(g1)
% b + c(h2/h1) + d(h3/h1) = (h4)/(h1)
% simplify variables
% p1 = f2/f1
% p2 = f3/f1
% p3 = f4/f1
/p1 { f2 f1 div} store
/p2 { f3 f1 div} store
/p3 { f4 f1 div} store
% q1 = g2/g1
% q2 = g3/g1
% q3 = g4/g1
/q1 { g2 g1 div} store
/q2 { g3 g1 div} store
/q3 { g4 g1 div} store
% r1 = h2/h1
% r2 = h3/h1
% r3 = h4/h1
/r1 { h2 h1 div} store
/r2 { h3 h1 div} store
/r3 { h4 h1 div} store
% b + c(p1) + d(p2) = (p3) % eqn II
% b + c(q1) + d(q2) = (q3) % eqn III
% b + c(r1) + d(r2) = (r3) % eqn IV
% subtract eqn III from II and eqn IV from III
% c((p1)-(q1)) + d((p2)-(q2) = (p3)-(q3)
% c((q1)-(r1)) + d((q2)-(r2) = (q3)-(r3)
% simplify variables
% s1 = (p1) - (q1)
% s2 = (p2) - (q2)
% s3 = (p3) - (q3)
/s1 { p1 q1 sub} store
/s2 { p2 q2 sub} store
/s3 { p3 q3 sub} store
% t1 = (q1) - (r1)
% t2 = (q2) - (r2)
% t3 = (q3) - (r3)
/t1 { q1 r1 sub} store
/t2 { q2 r2 sub} store
/t3 { q3 r3 sub} store
% c(s1) + d(s2) = (s3)
% c(t1) + d(t2) = (t3)
% unify c
% c + d (s2)/(s1) = (s3)/(s1)
% c + d (t2)/(t1) = (t3)/(t1)
% simplify variables
% u1 = (s2)/(s1)
% u2 = (s3)/(s1)
/u1 { s2 s1 div} store
/u2 { s3 s1 div} store
% v1 = (t2)/(t1)
% v2 = (t3)/(t1)
/v1 { t2 t1 div} store
/v2 { t3 t1 div} store
% c + d(u1) = (u2) % eqn V
% c + d(v1) = (v2) % eqn VI
% subtract VI from V
% d ((u1) - (v1)) = (u2) - (v2)
% AND SOLVE FOR DDD
% d = (u2) - (v2)/((u1) - (v1)) % DDDDDDDDD
/dd { u2 v2 sub u1 v1 sub div } store
% find c from eqn V
% c + d(u1) = (u2)
% c = (u2) - d(u1) <------- CCCCC
/cc {u2 dd u1 mul sub} store
% find b from eqn II
% b + c(p1) + d(p2) = (p3)
% b = (p3) - c(p1) - d(p2) <------- BBBBB
/bb { p3 cc p1 mul sub dd p2 mul sub } store
% and a from eqn I
% a = 1 - b - c - d
/aa {1 bb sub cc sub dd sub} store
% ==========
/x1 0.2 store
/y1 0.1 store
/x2 0.4 store
/y2 0.3 store
/x3 0.7 store
/y3 0.8 store
100 100 10 setgrid
20 20 showgrid
0 0 mt dot
1 20 mul 1 20 mul mt dot
x1 20 mul y1 20 mul mt dot
x2 20 mul y2 20 mul mt dot
x3 20 mul y3 20 mul mt dot
line1 0 0 mt 20 20 lineto stroke
line2
0 0 moveto
0 0.01 1 { /xx exch store
xx 20 mul
xx dup mul dup mul aa mul
xx dup dup mul mul bb mul add
xx dup mul cc mul add
xx dd mul add
20 mul lineto
} for
stroke
showpage
snap3 restore
%%%%%%%%%%% QUINTIC SIX-POINT FIT %%%%%%%%%%%%%%
save /snap4 exch store
% fit a quintic to six data points
% points are 0,0 xx1,yy1, xx2,yy2, xx3,yy3 xx4,yy4, and 1,1
% seek general solution.
% y = ax5 + bx4 + cx3 + dx2 + ex + f no f since y=0 at x=0
% 1 = a + b + c + d + e at x = 1 EQUATION AA
% a = 1 - b - c - d - e % egn I
% a(x1)5 + b(x1)4 + c(x1)3 + d(x1)2 + e(x1) = (y1) % evaluating first point
% a(x2)5 + b(x2)4 + c(x2)3 + d(x2)2 + e(x2) = (y2) % second point
% a(x3)5 + b(x3)4 + c(x3)3 + d(x3)2 + e(x3) = (y3) % third point
% a(x4)5 + b(x4)4 + c(x4)3 + d(x4)2 + e(x4) = (y4) % third point
% .00032 a + .0016 b + .008 c + .04 d + .2000 e = .04 at x = .2
% .01024 a + .0256 b + .064 c + .16 d + .4000 e = .16 at x = .4
% .07776 a + .1296 b + .096 c + .36 d + .6000 e = .36 at x = .6
% .32768 a + .4096 b + .512 c + .64 d + .8000 e = .64 at x = 0.8
% scale to clean a
% a + b/(x1) + c/(x1)2 + d/(x1)3 + e/(x1)4 = (y1)/(x1)5 % divide by (x1)4
% a + b/(x2) + c/(x2)2 + d/(x2)3 + e/(x2)4 = (y2)/(x2)5 % divide by (x2)4
% a + b/(x3) + c/(x3)2 + d/(x3)3 + e/(x3)4 = (y3)/(x3)5 % divide by (x3)4
% a + b/(x4) + c/(x4)2 + d/(x4)3 + e/(x4)4 = (y4)/(x4)5 % divide by (x4)4
% reduce a with eqn I
% 1 - b - c - d - e + b/(x1) + c/(x1)2 + d/(x1)3 + e/(x1)4 = (y1)/(x1)5
% 1 - b - c - d - e + b/(x2) + c/(x2)2 + d/(x2)3 + e/(x2)4 = (y2)/(x2)5
% 1 - b - c - d - e + b/(x3) + c/(x3)2 + d/(x3)3 + e/(x3)4 = (y3)/(x3)5
% 1 - b - c - d - e + b/(x4) + c/(x4)2 + d/(x4)3 + e/(x4)4 = (y4)/(x4)5
% rearranging
% b (1-(x1))/(x1) + c (1-(x1)2)/(x1)2 + d (1-(x1)3)/(x1)3 + e(1-(x1)4)/(x1)4
% = (y1)/(x1)5 - 1
% b (1-(x2))/(x2) + c (1-(x2)2)/(x2)2 + d (1-(x2)3)/(x2)3 + e(1-(x2)4)/(x2)4
% = (y2)/(x2)5 - 1
% b (1-(x3))/(x3) + c (1-(x3)2)/(x3)2 + d (1-(x3)3)/(x3)3 + e(1-(x3)4)/(x3)4
% = (y3)/(x3)5 - 1
% b (1-(x4))/(x4) + c (1-(x4)2)/(x4)2 + d (1-(x4)3)/(x4)3 + e(1-(x4)4)/(x4)4
% = (y4)/(x4)5 - 1
% simplifying variables
% f1 = (1-(x1))/(x1)
% f2 = (1-(x1)2)/(x1)2
% f3 = (1-(x1)3)/(x1)3
% f4 = (1-(x1)4)/(x1)4
% f5 = (y1)/(x1)5 - 1
/f1 { 1 x1 sub x1 div } store
/f2 { 1 x1 dup mul sub x1 dup mul div } store
/f3 { 1 x1 dup dup mul mul sub x1 dup dup mul mul div } store
/f4 { 1 x1 dup mul dup mul sub x1 dup mul dup mul div } store
/f5 { y1 x1 dup dup mul dup mul mul div 1 sub } store
% g1 = (1-(x2))/(x2)
% g2 = (1-(x2)2)/(x2)2
% g3 = (1-(x2)3)/(x2)3
% g4 = (1-(x2)4)/(x2)4
% g5 = (y2)/(x2)5 - 1
/g1 { 1 x2 sub x2 div } store
/g2 { 1 x2 dup mul sub x2 dup mul div } store
/g3 { 1 x2 dup dup mul mul sub x2 dup dup mul mul div } store
/g4 { 1 x2 dup mul dup mul sub x2 dup mul dup mul div } store
/g5 { y2 x2 dup dup mul dup mul mul div 1 sub } store
% h1 = (1-(x3))/(x3)
% h2 = (1-(x3)2)/(x3)2
% h3 = (1-(x3)3)/(x3)3
% h4 = (1-(x3)4)/(x3)4
% h5 = (y3)/(x3)5 - 1
/h1 { 1 x3 sub x3 div } store
/h2 { 1 x3 dup mul sub x3 dup mul div } store
/h3 { 1 x3 dup dup mul mul sub x3 dup dup mul mul div } store
/h4 { 1 x3 dup mul dup mul sub x3 dup mul dup mul div } store
/h5 { y3 x3 dup dup mul dup mul mul div 1 sub } store
% i1 = (1-(x4))/(x4)
% i2 = (1-(x4)2)/(x4)2
% i3 = (1-(x4)3)/(x4)3
% i4 = (1-(x4)4)/(x4)4
% i5 = (y4)/(x4)5 - 1
/i1 { 1 x4 sub x4 div } store
/i2 { 1 x4 dup mul sub x4 dup mul div } store
/i3 { 1 x4 dup dup mul mul sub x4 dup dup mul mul div } store
/i4 { 1 x4 dup mul dup mul sub x4 dup mul dup mul div } store
/i5 { y4 x4 dup dup mul dup mul mul div 1 sub } store
% b(f1) + c(f2) + d(f3) + e(f4) = (f5)
% b(g1) + c(g2) + d(g3) + e(g4) = (g5)
% b(h1) + c(h2) + d(h3) + e(h4) = (h5)
% b(i1) + c(i2) + d(i3) + e(i4) = (i5)
% unify b
% b + c(f2/f1) + d(f3/f1) + e(f4/f1) = (f5)/(f1)
% b + c(g2/g1) + d(g3/g1) + e(g4/g1) = (g5)/(g1)
% b + c(h2/h1) + d(h3/h1) + e(h4/h1) = (h5)/(h1)
% b + c(i2/i1) + d(i3/i1) + e(i4/i1) = (i5)/(i1)
% simplify variables
% p1 = f2/f1
% p2 = f3/f1
% p3 = f4/f1
% p4 = f5/f1
/p1 { f2 f1 div} store
/p2 { f3 f1 div} store
/p3 { f4 f1 div} store
/p4 { f5 f1 div} store
% q1 = g2/g1
% q2 = g3/g1
% q3 = g4/g1
% q4 = g5/g1
/q1 { g2 g1 div} store
/q2 { g3 g1 div} store
/q3 { g4 g1 div} store
/q4 { g5 g1 div} store
% r1 = h2/h1
% r2 = h3/h1
% r3 = h4/h1
% r4 = h5/h1
/r1 { h2 h1 div} store
/r2 { h3 h1 div} store
/r3 { h4 h1 div} store
/r4 { h5 h1 div} store
% s1 = i2/i1
% s2 = i3/i1
% s3 = i4/i1
% s4 = i5/i1
/s1 { i2 i1 div} store
/s2 { i3 i1 div} store
/s3 { i4 i1 div} store
/s4 { i5 i1 div} store
% b + c(p1) + d(p2) + e(p3) = (p4) % eqn II
% b + c(q1) + d(q2) + e(q3) = (q4) % eqn III
% b + c(r1) + d(r2) + e(r3) = (r4) % eqn IV
% b + c(s1) + d(s2) + e(s3) = (s4) % eqn V
% subtract eqn III from II and eqn IV from III and eqn V from IV
% c((p1)-(q1)) + d((p2)-(q2) + e ((p3)-(q3)) = (p4)-(q4)
% c((q1)-(r1)) + d((q2)-(r2) + e ((q3)-(r3)) = (q4)-(r4)
% c((r1)-(s1)) + d((r2)-(s2) + e ((r3)-(s3)) = (r4)-(s4)
% simplify variables
% t1 = (p1) - (q1)
% t2 = (p2) - (q2)
% t3 = (p3) - (q3)
% t4 = (p4) - (q4)
/t1 { p1 q1 sub} store
/t2 { p2 q2 sub} store
/t3 { p3 q3 sub} store
/t4 { p4 q4 sub} store
% u1 = (q1) - (r1)
% u2 = (q2) - (r2)
% u3 = (q3) - (r3)
% u4 = (q4) - (r4)
/u1 { q1 r1 sub} store
/u2 { q2 r2 sub} store
/u3 { q3 r3 sub} store
/u4 { q4 r4 sub} store
% v1 = (r1) - (s1)
% v2 = (r2) - (s2)
% v3 = (r3) - (s3)
% v4 = (r4) - (s4)
/v1 { r1 s1 sub} store
/v2 { r2 s2 sub} store
/v3 { r3 s3 sub} store
/v4 { r4 s4 sub} store
% c(t1) + d(t2) + e(t3) = (t4)
% c(u1) + d(u2) + e(u3) = (u4)
% c(v1) + d(v2) + e(v3) = (v4)
% unify c
% c + d(t2)/(t1) + e(t3)/(t1) = (t4/(t1)
% c + d(u2)/(u1) + e(u3)/(u1) = (u4/(u1)
% c + d(v2)/(v1) + e(v3)/(v1) = (v4/(v1)
% simplify variables
% w1 = (t2)/(t1)
% w2 = (t3)/(t1)
% w3 = (t4)/(t1)
/w1 {t2 t1 div} store
/w2 {t3 t1 div} store
/w3 {t4 t1 div} store
% x1 = (u2)/(u1)
% x2 = (u3)/(u1)
% x3 = (u4)/(u1)
/xx1 {u2 u1 div} store
/xx2 {u3 u1 div} store
/xx3 {u4 u1 div} store
% y1 = (v2)/(v1)
% y2 = (v3)/(v1)
% y3 = (v4)/(v1)
/yy1 {v2 v1 div} store
/yy2 {v3 v1 div} store
/yy3 {v4 v1 div} store
% c + d(w1) + e(w2) = (w3) % eqn VI
% c + d(x1) + e(x2) = (x3) % eqn VII
% c + d(y1) + e(y2) = (y3) % eqn VIII
% subtract VII from VI and VIII from VII
% d ((w1)-(xx1)) + e ((w2)-(xx2)) = ((w3)-(xx3))
% d ((xx1)-(yy1)) + e ((xx2)-(yy2)) = ((xx3)-(yy3))
% unitiize d
% d + e ((w2)-(xx2))/((w1)-(xx1)) = ((w3)-(xx3))/((w1)-(xx1))
% d + e ((xx2)-(yy2))/((xx1)-(yy1)) = ((xx3)-(yy3))/((xx1)-(yy1))
% simplify variables
% /aa1 = ((w2)-(xx2))/((w1)-(xx1))
% /aa2 = ((w3)-(xx3))/((w1)-(xx1))
% /bb1 = ((xx2)-(yy2))/((xx1)-(yy1))
% /bb2 = ((xx3)-(yy3))/((xx1)-(yy1))
/aa1 { w2 xx2 sub w1 xx1 sub div } store
/aa2 { w3 xx3 sub w1 xx1 sub div } store
/bb1 { xx2 yy2 sub xx1 yy1 sub div } store
/bb2 { xx3 yy3 sub xx1 yy1 sub div } store
% d + e(aa1) = (aa2) eqn IX
% d + e(bb1) = (bb2) eqn X
% subtract and solve for EEEE
% e((aa1)-(bb1)) = ((aa2)-(bb2))
% e = ((aa2)-(bb2))/((aa1)-(bb1))
/ee {aa2 bb2 sub aa1 bb1 sub div
} store % <----EEEE
% d = (aa2) - e(aa1) from IX
/dd {aa2 aa1 ee mul sub
} store % <----DDDD
% c = (w3) - d(w1) - e(w2) % eqn VI
/cc {w3 w1 dd mul sub w2 ee mul sub
} store % <---- CCCC
% b + c(p1) + d(p2) + e(p3) = (p4) % eqn II
% b = (p4) - c(p1) -d(p2) -e(p3)
/bb { p4 p1 cc mul sub p2 dd mul sub p3 ee mul sub
} store % <----BBBB
/aa {1
bb sub
cc sub dd sub ee sub } store % <----- AAAA
% ==========
/x1 0.2 store
/y1 0.1 store
/x2 0.4 store
/y2 0.3 store
/x3 0.5 store
/y3 0.35 store
/x4 0.7 store
/y4 0.65 store
100 100 10 setgrid
20 20 showgrid
0 0 mt dot
1 20 mul 1 20 mul mt dot
x1 20 mul y1 20 mul mt dot
x2 20 mul y2 20 mul mt dot
x3 20 mul y3 20 mul mt dot
x4 20 mul y4 20 mul mt dot
line1 0 0 mt 20 20 lineto stroke
line2
0 0 moveto
0 0.01 1 { /xx exch store
xx 20 mul
xx dup dup mul dup mul mul aa mul
xx dup mul dup mul bb mul add
xx dup dup mul mul cc mul add
xx dup mul dd mul add
xx ee mul add
20 mul lineto
} for
stroke
showpage
snap4 restore
%%%%%%%%%%% HEXIC SEVEN-POINT FIT %%%%%%%%%%%%%%
save /snap5 exch store
% fit a hexic to seven data points
% points are 0,0 xx1,yy1, xx2,yy2, xx3,yy3 xx4,yy4, ss5,yy5, and 1,1
% seek general solution.
% y = ax6 + bx5 + cx4 + dx3 + ex2 + fx + g no g since y=0 at x=0
% 1 = a + b + c + d + e + f at x = 1 EQUATION AA
% a = 1 - b - c - d - e - f % egn I
% write and scale to clean a
% a + b/(x1) + c/(x1)2 + d/(x1)3 + e/(x1)4 + f/(x1)5 = (y1)/(x1)6 % first point
% a + b/(x2) + c/(x2)2 + d/(x2)3 + e/(x2)4 + f/(x2)5 = (y2)/(x2)6 % second point
% a + b/(x3) + c/(x3)2 + d/(x3)3 + e/(x3)4 + f/(x3)5 = (y3)/(x3)6 % third point
% a + b/(x4) + c/(x4)2 + d/(x4)3 + e/(x4)4 + f/(x4)5 = (y4)/(x4)6 % fourth point
% a + b/(x5) + c/(x5)2 + d/(x5)3 + e/(x5)4 + f/(x5)5 = (y5)/(x5)6 % fifth point
% substitute and rearrange
% b (1-(x1))/(x1) + c (1-(x1)2)/(x1)2 + d (1-(x1)3)/(x1)3 + e(1-(x1)4)/(x1)4 +
% f (1-(x1)5/(x1)5 = (y1)/(x1)6 - 1
% b (1-(x2))/(x2) + c (1-(x2)2)/(x2)2 + d (1-(x2)3)/(x2)3 + e(1-(x2)4)/(x2)4 +
% f (1-(x2)5/(x2)5 = (y2)/(x2)6 - 1
% b (1-(x3))/(x3) + c (1-(x3)2)/(x3)2 + d (1-(x3)3)/(x)3 + e(1-(x3)4)/(x3)4 +
% f (1-(x3)5/(x3)5 = (y3)/(x3)6 - 1
% b (1-(x4))/(x4) + c (1-(x4)2)/(x4)2 + d (1-(x4)3)/(x4)3 + e(1-(x4)4)/(x4)4 +
% f (1-(x4)5/(x4)5 = (y4)/(x4)6 - 1
% b (1-(x5))/(x5) + c (1-(x5)2)/(x5)2 + d (1-(x5)3)/(x5)3 + e(1-(x5)4)/(x5)4 +
% f (1-(x5)5/(x5)5 = (y5)/(x5)6 - 1
% internalize p00 here
/p00 { 1 x1 sub x1 div } store
/p01 { 1 x1 2 exp sub x1 2 exp div p00 div} store
/p02 { 1 x1 3 exp sub x1 3 exp div p00 div} store
/p03 { 1 x1 4 exp sub x1 4 exp div p00 div} store
/p04 { 1 x1 5 exp sub x1 5 exp div p00 div} store
/p05 { y1 x1 6 exp div 1 sub p00 div} store
/p10 {1 x2 sub x2 div } store
/p11 { 1 x2 2 exp sub x2 2 exp div p10 div} store
/p12 { 1 x2 3 exp sub x2 3 exp div p10 div} store
/p13 { 1 x2 4 exp sub x2 4 exp div p10 div} store
/p14 { 1 x2 5 exp sub x2 5 exp div p10 div} store
/p15 { y2 x2 6 exp div 1 sub p10 div} store
/p20 {1 x3 sub x3 div } store
/p21 { 1 x3 2 exp sub x3 2 exp div p20 div} store
/p22 { 1 x3 3 exp sub x3 3 exp div p20 div} store
/p23 { 1 x3 4 exp sub x3 4 exp div p20 div} store
/p24 { 1 x3 5 exp sub x3 5 exp div p20 div} store
/p25 { y3 x3 6 exp div 1 sub p20 div} store
/p30 {1 x4 sub x4 div } store
/p31 { 1 x4 2 exp sub x4 2 exp div p30 div} store
/p32 { 1 x4 3 exp sub x4 3 exp div p30 div} store
/p33 { 1 x4 4 exp sub x4 4 exp div p30 div} store
/p34 { 1 x4 5 exp sub x4 5 exp div p30 div} store
/p35 { y4 x4 6 exp div 1 sub p30 div } store
/p40 {1 x5 sub x5 div } store
/p41 { 1 x5 2 exp sub x5 2 exp div p40 div} store
/p42 { 1 x5 3 exp sub x5 3 exp div p40 div} store
/p43 { 1 x5 4 exp sub x5 4 exp div p40 div} store
/p44 { 1 x5 5 exp sub x5 5 exp div p40 div} store
/p45 { y5 x5 6 exp div 1 sub p40 div } store
% unifying b
% b + c(p01) + d(p02) + e(p03)+ f (p04) = (p05) eqn II
% b + c(p11) + d(p12) + e(p13)+ f (p14) = (p15) eqn III
% b + c(p21) + d(p22) + e(p23)+ f (p24) = (p25) eqn IV
% b + c(p31) + d(p32) + e(p33)+ f (p34) = (p35) eqn V
% b + c(p41) + d(p42) + e(p43)+ f (p44) = (p45) eqn VI
% subtracting equation pairs
% c((p01)-(p11)) + d((p02)-(p12)) + e((p03)-(p13)) + f((p04)-(p14)) = ((p05)-(p15))
% etc... i hope
% and normalizing
/q00 {p01 p11 sub} store
/q01 {p02 p12 sub q00 div} store
/q02 {p03 p13 sub q00 div} store
/q03 {p04 p14 sub q00 div} store
/q04 {p05 p15 sub q00 div} store
/q10 {p11 p21 sub} store
/q11 {p12 p22 sub q10 div} store
/q12 {p13 p23 sub q10 div} store
/q13 {p14 p24 sub q10 div} store
/q14 {p15 p25 sub q10 div} store
/q20 {p21 p31 sub} store
/q21 {p22 p32 sub q20 div} store
/q22 {p23 p33 sub q20 div} store
/q23 {p24 p34 sub q20 div} store
/q24 {p25 p35 sub q20 div} store
/q30 {p31 p41 sub} store
/q31 {p32 p42 sub q30 div} store
/q32 {p33 p43 sub q30 div} store
/q33 {p34 p44 sub q30 div} store
/q34 {p35 p45 sub q30 div} store
% c + d(q01) + e(q02) + f(q03) = (q04) eqn VII
% c + d(q11) + e(q12) + f(q13) = (q14) eqn VIII
% c + d(q21) + e(q22) + f(q23) = (q24) eqn IX
% c + d(q31) + e(q32) + f(q33) = (q34) eqn X
% subtract equation pairs
% d((q01)-(q11)) + e((q02)-(q12)) + f((q03)-(q13)) = ((q04)-(q14))
% etc I hope
% and normalizing
/r00 {q01 q11 sub} store
/r01 {q02 q12 sub r00 div} store
/r02 {q03 q13 sub r00 div} store
/r03 {q04 q14 sub r00 div} store
/r10 {q11 q21 sub} store
/r11 {q12 q22 sub r10 div} store
/r12 {q13 q23 sub r10 div} store
/r13 {q14 q24 sub r10 div} store
/r20 {q21 q31 sub} store
/r21 {q22 q32 sub r20 div} store
/r22 {q23 q33 sub r20 div} store
/r23 {q24 q34 sub r20 div} store
% d + e(r01) + f(r02) = (r03) eqn XI
% d + e(r11) + f(r12) = (r13) eqn XII
% d + e(r21) + f(r22) = (r23) eqn XIII
% subtract equation pairs
% e((r01)-(r11)) + f ((r02)-r12)) = ((r03)-(r13))
% e((r11)-(r21)) + f ((r12)-r22)) = ((r13)-(r23))
/s00 {r01 r11 sub} store
/s01 {r02 r12 sub s00 div} store
/s02 {r03 r13 sub s00 div} store
/s10 {r11 r21 sub} store
/s11 {r12 r22 sub s10 div} store
/s12 {r13 r23 sub s10 div} store
% e + f(s01) = (s02)
% e + f(s11) = (s12)
/ff {s02 s12 sub s01 s11 sub div} store
/ee {s12 s11 ff mul sub} store
% d + e(r01) + f(r02) = (r03) eqn XI
/dd {r03 r01 ee mul sub r02 ff mul sub} store
% c + d(q01) + e(q02) + f(q03) = (q04) eqn VII
/cc { q04 dd q01 mul sub ee q02 mul sub ff q03 mul sub} store
% b + c(p01) + d(p02) + e(p03)+ f (p04) = (p05) eqn II
/bb {p05 cc p01 mul sub dd p02 mul sub ee p03 mul sub ff p04 mul sub } store
/aa {1 bb sub cc sub dd sub ee sub ff sub} store
% ========== demo plot ===========
/x1 0.2 store
/y1 0.1 store
/x2 0.4 store
/y2 0.3 store
/x3 0.5 store
/y3 0.35 store
/x4 0.7 store
/y4 0.65 store
/x5 0.8 store
/y5 0.85 store
100 100 10 setgrid
20 20 showgrid
0 0 mt dot
1 20 mul 1 20 mul mt dot
x1 20 mul y1 20 mul mt dot
x2 20 mul y2 20 mul mt dot
x3 20 mul y3 20 mul mt dot
x4 20 mul y4 20 mul mt dot
x5 20 mul y5 20 mul mt dot
line1 0 0 mt 20 20 lineto stroke
line2
0 0 moveto
0 0.01 1 { /xx exch store
xx 20 mul
xx 6 exp aa mul
xx 5 exp bb mul add
xx 4 exp cc mul add
xx 3 exp dd mul add
xx 2 exp ee mul add
xx ff mul add
20 mul lineto
} for
stroke
showpage
snap5 restore
%%%%%%%%%%%%%%%%
% Higher orders and consulting services available http://www.tinaja.com/info01.html
% EOF