text:
"""
rq_periodic_kernel function created by kernel_coder(). Requires 3 hyperparameters. Likely created using rq_periodic_kernel_base() as an input.
Use with include("src/kernels/rq_periodic_kernel.jl").
hyperparameters == ["α", "λ", "rq_P"]
"""
function rq_periodic_kernel(
    hyperparameters::Vector{<:Real},
    δ::Real;
    dorder::Vector{<:Integer}=zeros(Int64, length(hyperparameters) + 2))

    @assert length(hyperparameters)==3 "hyperparameters is the wrong length"
    @assert length(dorder)==(3 + 2) "dorder is the wrong length"
    even_time_derivative = powers_of_negative_one(dorder[2])
    @assert maximum(dorder) < 3 "No more than two time derivatives for either t1 or t2 can be calculated"

    dorder = append!([sum(dorder[1:2])], dorder[3:end])

    α = hyperparameters[1]
    λ = hyperparameters[2]
    rq_P = hyperparameters[3]

    if dorder==[4, 0, 0, 2]
        func = -320*pi^4*λ*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^6 + 320*pi^4*λ*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^6 + 64*pi^6*δ^2*λ*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^8 - 64*pi^6*δ^2*λ*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^8 - 960*pi^4*λ^2*(-1 - α)*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^6) - 960*pi^4*λ^2*(-1 - α)*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^6) - 640*pi^5*δ*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^7 + 7040*pi^4*λ^2*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^6) + 960*pi^6*δ^2*λ^2*(-1 - α)*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^8) + 960*pi^6*δ^2*λ^2*(-1 - α)*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^8) - 7680*pi^4*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^6) + 7680*pi^4*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^6) + 9600*pi^5*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^7) - 9600*pi^5*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^7) + 960*pi^6*δ^2*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^6*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^8) - 960*pi^6*δ^2*λ^3*(-2 - α)*(-1 - α)*cos(pi*δ/rq_P)^6*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^8) - 6016*pi^6*δ^2*λ^2*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^8) - 5120*pi^4*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^6) - 9600*pi^5*δ*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^7) + 44800*pi^5*δ*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^7) - 9600*pi^5*δ*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^7) + 22080*pi^6*δ^2*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^8) - 22080*pi^6*δ^2*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^8) - 25600*pi^5*δ*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^7) + 25600*pi^5*δ*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^7) - 11520*pi^6*δ^2*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^6*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^8) + 43520*pi^6*δ^2*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^8) - 11520*pi^6*δ^2*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^6*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^8) - 10240*pi^5*δ*λ^5*(-4 - α)*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-5 - α)/(α^4*rq_P^7) - 15360*pi^6*δ^2*λ^5*(-4 - α)*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^6*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-5 - α)/(α^4*rq_P^8) + 15360*pi^6*δ^2*λ^5*(-4 - α)*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^6*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-5 - α)/(α^4*rq_P^8) - 4096*pi^6*δ^2*λ^6*(-4 - α)*(-3 - α)*(-2 - α)*(-1 - α)*(-5 - α)*sin(pi*δ/rq_P)^6*cos(pi*δ/rq_P)^6*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-6 - α)/(α^5*rq_P^8)
    end

    if dorder==[3, 0, 0, 2]
        func = 192*pi^3*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^5 - 128*pi^4*δ*λ*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^6 + 128*pi^4*δ*λ*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^6 - 64*pi^5*δ^2*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^7 - 576*pi^3*λ^2*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) + 576*pi^3*λ^2*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) - 384*pi^4*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^6) - 384*pi^4*δ*λ^2*(-1 - α)*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^6) - 768*pi^3*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) + 2816*pi^4*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^6) + 960*pi^5*δ^2*λ^2*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^7) - 960*pi^5*δ^2*λ^2*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^7) - 3072*pi^4*δ*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^6) + 3072*pi^4*δ*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^6) - 960*pi^5*δ^2*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^7) + 4480*pi^5*δ^2*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^7) - 960*pi^5*δ^2*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^7) - 2048*pi^4*δ*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^6) - 2560*pi^5*δ^2*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^7) + 2560*pi^5*δ^2*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^7) - 1024*pi^5*δ^2*λ^5*(-4 - α)*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-5 - α)/(α^4*rq_P^7)
    end

    if dorder==[2, 0, 0, 2]
        func = 24*pi^2*λ*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 - 24*pi^2*λ*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 - 16*pi^4*δ^2*λ*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^6 + 16*pi^4*δ^2*λ*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^6 + 96*pi^3*δ*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^5 - 96*pi^2*λ^2*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 48*pi^4*δ^2*λ^2*(-1 - α)*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^6) - 48*pi^4*δ^2*λ^2*(-1 - α)*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^6) - 288*pi^3*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) + 288*pi^3*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) + 352*pi^4*δ^2*λ^2*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^6) - 384*pi^3*δ*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) - 384*pi^4*δ^2*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^6) + 384*pi^4*δ^2*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^6) - 256*pi^4*δ^2*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^6)
    end

    if dorder==[1, 0, 0, 2]
        func = -8*pi*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^3 + 16*pi^2*δ*λ*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 - 16*pi^2*δ*λ*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 + 16*pi^3*δ^2*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^5 - 64*pi^2*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 48*pi^3*δ^2*λ^2*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) + 48*pi^3*δ^2*λ^2*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) - 64*pi^3*δ^2*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5)
    end

    if dorder==[0, 0, 0, 2]
        func = 4*pi^2*δ^2*λ*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 - 4*pi^2*δ^2*λ*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 - 8*pi*δ*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^3 - 16*pi^2*δ^2*λ^2*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4)
    end

    if dorder==[4, 0, 1, 1]
        func = 64*pi^4*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^5 - 64*pi^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^5 + 64*pi^5*δ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^6 + 512*pi^4*λ*(-1 - α)*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) + 384*pi^4*λ*(-1 - α)*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) + 384*pi^4*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^6*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) - 2944*pi^4*λ*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) + 4992*pi^4*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) - 7424*pi^4*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) - 1920*pi^5*δ*λ*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^6) + 2048*pi^5*δ*λ*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^6) + 7168*pi^4*λ^3*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^5) - 3072*pi^4*λ^3*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^6*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^5) + 2880*pi^5*δ*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^6) - 15360*pi^5*δ*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^6) + 4800*pi^5*δ*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^6) + 2048*pi^4*λ^4*(-4 - α)*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^6*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-5 - α)/(α^4*rq_P^5) + 12160*pi^5*δ*λ^3*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^6) - 19200*pi^5*δ*λ^3*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^6) + 1920*pi^5*δ*λ^3*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^7*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^6) + 10240*pi^5*δ*λ^4*(-4 - α)*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-5 - α)/(α^4*rq_P^6) - 5120*pi^5*δ*λ^4*(-4 - α)*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^7*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-5 - α)/(α^4*rq_P^6) + 2048*pi^5*δ*λ^5*(-4 - α)*(-3 - α)*(-2 - α)*(-1 - α)*(-5 - α)*sin(pi*δ/rq_P)^7*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-6 - α)/(α^5*rq_P^6)
    end

    if dorder==[3, 0, 1, 1]
        func = -48*pi^3*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 + 16*pi^4*δ*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^5 - 16*pi^4*δ*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^5 + 288*pi^3*λ*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 384*pi^3*λ*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) + 128*pi^4*δ*λ*(-1 - α)*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) + 96*pi^4*δ*λ*(-1 - α)*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) + 864*pi^3*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 288*pi^3*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) + 96*pi^4*δ*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^6*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) - 736*pi^4*δ*λ*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) + 384*pi^3*λ^3*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4) + 1248*pi^4*δ*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) - 1856*pi^4*δ*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) + 1792*pi^4*δ*λ^3*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^5) - 768*pi^4*δ*λ^3*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^6*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^5) + 512*pi^4*δ*λ^4*(-4 - α)*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^6*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-5 - α)/(α^4*rq_P^5)
    end

    if dorder==[2, 0, 1, 1]
        func = -8*pi^2*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^3 + 8*pi^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^3 - 16*pi^3*δ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 - 16*pi^2*λ*(-1 - α)*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) + 80*pi^2*λ*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) + 64*pi^2*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^3) + 96*pi^3*δ*λ*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 128*pi^3*δ*λ*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) + 288*pi^3*δ*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 96*pi^3*δ*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) + 128*pi^3*δ*λ^3*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4)
    end

    if dorder==[1, 0, 1, 1]
        func = 4*pi*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^2 - 4*pi^2*δ*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^3 + 4*pi^2*δ*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^3 + 8*pi*λ*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^2) - 8*pi^2*δ*λ*(-1 - α)*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) + 40*pi^2*δ*λ*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) + 32*pi^2*δ*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^3)
    end

    if dorder==[0, 0, 1, 1]
        func = 4*pi*δ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^2 + 8*pi*δ*λ*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^2)
    end

    if dorder==[4, 1, 0, 1]
        func = -192*pi^4*λ^2*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) - 192*pi^4*λ^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) + 64*pi^4*λ*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^5 - 64*pi^4*λ*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^5 - 16*pi^4*λ*(4*pi*δ*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 8*pi*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)/(α^3*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*pi*δ*λ*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α^2*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 + 16*pi^4*λ*(4*pi*δ*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 8*pi*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)/(α^3*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*pi*δ*λ*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α^2*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 - 192*pi^4*λ^2*(-1 - α)*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^5) - 192*pi^4*λ^2*(-1 - α)*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^5) + 1408*pi^4*λ^2*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) + 1408*pi^4*λ^2*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^5) - 1536*pi^4*λ^3*(-2 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) + 1536*pi^4*λ^3*(-2 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) - 1536*pi^4*λ^3*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) + 1536*pi^4*λ^3*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) + 192*pi^4*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) + 192*pi^4*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) - 48*pi^4*λ^2*(-1 - α)*(4*pi*δ*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 8*pi*δ*λ^2*(-2 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)/(α^3*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*pi*δ*λ*(-2 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α^2*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 48*pi^4*λ^2*(-1 - α)*(4*pi*δ*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 8*pi*δ*λ^2*(-2 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)/(α^3*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*pi*δ*λ*(-2 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α^2*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) + 896*pi^5*δ*λ^2*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^6) - 896*pi^5*δ*λ^2*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^6) + 64*pi^5*δ*λ*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^6 - 3072*pi^4*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^5) + 3072*pi^4*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^5) - 1024*pi^4*λ^4*(-3 - α)*(-2 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^5) - 1024*pi^4*λ^4*(-3 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^5) - 1024*pi^4*λ^4*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^5) - 1408*pi^4*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) + 352*pi^4*λ^2*(-1 - α)*(4*pi*δ*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 8*pi*δ*λ^2*(-2 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)/(α^3*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*pi*δ*λ*(-2 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α^2*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) + 896*pi^5*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^6) - 896*pi^5*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^6) - 960*pi^5*δ*λ^3*(-2 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^6) + 4480*pi^5*δ*λ^3*(-2 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^6) - 960*pi^5*δ*λ^3*(-2 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^6) - 768*pi^5*δ*λ^3*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^6) + 3072*pi^5*δ*λ^3*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^6) - 768*pi^5*δ*λ^3*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^6) - 3072*pi^4*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^4*rq_P^5) + 1536*pi^4*λ^3*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) - 1536*pi^4*λ^3*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) - 384*pi^4*λ^3*(-2 - α)*(-1 - α)*(4*pi*δ*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 8*pi*δ*λ^2*(-3 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)/(α^3*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*pi*δ*λ*(-3 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α^2*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) + 384*pi^4*λ^3*(-2 - α)*(-1 - α)*(4*pi*δ*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 8*pi*δ*λ^2*(-3 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)/(α^3*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*pi*δ*λ*(-3 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α^2*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 1728*pi^5*δ*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^6) + 7552*pi^5*δ*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^6) - 1728*pi^5*δ*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^6) - 2560*pi^5*δ*λ^4*(-3 - α)*(-2 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^6) + 2560*pi^5*δ*λ^4*(-3 - α)*(-2 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^6) - 2560*pi^5*δ*λ^4*(-3 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^6) + 2560*pi^5*δ*λ^4*(-3 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^6) - 1024*pi^5*δ*λ^4*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^6) + 1024*pi^5*δ*λ^4*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^6) - 896*pi^5*δ*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^6) + 896*pi^5*δ*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^6) - 64*pi^5*δ*λ^2*(-1 - α)*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^6) + 64*pi^5*δ*λ^2*(-1 - α)*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^6) + 1024*pi^4*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*(-2*λ*(-4 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^5) - 256*pi^4*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)*(4*pi*δ*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 8*pi*δ*λ^2*(-4 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)/(α^3*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*pi*δ*λ*(-4 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α^2*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))/(α^3*rq_P^4) - 6144*pi^5*δ*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^4*rq_P^6) + 6144*pi^5*δ*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^4*rq_P^6) - 1024*pi^5*δ*λ^5*(-4 - α)*(-3 - α)*(-2 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-5 - α)/(α^4*rq_P^6) - 1024*pi^5*δ*λ^5*(-4 - α)*(-3 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-5 - α)/(α^4*rq_P^6) - 1024*pi^5*δ*λ^5*(-4 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-5 - α)/(α^4*rq_P^6) + 768*pi^5*δ*λ^3*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^6) - 3072*pi^5*δ*λ^3*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^6) + 768*pi^5*δ*λ^3*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^6) + 192*pi^5*δ*λ^3*(-2 - α)*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^6) - 1408*pi^5*δ*λ^3*(-2 - α)*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^6) + 192*pi^5*δ*λ^3*(-2 - α)*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^6) - 3072*pi^5*δ*λ^5*(-4 - α)*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-5 - α)/(α^5*rq_P^6) + 1024*pi^5*δ*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*(-2*λ*(-4 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^6) - 1024*pi^5*δ*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*(-2*λ*(-4 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^6) + 1536*pi^5*δ*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^6) - 1536*pi^5*δ*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^6) + 1024*pi^5*δ*λ^5*(-4 - α)*(-3 - α)*(-2 - α)*(-1 - α)*(-2*λ*(-4 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-5 - α)/(α^4*rq_P^6)
    end

    if dorder==[3, 1, 0, 1]
        func = -144*pi^3*λ^2*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) + 144*pi^3*λ^2*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 48*pi^3*λ*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 + 16*pi^3*λ*(4*pi*δ*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 8*pi*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)/(α^3*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*pi*δ*λ*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α^2*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^3 - 48*pi^4*δ*λ^2*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) - 48*pi^4*δ*λ^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) + 16*pi^4*δ*λ*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^5 - 16*pi^4*δ*λ*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^5 - 144*pi^3*λ^2*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^4) + 144*pi^3*λ^2*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^4) - 192*pi^3*λ^3*(-2 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 192*pi^3*λ^3*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 48*pi^4*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^5) - 48*pi^4*δ*λ^2*(-1 - α)*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^5) + 288*pi^4*δ*λ^2*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) - 384*pi^3*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^4) + 144*pi^3*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 144*pi^3*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 48*pi^3*λ^2*(-1 - α)*(4*pi*δ*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 8*pi*δ*λ^2*(-2 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)/(α^3*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*pi*δ*λ*(-2 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α^2*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) + 48*pi^3*λ^2*(-1 - α)*(4*pi*δ*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 8*pi*δ*λ^2*(-2 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)/(α^3*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*pi*δ*λ*(-2 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α^2*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) + 288*pi^4*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^5) - 384*pi^4*δ*λ^3*(-2 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) + 384*pi^4*δ*λ^3*(-2 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) - 192*pi^4*δ*λ^3*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) + 192*pi^4*δ*λ^3*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) + 48*pi^4*δ*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) + 48*pi^4*δ*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) + 192*pi^3*λ^3*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 64*pi^3*λ^3*(-2 - α)*(-1 - α)*(4*pi*δ*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 8*pi*δ*λ^2*(-3 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)/(α^3*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*pi*δ*λ*(-3 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α^2*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^3) - 576*pi^4*δ*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^5) + 576*pi^4*δ*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^5) - 256*pi^4*δ*λ^4*(-3 - α)*(-2 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^5) - 256*pi^4*δ*λ^4*(-3 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^5) - 288*pi^4*δ*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) - 64*pi^4*δ*λ^2*(-1 - α)*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) - 512*pi^4*δ*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^4*rq_P^5) + 192*pi^4*δ*λ^3*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) - 192*pi^4*δ*λ^3*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) + 192*pi^4*δ*λ^3*(-2 - α)*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) - 192*pi^4*δ*λ^3*(-2 - α)*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) + 256*pi^4*δ*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^5)
    end

    if dorder==[2, 1, 0, 1]
        func = -8*pi^2*λ*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^3 + 8*pi^2*λ*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^3 + 4*pi^2*λ*(4*pi*δ*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 8*pi*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)/(α^3*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*pi*δ*λ*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α^2*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^2 - 4*pi^2*λ*(4*pi*δ*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 8*pi*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)/(α^3*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*pi*δ*λ*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α^2*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^2 - 32*pi^2*λ^2*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) - 32*pi^2*λ^2*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^3) - 32*pi^3*δ*λ^2*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) + 32*pi^3*δ*λ^2*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 16*pi^3*δ*λ*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 + 32*pi^2*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) - 16*pi^2*λ^2*(-1 - α)*(4*pi*δ*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 8*pi*δ*λ^2*(-2 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)/(α^3*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*pi*δ*λ*(-2 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α^2*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^2) - 32*pi^3*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^4) + 32*pi^3*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^4) - 64*pi^3*δ*λ^3*(-2 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 64*pi^3*δ*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^4) + 32*pi^3*δ*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 32*pi^3*δ*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) + 16*pi^3*δ*λ^2*(-1 - α)*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 16*pi^3*δ*λ^2*(-1 - α)*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) + 64*pi^3*δ*λ^3*(-2 - α)*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4)
    end

    if dorder==[1, 1, 0, 1]
        func = 4*pi*λ*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^2 - 4*pi*λ*(4*pi*δ*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 8*pi*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)/(α^3*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*pi*δ*λ*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)/(α^2*rq_P^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P - 4*pi^2*δ*λ*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^3 + 4*pi^2*δ*λ*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^3 + 16*pi^2*δ*λ^2*(-1 - α)*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3)
    end

    if dorder==[0, 1, 0, 1]
        func = 8*pi*δ*λ^2*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^2) + 4*pi*δ*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(2*λ*sin(pi*δ/rq_P)^2/(α*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^2
    end

    if dorder==[4, 0, 0, 1]
        func = 64*pi^4*λ*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^5 - 64*pi^4*λ*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^5 + 192*pi^4*λ^2*(-1 - α)*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) + 192*pi^4*λ^2*(-1 - α)*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) + 64*pi^5*δ*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^6 - 1408*pi^4*λ^2*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) + 1536*pi^4*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) - 1536*pi^4*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) - 960*pi^5*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^6) + 960*pi^5*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^6) + 1024*pi^4*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^5) + 960*pi^5*δ*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^6) - 4480*pi^5*δ*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^6) + 960*pi^5*δ*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^6) + 2560*pi^5*δ*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^6) - 2560*pi^5*δ*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^6) + 1024*pi^5*δ*λ^5*(-4 - α)*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^5*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-5 - α)/(α^4*rq_P^6)
    end

    if dorder==[3, 0, 0, 1]
        func = -48*pi^3*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 + 16*pi^4*δ*λ*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^5 - 16*pi^4*δ*λ*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^5 + 144*pi^3*λ^2*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 144*pi^3*λ^2*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) + 48*pi^4*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) + 48*pi^4*δ*λ^2*(-1 - α)*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) + 192*pi^3*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 352*pi^4*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^5) + 384*pi^4*δ*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) - 384*pi^4*δ*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^5) + 256*pi^4*δ*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^5)
    end

    if dorder==[2, 0, 0, 1]
        func = -8*pi^2*λ*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^3 + 8*pi^2*λ*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^3 - 16*pi^3*δ*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 + 32*pi^2*λ^2*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) + 48*pi^3*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 48*pi^3*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) + 64*pi^3*δ*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4)
    end

    if dorder==[1, 0, 0, 1]
        func = 4*pi*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^2 - 4*pi^2*δ*λ*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^3 + 4*pi^2*δ*λ*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^3 + 16*pi^2*δ*λ^2*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3)
    end

    if dorder==[0, 0, 0, 1]
        func = 4*pi*δ*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^2
    end

    if dorder==[4, 0, 2, 0]
        func = -160*pi^4*(-1 - α)*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 96*pi^4*(-1 - α)*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) + 768*pi^4*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 448*pi^4*λ*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^6*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 192*pi^4*λ^2*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^8*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4) - 2688*pi^4*λ*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) + 5184*pi^4*λ*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 7872*pi^4*λ^2*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4) + 6016*pi^4*λ^2*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^6*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4) - 5632*pi^4*λ^3*(-4 - α)*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^6*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-5 - α)/(α^4*rq_P^4) + 1536*pi^4*λ^3*(-4 - α)*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^8*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-5 - α)/(α^4*rq_P^4) - 1024*pi^4*λ^4*(-4 - α)*(-3 - α)*(-2 - α)*(-1 - α)*(-5 - α)*sin(pi*δ/rq_P)^8*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-6 - α)/(α^5*rq_P^4)
    end

    if dorder==[3, 0, 2, 0]
        func = -96*pi^3*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) + 160*pi^3*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) - 768*pi^3*λ*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^3) + 448*pi^3*λ*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^3) - 960*pi^3*λ^2*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^3) + 192*pi^3*λ^2*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^7*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^3) - 256*pi^3*λ^3*(-4 - α)*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^7*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-5 - α)/(α^4*rq_P^3)
    end

    if dorder==[2, 0, 2, 0]
        func = 16*pi^2*(-1 - α)*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^2) - 48*pi^2*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^2) + 16*pi^2*λ*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^6*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^2) - 144*pi^2*λ*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^2) - 64*pi^2*λ^2*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^6*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^2)
    end

    if dorder==[1, 0, 2, 0]
        func = -16*pi*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P) - 16*pi*λ*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P)
    end

    if dorder==[0, 0, 2, 0]
        func = -4*(-1 - α)*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/α
    end

    if dorder==[4, 1, 1, 0]
        func = -16*pi^4*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 + 16*pi^4*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 + 96*pi^4*λ*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) + 96*pi^4*λ*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 16*pi^4*λ*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)*(-2*sin(pi*δ/rq_P)^2/(α*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 2*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) + 4*λ*(-1 - α)*sin(pi*δ/rq_P)^4/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2))/rq_P^4 + 16*pi^4*λ*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)*(-2*sin(pi*δ/rq_P)^2/(α*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 2*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) + 4*λ*(-1 - α)*sin(pi*δ/rq_P)^4/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2))/rq_P^4 + 96*pi^4*λ*(-1 - α)*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^4) + 96*pi^4*λ*(-1 - α)*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^4) + 96*pi^4*λ^2*(-2 - α)*sin(pi*δ/rq_P)^6*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 704*pi^4*λ*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) + 96*pi^4*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^6*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^4) - 704*pi^4*λ*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^4) + 1248*pi^4*λ^2*(-2 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 1856*pi^4*λ^2*(-2 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) + 1152*pi^4*λ^2*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 1152*pi^4*λ^2*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 96*pi^4*λ*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 96*pi^4*λ*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 32*pi^4*λ*(-1 - α)*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 48*pi^4*λ^2*(-1 - α)*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)*(-2*sin(pi*δ/rq_P)^2/(α*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 2*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) + 4*λ*(-2 - α)*sin(pi*δ/rq_P)^4/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2))/(α*rq_P^4) - 48*pi^4*λ^2*(-1 - α)*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)*(-2*sin(pi*δ/rq_P)^2/(α*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 2*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) + 4*λ*(-2 - α)*sin(pi*δ/rq_P)^4/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2))/(α*rq_P^4) + 2400*pi^4*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^4) - 3008*pi^4*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^4) + 1792*pi^4*λ^3*(-3 - α)*(-2 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4) - 768*pi^4*λ^3*(-3 - α)*(-2 - α)*sin(pi*δ/rq_P)^6*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4) + 1792*pi^4*λ^3*(-3 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4) - 768*pi^4*λ^3*(-3 - α)*(-1 - α)*sin(pi*δ/rq_P)^6*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4) + 1024*pi^4*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4) - 96*pi^4*λ^2*(-2 - α)*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^6*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) + 704*pi^4*λ*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) + 32*pi^4*λ*(-1 - α)*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) + 352*pi^4*λ^2*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)*(-2*sin(pi*δ/rq_P)^2/(α*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 2*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) + 4*λ*(-2 - α)*sin(pi*δ/rq_P)^4/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2))/(α*rq_P^4) + 4608*pi^4*λ^3*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^4*rq_P^4) - 1536*pi^4*λ^3*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^6*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^4*rq_P^4) + 512*pi^4*λ^4*(-4 - α)*(-3 - α)*(-2 - α)*sin(pi*δ/rq_P)^6*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-5 - α)/(α^4*rq_P^4) + 512*pi^4*λ^4*(-4 - α)*(-3 - α)*(-1 - α)*sin(pi*δ/rq_P)^6*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-5 - α)/(α^4*rq_P^4) + 512*pi^4*λ^4*(-4 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^6*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-5 - α)/(α^4*rq_P^4) - 1152*pi^4*λ^2*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) + 1152*pi^4*λ^2*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 96*pi^4*λ^2*(-2 - α)*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) + 704*pi^4*λ^2*(-2 - α)*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 384*pi^4*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)*(-2*sin(pi*δ/rq_P)^2/(α*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 2*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) + 4*λ*(-3 - α)*sin(pi*δ/rq_P)^4/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2))/(α^2*rq_P^4) + 384*pi^4*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)*(-2*sin(pi*δ/rq_P)^2/(α*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 2*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) + 4*λ*(-3 - α)*sin(pi*δ/rq_P)^4/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2))/(α^2*rq_P^4) + 1536*pi^4*λ^4*(-4 - α)*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^6*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-5 - α)/(α^5*rq_P^4) - 1024*pi^4*λ^3*(-3 - α)*(-2 - α)*(-1 - α)*(-2*λ*(-4 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4) - 768*pi^4*λ^3*(-3 - α)*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4) + 768*pi^4*λ^3*(-3 - α)*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^6*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4) - 256*pi^4*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)*(-2*sin(pi*δ/rq_P)^2/(α*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 2*(-4 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) + 4*λ*(-4 - α)*sin(pi*δ/rq_P)^4/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2))/(α^3*rq_P^4) - 512*pi^4*λ^4*(-4 - α)*(-3 - α)*(-2 - α)*(-1 - α)*(-2*λ*(-4 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^6*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-5 - α)/(α^4*rq_P^4)
    end

    if dorder==[3, 1, 1, 0]
        func = 16*pi^3*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^3 + 96*pi^3*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) - 96*pi^3*λ*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) + 16*pi^3*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)*(-2*sin(pi*δ/rq_P)^2/(α*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 2*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) + 4*λ*(-1 - α)*sin(pi*δ/rq_P)^4/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2))/rq_P^3 + 96*pi^3*λ*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^3) - 96*pi^3*λ*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^3) + 288*pi^3*λ^2*(-2 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^3) - 96*pi^3*λ^2*(-2 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^3) + 192*pi^3*λ^2*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^3) + 480*pi^3*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^3) - 96*pi^3*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^3) + 128*pi^3*λ^3*(-3 - α)*(-2 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^3) + 128*pi^3*λ^3*(-3 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^3) - 96*pi^3*λ*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) + 96*pi^3*λ*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) + 32*pi^3*λ*(-1 - α)*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) - 48*pi^3*λ^2*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)*(-2*sin(pi*δ/rq_P)^2/(α*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 2*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) + 4*λ*(-2 - α)*sin(pi*δ/rq_P)^4/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2))/(α*rq_P^3) + 48*pi^3*λ^2*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)*(-2*sin(pi*δ/rq_P)^2/(α*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 2*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) + 4*λ*(-2 - α)*sin(pi*δ/rq_P)^4/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2))/(α*rq_P^3) + 256*pi^3*λ^3*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^4*rq_P^3) - 192*pi^3*λ^2*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^3) - 96*pi^3*λ^2*(-2 - α)*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^3) + 96*pi^3*λ^2*(-2 - α)*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^3) - 64*pi^3*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)*(-2*sin(pi*δ/rq_P)^2/(α*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 2*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) + 4*λ*(-3 - α)*sin(pi*δ/rq_P)^4/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2))/(α^2*rq_P^3) - 128*pi^3*λ^3*(-3 - α)*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^3)
    end

    if dorder==[2, 1, 1, 0]
        func = 4*pi^2*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^2 - 4*pi^2*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^2 + 4*pi^2*λ*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)*(-2*sin(pi*δ/rq_P)^2/(α*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 2*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) + 4*λ*(-1 - α)*sin(pi*δ/rq_P)^4/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2))/rq_P^2 - 4*pi^2*λ*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)*(-2*sin(pi*δ/rq_P)^2/(α*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 2*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) + 4*λ*(-1 - α)*sin(pi*δ/rq_P)^4/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2))/rq_P^2 + 32*pi^2*λ*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^2) + 32*pi^2*λ*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^2) + 32*pi^2*λ^2*(-2 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^2) + 8*pi^2*λ*(-1 - α)*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^2) + 32*pi^2*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^2) - 32*pi^2*λ*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^2) - 8*pi^2*λ*(-1 - α)*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^2) - 16*pi^2*λ^2*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)*(-2*sin(pi*δ/rq_P)^2/(α*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 2*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) + 4*λ*(-2 - α)*sin(pi*δ/rq_P)^4/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2))/(α*rq_P^2) - 32*pi^2*λ^2*(-2 - α)*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^2)
    end

    if dorder==[1, 1, 1, 0]
        func = -4*pi*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P - 4*pi*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)*(-2*sin(pi*δ/rq_P)^2/(α*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 2*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) + 4*λ*(-1 - α)*sin(pi*δ/rq_P)^4/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2))/rq_P - 8*pi*λ*(-1 - α)*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P)
    end

    if dorder==[0, 1, 1, 0]
        func = -2*sin(pi*δ/rq_P)^2*(2*λ*sin(pi*δ/rq_P)^2/(α*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α) - 4*λ*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/α^2
    end

    if dorder==[4, 0, 1, 0]
        func = -16*pi^4*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 + 16*pi^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 - 128*pi^4*λ*(-1 - α)*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 96*pi^4*λ*(-1 - α)*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 96*pi^4*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^6*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) + 736*pi^4*λ*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 1248*pi^4*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) + 1856*pi^4*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 1792*pi^4*λ^3*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4) + 768*pi^4*λ^3*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^6*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4) - 512*pi^4*λ^4*(-4 - α)*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^6*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-5 - α)/(α^4*rq_P^4)
    end

    if dorder==[3, 0, 1, 0]
        func = 16*pi^3*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^3 - 96*pi^3*λ*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) + 128*pi^3*λ*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) - 288*pi^3*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^3) + 96*pi^3*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^3) - 128*pi^3*λ^3*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^5*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^3)
    end

    if dorder==[2, 0, 1, 0]
        func = 4*pi^2*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^2 - 4*pi^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^2 + 8*pi^2*λ*(-1 - α)*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^2) - 40*pi^2*λ*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^2) - 32*pi^2*λ^2*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^2)
    end

    if dorder==[1, 0, 1, 0]
        func = -4*pi*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P - 8*pi*λ*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P)
    end

    if dorder==[0, 0, 1, 0]
        func = -2*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)
    end

    if dorder==[4, 2, 0, 0]
        func = -96*pi^4*λ^2*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^4) - 96*pi^4*λ^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^4) - 16*pi^4*λ*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))^2*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 + 16*pi^4*λ*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 - 16*pi^4*λ*(4*λ*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 4*λ^2*(-1 - α)*sin(pi*δ/rq_P)^4/(α^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 + 16*pi^4*λ*(4*λ*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 4*λ^2*(-1 - α)*sin(pi*δ/rq_P)^4/(α^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 - 96*pi^4*λ^2*(-1 - α)*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^3*rq_P^4) - 96*pi^4*λ^2*(-1 - α)*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^3*rq_P^4) + 704*pi^4*λ^2*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^4) - 768*pi^4*λ^3*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) + 768*pi^4*λ^3*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) + 96*pi^4*λ^2*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) + 96*pi^4*λ^2*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) + 704*pi^4*λ^2*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^3*rq_P^4) - 1536*pi^4*λ^3*(-2 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^4) + 1536*pi^4*λ^3*(-2 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^4) - 1536*pi^4*λ^3*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^4) + 1536*pi^4*λ^3*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^4) - 512*pi^4*λ^4*(-3 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4) - 512*pi^4*λ^4*(-2 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4) - 512*pi^4*λ^4*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4) + 96*pi^4*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^4) + 96*pi^4*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^4) - 704*pi^4*λ^2*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 48*pi^4*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))^2*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 48*pi^4*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 48*pi^4*λ^2*(-1 - α)*(4*λ*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 4*λ^2*(-2 - α)*sin(pi*δ/rq_P)^4/(α^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 48*pi^4*λ^2*(-1 - α)*(4*λ*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 4*λ^2*(-2 - α)*sin(pi*δ/rq_P)^4/(α^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 2304*pi^4*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^4*rq_P^4) + 2304*pi^4*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^4*rq_P^4) - 1536*pi^4*λ^4*(-3 - α)*(-2 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^4*rq_P^4) - 1536*pi^4*λ^4*(-3 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^4*rq_P^4) - 1536*pi^4*λ^4*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^4*rq_P^4) - 704*pi^4*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^4) + 768*pi^4*λ^3*(-2 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 768*pi^4*λ^3*(-2 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) + 768*pi^4*λ^3*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 768*pi^4*λ^3*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) + 352*pi^4*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))^2*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) + 352*pi^4*λ^2*(-1 - α)*(4*λ*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 4*λ^2*(-2 - α)*sin(pi*δ/rq_P)^4/(α^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 3072*pi^4*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^5*rq_P^4) + 1536*pi^4*λ^3*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^4) - 1536*pi^4*λ^3*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^4) + 512*pi^4*λ^4*(-3 - α)*(-2 - α)*(-2*λ*(-4 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4) + 512*pi^4*λ^4*(-3 - α)*(-1 - α)*(-2*λ*(-4 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4) + 512*pi^4*λ^4*(-2 - α)*(-1 - α)*(-2*λ*(-4 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4) - 384*pi^4*λ^3*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))^2*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) + 384*pi^4*λ^3*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))^2*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 384*pi^4*λ^3*(-2 - α)*(-1 - α)*(4*λ*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 4*λ^2*(-3 - α)*sin(pi*δ/rq_P)^4/(α^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) + 384*pi^4*λ^3*(-2 - α)*(-1 - α)*(4*λ*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 4*λ^2*(-3 - α)*sin(pi*δ/rq_P)^4/(α^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) + 1536*pi^4*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*(-2*λ*(-4 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^4*rq_P^4) - 256*pi^4*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*(-2*λ*(-4 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))^2*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4) - 256*pi^4*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*(4*λ*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 4*λ^2*(-4 - α)*sin(pi*δ/rq_P)^4/(α^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*λ*(-4 - α)*sin(pi*δ/rq_P)^2/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4)
    end

    if dorder==[3, 2, 0, 0]
        func = -96*pi^3*λ^2*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^3) + 96*pi^3*λ^2*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^3) - 128*pi^3*λ^3*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^3) + 16*pi^3*λ*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))^2*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^3 + 16*pi^3*λ*(4*λ*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 4*λ^2*(-1 - α)*sin(pi*δ/rq_P)^4/(α^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^3 - 96*pi^3*λ^2*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^3*rq_P^3) + 96*pi^3*λ^2*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^3*rq_P^3) - 256*pi^3*λ^3*(-2 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^3) - 256*pi^3*λ^3*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^3) + 96*pi^3*λ^2*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) - 96*pi^3*λ^2*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) - 384*pi^3*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^4*rq_P^3) + 96*pi^3*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^3) - 96*pi^3*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^3) + 128*pi^3*λ^3*(-2 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^3) + 128*pi^3*λ^3*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^3) - 48*pi^3*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))^2*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) + 48*pi^3*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))^2*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) - 48*pi^3*λ^2*(-1 - α)*(4*λ*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 4*λ^2*(-2 - α)*sin(pi*δ/rq_P)^4/(α^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) + 48*pi^3*λ^2*(-1 - α)*(4*λ*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 4*λ^2*(-2 - α)*sin(pi*δ/rq_P)^4/(α^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) + 256*pi^3*λ^3*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^3) - 64*pi^3*λ^3*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))^2*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^3) - 64*pi^3*λ^3*(-2 - α)*(-1 - α)*(4*λ*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 4*λ^2*(-3 - α)*sin(pi*δ/rq_P)^4/(α^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^3)
    end

    if dorder==[2, 2, 0, 0]
        func = 4*pi^2*λ*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))^2*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^2 - 4*pi^2*λ*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^2 + 4*pi^2*λ*(4*λ*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 4*λ^2*(-1 - α)*sin(pi*δ/rq_P)^4/(α^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^2 - 4*pi^2*λ*(4*λ*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 4*λ^2*(-1 - α)*sin(pi*δ/rq_P)^4/(α^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^2 - 32*pi^2*λ^2*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^2) - 32*pi^2*λ^2*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^3*rq_P^2) + 32*pi^2*λ^2*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^2) + 32*pi^2*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^2) - 16*pi^2*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))^2*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^2) - 16*pi^2*λ^2*(-1 - α)*(4*λ*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 4*λ^2*(-2 - α)*sin(pi*δ/rq_P)^4/(α^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^2)
    end

    if dorder==[1, 2, 0, 0]
        func = -4*pi*λ*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))^2*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P - 4*pi*λ*(4*λ*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - 4*λ^2*(-1 - α)*sin(pi*δ/rq_P)^4/(α^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^2) + 4*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P
    end

    if dorder==[0, 2, 0, 0]
        func = (2*λ*sin(pi*δ/rq_P)^2/(α*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-α) + 4*λ^2*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/α^3
    end

    if dorder==[4, 1, 0, 0]
        func = 48*pi^4*λ^2*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) + 48*pi^4*λ^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 16*pi^4*λ*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 + 16*pi^4*λ*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 + 48*pi^4*λ^2*(-1 - α)*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^4) + 48*pi^4*λ^2*(-1 - α)*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^4) - 352*pi^4*λ^2*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 352*pi^4*λ^2*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^4) + 384*pi^4*λ^3*(-2 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 384*pi^4*λ^3*(-2 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) + 384*pi^4*λ^3*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 384*pi^4*λ^3*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 48*pi^4*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 48*pi^4*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) + 768*pi^4*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^4) - 768*pi^4*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^4) + 256*pi^4*λ^4*(-3 - α)*(-2 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4) + 256*pi^4*λ^4*(-3 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4) + 256*pi^4*λ^4*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4) + 352*pi^4*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) + 768*pi^4*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^4*rq_P^4) - 384*pi^4*λ^3*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) + 384*pi^4*λ^3*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 256*pi^4*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*(-2*λ*(-4 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4)
    end

    if dorder==[3, 1, 0, 0]
        func = 48*pi^3*λ^2*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) - 48*pi^3*λ^2*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) + 16*pi^3*λ*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^3 + 48*pi^3*λ^2*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^3) - 48*pi^3*λ^2*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^3) + 64*pi^3*λ^3*(-2 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^3) + 64*pi^3*λ^3*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^3) + 128*pi^3*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^3*rq_P^3) - 48*pi^3*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) + 48*pi^3*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) - 64*pi^3*λ^3*(-2 - α)*(-1 - α)*(-2*λ*(-3 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^3)
    end

    if dorder==[2, 1, 0, 0]
        func = 4*pi^2*λ*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^2 - 4*pi^2*λ*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^2 + 16*pi^2*λ^2*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^2) + 16*pi^2*λ^2*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α^2*rq_P^2) - 16*pi^2*λ^2*(-1 - α)*(-2*λ*(-2 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^2)
    end

    if dorder==[1, 1, 0, 0]
        func = -4*pi*λ*(-2*λ*(-1 - α)*sin(pi*δ/rq_P)^2/(α^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P
    end

    if dorder==[0, 1, 0, 0]
        func = (2*λ*sin(pi*δ/rq_P)^2/(α*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)) - log(1 + 2*λ*sin(pi*δ/rq_P)^2/α))*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-α)
    end

    if dorder==[4, 0, 0, 0]
        func = -16*pi^4*λ*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 + 16*pi^4*λ*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^4 - 48*pi^4*λ^2*(-1 - α)*sin(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 48*pi^4*λ^2*(-1 - α)*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) + 352*pi^4*λ^2*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^4) - 384*pi^4*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) + 384*pi^4*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^4) - 256*pi^4*λ^4*(-3 - α)*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^4*cos(pi*δ/rq_P)^4*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-4 - α)/(α^3*rq_P^4)
    end

    if dorder==[3, 0, 0, 0]
        func = 16*pi^3*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^3 - 48*pi^3*λ^2*(-1 - α)*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) + 48*pi^3*λ^2*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^3) - 64*pi^3*λ^3*(-2 - α)*(-1 - α)*sin(pi*δ/rq_P)^3*cos(pi*δ/rq_P)^3*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-3 - α)/(α^2*rq_P^3)
    end

    if dorder==[2, 0, 0, 0]
        func = 4*pi^2*λ*sin(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^2 - 4*pi^2*λ*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P^2 - 16*pi^2*λ^2*(-1 - α)*sin(pi*δ/rq_P)^2*cos(pi*δ/rq_P)^2*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-2 - α)/(α*rq_P^2)
    end

    if dorder==[1, 0, 0, 0]
        func = -4*pi*λ*sin(pi*δ/rq_P)*cos(pi*δ/rq_P)*(1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-1 - α)/rq_P
    end

    if dorder==[0, 0, 0, 0]
        func = (1 + 2*λ*sin(pi*δ/rq_P)^2/α)^(-α)
    end

    return even_time_derivative * float(func)  # correcting for amount of t2 derivatives

end

return rq_periodic_kernel, 3  # the function handle and the number of kernel hyperparameters


text:
# Define abstract interface
"""
$(TYPEDEF)

Abstract supertype for quadrature nodes.
"""
abstract type QuadratureNodes{T} end

"""
$(SIGNATURES)

Implements the julia iteration interface returning (node,weight) tuples
while iterating over `quadratureNodes`.
"""
function Base.iterate(quadratureNodes::QuadratureNodes{T}) where {T}
    error("Base.iterate is not implemented for quadratureNodes::$(typeof(quadratureNodes)).")
end

function Base.iterate(quadratureNodes::QuadratureNodes{T}, state) where {T}
    error("Base.iterate is not implemented for quadratureNodes::$(typeof(quadratureNodes)).")
end

# Define specific subtypes
## gauss legendre nodes provide very good results in case, where the magnetization is continous on the integration domain
struct GaussLegendre{T} <: QuadratureNodes{T}
    nodesx::Vector{T}
    weightsx::Vector{T}
    nodesy::Vector{T}
    weightsy::Vector{T}
    nodesz::Vector{T}
    weightsz::Vector{T}
    n::Tuple{Int,Int,Int}
end

function Base.iterate(quadratureNodes::GaussLegendre{T}, state = (CartesianIndices(quadratureNodes.n),1)) where {T}
    (cartesian,counter) = state

    if counter>length(cartesian)
        return nothing
    else
        i,j,k = Tuple(cartesian[counter])
        r = SVector{3,T}(quadratureNodes.nodesx[i],quadratureNodes.nodesy[j],quadratureNodes.nodesz[k])
        w = quadratureNodes.weightsx[i]*quadratureNodes.weightsy[j]*quadratureNodes.weightsz[k]
        return ((r,w),(cartesian,counter+1))
    end
end

## midpoint nodes are the fall back option, where the magnetization is non-continous on the integration domain
struct MidPoint{T,U<:AbstractVector{T},V<:AbstractVector{T},W<:AbstractVector{T}} <: QuadratureNodes{T}
    nodesx::U
    nodesy::V
    nodesz::W
    weight::T
    n::Tuple{Int,Int,Int}
end

function Base.iterate(quadratureNodes::MidPoint{T}, state = (CartesianIndices(quadratureNodes.n),1)) where {T}
    (cartesian,counter) = state

    if counter>length(cartesian)
        return nothing
    else
        i,j,k = Tuple(cartesian[counter])
        r = SVector{3,T}(quadratureNodes.nodesx[i],quadratureNodes.nodesy[j],quadratureNodes.nodesz[k])
        return ((r,quadratureNodes.weight),(cartesian,counter+1))
    end
end

function MidPoint(lfb::SVector{3,T},rrt::SVector{3,T},density=10/1e-3) where {T}
    nx,ny,nz = round.(Int,abs.(density*(rrt-lfb)))
    nx = max(nx,1)
    ny = max(ny,1)
    nz = max(nz,1)
    n = tuple(nx,ny,nz)

    center = (rrt+lfb)/2
    halfLength = abs.((rrt-lfb)/2)
    nodesx = range(center[1]-(1-1/nx)*halfLength[1], center[1]+(1-1/nx)*halfLength[1], length=nx)
    nodesy = range(center[2]-(1-1/ny)*halfLength[2], center[2]+(1-1/ny)*halfLength[2], length=ny)
    nodesz = range(center[3]-(1-1/nz)*halfLength[3], center[3]+(1-1/nz)*halfLength[3], length=nz)
    weight = 2*halfLength[1] * 2*halfLength[2] * 2*halfLength[3] / prod(n)

    return MidPoint{T,typeof(nodesx),typeof(nodesy),typeof(nodesz)}(nodesx,nodesy,nodesz,weight,n)
end

ext: jl


text:
using Flux, DiffEqFlux, DiffEqSensitivity, DelayDiffEq, Zygote, Test

## Setup DDE to optimize
function delay_lotka_volterra(du,u,h,p,t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = (α - β*y)*h(p,t-0.1)[1]
  du[2] = dy = (δ*x - γ)*y
end
h(p,t) = ones(eltype(p),2)
prob = DDEProblem(delay_lotka_volterra,[1.0,1.0],h,(0.0,10.0),constant_lags=[0.1])
p = [2.2, 1.0, 2.0, 0.4]
function predict_fd_dde(p)
  concrete_solve(prob,MethodOfSteps(Tsit5()),prob.u0,p,saveat=0.0:0.1:10.0,reltol=1e-4,sensealg=ForwardDiffSensitivity())[1,:]
end
loss_fd_dde(p) = sum(abs2,x-1 for x in predict_fd_dde(p))
loss_fd_dde(p)
@test !iszero(Zygote.gradient(loss_fd_dde,p)[1])

function predict_rd_dde(p)
  concrete_solve(prob,MethodOfSteps(Tsit5()),prob.u0,p,saveat=0.1,reltol=1e-4,sensealg=TrackerAdjoint())[1,:]
end
loss_rd_dde(p) = sum(abs2,x-1 for x in predict_rd_dde(p))
loss_rd_dde(p)
@test !iszero(Zygote.gradient(loss_rd_dde,p)[1])

@test Zygote.gradient(loss_fd_dde,p)[1] ≈ Zygote.gradient(loss_rd_dde,p)[1] rtol=1e-2


text:
# This implements the ScikitLearn.jl interface and makes GMM compatible with
# ScikitLearn pipelines

import ScikitLearnBase

n_components(gmm::GMM) = size(gmm.μ, 1)

# We set d=1, but that's OK, we'll change it in fit!
GMM(; n_components=1, kind=:diag) = GMM(n_components, 2; kind=kind)
ScikitLearnBase.is_classifier(::GMM) = false
ScikitLearnBase.clone(gmm::GMM) =
    GMM(; n_components=n_components(gmm),
        kind=kind(gmm))

function ScikitLearnBase.get_params(gmm::GMM)
    return Dict(:n_components=>n_components(gmm))
end

function ScikitLearnBase.set_params!(gmm::GMM; params...)
    for (param, val) in params
        if param==:n_components
            # We don't have to update Σ, it'll be wiped in fit! anyway
            gmm.μ = zeros(val, 1)
        else
            throw(ArgumentError("Bad parameter: $param"))
        end
    end
    gmm
end

function Base.copy!(gmm_dest::GMM, gmm_src::GMM)
    # shallow copy - used below
    for f in fieldnames(gmm_dest)
        setfield!(gmm_dest, f, getfield(gmm_src, f))
    end
end

function ScikitLearnBase.fit!(gmm::GMM, X::AbstractMatrix, y=nothing)
    n = n_components(gmm)
    # Creating a temporary is not great, but it's negligible in the grand
    # scheme of thing.  We'd just need a slight refactor in
    # GaussianMixtures/src/train.jl to avoid it
    gmm_temp = GMM(n, X; kind=kind(gmm))
    copy!(gmm, gmm_temp)
    gmm
end

ScikitLearnBase.predict_log_proba(gmm::GMM, X) = log(gmmposterior(gmm, X)[1])
ScikitLearnBase.predict_proba(gmm::GMM, X) = gmmposterior(gmm, X)[1]
ScikitLearnBase.predict(gmm::GMM, X) =
    # This is just `argmax(axis=2)`. It's very verbose in Julia.
    ind2sub(size(X),
            vec(findmax(ScikitLearnBase.predict_proba(gmm, X), 2)[2]))[2]

""" `density(gmm::GMM, X)` returns `log(P(X|μ, Σ))` """
function density(gmm::GMM, X)
    # Let mᵢ be "X came from mixture #i"
    # P(X|μ, Σ) = P(X|m₁, μ, Σ) * P(mᵢ) + P(X|m₂, μ, Σ) * P(m₂) + ...
    logPrior = reshape(log(gmm.w), 1, length(gmm.w))
    PX = logsumexp(broadcast(+, llpg(gmm, X), logPrior), 2)
    return squeeze(PX, 2)::Vector
end

# score_samples is underspecified by the scikit-learn API, so we're more or
# less free to return what we want
# http://scikit-learn.org/stable/modules/generated/sklearn.mixture.GMM.html#sklearn.mixture.GMM
# http://scikit-learn.org/stable/modules/generated/sklearn.neighbors.KernelDensity.html
"""    score_samples(gmm::GMM, X::Matrix)
Return the per-sample likelihood of the data under the model. """
ScikitLearnBase.score_samples(gmm::GMM, X) = density(gmm, X)

# Used for cross-validation. Higher is better.
ScikitLearnBase.score(gmm::GMM, X) = avll(gmm, X)