docs: def get_matrix(theta, phi, psi, dx, dy, dz,
                     matrix=np.zeros((4, 4), dtype=DTYPE),
                     angles=np.zeros(3, dtype=DTYPE)):
          """
          Build the rotation-translation matrix.
  
          It has the form:
  
          [         | dx ]
          [    R    | dy ]
          [         | dz ]
          [ 0  0  0 | 1  ]
          """
          # NB!: matrix and angles by default are being overwritten on each call
          # thus, only created once at compile time.
  
          angles[0] = theta
          angles[1] = phi
          angles[2] = psi
  
          cx, cy, cz = np.cos(angles)
          sx, sy, sz = np.sin(angles)
  
          rotation = matrix[:3, :3]
          rotation.flat = (cx * cz - sx * cy * sz,
                           cx * sz + sx * cy * cz, sx * sy,
                           -sx * cz - cx * cy * sz,
                           -sx * sz + cx * cy * cz, cx * sy,
                           sy * sz,
                           -sy * cz, cy)
  
          # Translation component
          matrix[:3, 3] = dx, dy, dz
          matrix[3, 3] = 1.
          return matrix