-- LSCurveImpl.mesa
-- Last edited by Michael Plass November 29, 1982 9:52 am
-- Michael Plass and Maureen Stone Oct-81
DIRECTORY
Cubic,
Complex,
LinearSystem,
LSCurve,
PiecewiseCubic,
Real USING [RealException],
Seq,
Vector;
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]];