MGDiskRollingDownInclinedPlaneEnergy.html   (MotionGenesis input/output).
```   (1) % MotionGenesis file:  MGDiskRollingDownInclinedPlaneEnergy.txt
(2) % Copyright (c) 2020 Motion Genesis LLC.
(3) %--------------------------------------------------------------------
(4) NewtonianFrame  N       % Inclined plane with Nx> downhill and Ny> perpendicular to the plane.
(5) RigidFrame      E       % Earth with Ex> horizontally-right, Ey> vertically-upward.
(6) RigidBody       B       % Cylindrical disk whose symmetry axis parallel to Nz>.
(7) Point           BN( B ) % Point of B in contact with inclined plane N.
(8) %--------------------------------------------------------------------
(9) Constant  r             % Cylinder's radius.
(10) Constant  g             % Earth's gravitational acceleration.
(11) Constant  psi           % Angle inclined plane makes with Earth's horizontal.
(12) Variable  theta'        % B's rotation angle from Ey> with sense -Nz>.
(13) Variable  x'            % Nx> measure of Bcm's position from No.
(14) Variable  Fx, Fy        % Friction and normal force on B at BN.
(15) B.SetMassInertia( m,  Ixx, Iyy, Izz )
(16) %--------------------------------------------------------------------
(17) %   Rotational kinematics.
(18) E.SetRotationMatrixZ( N,  psi )
-> (19) E_N = [cos(psi), sin(psi), 0;  -sin(psi), cos(psi), 0;  0, 0, 1]

(20) B.RotateNegativeZ( E,  theta )
-> (21) B_E = [cos(theta), -sin(theta), 0;  sin(theta), cos(theta), 0;  0, 0, 1]
-> (22) w_B_E> = -theta'*Bz>

(23) %--------------------------------------------------------------------
(24) %   Translational kinematics.
(25) Bcm.SetPositionVelocity( No,  x*Nx> + r*Ny> )
-> (26) p_No_Bcm> = x*Nx> + r*Ny>
-> (27) v_Bcm_N> = x'*Nx>

(28) BN.SetPositionVelocity( Bcm,  -r*Ny> )
-> (29) p_Bcm_BN> = -r*Ny>
-> (30) v_BN_N> = (x'-r*theta')*Nx>

(31) %--------------------------------------------------------------------
(32) %   Motion constraints (rolling at BN).
(33) RollingConstraint = Dot( BN.GetVelocity(N),  Nx> )
-> (34) RollingConstraint = x' - r*theta'

(35) Solve( RollingConstraint = 0,  x' )
-> (36) x' = r*theta'

(37) %--------------------------------------------------------------------
(38) %   Add relevant contact/distance forces on body B.
(39) B.AddForceGravity( -g*Ey> )      % Gravity force on B.
-> (40) Force_Bcm> = -m*g*Ey>

(41) BN.AddForce( Fx*Nx> + Fy*Ny> )   % Friction and normal force on B at BN.
-> (42) Force_BN> = Fx*Nx> + Fy*Ny>

(43) %--------------------------------------------------------------------
(44) %   Form B's kinetic energy in N.
(45) KE = B.GetKineticEnergy()        % Valid when B slides or rolls on N.
-> (46) KE = 0.5*Izz*theta'^2 + 0.5*m*x'^2

(47) Explicit( KE,  m, r, theta' )    % Valid when B rolls on N.
-> (48) KE = 0.5*(Izz+m*r^2)*theta'^2

(49) %--------------------------------------------------------------------
(50) %   Calculate a gravitational potential energy PEgravity.
(51) PEgravity = System.GetForceGravityPotentialEnergy( -g*Ey>, No )
-> (52) PEgravity = m*g*(r*cos(psi)-sin(psi)*x)

(53) %--------------------------------------------------------------------
(54) %   For this system, mechanical energy (kinetic energy plus potential energy) is conserved.
(55) MechanicalEnergy = Explicit( KE + PEgravity,   m, r, theta' )
-> (56) MechanicalEnergy = m*g*(r*cos(psi)-sin(psi)*x) + 0.5*(Izz+m*r^2)*theta'^2

(57) %--------------------------------------------------------------------
(58) %   Determine the initial value of MechanicalEnergy.
(59) InitialMechanicalEnergy = Evaluate( MechanicalEnergy,  x = 0,  theta' = 0 )
-> (60) InitialMechanicalEnergy = m*g*r*cos(psi)

(61) %--------------------------------------------------------------------
(62) %   Use the fact that MechanicalEnergy is constant to solve for theta'.
(63) eqnToSolve = Explicit( MechanicalEnergy - InitialMechanicalEnergy,  m, r, theta' )
-> (64) eqnToSolve = 0.5*(Izz+m*r^2)*theta'^2 - m*g*sin(psi)*x

(65) SolveQuadraticPositiveRoot( eqnToSolve = 0,  theta' )
-> (66) theta' = 1.414214*sqrt(m*g*sin(psi)*x/(Izz+m*r^2))

(67) %--------------------------------------------------------------------
(68) %   Determine w = theta' for various shapes of B (solid/hollow disk or sphere).
(69) wSolidDisk = Evaluate( theta',  Izz = 1/2*m*r^2 )
-> (70) wSolidDisk = 1.154701*sqrt(g*sin(psi)*x/r^2)

(71) wHollowDisk = Evaluate( theta',  Izz = m*r^2 )
-> (72) wHollowDisk = 1*sqrt(g*sin(psi)*x/r^2)

(73) wSolidSphere = Evaluate( theta',  Izz = 2/5*m*r^2 )
-> (74) wSolidSphere = 1.195229*sqrt(g*sin(psi)*x/r^2)

(75) wHollowSphere = Evaluate( theta',  Izz = 2/3*m*r^2 )
-> (76) wHollowSphere = 1.095445*sqrt(g*sin(psi)*x/r^2)

(77) %--------------------------------------------------------------------
(78) %   Form power on B (consider all forces -- e.g., gravity and normal/friction force).
(79) %   Note: Due to contact, Fy does not show up in power (normal force is a "workless" force).
(80) %   Note: When B rolls,   Fx does not show up in power (friction force is "workless" force when rolling).
(81) PowerGravity = Dot( Bcm.GetResultantForce(), Bcm.GetVelocity(N) )
-> (82) PowerGravity = m*g*sin(psi)*x'

(83) PowerContact = Dot(  BN.GetResultantForce(),  BN.GetVelocity(N) )
-> (84) PowerContact = Fx*(x'-r*theta')

(85) Explicit( PowerContact )
-> (86) PowerContact = 0

(87) Power = Explicit( PowerGravity + PowerContact,  x' )
-> (88) Power = m*g*sin(psi)*x'

(89) %--------------------------------------------------------------------
(90) %   Optional: Verify that Power = -Dt( PEgravity ).
(91) isEqualPower = IsSimplifyEqual( Power,  -Dt(PEgravity) )
-> (92) isEqualPower = true

(93) %--------------------------------------------------------------------
(94) %   Record input together with responses.
```