% given points (x0,y0) (x1,y1) (x2,y2) (x3,y3) on stack, % draw the shortest Bezier cubic spline through them... /Bez4ptsviaLagrange { /y3 exch store /x3 exch store /y2 exch store /x2 exch store /y1 exch store /y1 exch store /y0 exch store /x0 exch store % often the currentpoint pair % Finding the Lagrange Polynomial... /y0d0 y0 x0 x1 sub x0 x2 sub mul x0 x3 sub mul div def /y1d1 y1 x1 x0 sub x1 x2 sub mul x1 x3 sub mul div def /y2d2 y2 x2 x0 sub x2 x1 sub mul x2 x3 sub mul div def /y3d3 y3 x3 x0 sub x3 x1 sub mul x3 x2 sub mul div def /a3 y0d0 y1d1 add y2d2 add y3d3 add def /a2 y0d0 x1 x2 add x3 add mul neg y1d1 x0 x2 add x3 add mul sub y2d2 x0 x1 add x3 add mul sub y3d3 x0 x1 add x2 add mul sub def /a1 y0d0 x1 x2 mul x1 x3 mul add x2 x3 mul add mul y1d1 x0 x2 mul x0 x3 mul add x2 x3 mul add mul add y2d2 x0 x1 mul x0 x3 mul add x1 x3 mul add mul add y3d3 x0 x1 mul x0 x2 mul add x1 x2 mul add mul add def /a0 y0d0 x1 x2 mul x3 mul mul neg y1d1 x0 x2 mul x3 mul mul sub y2d2 x0 x1 mul x3 mul mul sub y3d3 x0 x1 mul x2 mul mul sub def % % L(x) = a3 * x^3 + a2 * x^2 + a1 * x + a0 % % L is called with x on stack, returns x and y = L(x) on stack % /L { dup dup dup mul mul a3 mul 1 index dup mul a2 mul add 1 index a1 mul add a0 add } store % Converting the Lagrange Polynomial to a cubic spline curveto... /L' {/xx exch store % find slope of x on stack xx dup mul a3 mul 3 mul xx a2 mul 2 mul add a1 add } store /cx1 {x0 x3 x0 sub 1 3 div mul add } store /cx2 {x3 x3 x0 sub 1 3 div mul sub } store /m1 {x0 L' } store /b1 {y0 x0 ml mul sub } store /cy1 {mu cx1 mul b1 add} store /m2 {x3 L' mul } store /b2 {y3 m2 x3 mul sub } store /cy2 {m2 cx2 mul b2 add } store % Creating or appending the curve x0 y0 moveto % or use currentpoint cx1 xy1 cx2 cy2 x3 y3 curveto } store