PiecewiseCubicImpl.mesa
Copyright Ó 1985, 1992 by Xerox Corporation. All rights reserved.
Michael Plass 15-Dec-81 8:38:56
Doug Wyatt, September 5, 1985 1:36:00 pm PDT
DIRECTORY
PiecewiseCubic USING [Handle, PieceProc, PieceRec];
PiecewiseCubicImpl: CEDAR PROGRAM
EXPORTS PiecewiseCubic
= BEGIN
Handle: TYPE = PiecewiseCubic.Handle;
PieceRec: TYPE = PiecewiseCubic.PieceRec;
Zero: PUBLIC PROCEDURE RETURNS [Handle] = {RETURN [NIL]};
Piece: PUBLIC PROCEDURE
[domainStart, domainEnd: REAL,
initValue, initSlope, finalSlope, finalValue: REAL]
RETURNS [Handle] =
BEGIN
IF domainStart<domainEnd AND
(initValue#0 OR initSlope#0 OR finalSlope#0 OR finalValue#0) THEN
BEGIN
piece: PieceRec ¬ [domainStart:domainStart, domainEnd:domainEnd, initValue:initValue, initSlope:initSlope,finalSlope:finalSlope, finalValue:finalValue];
RETURN[LIST[piece]]
END
ELSE RETURN[NIL]
END;
ScalePiece: PROCEDURE [p: PieceRec, a: REAL] RETURNS [PieceRec] =
BEGIN OPEN p;
RETURN [[domainStart:domainStart, domainEnd:domainEnd,
initValue:initValue*a,
initSlope:initSlope*a,
finalSlope:finalSlope*a,
finalValue:finalValue*a]]
END;
BreakPiece: PROCEDURE [f: Handle, d: REAL] RETURNS [Handle] =
BEGIN
IF d<=f.first.domainStart OR d>=f.first.domainEnd THEN RETURN [f]
ELSE BEGIN OPEN f.first;
jointValue,jointSlope: REAL;
[value:jointValue, deriv:jointSlope] ¬ EvalAll[f,d];
RETURN [CONS[[domainStart:domainStart, domainEnd:d,
initValue:initValue,
initSlope:initSlope,
finalSlope:jointSlope,
finalValue:jointValue],
CONS[[domainStart:d, domainEnd:domainEnd,
initValue:jointValue,
initSlope:jointSlope,
finalSlope:finalSlope,
finalValue:finalValue],f.rest]]]
END
END;
CombinePieces: PROCEDURE [a: REAL, p: PieceRec, b: REAL, q: PieceRec]
RETURNS [PieceRec] =
BEGIN -- p and q should have the same domain
RETURN [[domainStart:p.domainStart, domainEnd:p.domainEnd,
initValue:a*p.initValue+b*q.initValue,
initSlope:a*p.initSlope+b*q.initSlope,
finalSlope:a*p.finalSlope+b*q.finalSlope,
finalValue:a*p.finalValue+b*q.finalValue]]
END;
Scale: PUBLIC PROCEDURE [f: Handle, a: REAL] RETURNS [Handle] =
BEGIN
IF a=0 OR f=NIL THEN RETURN [NIL]
ELSE IF a=1 THEN RETURN [f]
ELSE RETURN [CONS[ScalePiece[f.first,a],Scale[f.rest,a]]]
END;
EnumerateCommonPieces: PUBLIC PROCEDURE [f: Handle, g: Handle, P: PiecewiseCubic.PieceProc] =
BEGIN
WHILE f#NIL OR g#NIL DO
WHILE f#NIL AND (g=NIL OR f.first.domainStart<g.first.domainStart) DO
f ¬ BreakPiece[f,f.first.domainStart];
P[f.first, [f.first.domainStart, f.first.domainEnd, 0, 0, 0, 0]];
f ¬ f.rest;
ENDLOOP;
WHILE g#NIL AND (f=NIL OR g.first.domainStart<f.first.domainStart) DO
g ¬ BreakPiece[g,g.first.domainStart];
P[[g.first.domainStart, g.first.domainEnd, 0, 0, 0, 0], g.first];
g ¬ g.rest;
ENDLOOP;
WHILE f#NIL AND g#NIL AND (f.first.domainStart=g.first.domainStart) DO
f ¬ BreakPiece[f,g.first.domainEnd];
g ¬ BreakPiece[g,f.first.domainEnd];
P[f.first,g.first];
f ¬ f.rest;
g ¬ g.rest;
ENDLOOP;
ENDLOOP;
END;
Combine: PUBLIC PROCEDURE [a: REAL, f: Handle, b: REAL, g: Handle]
RETURNS [Handle] =
BEGIN
IF f=NIL THEN RETURN[Scale[g,b]]
ELSE IF g=NIL THEN RETURN[Scale[f,a]]
ELSE IF g.first.domainStart>f.first.domainStart THEN
RETURN[Combine[b,g,a,f]]
ELSE IF g.first.domainStart=f.first.domainStart THEN
BEGIN
IF g.first.domainEnd=f.first.domainEnd THEN
RETURN[CONS[CombinePieces[a,f.first,b,g.first],
Combine[a,f.rest,b,g.rest]]]
ELSE IF g.first.domainEnd<f.first.domainEnd THEN
RETURN[Combine[a,BreakPiece[f,g.first.domainEnd],b,g]]
ELSE RETURN[Combine[a,f,b,BreakPiece[g,f.first.domainEnd]]]
END
ELSE IF g.first.domainStart<f.first.domainStart THEN
BEGIN
IF g.first.domainEnd<=f.first.domainStart THEN
RETURN[CONS[ScalePiece[g.first,b],Combine[a,f,b,g.rest]]]
ELSE RETURN[Combine[a,f,b,BreakPiece[g,f.first.domainStart]]]
END
ELSE ERROR
END;
Eval: PUBLIC PROCEDURE [f: Handle, t: REAL] RETURNS [REAL] =
BEGIN
WHILE f#NIL AND t>f.first.domainEnd DO f ¬ f.rest ENDLOOP;
IF f=NIL OR t<f.first.domainStart THEN RETURN [0]
ELSE
BEGIN OPEN f.first;
interval: REAL ¬ domainEnd-domainStart;
c0: REAL ¬ initValue;
c1: REAL ¬ initSlope;
c2: REAL ¬ ((-2*initSlope-finalSlope)*interval+3*(finalValue-initValue))/(interval*interval);
c3: REAL ¬ ((initSlope+finalSlope)*interval+2*(initValue-finalValue))/(interval*interval*interval);
t ¬ (t-domainStart);
RETURN[((c3*t+c2)*t+c1)*t+c0]
END
END;
EvalAll: PUBLIC PROCEDURE [f: Handle, t: REAL] RETURNS [value,deriv,derivDeriv,derivDerivDeriv,domainStart,domainEnd: REAL] =
BEGIN
WHILE f#NIL AND t>f.first.domainEnd DO f ¬ f.rest ENDLOOP;
IF f=NIL OR t<f.first.domainStart THEN RETURN [0,0,0,0,0,0]
ELSE
BEGIN OPEN f.first;
interval: REAL ¬ domainEnd-domainStart;
c0: REAL ¬ initValue;
c1: REAL ¬ initSlope;
c2: REAL ¬ ((-2*initSlope-finalSlope)*interval+3*(finalValue-initValue))/(interval*interval);
c3: REAL ¬ ((initSlope+finalSlope)*interval+2*(initValue-finalValue))/(interval*interval*interval);
t ¬ (t-domainStart);
RETURN[((c3*t+c2)*t+c1)*t+c0, (3*c3*t+2*c2)*t+c1, 6*c3*t+2*c2, 6*c3, f.first.domainStart, f.first.domainEnd]
END
END;
END.