From 3eaad76928dd1484c93c196702bc12782ecb25e7 Mon Sep 17 00:00:00 2001 From: Matias Bundgaard-Nielsen <56378038+mabuni1998@users.noreply.github.com> Date: Sun, 18 Aug 2024 16:52:20 +0200 Subject: [PATCH] WaveguideInteraction loop bugfix and test update Looping now also works for WaveguideInteraction operators and more test for WaveguideInteraction with delayed operators has been hadded --- docs/src/index.md | 2 +- src/WaveguideInteraction.jl | 143 ++++++++++++++++++++++++++++-------- test/test_operators.jl | 44 +++++++++++ test/test_singlecavity.jl | 2 +- 4 files changed, 159 insertions(+), 32 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index 75290b1..bd491f9 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -7,7 +7,7 @@ Added functionalities: * [`WaveguideBasis`](@ref) for representing the waveguide Hilbert space and the related functions for generating states in this Hilbert space: [`zerophoton`](@ref), [`onephoton`](@ref), and [`twophoton`](@ref). Also see [`OnePhotonView`](@ref), [`TwoPhotonView`](@ref), and [`plot_twophoton!`](@ref) for viewing the waveguide states and plotting them. Note that [`WaveguideBasis`](@ref) can contain multiple waveguides. * [`WaveguideOperator`](@ref) are specialized operators allowing efficient annihilation and creation operators at each time-bin in the waveguide. They are created by giving a basis to [`WaveguideQED.destroy`](@ref) and [`WaveguideQED.create`](@ref) * Since the interaction between the waveguide time-bin mode $k$ and cavity/emitter is given as: $a^\dagger w_k - a w_k^\dagger$ there are specially optimized functions for doing these operations called [`CavityWaveguideOperator`](@ref) which are created using a fockbasis and a waveguide basis and the functions [`emission`](@ref) and [`absorption`](@ref). -* (OBSOLETE. SEE [Beamsplitter](@ref) INSTEAD). [`Detector`](@ref), [`LazyTensorKet`](@ref), and [`LazySumKet`](@ref), together with [`detect_single_click`](@ref) and [`detect_double_click`](@ref) allow one to do a beamsplitter interference and subsequent detection of photons coming from two waveguides. +* (OBSOLETE. SEE [Beamsplitter](@ref Beamsplitter) INSTEAD). [`Detector`](@ref), [`LazyTensorKet`](@ref), and [`LazySumKet`](@ref), together with [`detect_single_click`](@ref) and [`detect_double_click`](@ref) allow one to do a beamsplitter interference and subsequent detection of photons coming from two waveguides. ![alt text](./animations/firstgif.gif) diff --git a/src/WaveguideInteraction.jl b/src/WaveguideInteraction.jl index 79dce85..d140cf7 100644 --- a/src/WaveguideInteraction.jl +++ b/src/WaveguideInteraction.jl @@ -208,11 +208,26 @@ function waveguide_interaction_mul!(result,a::WaveguideCreate{B1,B2,1,idx1},bop: if !isone(beta) rmul!(result,beta) end - timeindex_a = a.timeindex - timeindex_b = bop.timeindex + + nsteps = a.basis_l.nsteps + timeindex_a = (a.timeindex + a.delay -1) % (nsteps) +1 + timeindex_b = (bop.timeindex + bop.delay -1) % (nsteps) +1 + nsteps = a.basis_l.nsteps - if timeindex_a>0 && timeindex_b>0 - @inbounds result[1+(idx1-1)*nsteps+timeindex_a] += alpha*a.factor*bop.factor*b[1+(idx2-1)*nsteps+timeindex_b] + @inbounds result[1+(idx1-1)*nsteps+timeindex_a] += alpha*a.factor*bop.factor*b[1+(idx2-1)*nsteps+timeindex_b] + result +end + +function waveguide_interaction_mul!(result,a::WaveguideDestroy{B1,B2,1,idx1},bop::WaveguideDestroy{B1,B2,1,idx2},b,alpha,beta) where {B1,B2,idx1,idx2} + if !isone(beta) + rmul!(result,beta) + end + result +end + +function waveguide_interaction_mul!(result,a::WaveguideCreate{B1,B2,1,idx1},bop::WaveguideCreate{B1,B2,1,idx2},b,alpha,beta) where {B1,B2,idx1,idx2} + if !isone(beta) + rmul!(result,beta) end result end @@ -221,32 +236,35 @@ function waveguide_interaction_mul!(result,a::WaveguideCreate{B1,B2,2,idx1},bop: if !isone(beta) rmul!(result,beta) end - timeindex = a.timeindex nsteps = a.basis_l.nsteps + + timeindex_a = (a.timeindex + a.delay -1) % (nsteps) +1 + timeindex_b = (bop.timeindex + bop.delay -1) % (nsteps) +1 + Nw = get_number_of_waveguides(a.basis_l) - @inbounds result[1+(idx1-1)*nsteps+timeindex] += alpha*a.factor*bop.factor*b[1+(idx2-1)*nsteps+timeindex] + @inbounds result[1+(idx1-1)*nsteps+timeindex_a] += alpha*a.factor*bop.factor*b[1+(idx2-1)*nsteps+timeindex_b] - twophotonview_b = TwoPhotonTimestepView(b,timeindex,nsteps,Nw*nsteps+1+(idx2-1)*(nsteps*(nsteps+1))÷2) + twophotonview_b = TwoPhotonTimestepView(b,timeindex_b,nsteps,Nw*nsteps+1+(idx2-1)*(nsteps*(nsteps+1))÷2) i,j = min(idx1,idx2),max(idx1,idx2) index = (i-1)*Nw + j - (i*(i+1))÷2 - twophotonview_r = TwoWaveguideTimestepView(result,timeindex,nsteps,1+Nw*nsteps+Nw*(nsteps*(nsteps+1))÷2+(index-1)*nsteps^2,i==idx1) + twophotonview_r = TwoWaveguideTimestepView(result,timeindex_a,nsteps,1+Nw*nsteps+Nw*(nsteps*(nsteps+1))÷2+(index-1)*nsteps^2,i==idx1) axpy!(alpha*a.factor*bop.factor,twophotonview_b,twophotonview_r) i,j = min(idx1,idx2),max(idx1,idx2) index = (i-1)*Nw + j - (i*(i+1))÷2 - twophotonview_b = TwoWaveguideTimestepView(b,timeindex,nsteps,1+Nw*nsteps+Nw*(nsteps*(nsteps+1))÷2+(index-1)*nsteps^2,i==idx2) - twophotonview_r = TwoPhotonTimestepView(result,timeindex,nsteps,1+Nw*nsteps+(idx1-1)*(nsteps*(nsteps+1))÷2) + twophotonview_b = TwoWaveguideTimestepView(b,timeindex_b,nsteps,1+Nw*nsteps+Nw*(nsteps*(nsteps+1))÷2+(index-1)*nsteps^2,i==idx2) + twophotonview_r = TwoPhotonTimestepView(result,timeindex_a,nsteps,1+Nw*nsteps+(idx1-1)*(nsteps*(nsteps+1))÷2) axpy!(alpha*a.factor*bop.factor,twophotonview_b,twophotonview_r) @simd for k in filter(x -> x != idx1 && x != idx2, 1:Nw) m,l = min(k,idx1),max(k,idx1) index_r = (m-1)*Nw + l - (m*(m+1))÷2 - twophotonview_r = TwoWaveguideTimestepView(result,timeindex,nsteps,1+Nw*nsteps+Nw*(nsteps*(nsteps+1))÷2+(index_r-1)*nsteps^2,m==idx1) + twophotonview_r = TwoWaveguideTimestepView(result,timeindex_a,nsteps,1+Nw*nsteps+Nw*(nsteps*(nsteps+1))÷2+(index_r-1)*nsteps^2,m==idx1) i,j = min(k,idx2),max(k,idx2) index = (i-1)*Nw + j - (i*(i+1))÷2 - twophotonview_b = TwoWaveguideTimestepView(b,timeindex,nsteps,1+Nw*nsteps+Nw*(nsteps*(nsteps+1))÷2+(index-1)*nsteps^2,i==idx2) + twophotonview_b = TwoWaveguideTimestepView(b,timeindex_b,nsteps,1+Nw*nsteps+Nw*(nsteps*(nsteps+1))÷2+(index-1)*nsteps^2,i==idx2) axpy!(alpha*a.factor*bop.factor,twophotonview_b,twophotonview_r) end result @@ -256,21 +274,24 @@ function waveguide_interaction_mul!(result,a::WaveguideCreate{B1,B2,2,idx},bop:: if !isone(beta) rmul!(result,beta) end - timeindex = a.timeindex nsteps = a.basis_l.nsteps + timeindex_a = (a.timeindex + a.delay -1) % (nsteps) +1 + timeindex_b = (bop.timeindex + bop.delay -1) % (nsteps) +1 + + Nw = get_number_of_waveguides(a.basis_l) - @inbounds result[1+(idx-1)*nsteps+timeindex] += alpha*a.factor*bop.factor*b[1+(idx-1)*nsteps+timeindex] + @inbounds result[1+(idx-1)*nsteps+timeindex_a] += alpha*a.factor*bop.factor*b[1+(idx-1)*nsteps+timeindex_b] - twophotonview_b = TwoPhotonTimestepView(b,timeindex,nsteps,Nw*nsteps+1+(idx-1)*(nsteps*(nsteps+1))÷2) - twophotonview_r = TwoPhotonTimestepView(result,timeindex,nsteps,Nw*nsteps+1+(idx-1)*(nsteps*(nsteps+1))÷2) + twophotonview_b = TwoPhotonTimestepView(b,timeindex_b,nsteps,Nw*nsteps+1+(idx-1)*(nsteps*(nsteps+1))÷2) + twophotonview_r = TwoPhotonTimestepView(result,timeindex_a,nsteps,Nw*nsteps+1+(idx-1)*(nsteps*(nsteps+1))÷2) axpy!(alpha*a.factor*bop.factor,twophotonview_b,twophotonview_r) @simd for k in filter(x -> x != idx, 1:Nw) m,l = min(k,idx),max(k,idx) index_r = (m-1)*Nw + l - (m*(m+1))÷2 - twophotonview_r = TwoWaveguideTimestepView(result,timeindex,nsteps,1+Nw*nsteps+Nw*(nsteps*(nsteps+1))÷2+(index_r-1)*nsteps^2,m==idx) - twophotonview_b = TwoWaveguideTimestepView(b,timeindex,nsteps,1+Nw*nsteps+Nw*(nsteps*(nsteps+1))÷2+(index_r-1)*nsteps^2,m==idx) + twophotonview_r = TwoWaveguideTimestepView(result,timeindex_a,nsteps,1+Nw*nsteps+Nw*(nsteps*(nsteps+1))÷2+(index_r-1)*nsteps^2,m==idx) + twophotonview_b = TwoWaveguideTimestepView(b,timeindex_b,nsteps,1+Nw*nsteps+Nw*(nsteps*(nsteps+1))÷2+(index_r-1)*nsteps^2,m==idx) axpy!(alpha*a.factor*bop.factor,twophotonview_b,twophotonview_r) end result @@ -280,21 +301,36 @@ function waveguide_interaction_mul!(result,a::WaveguideDestroy{B1,B2,2,idx},bop: if !isone(beta) rmul!(result,beta) end - timeindex = a.timeindex nsteps = a.basis_l.nsteps + timeindex_a = (a.timeindex + a.delay -1) % (nsteps) +1 + timeindex_b = (bop.timeindex + bop.delay -1) % (nsteps) +1 + + t2 = timeindex_a >= timeindex_b ? timeindex_a : timeindex_b + t1 = timeindex_a >= timeindex_b ? timeindex_b : timeindex_a + + Nw = get_number_of_waveguides(a.basis_l) - result[1] += sqrt(2)*alpha*a.factor*bop.factor*b[1+Nw*nsteps+(idx-1)*(nsteps*(nsteps+1))÷2+twophoton_index(timeindex,nsteps,timeindex)] + factor = timeindex_a==timeindex_b ? sqrt(2) : 1 + + result[1] += factor*alpha*a.factor*bop.factor*b[1+Nw*nsteps+(idx-1)*(nsteps*(nsteps+1))÷2+twophoton_index(t1,nsteps,t2)] result end function waveguide_interaction_mul!(result,a::WaveguideCreate{B1,B2,2,idx},bop::WaveguideCreate{B1,B2,2,idx},b,alpha,beta) where {B1,B2,idx} if !isone(beta) rmul!(result,beta) end - timeindex = a.timeindex nsteps = a.basis_l.nsteps + timeindex_a = (a.timeindex + a.delay -1) % (nsteps) +1 + timeindex_b = (bop.timeindex + bop.delay -1) % (nsteps) +1 + + t2 = timeindex_a >= timeindex_b ? timeindex_a : timeindex_b + t1 = timeindex_a >= timeindex_b ? timeindex_b : timeindex_a + Nw = get_number_of_waveguides(a.basis_l) - result[1+Nw*nsteps+(idx-1)*(nsteps*(nsteps+1))÷2+twophoton_index(timeindex,nsteps,timeindex)] += sqrt(2)*alpha*a.factor*bop.factor*b[1] + factor = timeindex_a==timeindex_b ? sqrt(2) : 1 + + result[1+Nw*nsteps+(idx-1)*(nsteps*(nsteps+1))÷2+twophoton_index(t1,nsteps,t2)] += factor*alpha*a.factor*bop.factor*b[1] result end @@ -302,37 +338,84 @@ function waveguide_interaction_mul!(result,a::WaveguideDestroy{B1,B2,2,idx1},bop if !isone(beta) rmul!(result,beta) end - timeindex = a.timeindex nsteps = a.basis_l.nsteps + timeindex_a = (a.timeindex + a.delay -1) % (nsteps) +1 + timeindex_b = (bop.timeindex + bop.delay -1) % (nsteps) +1 + + t1 = idx1 >= idx2 ? timeindex_a : timeindex_b + t2 = idx1 >= idx2 ? timeindex_b : timeindex_a + + Nw = get_number_of_waveguides(a.basis_l) m,l = min(idx1,idx2),max(idx1,idx2) index_r = (m-1)*Nw + l - (m*(m+1))÷2 - result[1] += alpha*a.factor*bop.factor*b[1+Nw*nsteps+Nw*(nsteps*(nsteps+1))÷2+(index_r-1)*nsteps^2+(timeindex-1)*nsteps + timeindex] + result[1] += alpha*a.factor*bop.factor*b[1+Nw*nsteps+Nw*(nsteps*(nsteps+1))÷2+(index_r-1)*nsteps^2+(t1-1)*nsteps + t2] result end function waveguide_interaction_mul!(result,a::WaveguideCreate{B1,B2,2,idx1},bop::WaveguideCreate{B1,B2,2,idx2},b,alpha,beta) where {B1,B2,idx1,idx2} if !isone(beta) rmul!(result,beta) end - timeindex = a.timeindex nsteps = a.basis_l.nsteps + timeindex_a = (a.timeindex + a.delay -1) % (nsteps) +1 + timeindex_b = (bop.timeindex + bop.delay -1) % (nsteps) +1 + t1 = idx1 >= idx2 ? timeindex_a : timeindex_b + t2 = idx1 >= idx2 ? timeindex_b : timeindex_a + #a >= timeindex_b ? timeindex_a : timeindex_b + #t2 = timeindex_b + #>= timeindex_b ? timeindex_b : timeindex_a + Nw = get_number_of_waveguides(a.basis_l) m,l = min(idx1,idx2),max(idx1,idx2) index_r = (m-1)*Nw + l - (m*(m+1))÷2 - result[1+Nw*nsteps+Nw*(nsteps*(nsteps+1))÷2+(index_r-1)*nsteps^2+(timeindex-1)*nsteps + timeindex] += alpha*a.factor*bop.factor*b[1] + + result[1+Nw*nsteps+Nw*(nsteps*(nsteps+1))÷2+(index_r-1)*nsteps^2+(t1-1)*nsteps + t2] += alpha*a.factor*bop.factor*b[1] end -function waveguide_interaction_mul!(result,a::WaveguideDestroy{B1,B2,Np,idx1},bop::WaveguideCreate{B1,B2,1,idx2},b,alpha,beta) where {B1,B2,Np,idx1,idx2} +function waveguide_interaction_mul!(result,a::WaveguideDestroy{B1,B2,1,idx1},bop::WaveguideCreate{B1,B2,1,idx2},b,alpha,beta) where {B1,B2,idx1,idx2} + if !isone(beta) + rmul!(result,beta) + end + nsteps = a.basis_l.nsteps + timeindex_a = (a.timeindex + a.delay -1) % (nsteps) +1 + timeindex_b = (bop.timeindex + bop.delay -1) % (nsteps) +1 + + if timeindex_a == timeindex_b && idx1 == idx2 + @inbounds result[1] += alpha*a.factor*bop.factor*b[1] + end + + #result[1] = alpha*a.factor*bop.factor*b[1] + #@inbounds result[1+(idx1-1)*nsteps+timeindex_a] += alpha*a.factor*bop.factor*b[1+(idx2-1)*nsteps+timeindex_b] + result +end + + +function waveguide_interaction_mul!(result,a::WaveguideDestroy{B1,B2,2,idx1},bop::WaveguideCreate{B1,B2,2,idx2},b,alpha,beta) where {B1,B2,idx1,idx2} if !isone(beta) rmul!(result,beta) end - timeindex = a.timeindex nsteps = a.basis_l.nsteps + timeindex_a = (a.timeindex + a.delay -1) % (nsteps) +1 + timeindex_b = (bop.timeindex + bop.delay -1) % (nsteps) +1 + + Nw = get_number_of_waveguides(a.basis_l) + + if timeindex_a == timeindex_b && idx1==idx2 + @inbounds result[1] += alpha*a.factor*bop.factor*b[1] + for j in 1:Nw + @inbounds result[1+(j-1)*nsteps+1:1+(j-1)*nsteps+timeindex_a - 1] .+= b[1+(j-1)*nsteps+1:1+(j-1)*nsteps+timeindex_a - 1] + @inbounds result[1+(j-1)*nsteps+timeindex_a + 1:1+(j)*nsteps] .+= b[1+(j-1)*nsteps+timeindex_a + 1:1+(j)*nsteps] + factor = j==idx1 ? 2 : 1 + @inbounds result[1+(j-1)*nsteps+timeindex_a] += factor*b[1+(j-1)*nsteps+timeindex_a] + end + else + @inbounds result[1+(idx2-1)*nsteps+timeindex_b] += alpha*a.factor*bop.factor*b[1+(idx1-1)*nsteps+timeindex_a] + end + + #result[1] = alpha*a.factor*bop.factor*b[1] - @inbounds result[1] = alpha*a.factor*bop.factor*b[1] - @inbounds result[1+(idx1-1)*nsteps+timeindex] += alpha*a.factor*bop.factor*b[1+(idx2-1)*nsteps+timeindex] result end diff --git a/test/test_operators.jl b/test/test_operators.jl index b5af45e..90ae07b 100644 --- a/test/test_operators.jl +++ b/test/test_operators.jl @@ -312,3 +312,47 @@ end @test test_interaction(right_1photon[1],3,c[1],c[2],[3,1,2],right_1photon[2][1]) end end + +@testset "Waveguide Delay Interaction" begin + times = 0:1:10 + dt = times[2] - times[1] + + for i in 1:2 + bw = WaveguideBasis(i,2,times) + + tau = 5 + + w1 = destroy(bw,1) + wd1 = create(bw,1); + w1_delay = destroy(bw,1;delay=tau/dt) + wd1_delay = create(bw,1;delay=tau/dt); + + w2 = destroy(bw,2) + wd2 = create(bw,2); + w2_delay = destroy(bw,2;delay=tau/dt) + wd2_delay = create(bw,2;delay=tau/dt); + + ξfun(t1) = 1 + input = Ket(bw) + input.data .= 1 + normalize!(input) + temp1 = copy(input) + temp2 = copy(input) + + + operators = [w1, wd1, w1_delay, wd1_delay, w2, wd2,w2_delay, wd2_delay] + + for o1 in operators + for o2 in operators + c1 = o1*o2 + c2 = LazyProduct([o1,o2]) + + QuantumOptics.mul!(temp1,c1,input) + QuantumOptics.mul!(temp2,c2,input) + + @test isapprox(temp1.data,temp2.data) + end + end + end + +end diff --git a/test/test_singlecavity.jl b/test/test_singlecavity.jl index 68b9414..d1e04a9 100644 --- a/test/test_singlecavity.jl +++ b/test/test_singlecavity.jl @@ -35,7 +35,7 @@ include("helper_functions.jl") #REFERENCE SOLUTION sol1 = solve_differentialeq(param,ξfun) - ref_sol = ξfun.(sol1.t,param.σ,param.t0)-sqrt(param.γ)*sol1 + ref_sol = ξfun.(sol1.t,param.σ,param.t0) .- sqrt(param.γ)*sol1 ψ_single = OnePhotonView(ψ)/sqrt(dt) @test isapprox(ψ_single,ref_sol,rtol=0.05)