Use FastGaussQuadrature for transform between Position and Fock bases (…

…#142)

Project.toml

@@ -6,6 +6,7 @@ version = "0.4.15"
 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"

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

test/runtests.jl

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