Use FastGaussQuadrature for transform between Position and Fock bases

qojulia · Sep 14, 2023 · 2e6b9da · 2e6b9da
1 parent 032dd82
commit 2e6b9da
Show file tree

Hide file tree

Showing 7 changed files with 12 additions and 170 deletions.
diff --git 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
@@ -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
diff --git 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,16 @@ 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) where S
+    T = Matrix{S}(undef, length(bp), length(bf))
+    xvec = samplepoints(bp)
+    C = pi^(-1/4)*sqrt(spacing(bp))
+    for i in 1:length(bp)
+        T[i,:] = C*hermpoly_rec(bf.offset:bf.N, sqrt(2)*xvec[i])
     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) 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
@@ -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
diff --git 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)