LSCurveImpl.mesa
Copyright © 1985 by Xerox Corporation. All rights reserved.
Michael Plass and Maureen Stone Oct-81
Last edited by Michael Plass November 29, 1982 9:52 am
Doug Wyatt, September 5, 1985 1:34:22 pm PDT
DIRECTORY
Complex USING [Abs, Add, Arg, Sub, VEC],
Cubic USING [Bezier, BezierToCoeffs, Coeffs],
LinearSystem USING [ColumnN, MatrixN, MatrixSeq, SolveN, VecSeq],
LSCurve USING [Handle, Patch, PatchSequence, PatchSequenceRec, SplineBasis, SplineBasisRec, SplineFunction, StateRec],
PiecewiseCubic USING [Combine, EnumerateCommonPieces, Eval, EvalAll, Handle, Piece, PieceProc, Zero],
Real USING [RealException],
Seq USING [BooleanSequence, BooleanSequenceRec, ComplexSequence, RealSequence, RealSequenceRec],
Vector USING [Div, Dot];
LSCurveImpl: CEDAR PROGRAM
IMPORTS Complex, Cubic, LinearSystem, PiecewiseCubic, Real, Vector
EXPORTS LSCurve
= BEGIN OPEN Seq;
RealProc: TYPE = PROCEDURE[i:NAT] RETURNS [REAL];
PointNumber: TYPE = NAT;
Patch: TYPE = LSCurve.Patch;
PatchSequence: TYPE = LSCurve.PatchSequence;
PatchSequenceRec: TYPE = LSCurve.PatchSequenceRec;
SplineBasis: TYPE = LSCurve.SplineBasis;
SplineBasisRec: TYPE = LSCurve.SplineBasisRec;
SplineFunction: TYPE = LSCurve.SplineFunction;
Handle: TYPE = LSCurve.Handle;
StateRec: TYPE = LSCurve.StateRec;
Create:
PUBLIC
PROCEDURE [sa: Seq.ComplexSequence]
RETURNS [h: Handle] =
BEGIN
h ← SamplesToHandle[sa];
h.weight ← NEW[RealSequenceRec[h.z.length]];
FOR i:NAT IN [0..h.z.length) DO h.weight[i] ← 1.0 ENDLOOP;
h.t ← NEW[RealSequenceRec[h.z.length]];
h.oldt ← NEW[RealSequenceRec[h.z.length]];
h.oldoldt ← NEW[RealSequenceRec[h.z.length]];
END;
XYat:
PUBLIC
PROCEDURE [h: Handle, t:
REAL]
RETURNS [f: Complex.
VEC] =
BEGIN
patchNumber: NAT ← 0;
ValidateCache[h];
f.x ← PiecewiseCubic.Eval[h.patchCache.x,t];
f.y ← PiecewiseCubic.Eval[h.patchCache.y,t];
END;
SamplesToHandle:
PROCEDURE [s: Seq.ComplexSequence]
RETURNS [h: Handle] =
{h ← NEW[StateRec]; BEGIN OPEN h^;
z ← s;
l ← 0;
n ← s.length-1;
END};
CountTrues:
PROCEDURE [b: BooleanSequence]
RETURNS [count:
NAT] =
BEGIN
count ← 0;
FOR i:NAT IN [0..b.length) DO IF b[i] THEN count ← count + 1 ENDLOOP;
END;
EvalBasis:
PUBLIC
PROCEDURE [h: Handle, i:
NAT, t:
REAL]
RETURNS [z: Complex.VEC] = INLINE
{OPEN h.splineBasis[i]; RETURN [[PiecewiseCubic.Eval[x,t],PiecewiseCubic.Eval[y,t]]]};
SolveCurve:
PUBLIC
PROCEDURE [h: Handle] =
BEGIN OPEN h,LinearSystem;
nFree: NAT ← CountTrues[free];
b: ColumnN ← NEW[VecSeq[nFree]];
c: ColumnN ← NIL;
A: MatrixN ← NEW[MatrixSeq[nFree]];
FOR i:
NAT
IN [0..nFree)
DO
A[i] ← NEW[VecSeq[nFree]];
FOR j: NAT IN [0..nFree) DO A[i][j] ← 0 ENDLOOP;
b[i] ← 0;
ENDLOOP;
FOR k:
NAT
IN [l..n]
DO
tk:REAL = t[k];
row: NAT ← 0;
FOR i:
NAT
IN [0..a.length)
DO
IF free[i]
THEN
BEGIN
gi:Complex.VEC ← EvalBasis[h,i,tk];
IF gi#[0,0]
THEN
BEGIN
column: NAT ← 0;
FOR j:
NAT
IN [0..a.length)
DO
gj: Complex.VEC ← EvalBasis[h,j,tk];
IF free[j]
THEN
A[row][column] ← A[row][column] + weight[k]*Vector.Dot[gi,gj]
ELSE
b[row] ← b[row] - weight[k]*a[j]*Vector.Dot[gi,gj];
IF free[j] THEN column ← column + 1;
ENDLOOP;
END;
b[row] ← b[row] + weight[k]*Vector.Dot[z[k],gi];
row ← row + 1;
END ENDLOOP;
ENDLOOP;
c ← SolveN[A,b,nFree!
Real.RealException => CHECKED {CONTINUE}];
IF c#
NIL
THEN
BEGIN
row: NAT ← 0;
FOR i:
NAT
IN [0..a.length)
DO
IF free[i] THEN a[i] ← c[row];
IF free[i] THEN row ← row + 1;
ENDLOOP;
patchCacheValid ← FALSE;
END;
END;
ArcLengthInitialTValues:
PUBLIC
PROCEDURE [h: Handle] =
BEGIN OPEN h;
arcLen: REAL ← 0;
t[0] ← 0;
FOR i:
NAT
IN [1..z.length)
DO
arcLen ← arcLen + Complex.Abs[Complex.Sub[z[i],z[i-1]]];
t[i] ← arcLen;
ENDLOOP;
FOR i:
NAT
IN [0..z.length)
DO
t[i] ← t[i]/arcLen
ENDLOOP;
END;
UnitInitialTValues:
PUBLIC
PROCEDURE [h: Handle] =
BEGIN OPEN h;
IF z.length < 2 THEN RETURN;
FOR i:
NAT
IN [0..z.length)
DO
t[i] ← i*1.0/(z.length-1);
ENDLOOP;
END;
AngleInitialTValues:
PUBLIC
PROCEDURE [h: Handle] =
BEGIN OPEN h;
twopi: REAL = 2*3.1415926;
m: Complex.VEC ← [0,0];
FOR i:
NAT
IN [0..z.length)
DO
m ← Complex.Add[m,z[i]];
ENDLOOP;
m ← Vector.Div[m,z.length];
FOR i:
NAT
IN [0..z.length)
DO
t[i] ← Complex.Arg[Complex.Sub[z[i],m]]/twopi;
WHILE t[i] < 0 DO t[i] ← t[i] + 1 ENDLOOP;
WHILE t[i] > 1 DO t[i] ← t[i] - 1 ENDLOOP;
ENDLOOP;
END;
AdjustTValues:
PUBLIC
PROCEDURE [h: Handle]
RETURNS [maxChange:
REAL] =
BEGIN OPEN h;
delta: REAL;
{newt: RealSequence ← oldoldt; oldoldt ← oldt; oldt ← t; t ← newt};
ValidateCache[h];
maxChange ← 0;
FOR i:
NAT
IN [l..n]
DO
delta ← Adjustment[z[i],patchCache,oldt[i]];
IF ABS[delta] > maxChange THEN maxChange ← ABS[delta];
t[i] ← oldt[i]+delta;
IF closedCurve
THEN
BEGIN
WHILE t[i] < 0 DO t[i] ← t[i] + 1 ENDLOOP;
WHILE t[i] > 1 DO t[i] ← t[i] - 1 ENDLOOP;
END
ENDLOOP;
IF initialEndFree
THEN
BEGIN
t0: REAL ← t[l];
FOR i:
NAT
IN [l..n]
DO
t[i] ← t[i] - t0;
ENDLOOP;
END;
IF finalEndFree
THEN
BEGIN
interval: REAL ← t[n];
FOR i:
NAT
IN [l..n]
DO
t[i] ← t[i]/interval;
ENDLOOP;
END;
END;
Adjustment:
PROCEDURE [z: Complex.
VEC, f: SplineFunction, t:
REAL]
RETURNS [delta:REAL] =
BEGIN
epsilon: REAL ← 1.0/8388608;
x,dx,ddx,dddx: REAL;
y,dy,ddy,dddy: REAL;
lo,hi: REAL;
derivSqrDist: REAL;
derivDerivSqrDist: REAL;
[x,dx,ddx,dddx,lo,hi] ← PiecewiseCubic.EvalAll[f.x, t];
[y,dy,ddy,dddy,lo,hi] ← PiecewiseCubic.EvalAll[f.y, t];
derivSqrDist ← (x-z.x)*dx + (y-z.y)*dy;
derivDerivSqrDist ← dx*dx + dy*dy + (x-z.x)*ddx + (y-z.y)*ddy;
delta ← IF derivDerivSqrDist<=0 THEN 0 ELSE -derivSqrDist/derivDerivSqrDist;
IF t+delta>hi THEN delta ← hi+epsilon-t;
IF t+delta<lo THEN delta ← lo-epsilon-t;
END;
Sqr: PROC [x:REAL] RETURNS [REAL] = INLINE{RETURN[x*x]};
AccelTValues:
PUBLIC
PROCEDURE [h: Handle]
RETURNS [r, maxChange:
REAL] =
BEGIN OPEN h;
sqrOldDelta, deltaDotOldDelta: REAL ← 0;
r ← 0;
maxChange ← AdjustTValues[h];
SolveCurve[h]; maxChange ← AdjustTValues[h];
FOR i:
NAT
IN [l..n]
DO
sqrOldDelta ← sqrOldDelta + Sqr[oldt[i]-oldoldt[i]];
deltaDotOldDelta ← deltaDotOldDelta + (t[i]-oldt[i])*(oldt[i]-oldoldt[i]);
ENDLOOP;
IF sqrOldDelta=0 THEN RETURN;
r ← deltaDotOldDelta/sqrOldDelta;
IF ABS[r]>=1 THEN RETURN;
maxChange ← 0;
FOR i:
NAT
IN [l..n]
DO
delta: REAL ← (t[i] - oldt[i])/(1-r);
t[i] ← oldt[i] + delta;
IF ABS[delta] > maxChange THEN maxChange ← ABS[delta];
ENDLOOP;
END;
PatchesOf:
PUBLIC
PROCEDURE [h: Handle]
RETURNS [x, y: PatchSequence, knots: RealSequence] =
BEGIN
i,nPatches: NAT ← 0;
CountPatches: PiecewiseCubic.PieceProc = {nPatches ← nPatches+1};
StorePatch: PiecewiseCubic.PieceProc =
BEGIN
interval: REAL ← p.domainEnd-p.domainStart;
b: Cubic.Bezier ←
[b0: [p.initValue, q.initValue],
b1: [p.initValue + p.initSlope*interval/3, q.initValue + q.initSlope*interval/3],
b2: [p.finalValue - p.finalSlope*interval/3, q.finalValue - q.finalSlope*interval/3],
b3: [p.finalValue, q.finalValue]];
c: Cubic.Coeffs ← Cubic.BezierToCoeffs[b];
x[i] ← [c0:c.c0.x, c1:c.c1.x, c2:c.c2.x, c3:c.c3.x];
y[i] ← [c0:c.c0.y, c1:c.c1.y, c2:c.c2.y, c3:c.c3.y];
knots[i] ← p.domainStart;
knots[i+1] ← p.domainEnd;
i ← i + 1;
END;
ValidateCache[h];
PiecewiseCubic.EnumerateCommonPieces[h.patchCache.x,h.patchCache.y,CountPatches];
x ← NEW[PatchSequenceRec[nPatches]];
y ← NEW[PatchSequenceRec[nPatches]];
knots ← NEW[RealSequenceRec[nPatches+1]];
PiecewiseCubic.EnumerateCommonPieces[h.patchCache.x,h.patchCache.y,StorePatch];
END;
ValidateCache:
PROCEDURE [h: Handle] =
BEGIN OPEN h;
IF patchCacheValid THEN RETURN;
patchCache ← [PiecewiseCubic.Zero[],PiecewiseCubic.Zero[]];
FOR i:
NAT
DECREASING
IN [0..a.length)
DO
patchCache ← [PiecewiseCubic.Combine[1,patchCache.x,a[i],splineBasis[i].x],
PiecewiseCubic.Combine[1,patchCache.y,a[i],splineBasis[i].y]]
ENDLOOP;
patchCacheValid←TRUE;
END;
PointSlopeBasis:
PUBLIC
PROCEDURE [h: Handle, b: Cubic.Bezier] =
BEGIN OPEN h,PiecewiseCubic;
vx,vy: REAL;
closedCurve ← FALSE;
initialEndFree ← FALSE;
finalEndFree ← FALSE;
patchCacheValid ← FALSE;
splineBasis ← NEW[SplineBasisRec[8]];
a ← NEW[RealSequenceRec[8]];
free ← NEW[BooleanSequenceRec[8]];
(a0,a1) and (a6,a7) are the positions of the endpoints
a2 is the velocity component parallel to b1-b0
a3 is the velocity component perpendicular to b1-b0
a4 is the velocity component parallel to b2-b3
a5 is the velocity component perpendicular to b2-b3
splineBasis[0] ← [Piece[0,1,1,0,0,0],Zero[]]; a[0] ← b.b0.x;
splineBasis[1] ← [Zero[],Piece[0,1,1,0,0,0]]; a[1] ← b.b0.y;
[[vx,vy]] ← Complex.Sub[b.b1,b.b0];
splineBasis[2] ← [Combine[3*vx,Piece[0,1,0,1,0,0],0,Zero[]],
Combine[3*vy,Piece[0,1,0,1,0,0],0,Zero[]]]; a[2] ← 1;
splineBasis[3] ← [Combine[-3*vy,Piece[0,1,0,1,0,0],0,Zero[]],
Combine[3*vx,Piece[0,1,0,1,0,0],0,Zero[]]]; a[3] ← 0;
[[vx,vy]] ← Complex.Sub[b.b3,b.b2];
splineBasis[4] ← [Combine[3*vx,Piece[0,1,0,0,1,0],0,Zero[]],
Combine[3*vy,Piece[0,1,0,0,1,0],0,Zero[]]]; a[4] ← 1;
splineBasis[5] ← [Combine[-3*vy,Piece[0,1,0,0,1,0],0,Zero[]],
Combine[3*vx,Piece[0,1,0,0,1,0],0,Zero[]]]; a[5] ← 0;
splineBasis[6] ← [Piece[0,1,0,0,0,1],Zero[]]; a[6] ← b.b3.x;
splineBasis[7] ← [Zero[],Piece[0,1,0,0,0,1]]; a[7] ← b.b3.y;
FOR i:
NAT
IN [0..8)
DO
free[i] ← FALSE;
ENDLOOP;
END;
SmoothBasis:
PUBLIC
PROCEDURE [h: Handle, nKnots:
NAT] =
BEGIN OPEN h;
closedCurve ← TRUE;
initialEndFree ← FALSE;
finalEndFree ← FALSE;
patchCacheValid ← FALSE;
splineBasis ← NEW[SplineBasisRec[4*nKnots]];
a ← NEW[RealSequenceRec[4*nKnots]];
free ← NEW[BooleanSequenceRec[4*nKnots]];
FOR i:
NAT
IN [0..nKnots)
DO
even: PiecewiseCubic.Handle ← EvenBasis[i*1.0/nKnots,1.0/nKnots];
odd: PiecewiseCubic.Handle ← OddBasis[i*1.0/nKnots,1.0/nKnots];
splineBasis[4*i] ← [even,PiecewiseCubic.Zero[]];
splineBasis[4*i+1] ← [PiecewiseCubic.Zero[],even];
splineBasis[4*i+2] ← [odd,PiecewiseCubic.Zero[]];
splineBasis[4*i+3] ← [PiecewiseCubic.Zero[],odd];
ENDLOOP;
FOR i:
NAT
IN [0..4*nKnots)
DO
a[i] ← 0;
free[i] ← TRUE;
ENDLOOP;
END;
EvenBasis:
PROCEDURE [center, radius:
REAL]
RETURNS [PiecewiseCubic.Handle] =
BEGIN OPEN PiecewiseCubic;
l: REAL ← center-radius;
IF l<0 THEN l←l+1;
RETURN [Combine[1,Piece[l,l+radius,0,0,0,1],
1,Piece[center,center+radius,1,0,0,0]]]
END;
OddBasis:
PROCEDURE [center, radius:
REAL]
RETURNS [PiecewiseCubic.Handle] =
BEGIN OPEN PiecewiseCubic;
l: REAL ← center-radius;
IF l<0 THEN l←l+1;
RETURN [Combine[1,Piece[l,l+radius,0,0,1,0],
1,Piece[center,center+radius,0,1,0,0]]]
END;
BSplineBasis:
PUBLIC
PROCEDURE [h: Handle, nKnots:
NAT] =
BEGIN OPEN h;
closedCurve ← TRUE;
initialEndFree ← FALSE;
finalEndFree ← FALSE;
patchCacheValid ← FALSE;
splineBasis ← NEW[SplineBasisRec[2*nKnots]];
a ← NEW[RealSequenceRec[2*nKnots]];
free ← NEW[BooleanSequenceRec[2*nKnots]];
FOR i:
NAT
IN [0..nKnots)
DO
elem: PiecewiseCubic.Handle ← BSplineBasisElement[i,nKnots];
splineBasis[2*i] ← [elem,PiecewiseCubic.Zero[]];
splineBasis[2*i+1] ← [PiecewiseCubic.Zero[],elem];
ENDLOOP;
FOR i:
NAT
IN [0..2*nKnots)
DO
a[i] ← 0;
free[i] ← TRUE;
ENDLOOP;
END;
BSplineBasisElement:
PROCEDURE [i, n:
NAT]
RETURNS [f: PiecewiseCubic.Handle] =
BEGIN
k: REAL ← i-2;
IF k<0 THEN k←k+n;
f ← PiecewiseCubic.Piece[k/n,(k+1)/n,0,0,0.75*n,0.25];
k←k+1; IF k>n-1 THEN k←k-n;
f ← PiecewiseCubic.Combine[1,f,1,PiecewiseCubic.Piece[k/n,(k+1)/n,0.25,0.75*n,0,1]];
k←k+1; IF k>n-1 THEN k←k-n;
f ← PiecewiseCubic.Combine[1,f,1,PiecewiseCubic.Piece[k/n,(k+1)/n,1,0,-0.75*n,0.25]];
k←k+1; IF k>n-1 THEN k←k-n;
f ← PiecewiseCubic.Combine[1,f,1,PiecewiseCubic.Piece[k/n,(k+1)/n,0.25,-0.75*n,0,0]];
END;
CubicBasis:
PUBLIC
PROCEDURE [h: Handle] =
BEGIN OPEN h;
zero: PiecewiseCubic.Handle ← PiecewiseCubic.Zero[];
one: PiecewiseCubic.Handle ← PiecewiseCubic.Piece[-10,10,1,0,0,1];
iden: PiecewiseCubic.Handle ← PiecewiseCubic.Piece[-10,10,-10,1,1,10];
sqr: PiecewiseCubic.Handle ← PiecewiseCubic.Piece[-10,10,100,-20,20,100];
cube: PiecewiseCubic.Handle ← PiecewiseCubic.Piece[-10,10,-1000,300,300,1000];
closedCurve ← FALSE;
initialEndFree ← TRUE;
finalEndFree ← TRUE;
patchCacheValid ← FALSE;
splineBasis ← NEW[SplineBasisRec[8]];
a ← NEW[RealSequenceRec[8]];
free ← NEW[BooleanSequenceRec[8]];
splineBasis[0] ← [one,zero];
splineBasis[1] ← [zero,one];
splineBasis[2] ← [iden,zero];
splineBasis[3] ← [zero,iden];
splineBasis[4] ← [sqr,zero];
splineBasis[5] ← [zero,sqr];
splineBasis[6] ← [cube,zero];
splineBasis[7] ← [zero,cube];
FOR i:
NAT
IN [0..8)
DO
a[i] ← 0;
free[i] ← TRUE;
ENDLOOP;
END;
CubicPieceBasis:
PUBLIC
PROCEDURE [h: Handle] =
BEGIN OPEN h;
zero: PiecewiseCubic.Handle ← PiecewiseCubic.Zero[];
one: PiecewiseCubic.Handle ← PiecewiseCubic.Piece[0,1,1,0,0,1];
iden: PiecewiseCubic.Handle ← PiecewiseCubic.Piece[0,1,0,1,1,1];
sqr: PiecewiseCubic.Handle ← PiecewiseCubic.Piece[0,1,0,0,2,1];
cube: PiecewiseCubic.Handle ← PiecewiseCubic.Piece[0,1,0,0,3,1];
closedCurve ← FALSE;
initialEndFree ← TRUE;
finalEndFree ← TRUE;
patchCacheValid ← FALSE;
splineBasis ← NEW[SplineBasisRec[8]];
a ← NEW[RealSequenceRec[8]];
free ← NEW[BooleanSequenceRec[8]];
splineBasis[0] ← [one,zero];
splineBasis[1] ← [zero,one];
splineBasis[2] ← [iden,zero];
splineBasis[3] ← [zero,iden];
splineBasis[4] ← [sqr,zero];
splineBasis[5] ← [zero,sqr];
splineBasis[6] ← [cube,zero];
splineBasis[7] ← [zero,cube];
FOR i:
NAT
IN [0..8)
DO
a[i] ← 0;
free[i] ← TRUE;
ENDLOOP;
END;
ParabolaBasis:
PUBLIC
PROCEDURE [h: Handle] =
BEGIN OPEN h;
zero: PiecewiseCubic.Handle ← PiecewiseCubic.Zero[];
one: PiecewiseCubic.Handle ← PiecewiseCubic.Piece[-10,10,1,0,0,1];
iden: PiecewiseCubic.Handle ← PiecewiseCubic.Piece[-10,10,-10,1,1,10];
sqr: PiecewiseCubic.Handle ← PiecewiseCubic.Piece[-10,10,100,-20,20,100];
closedCurve ← FALSE;
initialEndFree ← TRUE;
finalEndFree ← TRUE;
patchCacheValid ← FALSE;
splineBasis ← NEW[SplineBasisRec[6]];
a ← NEW[RealSequenceRec[6]];
free ← NEW[BooleanSequenceRec[6]];
splineBasis[0] ← [one,zero];
splineBasis[1] ← [zero,one];
splineBasis[2] ← [iden,zero];
splineBasis[3] ← [zero,iden];
splineBasis[4] ← [sqr,zero];
splineBasis[5] ← [zero,sqr];
FOR i:
NAT
IN [0..6)
DO
a[i] ← 0;
free[i] ← TRUE;
ENDLOOP;
END;
CircleBasis:
PUBLIC
PROCEDURE [h: Handle] =
BEGIN OPEN h;
twopi: REAL ← 2*3.1415926;
zero: PiecewiseCubic.Handle ← PiecewiseCubic.Zero[];
one: PiecewiseCubic.Handle ← PiecewiseCubic.Piece[0,1,1,0,0,1];
cos: PiecewiseCubic.Handle ← CosBasis[];
sin: PiecewiseCubic.Handle ← SinBasis[];
closedCurve ← TRUE;
initialEndFree ← FALSE;
finalEndFree ← FALSE;
patchCacheValid ← FALSE;
splineBasis ← NEW[SplineBasisRec[3]];
a ← NEW[RealSequenceRec[3]];
free ← NEW[BooleanSequenceRec[3]];
splineBasis[0] ← [one,zero];
splineBasis[1] ← [zero,one];
splineBasis[2] ← [cos,sin];
FOR i:
NAT
IN [0..3)
DO
a[i] ← 0;
free[i] ← TRUE;
ENDLOOP;
END;
EllipseBasis:
PUBLIC
PROCEDURE [h: Handle] =
BEGIN OPEN h;
twopi: REAL ← 2*3.1415926;
zero: PiecewiseCubic.Handle ← PiecewiseCubic.Zero[];
one: PiecewiseCubic.Handle ← PiecewiseCubic.Piece[0,1,1,0,0,1];
cos: PiecewiseCubic.Handle ← CosBasis[];
sin: PiecewiseCubic.Handle ← SinBasis[];
closedCurve ← TRUE;
initialEndFree ← FALSE;
finalEndFree ← FALSE;
patchCacheValid ← FALSE;
splineBasis ← NEW[SplineBasisRec[5]];
a ← NEW[RealSequenceRec[5]];
free ← NEW[BooleanSequenceRec[5]];
splineBasis[0] ← [one,zero];
splineBasis[1] ← [zero,one];
splineBasis[2] ← [cos,sin];
splineBasis[3] ← [sin,zero];
splineBasis[4] ← [zero,cos];
FOR i:
NAT
IN [0..5)
DO
a[i] ← 0;
free[i] ← TRUE;
ENDLOOP;
END;
CosBasis:
PROC
RETURNS[PiecewiseCubic.Handle] =
{RETURN[LIST[[domainStart: 0, domainEnd: 0.25, initValue: 99.98286, initSlope: 1.302177e-2, finalSlope: -662.7453, finalValue: -1.397745e-4], [domainStart: 0.25, domainEnd: 0.5, initValue: -1.397745e-4, initSlope: -662.7453, finalSlope: 2.099349e-3, finalValue: -99.9828], [domainStart: 0.5, domainEnd: 0.75, initValue: -99.9828, initSlope: 2.099349e-3, finalSlope: 662.7528, finalValue: 6.351456e-6], [domainStart: 0.75, domainEnd: 1, initValue: 6.351456e-6, initSlope: 662.7528, finalSlope: 1.302177e-2, finalValue: 99.98286]]]};
SinBasis:
PROC
RETURNS[PiecewiseCubic.Handle] =
{RETURN[LIST[[domainStart: 0, domainEnd: 0.25, initValue: 1.728265e-4, initSlope: 662.7498, finalSlope: -4.030657e-4, finalValue: 99.98302], [domainStart: 0.25, domainEnd: 0.5, initValue: 99.98302, initSlope: -4.030657e-4, finalSlope: -662.7458, finalValue: 1.005444e-4], [domainStart: 0.5, domainEnd: 0.75, initValue: 1.005444e-4, initSlope: -662.7458, finalSlope: 1.04512e-3, finalValue: -99.98254], [domainStart: 0.75, domainEnd: 1, initValue: -99.98254, initSlope: 1.04512e-3, finalSlope: 662.7498, finalValue: 1.728265e-4]]]};
END.
Quick approximations
cos ← PiecewiseCubic.Combine[1,PiecewiseCubic.Combine[1,PiecewiseCubic.Combine[1,
PiecewiseCubic.Piece[0,0.25,1,0,-twopi,0],1,
PiecewiseCubic.Piece[0.25,0.5,0,-twopi,0,-1]],1,
PiecewiseCubic.Piece[0.5,0.75,-1,0,twopi,0]],1,
PiecewiseCubic.Piece[0.75,1,0,twopi,0,1]];
sin ← PiecewiseCubic.Combine[1,PiecewiseCubic.Combine[1,PiecewiseCubic.Combine[1,
PiecewiseCubic.Piece[0,0.25,0,twopi,0,1],1,
PiecewiseCubic.Piece[0.25,0.5,1,0,-twopi,0]],1,
PiecewiseCubic.Piece[0.5,0.75,0,-twopi,0,-1]],1,
PiecewiseCubic.Piece[0.75,1,-1,0,twopi,0]];