Skip to content

Instantly share code, notes, and snippets.

@ikirill

ikirill/Q.jl Secret

Created Jun 23, 2018
Embed
What would you like to do?
module Q
using Cubature
using QuadGK
using GSL
function fi::Real, α::Real, E::Number, k4::Real, kp::Real)
p1 = (k4^2 + E/2)^2 + kp^2 + α^2
p2 = (k4^2 - E/2)^2 + kp^2 + (1-α)^2
2 * 2kp *
* (1+p1)^1.5ξ * (1+p2)^0.5ξ) /
((1+p1*(1+p1)^2ξ) * (1+p2*(1+p2)^2ξ))
end
function f::Real, α::Real, E::Real)
function ii(x)
k4, kp = x[1]/(1-x[1]), x[2]/(1-x[2])
# @show x, k4, kp
fi(ξ, α, E, k4, kp) / ((1-x[1]) * (1-x[2]))^2
end
val, err = hcubature(ii, (0.0, 0.0), (1.0, 1.0), reltol=0, abstol=1e-8, maxevals=10^6)
# @show val, err
val
end
function f::Real, α::Real, E::Complex)
function ii(x, out)
k4, kp = x[1]/(1-x[1]), x[2]/(1-x[2])
# @show x, k4, kp
w = fi(ξ, α, E, k4, kp) / ((1-x[1]) * (1-x[2]))^2
out[1] = real(w)
out[2] = imag(w)
end
val, err = hcubature(2, ii, (0.0, 0.0), (1.0, 1.0), reltol=0, abstol=1e-8, maxevals=10^6)
# @show val, err
val[1] + val[2] * im
end
function g1::Real, α::Real, t::Real, e0::Real)
const θ = 0
γ, δ = cis(θ), cis- θ)
quadgk(u ->
γ * cis(t * γ * u) / ((γ*u)^2 + e0^2) * f(ξ, α, γ*u) -
δ * cis(t * δ * u) / ((δ*u)^2 + e0^2) * f(ξ, α, δ*u),
0, Inf)
end
function gslfun(x::Cdouble, params::Ptr{Void})
(unsafe_pointer_to_objref(params)::Function)(x)::Cdouble
end
gslfun_c = cfunction(gslfun, Cdouble, (Cdouble, Ptr{Void}))
Base.convert(::Type{GSL.gsl_function}, f::Function) = GSL.gsl_function(gslfun_c, pointer_from_objref(f))
function g2::Real, α::Real, t::Real, e0::Real)
limit = 1024
ws1 = GSL.integration_workspace_alloc(limit)
ws2 = GSL.integration_workspace_alloc(1024)
qawo_cos = ccall((:gsl_integration_qawo_table_alloc, GSL.libgsl),
Ptr{GSL.gsl_integration_qawo_table},
(Cdouble, Cdouble, Cint, Csize_t),
t, 1.0, 0, 10)
qawo_sin = ccall((:gsl_integration_qawo_table_alloc, GSL.libgsl),
Ptr{GSL.gsl_integration_qawo_table},
(Cdouble, Cdouble, Cint, Csize_t),
t, 1.0, 1, 10)
f1 = E -> (f(ξ, α, E) + f(ξ, α, -E))/(E^2+e0^2)
f2 = E -> (f(ξ, α, E) - f(ξ, α, -E))/(E^2+e0^2)
result1, result2, abserr = Ref(NaN), Ref(NaN), Ref(NaN)
ret = ccall((:gsl_integration_qawf, GSL.libgsl), Cint,
(Ref{gsl_function}, Cdouble, Cdouble, Csize_t, Ptr{GSL.gsl_integration_workspace}, Ptr{GSL.gsl_integration_workspace}, Ptr{GSL.gsl_integration_qawo_table}, Ref{Cdouble}, Ref{Cdouble}),
f1, 0.0, 1e-8, limit, ws1, ws2, qawo_cos, result1, abserr)
@show ret
@show result1, abserr
ret == 0 || warn("FAILED: gsl_integration_qawf")
ret = ccall((:gsl_integration_qawf, GSL.libgsl), Cint,
(Ref{gsl_function}, Cdouble, Cdouble, Csize_t, Ptr{GSL.gsl_integration_workspace}, Ptr{GSL.gsl_integration_workspace}, Ptr{GSL.gsl_integration_qawo_table}, Ref{Cdouble}, Ref{Cdouble}),
f2, 0.0, 1e-8, limit, ws1, ws2, qawo_sin, result2, abserr)
@show ret
@show result2, abserr
ret == 0 || warn("FAILED: gsl_integration_qawf")
result1[] + im * result2[]
end
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment