function A = ConstructRotArbLine(L, theta)
  % -----------------------------------------------------------------------
  % ConstructRotArbLine - construct the 4 x 4 transformation matrix which
  %                       rotates objects about line L by angle theta.
  %
  % Args: L - the 4 x 2 matrix of endpoints of the line to rotate about
  %       theta - the angle by which to rotate about L
  %
  % Rtns: A - the 4 x 4 matrix of the transform
  % -----------------------------------------------------------------------

  d = L(:,2) - L(:,1); x = d(1); y = d(2); z = d(3);

  To = [ % translate first of L point to origin
     1    0    0  -L(1,1)
     0    1    0  -L(2,1)
     0    0    1  -L(3,1)
     0    0    0    1
  ];
  Toi = [ % inverse of To
     1    0    0   L(1,1)
     0    1    0   L(2,1)
     0    0    1   L(3,1)
     0    0    0    1
  ];

  h = sqrt(y*y + z*z);
  Rx = [ % rotate L about x into the zx plane
     1    0    0    0
     0   z/h -y/h   0
     0   y/h  z/h   0 
     0    0    0    1
  ];
  Rxi = Rx'; % inverse of Rx

  r = sqrt(x*x + y*y + z*z);
  Ry = [ % rotate L about y into the z axis
    h/r   0  -x/r   0
     0    1    0    0
    x/r   0   h/r   0 
     0    0    0    1
  ];
  Ryi = Ry'; % inverse of Ry

  c = cos(theta); s = sin(theta);
  Rz = [ % rotate about z by theta
     c   -s    0    0
     s    c    0    0 
     0    0    1    0
     0    0    0    1
  ];

  A = Toi * Rxi * Ryi * Rz * Ry * Rx * To; % the composite transform
end
