<> <> <> DIRECTORY BasicObject3d, CastRays, CSG, ObjectCast, Polynomial, RealFns, SV2d, Roots; ObjectCastImplC: PROGRAM IMPORTS CastRays, Polynomial, RealFns EXPORTS ObjectCast = BEGIN <> Classification: TYPE = CSG.Classification; Point2d: TYPE = SV2d.Point2d; Primitive: TYPE = CSG.Primitive; Ray: TYPE = CSG.Ray; Surface: TYPE = CSG.Surface; TorusRec: TYPE = BasicObject3d.TorusRec; globalPoly: Polynomial.Ref; 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 _ Polynomial.CheapRealRoots[globalPoly]; -- optional <<[realRootsOneOrigin, rootCount] _ Roots.RealQuartic[a4, a3, a2, a1, a0];>> <> <> 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] = { 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; class _ CastRays.GetClassFromPool[]; torusData _ NARROW[torusRec]; B _ torusData.bigRadius; C _ torusData.sectionRadius; n _ localRay.basePt[1]; N _ localRay.direction[1]; q _ localRay.basePt[2]; Q _ localRay.direction[2]; l _ localRay.basePt[3]; L _ localRay.direction[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[a4, a3, a2, a1, a0]; <<>> 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, theta: 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; theta _ RealFns.ArcTan[-z,x]; class.normals[1] _ [-B*RealFns.Cos[theta] + x, y, B*RealFns.Sin[theta] + z]; <> x _ n + class.params[2]*N; y _ q + class.params[2]*Q; z _ l + class.params[2]*L; theta _ RealFns.ArcTan[-z,x]; class.normals[2] _ [-B*RealFns.Cos[theta] + x, y, B*RealFns.Sin[theta] + z]; }; 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, theta: 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[1] _ TRUE; class.classifs[2] _ 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; theta _ RealFns.ArcTan[-z,x]; class.normals[1] _ [-B*RealFns.Cos[theta] + x, y, B*RealFns.Sin[theta] + z]; <> x _ n + class.params[2]*N; y _ q + class.params[2]*Q; z _ l + class.params[2]*L; theta _ RealFns.ArcTan[-z,x]; class.normals[2] _ [-B*RealFns.Cos[theta] + x, y, B*RealFns.Sin[theta] + z]; x _ n + class.params[3]*N; y _ q + class.params[3]*Q; z _ l + class.params[3]*L; theta _ RealFns.ArcTan[-z,x]; class.normals[3] _ [-B*RealFns.Cos[theta] + x, y, B*RealFns.Sin[theta] + z]; x _ n + class.params[4]*N; y _ q + class.params[4]*Q; z _ l + class.params[4]*L; theta _ RealFns.ArcTan[-z,x]; class.normals[4] _ [-B*RealFns.Cos[theta] + x, y, B*RealFns.Sin[theta] + z]; }; 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; }; <> <> <> <> < roots[j+1] THEN {>> <> <> <> <> <> <<}; >> Init: PROC = { globalPoly _ Polynomial.Quartic[[0,0,0,0,0]]; }; Init[]; END.