-- File: SweepCastImpl.mesa
-- Last edited by Bier on December 18, 1982 1:41 am
-- Author: Eric Bier on July 4, 1983 11:05 pm
-- Contents: Ray Tracing Code for linear and revolute sweep shapes.

DIRECTORY
 CastRays,
CSG,
 Matrix3d,
 SV2d,
 SVBoundBox,
 SVFaces,
 SVFacesCast,
 SVLines2d,
 SVPolygon2d,
 SVPolygon3d,
 SweepCast,
 SweepGeometry,
 SVVector3d;


SweepCastImpl: PROGRAM
IMPORTS CastRays, SVBoundBox, SVFaces, SVFacesCast, SVLines2d, SVPolygon2d, SVPolygon3d
EXPORTS SweepCast =
BEGIN

Point2d: TYPE = SV2d.Point2d;
Point3d: TYPE = Matrix3d.Point3d;
Polygon: TYPE = SV2d.Polygon;
Vector: TYPE = SVVector3d.Vector;


Cone: TYPE = REF ConeObj;
ConeObj: TYPE = SVFaces.ConeObj;
DiskRing: TYPE = REF DiskRingObj;
DiskRingObj: TYPE = SVFaces.DiskRingObj;
Cylinder: TYPE = REF CylinderObj;
CylinderObj: TYPE = SVFaces.CylinderObj;
FaceHitRec: TYPE = REF FaceHitRecObj;
FaceHitRecObj: TYPE = SVFacesCast.FaceHitRecObj;
Plane: TYPE = SVPolygon3d.Plane;

Primitive: TYPE = REF PrimitiveObj;
PrimitiveObj: TYPE = CSG.PrimitiveObj;

Ray: TYPE = REF RayObj;
RayObj: TYPE = CastRays.RayObj;

Classification: TYPE = REF ClassificationObj;
ClassificationObj: TYPE = CastRays.ClassificationObj;

LinearMesh: TYPE = REF LinearMeshRecord;
LinearMeshRecord: TYPE = SweepGeometry.LinearMeshRecord;
RevoluteMesh: TYPE = REF RevoluteMeshRecord;
RevoluteMeshRecord: TYPE = SweepGeometry.RevoluteMeshRecord;

RevoSweepFaces: TYPE = REF RevoSweepFacesObj;
RevoSweepFacesObj: TYPE = SweepCast.RevoSweepFacesObj;
SubBoxesBody: TYPE = REF SubBoxesBodyObj;
SubBoxesBodyObj: TYPE = SweepCast.SubBoxesBodyObj;
RevoHitArray: TYPE = REF RevoHitArrayObj;
RevoHitArrayObj: TYPE = SweepCast.RevoHitArrayObj;
RevoHitRec: TYPE = REF RevoHitRecObj;
RevoHitRecObj: TYPE = SweepCast.RevoHitRecObj;

LinSweepFaces: TYPE = REF LinSweepFacesObj;
LinSweepFacesObj: TYPE = SweepCast.LinSweepFacesObj;
LinHitArray: TYPE = REF LinHitArrayObj;
LinHitArrayObj: TYPE = SweepCast.LinHitArrayObj;
LinHitRec: TYPE = REF LinHitRecObj;
LinHitRecObj: TYPE = SweepCast.LinHitRecObj;

globalRevoHits: RevoHitArray;
globalLinHits: LinHitArray;

almostZero: REAL = 1.0e-12;

LinCast: PUBLIC PROC [localRay: Ray, prim: Primitive, linMesh: LinearMesh, faces: LinSweepFaces] RETURNS [class: Classification] = {
-- A linear sweep in instance coordinates has its two complicated polygons parallel to the xy plane. Hence, all of the other (rectangular) faces are edge-on when projected on to the xy plane. We can use this fact to simplify their plane equations to line equations for ray casting.
-- We must do 2 ray-xyPolygon tests and "len" ray-edgeOnRectangle tests.
-- the shape may be convex so the ray may hit any even number of faces (in a perfect world), and any number of faces at all in our floating point world. We ignore these for now.

  poly: Polygon;
  faceHit: FaceHitRec;

 class ← CastRays.GetClassFromPool[];
 globalLinHits.count ← 0;

-- do 2 ray-xyPolygon tests.


-- front face;
 poly ← faces.poly;
 faceHit ← RayXYPolygon[localRay, poly, linMesh.array[1][1][3], [0, 0, 1]];
IF faceHit.count = 1 THEN {globalLinHits.count ← globalLinHits.count + 1;
         globalLinHits.array[globalLinHits.count]^ ←
         [faceHit.params[1], faceHit.normals[1], 0]};
 SVFacesCast.ReturnFaceHitToPool[faceHit];

-- back face
 faceHit ← RayXYPolygon[localRay, poly, linMesh.array[1][2][3], [0, 0, -1]];
IF faceHit.count = 1 THEN {globalLinHits.count ← globalLinHits.count + 1;
         globalLinHits.array[globalLinHits.count]^ ←
         [faceHit.params[1], faceHit.normals[1], linMesh.len+1]};
 SVFacesCast.ReturnFaceHitToPool[faceHit];

-- do ray-edgeOnRectangle tests
FOR i: NAT IN [0..faces.len) DO
  faceHit ← SVFacesCast.EdgeOnRectCast[faces[i], localRay];
  IF faceHit.count = 1 THEN {globalLinHits.count ← globalLinHits.count + 1;
         globalLinHits.array[globalLinHits.count]^ ←
         [faceHit.params[1], faceHit.normals[1], i+1]};
  SVFacesCast.ReturnFaceHitToPool[faceHit];
ENDLOOP;

-- now combine the results

SELECT globalLinHits.count FROM
  0 => {-- miss.
    CastRays.MakeClassAMiss[class]};
  1 => {-- ray just skims one face. Treat as a miss.
    CastRays.MakeClassAMiss[class]};
  ENDCASE => {-- Sort the hits in order of depth along ray. Assume for now that the ray doesn't skim edges or miss surfaces or go through cracks. Hence every other hit is an enter followed by an exit. If there is an odd number of hits, ignore the last one.
    iPlusOne: NAT;
    SortLinHits[globalLinHits];
    IF (globalLinHits.count/2)*2 # globalLinHits.count
     THEN globalRevoHits.count ← globalLinHits.count - 1;
    class.count ← globalLinHits.count;
    FOR i: NAT ← 1, i+2 UNTIL i > globalLinHits.count DO
     iPlusOne ← i+1;
     class.params[i] ← globalLinHits.array[i].param;
     class.params[iPlusOne] ← globalLinHits.array[iPlusOne].param;
      class.normals[i] ← globalLinHits.array[i].normal;
      class.normals[iPlusOne] ← globalLinHits.array[iPlusOne].normal;
      class.surfaces[i] ← NIL;
      class.surfaces[iPlusOne] ← NIL;
      class.classifs[i] ← FALSE;
      class.classifs[iPlusOne] ← TRUE;
      class.primitives[i] ← prim;
      class.primitives[iPlusOne] ← prim;
     ENDLOOP;
     class.classifs[globalLinHits.count+1] ← FALSE;
    }; -- end of 2 or more hits 
  
 }; -- end of LinCast

RevoCast: PUBLIC PROC [cameraPoint: Point2d, localRay: Ray, prim: Primitive, faces: RevoSweepFaces, boxes: SubBoxesBody] RETURNS [class: Classification] = {
-- A revolute sweep in instance coordinates has rotational symmetry about the y axis. Since it consists of straight line segments rotated about the y axis, its faces are cylindrical tubes, planar rings, and conical shouds, depending on whether the line segment in question was vertical, horizontal or neither respectively.
-- given a list of these "faces", we test the ray against each one for intersections. We may get as many as 4 intersections if the ray skims the cracks between two pairs of faces. However, we should only come up with two unique ray parameters t. We will award the intersection to the face which is higher.
 class ← CastRays.GetClassFromPool[];
 globalRevoHits.count ← 0;
FOR i: NAT IN[0..faces.len) DO
  WITH faces[i] SELECT FROM
  cone: Cone => {
   hits: FaceHitRec;
   IF NOT SVBoundBox.PointInBoundBox[cameraPoint, boxes[i]] THEN LOOP;
   hits ← SVFacesCast.ConeCast[cone, localRay];
   FOR j: NAT IN[1..hits.count] DO
    globalRevoHits.count ← globalRevoHits.count + 1;
    globalRevoHits.array[globalRevoHits.count]^ ← [hits.params[j], hits.normals[j], i];
   ENDLOOP;
   SVFacesCast.ReturnFaceHitToPool[hits];
   };
  disk: DiskRing => {
   hit: FaceHitRec;
   IF NOT SVBoundBox.PointInBoundBox[cameraPoint, boxes[i]] THEN LOOP;
   hit ← SVFacesCast.DiskRingCast[disk, localRay];
   IF hit.count = 1 THEN {
    globalRevoHits.count ← globalRevoHits.count + 1;
    globalRevoHits.array[globalRevoHits.count]^ ← [hit.params[1], hit.normals[1], i];
    };
   SVFacesCast.ReturnFaceHitToPool[hit];
   };
  cyl: Cylinder => {
   hits: FaceHitRec;
   IF NOT SVBoundBox.PointInBoundBox[cameraPoint, boxes[i]] THEN LOOP;
   hits ← SVFacesCast.CylinderCast[cyl, localRay];
   FOR j: NAT IN[1..hits.count] DO
    globalRevoHits.count ← globalRevoHits.count + 1;
    globalRevoHits.array[globalRevoHits.count]^ ← [hits.params[j], hits.normals[j], i];
   ENDLOOP;
   SVFacesCast.ReturnFaceHitToPool[hits];
   };
  ENDCASE => ERROR;
ENDLOOP; -- next face
  -- now see how many hits we have. We can have any number up to linesOfLatitude + 1 since the outline doesn't have to be convex. The number can be odd if the ray hits a crack between faces and both notice it. We proceed as follows: If there are 2 or more hits then we sort all of the hits by t. If 2 adjacent hits have similar y values and similar t values then we count them as the same hit and throw out the lower one. We we are done with this processing, if we have an odd number of hits, it is an error. Next we use the surface normals to see if the ray is entering or exitting at this point. If this does not alternate (starting with in), we have an error. Finally, we make a classification with the sorted params.
SELECT globalRevoHits.count FROM
  0 => {-- miss.
    CastRays.MakeClassAMiss[class]};
  1 => {-- ray just skims one face. Treat as a miss.
    CastRays.MakeClassAMiss[class]};
  ENDCASE => {-- Sort the hits in order of depth along ray. Assume for now that the ray doesn't skim edges or miss surfaces or go through cracks. Hence every other hit is an enter followed by an exit. If there is an odd number of hits, ignore the last one.
    iPlusOne: NAT;
    SortRevoHits[globalRevoHits];
    IF (globalRevoHits.count/2)*2 # globalRevoHits.count
     THEN globalRevoHits.count ← globalRevoHits.count - 1;
    class.count ← globalRevoHits.count;
    FOR i: NAT ← 1, i+2 UNTIL i > globalRevoHits.count DO
     iPlusOne ← i+1;
     class.params[i] ← globalRevoHits.array[i].param;
     class.params[iPlusOne] ← globalRevoHits.array[iPlusOne].param;
      class.normals[i] ← globalRevoHits.array[i].normal;
      class.normals[iPlusOne] ← globalRevoHits.array[iPlusOne].normal;
      class.surfaces[i] ← NIL;
      class.surfaces[iPlusOne] ← NIL;
      class.classifs[i] ← FALSE;
      class.classifs[iPlusOne] ← TRUE;
      class.primitives[i] ← prim;
      class.primitives[iPlusOne] ← prim;
     ENDLOOP;
     class.classifs[globalRevoHits.count+1] ← FALSE;
    }; -- end of 2 or more hits 
  
 }; -- end of RevoCast

RevoCastNoBBoxes: PUBLIC PROC [localRay: Ray, prim: Primitive, faces: RevoSweepFaces] RETURNS [class: Classification] = {
-- A revolute sweep in instance coordinates has rotational symmetry about the y axis. Since it consists of straight line segments rotated about the y axis, its faces are cylindrical tubes, planar rings, and conical shouds, depending on whether the line segment in question was vertical, horizontal or neither respectively.
-- given a list of these "faces", we test the ray against each one for intersections. We may get as many as 4 intersections if the ray skims the cracks between two pairs of faces. However, we should only come up with two unique ray parameters t. We will award the intersection to the face which is higher.
 class ← CastRays.GetClassFromPool[];
 globalRevoHits.count ← 0;
FOR i: NAT IN[0..faces.len) DO
  WITH faces[i] SELECT FROM
  cone: Cone => {
   hits: FaceHitRec;
   hits ← SVFacesCast.ConeCast[cone, localRay];
   FOR j: NAT IN[1..hits.count] DO
    globalRevoHits.count ← globalRevoHits.count + 1;
    globalRevoHits.array[globalRevoHits.count]^ ← [hits.params[j], hits.normals[j], i];
   ENDLOOP;
   SVFacesCast.ReturnFaceHitToPool[hits];
   };
  disk: DiskRing => {
   hit: FaceHitRec;
   hit ← SVFacesCast.DiskRingCast[disk, localRay];
   IF hit.count = 1 THEN {
    globalRevoHits.count ← globalRevoHits.count + 1;
    globalRevoHits.array[globalRevoHits.count]^ ← [hit.params[1], hit.normals[1], i];
    };
   SVFacesCast.ReturnFaceHitToPool[hit];
   };
  cyl: Cylinder => {
   hits: FaceHitRec;
   hits ← SVFacesCast.CylinderCast[cyl, localRay];
   FOR j: NAT IN[1..hits.count] DO
    globalRevoHits.count ← globalRevoHits.count + 1;
    globalRevoHits.array[globalRevoHits.count]^ ← [hits.params[j], hits.normals[j], i];
   ENDLOOP;
   SVFacesCast.ReturnFaceHitToPool[hits];
   };
  ENDCASE => ERROR;
ENDLOOP; -- next face
  -- now see how many hits we have. We can have any number up to linesOfLatitude + 1 since the outline doesn't have to be convex. The number can be odd if the ray hits a crack between faces and both notice it. We proceed as follows: If there are 2 or more hits then we sort all of the hits by t. If 2 adjacent hits have similar y values and similar t values then we count them as the same hit and throw out the lower one. We we are done with this processing, if we have an odd number of hits, it is an error. Next we use the surface normals to see if the ray is entering or exitting at this point. If this does not alternate (starting with in), we have an error. Finally, we make a classification with the sorted params.
SELECT globalRevoHits.count FROM
  0 => {-- miss.
    CastRays.MakeClassAMiss[class]};
  1 => {-- ray just skims one face. Treat as a miss.
    CastRays.MakeClassAMiss[class]};
  ENDCASE => {-- Sort the hits in order of depth along ray. Assume for now that the ray doesn't skim edges or miss surfaces or go through cracks. Hence every other hit is an enter followed by an exit. If there is an odd number of hits, ignore the last one.
    iPlusOne: NAT;
    SortRevoHits[globalRevoHits];
    IF (globalRevoHits.count/2)*2 # globalRevoHits.count
     THEN globalRevoHits.count ← globalRevoHits.count - 1;
    class.count ← globalRevoHits.count;
    FOR i: NAT ← 1, i+2 UNTIL i > globalRevoHits.count DO
     iPlusOne ← i+1;
     class.params[i] ← globalRevoHits.array[i].param;
     class.params[iPlusOne] ← globalRevoHits.array[iPlusOne].param;
      class.normals[i] ← globalRevoHits.array[i].normal;
      class.normals[iPlusOne] ← globalRevoHits.array[iPlusOne].normal;
      class.surfaces[i] ← NIL;
      class.surfaces[iPlusOne] ← NIL;
      class.classifs[i] ← FALSE;
      class.classifs[iPlusOne] ← TRUE;
      class.primitives[i] ← prim;
      class.primitives[iPlusOne] ← prim;
     ENDLOOP;
     class.classifs[globalRevoHits.count+1] ← FALSE;
    }; -- end of 2 or more hits 
  
 }; -- end of RevoCast
OddHits: ERROR = CODE;
InOutsDontAlternate: ERROR = CODE;

HitsACrack: PRIVATE PROC [localRay: Ray, t, y: REAL, meshRecord: RevoluteMesh] RETURNS [matches: BOOL, index: NAT] = {
 x, z, r: REAL;
FOR i: NAT IN[1..meshRecord.linesOfLatitude] DO
  IF AlmostEqualHeights[y, meshRecord.array[i][1][2]] THEN {
   x ← localRay.basePt[1] + t*localRay.direction[1];
   z ← localRay.basePt[3] + t*localRay.direction[3];
   r ← meshRecord.array[i][meshRecord.linesOfLongitude][1];
   IF x*x + z*z = r*r THEN
    {matches ← TRUE; index ← i; RETURN}
   ELSE {matches ← FALSE; RETURN}};
ENDLOOP;
 matches ← FALSE;
 }; -- end of HitsACrack

HitsNeighbors: PRIVATE PROC [localRay: Ray, index: NAT, mesh: RevoluteMesh] RETURNS [BOOL] = {
IF index = 1 THEN RETURN[HitsDisk[localRay, mesh.array[2][1][2], mesh.array[2][mesh.linesOfLongitude][1]]];
IF index = mesh.linesOfLatitude THEN RETURN[
  HitsDisk[localRay, mesh.array[mesh.linesOfLatitude-1][1][2],
   mesh.array[mesh.linesOfLatitude-1][mesh.linesOfLongitude][1]]];
RETURN[HitsDisk[localRay, mesh.array[index-1][1][2],
   mesh.array[index-1][mesh.linesOfLongitude][1]]
   OR
   HitsDisk[localRay, mesh.array[index+1][1][2],
   mesh.array[index+1][mesh.linesOfLongitude][1]]];
   };
   
HitsDisk: PRIVATE PROC [localRay: Ray, y: REAL, r: REAL] RETURNS [success: BOOL] = {
-- checks ray for intersections with the disk in plane Y=y of radius r.
-- do 1 ray-plane intersection and test to see if it is inside the circle
-- Plane is y = y;
-- For ray, y = p2 + t*d2;
-- Solving, y = p2 + t*d2; (y - p2)/d2 = t;
 x, z, t: REAL;
 p: Point3d ← localRay.basePt;
 d: Vector ← localRay.direction;
 success ← FALSE;
IF Abs[d[2]] < almostZero THEN RETURN; -- ray parallel to plane of disk
 t ← (y - p[2])/d[2];
IF t > 0 THEN {
  x ← p[1] + t*d[1];
  z ← p[3] + t*d[3];
  IF x*x + z*z <= r*r THEN success ← TRUE;
  };
 }; -- end of HitsDisk
  

AlmostEqualParams: PRIVATE PROC [t1, t2: REAL] RETURNS [BOOL] = {
RETURN[Abs[t2-t1] < 1.0e-2];
 };

AlmostEqualHeights: PRIVATE PROC [h1, h2: REAL] RETURNS [BOOL] = {
RETURN[Abs[h2-h1] < 1.0];
 };

DeleteFromHits: PRIVATE PROC [hits: RevoHitArray, index: NAT] = {
-- shift all elements above this one down and reduce count by 1
 temp: RevoHitRec ← hits.array[index];
FOR i: NAT IN[index..hits.count-1] DO
  hits.array[i] ← hits.array[i+1];
ENDLOOP;
 hits.array[hits.count] ← temp; -- replace the storage for later use
 hits.count ← hits.count - 1;
 }; -- end of DeleteFromHits


RayGoesIn: PRIVATE PROC [ray: Ray, hit: RevoHitRec] RETURNS [BOOL] = {
-- TRUE means the ray is entering the plane from the normal side.
 hitPoint: Point3d;
 plane: Plane;
 basePtOnNormalSide: BOOL;
 hitPoint[1] ← ray.basePt[1] + hit.param*ray.direction[1];
 hitPoint[2] ← ray.basePt[2] + hit.param*ray.direction[2];
 hitPoint[3] ← ray.basePt[3] + hit.param*ray.direction[3];
 plane ← SVPolygon3d.PlaneFromPointAndNormal[hitPoint, hit.normal];
 basePtOnNormalSide ← SVPolygon3d.PointOnNormalSideOfPlane[ray.basePt, plane];
-- if basePt is on the normal side then we are going in.
RETURN[basePtOnNormalSide];
 }; -- end of RayGoesIn

SortLinHits: PROC [hits: LinHitArray] = {
 temp: LinHitRec;
FOR i: NAT IN[1..hits.count-1] DO
  FOR j: NAT IN[1..hits.count-i] DO
   IF hits.array[j].param > hits.array[j+1].param THEN {
    temp ← hits.array[j];
    hits.array[j] ← hits.array[j+1];
    hits.array[j+1] ← temp;};
  ENDLOOP;
ENDLOOP;
 };

SortRevoHits: PROC [hits: RevoHitArray] = {
 temp: RevoHitRec;
FOR i: NAT IN[1..hits.count-1] DO
  FOR j: NAT IN[1..hits.count-i] DO
   IF hits.array[j].param > hits.array[j+1].param THEN {
    temp ← hits.array[j];
    hits.array[j] ← hits.array[j+1];
    hits.array[j+1] ← temp;};
  ENDLOOP;
ENDLOOP;
 };

Abs: PROC [r: REAL] RETURNS [REAL] = {
RETURN[IF r >= 0 THEN r ELSE -r];
 };

Point3dToPoint2d: PRIVATE PROC [p3: Point3d] RETURNS [p2: Point2d] = {
 p2[1] ← p3[1];
 p2[2] ← p3[2];
 };

MakeLinSweepFaces: PUBLIC PROC [linMesh: LinearMesh] RETURNS [faces: LinSweepFaces] = {
-- interate through consequtive pairs of vertices of the front polygon of the linear mesh. Make linMesh.len faces using SVFaces.EdgeOnRectFromTrigLineSeg
 lastPoint, thisPoint: Point2d;
 seg: SV2d.TrigLineSeg;
 poly: Polygon ← SVPolygon2d.CreatePoly[linMesh.len];
-- extract the defining polygon
FOR i: NAT IN[1..linMesh.len] DO
  poly ← SVPolygon2d.AddPolyPoint[poly, [linMesh.array[i][1][1], linMesh.array[i][1][2]] ];
ENDLOOP;
-- make the edge-on-rectangles
 faces ← NEW[LinSweepFacesObj[linMesh.len]];
 faces.len ← linMesh.len;
 faces.poly ← poly;
 lastPoint ← Point3dToPoint2d[linMesh.array[linMesh.len][1]];
FOR i: NAT IN[1..linMesh.len] DO
  thisPoint ← Point3dToPoint2d[linMesh.array[i][1]];
  seg ← SVLines2d.CreateTrigLineSeg[lastPoint, thisPoint];
  faces[i-1] ← SVFaces.EdgeOnRectFromTrigLineSeg[seg: seg, frontZ: linMesh.array[1][1][3],
   backZ: linMesh.array[1][2][3]];
  lastPoint ← thisPoint;
ENDLOOP;
 };

MakeRevoSweepFaces: PUBLIC PROC [revMesh: RevoluteMesh] RETURNS [faces: RevoSweepFaces] = {
-- currently the 2d polygon from which a revMesh is derived just happens to be in the last line of longitude (ignoring the z component which is zero); I assume the revolute mesh is defined clockwise (in the positive x half-plane seen from positive z). A trig line segment is oriented, so we don't lose the clockwise information when we make one. This allows SVFaces to orient normals properly.
 lastPoint, thisPoint: Point2d;
 lastLong: NAT ← revMesh.linesOfLongitude;
 seg: SV2d.TrigLineSeg;
 faces ← NEW[RevoSweepFacesObj[revMesh.linesOfLatitude+1]];
 faces.len ← revMesh.linesOfLatitude+1;
-- the top of a sweep is a disk
 lastPoint ← Point3dToPoint2d[revMesh.array[1][lastLong]];
 seg ← SVLines2d.CreateTrigLineSeg[[0,lastPoint[2]], lastPoint];
 faces[0] ← SVFaces.DiskRingFromTrigLineSeg[seg];
FOR i: NAT IN[2..revMesh.linesOfLatitude] DO -- the intermediate faces
  thisPoint ← Point3dToPoint2d[revMesh.array[i][lastLong]];
  seg ← SVLines2d.CreateTrigLineSeg[lastPoint, thisPoint];
  SELECT seg.line.theta FROM
  0, 180 => -- make a disk
     faces[i-1] ← SVFaces.DiskRingFromTrigLineSeg[seg];
  IN (0..90), IN (90..180), IN (-90..0), IN (-180..-90) => -- make a cone
     faces[i-1] ← SVFaces.ConeFromTrigLineSeg[seg];
  90, -90 => -- make a cylinder
     faces[i-1] ← SVFaces.CylinderFromTrigLineSeg[seg];
  ENDCASE => ERROR;
  lastPoint ← thisPoint;
ENDLOOP;
-- the bottom face is a disk
 seg ← SVLines2d.CreateTrigLineSeg[lastPoint, [0,lastPoint[2]]];
 faces[revMesh.linesOfLatitude] ← SVFaces.DiskRingFromTrigLineSeg[seg];
 faces.mesh ← revMesh;
 };

RayXYPolygon: PRIVATE PROC [localRay: Ray, poly: Polygon, depth: REAL, normal: Vector] RETURNS [linHit: FaceHitRec] = {
-- Find intersection of localRay with the plane of the polygon. Use PointInPoly test to see if this point is inside of the polygon. Return a filled LinHitRec if so.
-- plane equation Z=z. Ray equation is, as usual, z = p3 + t*d3. Solve for t.
-- (z-p3)/d3 = t;
 x, y, z, t: REAL;
 pointClass: SVPolygon2d.PointPolyClass;
 linHit ← SVFacesCast.GetFaceHitFromPool[];
 z ← depth;
IF Abs[localRay.direction[3]] < almostZero THEN {linHit.count ← 0; RETURN};
 t ← (z - localRay.basePt[3])/localRay.direction[3];
IF t <= 0 THEN {linHit.count ← 0; RETURN};
 x ← localRay.basePt[1] + t*localRay.direction[1];
 y ← localRay.basePt[2] + t*localRay.direction[2];
 pointClass ← SVPolygon2d.BoysePointInPoly[[x,y], poly];
IF pointClass = in OR pointClass = on THEN {linHit.count ← 1;
  linHit.params[1] ← t;
  linHit.normals[1] ← normal;
  }
ELSE linHit.count ← 0;
 }; -- end of RayXYPolygon

  
Init: PROC = {
 globalRevoHits ← NEW[RevoHitArrayObj];
 globalRevoHits.count ← 0;
FOR i: NAT IN[1..SweepGeometry.maxLinesOfLat+1] DO
  globalRevoHits.array[i] ← NEW[RevoHitRecObj];
ENDLOOP;
 globalLinHits ← NEW[LinHitArrayObj];
 globalLinHits.count ← 0;
FOR i: NAT IN[1..SweepGeometry.maxMeshLen-1] DO
  globalLinHits.array[i] ← NEW[LinHitRecObj];
ENDLOOP;
 };

Init[];

END.