docs: def polarization_vector(phi, theta, alpha, beta, p,
                          numeric=False, abstract=False):
      """This function returns a unitary vector describing the polarization
      of plane waves.:
  
      INPUT:
  
      -  ``phi`` - The spherical coordinates azimuthal angle of the wave vector\
   k.
      -  ``theta`` - The spherical coordinates polar angle of the wave vector k.
      -  ``alpha`` - The rotation of a half-wave plate.
      -  ``beta`` - The rotation of a quarter-wave plate.
      -  ``p`` - either 1 or -1 to indicate whether to return epsilon^(+) or\
   epsilon^(-) respectively.
  
      If alpha and beta are zero, the result will be linearly polarized light
      along some fast axis. alpha and beta are measured from that fast axis.
  
      Propagation towards y, linear polarization (for pi transitions):
  
      >>> from sympy import pi
      >>> polarization_vector(phi=pi/2, theta=pi/2, alpha=pi/2, beta= 0,p=1)
      Matrix([
      [0],
      [0],
      [1]])
  
      Propagation towards +z, circular polarization (for sigma + transitions):
  
      >>> polarization_vector(phi=0, theta= 0, alpha=pi/2, beta= pi/8,p=1)
      Matrix([
      [  -sqrt(2)/2],
      [-sqrt(2)*I/2],
      [           0]])
  
      Propagation towards -z, circular polarization for sigma + transitions:
  
      >>> polarization_vector(phi=0, theta=pi, alpha=   0, beta=-pi/8,p=1)
      Matrix([
      [  -sqrt(2)/2],
      [-sqrt(2)*I/2],
      [           0]])
  
      Components + and - are complex conjugates of each other
  
      >>> from sympy import symbols
      >>> phi, theta, alpha, beta = symbols("phi theta alpha beta", real=True)
      >>> ep = polarization_vector(phi,theta,alpha,beta, 1)
      >>> em = polarization_vector(phi,theta,alpha,beta,-1)
      >>> ep-em.conjugate()
      Matrix([
      [0],
      [0],
      [0]])
  
      We can also define abstract polarization vectors without explicit \
      components
  
      >>> polarization_vector(0, 0, 0, 0, 1, abstract=True)
      epsilonp
      >>> polarization_vector(0, 0, 0, 0, -1, abstract=True)
      epsilonm
  
      """
      if abstract:
          Nl = symbols("N_l", integer=True)
          if p == 1:
              epsilon = Vector3D(IndexedBase("epsilonp", shape=(Nl,)))
          else:
              epsilon = Vector3D(IndexedBase("epsilonm", shape=(Nl,)))
          return epsilon
  
      epsilon = Matrix([cos(2*beta), p*I*sin(2*beta), 0])
  
      R1 = Matrix([[cos(2*alpha), -sin(2*alpha), 0],
                  [sin(2*alpha), cos(2*alpha), 0],
                  [0, 0, 1]])
  
      R2 = Matrix([[cos(theta), 0, sin(theta)],
                   [0, 1, 0],
                   [-sin(theta), 0, cos(theta)]])
  
      R3 = Matrix([[cos(phi), -sin(phi), 0],
                   [sin(phi), cos(phi), 0],
                   [0, 0, 1]])
  
      epsilon = R3*R2*R1*epsilon
      if numeric:
          epsilon = nparray([complex(epsilon[i]) for i in range(3)])
      return epsilon