module BZInt2D
using DoubleFloats
using StaticArrays
using Base.Threads
export Lin2DRuleÎ˜,Lin2DRuleÎ´,Quad2DRuleÎ˜,Quad2DRuleÎ´,
       Quad2DRuleÎ˜ð�”‡,Quad2DRuleÎ˜Î´,Quad2DRuleÎ˜Î˜,Quad2DRuleÎ˜Î˜ð�”‡,
       Quad2DRuleÎ´Î´,Quad2DRuleÎ˜Î˜Î´,Quad2DRuleð�’²,Quad2DRuleð�’²ð�’²,Quad2DRuleð�’²ð�’²ð�’²,Quad2DRuleð�’²ð�”‡,Quad2DRuleð�’²ð�’²ð�”‡
       
include("SplitMesh.jl")
include("QuadTrig.jl")
include("TrigSupply.jl")


"""
Linear triangle method for weight function W(k) = Î˜(eF-E(k))
with BlÃ¶chl correction
"""
function Lin2DRuleÎ˜(Emesh,eF)
    Trigs = Mesh2Trig(size(Emesh)...)
    ETrigs = Emesh[Trigs]
    WTrigs = zeros(size(ETrigs)...)

    @views @threads for i in 1:size(Trigs,1)
        WTrigs[i,:] = LinTrigÎ˜_BlÃ¶chl(SVector{3}(ETrigs[i,:]),eF)
    end

    Wmesh = zeros(typeof(float(eF)),size(Emesh)...)
    @views for i in 1:size(Trigs,1)
        Wmesh[Trigs[i,:]] += WTrigs[i,:]
    end

    return Wmesh/size(Trigs,1)*2
end
"""
Linear triangle method for weight function W(k) = Î´(eF-E(k))
"""
function Lin2DRuleÎ´(Emesh,eF)
    Trigs = Mesh2Trig(size(Emesh)...)
    ETrigs = Emesh[Trigs]
    WTrigs = zeros(size(ETrigs)...)

    @views @threads for i in 1:size(Trigs,1)
        WTrigs[i,:] = LinTrigÎ´(SVector{3}(ETrigs[i,:]),eF)
    end

    Wmesh = zeros(typeof(float(eF)),size(Emesh)...)
    @views for i in 1:size(Trigs,1)
        Wmesh[Trigs[i,:]] += WTrigs[i,:]
    end

    return Wmesh/size(Trigs,1)*2
end
"""
Recursive triangle method for weight function W(k) = Î˜(eF-E(k))
"""
function Quad2DRuleÎ˜(Emesh,eF,iter=2)
    QTrigs = Mesh2QTrig(size(Emesh)...)
    EQTrigs = Emesh[QTrigs]
    WQTrigs = zeros(size(EQTrigs)...)
    @views @threads for i in 1:size(QTrigs,1)
        WQTrigs[i,:] = QuadTrigÎ˜(SVector{6}(EQTrigs[i,:]),eF,iter)
    end

    Wmesh = zeros(typeof(float(eF)),size(Emesh)...)
    @views for i in 1:size(QTrigs,1)
        Wmesh[QTrigs[i,:]] += WQTrigs[i,:]
    end
    return Wmesh/size(QTrigs,1)*2
end
"""
Recursive triangle method for weight function W(k) = W(k) = 1/D(k) Î˜(eF-E(k))
"""
function Quad2DRuleÎ˜ð�”‡(Emesh,eF,Dmesh,iter=2)
    QTrigs = Mesh2QTrig(size(Emesh)...)
    EQTrigs = Emesh[QTrigs]
    DQTrigs = Dmesh[QTrigs]
    WQTrigs = zeros(size(EQTrigs)...)
    @views @threads for i in 1:size(QTrigs,1)
        WQTrigs[i,:] = QuadTrigÎ˜ð�”‡(SVector{6}(EQTrigs[i,:]),eF,SVector{6}(DQTrigs[i,:]),iter)
    end
    
    Wmesh = zeros(typeof(float(eF)),size(Emesh)...)
    @views for i in 1:size(QTrigs,1)
        Wmesh[QTrigs[i,:]] += WQTrigs[i,:]
    end
    return Wmesh/size(QTrigs,1)*2
end
"""
Recursive triangle method for weight function W(k) = Î´(eF-E(k))
"""
function Quad2DRuleÎ´(Emesh,eF,iter=2)
    QTrigs = Mesh2QTrig(size(Emesh)...)
    EQTrigs = Emesh[QTrigs]
    WQTrigs = zeros(size(EQTrigs)...)
    @views @threads for i in 1:size(QTrigs,1)
        WQTrigs[i,:] = QuadTrigÎ´(SVector{6}(EQTrigs[i,:]),eF,iter)
    end
    
    Wmesh = zeros(typeof(float(eF)),size(Emesh)...)
    @views for i in 1:size(QTrigs,1)
        Wmesh[QTrigs[i,:]] += WQTrigs[i,:]
    end
    return Wmesh/size(QTrigs,1)*2
end
"""
Recursive triangle method for weight function W(k) = Î´(D(k)) Î˜(eF-E(k))
"""
function Quad2DRuleÎ˜Î´(Emesh,eF,Dmesh,iter=2)
    QTrigs = Mesh2QTrig(size(Emesh)...)
    EQTrigs = Emesh[QTrigs]
    DQTrigs = Dmesh[QTrigs]
    WQTrigs = zeros(size(EQTrigs)...)
    @views @threads for i in 1:size(QTrigs,1)
        WQTrigs[i,:] = QuadTrigÎ˜Î´(SVector{6}(EQTrigs[i,:]),eF,SVector{6}(DQTrigs[i,:]),iter)
    end
    
    Wmesh = zeros(typeof(float(eF)),size(Emesh)...)
    @views for i in 1:size(QTrigs,1)
        Wmesh[QTrigs[i,:]] += WQTrigs[i,:]
    end
    return Wmesh/size(QTrigs,1)*2
end

"""
Recursive triangle method for weight function W(k) = Î˜(x1(k))*Î˜(x2(k))
"""
function Quad2DRuleÎ˜Î˜(X1mesh,X2mesh,iter=2)
    QTrigs = Mesh2QTrig(size(X1mesh)...)
    X1QTrigs = X1mesh[QTrigs]
    X2QTrigs = X2mesh[QTrigs]
    WQTrigs = zeros(size(X1QTrigs)...)
    @views @threads for i in 1:size(QTrigs,1)
        WQTrigs[i,:] = QuadTrigÎ˜Î˜(SVector{6}(X1QTrigs[i,:]),SVector{6}(X2QTrigs[i,:]),iter)
    end

    Wmesh = zeros(typeof(float(X1mesh[1])),size(X1mesh)...)
    @views for i in 1:size(QTrigs,1)
        Wmesh[QTrigs[i,:]] += WQTrigs[i,:]
    end
    return Wmesh/size(QTrigs,1)*2
end

"""
Recursive triangle method for weight function W(k) = 1/D(k) Î˜(x1(k))*Î˜(x2(k))
"""
function Quad2DRuleÎ˜Î˜ð�”‡(X1mesh,X2mesh,Dmesh,iter=2)
    QTrigs = Mesh2QTrig(size(X1mesh)...)
    X1QTrigs = X1mesh[QTrigs]
    X2QTrigs = X2mesh[QTrigs]
    DQTrigs = Dmesh[QTrigs]
    WQTrigs = zeros(size(X1QTrigs)...)
    @views @threads for i in 1:size(QTrigs,1)
        WQTrigs[i,:] = QuadTrigÎ˜Î˜ð�”‡(SVector{6}(X1QTrigs[i,:]),SVector{6}(X2QTrigs[i,:]),SVector{6}(DQTrigs[i,:]),iter)
    end

    Wmesh = zeros(typeof(float(X1mesh[1])),size(X1mesh)...)
    @views for i in 1:size(QTrigs,1)
        Wmesh[QTrigs[i,:]] += WQTrigs[i,:]
    end
    return Wmesh/size(QTrigs,1)*2
end
"""
Recursive triangle method for weight function W(k) = Î´(D(k)) Î˜(x1(k))*Î˜(x2(k))
"""
function Quad2DRuleÎ˜Î˜Î´(X1mesh,X2mesh,Dmesh,iter=2)
    QTrigs = Mesh2QTrig(size(X1mesh)...)
    X1QTrigs = X1mesh[QTrigs]
    X2QTrigs = X2mesh[QTrigs]
    DQTrigs = Dmesh[QTrigs]
    WQTrigs = zeros(size(X1QTrigs)...)
    @views @threads for i in 1:size(QTrigs,1)
        WQTrigs[i,:] = QuadTrigÎ˜Î˜Î´(SVector{6}(X1QTrigs[i,:]),SVector{6}(X2QTrigs[i,:]),SVector{6}(DQTrigs[i,:]),iter)
    end

    Wmesh = zeros(typeof(float(X1mesh[1])),size(X1mesh)...)
    @views for i in 1:size(QTrigs,1)
        Wmesh[QTrigs[i,:]] += WQTrigs[i,:]
    end
    return Wmesh/size(QTrigs,1)*2
end

"""
Recursive triangle method for weight function W(k) = Î´(x1(k))*Î´(x2(k))
"""
function Quad2DRuleÎ´Î´(X1mesh,X2mesh,iter=2)
    dx1 = maximum(abs.(X1mesh))*5.0e-4
    xF = zero(X1mesh[1])
    Wmesh0 = Quad2DRuleÎ˜Î´(X1mesh,xF,X2mesh,iter)
    Wmeshd = Quad2DRuleÎ˜Î´(X1mesh,dx1+xF,X2mesh,iter)
    Wmesh = (Wmeshd-Wmesh0)/dx1
    return Wmesh
end


"""
Recursive triangle method for weight function W(k) = ð�’²(x1(k))
default integral type is wtype=Float64
"""
function Quad2DRuleð�’²(ð�’²,X1mesh,iter=2,wtype=Float64)
    QTrigs = Mesh2QTrig(size(X1mesh)...)
    X1QTrigs = X1mesh[QTrigs] 
    WQTrigs = zeros(size(X1QTrigs)...) 
    @views @threads for i in 1:size(QTrigs,1)
        WQTrigs[i,:] = QuadTrigð�’²(ð�’²,SVector{6}(X1QTrigs[i,:]),iter,wtype)
    end

    Wmesh = zeros(size(X1mesh)...)
    @views for i in 1:size(QTrigs,1)
        Wmesh[QTrigs[i,:]] += WQTrigs[i,:]
    end
    return Wmesh/size(QTrigs,1)*2
end

"""
Recursive triangle method for weight function W(k) = ð�’²1(x1(k))*ð�’²2(x2(k))
default integral type is wtype=Float64
"""
function Quad2DRuleð�’²ð�’²(ð�’²1,ð�’²2,X1mesh,X2mesh,iter=2,wtype=Float64)
    QTrigs = Mesh2QTrig(size(X1mesh)...)
    X1QTrigs = X1mesh[QTrigs] 
    X2QTrigs = X2mesh[QTrigs] 
    WQTrigs = zeros(size(X1QTrigs)...) 
    @views @threads for i in 1:size(QTrigs,1)
        WQTrigs[i,:] = QuadTrigð�’²ð�’²(ð�’²1,ð�’²2,SVector{6}(X1QTrigs[i,:]),SVector{6}(X2QTrigs[i,:]),iter,wtype)
    end

    Wmesh = zeros(size(X1mesh)...)
    @views for i in 1:size(QTrigs,1)
        Wmesh[QTrigs[i,:]] += WQTrigs[i,:]
    end
    return Wmesh/size(QTrigs,1)*2
end

"""
Recursive triangle method for weight function W(k) = ð�’²1(x1(k))*ð�’²2(x2(k))*ð�’²3(x3(k))
default integral type is wtype=Float64
"""
function Quad2DRuleð�’²ð�’²ð�’²(ð�’²1,ð�’²2,ð�’²3,X1mesh,X2mesh,X3mesh,iter=2,wtype=Float64)
    QTrigs = Mesh2QTrig(size(X1mesh)...)
    X1QTrigs = X1mesh[QTrigs] 
    X2QTrigs = X2mesh[QTrigs] 
    X3QTrigs = X3mesh[QTrigs] 
    WQTrigs = zeros(size(X1QTrigs)...) 
    @views @threads for i in 1:size(QTrigs,1)
        WQTrigs[i,:] = QuadTrigð�’²ð�’²(ð�’²1,ð�’²2,ð�’²3,SVector{6}(X1QTrigs[i,:]),SVector{6}(X2QTrigs[i,:]),SVector{6}(X3QTrigs[i,:]),iter,wtype)
    end

    Wmesh = zeros(size(X1mesh)...)
    @views for i in 1:size(QTrigs,1)
        Wmesh[QTrigs[i,:]] += WQTrigs[i,:]
    end
    return Wmesh/size(QTrigs,1)*2
end

"""
Recursive triangle method for weight function W(k) = ð�’²(x1(k))/ D(k)
default integral type is wtype=Float64
"""
function Quad2DRuleð�’²ð�”‡(ð�’²,X1mesh,Dmesh,iter=2,wtype=Float64)
    QTrigs = Mesh2QTrig(size(X1mesh)...)
    X1QTrigs = X1mesh[QTrigs] 
    DQTrigs = Dmesh[QTrigs] 
    WQTrigs = zeros(size(X1QTrigs)...) 
    @views @threads for i in 1:size(QTrigs,1)
        WQTrigs[i,:] = QuadTrigð�’²ð�”‡(ð�’²,SVector{6}(X1QTrigs[i,:]),SVector{6}(DQTrigs[i,:]),iter,wtype)
    end

    Wmesh = zeros(size(X1mesh)...)
    @views for i in 1:size(QTrigs,1)
        Wmesh[QTrigs[i,:]] += WQTrigs[i,:]
    end
    return Wmesh/size(QTrigs,1)*2
end

"""
Recursive triangle method for weight function W(k) = ð�’²1(x1(k))*ð�’²2(x2(k))/ D(k)
default integral type is wtype=Float64
"""
function Quad2DRuleð�’²ð�’²ð�”‡(ð�’²1,ð�’²2,X1mesh,X2mesh,Dmesh,iter=2,wtype=Float64)
    QTrigs = Mesh2QTrig(size(X1mesh)...)
    X1QTrigs = X1mesh[QTrigs] 
    X2QTrigs = X2mesh[QTrigs] 
    DQTrigs = Dmesh[QTrigs] 
    WQTrigs = zeros(size(X1QTrigs)...) 
    @views @threads for i in 1:size(QTrigs,1)
        WQTrigs[i,:] = QuadTrigð�’²ð�’²ð�”‡(ð�’²1,ð�’²2,SVector{6}(X1QTrigs[i,:]),SVector{6}(X2QTrigs[i,:]),SVector{6}(DQTrigs[i,:]),iter,wtype)
    end

    Wmesh = zeros(size(X1mesh)...)
    @views for i in 1:size(QTrigs,1)
        Wmesh[QTrigs[i,:]] += WQTrigs[i,:]
    end
    return Wmesh/size(QTrigs,1)*2
end

function fÏµk(Ïµ_Î¼,Î²)
    x = Ïµ_Î¼*Î²
    if x>20 
        res = 0.0
    elseif x<-20
        res = 1.0
    else
        res = 1.0/(1+exp(x))
    end
    return res
end

function dfÏµk_dÏµ(Ïµ_Î¼,Î²)
    x = Ïµ_Î¼*Î²
    if x>20 || x<-20 
        res = 0.0
    else
        e =exp(x) 
        res = -Î²*e/(1+e)^2
    end
    return res
end

# Quad2DRuleð�‘“(Emesh,Î¼,Î²,iter=2)=Quad2DRuleð�’²(e->fÏµk(e-Î¼,Î²),Emesh,iter)

# Quad2DRuleð�‘‘ð�‘“(Emesh,Î¼,Î²,iter=2)=Quad2DRuleð�’²(e->dfÏµk_dÏµ(e-Î¼,Î²),Emesh,iter)

# Quad2DRuleð�‘‘ð�‘“ð�‘‘ð�‘“(E1mesh,E2mesh,Î¼,Î²,iter=2)=Quad2DRuleð�’²ð�’²(e->dfÏµk_dÏµ(e-Î¼,Î²),e->dfÏµk_dÏµ(e-Î¼,Î²),E1mesh,E2mesh,iter)

# Quad2DRuleð�‘“ð�”“(Emesh,Dmesh,Î¼,Î²,Î·,iter=2)=Quad2DRuleð�’²ð�’²(e->dfÏµk_dÏµ(e-Î¼,Î²),d->1/(d+1im*Î·),Emesh,Dmesh,iter,ComplexF64)

# Quad2DRuleð�‘“ð�‘“ð�”“(E1mesh,E2mesh,Dmesh,Î¼,Î²,Î·,iter=2)=Quad2DRuleð�’²ð�’²ð�’²(e->dfÏµk_dÏµ(e-Î¼,Î²),e->dfÏµk_dÏµ(e-Î¼,Î²),d->1/(d+1im*Î·),E1mesh,E2mesh,Dmesh,iter,ComplexF64)

end


