MGBeamOnTwoCablesDynamicsViaKaneEmbedConstraintsAndLinearize.html  (MotionGenesis input/output).
```   (1) %    File: MGBeamOnTwoCablesDynamicsViaKaneEmbedConstraintsAndLinearize.txt
(2) % Problem: Dynamics of a beam supported by two cables.
(3) %--------------------------------------------------------------------
(4) NewtonianFrame N         % Ceiling with Nx> horizontally right, Ny> vertically down
(5) RigidBody      B         % Beam with Bx> pointing from Bo to Bc and Bz> = Nz>.
(6) Point          Nc( N )   % Point of N attached to cable C.
(7) Point          Bc( B )   % Point of B attached to cable C.
(8) %--------------------------------------------------------------------
(9) Constant  LN = 6 m       % Distance between No and NC.
(10) Constant  LB = 4 m       % Distance between Bo and BC.
(11) Constant  LA = 2.7 m     % Length of cable A (connects Bo to No).
(12) Constant  LC = 3.7 m     % Length of cable C (connects BC to NC).
(13) Constant  g = 9.8 m/s^2  % Earth's gravitational acceleration.
(14) Constant  b = 230 N*s/m  % Aerodynamic damping constant.
(15) B.SetMass( m = 120 kg )
(16) B.SetInertia( Bcm, Izz = m*LB^2/12, 0, Izz )
-> (17) Izz = 0.08333333*m*LB^2

(18) %--------------------------------------------------------------------
(19) Variable  x'', y''       % Nx> and Ny> measures of Bo's position vector from No.
(20) Variable  q''            % Angle from Nx> to Bx> with positive sense +Nz>.
(21) %--------------------------------------------------------------------
(22) %   Rotational kinematics relating  Bx>, By>, Bz>  to  Nx>, Ny>, Nz>.
(23) B.RotateZ( N,  q )
-> (24) B_N = [cos(q), sin(q), 0;  -sin(q), cos(q), 0;  0, 0, 1]
-> (25) w_B_N> = q'*Bz>
-> (26) alf_B_N> = q''*Bz>

(27) %--------------------------------------------------------------------
(28) %   Translational kinematics.
(29) Bo.Translate( No,  x*Nx> + y*Ny> )
-> (30) p_No_Bo> = x*Nx> + y*Ny>
-> (31) v_Bo_N> = x'*Nx> + y'*Ny>
-> (32) a_Bo_N> = x''*Nx> + y''*Ny>

(33) Bcm.Translate( Bo,  0.5*LB*Bx> )
-> (34) p_Bo_Bcm> = 0.5*LB*Bx>
-> (35) v_Bcm_N> = 0.5*LB*q'*By> + x'*Nx> + y'*Ny>
-> (36) a_Bcm_N> = -0.5*LB*q'^2*Bx> + 0.5*LB*q''*By> + x''*Nx> + y''*Ny>

(37) BC.SetPosition( Bo,  LB*Bx> )
-> (38) p_Bo_Bc> = LB*Bx>

(39) NC.SetPosition( No,  LN*Nx> )
-> (40) p_No_Nc> = LN*Nx>

(41) %--------------------------------------------------------------------
(42) %   Constraint relating cable lengths to magnitudes of position vectors.
(43) LengthConstraint = LA^2 - Bo.GetDistanceSquared( No )
-> (44) LengthConstraint = LA^2 - x^2 - y^2

(45) LengthConstraint = LC^2 - Bc.GetDistanceSquared( Nc )
-> (46) LengthConstraint = LC^2 + 2*LB*cos(q)*(LN-x) - LB^2 - y^2 - 2*LB*y*sin(q)
- (LN-x)^2

(47) %--------------------------------------------------------------------
(48) %   Relevant forces on beam (to calculate generalized forces).
(49) %   Note: Cable tensions TA and TC do not contributes to generalized force
(50) %   because both LengthConstraint and LengthConstraint are embedded.
(51) Bcm.AddForce( m*g*Ny> - b * Bcm.GetVelocity(N) )
-> (52) Force_Bcm> = -0.5*b*LB*q'*By> - b*x'*Nx> + (m*g-b*y')*Ny>

(53) %--------------------------------------------------------------------
(54) %   Provide initial value for x here (to solve for initial values of y, q).
(55) Input  x = 2 m,  x' = 0 m/s
(56) %--------------------------------------------------------------------
(57) %   Use length constraints to set initial values for y and q (with guess).
(58) %   Differentiate length constraints to solve for y' and q' in terms of x'.
(59) %   Differentiate again and solve for y'' and q'' in terms of x''.
(60) SolveSetInputDt( LengthConstraint = 0,   y = 3 m,  q = 15 deg )

->   %  INPUT has been assigned as follows:
->   %   y                         1.813835714721706       m
->   %   q                         27.67399398437716       deg

-> (61) y' = -x*x'/y
-> (62) q' = -(LB*cos(q)-LN-LB*x*sin(q)/y)*x'/(LB*(y*cos(q)+sin(q)*(LN-x)))
-> (63) y'' = -(x'^2+y'^2+x*x'')/y
-> (64) q'' = (2*LB*sin(q)*q'*x'+LB*q'*(y*sin(q)*q'-2*cos(q)*y'-cos(q)*(LN-x)*
q')+(LN-LB*cos(q))*x''+LB*sin(q)*(x'^2+y'^2+x*x'')/y)/(LB*(y*cos(q)+sin
(q)*(LN-x)))

(65) %--------------------------------------------------------------------
(66) %   Dynamics via Kane's embedded method (y' and q' are in terms of x').
(67) %   This method embeds both LengthConstraint and LengthConstraint.
(68) SetGeneralizedSpeeds( x' )
(69) Dynamics = System.GetDynamicsKane()
-> (70) Dynamics = b*x' + 0.5*x*(2*m*g+2*b*x*x'/y-b*LB*cos(q)*q')/y + 0.25*(
LB*cos(q)-LN-LB*x*sin(q)/y)*(2*b*sin(q)*x'+2*cos(q)*(m*g+b*x*x'/y)-b*
LB*q')/(y*cos(q)+sin(q)*(LN-x)) + 0.25*m*(2*x*(2*(1+x^2/y^2)*x'^2/y+LB*
sin(q)*q'^2-LB*cos(q)*(sin(q)*(1+x^2/y^2)*x'^2/y+4*sin(q)*x'*q'+q'*(2*x
*cos(q)*x'/y+2*y*sin(q)*q'-cos(q)*(LN-x)*q'))/(y*cos(q)+sin(q)*(LN-x)))
/y-2*LB*cos(q)*q'^2-(2*LB*sin(q)*(sin(q)*(1+x^2/y^2)*x'^2/y+2*sin(q)*
x'*q'-(LN*cos(q)-y*sin(q))*q'^2)-(LB*cos(q)-LN-LB*x*sin(q)/y)*(2*cos(q)
*(1+x^2/y^2)*x'^2/y-LB*(sin(q)*(1+x^2/y^2)*x'^2/y+2*sin(q)*x'*q'+q'*(2*
x*cos(q)*x'/y+y*sin(q)*q'-cos(q)*(LN-x)*q'))/(y*cos(q)+sin(q)*(LN-x))))
/(y*cos(q)+sin(q)*(LN-x))) + 0.25*(4*Izz*(LB*cos(q)-LN-LB*x*sin(q)/y)^2
/(LB^2*(y*cos(q)+sin(q)*(LN-x))^2)+m*(4+(LB*cos(q)-LN-LB*x*sin(q)/y)^2/
(y*cos(q)+sin(q)*(LN-x))^2+4*sin(q)*(LB*cos(q)-LN-LB*x*sin(q)/y)/(y*cos
(q)+sin(q)*(LN-x))+4*x*(x/y+cos(q)*(LB*cos(q)-LN-LB*x*sin(q)/y)/(y*cos(
q)+sin(q)*(LN-x)))/y))*x'' - 0.5*b*LB*sin(q)*q' - Izz*(LB*cos(q)-LN-LB*
x*sin(q)/y)*(sin(q)*(1+x^2/y^2)*x'^2/y+2*sin(q)*x'*q'+q'*(2*x*cos(q)*
x'/y+y*sin(q)*q'-cos(q)*(LN-x)*q'))/(LB*(y*cos(q)+sin(q)*(LN-x))^2)

(71) %********************************************************************
(72) %   LINEARIZATION AND FREQUENCY ANALYSIS
(73) %--------------------------------------------------------------------
(74) Variable  dx'', dy'', dq''            % Perturbations of x, y, q.
(75) Constant  xStatic, yStatic, qStatic   % Static (nominal) values of x, y, q.
(76) %--------------------------------------------------------------------
(77) %   Find static (nominal) solution for x, y, q.
(78) Nominal = Evaluate( [Dynamics; LengthConstraint],  x = xStatic, y = yStatic, q = qStatic, x' = 0,  x'' = 0 )
-> (79) Nominal = 0.5*m*g*(2*xStatic/yStatic+cos(qStatic)*(LB*cos(qStatic)-
LN-LB*xStatic*sin(qStatic)/yStatic)/(yStatic*cos(qStatic)+sin(qStatic)*
(LN-xStatic)))

-> (80) Nominal = LA^2 - xStatic^2 - yStatic^2
-> (81) Nominal = LC^2 + 2*LB*cos(qStatic)*(LN-xStatic) - LB^2 - yStatic^2
- 2*LB*yStatic*sin(qStatic) - (LN-xStatic)^2

(82) Solve( Nominal = 0,   xStatic = 1 m,  yStatic = 3 m,  qStatic = 15 deg )
-> (83) xStatic = 0.8153669
-> (84) yStatic = 2.573942
-> (85) qStatic = 0.2257464       %  or  qStatic = 12.93432 deg.

(86) %--------------------------------------------------------------------
(87) %   Linearize constraints about static (nominal) solution.
(88) %   Solve for dy and dq in terms of dx.
(89) LinearizedConstraint = Linearize1( LengthConstraint,  x = xStatic : dx,  y = yStatic : dy,  q = qStatic : dq )
-> (90) LinearizedConstraint = -2*xStatic*dx - 2*yStatic*dy
-> (91) LinearizedConstraint = 2*(LN-xStatic-LB*cos(qStatic))*dx - 2*(yStat
ic+LB*sin(qStatic))*dy - 2*LB*(yStatic*cos(qStatic)+sin(qStatic)*(LN-
xStatic))*dq

(92) SolveDt( LinearizedConstraint = 0,  dy,  dq )
-> (93) dy = -xStatic*dx/yStatic
-> (94) dq = -(LB*cos(qStatic)-LN-LB*xStatic*sin(qStatic)/yStatic)*dx/(LB*(ySt
atic*cos(qStatic)+sin(qStatic)*(LN-xStatic)))

-> (95) dy' = -xStatic*dx'/yStatic
-> (96) dq' = -(LB*cos(qStatic)-LN-LB*xStatic*sin(qStatic)/yStatic)*dx'/(LB*(
yStatic*cos(qStatic)+sin(qStatic)*(LN-xStatic)))

-> (97) dy'' = -xStatic*dx''/yStatic
-> (98) dq'' = -(LB*cos(qStatic)-LN-LB*xStatic*sin(qStatic)/yStatic)*dx''/(LB*(
yStatic*cos(qStatic)+sin(qStatic)*(LN-xStatic)))

(99) %--------------------------------------------------------------------
(100) %   Linearize dynamics about nominal solution.
(101) SetAutoZee( OFF )
(102) LinearizedDynamics = Linearize1( Dynamics,  x = xStatic : dx,  x' = 0 : dx',  x'' = 0 : dx'',  &
y = yStatic : dy,  q = qStatic : dq )
-> (103) LinearizedDynamics = 0.5*m*g*(LB*xStatic*sin(qStatic)*cos(qStatic)/(yStatic^2
*(yStatic*cos(qStatic)+sin(qStatic)*(LN-xStatic)))-2*xStatic/yStatic^2
-cos(qStatic)^2*(LB*cos(qStatic)-LN-LB*xStatic*sin(qStatic)/yStatic)/(
yStatic*cos(qStatic)+sin(qStatic)*(LN-xStatic))^2)*dy + 0.5*m*g*(2/yStatic
+2*xStatic^2/yStatic^3+sin(qStatic)*(LB*cos(qStatic)-LN-LB*xStatic*sin
(qStatic)/yStatic)^2/(LB*(yStatic*cos(qStatic)+sin(qStatic)*(LN-xStat
ic))^2)+cos(qStatic)*(LB*cos(qStatic)-LN-LB*xStatic*sin(qStatic)/yStatic)
*(sin(qStatic)*(LN*yStatic+LB*xStatic*sin(qStatic)-LB*yStatic*cos(qSt
atic))+cos(qStatic)*(LN-xStatic)*(LB*cos(qStatic)-LN-LB*xStatic*sin(
qStatic)/yStatic))/(LB*(yStatic*cos(qStatic)+sin(qStatic)*(LN-xStatic))^3)
-cos(qStatic)*(LB*sin(qStatic)*(xStatic^2+yStatic^2)/yStatic^3-2*sin(
qStatic)*(LB*cos(qStatic)-LN-LB*xStatic*sin(qStatic)/yStatic)/(yStatic
*cos(qStatic)+sin(qStatic)*(LN-xStatic))-2*xStatic*cos(qStatic)*(LB*
cos(qStatic)-LN-LB*xStatic*sin(qStatic)/yStatic)/(yStatic*(yStatic*cos
(qStatic)+sin(qStatic)*(LN-xStatic))))/(yStatic*cos(qStatic)+sin(qSta
tic)*(LN-xStatic)))*dx + 0.25*b*(4+4*xStatic*(xStatic/yStatic+cos(qSt
atic)*(LB*cos(qStatic)-LN-LB*xStatic*sin(qStatic)/yStatic)/(yStatic*
cos(qStatic)+sin(qStatic)*(LN-xStatic)))/yStatic+(LB*cos(qStatic)-LN-
LB*xStatic*sin(qStatic)/yStatic)*(4*sin(qStatic)+(LB*cos(qStatic)-LN-
LB*xStatic*sin(qStatic)/yStatic)/(yStatic*cos(qStatic)+sin(qStatic)*(
LN-xStatic)))/(yStatic*cos(qStatic)+sin(qStatic)*(LN-xStatic)))*dx'
+ 0.25*(m*(4+(LB*cos(qStatic)-LN-LB*xStatic*sin(qStatic)/yStatic)^2/(
yStatic*cos(qStatic)+sin(qStatic)*(LN-xStatic))^2+4*sin(qStatic)*(LB*
cos(qStatic)-LN-LB*xStatic*sin(qStatic)/yStatic)/(yStatic*cos(qStatic)
+sin(qStatic)*(LN-xStatic))+4*xStatic*(xStatic/yStatic+cos(qStatic)*(
LB*cos(qStatic)-LN-LB*xStatic*sin(qStatic)/yStatic)/(yStatic*cos(qSta
tic)+sin(qStatic)*(LN-xStatic)))/yStatic)+4*(LB*cos(qStatic)-LN-LB*xS
tatic*sin(qStatic)/yStatic)^2*Izz/(LB^2*(yStatic*cos(qStatic)+sin(qSt
atic)*(LN-xStatic))^2))*dx'' - 0.5*m*g*(LB*cos(qStatic)*(2*sin(qStatic)
+xStatic*cos(qStatic)/yStatic)*(yStatic*cos(qStatic)+sin(qStatic)*(LN-
xStatic))-sin(qStatic)*(LN+LB*xStatic*sin(qStatic)/yStatic)*(yStatic*
cos(qStatic)+sin(qStatic)*(LN-xStatic))-cos(qStatic)*(yStatic*sin(qSt
atic)-cos(qStatic)*(LN-xStatic))*(LB*cos(qStatic)-LN-LB*xStatic*sin(
qStatic)/yStatic))*dq/(yStatic*cos(qStatic)+sin(qStatic)*(LN-xStatic))^2

(104) %--------------------------------------------------------------------
(105) %   Evaluate coefficients to numerical values (except keep explicit in b)
(106) LinearizedDynamics := EvaluateToNumber( LinearizedDynamics, b = b )
-> (107) LinearizedDynamics = 465.3305*dx + 0.8597899*b*dx' + 107.4004*dx''

(108) mEffective = GetCoefficient( LinearizedDynamics,  dx'' )  % Effective mass.
-> (109) mEffective = 107.4004

(110) bEffective = GetCoefficient( LinearizedDynamics,  dx'  )  % Effective damping.
-> (111) bEffective = 0.8597899*b

(112) kEffective = GetCoefficient( LinearizedDynamics,  dx   )  % Effective stiffness.
-> (113) kEffective = 465.3305

(114) %--------------------------------------------------------------------
(115) %   Calculated undamped and damped natural frequencies (explicit in b).
(116) wn = sqrt( kEffective / mEffective )                      % Undamped natural frequency.
-> (117) wn = 2.081507

(118) TauPeriod = 2 * pi / wn                                   % Undamped period.
-> (119) TauPeriod = 3.018576

(120) zeta =  bEffective / (2 * sqrt(mEffective * kEffective))  % Damping ratio.
-> (121) zeta = 0.002236589*bEffective

(122) wd = wn * sqrt(1 - zeta^2)                                % Damped frequency.
-> (123) wd = 2.081507*sqrt(1-zeta^2)

(124) TauPeriodDamped = 2 * pi / wd                             % Damped period.
-> (125) TauPeriodDamped = 6.283185/wd

(126) %--------------------------------------------------------------------
(127) %   Save input and output responses.
```