include("LinTrig.jl")
include("QTrigInterp.jl")

"""
Recursive triangle rule in a single quadratic triangle, weight function is W(k) = Î˜(eF-E(k)).
"""
function QuadTrigÎ˜(Eqtrig,eF,iter)
    FloatType = typeof(float(eF))
    if iter==0 
        if maximum(Eqtrig) <= eF
            # qw = FloatType.([0,0,0,1/6,1/6,1/6])
            qw = FloatType.(@SArray[0.041666666666666664, 0.041666666666666664, 
            0.041666666666666664, 0.125, 0.125, 0.125])
        else
            tmp = @MArray zeros(FloatType,6)
            if minimum(Eqtrig) < eF
                @views for i=1:4
                    res = LinTrigÎ˜( (SVector{3}(Eqtrig[qtrig2trig[i,:]])) ,eF)/4
                    for j=1:3
                        @inbounds tmp[qtrig2trig[i,j]] += res[j]
                    end
                end
            end
            qw = SArray(tmp)
        end
    else 
        eqtrigs = QTrigInterpolation(Eqtrig)
        qweights = @MArray zeros(FloatType,4,6)
        for i =1:4
            qweights[i,:]=QuadTrigÎ˜(eqtrigs[i,:],eF,iter-1)
        end
        qw = CollectQWeights(SArray(qweights))
    end
    return qw
end
"""
Recursive triangle rule in a single quadratic triangle, weight function is W(k) = Î˜(X(k)).
"""
QuadTrigÎ˜(Xqtrig,iter) = QuadTrigÎ˜(-Xqtrig,zero(Xqtrig[1]),iter)

"""
Recursive triangle rule in a single quadratic triangle, weight function is W(k) = Î´(eF-E(k)).
"""
function QuadTrigÎ´(Eqtrig,eF,iter)
    FloatType = typeof(float(eF))
        if iter ==0
            if minimum(Eqtrig) < eF < maximum(Eqtrig)
                tmp = @MArray zeros(FloatType,6)
                @views for i=1:4
                    res = LinTrigÎ´( (SVector{3}(Eqtrig[qtrig2trig[i,:]])) ,eF)/4
                    for j=1:3
                        @inbounds tmp[qtrig2trig[i,j]] += res[j]
                    end
                end
                qw = SArray(tmp)
            else
                qw = @SArray zeros(FloatType,6)
            end
        else
            eqtrigs = QTrigInterpolation(Eqtrig)
            qweights = @MArray zeros(FloatType,4,6)
            for i =1:4
                qweights[i,:]=QuadTrigÎ´(eqtrigs[i,:],eF,iter-1)
            end
            qw = CollectQWeights(SArray(qweights))
        end
    
    return qw
end
"""
Recursive triangle rule in a single quadratic triangle, weight function is W(k) = Î´(X(k)).
"""
QuadTrigÎ´(Xqtrig,iter)=QuadTrigÎ´(-Xqtrig,zero(Xqtrig[1]),iter)

"""
Recursive triangle rule in a single quadratic triangle, weight function is W(k) = Î˜(eF-E(k))Î´(D(k)).
"""
function QuadTrigÎ˜Î´(Eqtrig,eF,Dqtrig,iter)
    FloatType = typeof(float(eF))
    if iter ==0
        if minimum(Eqtrig) < eF 
            tmp = @MArray zeros(FloatType,6)
            @views for i=1:4
                res = LinTrigÎ˜Î´(SVector{3}(Eqtrig[qtrig2trig[i,:]]) ,eF,SVector{3}(Dqtrig[qtrig2trig[i,:]]))/4
                for j=1:3
                    @inbounds tmp[qtrig2trig[i,j]] +=res[j]
                end
            end
            qw = SArray(tmp)
        else
            qw = @SArray zeros(FloatType,6)
        end
    else
        eqtrigs = QTrigInterpolation(Eqtrig)
        dqtrigs = QTrigInterpolation(Dqtrig)
        qweights = @MArray zeros(FloatType,4,6)
        for i =1:4
        # for i=1:4    
            qweights[i,:]=QuadTrigÎ˜Î´(eqtrigs[i,:],eF,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)).
"""
QuadTrigÎ˜Î´(X1qtrig,X2qtrig,iter) = QuadTrigÎ˜Î´(-X1qtrig,zero(X1qtrig[1]),X2qtrig,iter) 

"""
Recursive triangle rule in a single quadratic triangle, weight function is W(k) = Î˜(eF-E(k))â‹…1/D(k)
"""
function QuadTrigÎ˜ð�”‡(Eqtrig,eF,Dqtrig,iter)
    FloatType = typeof(float(eF))
    dmax = maximum(Dqtrig)
    dmin = minimum(Dqtrig)
    if dmax*dmin>0 && isapprox(dmax,dmin,rtol = 0.25)
        qw = QuadTrigÎ˜(Eqtrig,eF,iter) ./ Dqtrig
    else 
        if iter ==0
            tmp = @MArray zeros(FloatType,6)
            if minimum(Eqtrig) < eF 
                @views for i=1:4
                    res = 
                    LinTrigÎ˜ð�”‡(SVector{3}(Eqtrig[qtrig2trig[i,:]]),eF,SVector{3}(Dqtrig[qtrig2trig[i,:]]))/4
                    for j=1:3
                        @inbounds tmp[qtrig2trig[i,j]] +=res[j]
                    end
                end
            end
            qw = SArray(tmp)
        else 
            eqtrigs = QTrigInterpolation(Eqtrig)
            dqtrigs = QTrigInterpolation(Dqtrig)
            qweights = @MArray zeros(FloatType,4,6)
            for i =1:4
                qweights[i,:]=QuadTrigÎ˜ð�”‡(eqtrigs[i,:],eF,dqtrigs[i,:],iter-1)
            end
            qw = CollectQWeights(SArray(qweights))
        end
    end
    return qw
end
"""
Recursive triangle rule in a single quadratic triangle, weight function is W(k) = Î˜(X(k))â‹…1/D(k)
"""
QuadTrigÎ˜ð�”‡(Xqtrig,Dqtrig,iter) = QuadTrigÎ˜ð�”‡(-Xqtrig,zero(Xqtrig[1]),Dqtrig,iter)