%!PS

% Bezier through four points via Lagrange demo
% =============================================
% by  Don Lancaster

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Copyright c 2005, 2010 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 2.1 

%  Creates the shortest, simplest and most attractive Bezier Cubic Spline through
%  four data points. But does not preserve initial slope continuity.

%  Certain XY mnipulations may be 5X slower and increase with the resolution square.
%  A minimum of 250 megabytes of system RAM is recommended with no other programs active.

%  NOTE THAT ALL PS FILENAME STRINGS !!!DEMAND!!! DOUBLE REVERSE SLASHES.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% (C:\\Documents and Settings\\don\\Desktop\\gonzo\\gonzo.ps) run  % use internal 

(C:\\Users\\don\\Desktop\\Ghost\\gonzo.ps) run  % modify for current machine


% ========

/guru { gonzo begin
ps.util.1 begin printerror nuisance begin} def   

guru                                           % optionally activate gonzo utilities

% =========


100 200 20  setgrid     % draw a demo grid
20 20 showgrid

/dotsrus [ 3 5 7 14 11 2 18 5 ] store     % Bezier control points x0 y0 x1 ....



/Bez4ptsviaLagrange {                     % the routine

aload pop   

/y3 exch store /x3 exch store
/y2 exch store /x2 exch store
/y1 exch store /x1 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 m1 mul sub } store
/cy1 {m1 cx1 mul b1 add} store

/m2 {x3 L' } 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 

(x0 is now ) print x0 ==
(y0 is now ) print x1 ==

(cx1 is now ) print cx1 ==
(cy1 is now ) print cy1 ==

(cx2 is now ) print cx2 ==
(cy2 is now ) print cy2 ==



cx1 cy1 cx2 cy2 x3 y3 curveto

} store


dotsrus Bez4ptsviaLagrange 

.6 .6 1 setrgbcolor line3 stroke

black dotsrus aload pop  4 {mt dot} repeat      % Plot the dots

1 0 0 setrgbcolor cx1 cy1 mt dot  cx2 cy2 mt dot 




 showpage

% EOF