MGLaserBeamShortestDistanceToInclinedEllipse.html (MotionGenesis input/output).
(1) % MotionGenesis file: MGLaserBeamShortestDistanceToInclinedEllipse.txt
(2) % Copyright (c) 2019 Motion Genesis LLC. All rights reserved.
(3) % Note: Ellipse equation: x^2/a^2 + y^2/b^2 + z^2/c^2 = 1
(4) % For this ellipse: a = 2, b = 1, c = 1.
(5) %------------------------------------------------------------
(6) RigidFrame N % Horizontal plane.
(7) RigidFrame A % Laser (also creates point Ao).
(8) RigidFrame B % Inclined plane (also creates point Bo).
(9) RigidFrame C % Inclined ellipse (also creates point Co).
(10) Point CB( C ) % Point of C in contact with B.
(11) Point Q( C ) % Point of C closest to Ao.
(12) %------------------------------------------------------------
(13) % Rotation matrix relating inclined plane and horizontal plane.
(14) B.SetRotationMatrixZ( N, 30 * ConvertUnits(deg, rad) )
-> (15) B_N = [0.8660254, 0.5, 0; -0.5, 0.8660254, 0; 0, 0, 1]
(16) %------------------------------------------------------------
(17) % Known position vectors.
(18) Ao.SetPosition( No, Nx> + 2*Ny> )
-> (19) p_No_Ao> = Nx> + 2*Ny>
(20) Bo.SetPosition( No, 4*Nx> + Nz> )
-> (21) p_No_Bo> = 4*Nx> + Nz>
(22) CB.SetPosition( Bo, 2*Bx> )
-> (23) p_Bo_CB> = 2*Bx>
(24) Co.SetPosition( CB, 1*By> )
-> (25) p_CB_Co> = By>
(26) %------------------------------------------------------------
(27) % Q's position vector from Co unknown.
(28) Variable Lambda, x, y, z
(29) Q.SetPosition( Co, x*Bx> + y*By> + z*Bz> )
-> (30) p_Co_Q> = x*Bx> + y*By> + z*Bz>
(31) %------------------------------------------------------------
(32) % Form the equation for the ellipse's surface and its gradient.
(33) % Note: The gradient's direction is outward normal to the ellipse's surface.
(34) EllipseSurface = x^2/2^2 + y^2 + z^2 - 1
-> (35) EllipseSurface = -1 + y^2 + z^2 + 0.25*x^2
(36) Gradient_Surface> = D(EllipseSurface, x) * Bx> + D(EllipseSurface, y) * By> + D(EllipseSurface, z) * Bz>
-> (37) Gradient_Surface> = 0.5*x*Bx> + 2*y*By> + 2*z*Bz>
(38) %------------------------------------------------------------
(39) % Q's position vector from Ao is opposite the ellipse's outward normal at Q.
(40) % This can be proved by Calculus/minimization.
(41) Q.SetPosition( Ao, -Lambda * Gradient_Surface> )
-> (42) p_Ao_Q> = -0.5*Lambda*x*Bx> - 2*Lambda*y*By> - 2*Lambda*z*Bz>
(43) %------------------------------------------------------------
(44) % Solution must satisfy position vector loop equation.
(45) Loop> = P_No_Ao> + P_Ao_Q> + P_Q_Co> + P_Co_CB> + P_CB_Bo> + P_Bo_No>
-> (46) Loop> = (-2-x-0.5*Lambda*x)*Bx> + (-1-y-2*Lambda*y)*By> - z*(1+2*Lambda)*Bz>
- 3*Nx> + 2*Ny> - Nz>
(47) Loop = Matrix( N, Loop> ) % Dot-product Nx>, Ny>, Nz> with Loop>.
-> (48) Loop[1] = -4.232051 + 0.5*y + Lambda*y - 0.8660254*x - 0.4330127*Lambda*x
-> (49) Loop[2] = 0.1339746 - 0.8660254*y - 0.5*x - 1.732051*Lambda*y - 0.25*
Lambda*x
-> (50) Loop[3] = -1 - z*(1+2*Lambda)
(51) %------------------------------------------------------------
(52) % There are 4 unknowns: x, y, z, Lambda.
(53) % There are 4 equations: EllipseSurface = 0, Loop[i] = 0 (i=1,2,3).
(54) % Step 1: Solve for x, y, z in terms of Lambda.
(55) % Step 2 (optional): Make EllipseSurface explicit in Lambda.
(56) Solve( Loop = 0, x, y, z )
-> (57) x = -7.196152/(2+Lambda)
-> (58) y = 2.232051/(1+2*Lambda)
-> (59) z = -1/(1+2*Lambda)
(60) Explicit( EllipseSurface, Lambda )
-> (61) EllipseSurface = -1 + 12.94615/(2+Lambda)^2 + 5.982051/(1+2*Lambda)^2
(62) %------------------------------------------------------------
(63) % Math digression: One way to solve for Lambda is to rearrange
(64) % EllipseSurface into a 4th-order polynomial equation and solve.
(65) PolynomialEqn = Expand( Rhs(EllipseSurface) * (2+Lambda)^2 * (1+2*Lambda)^2 )
-> (66) PolynomialEqn = 5.982051*(2+Lambda)^2 + 12.94615*(1+2*Lambda)^2 - (2+
Lambda)^2*(1+2*Lambda)^2
(67) Expand( PolynomialEqn, 0:2 )
-> (68) PolynomialEqn = 32.87436 + 55.71281*Lambda + 24.76666*Lambda^2 - 20*Lambda^3
- 4*Lambda^4
(69) roots = GetRoots( PolynomialEqn = 0, Lambda, 4 )
-> (70) roots = [-5.701827; -0.6912399 - 0.4623138*imaginary; -0.6912399 + 0.4623138
*imaginary; 2.084307]
(71) DistanceAoToQmin = EvaluateToNumber( Q.GetDistance(Ao), Lambda = roots[1] )
-> (72) DistanceAoToQmin = 6.156405
(73) DistanceAoToQmax = EvaluateToNumber( Q.GetDistance(Ao), Lambda = roots[4] )
-> (74) DistanceAoToQmax = 2.694949
(75) %------------------------------------------------------------
(76) % Another way to solve for Lambda is with nonlinear solver which needs a guess for Lambda.
(77) % The laser hits the ellipses twice (near side which is minimum distance and far side
(78) % which is maximum distance). To solve for minimum distance, guess a small positive Lambda.
(79) Solve( EllipseSurface = 0, Lambda = 0.1 )
-> (80) Lambda = 2.084307
(81) DistanceBetweenQAndAo = EvaluateToNumber( Q.GetDistance(Ao) )
-> (82) DistanceBetweenQAndAo = 2.694949
(83) %------------------------------------------------------------
(84) % Determine the unit vector u> in the direction of Q's position from Ao.
(85) % Express it in terms of both Bx>, By>, Bz> and Nx>, Ny>, Nz>.
(86) uB> = EvaluateToNumber( Q.GetUnitPositionVector(Ao) )
-> (87) uB> = 0.6813389*Bx> - 0.6679919*By> + 0.2992727*Bz>
(88) uN> = Express( uB>, N )
-> (89) uN> = 0.9240527*Nx> - 0.2378285*Ny> + 0.2992727*Nz>
(90) %------------------------------------------------------------
(91) % Form the gradient of s (derivative of s with respect to Q's position vector from Ao).
(92) % Use the definition, and express the results in N. Verify that it is opposite uN>.
(93) Variable xN, yN, zN
(94) Ao.SetPosition( No, xN * Nx> + yN * Ny> + zN * Nz> )
-> (95) p_No_Ao> = xN*Nx> + yN*Ny> + zN*Nz>
(96) r> = P_Ao_No> + P_No_Bo> + P_Bo_CB> + P_CB_Co> + P_Co_Q>
-> (97) r> = (2+x)*Bx> + (1+y)*By> + z*Bz> + (4-xN)*Nx> - yN*Ny> + (1-zN)*Nz>
(98) s = GetMagnitude( r> )
-> (99) s = sqrt(yN^2+z^2+(1+y)^2+(2+x)^2+(1-zN)^2+(4-xN)^2+2*z*(1-zN)+1.732051
*(2+x)*(4-xN)-x*yN-1.732051*yN*(2.154701+y)-(1+y)*(4-xN))
(100) Gradient_s_N> = EvaluateToNumber( D(s, xN) * Nx> + D(s, yN) * Ny> + D(s, zN) * Nz>, &
xN = 1, yN = 2, zN = 0 )
-> (101) Gradient_s_N> = -0.9240527*Nx> + 0.2378285*Ny> - 0.2992727*Nz>
(102) unitGradient_s_N> = GetUnitVector( Gradient_s_N> )
-> (103) unitGradient_s_N> = -0.9240527*Nx> + 0.2378285*Ny> - 0.2992727*Nz>
(104) ShouldBeZeroN> = uN> + unitGradient_s_N>
-> (105) ShouldBeZeroN> = 0>
(106) %------------------------------------------------------------
(107) % To show the previous result is basis-independent, again calculate the
(108) % gradient of s, but instead use Bx>, By>, Bz>. Verify it is opposite uB>.
(109) Variable xB, yB, zB
(110) Ao.SetPosition( No, xB * Bx> + yB * By> + zB * Bz> )
-> (111) p_No_Ao> = xB*Bx> + yB*By> + zB*Bz>
(112) r> := P_Ao_No> + P_No_Bo> + P_Bo_CB> + P_CB_Co> + P_Co_Q>
-> (113) r> = (2+x-xB)*Bx> + (1+y-yB)*By> + (z-zB)*Bz> + 4*Nx> + Nz>
(114) s := GetMagnitude( r> )
-> (115) s = sqrt(26.85641+2*z+4*yB+6.928203*x+(z-zB)^2+(1+y-yB)^2+(2+x-xB)^2-6.928203
*xB-4*y-2*zB)
(116) Gradient_s_B> = EvaluateToNumber( D(s, xB) * Bx> + D(s, yB) * By> + D(s, zB) * Bz>, &
xB = 1*cosDegrees(30) + 2*sinDegrees(30), yB = -1*sinDegrees(30) + 2*cosDegrees(30), zB = 0 )
-> (117) Gradient_s_B> = -0.6813389*Bx> + 0.6679919*By> - 0.2992727*Bz>
(118) unitGradient_s_B> = GetUnitVector( Gradient_s_B> )
-> (119) unitGradient_s_B> = -0.6813389*Bx> + 0.6679919*By> - 0.2992727*Bz>
(120) ShouldBeZeroB> = uB> + unitGradient_s_B>
-> (121) ShouldBeZeroB> = 0>
(122) %------------------------------------------------------------
(123) % Save input and program responses.
Saved by Motion Genesis LLC. Command names and syntax: Copyright (c) 2009-2021 Motion Genesis LLC. All rights reserved.