update 3D PB and OU random field

Project.toml

@@ -1,7 +1,7 @@
 name = "SpinShuttling"
 uuid = "eff49e45-ae6d-44ae-9170-2124a87971c9"
 authors = ["Yuning Zhang"]
-version = "0.2.4"
+version = "0.2.5"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/SpinShuttling.jl

@@ -131,7 +131,8 @@ function TwoSpinModel(Ψ::Vector{<:Complex}, T::Real, N::Int,
 
     X = [x₁, x₂]
     t = range(0, T, N)
-    P = vcat(collect(zip(t, x₁.(t))), collect(zip(t, x₂.(t))))
+    f(x::Function, t::Real) = (t, x(t)...)
+    P = vcat(f.(x₁, t), f.(x₂, t))
     R = RandomFunction(P, B)
     model = ShuttlingModel(2, Ψ, T, N, B, X, R)
     return model
@@ -177,7 +178,6 @@ The total shuttling time is `T` and the length of the path is `L` in `μm`.
 """
 function TwoSpinParallelModel(T::Real, D::Real, L::Real, N::Int,
     B::GaussianRandomField)
-    @assert length(B.θ) >= 3
     x₁(t::Real)::Tuple{Real,Real} = (L / T * t, 0)
     x₂(t::Real)::Tuple{Real,Real} = (L / T * t, D)
     Ψ = 1 / √2 .* [0.0, 1+0im, -1+0im, 0.0]
@@ -373,7 +373,7 @@ function W(T0::Real, T1::Real, L::Real, B::OrnsteinUhlenbeckField; path=:sequenc
     if path == :sequenced
         return exp(-σ^2 / (4 * κₜ * κₓ) / κₜ^2 * (F1(β, γ, τ) - F2(β, γ, τ)))
     elseif path == :parallel
-        missing("Parallel path not implemented yet.")
+        return exp(-σ^2 / (8 *κₜ*κₓ*κₓ) / κₜ^2 *(1-exp(-κₓ*T1)) * P1(κₜ*T0, κₓ*L))
     else
         error("Path not recognized. Use :sequenced or :parallel for two-spin EPR pair shuttling model.")
     end

test/runtests.jl

@@ -4,8 +4,8 @@ using SpinShuttling
 using UnicodePlots
 using ProgressMeter
 
-# include("testintegration.jl")
-# include("testfidelity.jl")
-# include("teststochastics.jl")
-# include("testspectrum.jl")
+include("testintegration.jl")
+include("testfidelity.jl")
+include("teststochastics.jl")
+include("testspectrum.jl")
 include("testdephasingfactor.jl")

test/testfidelity.jl

@@ -54,23 +54,6 @@ end
     println("TH:", f3)
 end
 
-##
-@testset begin "test two spin sequenced shuttling fidelity"
-    L=10; σ =sqrt(2)/20; M=5000; N=501; T1=200; T0=25*0.05; κₜ=1/20; κₓ=1/0.1;
-    B=OrnsteinUhlenbeckField(0,[κₜ,κₓ],σ)
-    model=TwoSpinSequentialModel(T0, T1, L, N, B)
-    if visualize
-        display(heatmap(collect(model.R.Σ), title="cross covariance matrix, two spin EPR"))
-    end
-    f1=statefidelity(model)
-    f2, f2_err=sampling(model, statefidelity, M)
-    f3=1/2*(1+W(T0, T1, L,B))
-    @test isapprox(f1, f3,rtol=3e-2)
-    @test isapprox(f2, f3, rtol=3e-2) 
-    println("NI:", f1)
-    println("MC:", f2)
-    println("TH:", f3)
-end
 
 #
 @testset "1/f noise chacacteristics" begin
@@ -110,4 +93,64 @@ end
             lineplot!(fig, t, f_th, name="theoretical fidelity")
         display(fig)
     end
-end
+end
+
+##
+@testset begin "test two spin sequenced shuttling fidelity"
+    L=10; σ =sqrt(2)/20; M=5000; N=501; T1=200; T0=25*0.05; κₜ=1/20; κₓ=1/0.1;
+    B=OrnsteinUhlenbeckField(0,[κₜ,κₓ],σ)
+    model=TwoSpinSequentialModel(T0, T1, L, N, B)
+    if visualize
+        display(heatmap(collect(model.R.Σ), title="cross covariance matrix, two spin EPR"))
+    end
+    f1=statefidelity(model)
+    f2, f2_err=sampling(model, statefidelity, M)
+    f3=1/2*(1+W(T0, T1, L,B))
+    @test isapprox(f1, f3,rtol=3e-2)
+    @test isapprox(f2, f3, rtol=3e-2) 
+    println("NI:", f1)
+    println("MC:", f2)
+    println("TH:", f3)
+end
+
+##
+@testset begin "test two spin parallel shuttling fidelity"
+    L=10; σ =sqrt(2)/20; M=5000; N=501; T=200; κₜ=1/20; κₓ=1/0.1;
+    D=0.3;
+    B=OrnsteinUhlenbeckField(0,[κₜ,κₓ,κₓ],σ)
+    model=TwoSpinParallelModel(T, D, L, N, B)
+    if visualize
+        display(heatmap(collect(model.R.Σ), title="cross covariance matrix, two spin EPR"))
+    end
+    f1=statefidelity(model)
+    f2, f2_err=sampling(model, statefidelity, M)
+    w=exp(-σ^2 / (8 *κₜ*κₓ*κₓ) / κₜ^2 *(1-exp(-κₓ*D)) * SpinShuttling.P1(κₜ*T, κₓ*L))
+    f3=1/2*(1+w)
+    @test isapprox(f1, f3,rtol=3e-2)
+    @test isapprox(f2, f3, rtol=3e-2) 
+    println("NI:", f1)
+    println("MC:", f2)
+    println("TH:", f3)
+end
+
+@testset "test two spin parallel shuttling fidelity with 1/f noise" begin
+    σ = sqrt(2)/20; M = 5000; N=501; L=10; γ=(1e-9,1e3); # MHz
+    # 0.01 ~ 100 μs
+    # v = 0.1 ~ 1000 m/s
+    v=1; T=L/v; κₓ=10;
+    # T=10 
+    # 1/T=0.1 N/T=20 
+    D=0.3;
+    B=PinkBrownianField(0,[κₓ,κₓ],σ, γ)
+    model=TwoSpinParallelModel(T, D, L, N, B)
+
+    f1=statefidelity(model)
+    f2, f2_err=sampling(model, statefidelity, M)
+    w=exp(-B.σ^2 * T^2 * 2(1-exp(-κₓ*D)) * SpinShuttling.F3(γ.*T, κₓ*L))
+    f3=1/2*(1+w)
+    @test isapprox(f1, f3,rtol=3e-2)
+    @test isapprox(f2, f3, rtol=3e-2)
+    println("NI:", f1)
+    println("MC:", f2)
+    println("TH:", f3)
+end