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

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

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,59 +301,121 @@ 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
 
 function waveguide_interaction_mul!(result,a::WaveguideDestroy{B1,B2,2,idx1},bop::WaveguideDestroy{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
+
+
  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

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

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)