diff --git a/Project.toml b/Project.toml index 4ca150fc..9aa04bf8 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "0.4.14" Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" FastExpm = "7868e603-8603-432e-a1a1-694bd70b01f2" +FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -19,6 +20,7 @@ UnsafeArrays = "c4a57d5a-5b31-53a6-b365-19f8c011fbd6" Adapt = "1, 2, 3.3" FFTW = "1.2" FastExpm = "1.1.0" +FastGaussQuadrature = "0.5" FillArrays = "0.13, 1" LRUCache = "1" QuantumInterface = "0.3.0" diff --git a/src/QuantumOpticsBase.jl b/src/QuantumOpticsBase.jl index 1ee5694e..14608f2e 100644 --- a/src/QuantumOpticsBase.jl +++ b/src/QuantumOpticsBase.jl @@ -65,7 +65,6 @@ export Basis, GenericBasis, CompositeBasis, basis, #apply apply! -include("polynomials.jl") include("bases.jl") include("states.jl") include("operators.jl") diff --git a/src/polynomials.jl b/src/polynomials.jl deleted file mode 100644 index ef5615c5..00000000 --- a/src/polynomials.jl +++ /dev/null @@ -1,96 +0,0 @@ -""" - horner(coefficients, x) - -Evaluate the given polynomial at position x using the Horner scheme. - -```math -p(x) = \\sum_{n=0}^N c_n x^n -``` -""" -function horner(coefficients, x) - bn = coefficients[end] - for n=length(coefficients)-1:-1:1 - bn = coefficients[n] + bn*x - end - bn -end - - -module hermite - -""" - a_nk(N) - -Calculate the all coefficients for all Hermite polynomials up to order N. - -```math -H_n(x) = \\sum_{k=0}^n a_{n,k} x^k -``` - -Returns a vector of length N+1 where the n-th entry contains all coefficients -for the n-th Hermite polynomial. -""" -function a(N::T) where T - a = Vector{Vector{T}}(undef, N+1) - a[1] = [1] - a[2] = [0,2] - am = a[2] - for n=2:N - an = zeros(Int, n+1) - a[n+1] = an - if iseven(n) - an[1] = -am[2] - end - an[n+1] = 2*am[n] - if iseven(n) - for k=3:2:n-1 - an[k] = 2*am[k-1] - am[k+1]*k - end - else - for k=2:2:n-1 - an[k] = 2*am[k-1] - am[k+1]*k - end - end - am = an - end - a -end - -""" - A_nk(N) - -Calculate the all scaled coefficients for all Hermite polynomials up to order N. - -The scaled coefficients `A` are connected to the unscaled coefficients `a` by -the relation ``A_{n,k} = \\frac{a_{n,k}}{\\sqrt{2^n n!}}`` - -Returns a vector of length N+1 where the n-th entry contains all scaled -coefficients for the n-th Hermite polynomial. -""" -function A(N::T) where T - A = Vector{Vector{float(T)}}(undef, N+1) - A[1] = [1.] - A[2] = [0., sqrt(2)] - Am = A[2] - for n=2:N - An = zeros(float(T), n+1) - A[n+1] = An - if iseven(n) - An[1] = -Am[2]/sqrt(2*n) - end - An[n+1] = Am[n]*sqrt(2/n) - if iseven(n) - for k=3:2:n-1 - An[k] = Am[k-1]*sqrt(2/n) - Am[k+1]*k/sqrt(2*n) - end - else - for k=2:2:n-1 - An[k] = Am[k-1]*sqrt(2/n) - Am[k+1]*k/sqrt(2*n) - end - end - Am = An - end - A -end - -end # module hermite diff --git a/src/transformations.jl b/src/transformations.jl index 9f0b86c5..bbba8813 100644 --- a/src/transformations.jl +++ b/src/transformations.jl @@ -1,3 +1,5 @@ +import FastGaussQuadrature: hermpoly_rec + """ transform([S=ComplexF64, ]b1::PositionBasis, b2::FockBasis; x0=1) transform([S=ComplexF64, ]b1::FockBasis, b2::PositionBasis; x0=1) @@ -15,38 +17,17 @@ a harmonic trap potential at position ``x``, i.e.: \\frac{1}{\\sqrt{2^n n!}} H_n\\left(\\frac{x}{x_0}\\right) ``` """ -function transform(::Type{S}, b1::PositionBasis, b2::FockBasis; x0=1) where S - T = Matrix{S}(undef, length(b1), length(b2)) - xvec = samplepoints(b1) - A = hermite.A(b2.N) - delta_x = spacing(b1) - c = 1.0/sqrt(x0*sqrt(pi))*sqrt(delta_x) - for i in 1:length(b1) - u = xvec[i]/x0 - C = c*exp(-u^2/2) - for n=b2.offset:b2.N - idx = n-b2.offset+1 - T[i,idx] = C*horner(A[n+1], u) - end +function transform(::Type{S}, bp::PositionBasis, bf::FockBasis; x0=1) where S + T = Matrix{S}(undef, length(bp), length(bf)) + xvec = samplepoints(bp) + C = pi^(-1/4)*sqrt(spacing(bp)/x0) + for i in 1:length(bp) + T[i,:] = C*hermpoly_rec(bf.offset:bf.N, sqrt(2)*xvec[i]/x0) end - DenseOperator(b1, b2, T) + DenseOperator(bp, bf, T) end -function transform(::Type{S}, b1::FockBasis, b2::PositionBasis; x0=1) where S - T = Matrix{S}(undef, length(b1), length(b2)) - xvec = samplepoints(b2) - A = hermite.A(b1.N) - delta_x = spacing(b2) - c = 1.0/sqrt(x0*sqrt(pi))*sqrt(delta_x) - for i in 1:length(b2) - u = xvec[i]/x0 - C = c*exp(-u^2/2) - for n=b1.offset:b1.N - idx = n-b1.offset+1 - T[idx,i] = C*horner(A[n+1], u) - end - end - DenseOperator(b1, b2, T) -end +transform(::Type{S}, bf::FockBasis, bp::PositionBasis; x0=1) where S = + dagger(transform(S, bp, bf; x0=x0)) transform(b1::Basis,b2::Basis;kwargs...) = transform(ComplexF64,b1,b2;kwargs...) diff --git a/test/runtests.jl b/test/runtests.jl index 5a719baf..634a8bc9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,5 @@ names = [ "test_sortedindices.jl", - "test_polynomials.jl", "test_bases.jl", "test_states.jl", diff --git a/test/test_polynomials.jl b/test/test_polynomials.jl deleted file mode 100644 index b9299a37..00000000 --- a/test/test_polynomials.jl +++ /dev/null @@ -1,28 +0,0 @@ -using Test -using QuantumOpticsBase: horner, hermite - -@testset "polynomials" begin - -# Test Horner scheme -c = [0.2, 0.6, 1.7] -x0 = 1.3 -@test horner(c, x0) == c[1] + c[2]*x0 + c[3]*x0^2 - -# Test Hermite polynomials -an = Vector{Vector{Int}}(undef, 8) -an[1] = [1] -an[2] = [0,2] -an[3] = [-2, 0, 4] -an[4] = [0, -12, 0, 8] -an[5] = [12, 0, -48, 0, 16] -an[6] = [0, 120, 0, -160, 0, 32] -an[7] = [-120, 0, 720, 0, -480, 0, 64] -an[8] = [0, -1680, 0, 3360, 0, -1344, 0, 128] -@test hermite.a(7) == an - -A = hermite.A(7) -for n=0:7 - @test A[n+1] ≈ an[n+1]/sqrt(2^n*factorial(n)) -end - -end # testset