# Example 21.1
ixx = (1/12) * 2 kg * (0.2 m)^2 + 2 kg * (0.1 m)^2 \
    + 2 kg * (0.2 m)^2 \
    + (1/12) * (4 kg) * (0.4 m)^2 + 4 kg * ((0.2 m)^2 + (0.2 m)^2)
iyy = (1/12) * (2 kg) * (0.2 m)^2 + 2 kg * (0.1 m)^2 \
    + (1/12) * (2 kg) * (0.2 m)^2 + 2 kg * ((-0.1 m)^2 + (0.2 m)^2) \
    + 4 kg * ((-0.2 m)^2 + (0.2 m)^2)
izz = (1/12) * 2 kg * (0.2 m)^2 + 2 kg * (0.1 m)^2 \
    + (1/12) * 4 kg * (0.4 m)^2 + 4 kg * ((-0.2 m)^2 + (0.2 m)^2)
ixy = 4 kg * (-0.2 m) * (0.2 m)
iyz = 4 kg * (0.2 m) * (0.2 m)
izx = 2 kg * (0.2 m) * (-0.1 m) + 4 kg * (0.2 m) * (-0.2 m)
rd = [-0.2, 0.4, 0.2]
u = rd / mag(rd)
iaa = ixx * u[1]^2 + iyy * u[2]^2 + izz * u[3]^2 \
        - 2 * ixy * u[1] * u[2] \
        - 2 * iyz * u[2] * u[3] \
        - 2 * izx * u[3] * u[1] 

# Example 21.2
FPS
w = 1.5 lb/ft
len = 3 ft
w0 = len * w
m0 = w0 / grav()
vg1[3] = -10 ft/s
rga = [-0.667, 0.5] m
ixx = 0.0272 slug ft^2
iyy = 0.0155 slug ft^2
izz = 0.0427 slug ft^2
begin
 v = 1 m/s
 omx = 1 rad/s
 omy = 1 rad/s
 omz = 1 rad/s
 om = [omx, omy, omz]
 vg2[3] = -v
 hg2 = [ixx * omx, iyy * omy, izz * omz]
 MM = cross(rga, m0 * vg1) - (cross(rga, m0 * vg2) + hg2)
 ee = vg2 - cross(om, rga)
end
solve(MM[1], MM[2], MM[3], ee[3], v, omx, omy, omz)

# Example 21.3
MKS(N,  J)
tt = 5 N m
m0 = 10 kg
work = tt * 4 * pi() rad
r1 = 0.1 m
r2 = 0.3 m
ix = 0.5 * m0 * r1^2
iy = 0.25 * m0 * r1^2 + m0 * r2^2
iz = iy
begin
  omcd = 1 rad/s
  omce = omcd * (r2 / r1)
  oma = [-omce, 0, omcd]
  KE = 0.5 * (ix * oma[1]^2 + iy * oma[2]^2 + iz * oma[3]^2)
  err = KE - work
end
solve(err, omcd)

# Example 21.4
m0 = 10 kg
q = 10 deg
iz = 0.1 kg m^2
ix = 0.05 kg m^2
iy = ix
om = 30 rad/s
omx = 0 rad/s
omy = - om * sin(q)
omz = om * cos(q)
begin
  Ax = 0 N
  Ay = 70 N
  Bx = 0 N
  By = 26 N
  Mx = -Ay * (0.2 m) + By * (0.25 m) -( -(iy - iz) * omy * omz)
  My =  Ax * 0.2 m * cos(q) - Bx * 0.25 m * cos(q)
  Mz = Ax * 0.2 m * sin(q) - Bx * 0.25 m * sin(q) 
  ee = Ax
end
solve(Mx, My, Mz, ee, Ax, Ay, Bx, By)

# Example 21.6
m0 = 10 kg
oms = 6 rad/s
omp = 3 rad/s
OM[3] = omp
r = 0.2 m
ix = 0.25 * m0 * r^2
iz = ix
iy = 0.5 * m0 * r^2
omx = 0 rad/s
omy = oms
omz = omp
agy = -0.5 m * omp^2
begin
ax = 1 N
ay = 1 N
az = 1 N
bx = 1 N
bz = 1 N
Mx = -az * 0.5 m + bz * 0.5 * m -( iy * OM[3] * omy)
Mz = ax * 0.5 m - bx * 0.5 m
fx = ax + bx
fy = ay - m0 * agy
fz = az + bz - m0 * grav()
end
solve(Mx, Mz, fx, fy, fz, ax, ay, az, bx, bz)

# Example 21.7
m0 = 0.5 kg
q = 60 deg
ws = 100 rad/s
iz = 0.45e-3 kg m^2
ii = 1.20e-3 kg m^2
x = 0.05 m
begin
  phidot = 100 rad/s
  mx = m0 * grav() * x * sin(q) - (-ii * phidot^2 * sin(q) * cos(q) \
   + iz * phidot * sin(q) * (phidot * cos(q) + ws))
end
solve(mx, phidot)

# Example 21.8
m0 = 1 kg
omd = 70 rad/s
mb = 2 kg
omp = 0.5 rad/s
r = 0.05 m
iz = 0.5 * m0 * r^2
x = 0.2 m
begin
ss = 1 mm
f = 1 N
mx = m0 * grav() * ss - f * ss - (omp * iz * (-omd))
my = mb * grav() * ss - f * x
end
solve(mx, my, ss, f)

# Example 21.9
q = 30 deg
phidot = 3 rad/s
ii = 1 kg/m^2
iz = (1/3) * ii
psidot = (ii - iz) / iz * phidot * cos(q)
omx = 0 rad/s
omy = phidot * ii * sin(q) / ii
omz = phidot * ii * cos(q) / iz
sqrt(omx^2 + omy^2 + omz^2)

# Homework Problems:

TBD