DIRECTORY
SVMasterObject, SVCastRays, SVRay, SVMasterObjectTypes, SVObjectCast, Polynomial, RealFns, SV2d, SV3d, SVRayTypes, Roots;


SVObjectCastImplC: CEDAR PROGRAM
IMPORTS SVCastRays, SVRay, Polynomial, RealFns
EXPORTS SVObjectCast =

BEGIN

Classification: TYPE = SVRayTypes.Classification;
Point2d: TYPE = SV2d.Point2d;
Point3d: TYPE = SV3d.Point3d;
Primitive: TYPE = SVRayTypes.Primitive;
Ray: TYPE = SVRayTypes.Ray;
Surface: TYPE = SVRayTypes.Surface;
TorusRec: TYPE = SVMasterObjectTypes.TorusRec;
Vector3d: TYPE = SV3d.Vector3d;

ShortRealRootRec: TYPE = Polynomial.ShortRealRootRec;
Ref: TYPE = Polynomial.Ref;

globalPoly: Polynomial.Ref;
newWay: BOOL _ TRUE;


PositiveRoots: PROCEDURE [a4, a3, a2, a1, a0: REAL] RETURNS [rootArray: Polynomial.ShortRealRootRec] = {
i: NAT _ 1;
globalPoly[4] _ a4;  -- optional
globalPoly[3] _ a3;  -- optional
globalPoly[2] _ a2;  -- optional
globalPoly[1] _ a1;  -- optional
globalPoly[0] _ a0;  -- optional
rootArray _ CheapRealRoots[globalPoly]; -- optional

IF rootArray.nRoots = 0 THEN RETURN;
WHILE i <= rootArray.nRoots DO
IF rootArray.realRoot[i-1] <= 0 THEN {
FOR j: NAT IN [i..rootArray.nRoots-1] DO
rootArray.realRoot[j-1] _ rootArray.realRoot[j];
ENDLOOP;
rootArray.nRoots _ rootArray.nRoots - 1;
}
ELSE {i _ i + 1};
ENDLOOP;
};

ToroidCast: PUBLIC PROC [localRay: Ray,  prim: Primitive, torusRec: REF ANY, surface: Surface]
RETURNS [class: Classification] = {
IF newWay THEN RETURN[NewToroidCast[localRay, prim, torusRec, surface]]
ELSE RETURN[OldToroidCast[localRay, prim, torusRec, surface]];
};

OldToroidCast: PUBLIC PROC [localRay: Ray,  prim: Primitive, torusRec: REF ANY, surface: Surface]
RETURNS [class: Classification] = {
torusData: TorusRec;
realRoots: Polynomial.ShortRealRootRec;  -- optional
B,C: REAL;
n,N,q,Q,l,L: REAL;
BB, CC: REAL;
D,E,F: REAL;
nn, NN, qq, QQ, ll, LL, NNNN, QQQQ, LLLL, nN, qQ, lL: REAL;
a4,a3,a2,a1,a0: REAL;
p: Point3d;
d: Vector3d;
class _ SVCastRays.GetClassFromPool[];
torusData _ NARROW[torusRec];
[p, d] _ SVRay.GetLocalRay[localRay];
B _ torusData.bigRadius;
C _ torusData.sectionRadius;
n _ p[1];	N _ d[1];
q _ p[2];	Q _ d[2];
l _ p[3];	L _ d[3];
BB _ B*B;
CC _ C*C;
D _ BB + CC;
E _ BB - CC;
F _ BB*CC;
nn _ n*n;  NN _ N*N;  qq _ q*q;  QQ _ Q*Q;  ll _ l*l;  LL _ L*L;
NNNN _ NN*NN;  QQQQ _ QQ*QQ;  LLLL _ LL*LL;
nN _ n*N;  qQ _ q*Q;  lL _ l*L;
a4 _ NNNN + QQQQ + LLLL + 2*(NN*(QQ + LL) + QQ*LL);
a3 _ 4*(nN*NN + qQ*QQ + lL*LL + nN*(QQ + LL) + NN*(qQ + lL) + lL*QQ + qQ*LL);
a2 _ 2*(3*(nn*NN + qq*QQ + ll*LL) + NN*(qq + ll) + nn*(QQ+LL) + 4*nN*(qQ+lL)
		- D*NN + E*QQ - D*LL + ll*QQ + LL*qq + 4*qQ*lL);
a1 _ 4*(nn*nN + qq*qQ + ll*lL + nN*(qq + ll) + (qQ + lL)*nn - nN*D + qQ*E - lL*D + qQ*ll + lL*qq);
a0 _ nn*nn + qq*qq + ll*ll + BB*BB + CC*CC +
	2*(nn*(qq + ll) - D*nn + E*qq - D*ll + qq*ll - F);

realRoots _ PositiveRoots[1.0, a3/a4, a2/a4, a1/a4, a0/a4];
SELECT realRoots.nRoots FROM
0 =>  {-- complete miss.
class.count _ 0;
class.classifs[1] _ FALSE;
};
1 => {-- a skim.  Count as a miss.
class.count _ 0;
class.classifs[1] _ FALSE;
};
2 => {-- cuts one cross section of the doughnut.  Sort roots by size.
x, y, z, r, OneMinusBoverR: REAL;
class.count _ 2;
IF realRoots.realRoot[0] < realRoots.realRoot[1] THEN
{class.params[1] _ realRoots.realRoot[0];
 class.params[2] _ realRoots.realRoot[1]}
ELSE
{class.params[2] _ realRoots.realRoot[0];
 class.params[1] _ realRoots.realRoot[1]};
class.surfaces[1] _ surface; class.surfaces[2] _ surface;
class.classifs[1] _ FALSE;  class.classifs[2] _ TRUE;  class.classifs[3] _ FALSE;
class.primitives[1] _ class.primitives[2] _ prim;
x _ n + class.params[1]*N;  y _ q + class.params[1]*Q;  z _ l + class.params[1]*L;
r _ RealFns.SqRt[x*x + z*z];
OneMinusBoverR _ 1 - (B/r);
class.normals[1] _ [x*OneMinusBoverR, y, z*OneMinusBoverR];
x _ n + class.params[2]*N;  y _ q + class.params[2]*Q;  z _ l + class.params[2]*L;
r _ RealFns.SqRt[x*x + z*z];
OneMinusBoverR _ 1 - (B/r);
class.normals[2] _ [x*OneMinusBoverR, y, z*OneMinusBoverR];
};
3 => {-- cuts one cross section.  Skims the other.  For now, count as a miss.
class.count _ 0;
class.classifs[1] _ FALSE;
};
4 => {-- cuts both cross sections.  Sort the roots.  Pair them in order.
x, y, z, r, OneMinusBoverR: REAL;
SortHitRec[realRoots];
class.count _ 4;
class.params[1] _ realRoots.realRoot[0];
class.params[2] _ realRoots.realRoot[1];
class.params[3] _ realRoots.realRoot[2];
class.params[4] _ realRoots.realRoot[3];
class.surfaces[1] _ class.surfaces[2] _ class.surfaces[3] _ class.surfaces[4] _ surface;
class.classifs[1] _ FALSE;  class.classifs[2] _ TRUE;  class.classifs[3] _ FALSE;
class.classifs[4] _ TRUE;  class.classifs[5] _ FALSE;
class.primitives[1] _ class.primitives[2] _ class.primitives[3] _ class.primitives[4] _ prim;
x _ n + class.params[1]*N;  y _ q + class.params[1]*Q;  z _ l + class.params[1]*L;
r _ RealFns.SqRt[x*x + z*z];
OneMinusBoverR _ 1 - (B/r);
class.normals[1] _ [x*OneMinusBoverR, y, z*OneMinusBoverR];
x _ n + class.params[2]*N;  y _ q + class.params[2]*Q;  z _ l + class.params[2]*L;
r _ RealFns.SqRt[x*x + z*z];
OneMinusBoverR _ 1 - (B/r);
class.normals[2] _ [x*OneMinusBoverR, y, z*OneMinusBoverR];
x _ n + class.params[3]*N;  y _ q + class.params[3]*Q;  z _ l + class.params[3]*L;
r _ RealFns.SqRt[x*x + z*z];
OneMinusBoverR _ 1 - (B/r);
class.normals[3] _ [x*OneMinusBoverR, y, z*OneMinusBoverR];
x _ n + class.params[4]*N;  y _ q + class.params[4]*Q;  z _ l + class.params[4]*L;
r _ RealFns.SqRt[x*x + z*z];
OneMinusBoverR _ 1 - (B/r);
class.normals[4] _ [x*OneMinusBoverR, y, z*OneMinusBoverR];
};
ENDCASE => ERROR;
};


SortHitRec: PROC [roots: Polynomial.ShortRealRootRec] = {
temp: REAL;
FOR i: NAT IN[0..roots.nRoots-1) DO
FOR j: NAT IN[0..roots.nRoots-i-1) DO
IF roots.realRoot[j] > roots.realRoot[j+1] THEN {
temp _ roots.realRoot[j];
roots.realRoot[j] _ roots.realRoot[j+1];
roots.realRoot[j+1] _ temp;};
ENDLOOP;
ENDLOOP;
};

EvalTorus: PRIVATE PROC [RR, RRminusSS: REAL, degree: NAT, t: REAL, localRay: Ray] RETURNS [f: REAL] = {
HH, x, y, z: REAL;
p1,p2,p3,d1,d2,d3: REAL;
basePt: Point3d;
direction: Vector3d;
[basePt, direction] _ SVRay.GetLocalRay[localRay];
p1 _ basePt[1];
p2 _ basePt[2];
p3 _ basePt[3];
d1 _ direction[1];
d2 _ direction[2];
d3 _ direction[3];
x _ p1 + t*d1;
y _ p2 + t*d2;
z _ p3 + t*d3;
SELECT degree FROM
4=>{
xxPLUSzz: REAL;
xxPLUSzz _ x*x + z*z;
HH _ xxPLUSzz + y*y;
f _ RRminusSS*RRminusSS + HH*(HH+2.0*RRminusSS) - 4.0*RR*xxPLUSzz;
};
3=>{
TwoHHdot: REAL;
HH _ x*x + y*y + z*z;
TwoHHdot _ 2.0*(x*d1 + y*d2 + z*d3);
f _ TwoHHdot*(HH + 2.0*RRminusSS) + HH*TwoHHdot - 8.0*RR*(x*d1 + z*d3);
};
2=>{
HHdot, HHdotPrime: REAL;
HH _ x*x + y*y + z*z;
HHdot _ (x*d1 + y*d2 + z*d3);
HHdotPrime _ (d1*d1+d2*d2+d3*d3);
f _ 4.0*HHdotPrime*(HH + RRminusSS) + 8.0*(HHdot*HHdot - RR*(d1*d1+d3*d3));
};
1=>{
f _ 24.0*(x*d1 + y*d2 + z*d3)*(d1*d1+d2*d2+d3*d3);
};
ENDCASE => ERROR;
};

TorusRoots: PRIVATE PROC [RR, RRminusSS: REAL, localRay: Ray, degree: NAT, bound: REAL] RETURNS [roots: Polynomial.ShortRealRootRec] =
BEGIN
roots.nRoots _ 0;
IF degree=1 THEN
BEGIN
p: Point3d;
d: Vector3d;
[p,d] _ SVRay.GetLocalRay[localRay];
roots.nRoots _ 1;
roots.realRoot[0] _ -(p[1]*d[1] + p[2]*d[2] + p[3]*d[3])/(d[1]*d[1] + d[2]*d[2] + d[3]*d[3])
END
ELSE
BEGIN
BEGIN
extrema: Polynomial.ShortRealRootRec _ TorusRoots[RR, RRminusSS, localRay, degree-1, bound];
x: REAL;
IF extrema.nRoots>0 THEN
BEGIN
x _ TRootBetween[RR, RRminusSS, degree, -bound, extrema.realRoot[0], localRay];
IF x<=extrema.realRoot[0] THEN {roots.nRoots _ 1; roots.realRoot[0] _ x};
FOR i:NAT IN [0..extrema.nRoots-1) DO
x _ TRootBetween[RR, RRminusSS, degree, extrema.realRoot[i], extrema.realRoot[i+1], localRay];
IF x<=extrema.realRoot[i+1] THEN {roots.realRoot[roots.nRoots] _ x; roots.nRoots _ roots.nRoots + 1};
ENDLOOP;
x _ TRootBetween[RR, RRminusSS, degree, extrema.realRoot[extrema.nRoots-1],bound, localRay];
IF x<=bound THEN {roots.realRoot[roots.nRoots] _ x; roots.nRoots _ roots.nRoots + 1}
END
ELSE
BEGIN
x _ TRootBetween[RR, RRminusSS, degree, -bound, bound, localRay];
IF x<=bound THEN {roots.realRoot[0] _ x; roots.nRoots _ 1}
END
END
END
END;


TRootBound: PROC [R, S: REAL, p: Point3d, d: Vector3d] RETURNS [REAL] =
BEGIN -- finds a bound for the absolute values of the roots
p1, dSup: REAL;
p1 _ ABS[p[1]]+ABS[p[2]]+ABS[p[3]];
dSup _ MAX[ABS[d[1]], MAX[ABS[d[2]], ABS[d[3]]]];
RETURN[(p1 + R + S)/dSup];
END;

TRootBetween: PROC [RR, RRminusSS: REAL, degree: NAT, x0, x1: REAL, localRay: Ray]
RETURNS [x: REAL] =
BEGIN
y0: REAL _ EvalTorus[RR, RRminusSS, degree, x0, localRay];
y1: REAL _ EvalTorus[RR, RRminusSS, degree, x1, localRay];
xx: REAL;
xx0: REAL _ IF y0<y1 THEN x0 ELSE x1;
xx1: REAL _ IF y0<y1 THEN x1 ELSE x0;
yy0: REAL _ MIN[y0,y1];
yy1: REAL _ MAX[y0,y1];
NewtonStep: PROC [x: REAL] RETURNS [newx:REAL] =
BEGIN
y: REAL _ EvalTorus[RR, RRminusSS, degree, x, localRay];
deriv: REAL _ EvalTorus[RR, RRminusSS, degree-1, x, localRay];
IF y<0 THEN {IF y>yy0 THEN {xx0 _ x; yy0 _ y}};
IF y>0 THEN {IF y<yy1 THEN {xx1 _ x; yy1 _ y}};
newx _ IF deriv=0 THEN x+y ELSE x-y/deriv;
END;
IF y0=0 THEN RETURN[x0];
IF y1=0 THEN RETURN[x1];
IF (y0 > 0.0 AND y1 > 0.0) OR (y0 < 0.0 AND y1 < 0.0) THEN RETURN[x1+100]; -- return a non-sense value which the caller will check for (there is no root in this interval).
xx _ x _ (x0+x1)/2;
WHILE x IN [x0..x1] DO
newx:REAL _ NewtonStep[x];
IF x=newx THEN {
RETURN;
};
x _ NewtonStep[newx];
xx _ NewtonStep[xx];
IF xx=x THEN EXIT;
ENDLOOP;
IF xx0 IN [x0..x1] AND xx1 IN [x0..x1] THEN
BEGIN
x0 _ MIN[xx0,xx1];
x1 _ MAX[xx0,xx1];
END;
y0 _ EvalTorus[RR, RRminusSS, degree, x0, localRay];
y1 _ EvalTorus[RR, RRminusSS, degree, x1, localRay];
THROUGH [0..500) DO
y: REAL _ EvalTorus[RR, RRminusSS, degree, x _ (x0+x1)/2, localRay];
IF x=x0 OR x=x1 THEN
{IF ABS[y0] < ABS[y1] THEN RETURN[x0] ELSE RETURN[x1]};
IF y*y0 < 0 THEN {x1 _ x; y1 _ y}
ELSE {x0 _ x; y0 _ y};
IF y0*y1>=0 THEN {IF ABS[y0] < ABS[y1] THEN RETURN[x0] ELSE RETURN[x1]}
ENDLOOP;
ERROR;
END;


PositiveTorusRoots: PRIVATE PROC [RR, RRminusSS: REAL, localRay: Ray, bound: REAL] RETURNS [roots: Polynomial.ShortRealRootRec] = {
i: NAT _ 1;
roots _ TorusRoots[RR, RRminusSS, localRay, 4, bound];
IF roots.nRoots = 0 THEN RETURN;
WHILE i <= roots.nRoots DO
IF roots.realRoot[i-1] <= 0 THEN {
FOR j: NAT IN [i..roots.nRoots-1] DO
roots.realRoot[j-1] _ roots.realRoot[j];
ENDLOOP;
roots.nRoots _ roots.nRoots - 1;
}
ELSE {i _ i + 1};
ENDLOOP;
};

InsideTorus: PRIVATE PROC [R, S: REAL, testPt: Point3d] RETURNS [BOOL] = {
r, rMinusR: REAL;
r _ RealFns.SqRt[testPt[1]*testPt[1] + testPt[3]*testPt[3]];
rMinusR _ r-R;
RETURN[rMinusR*rMinusR + testPt[2]*testPt[2] <= S*S];
};

NewToroidCast: PUBLIC PROC [localRay: Ray, prim: Primitive, torusRec: REF ANY, surface: Surface]
RETURNS [class: Classification] = {
torusData: TorusRec;
realRoots: Polynomial.ShortRealRootRec;  -- optional
B, R, S, RR, RRminusSS: REAL; -- big and little torus radii
p: Point3d;
d: Vector3d;
bound: REAL;
class _ SVCastRays.GetClassFromPool[];
torusData _ NARROW[torusRec];
R _ B _ torusData.bigRadius;
S _ torusData.sectionRadius;
RR _ R*R;
RRminusSS _ RR - S*S;
[p, d] _ SVRay.GetLocalRay[localRay];

bound _ TRootBound[R, S, p, d];
realRoots _ PositiveTorusRoots[RR, RRminusSS, localRay, bound];

SELECT realRoots.nRoots FROM
0 =>  {-- complete miss.
class.count _ 0;
class.classifs[1] _ FALSE;
};
1 => {-- a skim.  Count as a miss.
IF InsideTorus[R, S, p] THEN {
x, y, z, r, OneMinusBoverR: REAL;
class.count _ 1;
class.classifs[1] _ TRUE; class.classifs[2] _ FALSE;
class.params[1] _ realRoots.realRoot[0];
[x, y, z] _ SVRay.EvaluateLocalRay2[localRay, class.params[1]];
r _ RealFns.SqRt[x*x + z*z];
OneMinusBoverR _ 1 - (B/r);
class.normals[1] _ [x*OneMinusBoverR, y, z*OneMinusBoverR];
class.surfaces[1] _ surface;
class.primitives[1] _ prim;
}
ELSE SVCastRays.MakeClassAMiss[class];
};
2 => {-- cuts one cross section of the doughnut.  Sort roots by size.
x, y, z, r, OneMinusBoverR: REAL;
class.count _ 2;
IF realRoots.realRoot[0] < realRoots.realRoot[1] THEN
{class.params[1] _ realRoots.realRoot[0];
 class.params[2] _ realRoots.realRoot[1]}
ELSE
{class.params[2] _ realRoots.realRoot[0];
 class.params[1] _ realRoots.realRoot[1]};
class.surfaces[1] _ surface; class.surfaces[2] _ surface;
class.classifs[1] _ FALSE;  class.classifs[2] _ TRUE;  class.classifs[3] _ FALSE;
class.primitives[1] _ class.primitives[2] _ prim;
[x, y, z] _ SVRay.EvaluateLocalRay2[localRay, class.params[1]];
r _ RealFns.SqRt[x*x + z*z];
OneMinusBoverR _ 1 - (B/r);
class.normals[1] _ [x*OneMinusBoverR, y, z*OneMinusBoverR];
[x, y, z] _ SVRay.EvaluateLocalRay2[localRay, class.params[2]];
r _ RealFns.SqRt[x*x + z*z];
OneMinusBoverR _ 1 - (B/r);
class.normals[2] _ [x*OneMinusBoverR, y, z*OneMinusBoverR];
};
3 => {-- cuts one cross section.  Skims the other (or starts inside of the torus).  For now, count as a miss.
class.count _ 0;
class.classifs[1] _ FALSE;
};
4 => {-- cuts both cross sections.  Sort the roots.  Pair them in order.
x, y, z, r, OneMinusBoverR: REAL;
SortHitRec[realRoots];
class.count _ 4;
class.params[1] _ realRoots.realRoot[0];
class.params[2] _ realRoots.realRoot[1];
class.params[3] _ realRoots.realRoot[2];
class.params[4] _ realRoots.realRoot[3];
class.surfaces[1] _ class.surfaces[2] _ class.surfaces[3] _ class.surfaces[4] _ surface;
class.classifs[1] _ FALSE;  class.classifs[2] _ TRUE;  class.classifs[3] _ FALSE;
class.classifs[4] _ TRUE;  class.classifs[5] _ FALSE;
class.primitives[1] _ class.primitives[2] _ class.primitives[3] _ class.primitives[4] _ prim;
[x, y, z] _ SVRay.EvaluateLocalRay2[localRay, class.params[1]];
r _ RealFns.SqRt[x*x + z*z];
OneMinusBoverR _ 1 - (B/r);
class.normals[1] _ [x*OneMinusBoverR, y, z*OneMinusBoverR];
[x, y, z] _ SVRay.EvaluateLocalRay2[localRay, class.params[2]];
r _ RealFns.SqRt[x*x + z*z];
OneMinusBoverR _ 1 - (B/r);
class.normals[2] _ [x*OneMinusBoverR, y, z*OneMinusBoverR];
[x, y, z] _ SVRay.EvaluateLocalRay2[localRay, class.params[3]];
r _ RealFns.SqRt[x*x + z*z];
OneMinusBoverR _ 1 - (B/r);
class.normals[3] _ [x*OneMinusBoverR, y, z*OneMinusBoverR];
[x, y, z] _ SVRay.EvaluateLocalRay2[localRay, class.params[4]];
r _ RealFns.SqRt[x*x + z*z];
OneMinusBoverR _ 1 - (B/r);
class.normals[4] _ [x*OneMinusBoverR, y, z*OneMinusBoverR];
};
ENDCASE => ERROR;
};

RootBound: PROC [poly: Ref] RETURNS [REAL] =
BEGIN -- finds a bound for the absolute values of the roots
d: NAT _ Degree[poly];
b: REAL _ 0;
FOR i:NAT IN [0..d) DO b _ b + ABS[poly[i]] ENDLOOP;
RETURN [b*1.01/ABS[poly[d]]+1];
END;

RootBetween: PROC [poly: Ref, x0, x1: REAL]
RETURNS [x: REAL] =
BEGIN
y0: REAL _ Eval[poly,x0];
y1: REAL _ Eval[poly,x1];
xx: REAL;
xx0: REAL _ IF y0<y1 THEN x0 ELSE x1;
xx1: REAL _ IF y0<y1 THEN x1 ELSE x0;
yy0: REAL _ MIN[y0,y1];
yy1: REAL _ MAX[y0,y1];
NewtonStep: PROC [x: REAL] RETURNS [newx:REAL] =
BEGIN
y: REAL _ Eval[poly, x];
deriv: REAL _ EvalDeriv[poly, x];
IF y<0 THEN {IF y>yy0 THEN {xx0 _ x; yy0 _ y}};
IF y>0 THEN {IF y<yy1 THEN {xx1 _ x; yy1 _ y}};
newx _ IF deriv=0 THEN x+y ELSE x-y/deriv;
END;
IF y0=0 THEN RETURN[x0];
IF y1=0 THEN RETURN[x1];
IF (y0 > 0.0 AND y1 > 0.0) OR (y0 < 0.0 AND y1 < 0.0) THEN RETURN[x1+100]; -- return a non-sense value which the caller will check for (there is no root in this interval).
xx _ x _ (x0+x1)/2;
WHILE x IN [x0..x1] DO
newx:REAL _ NewtonStep[x];
IF x=newx THEN {
RETURN;
};
x _ NewtonStep[newx];
xx _ NewtonStep[xx];
IF xx=x THEN EXIT;
ENDLOOP;
IF xx0 IN [x0..x1] AND xx1 IN [x0..x1] THEN
BEGIN
x0 _ MIN[xx0,xx1];
x1 _ MAX[xx0,xx1];
END;
y0 _ Eval[poly,x0];
y1 _ Eval[poly,x1];
THROUGH [0..500) DO
y: REAL _ Eval[poly, x _ (x0+x1)/2];
IF x=x0 OR x=x1 THEN
{IF ABS[y0] < ABS[y1] THEN RETURN[x0] ELSE RETURN[x1]};
IF y*y0 < 0 THEN {x1 _ x; y1 _ y}
ELSE {x0 _ x; y0 _ y};
IF y0*y1>=0 THEN {IF ABS[y0] < ABS[y1] THEN RETURN[x0] ELSE RETURN[x1]}
ENDLOOP;
ERROR;
END;

Degree: PUBLIC PROC [poly: Ref] RETURNS [NAT] =
BEGIN
FOR i:NAT DECREASING IN [0..poly.maxDegreePlusOne) DO
IF poly[i] # 0 THEN RETURN[i];
ENDLOOP;
RETURN[0];
END;

Differentiate: PUBLIC PROC [poly: Ref] =
BEGIN
d:NAT _ Degree[poly];
FOR i: NAT IN (0..d] DO poly[i-1] _ i*poly[i] ENDLOOP;
poly[d] _ 0.0;
END;

CheapRealRoots: PUBLIC PROC [poly: Ref] RETURNS [roots: ShortRealRootRec] =
BEGIN
d: NAT _ Degree[poly];
roots.nRoots _ 0;
IF d<=1 THEN
BEGIN
IF d=1 THEN {roots.nRoots _ 1; roots.realRoot[0] _ -poly[0]/poly[1]}
END
ELSE
BEGIN
bound: REAL _ RootBound[poly];
savedCoeff: ARRAY [0..5] OF REAL;
FOR i: NAT IN [0..d] DO savedCoeff[i] _ poly[i] ENDLOOP;
Differentiate[poly];
BEGIN
extrema: ShortRealRootRec _ CheapRealRoots[poly];
x: REAL;
FOR i: NAT IN [0..d] DO poly[i] _ savedCoeff[i] ENDLOOP;
IF extrema.nRoots>0 THEN
BEGIN
x _ RootBetween[poly,-bound,extrema.realRoot[0]];
IF x<=extrema.realRoot[0] THEN {roots.nRoots _ 1; roots.realRoot[0] _ x};
FOR i:NAT IN [0..extrema.nRoots-1) DO
x _ RootBetween[poly,extrema.realRoot[i],extrema.realRoot[i+1]];
IF x<=extrema.realRoot[i+1] THEN {roots.realRoot[roots.nRoots] _ x; roots.nRoots _ roots.nRoots + 1};
ENDLOOP;
x _ RootBetween[poly,extrema.realRoot[extrema.nRoots-1],bound];
IF x<=bound THEN {roots.realRoot[roots.nRoots] _ x; roots.nRoots _ roots.nRoots + 1}
END
ELSE
BEGIN
x _ RootBetween[poly,-bound,bound];
IF x<=bound THEN {roots.realRoot[0] _ x; roots.nRoots _ 1}
END
END
END
END;

EvalDeriv: PUBLIC PROC [f: Ref, x: REAL] RETURNS [y: REAL] =
BEGIN
y _ 0.0;
FOR i:NAT DECREASING IN (0..Degree[f]] DO
y _ i*f[i] + y*x
ENDLOOP
END;

Eval: PUBLIC PROC [f: Ref, x: REAL] RETURNS [y: REAL] =
BEGIN
y _ 0.0;
FOR i:NAT DECREASING IN [0..Degree[f]] DO
y _ y*x + f[i]
ENDLOOP
END;


Init: PROC = {
globalPoly _ Polynomial.Quartic[[0,0,0,0,0]];
};

Init[];

END.
File: SVObjectCastImplC.mesa
Last edited by Bier on January 2, 1983 8:31 pm
Copyright c 1984 by Xerox Corporation.  All rights reserved.
Author: Bier on March 29, 1987 5:06:28 pm PST

Debugging
A torus has a single surface.  Classifying a ray with respect to this surface involves solving a single quadric equation.
Let the radius from the center of the doughnut hole to the center of the doughnut be C.
Let the cross-sectional radius be B.
Let the doughnut be a circle revolved around the y axis.
Consider a cross section of the doughnut.  Call the line parallel to the y axis and running through the center of the cross section L.  Let r(y) be the distance perpendicular to L from L to the circular boundary of the cross section.  Then:
(1)  r^2 + y^2 = C^2
From which we find:
(2)  r = +- SqRt[C^2 - y^2]
Thus the two concentric circles cut by the plane y = y0 have equations:
x^2 + z^2  = (B +- SqRt[C^2 - y^2])^2 or compactly:
(3) x^2 + z^2  = (B +- r)^2 = B^2 +- 2Br + r^2
We can get rid of the square root by regrouping:
(4) x^2 + z^2 - B^2 - r^2 = +- 2Br
and squaring:
(5) (x^2 + z^2 - B^2 - r^2)^2 = 4B^2r^2
Expanding, we have:
x^4 + x^2z^2 - x^2*B^2 - x^2*r^2 +
z^2*x^2 + z^4 - z^2*B^2 - z^2*r^2
-B^2*x^2 -B^2*z^2 + B^4 + B^2*r^2
-r^2*x^2 - r^2*z^2 + B^2*r^2 + r^4	=	4*B^2*r^2;		OR

x^4  + z^4 + B^4 + r^4 +
2*x^2z^2 - 2*x^2*B^2 - 2*x^2*r^2 - 2*z^2*B^2 - 2*z^2*r^2 + 2*B^2*r^2
	= 4*B^2*r^2
	
Plugging in r^2 = (C^2 - y^2), r^4 = (C^4 - 2*C^2*y^2 + y^4):

x^4  + z^4 + B^4 + C^4 - 2*C^2*y^2 + y^4 +
2*x^2z^2 - 2*x^2*B^2 - 2*x^2*(C^2 - y^2) - 2*z^2*B^2 - 2*z^2*(C^2 - y^2) + 2*B^2*(C^2 - y^2)
	= 4*B^2*(C^2 - y^2) = 4*B^2*C^2 - 4*B^2*y^2	OR
	
 x^4  + y^4 + z^4 + B^4 + C^4 +
2*x^2z^2 + 2*x^2*y^2 + 2*z^2*y^2
x^2 (- 2*B^2 - 2*C^2) +
y^2 (- 2*B^2- 2*C^2) +
z^2 (- 2*B^2 - 2*C^2) +
- 2*C^2 + 2*B^2*C^2	=	4*B^2*C^2 - 4*B^2*y^2	OR

x^4  + y^4 + z^4 + B^4 + C^4 +
2*x^2z^2 + 2*x^2*y^2 + 2*z^2*y^2
x^2 (- 2*B^2 - 2*C^2) +
y^2 (2*B^2 - 2*C^2) +
z^2 (- 2*B^2 - 2*C^2) +
- 2*B^2*C^2	= 0

Let D = B^2 + C^2,  E = B^2 - C^2, F = C^2B^2  Then:

(6)  x^4  + y^4 + z^4 + B^4 + C^4 +
		+ 2(x^2*y^2 + x^2z^2 + y^2z^2
		  - Dx^2 + Ey^2 - Dz^2 - F) = 0

Now recall the ray equations (using single letters for all variables for compactness):
(7)  x = n + t*N
(8)  y = q + t*Q
(9)  z = l + t*L
We wish to solve for t by plugging (7) - (9) into (6).  We get:

(n + t*N)^4 + (q + t*Q)^4 + (l + t*L)^4 + B^4 + C^4
+ 2( (n + t*N)^2*(q + t*Q)^2 + (n + t*N)^2*(l + t*L)^2 + (q + t*Q)^2*(l + t*L)^2
	-D(n + t*N)^2 + E(q + t*Q)^2 - D(l + t*L)^2 - F) = 0
	
If I've done my algebra correctly, this boils down to:
a4*t^4 + a3*t^3 + a2*t^2 + a1*t + a0
Where:
a4 = N^4 + Q^4 + L^4 + 2(N^2(Q^2+L^2) + Q^2L^2);
a3 = 4(nN^3 + qQ^3 + lL^3 + nN(Q^2 + L^2) + N^2(qQ + lL) + lLQ^2 + qQL^2)
a2 = 2(3n^2N^2 + 3q^2Q^2 + 3l^2L^2 + N^2(q^2 + l^2) + n^2(Q^2+L^2) + 4nN(qQ + lL)
- DN^2 + EQ^2 - DL^2 + l^2Q^2 + L^2q^2 + 4qQlL)
a1 = 4(n^3N + q^3Q + l^3L + nN(q^2 + l^2) + (qQ + lL)n^2 - nND + qQE - lLD + qQl^2 + lLq^2)
a0 = n^4 + q^4 + l^4 + B^4 + C^4 +
2(n^2(q^2 + l^2) - Dn^2 + Eq^2 - Dl^2+ q^2l^2  - F)
Here goes:

Use only the positive roots.
[realRootsOneOrigin, rootCount] _ Roots.RealQuartic[a4, a3, a2, a1, a0];
realRoots[0] _ realRootsOneOrigin[1]; realRoots[1] _ realRootsOneOrigin[2];
realRoots[2] _ realRootsOneOrigin[3]; realRoots[3] _ realRootsOneOrigin[4];

realRoots: ARRAY[0..3] OF REAL;
realRootsOneOrigin: ARRAY[1..4] OF REAL;
rootCount: NAT;

SELECT rootCount FROM
Calculate the surface normals at the two hits points.
To normalize, divide by C.  (I won't do this since the lighting model normalizes anyway)
SortHitRec[realRoots, 4];
Calculate the surface normals at the four hits points.
To normalize, divide by C.  (I won't do this since the lighting model normalizes anyway)
A simple bubble sort. Sorts in increasing order (depth).
SortHitRec: PROC [roots: ARRAY[0..3] OF REAL, rootCount: NAT] = {
temp: REAL;
FOR i: NAT IN[0..rootCount-1) DO
FOR j: NAT IN[0..rootCount-i-1) DO
IF roots[j] > roots[j+1] THEN {
temp _ roots[j];
roots[j] _ roots[j+1];
roots[j+1] _ temp;};
ENDLOOP;
ENDLOOP;
};
	
debug Feedback.Append[IO.PutFR["%d degree, %d roots: ", IO.int[degree - 1], IO.int[extrema.nRoots]], TRUE, FALSE];
debug FOR i: NAT IN [0..extrema.nRoots) DO
debug Feedback.Append[IO.PutFR["%g ", IO.real[extrema.realRoot[i]]], TRUE, FALSE];
debug ENDLOOP;
debug Feedback.AppendTypescript[""];
We do this for the torus by reasoning geometrically.  In local coordinates, the torus is
centered on [0,0,0].  The ray base-point is at point p.  Hence (|p| + R + S)/|d| will bound
the value of d.  We can use the 1-norm to overestimate |p| and the sup-norm to underestimate |d| to get a conservative bound on t.
debug  xi: ARRAY [0..50) OF REAL;
debug  i: NAT _ 0;
IF y0*y1 > 0 THEN RETURN[x1+100]; -- return a non-sense value which the caller will check for (there is no root in this interval).
first try Newton's method
debug  Feedback.Append[IO.PutFR["TNewton: %d steps: ", IO.int[i]], TRUE,FALSE];
debug  FOR j: NAT IN [0..i) DO
debug  Feedback.Append[IO.PutFR["%g ", IO.real[xi[j]]], TRUE, FALSE];
debug  ENDLOOP;
debug  Feedback.AppendTypescript[""];
debug  xi[i] _ xx; i _ i+1;
If it oscillated or went out of range, try bisection
debug  SIGNAL Bisection;
debug  Feedback.Append[IO.PutFR["%d degree, %d roots: ", IO.int[4], IO.int[roots.nRoots]], TRUE, FALSE];
debug  FOR i: NAT IN [0..roots.nRoots) DO
debug  Feedback.Append[IO.PutFR["%g ", IO.real[roots.realRoot[i]]], TRUE, FALSE];
debug  ENDLOOP;
debug  Feedback.AppendTypescript[""];
SELECT rootCount FROM
Calculate the surface normals at the two hits points.
Calculate the surface normals at the two hits points.
To normalize, divide by C.  (I won't do this since the lighting model normalizes anyway)
SortHitRec[realRoots, 4];
Calculate the surface normals at the four hits points.
To normalize, divide by C.  (I won't do this since the lighting model normalizes anyway)
debug  xi: ARRAY [0..50) OF REAL;
debug  i: NAT _ 0;
IF y0*y1 > 0 THEN RETURN[x1+100];
first try Newton's method
debug  Feedback.Append[IO.PutFR["Newton: %d steps: ", IO.int[i]], TRUE,FALSE];
debug  FOR j: NAT IN [0..i) DO
debug  Feedback.Append[IO.PutFR["%g ", IO.real[xi[j]]], TRUE, FALSE];
debug  ENDLOOP;
debug  Feedback.AppendTypescript[""];
debug  xi[i] _ xx; i _ i+1;
If it oscillated or went out of range, try bisection
debug  SIGNAL Bisection;
debug  Feedback.Append[IO.PutFR["%g bound, %d degree %d roots: ", IO.real[bound], IO.int[d-1], IO.int[extrema.nRoots]], TRUE, FALSE];
debug  FOR i: NAT IN [0..extrema.nRoots) DO
debug  Feedback.Append[IO.PutFR["%g ", IO.real[extrema.realRoot[i]]], TRUE, FALSE];
debug  ENDLOOP;
debug  Feedback.AppendTypescript[""];
Κ#–

"Cedar" style˜Ihead1š
œ™Jš
œ.™.Jšœ
Οm
œ1™<Jš
œ.™.š
Οk	˜	JšœžœY˜z—J˜šœž
˜ Jšžœ
žœ˜.Jšžœ˜—šžœ
˜Jšœžœ˜1Jšœ	žœ˜Jšœ	žœ˜Jšœžœ˜'Jšœžœ˜Jšœ	žœ˜#Jšœ
žœ ˜.Jšœ
žœ˜˜J™	—Jšœžœ˜5Jšœžœ˜J˜Jš
œ˜Jšœžœžœ
˜—J˜Jšœ€Οb4œŸeœΑŸ$œΓ™ΓJ˜šΟn
œž	œžœžœ-˜hJ™Iprocšœžœ˜LšœΟc˜ Jšœ‘˜ Jšœ‘˜ Jšœ‘˜ Jšœ‘˜ Lšœ(‘˜4Lš
œH™HJš
œK™KJš
œL™LLšžœžœ
žœ
˜$šžœž˜šžœžœ˜&šžœžœ
žœž˜(Lš
œ0˜0—Lšžœ
˜Lš
œ(˜(L˜
—Lšžœ
˜—Lšžœ
˜Jš
œ˜—J˜š
 
œžœ
žœ-žœ
žœ˜^šžœ˜#Jšžœžœ
žœ2˜GJšžœ
žœ3˜>J˜J˜——š
 
œžœ
žœ-žœ
žœ˜aJšžœ˜#J˜Jšœ)‘˜4J™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šžœžœ
žœžœžœ
žœžœžœ
žœ
˜+Jšœž
œ
ž
œ
ž
œ
˜Jšœžœžœžœžœžœžœžœ
žœ˜3Jšœžœžœžœžœžœžœžœžœ˜MJšœžœžœžœžœžœ
žœž
œ
žœž
œ
žœž
œ
žœžœžœ˜JšœAž
œž
œž
œ˜bJšœžœ
žœžœ
žœž
œž
œž
œž
œ˜`J˜Jš
œ;˜;J™Jšžœž˜J™šœ‘˜Jš
œ˜Jšœžœ
˜Jš
œ˜—šœ‘˜"Jš
œ˜Jšœžœ
˜Jš
œ˜—šœ‘?˜EJšœžœ
˜!J˜šžœ/ž˜5Jš
œ)˜)Jš
œ)˜)—š
ž˜Jš
œ)˜)Jš
œ*˜*—Jš
œ9˜9Jšœžœžœžœ
˜QJš
œ1˜1J™5Jšœž
œž
œž
œ
˜RJ˜Jš
œ˜Jš
œ;˜;J™XJšœž
œž
œž
œ
˜RJ˜Jš
œ˜Jš
œ;˜;Jš
œ˜—šœ‘G˜MJš
œ˜Jšœžœ
˜Jš
œ˜—šœ‘B˜HJšœžœ
˜!Jš
œ™Jš
œ˜J˜Jš
œ(˜(Jš
œ(˜(Jš
œ(˜(Jš
œ(˜(Jš
œX˜XJšœžœžœžœ
˜QJšœžœžœ
˜5Jš
œ]˜]J™6Jšœž
œž
œž
œ
˜RJ˜Jš
œ˜Jš
œ;˜;J™XJšœž
œž
œž
œ
˜RJ˜Jš
œ˜Jš
œ;˜;Jšœž
œž
œž
œ
˜RJ˜Jš
œ˜Jš
œ;˜;Jšœž
œž
œž
œ
˜RJ˜Jš
œ˜Jš
œ;˜;J˜—Jšžœžœ
˜Jš
œ˜—J˜š 
œžœ)˜9Jš
œ8™8Jšœžœ
˜šžœžœ
žœž˜#šžœžœ
žœž˜%šžœ)žœ˜1Jš
œ˜Jš
œ(˜(Jš
œ˜——Jšžœ
˜—Jšžœ
˜Jš
œ˜—š 
œžœ1™AJšœžœ
™šžœžœ
žœž™ šžœžœ
žœž™"šžœžœ™Jš
œ™Jš
œ™Jš
œ™——Jšžœ
™—Jšžœ
™Jš
œ™—š 	œžœ
žœžœ
žœ
žœžœžœžœ˜hJšžœžœ
˜J˜J˜J˜Jšœ œ˜2Lš
œ˜Lš
œ˜Lš
œ˜Lš
œ˜Lš
œ˜Lš
œ˜Lš
œ˜L˜L˜Jšžœž˜š
œ˜Jšœ
žœ
˜J˜Jšžœ˜Jšœžœžœžœ
˜BJš
œ˜—š
œ˜Jšœ
žœ
˜Jšžœ˜J˜$Jšœžœžœžœ˜GJ˜—š
œ˜Jšœžœ
˜Jšžœ˜J˜Jš
œ!˜!Jšœžœ#žœ˜KJš
œ˜—š
œ˜Jš
œ2˜2J˜—Jšžœžœ
˜J˜J˜—š 
œžœ
žœžœ
žœžœ	žœžœ'˜†
Jš
ž˜J˜šžœ
ž˜š
ž˜J˜Jš
œ˜Jšœžœ˜$Jš
œ˜Jš
œ\˜\—Jš
ž˜—š
ž˜š
ž˜Jš
ž˜Jšœ2žœ(˜\Jšœžœ
˜Jšœžœ žœžœžœžœ™ršœžœžœ
žœž™*Jš	œžœžœžœžœ™R—Jšœžœ
™Jš
œ$™$šžœž˜Jš
ž˜Jšœžœ<˜OJšžœžœ+˜Išžœžœ
žœž˜%JšœžœK˜^JšžœžœE˜eJšžœ
˜—JšœžœI˜\Jšžœ
žœD˜TJš
ž˜—š
ž˜Jš
ž˜Jšœžœ.˜AJšžœ
žœ*˜:Jš
ž˜—Jš
ž˜—Jš
ž˜—Jšžœ
˜J˜—J˜š 
œžœž
œž
œžœžœžœ˜Gšžœ
‘5˜;J™·Jšœ
žœ
˜J˜#Jšœžœžœ˜1Jšžœž
œž
œ˜—Jšžœ
˜J˜—š œžœžœ
žœ
žœ
žœ˜RJšžœžœ˜Jš
ž˜Jšœžœ
žœ#˜:Jšœžœ
žœ#˜:Jš
œ!™!Jš
œ™Jšœžœ
˜	Jš	œžœžœžœžœ˜%Jš	œžœžœžœžœ˜%Jšœžœžœ˜Jšœžœžœ˜š
 
œžœžœžœžœ˜0Jš
ž˜Jšœžœ
žœ"˜8Jšœžœ
žœ$˜>Jšžœžœžœžœ˜/Jšžœžœžœžœ˜/Jšœžœ	žœžœ˜*Jšžœ
˜—Jšžœžœ
žœ˜Jšžœžœ
žœ˜Jš
žœžœžœžœžœ
žœ
‘`˜«
Jšžœžœ
žœ
‘`™‚
J˜Jš
œ™šžœžœ
ž˜Jšœžœ˜šžœž˜Jš	œžœžœ
žœ
žœ™Ošœžœžœ
žœž™Jš	œžœžœžœžœ™E—Jšœžœ
™J™%Jšžœ
˜J˜—J˜J˜Jš
œ™Jšžœžœ
žœ
˜Jšžœ
˜—Jš
œ4™4Jš
œ™š	žœžœ
žœžœ
ž˜+Jš
ž˜Jšœžœ
˜Jšœžœ
˜Jšžœ
˜—Jšœžœ#˜4Jšœžœ#˜4šžœ
ž˜Jšœžœ
žœ.˜Dšžœžœž˜Jšœ
žœ
žœžœžœ
žœžœ
žœ˜7—Jšžœ
žœ˜!Jšžœ˜Jšžœ
žœžœ
žœžœžœ
žœžœ
žœ˜GJšžœ
˜—Jšžœ
˜Jšžœ
˜J˜J˜—š œžœ
žœžœ
žœžœžœ)˜ƒ
Jšœžœ˜Jšœžœ!˜6Jšœžœ žœ	žœžœžœ™hšœžœžœ
žœž™)Jš	œžœžœžœžœ™Q—Jšœžœ
™Jš
œ%™%Jšžœžœ
žœ
˜ šžœž˜šžœžœ˜"šžœžœ
žœž˜$Lš
œ(˜(—Lšžœ
˜Lš
œ ˜ L˜
—Lšžœ
˜—Lšžœ
˜L˜L˜—š œžœ
žœž
œž
œžœžœžœ˜JJšœžœ
˜J˜<Jšœž
œ
˜Jšžœ*ž
œ
ž
œ˜5J˜J˜—š
 
œžœ
žœ,žœ
žœ˜`Jšžœ˜#J˜Jšœ)‘˜4Jšž
œž
œž
œžœ
žœ‘˜;J˜J˜Jšœžœ
˜Jš
œ&˜&Jšœžœ˜Jšž
œž
œ˜Jšž
œ˜Jšžœž
œ
ž
œ
˜	Jšœžœž
œ
ž
œ
˜Jšœ	žœ˜%J˜Jš
œ˜Jšœžœ˜?J˜Jšžœž˜J™šœ‘˜Jš
œ˜Jšœžœ
˜Jš
œ˜—šœ‘˜"šžœ
ž
œž
œžœ˜Jšœžœ
˜!Lš
œ˜Lšœžœžœ
˜4Lš
œ(˜(J™5Jšœžœ.˜?J˜Jšœž
œ˜Jš
œ;˜;Lš
œ˜Lš
œ˜L˜
—Lšžœ"˜&Jš
œ˜—šœ‘?˜EJšœžœ
˜!J˜šžœ/ž˜5Jš
œ)˜)Jš
œ)˜)—š
ž˜Jš
œ)˜)Jš
œ*˜*—Jš
œ9˜9Jšœžœžœžœ
˜QJš
œ1˜1J™5Jšœžœ.˜?J˜Jšœž
œ˜Jš
œ;˜;J™XJšœžœ.˜?J˜Jšœž
œ˜Jš
œ;˜;Jš
œ˜—šœ‘g˜mJš
œ˜Jšœžœ
˜Jš
œ˜—šœ‘B˜HJšœžœ
˜!Jš
œ™Jš
œ˜J˜Jš
œ(˜(Jš
œ(˜(Jš
œ(˜(Jš
œ(˜(Jš
œX˜XJšœžœžœžœ
˜QJšœžœžœ
˜5Jš
œ]˜]J™6Jšœžœ.˜?J˜Jšœž
œ˜Jš
œ;˜;J™XJšœžœ.˜?J˜Jšœž
œ˜Jš
œ;˜;Jšœžœ.˜?J˜Jšœž
œ˜Jš
œ;˜;Jšœžœ.˜?J˜Jšœž
œ˜Jš
œ;˜;J˜—Jšžœžœ
˜Jš
œ˜J˜—š 	œžœ
žœžœ˜,Jšžœ
‘5˜;Jšœžœ˜Jšœžœ˜Jšžœžœ
žœžœ	žœ
žœ
˜4Jšžœ	žœ
˜Jšžœ
˜J˜—š œžœžœ
˜+Jšžœžœ˜Jš
ž˜Jšœžœ	žœ
žœ
™!Jšœ
žœ™Jšœžœ˜Jšœžœ˜Jšœžœ
˜	Jš	œžœžœžœžœ˜%Jš	œžœžœžœžœ˜%Jšœžœžœ˜Jšœžœžœ˜š
 
œžœžœžœžœ˜0Jš
ž˜Jšœžœ˜Jšœžœ˜!Jšžœžœžœžœ˜/Jšžœžœžœžœ˜/Jšœžœ	žœžœ˜*Jšžœ
˜—Jšžœžœ
žœ˜Jšžœžœ
žœ˜Jš
žœžœžœžœžœ
žœ
‘`˜«
Jšžœžœ
žœ	™!J˜Jš
œ™šžœžœ
ž˜Jšœžœ˜šžœžœ
ž
˜Jš	œžœžœ
žœ
žœ™Nšœžœžœ
žœž™Jš	œžœžœžœžœ™E—Jšœžœ
™J™%Jšžœ
˜J˜—J˜J˜Jš
œ™Jšžœžœ
žœ
˜Jšžœ
˜—Jš
œ4™4Jš
œ™š	žœžœ
žœžœ
ž˜+Jš
ž˜Jšœžœ
˜Jšœžœ
˜Jšžœ
˜—J˜J˜šžœ
ž˜Jšœžœ˜$šžœžœž˜Jšœ
žœ
žœžœžœ
žœžœ
žœ˜7—Jšžœ
žœ˜!Jšžœ˜Jšžœ
žœžœ
žœžœžœ
žœžœ
žœ˜GJšžœ
˜—Jšžœ
˜Jšžœ
˜J˜—š
 œžœ
žœ
žœžœ˜/Jš
ž˜š	žœžœ
ž
œ
žœž˜5Jšžœ
žœ
žœ˜Jšžœ
˜—Jšžœ˜
Jšžœ
˜J˜—š 
œžœ
žœ˜(Jš
ž˜Jšœžœ˜Jš
žœžœ
žœžœžœ
˜6J˜Jšžœ
˜J˜—š œžœ
žœ
žœ˜KJš
ž˜Jšœžœ˜J˜šžœž˜Jš
ž˜Jšžœžœ9˜DJš
ž˜—š
ž˜Jš
ž˜Jšœžœ˜Jšœžœžœ
žœ
˜!Jš
žœžœ
žœžœžœ
˜8š
œ˜Jš
ž˜Jš
œ1˜1Jšœžœ
˜Jš
œžœ)žœžœžœžœžœ™…
šœžœžœ
žœž™+Jš	œžœžœžœžœ™S—Jšœžœ
™Jš
œ%™%Jš
žœžœ
žœžœžœ
˜8šžœž˜Jš
ž˜Jš
œ1˜1Jšžœžœ+˜Išžœžœ
žœž˜%Jš
œ@˜@JšžœžœE˜eJšžœ
˜—Jš
œ?˜?Jšžœ
žœD˜TJš
ž˜—š
ž˜Jš
ž˜J˜#Jšžœ
žœ*˜:Jš
ž˜—Jš
ž˜—Jš
ž˜—Jšžœ
˜J˜—š 	œžœ
žœ
žœžœžœ˜<Jš
ž˜J˜š	žœžœ
ž
œ
žœž˜)J˜Jš
ž˜—Jšžœ
˜J˜—š œžœ
žœ
žœžœžœ˜7Jš
ž˜J˜š	žœžœ
ž
œ
žœž˜)J˜Jš
ž˜—Jšžœ
˜J˜—J˜š œžœ˜Jš
œ-˜-Jš
œ˜—J˜Jš
œ˜Jšžœ
˜—…—Dδz