
"""
`module JAC.SphericalTensor`  
    ... a submodel of JAC that help deal with and simplify composed (single- and many-electron) spherical tensor operators,
        matrix elements and spherical amplitudes into a reduced matrix-element representation.
"""
module  SphericalTensor
  
    using   SymEngine,  ..AngularMomentum, ..Basics,  ..Defaults,  ..RacahAlgebra
    
    export  SphericalState, SphericalMatrixElement
    
    
    """
    `abstract type SphericalTensor.AbstractSphericalTensor` 
        ... defines an abstract type to comprise various spherical (electron) operators and spherical fields.
    """
    abstract type  AbstractSphericalTensor                          end
    
    
    """
    `abstract type SphericalTensor.AbstractSphericalOperator` 
        ... defines an abstract type to comprise various spherical (electron) operators.

        + CkOperator                ... to represent a C^(K) tensor operator.
        + CoulombOperator           ... to represent a V^(0, Coulomb) tensor operator.
        + DipoleOperator            ... to represent a D^(1) dipole operator.
        + MultipoleOperator         ... to represent a O^(L, multipole) operator.
        + TkOperator                ... to represent a general T^(K) tensor operator.
    """
    abstract type  AbstractSphericalOperator  <:  AbstractSphericalTensor  end


    """
    `struct  SphericalTensor.CkOperator  <: AbstractSphericalOperator`   ... to represent a C^(K) spherical tensor operator.

        + star          ::Bool      ... defines the complex-conjugate form, if true. 
        + k             ::Basic     ... Tensor rank. 
    """
    struct   CkOperator  <: AbstractSphericalOperator
        star            ::Bool 
        k               ::Basic
    end

    
    """
    `SphericalTensor.CkOperator(k::Basic)`  ... constructor for just setting the rank (star = false).
    """
    function CkOperator(k::Basic)
    	CkOperator(false, k)
    end


    function Base.string(op::CkOperator)
        sa = "C^(" * string(op.k);        if op.star   sa = sa * "*"   end;      sa = sa * ")"
        return( sa )
    end

    function Base.show(io::IO, op::CkOperator)
        sa = string(op) * " spherical tensor operator.";       print(io, sa, "\n")
    end


    """
    `struct  SphericalTensor.CoulombOperator  <: AbstractSphericalOperator`   ... to represent a V^(0, Coulomb) spherical tensor operator.

        + star          ::Bool      ... defines the complex-conjugate form, if true. 
    """
    struct   CoulombOperator  <: AbstractSphericalOperator
        star            ::Bool 
    end

    
    """
    `SphericalTensor.CoulombOperator()`  ... constructor for just setting the dipole operator (star = false).
    """
    function CoulombOperator()
    	CoulombOperator(false)
    end


    function Base.string(op::CoulombOperator)
        sa = "V^(0, Coulomb";        if op.star   sa = sa * "*"   end;      sa = sa * ")" 
        return( sa )
    end

    function Base.show(io::IO, op::CoulombOperator)
        sa = string(op) * " rank-0 Coulomb operator.";       print(io, sa, "\n")
    end


    """
    `struct  SphericalTensor.DipoleOperator  <: AbstractSphericalOperator`   ... to represent a D^(1) spherical tensor operator.

        + star          ::Bool      ... defines the complex-conjugate form, if true. 
    """
    struct   DipoleOperator  <: AbstractSphericalOperator
        star            ::Bool 
    end

    
    """
    `SphericalTensor.DipoleOperator()`  ... constructor for just setting the dipole operator (star = false).
    """
    function DipoleOperator()
    	DipoleOperator(false)
    end


    function Base.string(op::DipoleOperator)
        sa = "D^(1";        if op.star   sa = sa * "*"   end;      sa = sa * ")"
        return( sa )
    end

    function Base.show(io::IO, op::DipoleOperator)
        sa = string(op) * " spherical dipole operator.";       print(io, sa, "\n")
    end


    """
    `struct  SphericalTensor.TkOperator  <: AbstractSphericalOperator`   ... to represent a general rank-K T^(K) spherical tensor operator.

        + star          ::Bool      ... defines the complex-conjugate form, if true. 
        + k             ::Basic     ... Tensor rank. 
    """
    struct   TkOperator  <: AbstractSphericalOperator
        star            ::Bool 
        k               ::Basic
    end

    
    """
    `SphericalTensor.TkOperator(k::Basic)`  ... constructor for just setting the rank (star = false).
    """
    function TkOperator(k::Basic)
    	TkOperator(false, k)
    end


    function Base.string(op::TkOperator)
        sa = "T^(" * string(op.k);        if op.star   sa = sa * "*"   end;      sa = sa * ")"
        return( sa )
    end

    function Base.show(io::IO, op::TkOperator)
        sa = string(op) * " spherical tensor operator.";       print(io, sa, "\n")
    end
    
    
    """
    `abstract type SphericalTensor.AbstractSphericalConstant  <:  AbstractSphericalTensor` 
        ... defines an abstract type to comprise various (parametric-dependent) scalar constants, such as field amplitudes, etc.

        + ScalarConstant            ... to represent a scalar constant to comprise all all factors that do not
                                        affect the tensorial structure of matrix elements and amplitues.
    """
    abstract type  AbstractSphericalConstant   <:  AbstractSphericalTensor  end


    """
    `struct  SphericalTensor.ScalarConstant  <: AbstractSphericalConstant`   
        ... to represent a scalar constant which does not affect the tensorial evaluation of matrix elements.

        + star          ::Bool      ... defines the complex-conjugate form, if true. 
        + symbol        ::Basic     ... symbol of the constant.
    """
    struct   ScalarConstant  <: AbstractSphericalConstant
        star            ::Bool 
        symbol          ::Basic
    end


    function Base.string(op::ScalarConstant)
        sa = string(op.symbol);        if op.star   sa = sa * "^(*)"   end
        return( sa )
    end

    function Base.show(io::IO, op::ScalarConstant)
        sa = string(op) * " scalar constant.";       print(io, sa, "\n")
    end
    
    
    """
    `abstract type SphericalTensor.AbstractSphericalField  <:  AbstractSphericalTensor` 
        ... defines an abstract type to comprise various (external and electron-independent) field operators.

        + UVector                   ... to represent a u^(1) polarization (unit) vector.
        + TwistedU                  ... to represent a u^(1, twisted) unit vector.
    """
    abstract type  AbstractSphericalField   <:  AbstractSphericalTensor  end


    """
    `struct  SphericalTensor.UVector  <: AbstractSphericalField`   ... to represent a u^(1) polarization (unit) vector.

        + star          ::Bool      ... defines the complex-conjugate form, if true. 
    """
    struct   UVector  <: AbstractSphericalField
        star            ::Bool 
    end

    
    """
    `SphericalTensor.UVector()`  ... constructor for just setting the vector (star = false).
    """
    function UVector()
    	UVector(false)
    end


    function Base.string(op::UVector)
        sa = "u^(1";        if op.star   sa = sa * "*"   end;      sa = sa * ")"
        return( sa )
    end

    function Base.show(io::IO, op::UVector)
        sa = string(op) * " spherical polarization (unit) vector.";       print(io, sa, "\n")
    end


    
    """
    `struct  SphericalTensor.TensorProduct  <:  AbstractSphericalTensor` 
        ... defines the coupled tensor product of two spherical tensors [T1^(k1) * T2^(k2)]^(k).

        + star          ::Bool      ... defines the complex-conjugate form, if true. 
        + k             ::Basic     ... (total) tensor rank k.
        + tensor1       ::AbstractSphericalTensor  ... spherical tensor T1.
        + tensor2       ::AbstractSphericalTensor  ... spherical tensor T2.
    """
    struct   TensorProduct  <:  AbstractSphericalTensor
        star            ::Bool 
        k               ::Basic
        tensor1         ::AbstractSphericalTensor
        tensor2         ::AbstractSphericalTensor
    end

    
    """
    `SphericalTensor.TensorProduct(k::Basic, tensor1::AbstractSphericalTensor, tensor2::AbstractSphericalTensor)`  
        ... constructor for just setting the tensor product (star = false).
    """
    function TensorProduct(k::Basic, tensor1::AbstractSphericalTensor, tensor2::AbstractSphericalTensor)
    	TensorProduct(false, k, tensor1, tensor2)
    end


    function Base.string(op::TensorProduct)
        sa = "[" * string(op.tensor1) * " * " * string(op.tensor2) * "]^(" * string(op.k);        if op.star   sa = sa * "*"   end;      sa = sa * ")"
        return( sa )
    end

    function Base.show(io::IO, op::TensorProduct)
        sa = string(op) * " spherical tensor product.";       print(io, sa, "\n")
    end
    
    
    """
    `abstract type SphericalTensor.AbstractSphericalTensorComp` 
        ... defines an abstract type to comprise various spherical (electron) operators and spherical field tensor components.
    """
    abstract type  AbstractSphericalTensorComp     end


    """
    `struct  SphericalTensor.SphericalOperatorComp  <: AbstractSphericalTensorComp`   ... to represent a spherical tensor component.

        + tensor        ::AbstractSphericalOperator   ... tensor
        + q             ::Basic                       ... q-component
    """
    struct   SphericalOperatorComp  <: AbstractSphericalTensorComp
        tensor          ::AbstractSphericalTensor
        q               ::Basic
    end


    function Base.string(comp::SphericalOperatorComp)
        sa = string(comp.tensor) * "_" * string(comp.q)
        return( sa )
    end

    function Base.show(io::IO, comp::SphericalOperatorComp)
        sa = string(comp) * " spherical operator component.";       print(io, sa, "\n")
    end


    """
    `struct  SphericalTensor.SphericalFieldComp  <: AbstractSphericalTensorComp`   ... to represent a spherical field component.

        + tensor        ::AbstractSphericalField    ... exernal tensor field
        + q             ::Basic                     ... q-component
    """
    struct   SphericalFieldComp  <: AbstractSphericalTensorComp
        tensor          ::AbstractSphericalField
        q               ::Basic
    end


    function Base.string(comp::SphericalFieldComp)
        sa = string(comp.tensor) * "_" * string(comp.q)
        return( sa )
    end

    function Base.show(io::IO, comp::SphericalFieldComp)
        sa = string(comp) * " spherical field component.";       print(io, sa, "\n")
    end

    

    """
    `struct  SphericalTensor.SphericalState`   ... to represent a  |jm > state (spherical representation).

        + star          ::Bool      ... defines the complex-conjugate form, if true. 
        + j             ::Basic     ... total angular momentum
        + m             ::Basic     ... projection of total angular momentum
    """
    struct   SphericalState
        star            ::Bool 
        j               ::Basic
        m               ::Basic
    end

    
    """
    `SphericalTensor.SphericalState(j::Basic, m::Basic)`  ... constructor for just setting the ket-vector (star = false).
    """
    function SphericalState(j::Basic, m::Basic)
    	SphericalState(false, j, m)
    end


    function Base.string(state::SphericalState; plain=false)
        if  plain   sa =       string(state.j) * " " * string(state.m);        if state.star   sa = sa * " *"   end
        else        sa = "|" * string(state.j) * " " * string(state.m);        if state.star   sa = sa * " *"   end;      sa = sa * ">"
        end
        return( sa )
    end

    function Base.show(io::IO, state::SphericalState)
        sa = string(state) * " spherical state.";       print(io, sa, "\n")
    end

    

    """
    `struct  SphericalTensor.SphericalMatrixElement`   
        ... to represent a (full)  <ja ma|| T^(k)_q |jb mb>  matrix element of the spherical operator T^(k).

        + star             ::Bool                         ... defines the complex-conjugate form, if true. 
        + leftState        ::SphericalState               ... (reduced) bra-state.
        + constant         ::ScalarConstant               ... a scalar constant that is treated independently.
        + tensor           ::AbstractSphericalTensor      ... spherical tensor.
        + q                ::Basic                        ... q-quantum number of tensor
        + rightState       ::SphericalState               ... (reduced) ket-state.
    """
    struct   SphericalMatrixElement
        star               ::Bool
        leftState          ::SphericalState
        constant           ::ScalarConstant
        tensor             ::AbstractSphericalTensor
        q                  ::Basic 
        rightState         ::SphericalState
    end

    
    """
    `SphericalTensor.SphericalMatrixElement(leftState::SphericalState, tensor::AbstractSphericalTensor, q::Basic, rightState::SphericalState)`  
        ... constructor for just setting the matrix element without constant and with star = false.
    """
    function SphericalMatrixElement(leftState::SphericalState, tensor::AbstractSphericalTensor, q::Basic, rightState::SphericalState)
    	SphericalMatrixElement(false, leftState, ScalarConstant(false, Basic(1)), tensor, q, rightState)
    end


    function Base.string(me::SphericalMatrixElement)
        sa = "<" * string(me.leftState, plain=true) * " | " 
        if  me.constant.symbol != Basic(1)    sa = sa * string(me.constant)    end
        sa = sa * string(me.tensor) * "_" * string(me.q) * " | " * string(me.rightState, plain=true) * ">" 
        if  me.star                    sa = sa * "^(*)"                 end 
        return( sa )
    end

    function Base.show(io::IO, me::SphericalMatrixElement)
        sa = string(me) * " reduced matrix element.";       print(io, sa, "\n")
    end


    
    """
    `struct  SphericalTensor.ProductMatrixElement`   
        ... to represent a (reduced) <ja|| T^(k) |jb> matrix element of the spherical operator T^(k).

        + leftState        ::SphericalState                    ... (spherical) bra-state.
        + components       ::Array{SphericalOperatorComp,1}    ... (product of) spherical tensors.
        + rightState       ::SphericalState                    ... (spherical) ket-state.
    """
    struct   ProductMatrixElement
        leftState          ::SphericalState
        components         ::Array{SphericalOperatorComp,1}
        rightState         ::SphericalState
    end


    function Base.string(me::ProductMatrixElement)
        sa = "<" * string(me.leftState, plain=true) * " | " 
        for  comp in me.components    sa = sa * "  " * string(comp)      end
        sa = sa * " | " * string(me.rightState, plain=true) * "> " 
        return( sa )
    end

    function Base.show(io::IO, me::ProductMatrixElement)
        sa = string(me) * " product matrix element.";       print(io, sa, "\n")
    end

    

    """
    `struct  SphericalTensor.SphericalAmplitude`   
        ... to represent a (non-reduced, many-electron) spherical amplitude that may include one or several spherical matrix elements, spherical 
            field operator components as well as  various symbols from Racah's algebra. In general, these spherical amplitudes need further to be 
            simplified to (fully) reduced amplitudes prior to the calculations.

        + rex              ::RacahAlgebra.RacahExpression       ... contains all factors from Racah's algebra including summations and delta-factors
        + constants        ::Array{AbstractSphericalConstant,1} ... product (list) of scalar constants.
        + fields           ::Array{SphericalFieldComp,1}        ... product (list) of external field tensors.
        + productMes       ::Array{ProductMatrixElement,1}      ... product (list) of product matrix elements.
    """
    struct   SphericalAmplitude
        rex                ::RacahAlgebra.RacahExpression 
        constants          ::Array{AbstractSphericalConstant,1}
        fields             ::Array{SphericalFieldComp,1}
        productMes         ::Array{ProductMatrixElement,1}
    end


    function Base.string(amp::SphericalAmplitude)
        sa =      "\nRacahExpression:    " * string(amp.rex )
        sa = sa * "\nConstants:        " 
        for cst in amp.constants       sa = sa * "  " * string(cst)      end
        sa = sa * "\nField components: " 
        for field in amp.fields        sa = sa * "  " * string(field)    end
        sa = sa * "\nproductMes:       " 
        for rme in amp.productMes      sa = sa * "  " * string(rme)      end
        return( sa )
    end

    function Base.show(io::IO, amp::SphericalAmplitude)
        sa = "Spherical amplitude with: " * string(amp);       print(io, sa, "\n")
    end

    

    """
    `struct  SphericalTensor.ReducedState`   ... to represent a (reduced) |j> state (spherical representation).

        + j             ::Basic     ... total angular momentum
    """
    struct   ReducedState
        j               ::Basic
    end


    function Base.string(state::ReducedState; plain=false)
        if  plain      sa = string(state.j)    else    sa = "|" * string(state.j) * ">"    end
        return( sa )
    end

    function Base.show(io::IO, state::ReducedState)
        sa = string(state) * " reduced state.";       print(io, sa, "\n")
    end

    

    """
    `struct  SphericalTensor.ReducedMatrixElement`   
        ... to represent a (reduced) <ja|| T^(k) |jb> matrix element of the spherical operator T^(k).

        + leftState        ::ReducedState                      ... (reduced) bra-state.
        + tensors          ::Array{AbstractSphericalTensor,1}  ... (product of) spherical tensors.
        + rightState       ::ReducedState                      ... (reduced) ket-state.
    """
    struct   ReducedMatrixElement
        leftState          ::ReducedState
        tensors            ::Array{AbstractSphericalTensor,1}
        rightState         ::ReducedState
    end


    function Base.string(me::ReducedMatrixElement)
        sa = "<" * string(me.leftState, plain=true) * " || " 
        for  tn in me.tensors   sa = sa * " " * string(tn)      end
        sa = sa * " || " * string(me.rightState, plain=true) * ">" 
        return( sa )
    end

    function Base.show(io::IO, me::ReducedMatrixElement)
        sa = string(me) * " reduced matrix element.";       print(io, sa, "\n")
    end

    

    """
    `struct  SphericalTensor.ReducedAmplitude`   
        ... to represent a (reduced, many-electron) amplitude that may include one or several reduced matrix elements, spherical field operators
            and various symbols from Racah's algebra. In general, these reduced amplitudes need to be simplified prior to the calculations.

        + rex              ::RacahAlgebra.RacahExpression       ... contains all factors from Racah's algebra including summations and delta-factors
        + fields           ::Array{AbstractSphericalField,1}    ... product (list) of external field tensors.
        + reducedMes       ::Array{ReducedMatrixElement,1}      ... product (list) of reduced matrix elements.
    """
    struct   ReducedAmplitude
        rex                ::RacahAlgebra.RacahExpression 
        fields             ::Array{AbstractSphericalField,1}
        reducedMes         ::Array{ReducedMatrixElement,1}
    end


    function Base.string(amp::ReducedAmplitude)
        sa =      "\nRacahExpression:   " * string(amp.rex )
        sa = sa * "\nFields:           " 
        for field in amp.fields    sa = sa * " " * string(field)    end
        sa = sa * "\nreducedMes:       " 
        for rme in amp.reducedMes  sa = sa * " " * string(rme)      end
        return( sa )
    end

    function Base.show(io::IO, amp::ReducedAmplitude)
        sa = "Reduced amplitude with: " * string(amp);       print(io, sa, "\n")
    end


    
    
    """
    `SphericalTensor.conjugate(me::SphericalTensor.SphericalMatrixElement)`  
        ... conjugates the spherical matrix element by conjugating all parts of the tensor operator and by interchanging the 
            left- and right-hand states. A conjme::SphericalTensor.SphericalMatrixElement is returned.
    """
    function conjugate(me::SphericalTensor.SphericalMatrixElement)
        starx       = !me.star
        leftStatex  = me.rightState
        constantx   = ScalarConstant( !me.constant.star,   me.constant.symbol)
        tensorx     = SphericalTensor.conjugate( me.tensor )
        qx          = me.q
        rightStatex = me.leftState
        conjme      = SphericalTensor.SphericalMatrixElement(starx, leftStatex, constantx, tensorx, qx, rightStatex)
        
        return( conjme )    
    end
    
    
    """
    `SphericalTensor.conjugate(tn::SphericalTensor.TensorProduct)`  
        ... conjugates the (abstract) spherical matrix element by conjugating all parts of the tensor operator and by interchanging 
            the individual tensor components of each tensor product.
    """
    function conjugate(tn::SphericalTensor.TensorProduct)
        
        starx       = tn.star
        kx          = tn.k
        tensor1x    = SphericalTensor.conjugate( tn.tensor2 )
        tensor2x    = SphericalTensor.conjugate( tn.tensor1 )
        conjtn      = SphericalTensor.TensorProduct(starx, kx, tensor1x, tensor2x)
        
        return( conjtn )    
    end

    conjugate(tn::SphericalTensor.CkOperator)      = SphericalTensor.CkOperator( !tn.star, tn.k)
    conjugate(tn::SphericalTensor.CoulombOperator) = SphericalTensor.CoulombOperator( !tn.star)
    conjugate(tn::SphericalTensor.DipoleOperator)  = SphericalTensor.DipoleOperator( !tn.star)
    conjugate(tn::SphericalTensor.TkOperator)      = SphericalTensor.TkOperator( !tn.star, tn.k)

    conjugate(tn::SphericalTensor.ScalarConstant)  = SphericalTensor.ScalarConstant( !tn.star, tn.symbol)
    conjugate(tn::SphericalTensor.UVector)         = SphericalTensor.UVector( !tn.star)
    
    
    """
    `SphericalTensor.expandSphericalMatrixElements(mes::Array{SphericalTensor.SphericalMatrixElement,1})`  
        ... expand the (product of) spherical matrix elements into a reduced spherical amplitude; it separates the scalar constants from the 
            field operators and the reduced (electron) matrix elements. The Wigner-Eckardt theorem is applied to reduce the spherical
            matrix elements and a Clebsch-Gordan expansion for each product of tensor operators. All additional magnetic quantum numbers,
            that are needed for this expansion, are derived from the corresponding tensor ranks. The Clebsch-Gordan/Wigner symbols are
            'collected' in a Racah expression for further simplification. A sphericalAmp::SphericalTensor.SphericalAmplitude is returned.
    """
    function expandSphericalMatrixElements(mes::Array{SphericalTensor.SphericalMatrixElement,1})
        rex         = RacahAlgebra.RacahExpression()
        constants   = AbstractSphericalConstant[]
        fieldComps  = SphericalFieldComp[]
        productMes  = ProductMatrixElement[]
        
        for  mex in mes
           # Always work with the non-conjugated form of spherical ME
           if  mex.star  me = SphericalTensor.conjugate(mex)    else   me = mex     end
           # 'Add' the constants
           if me.constant.symbol != Basic(1)    push!(constants, me.constant)   end
           # Decompose the spherical operator recusively
           wa, wb, wc, wd = SphericalTensor.expandSphericalTensorComponent(me.tensor, me.q)
           rex = rex * wa;    append!(constants, wb);    append!(fieldComps, wc)
           ##x # Apply Wigner-Eckardt theorem and append reduced matrix element
           ##x k     = SphericalTensor.getRank(me.tensor)
           ##x rex   = rex * SphericalTensor.getWignerEckardtFactor(me.leftState, k, me.q, me.rightState)
           if  me.leftState.star  ||  me.rightState.star   error("SphericalState() assumes non-conjugated spherical states.")     end
           @show wd
           prodMe = SphericalTensor.ProductMatrixElement( me.leftState, wd, me.rightState)
           push!(productMes, prodMe )
        end
        sphericalAmp = SphericalTensor.SphericalAmplitude(rex, constants, fieldComps, productMes)
    
        return( sphericalAmp )    
    end
    
    
    """
    `SphericalTensor.expandSphericalTensorComponent(tn::SphericalTensor.AbstractSphericalTensorComp, q::Basic)`  
        ... expand a coupled product of spherical tensors into its (uncoupled) tensor components by means of a Clebsch-Gordan expansion.
            All necessary magnetic quantum numbers are 'derived' from the corresponding tensor ranks: K -> qK, 2 -> q2,  0 -> 0.
            No expansion is made for elementary spherical tensors. The tensor components of this expansion is returned by (the tuple of)
            three variables: (rex::RacahExpression, fields::Array{SphericalTensor.AbstractSphericalField,1}, 
                              operators::Array{SphericalTensor.AbstractSphericalOperator,1}).
    """
    function expandSphericalTensorComponent(tn::SphericalTensor.AbstractSphericalTensor, q::Basic)
        rex           = RacahAlgebra.RacahExpression()
        constants     = SphericalTensor.AbstractSphericalConstant[]
        fieldComps    = SphericalTensor.SphericalFieldComp[]
        operatorComps = SphericalTensor.SphericalOperatorComp[]
        
        if  typeof(tn) in [SphericalTensor.CkOperator, SphericalTensor.CoulombOperator, SphericalTensor.DipoleOperator, 
                           SphericalTensor.TkOperator]
            push!(operatorComps, SphericalOperatorComp(tn,q) )
        elseif   typeof(tn) in [SphericalTensor.ScalarConstant]
            push!(constants, tn )
        elseif   typeof(tn) in [SphericalTensor.UVector]
            push!(fieldComps, SphericalFieldComp(tn,q) )
        elseif   typeof(tn) in [SphericalTensor.TensorProduct]
            B0  = Basic(0)
            k1  = getRank(tn.tensor1);   if  k1 == B0 q1 = B0   else   q1 = Basic("q" * string(k1) * "_" *string(floor(10000*rand()))[1:2]) end
            k2  = getRank(tn.tensor2);   if  k2 == B0 q2 = B0   else   q2 = Basic("q" * string(k2) * "_" *string(floor(10000*rand()))[1:2]) end
            kk  = getRank(tn);           if  kk == B0 qq = B0   else   qq = Basic("q" * string(kk) * "_" *string(floor(10000*rand()))[1:2]) end
            
            rx1, cx1, fx1, ox1 = SphericalTensor.expandSphericalTensorComponent(tn.tensor1, q1)
            rx2, cx2, fx2, ox2 = SphericalTensor.expandSphericalTensorComponent(tn.tensor2, q2)
            rex = rex * rx1 * rx2 * RacahAlgebra.ClebschGordanExpansion(k1, q1, k2, q2, kk, qq)
            append!(constants,  cx1);      append!(constants,  cx2)
            append!(fieldComps, fx1);      append!(fieldComps, fx2)
            append!(operatorComps, ox1);   append!(operatorComps, ox2)
        else     error("stop a")
        end
        
        return( rex, constants, fieldComps, operatorComps )    
    end
    

    getRank(tn::SphericalTensor.CkOperator)      = tn.k
    getRank(tn::SphericalTensor.CoulombOperator) = Basic(0)
    getRank(tn::SphericalTensor.DipoleOperator)  = Basic(1)
    getRank(tn::SphericalTensor.TkOperator)      = tn.k

    getRank(tn::SphericalTensor.ScalarConstant)  = Basic(0)
    getRank(tn::SphericalTensor.UVector)         = Basic(1)

    getRank(tn::SphericalTensor.TensorProduct)   = tn.k
    
    
    """
    `SphericalTensor.getWignerEckardtFactor(leftState::SphericalState, k::Basic, q::Basic, rightState::SphericalState)`  
        ... Returns the pre-factor from the Wigner-Eckardt theorem. A rex::RacahAlgebra.RacahExpression is returned.
    """
    function getWignerEckardtFactor(leftState::SphericalState, k::Basic, q::Basic, rightState::SphericalState)
        rex = RacahAlgebra.ClebschGordan(rightState.j, rightState.m, k, q, leftState.j, leftState.j)
        return( rex )    
    end
   
end # module
   
