DIRECTORY G2dMatrix, G3dBasic, G3dMatrix, G3dPlane, G3dVector, Real, RealFns;

G3dMatrixImpl: CEDAR MONITOR
IMPORTS G2dMatrix, G3dPlane, G3dVector, Real, RealFns
EXPORTS G3dMatrix

~ BEGIN

singular:	PUBLIC ERROR ~ CODE;
Matrix:				TYPE ~ G3dMatrix.Matrix;
MatrixRep:			TYPE ~ G3dMatrix.MatrixRep;
MatrixSequence:		TYPE ~ G3dMatrix.MatrixSequence;
MatrixSequenceRep:	TYPE ~ G3dMatrix.MatrixSequenceRep;
Viewport:				TYPE ~ G3dMatrix.Viewport;
Pair:					TYPE ~ G3dBasic.Pair;
Triple:				TYPE ~ G3dBasic.Triple;
Quad:					TYPE ~ G3dBasic.Quad;
Ray:					TYPE ~ G3dBasic.Ray;
Screen:				TYPE ~ G3dBasic.Screen;

origin:				Triple ~ G3dBasic.origin;
CopyMatrix: PUBLIC PROC [matrix: Matrix, out: Matrix ¬ NIL] RETURNS [Matrix] ~ {
IF matrix = NIL THEN RETURN[NIL];
IF out = NIL THEN out ¬ NEW[MatrixRep];
out­ ¬ matrix­;
RETURN[out];
};

Identity: PUBLIC PROC [out: Matrix ¬ NIL] RETURNS [Matrix] ~ {
IF out = NIL THEN out ¬ NEW[MatrixRep];
out­ ¬ [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]];
RETURN[out];
};

MakeFromTriad: PUBLIC PROC [
v1, v2, v3: Triple,
p: Triple ¬ origin,
unitize: BOOL ¬ FALSE,
out: Matrix ¬ NIL]
RETURNS [Matrix]
~ {
IF out = NIL THEN out ¬ NEW[MatrixRep];
IF unitize THEN {
v1 ¬ G3dVector.Unit[v1];
v2 ¬ G3dVector.Unit[v2];
v3 ¬ G3dVector.Unit[v3];
};
out­ ¬ [[v1.x,v1.y,v1.z,0.0], [v2.x,v2.y,v2.z,0.0], [v3.x,v3.y,v3.z,0.0], [p.x,p.y,p.z,1.0]];
RETURN[out];
};

MakeFromRows: PUBLIC PROC [r0, r1, r2, r3: Quad, out: Matrix ¬ NIL] RETURNS [Matrix] ~ {
IF out = NIL THEN out ¬ NEW[MatrixRep];
out[0] ¬ [r0.x, r0.y, r0.z, r0.w];
out[1] ¬ [r1.x, r1.y, r1.z, r1.w];
out[2] ¬ [r2.x, r2.y, r2.z, r2.w];
out[3] ¬ [r3.x, r3.y, r3.z, r3.w];
RETURN[out];
};

MakeScale: PUBLIC PROC [s: Triple, out: Matrix ¬ NIL] RETURNS [Matrix] ~ {
IF out = NIL THEN out ¬ NEW[MatrixRep];
out­ ¬ [[s.x, 0.0, 0.0, 0.0], [0.0, s.y, 0.0, 0.0], [0.0, 0.0, s.z, 0.0], [0.0, 0.0, 0.0, 1.0]];
RETURN[out];
};

MakeTranslate: PUBLIC PROC [p: Triple, out: Matrix ¬ NIL] RETURNS [Matrix] ~ {
IF out = NIL THEN out ¬ NEW[MatrixRep];
out­ ¬ [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [p.x, p.y, p.z, 1.0]];
RETURN[out];
};

MakeRotate: PUBLIC PROC [
axis: Triple,
theta: REAL,
degrees: BOOL ¬ TRUE,
out: Matrix ¬ NIL]
RETURNS [Matrix]
~ {
sin: REAL ¬ IF degrees THEN RealFns.SinDeg[theta] ELSE RealFns.Sin[theta];
cos: REAL ¬ IF degrees THEN RealFns.CosDeg[theta] ELSE RealFns.Cos[theta];
cos1, sqx, sqy, sqz, xycos1, zxcos1, yzcos1, xsin, ysin, zsin: REAL;
axis ¬ G3dVector.Unit[axis];
[[sqx, sqy, sqz]] ¬ G3dVector.MulVectors[axis, axis];
cos1 ¬ 1.0-cos;
xycos1 ¬ axis.x*axis.y*cos1;	xsin ¬ axis.x*sin;
yzcos1 ¬ axis.y*axis.z*cos1;	ysin ¬ axis.y*sin;
zxcos1 ¬ axis.x*axis.z*cos1;	zsin ¬ axis.z*sin;
IF out = NIL THEN out ¬ NEW[MatrixRep];
out­ ¬ [
[sqx+(1.0-sqx)*cos,	xycos1+zsin,			zxcos1-ysin,			0.0],
[xycos1-zsin,		sqy+(1.0-sqy)*cos,	yzcos1+xsin,			0.0],
[zxcos1+ysin,		yzcos1-xsin,			sqz+(1.0-sqz)*cos,	0.0],
[0.0,					0.0,					0.0,					1.0]];
RETURN[out];
};

MakeRotateAbout: PUBLIC PROC [
axis: Triple,
theta: REAL,
degrees: BOOL ¬ TRUE,
base: Triple ¬ origin,
out: Matrix ¬ NIL]
RETURNS [Matrix]
~ {
rot: Matrix ¬ ObtainMatrix[];
tran: Matrix ¬ ObtainMatrix[];
mul: Matrix ¬ ObtainMatrix[];
out ¬ Translate[
			Mul[
			MakeTranslate[[-base.x, -base.y, -base.z], tran],
			MakeRotate[axis, theta, degrees, rot], mul],
			base, out];
ReleaseMatrix[rot];
ReleaseMatrix[tran];
ReleaseMatrix[mul];
RETURN[out];
};

MakePerspective: PUBLIC PROC [nInv, fInv, fov: REAL, out: Matrix ¬ NIL]
RETURNS [Matrix]
~ {
IF fov = 0.0
THEN RETURN[Identity[out]]
ELSE {
d: REAL ¬ 1.0/nInv;
s: REAL ¬ d*RealFns.TanDeg[0.5*fov];
out ¬ Identity[out];
out[0][0] ¬ 1.0/s;
out[1][1] ¬ 1.0/s;
out[2][2] ¬ 1.0/(d*(1.0-d*fInv));
out[2][3] ¬ nInv;
out[3][2] ¬ -1.0/(1.0-d*fInv);
out[3][3] ¬ 0.0;
RETURN[out];
};
};

HasPerspective: PUBLIC PROC [mat: Matrix] RETURNS [BOOL] ~ {
IF mat = NIL THEN RETURN[FALSE];
RETURN[mat[0][3] # 0. OR mat[1][3] # 0. OR mat[2][3] # 0. OR mat[3][3] # 1.0];
};

MakeVectorTransform: PUBLIC PROC [pointTransform: Matrix, out: Matrix ¬ NIL]
RETURNS [n: Matrix]
~ {
IF out = NIL THEN out ¬ NEW[MatrixRep];
out ¬ Transpose[Adjoint[pointTransform], out];
FOR k: INT IN [0..4) DO out[3][k] ¬ out[k][3] ¬ 0.0;  ENDLOOP;
out[3][3] ¬ 1.0;
RETURN[out];
};
TransformH: PUBLIC PROC [p: Triple, mat: Matrix] RETURNS [Quad] ~ {
RETURN[[
p.x*mat[0][0]+p.y*mat[1][0]+p.z*mat[2][0]+mat[3][0],
p.x*mat[0][1]+p.y*mat[1][1]+p.z*mat[2][1]+mat[3][1],
p.x*mat[0][2]+p.y*mat[1][2]+p.z*mat[2][2]+mat[3][2],
p.x*mat[0][3]+p.y*mat[1][3]+p.z*mat[2][3]+mat[3][3]]];
};

Transform: PUBLIC PROC [p: Triple, mat: Matrix] RETURNS [Triple] ~ {
q: Quad;
IF NOT HasPerspective[mat] THEN RETURN[[
p.x*mat[0][0]+p.y*mat[1][0]+p.z*mat[2][0]+mat[3][0],
p.x*mat[0][1]+p.y*mat[1][1]+p.z*mat[2][1]+mat[3][1],
p.x*mat[0][2]+p.y*mat[1][2]+p.z*mat[2][2]+mat[3][2]]];
q ¬ TransformH[p, mat];
RETURN[IF q.w = 1.0 THEN [q.x, q.y, q.z] ELSE [q.x/q.w, q.y/q.w, q.z/q.w]];
};

TransformD: PUBLIC PROC [p: Triple, mat: Matrix] RETURNS [Pair] ~ {
w: REAL;
pp: Pair ¬ [
p.x*mat[0][0]+p.y*mat[1][0]+p.z*mat[2][0]+mat[3][0],
p.x*mat[0][1]+p.y*mat[1][1]+p.z*mat[2][1]+mat[3][1]];
IF NOT HasPerspective[mat] THEN RETURN[pp];
w ¬ p.x*mat[0][3]+p.y*mat[1][3]+p.z*mat[2][3]+mat[3][3];
RETURN[IF w = 0.0 THEN [0., 0.] ELSE IF w = 1.0 THEN pp ELSE [pp.x/w, pp.y/w]];
};

TransformPair: PUBLIC PROC [p: Pair, mat: Matrix] RETURNS [Triple] ~ {
w: REAL;
t: Triple ¬ [
p.x*mat[0][0]+p.y*mat[1][0]+mat[3][0],
p.x*mat[0][1]+p.y*mat[1][1]+mat[3][1],
p.x*mat[0][2]+p.y*mat[1][2]+mat[3][2]];
IF NOT HasPerspective[mat] THEN RETURN[t];
w ¬ p.x*mat[0][3]+p.y*mat[1][3]+mat[3][3];
RETURN[IF w = 1.0 THEN t ELSE [t.x/w, t.y/w, t.z/w]];
};

TransformPairD: PUBLIC PROC [p: Pair, mat: Matrix] RETURNS [Pair] ~ {
w: REAL;
pp: Pair ¬ [p.x*mat[0][0]+p.y*mat[1][0]+mat[3][0], p.x*mat[0][1]+p.y*mat[1][1]+mat[3][1]];
IF NOT HasPerspective[mat] THEN RETURN[pp];
w ¬ p.x*mat[0][3]+p.y*mat[1][3]+mat[3][3];
RETURN[IF w = 1.0 THEN pp ELSE [pp.x/w, pp.y/w]];
};

TransformQuad: PUBLIC PROC [q: Quad, mat: Matrix] RETURNS [Quad] ~ {
RETURN[[
q.x*mat[0][0]+q.y*mat[1][0]+q.z*mat[2][0]+q.w*mat[3][0],
q.x*mat[0][1]+q.y*mat[1][1]+q.z*mat[2][1]+q.w*mat[3][1],
q.x*mat[0][2]+q.y*mat[1][2]+q.z*mat[2][2]+q.w*mat[3][2],
q.x*mat[0][3]+q.y*mat[1][3]+q.z*mat[2][3]+q.w*mat[3][3]]];
};


TransformVec: PUBLIC PROC [vec: Triple, mat: Matrix] RETURNS [Triple] ~ {
RETURN[[
vec.x*mat[0][0]+vec.y*mat[1][0]+vec.z*mat[2][0],
vec.x*mat[0][1]+vec.y*mat[1][1]+vec.z*mat[2][1],
vec.x*mat[0][2]+vec.y*mat[1][2]+vec.z*mat[2][2]]];
};

TransformVecDiffS: PUBLIC PROC [vec: Triple, mat: Matrix] RETURNS [Triple] ~ {
RETURN[Transform[vec, Cofactors[mat]]]; -- or check if sum of squares of rows are all equal.
};

TransformPlane: PUBLIC PROC [plane: Quad, mat: Matrix, matInverse: Matrix ¬ NIL]
RETURNS [Quad]
~ {
m: Matrix ¬ IF matInverse # NIL THEN matInverse ELSE Adjoint[mat];
RETURN[PremultiplyPlaneByMatrix[plane, m]];
};

PremultiplyPlaneByMatrix: PUBLIC PROC [plane: Quad, m: Matrix] RETURNS [Quad] ~ {
RETURN[[
plane.x*m[0][0]+plane.y*m[0][1]+plane.z*m[0][2]+plane.w*m[0][3],
plane.x*m[1][0]+plane.y*m[1][1]+plane.z*m[1][2]+plane.w*m[1][3],
plane.x*m[2][0]+plane.y*m[2][1]+plane.z*m[2][2]+plane.w*m[2][3],
plane.x*m[3][0]+plane.y*m[3][1]+plane.z*m[3][2]+plane.w*m[3][3]]];
};

TransformByViewport: PUBLIC PROC [p: Pair, vp: Viewport] RETURNS [x: Pair] ~ {
x ¬ [vp.scale.y*p.x+vp.translate.x, vp.scale.y*p.y+vp.translate.y]; -- uniform scale
};

Scale: PUBLIC PROC [in: Matrix, s: REAL, out: Matrix ¬ NIL] RETURNS [Matrix] ~ {
scratch: Matrix ¬ ObtainMatrix[];
out ¬ Mul[in, MakeScale[[s, s, s], scratch], out];
ReleaseMatrix[scratch];
RETURN[out];
};
DiffScale: PUBLIC PROC [in: Matrix, s: Triple, out: Matrix ¬ NIL] RETURNS [Matrix] ~ {
scratch: Matrix ¬ ObtainMatrix[];
out ¬ Mul[in, MakeScale[s, scratch], out];
ReleaseMatrix[scratch];
RETURN[out];
};
Rotate: PUBLIC PROC [
in: Matrix,
axis: Triple,
theta: REAL,
degrees: BOOL ¬ TRUE,
base: Triple ¬ origin,
out: Matrix ¬ NIL]
RETURNS [Matrix]
~ {
scratch: Matrix ¬ ObtainMatrix[];
out ¬ Mul[in, MakeRotateAbout[axis, theta, degrees, base, scratch], out];
ReleaseMatrix[scratch];
RETURN[out];
};

Translate: PUBLIC PROC [in: Matrix, p: Triple, out: Matrix ¬ NIL] RETURNS [Matrix] ~ {
scratch: Matrix ¬ ObtainMatrix[];
out ¬ Mul[in, MakeTranslate[p, scratch], out];
ReleaseMatrix[scratch];
RETURN[out];
};
LocalScale: PUBLIC PROC [in: Matrix, s: REAL, out: Matrix ¬ NIL] RETURNS [Matrix] ~ {
scratch: Matrix ¬ ObtainMatrix[];
out ¬ Mul[MakeScale[[s, s, s], scratch], in, out];
ReleaseMatrix[scratch];
RETURN[out];
};

LocalDiffScale: PUBLIC PROC [in: Matrix, s: Triple, out: Matrix ¬ NIL] RETURNS [Matrix] ~ {
scratch: Matrix ¬ ObtainMatrix[];
out ¬ Mul[MakeScale[s, scratch], in, out];
ReleaseMatrix[scratch];
RETURN[out];
};

LocalRotate: PUBLIC PROC [
in: Matrix,
axis: Triple,
theta: REAL,
degrees: BOOL ¬ TRUE,
base: Triple ¬ origin,
out: Matrix ¬ NIL]
RETURNS [Matrix]
~ {
scratch: Matrix ¬ ObtainMatrix[];
out ¬ Mul[MakeRotateAbout[axis, theta, degrees, base, scratch], in, out];
ReleaseMatrix[scratch];
RETURN[out];
};

LocalTranslate: PUBLIC PROC [in: Matrix, p: Triple, out: Matrix ¬ NIL] RETURNS [Matrix] ~ {
scratch: Matrix ¬ ObtainMatrix[];
out ¬ Mul[MakeTranslate[p, scratch], in];
ReleaseMatrix[scratch];
RETURN[out];
};
ScalarMul: PROC [in: Matrix, s: REAL, out: Matrix ¬ NIL] RETURNS [Matrix] ~ {
IF out = NIL THEN out ¬ NEW[MatrixRep];
FOR i: [0..4) IN [0..4) DO
FOR j: [0..4) IN [0..4) DO
out­[i][j] ¬ s*in[i][j];
ENDLOOP;
ENDLOOP;
RETURN[out];
};

Mul: PUBLIC PROC [left, rite: Matrix, out: Matrix ¬ NIL] RETURNS [Matrix] ~ {
ret: MatrixRep;
IF left = NIL OR rite = NIL THEN RETURN[out];

ret[0][0] ¬ left[0][0]*rite[0][0]+left[0][1]*rite[1][0]+left[0][2]*rite[2][0]+left[0][3]*rite[3][0];
ret[0][1] ¬ left[0][0]*rite[0][1]+left[0][1]*rite[1][1]+left[0][2]*rite[2][1]+left[0][3]*rite[3][1];
ret[0][2] ¬ left[0][0]*rite[0][2]+left[0][1]*rite[1][2]+left[0][2]*rite[2][2]+left[0][3]*rite[3][2];
ret[0][3] ¬ left[0][0]*rite[0][3]+left[0][1]*rite[1][3]+left[0][2]*rite[2][3]+left[0][3]*rite[3][3];

ret[1][0] ¬ left[1][0]*rite[0][0]+left[1][1]*rite[1][0]+left[1][2]*rite[2][0]+left[1][3]*rite[3][0];
ret[1][1] ¬ left[1][0]*rite[0][1]+left[1][1]*rite[1][1]+left[1][2]*rite[2][1]+left[1][3]*rite[3][1];
ret[1][2] ¬ left[1][0]*rite[0][2]+left[1][1]*rite[1][2]+left[1][2]*rite[2][2]+left[1][3]*rite[3][2];
ret[1][3] ¬ left[1][0]*rite[0][3]+left[1][1]*rite[1][3]+left[1][2]*rite[2][3]+left[1][3]*rite[3][3];

ret[2][0] ¬ left[2][0]*rite[0][0]+left[2][1]*rite[1][0]+left[2][2]*rite[2][0]+left[2][3]*rite[3][0];
ret[2][1] ¬ left[2][0]*rite[0][1]+left[2][1]*rite[1][1]+left[2][2]*rite[2][1]+left[2][3]*rite[3][1];
ret[2][2] ¬ left[2][0]*rite[0][2]+left[2][1]*rite[1][2]+left[2][2]*rite[2][2]+left[2][3]*rite[3][2];
ret[2][3] ¬ left[2][0]*rite[0][3]+left[2][1]*rite[1][3]+left[2][2]*rite[2][3]+left[2][3]*rite[3][3];

ret[3][0] ¬ left[3][0]*rite[0][0]+left[3][1]*rite[1][0]+left[3][2]*rite[2][0]+left[3][3]*rite[3][0];
ret[3][1] ¬ left[3][0]*rite[0][1]+left[3][1]*rite[1][1]+left[3][2]*rite[2][1]+left[3][3]*rite[3][1];
ret[3][2] ¬ left[3][0]*rite[0][2]+left[3][1]*rite[1][2]+left[3][2]*rite[2][2]+left[3][3]*rite[3][2];
ret[3][3] ¬ left[3][0]*rite[0][3]+left[3][1]*rite[1][3]+left[3][2]*rite[2][3]+left[3][3]*rite[3][3];

IF out = NIL THEN out ¬ NEW[MatrixRep];
out­ ¬ ret;
RETURN[out];
};

Cat: PUBLIC PROC [m1, m2: Matrix, m3, m4, m5, m6, out: Matrix ¬ NIL] RETURNS [Matrix] ~ {
IF out = NIL THEN out ¬ NEW[MatrixRep];
out ¬ Mul[m1, m2, out];
IF m3 # NIL THEN out ¬ Mul[out, m3, out];
IF m4 # NIL THEN out ¬ Mul[out, m4, out];
IF m5 # NIL THEN out ¬ Mul[out, m5, out];
IF m6 # NIL THEN out ¬ Mul[out, m6, out];
RETURN[out];
};

Invert: PUBLIC PROC [in: Matrix, out: Matrix ¬ NIL] RETURNS [ret: Matrix] ~ {
ret ¬ IF out = NIL THEN NEW[MatrixRep] ELSE out;
IF in[3][0] = 0.0 AND in[3][1] = 0.0 AND in[3][2] = 0.0 AND in[3][3] = 1.0
THEN {	-- a little trick from Martin Newell (or was it Jim Blinn?)
inv3x3: G2dMatrix.Matrix ¬ G2dMatrix.Invert[[
[in[0][0], in[0][1], in[0][2]],
[in[1][0], in[1][1], in[1][2]],
[in[2][0], in[2][1], in[2][2]]]];
q: Triple ¬ G2dMatrix.Transform[[-in[3][0], -in[3][1], -in[3][2]], inv3x3];
ret[0] ¬ [inv3x3.row1.x, inv3x3.row1.y, inv3x3.row1.z, 0.0];
ret[1] ¬ [inv3x3.row2.x, inv3x3.row2.y, inv3x3.row2.z, 0.0];
ret[2] ¬ [inv3x3.row3.x, inv3x3.row3.y, inv3x3.row3.z, 0.0];
ret[3] ¬ [q.x, q.y, q.z, 1.0];
}
ELSE {
adj: Matrix ¬ Adjoint[in, ObtainMatrix[]];
det: REAL ¬ in[0][0]*adj[0][0]+in[0][1]*adj[1][0]+in[0][2]*adj[2][0]+in[0][3]*adj[3][0];
IF det = 0.0 THEN ERROR singular;
ret ¬ ScalarMul[adj, 1.0/det, ret];
ReleaseMatrix[adj];
};
};

Adjoint: PUBLIC PROC [in: Matrix, out: Matrix ¬ NIL] RETURNS [Matrix] ~ {
scratch: Matrix ¬ ObtainMatrix[];
out ¬ Transpose[Cofactors[in, scratch], out];
ReleaseMatrix[scratch];
RETURN[out];
};

Det3x3: PROC [m: G2dMatrix.Matrix] RETURNS [REAL] ~ G2dMatrix.Determinant;

Cofactors: PUBLIC PROC [in: Matrix, out: Matrix ¬ NIL] RETURNS [Matrix] ~ {
m: Matrix ¬ in;
IF out = NIL THEN out ¬ NEW[MatrixRep];
out[0][0] ¬  Det3x3[[[m[1][1],m[1][2],m[1][3]],[m[2][1],m[2][2],m[2][3]],[m[3][1],m[3][2],m[3][3]]]]; 
out[0][1] ¬ -Det3x3[[[m[1][0],m[1][2],m[1][3]],[m[2][0],m[2][2],m[2][3]],[m[3][0],m[3][2],m[3][3]]]];
out[0][2] ¬  Det3x3[[[m[1][0],m[1][1],m[1][3]],[m[2][0],m[2][1],m[2][3]],[m[3][0],m[3][1],m[3][3]]]];
out[0][3] ¬ -Det3x3[[[m[1][0],m[1][1],m[1][2]],[m[2][0],m[2][1],m[2][2]],[m[3][0],m[3][1],m[3][2]]]];
out[1][0] ¬ -Det3x3[[[m[0][1],m[0][2],m[0][3]],[m[2][1],m[2][2],m[2][3]],[m[3][1],m[3][2],m[3][3]]]];
out[1][1] ¬  Det3x3[[[m[0][0],m[0][2],m[0][3]],[m[2][0],m[2][2],m[2][3]],[m[3][0],m[3][2],m[3][3]]]];
out[1][2] ¬ -Det3x3[[[m[0][0],m[0][1],m[0][3]],[m[2][0],m[2][1],m[2][3]],[m[3][0],m[3][1],m[3][3]]]];
out[1][3] ¬  Det3x3[[[m[0][0],m[0][1],m[0][2]],[m[2][0],m[2][1],m[2][2]],[m[3][0],m[3][1],m[3][2]]]];
out[2][0] ¬  Det3x3[[[m[0][1],m[0][2],m[0][3]],[m[1][1],m[1][2],m[1][3]],[m[3][1],m[3][2],m[3][3]]]];
out[2][1] ¬ -Det3x3[[[m[0][0],m[0][2],m[0][3]],[m[1][0],m[1][2],m[1][3]],[m[3][0],m[3][2],m[3][3]]]];
out[2][2] ¬  Det3x3[[[m[0][0],m[0][1],m[0][3]],[m[1][0],m[1][1],m[1][3]],[m[3][0],m[3][1],m[3][3]]]];
out[2][3] ¬ -Det3x3[[[m[0][0],m[0][1],m[0][2]],[m[1][0],m[1][1],m[1][2]],[m[3][0],m[3][1],m[3][2]]]];
out[3][0] ¬ -Det3x3[[[m[0][1],m[0][2],m[0][3]],[m[1][1],m[1][2],m[1][3]],[m[2][1],m[2][2],m[2][3]]]];
out[3][1] ¬  Det3x3[[[m[0][0],m[0][2],m[0][3]],[m[1][0],m[1][2],m[1][3]],[m[2][0],m[2][2],m[2][3]]]];
out[3][2] ¬ -Det3x3[[[m[0][0],m[0][1],m[0][3]],[m[1][0],m[1][1],m[1][3]],[m[2][0],m[2][1],m[2][3]]]];
out[3][3] ¬  Det3x3[[[m[0][0],m[0][1],m[0][2]],[m[1][0],m[1][1],m[1][2]],[m[2][0],m[2][1],m[2][2]]]];
RETURN[out];
};

Determinant: PUBLIC PROC [m: Matrix] RETURNS [REAL] ~ {
RETURN[
m[0][0]*
Det3x3[[[m[1][1],m[1][2],m[1][3]],[m[2][1],m[2][2],m[2][3]],[m[3][1],m[3][2],m[3][3]]]] -
m[0][1]*
Det3x3[[[m[1][0],m[1][2],m[1][3]],[m[2][0],m[2][2],m[2][3]],[m[3][0],m[3][2],m[3][3]]]] +
m[0][2]*
Det3x3[[[m[1][0],m[1][1],m[1][3]],[m[2][0],m[2][1],m[2][3]],[m[3][0],m[3][1],m[3][3]]]] -
m[0][3]*
Det3x3[[[m[1][0],m[1][1],m[1][2]],[m[2][0],m[2][1],m[2][2]],[m[3][0],m[3][1],m[3][2]]]]];
};

Transpose: PUBLIC PROC [in: Matrix, out: Matrix ¬ NIL] RETURNS [Matrix] ~ {
IF out = NIL THEN out ¬ NEW[MatrixRep];
FOR i: [0..4) IN [0..4) DO
out[i][i] ¬ in[i][i];
FOR j: [0..4) IN [0..i) DO
ij: REAL ¬ in[i][j];
ji: REAL ¬ in[j][i];
out[i][j] ¬ ji;
out[j][i] ¬ ij;
ENDLOOP;
ENDLOOP;
RETURN[out];
};

Atan: PROC [y, x: REAL] RETURNS [a: REAL] ~ INLINE {a ¬ RealFns.ArcTanDeg[y, x]};
Asin: PROC [x: REAL] RETURNS [a: REAL] ~ INLINE {a ¬ Atan[x, RealFns.SqRt[1.0-x*x]]};

ExtractRotate: PUBLIC PROC [m: Matrix] RETURNS [Triple] ~ {
RETURN[IF ABS[m[0][2]] # 1.0
THEN [Atan[m[1][2], m[2][2]], -Asin[m[0][2]], Atan[m[0][1], m[0][0]]]
ELSE  [Atan[m[1][0], m[2][0]], -Asin[m[0][2]], 0.]];
};

InTermsOf: PUBLIC PROC [A, B: Matrix, out: Matrix ¬ NIL] RETURNS [Matrix] ~ {
scratch: Matrix ¬ ObtainMatrix[];
out ¬ Invert[A, scratch];
out ¬ Mul[scratch, B, out];
ReleaseMatrix[scratch];
RETURN[out];
};
CopyMatrixSequence: PUBLIC PROC [matrices: MatrixSequence]
RETURNS [MatrixSequence] ~ {
copy: MatrixSequence ¬ NIL;
IF matrices # NIL THEN {
copy ¬ NEW[MatrixSequenceRep[matrices.length]];
copy.length ¬ matrices.length;
FOR n: NAT IN [0..matrices.length) DO copy[n] ¬ CopyMatrix[matrices[n]]; ENDLOOP;
};
RETURN[copy];
};

AddToMatrixSequence: PUBLIC PROC [matrices: MatrixSequence, matrix: Matrix]
RETURNS [MatrixSequence]
~ {
IF matrices = NIL THEN matrices ¬ NEW[MatrixSequenceRep[1]];
IF matrices.length = matrices.maxLength THEN matrices ¬ LengthenMatrixSequence[matrices];
matrices[matrices.length] ¬ matrix;
matrices.length ¬ matrices.length+1;
RETURN[matrices];
};

LengthenMatrixSequence: PUBLIC PROC [matrices: MatrixSequence, amount: REAL ¬ 1.3]
RETURNS [new: MatrixSequence] ~ {
newLength: NAT ¬ MAX[Real.Ceiling[amount*matrices.maxLength], 3];
new ¬ NEW[MatrixSequenceRep[newLength]];
FOR i: NAT IN [0..matrices.length) DO new[i] ¬ matrices[i]; ENDLOOP;
new.length ¬ matrices.length;
};
nScratchMatrices: NAT ~ 6;
scratchMatrices: ARRAY [0..nScratchMatrices) OF Matrix ¬ ALL[NIL];

ObtainMatrix: PUBLIC ENTRY PROC RETURNS [Matrix] ~ {
FOR i: NAT IN [0..nScratchMatrices) DO
matrix: Matrix ~ scratchMatrices[i];
IF matrix # NIL THEN {
scratchMatrices[i] ¬ NIL;
RETURN[matrix];
};
ENDLOOP;
RETURN[NEW[MatrixRep]];
};

ReleaseMatrix: PUBLIC ENTRY PROC [matrix: Matrix] ~ {
FOR i: NAT IN [0..nScratchMatrices) DO
IF scratchMatrices[i] = NIL THEN {
scratchMatrices[i] ¬ matrix;
RETURN;
};
ENDLOOP;
};
GetScreen: PUBLIC PROC [
point: Triple,
view: Matrix,
hasPerspective: BOOL,
viewport: Viewport ¬ []]
RETURNS [s: Screen]
~ {
IF hasPerspective
THEN {
q: Quad ¬ s.quad ¬ TransformH[point, view];
q.x ¬ s.quad.x ¬ q.x*viewport.aspectRecip; -- deform to square NDC clip-space
s.pos ¬ SELECT q.w FROM 0, 1 => [q.x, q.y], ENDCASE => [q.x/q.w, q.y/q.w];
s.l ¬ q.w+q.x < 0.0;		-- test left clipping
s.r ¬ q.w-q.x < 0.0;		-- test right clipping
s.b ¬ q.w+q.y < 0.0;		-- test bottom clipping
s.t ¬ q.w-q.y < 0.0;		-- test top clipping
s.n ¬ q.z < 0.0;			-- test near clipping
IF (s.visible ¬ NOT (s.l OR s.r OR s.b OR s.t OR s.n)) AND viewport # [] THEN {
s.pos.x ¬ viewport.scale.x*s.pos.x+viewport.translate.x; -- aspectRecip*scale.x=scale.y,
s.pos.y ¬ viewport.scale.y*s.pos.y+viewport.translate.y; -- so uniform scale
};
}
ELSE {
z: REAL ¬ point.x*view[0][2]+point.y*view[1][2]+point.z*view[2][2]+view[3][2];
s.pos ¬ [
point.x*view[0][0]+point.y*view[1][0]+point.z*view[2][0]+view[3][0],
point.x*view[0][1]+point.y*view[1][1]+point.z*view[2][1]+view[3][1]];
IF (s.visible ¬ NOT (s.n ¬ z < 0.0)) AND viewport # [] THEN {
s.pos.x ¬ viewport.scale.y*s.pos.x+viewport.translate.x; -- uniform scale
s.pos.y ¬ viewport.scale.y*s.pos.y+viewport.translate.y;
};
};
s.intPos ¬ [Real.Round[s.pos.x], Real.Round[s.pos.y]];
}; 
TripleFromScreen: PUBLIC PROC [p: Pair, view: Matrix, viewInverse: Matrix ¬ NIL]
RETURNS [t: Triple]
~ {
x: Matrix ¬ IF viewInverse # NIL THEN viewInverse ELSE Invert[view, ObtainMatrix[]];
t ¬ Transform[[p.x, p.y, 0.0], x];	-- z is irrelevant
IF viewInverse = NIL THEN ReleaseMatrix[x];
};

TripleFromScreenAndPlane: PUBLIC PROC [p: Pair, plane: Quad, view: Matrix]
RETURNS [t: Triple]
~ {
vScreen: Quad ¬ [1.0, 0.0, 0.0, -p.x];		-- vertical plane through the given screen point
hScreen: Quad ¬ [0.0, 1.0, 0.0, -p.y];		-- horizontal plane through the given screen point
vWorld: Quad ¬ PremultiplyPlaneByMatrix[vScreen, view];
hWorld: Quad ¬ PremultiplyPlaneByMatrix[hScreen, view];
t ¬ G3dPlane.IntersectionOf3[[vWorld.x, vWorld.y, vWorld.z, vWorld.w, FALSE], [hWorld.x, hWorld.y, hWorld.z, hWorld.w, FALSE], [plane.x, plane.y, plane.z, plane.w, FALSE]];
};

LineFromScreen: PUBLIC PROC [p: Pair, view: Matrix, viewInverse: Matrix ¬ NIL]
RETURNS [r: Ray]
~ {
x: Matrix ¬ IF viewInverse # NIL THEN viewInverse ELSE Invert[view, ObtainMatrix[]];
p0: Triple ¬ TripleFromScreen[p, view, x];
p1: Triple ¬ TripleFromScreen[[0.0, 0.0], view, x];
r ¬ [p0, G3dVector.Sub[p1, p0]];
IF viewInverse = NIL THEN ReleaseMatrix[x];
};

END.
..
by Crow:		/Cyan/Imaging/6.0/ThreeDWorld/Linear3DImpl.mesa
by Bier:		/Cedar/CedarChest6.0/SolidModeler/Matrix3dImpl.mesa
by Wyatt:	/Cedar/Cedar6.0/Imager/ImagerTransformationImpl.mesa
MakePerspective: PUBLIC PROC [dInv, fInv, fov: REAL, out: Matrix ¬ NIL] RETURNS [Matrix]
~ {
tan: REAL ¬ RealFns.TanDeg[0.5*fov];
out ¬ Identity[out];
out[2][2] ¬ tan*(1.0+fInv/MAX[0.0001, ABS[dInv]]);
out[2][3] ¬ tan;
RETURN[out];
};

	
G3dMatrixImpl.mesa
Copyright Σ 1984, 1992 by Xerox Corporation.  All rights reserved.
Bier, June 1982
Crow, October 31, 1985 12:31:44 pm PST
Bloomenthal, October 21, 1992 4:55 pm PDT
Ken Shoemake, August 30, 1989 8:17:57 pm PDT

See ``Principles of Interactive Computer Graphics,'' by Newman and Sproull.
Left-handed eyespace coordinate system (right, up, depth = into screen);
right-handed world coordinate system (longitude, latitude, altitude).
Ken Fishkin, November 6, 1991 1:30 pm PST

Type Declarations
Creation Operations
This performs correctly for right-handed coordinate system and right-handed rotations.
MakePerspective: PUBLIC PROC [nInv, fInv, fov: REAL, out: Matrix _ NIL]
RETURNS [Matrix]
~ {
tan: REAL _ RealFns.TanDeg[0.5*fov];
out _ Identity[out];
out[2][2] _ tan*(nInv+fInv)/(nInv-fInv);
out[2][3] _ tan;
out[3][2] _ -2.*tan/(nInv-fInv);
out[3][3] _ 0.;
RETURN[out];
};

RETURN[mat[2][3] # 0.0];  used to be
Transformation Operations
The following four transformation procedures test the transformation matrix for perspective (i.e.,
the matrix has perspective iff the last column is not [0, 0, 0, 1].  If a perspective element exists,
the return coordinates are first divided by w.

If no differential scaling exists within the transformation matrix, TransformVec may be used.
Otherwise, use TransformVecDiffS.
Concatenation Operations
The next 3 procs post-multiply the matrix, thus transforming in the reference coordinate space.



The next 3 procs pre-multiply the matrix, thus transforming in the local coordinate space.

Mathematical Operations
FOR i: INT IN [0..3] DO
FOR j: INT IN [0..3] DO
ret[i][j] _ left[i][0]*rite[0][j]+left[i][1]*rite[1][j]+left[i][2]*rite[2][j]+left[i][3]*rite[3][j];
ENDLOOP;
ENDLOOP;
Or: transpose; chk for no perspec; check constant sum of squares for rows
Return the ArcTangent of y and x.

Return the ArcSin of x.
Sequence Operations
Scratch Matrix Operations
Point Transformation and Visibility Testing
s.f _ q.w-q.z < 0.0;		-- test far clipping
Three-dimensional Information Derived from the Screen
 The screen planes must be transformed into world space; that is, transformed by the inverse of view.  To transform plane F by matrix M is to premultiply F by the inverse of M.  Thus, if M = view; transform plane F by M-1, yielding F'.   Thus, F' = MF.
References
The Wayside
dInverse = 0 for orthographic projection:
Basically, with an infinite far clipping plane, d is irrelevant.
Κz•NewlineDelimiter

–
"cedarcode" style™™šœ	Οeœ6™BJ™J™&J™)J™,J™—™KJ™HJ™E—™)J™—JšΟk	œD˜MJ˜—šΡbln
œžœ
ž˜Jšžœ.˜5Jšžœ
˜J˜—šœž˜J˜Jšœ
žœ
žœžœ
˜—
headš
Οl™Jšœžœ˜#Jšœ
žœ˜(Jšœžœ˜1Jšœžœ˜6Jšœ
žœ˜'Jšœ
žœ˜Jšœžœ˜"Jšœ
žœ˜Jšœ	žœ˜šœžœ˜"J˜—J˜$—š
 ™š
Οn
œžœ
žœ žœžœ
˜PJš
žœ
žœ
žœ
žœ
žœ˜!Jšžœžœ
žœžœ˜'J˜Jšžœ˜J˜J˜—š
‘œžœ
žœžœžœ
˜>Jšžœžœ
žœžœ˜'J˜`Jšžœ˜J˜J˜—š‘
œžœ
žœ˜Jš
œ˜J˜Jšœ	žœžœ
˜Jšœžœ
˜Jšžœ	˜J˜Jšžœžœ
žœžœ˜'šžœ	žœ˜J˜J˜J˜J˜—J˜]Jšžœ˜J˜J˜—š
‘œžœ
žœ&žœ
žœ
˜XJšžœžœ
žœžœ˜'J˜"J˜"J˜"J˜"Jšžœ˜J˜J˜—š
‘	œžœ
žœžœžœ
˜JJšžœžœ
žœžœ˜'J˜`Jšžœ˜J˜J˜—š
‘
œžœ
žœžœžœ
˜NJšžœžœ
žœžœ˜'J˜`Jšžœ˜J˜J˜—š‘
œžœ
žœ˜J˜
Jšœžœ
˜Jšœ	žœžœ
˜Jšœžœ
˜Jšžœ	˜J˜J™VJš	œžœžœ	žœžœ˜JJš	œžœžœ	žœžœ˜JJšœ?žœ
˜DJ˜J˜5J˜J˜/J˜/J˜/Jšžœžœ
žœžœ˜'˜J˜7J˜6J˜6J˜"—Jšžœ˜J˜J˜—š‘œžœ
žœ˜J˜
Jšœžœ
˜Jšœ	žœžœ
˜J˜Jšœžœ
˜Jšžœ	˜J˜J˜J˜J˜˜˜J˜4J˜/—J˜—J˜J˜J˜Jšžœ˜J˜J˜—š
‘œžœ
žœžœžœ
˜GJšžœ	˜J˜šžœ
˜Jšžœ
žœ˜šžœ˜Jšœžœ˜Jšœžœ˜$J˜J˜J˜J˜!J˜J˜J˜Jšžœ˜J˜——J˜J˜—š‘œžœžœžœ
™GJšžœ	™J™Jšœžœ™$J™J™(J™J™ J™Jšžœ™J™J™—š
‘œžœ
žœžœžœ˜<Jš
žœžœ
žœ
žœ
žœ˜ Jšžœžœ
™$Jšžœžœžœžœ˜NJ˜J˜—š‘œžœ
žœ(žœ
˜LJšžœ˜Jš
œ˜Jšžœžœ
žœžœ˜'J˜.Jš
žœžœ
žœžœžœ
˜>J˜Jšžœ˜J˜——š
 ™š‘
œžœ
žœžœ˜Cšžœ˜J˜4J˜4J˜4J˜6—J˜J˜—J™bJ™e™.J™—š‘	œžœ
žœžœ
˜DJ˜šžœ
žœžœ
žœ˜(J˜4J˜4J˜6—J˜Jšžœ
žœžœžœ˜KJ˜J˜—š‘
œžœ
žœžœ˜CJšœžœ
˜˜J˜4J˜5—Jšžœ
žœžœ
žœ˜+J˜8Jšžœ
žœ	žœ
žœ
žœ	žœžœ˜OJ˜J˜—š‘
œžœ
žœžœ
˜FJšœžœ
˜˜
J˜&J˜&J˜'—Jšžœ
žœžœ
žœ˜*J˜*Jšžœ
žœ	žœžœ˜5J˜J˜—š‘œžœ
žœžœ˜EJšœžœ
˜J˜ZJšžœ
žœžœ
žœ˜+J˜*Jšžœ
žœ	žœžœ˜1J˜J˜—š‘
œžœ
žœžœ˜Dšžœ˜J˜8J˜8J˜8J˜:—J˜J˜—J™]™!J˜—š‘œžœ
žœžœ
˜Išžœ˜J˜0J˜0J˜2—J˜J˜—š‘œžœ
žœžœ
˜NJšžœ"Οc4˜\J˜J˜—š‘œžœ
žœ1žœ
˜PJšžœ˜J˜Jš	œžœžœ
žœžœ˜BJšžœ%˜+J˜J˜—š‘œžœ
žœžœ˜Qšžœ˜J˜@J˜@J˜@J˜B—J˜J˜—š‘œžœžœ˜NJšœD’˜TJ˜——š
 ™™_J˜—š‘œžœ
žœžœžœžœ
˜PJ˜!J˜2J˜Jšžœ˜J˜J™—š
‘	œžœ
žœ'žœžœ
˜VJ˜!J˜*J˜Jšžœ˜J˜J™—š‘œžœ
žœ˜J˜J˜
Jšœžœ
˜Jšœ	žœžœ
˜J˜Jšœžœ
˜Jšžœ	˜J˜J˜!J˜IJ˜Jšžœ˜J˜J˜—š
‘	œžœ
žœ'žœžœ
˜VJ˜!J˜.J˜Jšžœ˜J˜J™—™ZJ™—š‘
œžœ
žœžœžœžœ
˜UJ˜!J˜2J˜Jšžœ˜J˜J˜—š
‘œžœ
žœ'žœžœ
˜[J˜!J˜*J˜Jšžœ˜J˜J˜—š‘œžœ
žœ˜J˜J˜
Jšœžœ
˜Jšœ	žœžœ
˜J˜Jšœžœ
˜Jšžœ	˜J˜J˜!J˜IJ˜Jšžœ˜J˜J˜—š
‘œžœ
žœ'žœžœ
˜[J˜!J˜)J˜Jšžœ˜J˜——š
 ™š
‘	œžœžœžœžœ
˜MJšžœžœ
žœžœ˜'šžœžœž˜šžœžœž˜J˜Jšžœ
˜—Jšžœ
˜—Jšžœ˜J˜J˜—š
‘œžœ
žœ$žœžœ
˜Mšžœžœ
žœž™šžœžœ
žœž™J™dJšžœ
™—Jšžœ
™—J˜Jšžœžœ
žœžœ
žœ
žœ˜-J˜J˜dJ˜dJ˜dJ˜dJ˜J˜dJ˜dJ˜dJ˜dJ˜J˜dJ˜dJ˜dJ˜dJ˜J˜dJ˜dJ˜dJ˜dJ˜Jšžœžœ
žœžœ˜'J˜Jšžœ˜J˜J˜—š
‘œžœ
žœ0žœžœ
˜YJšžœžœ
žœžœ˜'J˜Jšžœžœ
žœ˜)Jšžœžœ
žœ˜)Jšžœžœ
žœ˜)Jšžœžœ
žœ˜)Jšžœ˜J˜J˜—šΠbnΟb
œ
žœ
ž€
œžœžœ˜MJšœžœžœ
žœ
žœžœ˜0šžœžœžœžœ˜Jšžœ’;˜B˜-J˜J˜J˜!—J˜KJ˜<J˜<J˜<J˜J˜
—šžœ˜J˜*JšœžœO˜XJšžœžœ
žœ
˜!J˜#J˜J˜——J˜J˜—š
‘œžœ
žœžœžœ
˜IJ˜!J˜-J˜Jšžœ˜J™IJ˜J˜—š‘œžœžœžœ˜JJ˜—š
‘	œžœ
žœžœžœ
˜KJ˜Jšžœžœ
žœžœ˜'J˜fJ˜eJ˜eJ˜eJ˜eJ˜eJ˜eJ˜eJ˜eJ˜eJ˜eJ˜eJ˜eJ˜eJ˜eJ˜eJšžœ˜J˜J˜—š
‘œžœ
žœ
žœžœ˜7šžœ
˜˜J˜Y—˜J˜Y—˜J˜Y—˜J˜Y—J˜J˜——š
‘	œžœ
žœžœžœ
˜KJšžœžœ
žœžœ˜'šžœžœž˜J˜šžœžœž˜Jšœžœ˜Jšœžœ˜J˜J˜Jšžœ
˜—Jšžœ
˜—Jšžœ˜J˜J˜—š‘œžœžœžœžœžœ˜QJ™!J™—š‘œžœžœžœžœžœ&˜UJ™J˜—š‘
œžœ
žœ
žœ
˜;šžœ
žœ
žœ˜JšžœA˜EJšžœ0˜4—J˜J˜—š‘	œžœ
žœ‘
œžœžœ
˜MJ˜!J˜J˜J˜Jšžœ˜J˜——š
 ™š‘œžœ
žœ˜:Jšžœ˜Jšœžœ
˜šžœžœ
žœ˜Jšœžœ%˜/J˜Jš
žœžœ
žœžœ$žœ
˜QJ˜—Jšžœ˜
J˜J˜—š‘œžœ
žœ+˜KJšžœ˜J˜Jšžœžœ
žœžœ˜<Jšžœ&žœ-˜YJ˜#J˜$Jšžœ˜J˜J˜—š‘œžœ
žœ$žœ˜RJšžœ˜!Jšœžœžœ-˜AJšœžœ˜(Jš
žœžœ
žœžœžœ
˜DJ˜J˜——š
 ™Jšœžœ˜š	œžœžœ
žœ
žœ˜BJ˜—š
‘œžœ
žœ
žœ
žœ
˜4šžœžœ
žœž˜&J˜$šžœ
žœ
žœ˜Jšœžœ
˜Jšžœ	˜J˜—Jšžœ
˜—Jšžœ
žœ
˜J˜J˜—š‘
œžœ
žœ
žœ˜5šžœžœ
žœž˜&šžœžœ
žœ˜"J˜Jšžœ
˜J˜—Jšžœ
˜—J˜——š
 +™+š‘	œžœ
žœ˜Jš
œ˜Jš
œ
˜
Jšœžœ
˜J˜Jšžœ˜Jš
œ˜šžœ˜šžœ˜J˜+Jšœ+’"˜MJšœΟs
œ
₯
ž₯
œ₯
ž₯
œ₯
œ
₯
œ₯
œ₯
ž₯
œ₯
œ	₯
œ	˜JJšœ’˜+Jšœ’˜,Jšœ’˜-Jšœ’˜*Jšœ’˜(Jšœ’œ
™*šžœžœžœžœžœžœžœžœ˜OJšœ9’˜XJšœ9’˜LJ˜—J˜
—šžœ˜Jšœžœ₯
œ₯
œ₯
œ₯
œ₯
œ₯
œ₯
œ₯
œ˜N˜	Jšœ
₯
œ₯
œ₯
œ₯
œ₯
œ₯
œ₯
œ₯
œ˜DJšœ
₯
œ₯
œ₯
œ₯
œ₯
œ₯
œ₯
œ₯
œ˜E—šžœžœžœžœ˜=Jšœ9’˜IJ˜8J˜—J˜——J˜6Jš
œ˜——š
 5™5š‘œžœ
žœ/žœ
˜PJšžœ˜J˜Jš	œžœžœ
žœ
žœ˜TJšœ#’˜5Jšžœžœ
žœ˜+J˜J˜—š‘œžœ
žœ%˜JJšžœ˜J˜Jšœ(’0˜XJšœ(’2˜ZJš	œΫ
Οuœ¦
œ
¦
œ	¦
œ™ό
J˜7J˜7JšœFžœ,žœ(žœ˜¬
J˜J˜—š‘œžœ
žœ/žœ
˜NJšžœ	˜J˜Jš	œžœžœ
žœ
žœ˜TJ˜*J˜3J˜ Jšžœžœ
žœ˜+J˜—J˜—Jšžœ
˜J˜š
 
™
J˜9J˜=J˜>—š
 ™J™)š‘œžœ
žœžœžœžœ	˜XJ˜Jšœžœ˜$J˜Jšœžœ	žœ	˜2J™@J˜Jšžœ˜J˜J˜———…—PVo`