Use ClassicalOrthogonalPolynomials for transform between Position and…

… Fock bases

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"

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")

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

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...)

test/runtests.jl

@@ -1,6 +1,5 @@
 names = [
     "test_sortedindices.jl",
-    "test_polynomials.jl",
 
     "test_bases.jl",
     "test_states.jl",

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)