"""
Recursive triangle rule in a single quadratic triangle, weight function is W(k) = Î˜(X1(k))*Î˜(X2(k))
Î˜(x1)*Î˜(x2) = (Î˜(x1)+Î˜(x2)-Î˜(-x1*x2))/2
"""
function QuadTrigÎ˜Î˜(X1qtrig,X2qtrig,iter=2)
    FloatType = typeof(float(X1qtrig[1]))
    if iter==0
        if minimum(X1qtrig)>0 && minimum(X2qtrig)>0
            qw = QuadTrigÎ˜((@SArray zeros(FloatType,6)),one(FloatType),0)
        elseif maximum(X1qtrig)<0 && maximum(X2qtrig)<0
            qw = @SArray zeros(FloatType,6)
        else 
        Eqtrig = X1qtrig.*X2qtrig
        qw = (-QuadTrigÎ˜(-Eqtrig,0)+QuadTrigÎ˜(X1qtrig,0)+
              QuadTrigÎ˜(X2qtrig,0))/2
        end
    else 
        x1qtrigs = QTrigInterpolation(X1qtrig)
        x2qtrigs = QTrigInterpolation(X2qtrig)
        qweights = @MArray zeros(FloatType,4,6)
        for i =1:4
            qweights[i,:]=QuadTrigÎ˜Î˜(x1qtrigs[i,:],x2qtrigs[i,:],iter-1)
        end
        qw = CollectQWeights(SArray(qweights))
    end
    return qw
end

"""
Recursive triangle rule in a single quadratic triangle, weight function is W(k) = Î˜(X1(k))â‹…Î˜(X2(k))â‹…1/D(k)
Î˜(x1)*Î˜(x2) = (Î˜(x1)+Î˜(x2)-Î˜(-x1*x2))/2
"""
function QuadTrigÎ˜Î˜ğ�”‡(X1qtrig,X2qtrig,Dqtrig,iter=2)
    FloatType = typeof(float(X1qtrig[1]))
    if iter==0
        if minimum(X1qtrig)>0 && minimum(X2qtrig)>0
            qw = QuadTrigÎ˜ğ�”‡((@SArray zeros(FloatType,6)),one(FloatType),Dqtrig,0)
        elseif maximum(X1qtrig)<0 && maximum(X2qtrig)<0
            qw = @SArray zeros(FloatType,6)
        else
        Eqtrig = X1qtrig.*X2qtrig
        qw = (-QuadTrigÎ˜ğ�”‡(-Eqtrig,Dqtrig,0)+QuadTrigÎ˜ğ�”‡(X1qtrig,Dqtrig,0)+
              QuadTrigÎ˜ğ�”‡(X2qtrig,Dqtrig,0))/2
        end
    else 
        x1qtrigs = QTrigInterpolation(X1qtrig)
        x2qtrigs = QTrigInterpolation(X2qtrig)
        dqtrigs = QTrigInterpolation(Dqtrig)
        qweights = @MArray zeros(FloatType,4,6)
        for i =1:4
            qweights[i,:]=QuadTrigÎ˜Î˜ğ�”‡(x1qtrigs[i,:],x2qtrigs[i,:],dqtrigs[i,:],iter-1)
        end
        qw = CollectQWeights(SArray(qweights))
    end
    return qw
end

"""
Recursive triangle rule in a single quadratic triangle, weight function is W(k) = Î˜(X1(k))â‹…Î˜(X2(k))â‹…Î´(D(k))
Î˜(x1)*Î˜(x2) = (Î˜(x1)+Î˜(x2)-Î˜(-x1*x2))/2
"""
function QuadTrigÎ˜Î˜Î´(X1qtrig,X2qtrig,Dqtrig,iter=2)
    FloatType = typeof(float(X1qtrig[1]))
    if iter==0  
        if minimum(X1qtrig)>0 && minimum(X2qtrig)>0
            qw = QuadTrigÎ˜Î´((@SArray zeros(FloatType,6)),one(FloatType),Dqtrig,0)
        elseif maximum(X1qtrig)<0 && maximum(X2qtrig)<0
            qw = @SArray zeros(FloatType,6)
        else 
            Eqtrig = X1qtrig.*X2qtrig
            qw = (-QuadTrigÎ˜Î´(-Eqtrig,Dqtrig,0)+QuadTrigÎ˜Î´(X1qtrig,Dqtrig,0)+
                QuadTrigÎ˜Î´(X2qtrig,Dqtrig,0))/2
        end
    else 
        x1qtrigs = QTrigInterpolation(X1qtrig)
        x2qtrigs = QTrigInterpolation(X2qtrig)
        dqtrigs = QTrigInterpolation(Dqtrig)
        qweights = @MArray zeros(FloatType,4,6)
        for i =1:4
            qweights[i,:]=QuadTrigÎ˜Î˜Î´(x1qtrigs[i,:],x2qtrigs[i,:],dqtrigs[i,:],iter-1)
        end
        qw = CollectQWeights(SArray(qweights))
    end
    return qw
end

"""
Recursive triangle rule in a single quadratic triangle, 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 QuadTrigğ�’²(ğ�’²,X1qtrig,iter=2,wtype=Float64)
    FloatType = typeof(abs(X1qtrig[1]))
    WType = wtype==Float64 ? eltype(X1qtrig) : wtype
    if iter==0
        Wqtrig = @MArray zeros(WType,6)
        for i=1:6
            Wqtrig[i]= ğ�’²(X1qtrig[i])
        end
        if maximum(abs.(Wqtrig))==0
            qw = @SArray zeros(WType,6)
        else
            qw = QuadTrigÎ˜((@SArray zeros(FloatType,6)),one(FloatType),0).* SArray(Wqtrig)
        end
    else 
        x1qtrigs = QTrigInterpolation(X1qtrig)
        qweights = @MArray zeros(WType,4,6)
        for i =1:4
            qweights[i,:]=QuadTrigğ�’²(ğ�’²,x1qtrigs[i,:],iter-1,WType)
        end
        qw = CollectQWeights(SArray(qweights))
    end
    return qw
end

"""
Recursive triangle rule in a single quadratic triangle,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 QuadTrigğ�’²ğ�’²(ğ�’²1,ğ�’²2,X1qtrig,X2qtrig,iter=2,wtype=Float64)
    FloatType = typeof(abs(X1qtrig[1]))
    WType = wtype==Float64 ? eltype(X1qtrig) : wtype
    if iter==0
        W1qtrig = @MArray zeros(WType,6)
        W2qtrig = @MArray zeros(WType,6)
        for i=1:6
            W1qtrig[i]= ğ�’²1(X1qtrig[i])
            W2qtrig[i]= ğ�’²2(X2qtrig[i])
        end
        if maximum(abs.(W1qtrig))==0 || maximum(abs.(W2qtrig))==0
            qw = @SArray zeros(WType,6)
        else
            qw = QuadTrigÎ˜((@SArray zeros(FloatType,6)),one(FloatType),0).*SArray(W1qtrig.*W2qtrig)
        end
    else 
        x1qtrigs = QTrigInterpolation(X1qtrig)
        x2qtrigs = QTrigInterpolation(X2qtrig)
        qweights = @MArray zeros(WType,4,6)
         for i =1:4
            qweights[i,:]=QuadTrigğ�’²ğ�’²(ğ�’²1,ğ�’²2,x1qtrigs[i,:],x2qtrigs[i,:],iter-1,WType)
        end
        qw = CollectQWeights(SArray(qweights))
    end
    return qw
end

"""
Recursive triangle rule in a single quadratic triangle,weight function is W(k) = ğ�’²1(X1(k))*ğ�’²2(X2(k))*ğ�’²3(X3(k))
I think ğ�’²1,ğ�’²2 could be a complex function that takes in complex input
default integral type is wtype=Float64 
"""
function QuadTrigğ�’²ğ�’²ğ�’²(ğ�’²1,ğ�’²2,ğ�’²3,X1qtrig,X2qtrig,X3qtrig,iter=2,wtype=Float64)
    FloatType = typeof(abs(X1qtrig[1]))
    WType = wtype==Float64 ? eltype(X1qtrig) : wtype
    if iter==0
        W1qtrig = @MArray zeros(WType,6)
        W2qtrig = @MArray zeros(WType,6)
        W3qtrig = @MArray zeros(WType,6)
        for i=1:6
            W1qtrig[i]= ğ�’²1(X1qtrig[i])
            W2qtrig[i]= ğ�’²2(X2qtrig[i])
            W3qtrig[i]= ğ�’²2(X3qtrig[i])
        end
        if maximum(abs.(W1qtrig))==0 || maximum(abs.(W2qtrig))==0|| maximum(abs.(W3qtrig))==0
            qw = @SArray zeros(WType,6)
        else
            qw = QuadTrigÎ˜((@SArray zeros(FloatType,6)),one(FloatType),0).*SArray(W1qtrig.*W2qtrig.*W3qtrig)
        end
    else 
        x1qtrigs = QTrigInterpolation(X1qtrig)
        x2qtrigs = QTrigInterpolation(X2qtrig)
        x3qtrigs = QTrigInterpolation(X3qtrig)
        qweights = @MArray zeros(WType,4,6)
         for i =1:4
            qweights[i,:]=QuadTrigğ�’²ğ�’²ğ�’²(ğ�’²1,ğ�’²2,ğ�’²3,x1qtrigs[i,:],x2qtrigs[i,:],x3qtrigs[i,:],iter-1,WType)
        end
        qw = CollectQWeights(SArray(qweights))
    end
    return qw
end

"""
Recursive triangle rule in a single quadratic triangle,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 QuadTrigğ�’²ğ�”‡(ğ�’²,X1qtrig,Dqtrig,iter=2,wtype=Float64)
    FloatType = typeof(Dqtrig[1])
    WType = wtype==Float64 ? eltype(X1qtrig) : wtype
    if iter==0
        Wqtrig = @MArray zeros(WType,6)
        for i=1:6
            Wqtrig[i]= ğ�’²(X1qtrig[i])
        end
        if maximum(abs.(Wqtrig))==0
            qw = @SArray zeros(WType,6)
        else
            qw = QuadTrigÎ˜ğ�”‡((@SArray zeros(FloatType,6)),one(FloatType),Dqtrig,0).*SArray(Wqtrig)
        end
    else 
        x1qtrigs = QTrigInterpolation(X1qtrig)
        dqtrigs = QTrigInterpolation(Dqtrig)
        qweights =  @MArray zeros(WType,4,6)
        for i =1:4
            qweights[i,:]=QuadTrigğ�’²ğ�”‡(ğ�’²,x1qtrigs[i,:],dqtrigs[i,:],iter-1,WType)
        end
        qw = CollectQWeights(SArray(qweights))
    end
    return qw
end

"""
Recursive triangle rule in a single quadratic triangle,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 QuadTrigğ�’²ğ�’²ğ�”‡(ğ�’²1,ğ�’²2,X1qtrig,X2qtrig,Dqtrig,iter=2,wtype=Float64)
    FloatType = typeof(Dqtrig[1])
    WType = wtype==Float64 ? eltype(X1qtrig) : wtype
    if iter==0
        W1qtrig = @MArray zeros(WType,6)
        W2qtrig = @MArray zeros(WType,6)
        for i=1:6
            W1qtrig[i]= ğ�’²1(X1qtrig[i])
            W2qtrig[i]= ğ�’²2(X2qtrig[i])
        end
        if maximum(abs.(W1qtrig))==0 || maximum(abs.(W2qtrig))==0
            qw = @SArray zeros(WType,6)
        else
            qw = QuadTrigÎ˜ğ�”‡((@SArray zeros(FloatType,6)),one(FloatType),Dqtrig,0).*SArray(W1qtrig.*W2qtrig)
        end
    else 
        x1qtrigs = QTrigInterpolation(X1qtrig)
        x2qtrigs = QTrigInterpolation(X2qtrig)
        dqtrigs = QTrigInterpolation(Dqtrig)
        qweights = @MArray zeros(WType,4,6)
        for i =1:4
            qweights[i,:]=QuadTrigğ�’²ğ�’²ğ�”‡(ğ�’²1,ğ�’²2,x1qtrigs[i,:],x2qtrigs[i,:],dqtrigs[i,:],iter-1,WType)
        end
        qw = CollectQWeights(SArray(qweights))
    end
    return qw
end



