DIRECTORY RealFns, SV2d, SV3d, SVVector3d; SVVector3dImpl: PROGRAM IMPORTS RealFns EXPORTS SVVector3d = BEGIN Point3d: TYPE = SV3d.Point3d; Vector: TYPE = SV3d.Vector; Vector2d: TYPE = SV2d.Vector2d; galmostZero: REAL = 1.0e-10; -- Must be less than 1.0e-8 I've discovered. VectorToPlane: PRIVATE PROC [v: Vector, xaxis: Vector, normal: Vector] RETURNS [v2: Vector2d] = { dotProd, crossProdMag: REAL; crossProd: Vector; 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: Vector] RETURNS [BOOL] = { RETURN[TRUE]; }; VectorFromPlane: PUBLIC PROC [v2: Vector2d, xaxis: Vector, normal: Vector] RETURNS [v: Vector] = { }; VectorPlusAngleInPlane: PUBLIC PROC [v: Vector, degrees: REAL, normal: Vector] RETURNS [rotatedV: Vector2d] = { }; Add: PUBLIC PROC [v1: Vector, v2: Vector] RETURNS [sumV: Vector] = { sumV[1] _ v1[1] + v2[1]; sumV[2] _ v1[2] + v2[2]; sumV[3] _ v1[3] + v2[3]; }; Sub: PUBLIC PROC [v1: Vector, v2: Vector] RETURNS [v1Minusv2: Vector] = { v1Minusv2[1] _ v1[1] - v2[1]; v1Minusv2[2] _ v1[2] - v2[2]; v1Minusv2[3] _ v1[3] - v2[3]; }; Scale: PUBLIC PROC [v: Vector, scalar: REAL] RETURNS [scaledV: Vector] = { scaledV[1] _ v[1]*scalar; scaledV[2] _ v[2]*scalar; scaledV[3] _ v[3]*scalar; }; Normalize: PUBLIC PROC [v: Vector] RETURNS [normV: Vector] = { mag: REAL _ Magnitude[v]; normV[1] _ v[1]/mag; normV[2] _ v[2]/mag; normV[3] _ v[3]/mag; }; Negate: PUBLIC PROC [v: Vector] RETURNS [negV: Vector] = { negV[1] _ -v[1]; negV[2] _ -v[2]; negV[3] _ -v[3]; }; DotProduct: PUBLIC PROC [v1: Vector, v2: Vector] RETURNS [scalar: REAL] = { scalar _ v1[1]*v2[1] + v1[2]*v2[2] + v1[3]*v2[3]; }; ElementwiseProduct: PUBLIC PROC [v1, v2: Vector] RETURNS [prod: Vector] = { prod[1] _ v1[1] * v2[1]; prod[2] _ v1[2] * v2[2]; prod[3] _ v1[3] * v2[3]; }; CrossProduct: PUBLIC PROC [v1: Vector, v2: Vector] RETURNS [prodV: Vector] = { 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]; }; Magnitude: PUBLIC PROC [v: Vector] 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: Vector] 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: Vector] = { 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: Vector] RETURNS [BOOL] = { a,b,c: REAL; IF AlmostZero[v2[1]] THEN { IF NOT AlmostZero[v1[1]] THEN RETURN[FALSE] ELSE a _ 0} ELSE { -- v2[1] # 0 IF AlmostZero[v1[1]] THEN RETURN[FALSE] ELSE a _ v1[1]/v2[1]}; IF AlmostZero[v2[2]] THEN { IF NOT AlmostZero[v1[2]] THEN RETURN[FALSE] ELSE b _ 0} ELSE { -- v2[2] # 0 IF AlmostZero[v1[2]] THEN RETURN[FALSE] ELSE b _ v1[2]/v2[2]}; IF AlmostZero[v2[3]] THEN { IF NOT AlmostZero[v1[3]] THEN RETURN[FALSE] ELSE c _ 0} ELSE { -- v2[3] # 0 IF AlmostZero[v1[3]] THEN RETURN[FALSE] ELSE c _ v1[3]/v2[3]}; RETURN[AllAlmostEqual[a, b, c]]; }; AllAlmostEqual: PRIVATE PROC [a, b, c: REAL] RETURNS [BOOL] = { almostZero: REAL _ 1.0e-10; IF AlmostZero[a] THEN { IF AlmostZero[b] THEN RETURN[TRUE] -- both vectors have only a z component ELSE IF AlmostZero[c] THEN RETURN[TRUE] -- both vectors have only a y component ELSE RETURN[AlmostEqual[b, c]]; } ELSE { IF AlmostZero[b] THEN { IF AlmostZero[c] THEN RETURN[TRUE] -- both vectors have only an x component ELSE RETURN[AlmostEqual[a, c]]; } ELSE { IF AlmostZero[c] THEN RETURN[AlmostEqual[a, b]] ELSE RETURN[AlmostEqual[a, b] AND AlmostEqual[b, c]]; }; }; }; -- end of AllAlmostEqual Perpendicular: PUBLIC PROC [v1, v2: Vector] RETURNS [BOOL] = { RETURN[DotProduct[v1, v2] = 0]; }; Vector2DAsXYVector: PUBLIC PROC [vXY: Vector2d] RETURNS [vZeroZ: Vector] = { vZeroZ[1] _ vXY[1]; vZeroZ[2] _ vXY[2]; vZeroZ[3] _ 0; }; Vector2DAsYZVector: PUBLIC PROC [vYZ: Vector2d] RETURNS [vZeroX: Vector] = { vZeroX[1] _ 0; vZeroX[2] _ vYZ[1]; vZeroX[3] _ vYZ[2]; }; Vector2DAsZXVector: PUBLIC PROC [vZX: Vector2d] RETURNS [vZeroY: Vector] = { vZeroY[1] _ vZX[2]; vZeroY[2] _ 0; vZeroY[3] _ vZX[1]; }; ProjectOntoXYPlane: PUBLIC PROC [v: Vector] RETURNS [v2d: Vector2d] = { v2d[1] _ v[1]; v2d[2] _ v[2]; }; ProjectOntoYZPlane: PUBLIC PROC [v: Vector] RETURNS [v2d: Vector2d] = { v2d[1] _ v[2]; v2d[2] _ v[3]; }; ProjectOntoZXPlane: PUBLIC PROC [v: Vector] RETURNS [v2d: Vector2d] = { v2d[1] _ v[3]; v2d[2] _ v[1]; }; END. аFile: SVVector3dImpl.mesa Author: Eric Bier before January 12, 1983 1:35 pm Last edited by Bier on September 5, 1984 11:05:28 am PDT Contents: Vector addition and multiplication and stuff like that. Not yet tested. Not used anywhere yet. Find a two-dimensional vector v2 with the same length as v, such that the angle between v2 and the x-axis in 2-space is the same as the angle between v and xaxis in 3-space (where the positive sense of the angle is counter-clockwise around "normal"). Two methods to try: 1) v o xaxis = |v| |xaxis| cosq. |xaxis X v| = |v||xaxis| sinq. The ratio = tanq so we can easily find q and then use VectorFromAngle in two dimensions. 2) Normalize xaxis. Then v o xaxis = |v| cosq, |xaxis X v|= |v| sinq. These are the two components of v2. The second method seems the more efficient since no ArcTan is involved. Two square roots are involved to normalize xaxis and find the magnitude of |v X xaxis|. Does this have the same sense or opposite sense to normal? If opposite, then change the sign of the dot product. Not yet implemented. Are v1 and v2 more nearly parallel or anti-parallel? Taking the sign of the dot product would work but is more expensive then needed. Compare the signs of the components and let the majority vote win. Not yet implemented. The inverse of VectorToPlane. Given a vector in 2-space, and a plane p and x-axis vector in 3-space, find the vector v in plane p whose angle with xaxis is the same as the angle of the absolute angle of v2 in 2-space. Method: We know that v o normal = 0, we know the value of v o xaxis and we know the desired length of v. Let normal = [a,b,c], v = [x,y,z], xaxis = [p,q,r]. Let v2[1] (=|v| cosq ) be called A. Then we know: ax + by + cz = 0. px + qy + rz = A. Not yet implemented. Like SVVector2d.VectorPlusAngle, but we must specify a plane to rotate in. "normal" is a vector normal to such a plane (and must be normal to v as well for this operation to make sense. The sense of normal determines which rotational direction is counter-clockwise by the right-hand rule. We proceed as follows: Returns the unit vector with the same direction as v. | i j k | | v1x v1y v1z |=(v1y*v2z - v1z*v2y) i + (v1z*v2x - v1x*v2z) j | v2x v2y v2z |(v1x*v2y - v1y*v2x) k AllEqual: PRIVATE PROC [a, b, c: REAL] RETURNS [BOOL] = { IF a = 0 THEN { IF b = 0 THEN RETURN[TRUE]-- both vectors have only a z component ELSE IF c = 0 THEN RETURN[TRUE] -- both vectors have only a y component ELSE RETURN[b=c]; } ELSE { IF b = 0 THEN { IF c = 0 THEN RETURN[TRUE]-- both vectors have only an x component ELSE RETURN[a=c]; } ELSE { IF c = 0 THEN RETURN[a=b] ELSE RETURN[a=b AND b = c]; }; }; }; -- end of AllEqual Э╩ ╖Ц "cedar" styleШIheadЪЬЩIprocЪЬ1Щ1LЪЬ8Щ8LЪЬAЩALШЪ╧k Ш LЪЬШLШLЪЬШLЪЬШLШЧЪЬЭШLЪЭЬШLЪЭЬ ШЧLШLЪЭШШLЪЬ ЭЬШLЪЬЭЬШLЪЬ ЭЬШLШLЪЬ ЭЬ8ШIШLЩ'ЧЧЪ╧n ЬЭЬЭЬ,ЭЬШaLЪЬПЩПLЪЬ╧mЬ╧gЬЯЬаЬаЬаЬаЬ0ЩЬLЪЬЯЬаЬЯЬаЬаЬ%ЩmLЪЬФаЬЯЬаЬЩаLЪЬЭЬШLЪЬШLЪЬШLЪЬШЪЬ#Ш#LЪЬqЩqЧLЪЭЬЭЬЭЬШLЪЬ5Щ5LЪЬЭЬШLЪЬШLЪЬШLЪЬШLЪЬШЧЪЮЬЭЬЭЬ ЭЬШ:LЪЬ4Ш4LЪЬШЧЪ Ю ЬЭЬЭЬЭЬ ЭЬШKLЪЬ1Ш1LЪЬШЧЪЮЬЭЬЭЬЭЬШKLЪЬШLЪЬШLЪЬШLЪЬШLШЧЪЮЬЭЬЭЬЭЬШNLЪЬЩLЪЬ=Щ=LЪЬ$Щ$LЪЬ%Ш%LЪЬ%Ш%LЪЬ%Ш%LЪЬШЧLШЪ Ю ЬЭЬЭЬ ЭЬ ЭЬШ@LЪЬ:Ш:LЪЬШLШЧЪ ЮЬЭЬЭЬЭЬЭЬШ@LЪЬШLШЧLШЪЮЬЭЬ ЭЬЭЬШGLЪЬ-Ш-LЪЬШЧЪЮЬЭЬЭЬЭЬШNLЪЬ,Ш,LШЧLШЪЮЬЭЬЭЬШPLЪЬШLЪЬШLЪЬШLЪЬШLШЧLШЪЮ ЬЭЬЭЬЭЬЭЬЭЬШ5LЪЭЬЭЬШLШLШЧЪЮЬЭЬЭЬЭЬЭЬЭЬШ9LЪЭЬЭЬШLШLШЧЪ ЮЬЭЬЭЬЭЬЭЬШ9LЪЬЭЬШЪЭЬЭЬШLЪЭЬЭЬЭЬЭЬЭЬЭЬШ7ЧЪЭЬ╧cШLЪ ЭЬЭЬЭЬЭЬЭЬШ>LШЧЪЭЬЭЬШLЪЭЬЭЬЭЬЭЬЭЬЭЬШ7ЧЪЭЬбШLЪ ЭЬЭЬЭЬЭЬЭЬШ>LШЧЪЭЬЭЬШLЪЭЬЭЬЭЬЭЬЭЬЭЬШ7ЧЪЭЬбШLЪ ЭЬЭЬЭЬЭЬЭЬШ>ЧLЪЭЬШ LЪЬШLШЧЪЮЬЭЬЭЬЭЬЭЬЭЬЩ9ЪЭЬЭЬЩLЪ ЭЬЭЬЭЬЭЬб'ЩALЪЬЭЬЭЬЭЬЭЬЭЬб'ЩKLЪЬЭЬЭЬЩLЪЬЩЧЪЭЬЩЪЭЬЭЬЩLЪ ЭЬЭЬЭЬЭЬб(ЩBLЪЭЬЭЬЩLЪЬЩЧЪЭЬЩLЪЭЬЭЬЭЬЩLЪЭЬЭЬЭЬЩLЪЬЩЧLЪЬЩЧLЪЬбЩLЩЧЪЮЬЭЬЭЬЭЬЭЬЭЬШ?LЪЬЭЬШЪЭЬЭЬШLЪ ЭЬЭЬЭЬЭЬб(ШJLЪЬЭЬЭЬЭЬЭЬЭЬб'ШRLЪЬЭЬЭЬШ"LЪЬШЧЪЭЬШЪЭЬЭЬШLЪ ЭЬЭЬЭЬЭЬб)ШKLЪЭЬЭЬШLЪЬШЧЪЭЬШLЪЭЬЭЬЭЬШ/LЪЭЬЭЬЭЬШ5LЪЬШЧLЪЬШЧLЪЬб ЬШLШЧЪ Ю ЬЭЬЭЬЭЬЭЬШ>LЪЭЬШLЪЬШLШЧLШЪЮЬЭЬЭЬЭЬШLLЪЬШLЪЬШLЪЬШLЪЬШЧЪЮЬЭЬЭЬЭЬШLLЪЬШLЪЬШLЪЬШLЪЬШЧЪЮЬЭЬЭЬЭЬШLLЪЬШLЪЬШLЪЬШLЪЬШLШЧLШЪЮЬЭЬЭЬ ЭЬШGLЪЬШLЪЬШLЪЬШЧЪЮЬЭЬЭЬ ЭЬШGLЪЬШLЪЬШLЪЬШЧЪЮЬЭЬЭЬ ЭЬШGLЪЬШLЪЬШLЪЬШLШЧLШLЪЭЬШLШLШLШЧЕЧ╝(