code: def shear_matrix(angle, direction, point, normal):
      """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 = numpy.random.random(3) - 0.5
      >>> point = numpy.random.random(3) - 0.5
      >>> normal = numpy.cross(direct, numpy.random.random(3))
      >>> S = shear_matrix(angle, direct, point, normal)
      >>> numpy.allclose(1, numpy.linalg.det(S))
      True
  
      """
      normal = unit_vector(normal[:3])
      direction = unit_vector(direction[:3])
      if abs(numpy.dot(normal, direction)) > 1e-6:
          raise ValueError(
              "direction and normal vectors are not orthogonal"
          )
      angle = math.tan(angle)
      M = numpy.identity(4)
      M[:3, :3] += angle * numpy.outer(direction, normal)
      M[:3, 3] = -angle * numpy.dot(point[:3], normal) * direction
      return M
  
imports: import math
  import numpy
label_desc: Use image node_id 13 for calling a global function with example usage: shear_matrix(angle, direction, point, normal) and returns: M