module BZInt3D
using DoubleFloats
using StaticArrays
using Base.Threads
export Lin3DRuleÎ˜,Lin3DRuleÎ´,Quad3DRuleÎ˜,Quad3DRuleÎ´,
       Quad3DRuleÎ˜ð�”‡,Quad3DRuleÎ˜Î´,Quad3DRuleÎ˜Î˜,Quad3DRuleÎ˜Î˜ð�”‡,
       Quad3DRuleÎ´Î´,Quad3DRuleÎ˜Î˜Î´,Quad3DRuleð�’²,Quad3DRuleð�’²ð�’²,Quad3DRuleð�’²ð�’²ð�’²,Quad3DRuleð�’²ð�”‡,Quad3DRuleð�’²ð�’²ð�”‡
       
       
include("SplitMesh.jl")
include("QuadTetra.jl")
include("TetraSupply.jl")

"""
Linear tetrahedron method for weight function W(k) = Î˜(eF-E(k))
with BlÃ¶chl correction
"""
function Lin3DRuleÎ˜(Emesh,eF)
    Tetras = Mesh2Tetra(size(Emesh)...)
    ETetras = Emesh[Tetras]
    WTetras = zeros(size(ETetras)...)

    @views @threads for i in 1:size(Tetras,1)
        WTetras[i,:] = LinTetraÎ˜_BlÃ¶chl(SVector{4}(ETetras[i,:]),eF)
    end

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

    return Wmesh/size(Tetras,1)*6
end
"""
Linear tetrahedron method for weight function W(k) = Î´(eF-E(k))
"""
function Lin3DRuleÎ´(Emesh,eF)
    Tetras = Mesh2Tetra(size(Emesh)...)
    ETetras = Emesh[Tetras]
    WTetras = zeros(size(ETetras)...)

    @views @threads for i in 1:size(Tetras,1)
        WTetras[i,:] = LinTetraÎ´(SVector{4}(ETetras[i,:]),eF)
    end

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

    return Wmesh/size(Tetras,1)*6
end
"""
Recursive tetrahedron method for weight function W(k) = Î˜(eF-E(k))
"""
function Quad3DRuleÎ˜(Emesh,eF,iter=2)
    QTetras = Mesh2QTetra(size(Emesh)...)
    EQTetras = Emesh[QTetras] 
    WQTetras = zeros(size(EQTetras)...) 
    @views @threads for i in 1:size(QTetras,1)
        WQTetras[i,:] = QuadTetraÎ˜(SVector{10}(EQTetras[i,:]),eF,iter)
    end

    Wmesh = zeros(typeof(float(eF)),size(Emesh)...)
    @views for i in 1:size(QTetras,1)
        Wmesh[QTetras[i,:]] += WQTetras[i,:]
    end
    return Wmesh/size(QTetras,1)*6
end
"""
Recursive tetrahedron method for weight function W(k) = W(k) = 1/D(k) Î˜(eF-E(k))
"""
function Quad3DRuleÎ˜ð�”‡(Emesh,eF,Dmesh,iter=2)
    QTetras = Mesh2QTetra(size(Emesh)...)
    EQTetras = Emesh[QTetras]
    DQTetras = Dmesh[QTetras]
    WQTetras = zeros(size(EQTetras)...)
    @views @threads for i in 1:size(QTetras,1)
        WQTetras[i,:] = QuadTetraÎ˜ð�”‡(SVector{10}(EQTetras[i,:]),eF,SVector{10}(DQTetras[i,:]),iter)
    end
    
    Wmesh = zeros(typeof(float(eF)),size(Emesh)...)
    @views for i in 1:size(QTetras,1)
        Wmesh[QTetras[i,:]] += WQTetras[i,:]
    end
    return Wmesh/size(QTetras,1)*6
end
"""
Recursive tetrahedron method for weight function W(k) = Î´(eF-E(k))
"""
function Quad3DRuleÎ´(Emesh,eF,iter=2)
    QTetras = Mesh2QTetra(size(Emesh)...)
    EQTetras = Emesh[QTetras]
    WQTetras = zeros(size(EQTetras)...)
    @views @threads for i in 1:size(QTetras,1)
        WQTetras[i,:] = QuadTetraÎ´(SVector{10}(EQTetras[i,:]),eF,iter)
    end
    
    Wmesh = zeros(typeof(float(eF)),size(Emesh)...)
    @views for i in 1:size(QTetras,1)
        Wmesh[QTetras[i,:]] += WQTetras[i,:]
    end
    return Wmesh/size(QTetras,1)*6
end

"""
Recursive tetrahedron method for weight function W(k) = Î´(D(k)) Î˜(eF-E(k))
"""
function Quad3DRuleÎ˜Î´(Emesh,eF,Dmesh,iter=2)
    QTetras = Mesh2QTetra(size(Emesh)...)
    EQTetras = Emesh[QTetras]
    DQTetras = Dmesh[QTetras]
    WQTetras = zeros(size(EQTetras)...)
    @views @threads for i in 1:size(QTetras,1)
        WQTetras[i,:] = QuadTetraÎ˜Î´(SVector{10}(EQTetras[i,:]),eF,SVector{10}(DQTetras[i,:]),iter)
    end
    
    Wmesh = zeros(typeof(float(eF)),size(Emesh)...)
    @views for i in 1:size(QTetras,1)
        Wmesh[QTetras[i,:]] += WQTetras[i,:]
    end
    return Wmesh/size(QTetras,1)*6
end

"""
Recursive tetrahedron method for weight function W(k) = Î˜(x1(k))*Î˜(x2(k))
"""
function Quad3DRuleÎ˜Î˜(X1mesh,X2mesh,iter=2)
    QTetras = Mesh2QTetra(size(X1mesh)...)
    X1QTetras = X1mesh[QTetras]
    X2QTetras = X2mesh[QTetras]
    WQTetras = zeros(size(X1QTetras)...)
    @views @threads for i in 1:size(QTetras,1)
        WQTetras[i,:] = QuadTetraÎ˜Î˜(SVector{10}(X1QTetras[i,:]),SVector{10}(X2QTetras[i,:]),iter)
    end

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

"""
Recursive tetrahedron method for weight function W(k) = 1/D(k) Î˜(x1(k))*Î˜(x2(k))
"""
function Quad3DRuleÎ˜Î˜ð�”‡(X1mesh,X2mesh,Dmesh,iter=2)
    QTetras = Mesh2QTetra(size(X1mesh)...)
    X1QTetras = X1mesh[QTetras]
    X2QTetras = X2mesh[QTetras]
    DQTetras = Dmesh[QTetras]
    WQTetras = zeros(size(X1QTetras)...)
    @views @threads for i in 1:size(QTetras,1)
        WQTetras[i,:] = QuadTetraÎ˜Î˜ð�”‡(SVector{10}(X1QTetras[i,:]),SVector{10}(X2QTetras[i,:]),SVector{10}(DQTetras[i,:]),iter)
    end

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

"""
Recursive tetrahedron method for weight function W(k) = Î´(D(k)) Î˜(x1(k))*Î˜(x2(k))
"""
function Quad3DRuleÎ˜Î˜Î´(X1mesh,X2mesh,Dmesh,iter=2)
    QTetras = Mesh2QTetra(size(X1mesh)...)
    X1QTetras = X1mesh[QTetras]
    X2QTetras = X2mesh[QTetras]
    DQTetras = Dmesh[QTetras]
    WQTetras = zeros(size(X1QTetras)...)
    @views @threads for i in 1:size(QTetras,1)
        WQTetras[i,:] = QuadTetraÎ˜Î˜Î´(SVector{10}(X1QTetras[i,:]),SVector{10}(X2QTetras[i,:]),SVector{10}(DQTetras[i,:]),iter)
    end

    Wmesh = zeros(typeof(float(X1mesh[1])),size(X1mesh)...)
    @views for i in 1:size(QTetras,1)
        Wmesh[QTetras[i,:]] += WQTetras[i,:]
    end
    return Wmesh/size(QTetras,1)*6
end
"""
Recursive tetrahedron method for weight function W(k) = Î´(x1(k))*Î´(x2(k))
"""
function Quad3DRuleÎ´Î´(X1mesh,X2mesh,iter=2)
    dx1 = maximum(abs.(X1mesh))*5.0e-4
    xF = zero(X1mesh[1])
    Wmesh0 = Quad3DRuleÎ˜Î´(X1mesh,xF,X2mesh,iter)
    Wmeshd = Quad3DRuleÎ˜Î´(X1mesh,dx1+xF,X2mesh,iter)
    Wmesh = (Wmeshd-Wmesh0)/dx1
    return Wmesh
end


"""
Recursive tetrahedron method for weight function W(k) = ð�’²(x1(k))
default integral type is wtype=Float64
"""
function Quad3DRuleð�’²(ð�’²,X1mesh,iter=2,wtype=Float64)
    QTetras = Mesh2QTetra(size(X1mesh)...)
    X1QTetras = X1mesh[QTetras] 
    WQTetras = zeros(size(X1QTetras)...) 
    @views @threads for i in 1:size(QTetras,1)
        WQTetras[i,:] = QuadTetrað�’²(ð�’²,SVector{10}(X1QTetras[i,:]),iter,wtype)
    end

    Wmesh = zeros(size(X1mesh)...)
    @views for i in 1:size(QTetras,1)
        Wmesh[QTetras[i,:]] += WQTetras[i,:]
    end
    return Wmesh/size(QTetras,1)*6
end

"""
Recursive tetrahedron method for weight function W(k) = ð�’²1(x1(k))*ð�’²2(x2(k))
default integral type is wtype=Float64
"""
function Quad3DRuleð�’²ð�’²(ð�’²1,ð�’²2,X1mesh,X2mesh,iter=2,wtype=Float64)
    QTetras = Mesh2QTetra(size(X1mesh)...)
    X1QTetras = X1mesh[QTetras] 
    X2QTetras = X2mesh[QTetras] 
    WQTetras = zeros(size(X1QTetras)...) 
    @views @threads for i in 1:size(QTetras,1)
        WQTetras[i,:] = QuadTetrað�’²ð�’²(ð�’²1,ð�’²2,SVector{10}(X1QTetras[i,:]),SVector{10}(X2QTetras[i,:]),iter,wtype)
    end

    Wmesh = zeros(size(X1mesh)...)
    @views for i in 1:size(QTetras,1)
        Wmesh[QTetras[i,:]] += WQTetras[i,:]
    end
    return Wmesh/size(QTetras,1)*6
end

"""
Recursive tetrahedron method for weight function W(k) = ð�’²1(x1(k))*ð�’²2(x2(k))*ð�’²3(x3(k))
default integral type is wtype=Float64
"""
function Quad3DRuleð�’²ð�’²ð�’²(ð�’²1,ð�’²2,ð�’²3,X1mesh,X2mesh,X3mesh,iter=2,wtype=Float64)
    QTetras = Mesh2QTetra(size(X1mesh)...)
    X1QTetras = X1mesh[QTetras] 
    X2QTetras = X2mesh[QTetras] 
    X3QTetras = X3mesh[QTetras] 
    WQTetras = zeros(size(X1QTetras)...) 
    @views @threads for i in 1:size(QTetras,1)
        WQTetras[i,:] = QuadTetrað�’²ð�’²ð�’²(ð�’²1,ð�’²2,ð�’²3,SVector{10}(X1QTetras[i,:]),SVector{10}(X2QTetras[i,:]),SVector{10}(X3QTetras[i,:]),iter,wtype)
    end

    Wmesh = zeros(size(X1mesh)...)
    @views for i in 1:size(QTetras,1)
        Wmesh[QTetras[i,:]] += WQTetras[i,:]
    end
    return Wmesh/size(QTetras,1)*6
end

"""
Recursive tetrahedron method for weight function W(k) = ð�’²(x1(k))/D(k)
default integral type is wtype=Float64
"""
function Quad3DRuleð�’²ð�”‡(ð�’²,X1mesh,Dmesh,iter=2,wtype=Float64)
    QTetras = Mesh2QTetra(size(X1mesh)...)
    X1QTetras = X1mesh[QTetras] 
    DQTetras = Dmesh[QTetras] 
    WQTetras = zeros(size(X1QTetras)...) 
    @views @threads for i in 1:size(QTetras,1)
        WQTetras[i,:] = QuadTetrað�’²ð�”‡(ð�’²,SVector{10}(X1QTetras[i,:]),SVector{10}(DQTetras[i,:]),iter,wtype)
    end

    Wmesh = zeros(size(X1mesh)...)
    @views for i in 1:size(QTetras,1)
        Wmesh[QTetras[i,:]] += WQTetras[i,:]
    end
    return Wmesh/size(QTetras,1)*6
end

"""
Recursive tetrahedron method for weight function W(k) = ð�’²1(x1(k))*ð�’²2(x2(k))/D(k)
default integral type is wtype=Float64
"""
function Quad3DRuleð�’²ð�’²ð�”‡(ð�’²1,ð�’²2,X1mesh,X2mesh,Dmesh,iter=2,wtype=Float64)
    QTetras = Mesh2QTetra(size(X1mesh)...)
    X1QTetras = X1mesh[QTetras] 
    X2QTetras = X2mesh[QTetras] 
    DQTetras = Dmesh[QTetras] 
    WQTetras = zeros(size(X1QTetras)...) 
    @views @threads for i in 1:size(QTetras,1)
        WQTetras[i,:] = QuadTetrað�’²ð�’²ð�”‡(ð�’²1,ð�’²2,SVector{10}(X1QTetras[i,:]),SVector{10}(X2QTetras[i,:]),SVector{10}(DQTetras[i,:]),iter,wtype)
    end

    Wmesh = zeros(size(X1mesh)...)
    @views for i in 1:size(QTetras,1)
        Wmesh[QTetras[i,:]] += WQTetras[i,:]
    end
    return Wmesh/size(QTetras,1)*6
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
 
# Quad3DRuleð�‘“(Emesh,Î¼,Î²,iter=2)=Quad3DRuleð�’²(e->fÏµk(e-Î¼,Î²),Emesh,iter)

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

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

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

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


end

