interp.jl

#Function for numerical integration
function Kernel_Weight_Function(r::Float64)
    if r < 1/2
        return 8/π * (1 - 6r^2 + 6r^3)
    elseif r < 1
        return 8/π * 2 * (1-r)^3
    else
        return zero(r)
    end
end

#Numerical integration
function integr_ker(q_xy::Float64)
    dq = 0.0001
    F = zero(q_xy)
    qz = 0.0
    while qz < sqrt(1 - q_xy^2)
        F += Kernel_Weight_Function(sqrt(qz^2 + q_xy^2)) * dq
        qz += dq
    end
    return 2 * F
end


using Interpolations


unit_q = collect(0:0.0001:1)
F_qxy = zeros(length(unit_q))

#Obtaining Value for interpolation
for i in 1:length(unit_q)
    F_qxy[i] = integr_ker(unit_q[i])
end
knotsq = (unit_q,)
const Numerical_Kernel = interpolate(knotsq, F_qxy, Gridded(Linear()))

#For time test
function f(Times::Int64)
    A = rand(0:2,Times)
    for AA in A
        F = Numerical_Kernel[AA]
    end
end

@time f(1)
@time f(100)
@time f(1000)