"""
Recursive tetrahedron rule in a single quadratic tetrahedron, weight function is W(k) = Î˜(X1(k))*Î˜(X2(k))
Î˜(x1)*Î˜(x2) = (Î˜(x1)+Î˜(x2)-Î˜(-x1*x2))/2
"""
function QuadTetraÎ˜Î˜(X1qtetra,X2qtetra,iter=2)
    FloatType = typeof(float(X1qtetra[1]))
    if iter==0
        if minimum(X1qtetra)>0 && minimum(X2qtetra)>0
            qw = QuadTetraÎ˜((@SArray zeros(FloatType,10)),one(FloatType),0)
        elseif maximum(X1qtetra)<0 && maximum(X2qtetra)<0
            qw = @SArray zeros(FloatType,10)
        else 
        Eqtetra = X1qtetra.*X2qtetra
        qw = (-QuadTetraÎ˜(-Eqtetra,0)+QuadTetraÎ˜(X1qtetra,0)+
              QuadTetraÎ˜(X2qtetra,0))/2
        end
    else 
        x1qtetras = QTetraInterpolation(X1qtetra)
        x2qtetras = QTetraInterpolation(X2qtetra)
        qweights = @MArray zeros(FloatType,8,10)
        for i =1:8
            qweights[i,:]=QuadTetraÎ˜Î˜(x1qtetras[i,:],x2qtetras[i,:],iter-1)
        end
        qw = CollectQWeights(SArray(qweights))
    end
    return qw
end

"""
Recursive tetrahedron rule in a single quadratic tetrahedron, weight function is W(k) = Î˜(X1(k))â‹…Î˜(X2(k))â‹…1/D(k)
Î˜(x1)*Î˜(x2) = (Î˜(x1)+Î˜(x2)-Î˜(-x1*x2))/2
"""
function QuadTetraÎ˜Î˜ğ�”‡(X1qtetra,X2qtetra,Dqtetra,iter=2)
    FloatType = typeof(float(X1qtetra[1]))
    if iter==0
        if minimum(X1qtetra)>0 && minimum(X2qtetra)>0
            qw = QuadTetraÎ˜ğ�”‡((@SArray zeros(FloatType,10)),one(FloatType),Dqtetra,0)
        elseif maximum(X1qtetra)<0 && maximum(X2qtetra)<0
            qw = @SArray zeros(FloatType,10)
        else
        Eqtetra = X1qtetra.*X2qtetra
        qw = (-QuadTetraÎ˜ğ�”‡(-Eqtetra,Dqtetra,0)+QuadTetraÎ˜ğ�”‡(X1qtetra,Dqtetra,0)+
              QuadTetraÎ˜ğ�”‡(X2qtetra,Dqtetra,0))/2
        end
    else 
        x1qtetras = QTetraInterpolation(X1qtetra)
        x2qtetras = QTetraInterpolation(X2qtetra)
        dqtetras = QTetraInterpolation(Dqtetra)
        qweights = @MArray zeros(FloatType,8,10)
        for i =1:8
            qweights[i,:]=QuadTetraÎ˜Î˜ğ�”‡(x1qtetras[i,:],x2qtetras[i,:],dqtetras[i,:],iter-1)
        end
        qw = CollectQWeights(SArray(qweights))
    end
    return qw
end


"""
Recursive tetrahedron rule in a single quadratic tetrahedron, weight function is W(k) = Î˜(X1(k))â‹…Î˜(X2(k))â‹…Î´(D(k))
Î˜(x1)*Î˜(x2) = (Î˜(x1)+Î˜(x2)-Î˜(-x1*x2))/2
"""
function QuadTetraÎ˜Î˜Î´(X1qtetra,X2qtetra,Dqtetra,iter=2)
    FloatType = typeof(float(X1qtetra[1]))
    if iter==0
        if minimum(X1qtetra)>0 && minimum(X2qtetra)>0
            qw = QuadTetraÎ˜Î´((@SArray zeros(FloatType,10)),one(FloatType),Dqtetra,0)
        elseif maximum(X1qtetra)<0 && maximum(X2qtetra)<0
            qw = @SArray zeros(FloatType,10)
        else
        Eqtetra = X1qtetra.*X2qtetra
        qw = (-QuadTetraÎ˜Î´(-Eqtetra,Dqtetra,0)+QuadTetraÎ˜Î´(X1qtetra,Dqtetra,0)+
        QuadTetraÎ˜Î´(X2qtetra,Dqtetra,0))/2
        end
    else 
        x1qtetras = QTetraInterpolation(X1qtetra)
        x2qtetras = QTetraInterpolation(X2qtetra)
        dqtetras = QTetraInterpolation(Dqtetra)
        qweights = @MArray zeros(FloatType,8,10)
        for i =1:8
            qweights[i,:]=QuadTetraÎ˜Î˜Î´(x1qtetras[i,:],x2qtetras[i,:],dqtetras[i,:],iter-1)
        end
        qw = CollectQWeights(SArray(qweights))
    end
    return qw
end

"""
Recursive tetrahedron rule in a single quadratic tetrahedron, weight function is W(k) = ğ�’²(X1(k))
I think ğ�’² could be a complex function that takes in complex input
default integral type is wtype=Float64 
"""
function QuadTetrağ�’²(ğ�’²,X1qtetra,iter=2,wtype=Float64)
    FloatType = typeof(abs(X1qtetra[1]))
    WType = wtype==Float64 ? eltype(X1qtetra) : wtype
    if iter==0
        Wqtetra = @MArray zeros(WType,10)
        for i=1:10
            Wqtetra[i] = ğ�’²(X1qtetra[i])
        end
        if maximum(abs.(Wqtetra))==0
            qw = @SArray zeros(WType,10)
        else
            qw = QuadTetraÎ˜((@SArray zeros(FloatType,10)),one(FloatType),0).*SArray(Wqtetra)
        end
    else 
        x1qtetras = QTetraInterpolation(X1qtetra)
        qweights = @MArray zeros(WType,8,10)
        for i =1:8
            qweights[i,:]=QuadTetrağ�’²(ğ�’²,x1qtetras[i,:],iter-1,WType)
        end
        qw = CollectQWeights(SArray(qweights))
    end
    return qw
end

"""
Recursive tetrahedron rule in a single quadratic tetrahedron,weight function is W(k) = ğ�’²1(X1(k))*ğ�’²2(X2(k))
I think ğ�’²1,ğ�’²2 could be a complex function that takes in complex input
default integral type is wtype=Float64 
"""
function QuadTetrağ�’²ğ�’²(ğ�’²1,ğ�’²2,X1qtetra,X2qtetra,iter=2,wtype=Float64)
    FloatType = typeof(abs(X1qtetra[1]))
    WType = wtype==Float64 ? eltype(X1qtetra) : wtype
    if iter==0
        W1qtetra = @MArray zeros(WType,10)
        W2qtetra = @MArray zeros(WType,10)
        for i=1:10
            W1qtetra[i] = ğ�’²1(X1qtetra[i])
            W2qtetra[i] = ğ�’²2(X2qtetra[i])
        end
        if maximum(abs.(W1qtetra))==0 || maximum(abs.(W2qtetra))==0
            qw = @SArray zeros(WType,10)
        else
            qw = QuadTetraÎ˜((@SArray zeros(FloatType,10)),one(FloatType),0).*SArray(W1qtetra.*W2qtetra)
        end
    else 
        x1qtetras = QTetraInterpolation(X1qtetra)
        x2qtetras = QTetraInterpolation(X2qtetra)
        qweights = @MArray zeros(WType,8,10)
        for i =1:8
            qweights[i,:]=QuadTetrağ�’²ğ�’²(ğ�’²1,ğ�’²2,x1qtetras[i,:],x2qtetras[i,:],iter-1,WType)
        end
        qw = CollectQWeights(SArray(qweights))
    end
    return qw
end

"""
Recursive tetrahedron rule in a single quadratic tetrahedron,weight function is W(k) = ğ�’²1(X1(k))*ğ�’²2(X2(k))*ğ�’²3(X3(k))
I think ğ�’²1,ğ�’²2,ğ�’²3 could be a complex function that takes in complex input
default integral type is wtype=Float64 
"""
function QuadTetrağ�’²ğ�’²ğ�’²(ğ�’²1,ğ�’²2,ğ�’²3,X1qtetra,X2qtetra,X3qtetra,iter=2,wtype=Float64)
    FloatType = typeof(abs(X1qtetra[1]))
    WType = wtype==Float64 ? eltype(X1qtetra) : wtype
    if iter==0
        W1qtetra = @MArray zeros(WType,10)
        W2qtetra = @MArray zeros(WType,10)
        W3qtetra = @MArray zeros(WType,10)
        for i=1:10
            W1qtetra[i] = ğ�’²1(X1qtetra[i])
            W2qtetra[i] = ğ�’²2(X2qtetra[i])
            W3qtetra[i] = ğ�’²3(X3qtetra[i])
        end
        if maximum(abs.(W1qtetra))==0 || maximum(abs.(W2qtetra))==0|| maximum(abs.(W3qtetra))==0
            qw = @SArray zeros(WType,10)
        else
            qw = QuadTetraÎ˜((@SArray zeros(FloatType,10)),one(FloatType),0).*SArray(W1qtetra.*W2qtetra.*W3qtetra)
        end
    else 
        x1qtetras = QTetraInterpolation(X1qtetra)
        x2qtetras = QTetraInterpolation(X2qtetra)
        x3qtetras = QTetraInterpolation(X3qtetra)
        qweights = @MArray zeros(WType,8,10)
        for i =1:8
            qweights[i,:]=QuadTetrağ�’²ğ�’²ğ�’²(ğ�’²1,ğ�’²2,ğ�’²3,x1qtetras[i,:],x2qtetras[i,:],x3qtetras[i,:],iter-1,WType)
        end
        qw = CollectQWeights(SArray(qweights))
    end
    return qw
end
        
"""
Recursive tetrahedron rule in a single quadratic tetrahedron,weight function is  W(k) = ğ�’²(X1(k))/D(k)
I think ğ�’² could be a complex function that takes in complex input
default integral type is wtype=Float64 
"""
function QuadTetrağ�’²ğ�”‡(ğ�’²,X1qtetra,Dqtetra,iter=2,wtype=Float64)
    FloatType = typeof(Dqtetra[1])
    WType = wtype==Float64 ? eltype(X1qtetra) : wtype
    if iter==0
        Wqtetra = @MArray zeros(WType,10)
        for i=1:10
            Wqtetra[i] = ğ�’²(X1qtetra[i])
        end
        if maximum(abs.(Wqtetra))==0
            qw = @SArray zeros(WType,10)
        else
            qw = QuadTetraÎ˜ğ�”‡((@SArray zeros(FloatType,10)),one(FloatType),Dqtetra,0).*SArray(Wqtetra)
        end
    else 
        x1qtetras = QTetraInterpolation(X1qtetra)
        dqtetras = QTetraInterpolation(Dqtetra)
        qweights = @MArray zeros(WType,8,10)
        for i =1:8
            qweights[i,:]=QuadTetrağ�’²ğ�”‡(ğ�’²,x1qtetras[i,:],dqtetras[i,:],iter-1,WType)
        end
        qw = CollectQWeights(SArray(qweights))
    end
    return qw
end

"""
Recursive tetrahedron rule in a single quadratic tetrahedron,weight function is W(k) = ğ�’²1(X1(k))ğ�’²2(X2(k))/D(k)
I think ğ�’²1,ğ�’²2 could be a complex function that takes in complex input
default integral type is wtype=Float64 
"""
function QuadTetrağ�’²ğ�’²ğ�”‡(ğ�’²1,ğ�’²2,X1qtetra,X2qtetra,Dqtetra,iter=2,wtype=Float64)
    FloatType = typeof(Dqtetra[1])
    WType = wtype==Float64 ? eltype(X1qtetra) : wtype
    if iter==0
        W1qtetra = @MArray zeros(WType,10)
        W2qtetra = @MArray zeros(WType,10)
        for i=1:10
            W1qtetra[i] = ğ�’²1(X1qtetra[i])
            W2qtetra[i] = ğ�’²2(X2qtetra[i])
        end
        if maximum(abs.(W1qtetra))==0 || maximum(abs.(W2qtetra))==0
            qw = @SArray zeros(WType,10)
        else
            qw = QuadTetraÎ˜ğ�”‡((@SArray zeros(FloatType,10)),one(FloatType),Dqtetra,0).*SArray(W1qtetra.*W2qtetra)
        end
    else 
        x1qtetras = QTetraInterpolation(X1qtetra)
        x2qtetras = QTetraInterpolation(X2qtetra)
        dqtetras = QTetraInterpolation(Dqtetra)
        qweights = @MArray zeros(WType,8,10)
        for i =1:8
            qweights[i,:]=QuadTetrağ�’²ğ�’²ğ�”‡(ğ�’²1,ğ�’²2,x1qtetras[i,:],x2qtetras[i,:],dqtetras[i,:],iter-1,WType)
        end
        qw = CollectQWeights(SArray(qweights))
    end
    return qw
end




