def polarization_vector(phi, theta, alpha, beta, p):
    r"""This function returns a unitary vector describing the polarization
    of plane waves. It recieves as arguments:

    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]])

    """
    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]])

    return R3*R2*R1*epsilon
def polarization_vector ( phi , theta , alpha , beta , p ) : 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 ] ] ) return R3 * R2 * R1 * epsilon
r"""This function returns a unitary vector describing the polarization
    of plane waves. It recieves as arguments:

    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]])