<> <> <> <> DIRECTORY RealFns, SV2d, SV3d, SVVector3d; SVVector3dImpl: CEDAR PROGRAM IMPORTS RealFns EXPORTS SVVector3d = BEGIN Point3d: TYPE = SV3d.Point3d; Vector3d: TYPE = SV3d.Vector3d; Vector2d: TYPE = SV2d.Vector2d; galmostZero: REAL = 1.0e-10; -- Must be less than 1.0e-8 I've discovered. This is an old comment. VectorFromAngles: PUBLIC PROC [azimuth, slope: REAL] RETURNS [vector: Vector3d] = { <> <> sin sin cos sin cos vector[1] _ cos vector[3] _ -cos vector[2] _ sin }; VectorToPlane: PRIVATE PROC [v: Vector3d, xaxis: Vector3d, normal: Vector3d] RETURNS [v2: Vector2d] = { <> <> <<1) v VectorFromAngle in two dimensions.>> <<2) Normalize xaxis. Then v >> <> dotProd, crossProdMag: REAL; crossProd: Vector3d; xaxis _ Normalize[xaxis]; dotProd _ DotProduct[v, xaxis]; crossProd _ CrossProduct[xaxis, v]; <> IF NOT SameSense[crossProd, normal] THEN dotProd _ -dotProd; crossProdMag _ Magnitude[crossProd]; v2 _ [dotProd, crossProdMag]; }; <> SameSense: PUBLIC PROC [v1, v2: Vector3d] RETURNS [BOOL] = { <> RETURN[TRUE]; }; <> VectorFromPlane: PUBLIC PROC [v2: Vector2d, xaxis: Vector3d, normal: Vector3d] RETURNS [v: Vector3d] = { <> <> <> <> }; <<>> <> VectorPlusAngleInPlane: PUBLIC PROC [v: Vector3d, degrees: REAL, normal: Vector3d] RETURNS [rotatedV: Vector2d] = { <> }; Add: PUBLIC PROC [v1: Vector3d, v2: Vector3d] RETURNS [sumV: Vector3d] = { sumV[1] _ v1[1] + v2[1]; sumV[2] _ v1[2] + v2[2]; sumV[3] _ v1[3] + v2[3]; }; Sub: PUBLIC PROC [v1: Vector3d, v2: Vector3d] RETURNS [v1Minusv2: Vector3d] = { v1Minusv2[1] _ v1[1] - v2[1]; v1Minusv2[2] _ v1[2] - v2[2]; v1Minusv2[3] _ v1[3] - v2[3]; }; Scale: PUBLIC PROC [v: Vector3d, scalar: REAL] RETURNS [scaledV: Vector3d] = { scaledV[1] _ v[1]*scalar; scaledV[2] _ v[2]*scalar; scaledV[3] _ v[3]*scalar; }; Normalize: PUBLIC PROC [v: Vector3d] RETURNS [normV: Vector3d] = { <> mag: REAL _ Magnitude[v]; normV[1] _ v[1]/mag; normV[2] _ v[2]/mag; normV[3] _ v[3]/mag; }; Negate: PUBLIC PROC [v: Vector3d] RETURNS [negV: Vector3d] = { negV[1] _ -v[1]; negV[2] _ -v[2]; negV[3] _ -v[3]; }; DotProduct: PUBLIC PROC [v1: Vector3d, v2: Vector3d] RETURNS [scalar: REAL] = { scalar _ v1[1]*v2[1] + v1[2]*v2[2] + v1[3]*v2[3]; }; ElementwiseProduct: PUBLIC PROC [v1, v2: Vector3d] RETURNS [prod: Vector3d] = { prod[1] _ v1[1] * v2[1]; prod[2] _ v1[2] * v2[2]; prod[3] _ v1[3] * v2[3]; }; CrossProduct: PUBLIC PROC [v1: Vector3d, v2: Vector3d] RETURNS [prodV: Vector3d] = { <<| i j k |>> <<| v1x v1y v1z |=(v1y*v2z - v1z*v2y) i + (v1z*v2x - v1x*v2z) j>> <<| v2x v2y v2z |(v1x*v2y - v1y*v2x) k>> prodV[1] _ v1[2]*v2[3] - v1[3]*v2[2]; prodV[2] _ v1[3]*v2[1] - v1[1]*v2[3]; prodV[3] _ v1[1]*v2[2] - v1[2]*v2[1]; }; AngleCCWBetweenVectors: PUBLIC PROC [v1: Vector3d, v2: Vector3d] RETURNS [degrees: REAL] = { <> v1hat, v2hat, cross: Vector3d; sin, cos: REAL; v1hat _ Normalize[v1]; v2hat _ Normalize[v2]; cross _ CrossProduct[v1hat, v2hat]; sin _ Magnitude[cross]; cos _ DotProduct[v1hat, v2hat]; degrees _ RealFns.ArcTanDeg[sin, cos]; }; Magnitude: PUBLIC PROC [v: Vector3d] RETURNS [magnitude: REAL] = { magnitude _ RealFns.SqRt[(v[1]*v[1]+v[2]*v[2]+v[3]*v[3])]; }; Distance: PUBLIC PROC [p1, p2: Point3d] RETURNS [dist: REAL] = { dist _ Magnitude[Sub[p2, p1]]; }; MagnitudeSquared: PUBLIC PROC [v: Vector3d] RETURNS [magSquared: REAL]= { magSquared _ (v[1]*v[1]+v[2]*v[2]+v[3]*v[3]); }; DistanceSquared: PUBLIC PROC [p1, p2: Point3d] RETURNS [distSquared: REAL] = { distSquared _ MagnitudeSquared[Sub[p2, p1]]; }; VectorFromPoints: PUBLIC PROC [tail, head: Point3d] RETURNS [vector: Vector3d] = { vector[1] _ head[1] - tail[1]; vector[2] _ head[2] - tail[2]; vector[3] _ head[3] - tail[3]; }; AlmostZero: PRIVATE PROC [r: REAL] RETURNS [BOOL] = { RETURN[ABS[r] < galmostZero]; }; AlmostEqual: PRIVATE PROC [a, b: REAL] RETURNS [BOOL] = { RETURN[ABS[(a-b)/a] < 0.01]; }; <> <> <> <> <> <> <> <<>> <> <> <> <> <<>> <> <> <> <> <> <<};>> <<>> Parallel: PUBLIC PROC [v1, v2: Vector3d] RETURNS [BOOL] = { <> crossProd: REAL _ Magnitude[CrossProduct[Normalize[v1], Normalize[v2]]]; answer: BOOL; answer _ ABS[crossProd] < 1.0E-5; RETURN[answer]; }; <> <> <> << ELSE IF c = 0 THEN RETURN[TRUE] -- both vectors have only a y component>> << ELSE RETURN[b=c];>> <<}>> <> <> <> <> <<}>> <> <> <> <<};>> <<};>> <<}; -- end of AllEqual>> <<>> <> <> <> <> <> << ELSE IF AlmostZero[c] THEN RETURN[TRUE] -- both vectors have only a y component>> << ELSE RETURN[AlmostEqual[b, c]];>> <<}>> <> <> <> <> <<}>> <> <> <> <<};>> <<};>> <<}; -- end of AllAlmostEqual>> <<>> Perpendicular: PUBLIC PROC [v1, v2: Vector3d] RETURNS [BOOL] = { RETURN[DotProduct[v1, v2] = 0]; }; Vector2DAsXYVector: PUBLIC PROC [vXY: Vector2d] RETURNS [vZeroZ: Vector3d] = { vZeroZ[1] _ vXY[1]; vZeroZ[2] _ vXY[2]; vZeroZ[3] _ 0; }; Vector2DAsYZVector: PUBLIC PROC [vYZ: Vector2d] RETURNS [vZeroX: Vector3d] = { vZeroX[1] _ 0; vZeroX[2] _ vYZ[1]; vZeroX[3] _ vYZ[2]; }; Vector2DAsZXVector: PUBLIC PROC [vZX: Vector2d] RETURNS [vZeroY: Vector3d] = { vZeroY[1] _ vZX[2]; vZeroY[2] _ 0; vZeroY[3] _ vZX[1]; }; ProjectOntoXYPlane: PUBLIC PROC [v: Vector3d] RETURNS [v2d: Vector2d] = { v2d[1] _ v[1]; v2d[2] _ v[2]; }; ProjectOntoYZPlane: PUBLIC PROC [v: Vector3d] RETURNS [v2d: Vector2d] = { v2d[1] _ v[2]; v2d[2] _ v[3]; }; ProjectOntoZXPlane: PUBLIC PROC [v: Vector3d] RETURNS [v2d: Vector2d] = { v2d[1] _ v[3]; v2d[2] _ v[1]; }; END.