diff --git a/Project.toml b/Project.toml index 4ca150fc..921c5a60 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ version = "0.4.14" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" +ClassicalOrthogonalPolynomials = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" FastExpm = "7868e603-8603-432e-a1a1-694bd70b01f2" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" 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/metrics.jl b/src/metrics.jl index 721f8852..5ac00334 100644 --- a/src/metrics.jl +++ b/src/metrics.jl @@ -201,7 +201,7 @@ function ptranspose(rho::DenseOpType{B,B}, indices=1) where B<:CompositeBasis # works as long as QuantumOptics.jl doesn't change the implementation of `tensor`, i.e. tensor(a,b).data = kron(b.data,a.data) nsys = length(rho.basis_l.shape) mask = ones(Int, nsys) - mask[collect(indices)] .+= 1 + mask[vec(collect(indices))] .+= 1 pt_dims = reshape(1:2*nsys, (nsys,2)) # indices of the operator viewed as a tensor with 2nsys legs pt_idx = [[pt_dims[i,mask[i]] for i = 1 : nsys]; [pt_dims[i,3-mask[i]] for i = 1 : nsys] ] # permute the legs on the subsystem of `indices` # reshape the operator data into a 2nsys-legged tensor and shape it back with the legs permuted 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..04776985 100644 --- a/src/transformations.jl +++ b/src/transformations.jl @@ -1,3 +1,4 @@ +import ClassicalOrthogonalPolynomials: Hermite, Normalized, OrthonormalWeighted """ transform([S=ComplexF64, ]b1::PositionBasis, b2::FockBasis; x0=1) transform([S=ComplexF64, ]b1::FockBasis, b2::PositionBasis; x0=1) @@ -15,38 +16,14 @@ 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 - end - DenseOperator(b1, b2, T) +function transform(::Type{S}, bp::PositionBasis, bf::FockBasis) where S + xvec = samplepoints(bp) + C = sqrt(spacing(bp)) + H = OrthonormalWeighted(Normalized(Hermite())) + T = H[xvec, (bf.offset+1):(bf.N+1)] + DenseOperator(bp, bf, convert(Matrix{S}, C*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) where S = dagger(transform(S, bp, bf)) transform(b1::Basis,b2::Basis;kwargs...) = transform(ComplexF64,b1,b2;kwargs...) diff --git a/test/runtests.jl b/test/runtests.jl index 92456a83..8a9cbd94 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 diff --git a/test/test_transformations.jl b/test/test_transformations.jl index f05483e5..df873b66 100644 --- a/test/test_transformations.jl +++ b/test/test_transformations.jl @@ -24,20 +24,6 @@ psi_n = coherentstate(b_fock, α0) psi_x = gaussianstate(b_position, x0, p0, 1) @test 1e-10 > D(psi_x, Txn*psi_n) -# Test different characteristic length -x0 = 0.0 -p0 = 0.2 -α0 = (x0 + 1im*p0)/sqrt(2) -σ0 = 0.7 -Txn = transform(b_position, b_fock; x0=σ0) -Tnx = transform(b_fock, b_position; x0=σ0) -@test 1e-10 > D(dagger(Txn), Tnx) -@test 1e-10 > D(one(b_fock), Tnx*Txn) - -psi_n = coherentstate(b_fock, α0) -psi_x = gaussianstate(b_position, x0/σ0, p0/σ0, σ0) -@test 1e-10 > D(psi_x, Txn*psi_n) - # Test with offset in FockBasis b_fock = FockBasis(50,1) b_position = PositionBasis(-10, 10, 300)