code: def shear_matrix(angle, direction, point, normal):
      
      normal = unit_vector(normal[:3])
      direction = unit_vector(direction[:3])
      if abs(np.dot(normal, direction)) > 1e-6:
          raise ValueError("direction and normal vectors are not orthogonal")
      angle = math.tan(angle)
      M = np.identity(4)
      M[:3, :3] += angle * np.outer(direction, normal)
      M[:3, 3] = -angle * np.dot(point[:3], normal) * direction
      return M
docs: Return matrix to shear by angle along direction vector on shear plane.
  
      The shear plane is defined by a point and normal vector. The direction
      vector must be orthogonal to the plane's normal vector.
  
      A point P is transformed by the shear matrix into P" such that
      the vector P-P" is parallel to the direction vector and its extent is
      given by the angle of P-P'-P", where P' is the orthogonal projection
      of P onto the shear plane.
  
      >>> angle = (random.random() - 0.5) * 4*math.pi
      >>> direct = np.random.random(3) - 0.5
      >>> point = np.random.random(3) - 0.5
      >>> normal = np.cross(direct, np.random.random(3))
      >>> S = shear_matrix(angle, direct, point, normal)
      >>> np.allclose(1, np.linalg.det(S))
      True