code: def spherical_matrix(theta, phi, axes='sxyz'):
      
      result = euler_matrix(0.0, phi, theta, axes=axes)
      return result
docs: Give a spherical coordinate vector, find the rotation that will
      transform a [0,0,1] vector to those coordinates
  
      Parameters
      -----------
      theta: float, rotation angle in radians
      phi:   float, rotation angle in radians
  
      Returns
      ----------
      matrix: (4,4) rotation matrix where the following will
               be a cartesian vector in the direction of the
               input spherical coordinats:
                  np.dot(matrix, [0,0,1,0])