--DCubics.mesa
--m.stone September 26, 1980 5:37 PM
--constant step (64) forward difference, lines in between

DIRECTORY
PointDefs: FROM "PointDefs",
Real: FROM "Real",
SplineDefs: FROM "SplineDefs",
InlineDefs: FROM "InlineDefs";

DCubics: PROGRAM IMPORTS Real,InlineDefs EXPORTS SplineDefs =
BEGIN OPEN SplineDefs;

DisplayCubic: PUBLIC PROCEDURE
[coeffs: POINTER TO Coeffs,
MoveTo: PROCEDURE[PointDefs.ScrPt],
DrawTo: PROCEDURE[PointDefs.ScrPt]] =
BEGIN
i,j,nsteps: INTEGER;
pt: PointDefs.ScrPt;
ONE: REAL ← 1;
fd: ARRAY [0..3] OF ARRAY [1..NDIM] OF LONG INTEGER;

nsteps ← 64;
--make the forward differencing matrix
BEGIN
eps,eps2,eps3: REAL;
scaler: REAL ← 1000000B;--2↑18
eps ← ONE/nsteps;eps2 ← eps*eps;eps3 ← eps2*eps;
FOR i IN [1..NDIM] DO
fd[0][i] ← Real.RoundLI[scaler*(6*eps3*coeffs.t3[i])];
fd[1][i] ← Real.RoundLI[scaler*(6*eps3*coeffs.t3[i]+2*eps2*coeffs.t2[i])];
fd[2][i] ← Real.RoundLI[scaler*(eps3*coeffs.t3[i]+eps2*coeffs.t2[i]+eps*coeffs.t1[i])];
fd[3][i] ← Real.RoundLI[scaler*(coeffs.t0[i])];
ENDLOOP;
END;
pt[PointDefs.X] ← DoubleToInt[fd[3][1]];
pt[PointDefs.Y] ← DoubleToInt[fd[3][2]];
MoveTo[pt];
FOR i IN [1..nsteps] DO
--Turn the crank
FOR j IN [1..NDIM] DO
fd[3][j] ← fd[3][j] + fd[2][j];
fd[2][j] ← fd[2][j] + fd[1][j];
fd[1][j] ← fd[1][j] + fd[0][j];
ENDLOOP;

pt[PointDefs.X] ← DoubleToInt[fd[3][1]];
pt[PointDefs.Y] ← DoubleToInt[fd[3][2]];
DrawTo[pt];
ENDLOOP;
END;


--convert a REAL to an INTEGER

--round a Double to an INTEGER
DoubleToInt: PROCEDURE [double: LONG INTEGER] RETURNS[i: INTEGER] =
BEGIN
scaler: LONG INTEGER ← 1000000B;--2↑18
i ← InlineDefs.LowHalf[double/scaler]+(IF double MOD scaler >scaler/2 THEN 1 ELSE 0);
END;
END.