diff --git a/src/analyticfunctions.jl b/src/analyticfunctions.jl index a047d32..a09e8d3 100644 --- a/src/analyticfunctions.jl +++ b/src/analyticfunctions.jl @@ -46,7 +46,7 @@ where ``s`` is the order of the (blue) sideband that we are driving and ``L_{n}^ associated Laguerre polynomial. [ref](https://doi.org/10.1103/RevModPhys.75.281) """ -function rabiflop(tspan, Ω::Real, η::Real, n̄::Real; s::Int = 0) +function rabiflop(tspan, Ω::Real, η::Real, n̄::Real; s::Int=0) n̄ *= 1.0 Ω *= 2π p = Vector{Float64}(undef, 0) @@ -59,7 +59,7 @@ function rabiflop(tspan, Ω::Real, η::Real, n̄::Real; s::Int = 0) Ω / 2 * exp(-η^2 / 2) * η^s * - sqrt(1 / prod((n + 1):(n + s))) * + sqrt(1 / prod((n+1):(n+s))) * _alaguerre(η^2, n, s) * t )^2 @@ -72,7 +72,7 @@ end function _laguerre(x, n) L = 1.0, -x + 1 if n < 2 - return L[n + 1] + return L[n+1] end for i in 2:n L = L[2], ((2i - 1 - x) * L[2] - (i - 1) * L[1]) / i diff --git a/src/chambers.jl b/src/chambers.jl index e7979e8..07a723a 100644 --- a/src/chambers.jl +++ b/src/chambers.jl @@ -89,11 +89,11 @@ mutable struct Chamber _cnst_δB::Bool function Chamber(; iontrap::LinearChain, - B = 0, - Bhat = ẑ, - ∇B = 0, - δB::TδB = 0, - lasers = Laser[] + B=0, + Bhat=ẑ, + ∇B=0, + δB::TδB=0, + lasers=Laser[] ) where {TδB} warn = nothing for i in 1:length(lasers) @@ -102,7 +102,7 @@ mutable struct Chamber push!(lasers[i].pointing, (n, 1.0)) end end - for j in (i + 1):length(lasers) + for j in (i+1):length(lasers) if lasers[j] ≡ lasers[i] lasers[j] = copy(lasers[i]) if isnothing(warn) @@ -112,7 +112,7 @@ mutable struct Chamber end end end - @assert isapprox(norm(Bhat), 1, rtol = 1e-6) "!(|$Bhat| = 1)" + @assert isapprox(norm(Bhat), 1, rtol=1e-6) "!(|$Bhat| = 1)" for (li, l) in enumerate(lasers), p in l.pointing @assert p[1] <= length(iontrap.ions) ( """lasers[$li] points at iontrap.ions[$(p[1])], but there are only @@ -200,7 +200,7 @@ Sets `chamber.Bhat` to `Bhat` """ function bfield_unitvector!(chamber::Chamber, Bhat::NamedTuple{(:x, :y, :z)}) rtol = 1e-6 - @assert isapprox(norm(Bhat), 1, rtol = rtol) "!(|̂B| = 1)" + @assert isapprox(norm(Bhat), 1, rtol=rtol) "!(|̂B| = 1)" chamber.Bhat = Bhat end @@ -247,14 +247,14 @@ end basis(chamber::Chamber) Returns the composite basis describing the Hilbert space for `chamber`. This is the same as basis(iontrap(chain)). -""" +""" function basis(T::Chamber) - return tensor( - T.iontrap.ions..., - T.iontrap.selectedmodes.x..., - T.iontrap.selectedmodes.y..., - T.iontrap.selectedmodes.z..., - ) + return tensor( + T.iontrap.ions..., + T.iontrap.selectedmodes.x..., + T.iontrap.selectedmodes.y..., + T.iontrap.selectedmodes.z..., + ) end @@ -347,7 +347,9 @@ function intensity_from_pitime( s_indx = findall(x -> x[1] == ionnumber(ion), p) @assert length(s_indx) > 0 "This laser doesn't shine on this ion" s = p[s_indx[1]][2] - Ω = s * matrixelement(ion, transition, 1.0, polarization(laser), wavevector(laser), Bhat) + Ω = + s * + matrixelement(ion, transition, 1.0, polarization(laser), wavevector(laser), Bhat) if Ω < 3e-14 # After change from Efield to intensity: This inequality changed so that it serves the same function but now for intensity = 1.0 rather than efield = 1.0 (specificially, increased by a factor of √(2/(cϵ₀)) ∼ 30 in SI units) # even when coupling strength is zero, numerical error causes it to be finite # (on order 1e-14), this is a band-aid to prevent users from unknowingly setting @@ -373,7 +375,13 @@ function intensity_from_pitime( transition::Tuple, chamber::Chamber ) - return intensity_from_pitime(laser, pi_time, ion, transition, bfield_unitvector(chamber)) + return intensity_from_pitime( + laser, + pi_time, + ion, + transition, + bfield_unitvector(chamber) + ) end function intensity_from_pitime( laser_index::Int, @@ -387,7 +395,8 @@ function intensity_from_pitime( pi_time, ions(chamber)[ion_index], transition, - bfield_unitvector(chamber)) + bfield_unitvector(chamber) + ) end @@ -453,7 +462,7 @@ function intensity_from_rabifrequency( transition::Tuple, Bhat::NamedTuple{(:x, :y, :z)} ) - return intensity_from_pitime(laser, 1/(2*rabi_frequency), ion, transition, Bhat) + return intensity_from_pitime(laser, 1 / (2 * rabi_frequency), ion, transition, Bhat) end """ @@ -473,7 +482,13 @@ function intensity_from_rabifrequency( transition::Tuple, chamber::Chamber ) - return intensity_from_rabifrequency(laser, rabi_frequency, ion, transition, bfield_unitvector(chamber)) + return intensity_from_rabifrequency( + laser, + rabi_frequency, + ion, + transition, + bfield_unitvector(chamber) + ) end function intensity_from_rabifrequency( laser_index::Int, @@ -487,7 +502,8 @@ function intensity_from_rabifrequency( rabi_frequency, ions(chamber)[ion_index], transition, - bfield_unitvector(chamber)) + bfield_unitvector(chamber) + ) end @@ -519,7 +535,8 @@ function intensity_from_rabifrequency!( transition::Tuple, chamber::Chamber ) - I::Float64 = intensity_from_rabifrequency(laser, rabi_frequency, ion, transition, chamber) + I::Float64 = + intensity_from_rabifrequency(laser, rabi_frequency, ion, transition, chamber) intensity!(laser, I) return I end @@ -530,7 +547,13 @@ function intensity_from_rabifrequency!( transition::Tuple, chamber::Chamber ) - I::Float64 = intensity_from_rabifrequency(laser_index, rabi_frequency, ion_index, transition, chamber) + I::Float64 = intensity_from_rabifrequency( + laser_index, + rabi_frequency, + ion_index, + transition, + chamber + ) laser = lasers(chamber)[laser_index] intensity!(laser, I) return I @@ -555,43 +578,55 @@ end Returns The frequency of the transition `transition` including the Zeeman shift experienced by `ion` given its position in `T`. `ion` may be either an Ion or an Int indicating the desired ion's index within `chamber`. """ -transitionfrequency(ion::Ion, transition::Tuple, chamber::Chamber; ignore_manualshift = false) = - transitionfrequency( - ion, - transition; - B = bfield(chamber, ion), - ignore_manualshift = ignore_manualshift - ) -transitionfrequency(ion_index::Int, transition::Tuple, chamber::Chamber; ignore_manualshift = false) = - transitionfrequency( - ions(chamber)[ion_index], - transition, - chamber, - ignore_manualshift = ignore_manualshift - ) +transitionfrequency( + ion::Ion, + transition::Tuple, + chamber::Chamber; + ignore_manualshift=false +) = transitionfrequency( + ion, + transition; + B=bfield(chamber, ion), + ignore_manualshift=ignore_manualshift +) +transitionfrequency( + ion_index::Int, + transition::Tuple, + chamber::Chamber; + ignore_manualshift=false +) = transitionfrequency( + ions(chamber)[ion_index], + transition, + chamber, + ignore_manualshift=ignore_manualshift +) """ transitionwavelength(ion, transition::Tuple, chamber::Chamber; ignore_manualshift=false) Returns The wavelength of the transition `transition` including the Zeeman shift experienced by `ion` given its position in `T`. `ion` may be either an Ion or an Int indicating the desired ion's index within `chamber`. """ -transitionwavelength(ion::Ion, transition::Tuple, chamber::Chamber; ignore_manualshift = false) = - transitionwavelength( - ion, - transition; - B = bfield(chamber, ion), - ignore_manualshift = ignore_manualshift - ) +transitionwavelength( + ion::Ion, + transition::Tuple, + chamber::Chamber; + ignore_manualshift=false +) = transitionwavelength( + ion, + transition; + B=bfield(chamber, ion), + ignore_manualshift=ignore_manualshift +) transitionwavelength( ion_index::Int, transition::Tuple, chamber::Chamber; - ignore_manualshift = false + ignore_manualshift=false ) = transitionwavelength( ions(chamber)[ion_index], transition, chamber; - ignore_manualshift = ignore_manualshift + ignore_manualshift=ignore_manualshift ) """ @@ -611,12 +646,22 @@ Sets the wavelength of `laser` to the transition wavelength of `transition` in t at the magnetic field seen by `ion` in `chamber`. `ion` may be either an Ion or an Int indicating the desired ion's index within `chamber`. """ -function wavelength_from_transition!(laser::Laser, ion::Ion, transition::Tuple, chamber::Chamber) +function wavelength_from_transition!( + laser::Laser, + ion::Ion, + transition::Tuple, + chamber::Chamber +) λ = transitionwavelength(ion, transition, chamber) wavelength!(laser, λ) return λ end -function wavelength_from_transition!(laser::Laser, ion_index::Int, transition::Tuple, chamber::Chamber) +function wavelength_from_transition!( + laser::Laser, + ion_index::Int, + transition::Tuple, + chamber::Chamber +) ion = ions(chamber)[ion_index] λ = transitionwavelength(ion, transition, chamber) wavelength!(laser, λ) @@ -633,7 +678,8 @@ Calls `matrixelement(ion, transition, I, ϵhat, khat, Bhat)` with `I`, `ϵhat`, `ion` and `laser` must either both be indices or both their respective Structs. """ matrixelement(ion::Ion, transition::Tuple, laser::Laser, chamber::Chamber, time::Real) = - matrixelement(ion, + matrixelement( + ion, transition, intensity(laser)(time), polarization(laser), @@ -706,7 +752,7 @@ The Lamb-Dicke parameter: ``|k|cos(\\theta)\\sqrt{\\frac{\\hbar}{2m\\nu}}`` for a given vibrational mode, ion and laser. """ -function lambdicke(V::VibrationalMode, I::Ion, L::Laser; scaled = false) +function lambdicke(V::VibrationalMode, I::Ion, L::Laser; scaled=false) @fastmath begin k = 2π / L.λ scaled ? ν = 1 : ν = V.ν diff --git a/src/constants.jl b/src/constants.jl index cf7235b..a135e30 100644 --- a/src/constants.jl +++ b/src/constants.jl @@ -57,7 +57,7 @@ const c_rank2 = cat( 1 / 3 * [[-1, 0, 0] [0, -1, 0] [0, 0, 2]], #q=0 1 / sqrt(6) * [[0, 0, -1] [0, 0, im] [-1, im, 0]], #q=1 1 / sqrt(6) * [[1, -im, 0] [-im, -1, 0] [0, 0, 0]]; - dims = 3 + dims=3 ) #q=2 Base.print(pc::PhysicalConstant) = print("$(pc.x) [$(pc.units)]") @@ -83,9 +83,9 @@ end # module export x̂, ŷ, ẑ, ndot -const x̂ = (x = 1, y = 0, z = 0) -const ŷ = (x = 0, y = 1, z = 0) -const ẑ = (x = 0, y = 0, z = 1) +const x̂ = (x=1, y=0, z=0) +const ŷ = (x=0, y=1, z=0) +const ẑ = (x=0, y=0, z=1) function _print_axis(a::NamedTuple{(:x, :y, :z)}) if a == x̂ @@ -102,11 +102,11 @@ end ndot(a::NamedTuple{(:x, :y, :z)}, b::NamedTuple{(:x, :y, :z)}) = a.x * b.x + a.y * b.y + a.z * b.z function Base.:+(a::NamedTuple{(:x, :y, :z)}, b::NamedTuple{(:x, :y, :z)}) - return (x = a.x + b.x, y = a.y + b.y, z = a.z + b.z) + return (x=a.x + b.x, y=a.y + b.y, z=a.z + b.z) end function Base.:-(a::NamedTuple{(:x, :y, :z)}, b::NamedTuple{(:x, :y, :z)}) - return (x = a.x - b.x, y = a.y - b.y, z = a.z - b.z) + return (x=a.x - b.x, y=a.y - b.y, z=a.z - b.z) end -Base.:/(a::NamedTuple{(:x, :y, :z)}, b::Number) = (x = a.x / b, y = a.y / b, z = a.z / b) -Base.:*(a::NamedTuple{(:x, :y, :z)}, b::Number) = (x = a.x * b, y = a.y * b, z = a.z * b) -Base.:*(b::Number, a::NamedTuple{(:x, :y, :z)}) = (x = a.x * b, y = a.y * b, z = a.z * b) +Base.:/(a::NamedTuple{(:x, :y, :z)}, b::Number) = (x=a.x / b, y=a.y / b, z=a.z / b) +Base.:*(a::NamedTuple{(:x, :y, :z)}, b::Number) = (x=a.x * b, y=a.y * b, z=a.z * b) +Base.:*(b::Number, a::NamedTuple{(:x, :y, :z)}) = (x=a.x * b, y=a.y * b, z=a.z * b) diff --git a/src/hamiltonians.jl b/src/hamiltonians.jl index 31670dc..c7b2bb0 100644 --- a/src/hamiltonians.jl +++ b/src/hamiltonians.jl @@ -51,11 +51,11 @@ Constructs the Hamiltonian for `chamber` as a function of time. Return type is a """ function hamiltonian( chamber::Chamber; - timescale::Real = 1, - lamb_dicke_order::Union{Vector{Int}, Int} = 1, - rwa_cutoff::Real = Inf, - displacement::String = "truncated", - time_dependent_eta::Bool = false + timescale::Real=1, + lamb_dicke_order::Union{Vector{Int}, Int}=1, + rwa_cutoff::Real=Inf, + displacement::String="truncated", + time_dependent_eta::Bool=false ) return hamiltonian( chamber, @@ -104,7 +104,7 @@ function hamiltonian( S.data[i2, i1] = conj(bt_i) if length(cindxs[i]) != 0 flag = cindxs[i][1][1] - i3, i4 = cindxs[i][j + 1] + i3, i4 = cindxs[i][j+1] if flag == -1 S.data[i3, i4] = -conj_bt_i S.data[i4, i3] = -conj(conj_bt_i) @@ -194,7 +194,8 @@ function _setup_base_hamiltonian( Q = prod([shape(ion)[1] for ion in all_ions]) ion_arrays = [spdiagm(0 => [true for _ in 1:shape(ion)[1]]) for ion in all_ions] - ηm, Δm, Ωm = _ηmatrix(chamber), _Δmatrix(chamber, timescale), _Ωmatrix(chamber, timescale) + ηm, Δm, Ωm = + _ηmatrix(chamber), _Δmatrix(chamber, timescale), _Ωmatrix(chamber, timescale) lamb_dicke_order = _check_lamb_dicke_order(lamb_dicke_order, L) ld_array, rows, vals = _ld_array(mode_dims, lamb_dicke_order, νlist, timescale) if displacement == "truncated" && time_dependent_eta @@ -219,7 +220,7 @@ function _setup_base_hamiltonian( if length(all_ions) == 1 K = C else - K = kron(ion_arrays[1:(n - 1)]..., C, ion_arrays[(n + 1):length(all_ions)]...) + K = kron(ion_arrays[1:(n-1)]..., C, ion_arrays[(n+1):length(all_ions)]...) end ion_rows, ion_cols, ion_vals = findnz(K) ion_idxs = sortperm(real.(ion_vals)) @@ -245,7 +246,7 @@ function _setup_base_hamiltonian( Δ, Ω = Δm[rn, m][ti], Ωm[rn, m][ti] Δ_2π = Δ / 2π typeof(Ω) <: Number && continue # e.g. the laser doesn't shine on this ion - locs = view(ion_idxs, ((ti - 1) * ion_reps + 1):(ti * ion_reps)) + locs = view(ion_idxs, ((ti-1)*ion_reps+1):(ti*ion_reps)) for j in 1:prod(mode_dims) for i in nzrange(ld_array, j) ri = rows[i] @@ -426,7 +427,7 @@ function _setup_δν_hamiltonian(chamber, timescale) mode_op = number(mode) A = embed(basis(chamber), [N + l], [mode_op]).data mode_dim = mode.shape[1] - for i in 1:(mode_dim - 1) + for i in 1:(mode_dim-1) indices = [x[1] for x in getfield.(findall(x -> x .== complex(i, 0), A), :I)] push!(δν_indices_l, indices) end @@ -501,11 +502,11 @@ function _ηmatrix(T) for n in 1:N, m in 1:M, l in 1:L δν = frequency_fluctuation(vms[l]) ν = frequency(vms[l]) - eta = lambdicke(vms[l], all_ions[n], all_lasers[m], scaled = true) + eta = lambdicke(vms[l], all_ions[n], all_lasers[m], scaled=true) if eta == 0 - ηnml[n, m, L - l + 1] = 0 + ηnml[n, m, L-l+1] = 0 else - ηnml[n, m, L - l + 1] = + ηnml[n, m, L-l+1] = FunctionWrapper{Float64, Tuple{Float64}}(t -> eta / √(ν + δν(t))) end end @@ -527,8 +528,13 @@ function _Δmatrix(chamber, timescale) Btot = bfield(chamber, all_ions[n]) v = Vector{Float64}(undef, 0) for transition in subleveltransitions(all_ions[n]) - ωa = transitionfrequency(all_ions[n], transition, B = Btot) - push!(v, 2π * timescale * ((c / wavelength(all_lasers[m])) + detuning(all_lasers[m]) - ωa)) + ωa = transitionfrequency(all_ions[n], transition, B=Btot) + push!( + v, + 2π * + timescale * + ((c / wavelength(all_lasers[m])) + detuning(all_lasers[m]) - ωa) + ) end Δnmkj[n, m] = v end @@ -560,7 +566,14 @@ function _Ωmatrix(chamber, timescale) 2π * timescale * s * - matrixelement(all_ions[n], t, 1.0, polarization(all_lasers[m]), wavevector(all_lasers[m]), bfield_unitvector(chamber)) / 2.0 + matrixelement( + all_ions[n], + t, + 1.0, + polarization(all_lasers[m]), + wavevector(all_lasers[m]), + bfield_unitvector(chamber) + ) / 2.0 if Ω0 == 0 push!(v, 0) else @@ -607,11 +620,11 @@ end function _get_kron_indxs(indxs::Vector{Tuple{Int64, Int64}}, dims::Vector{Int64}) L = length(indxs) rowcol = Int64[0, 0] - @simd for i in 0:(L - 1) + @simd for i in 0:(L-1) if i == 0 - @inbounds rowcol .+= indxs[L - i] + @inbounds rowcol .+= indxs[L-i] else - @inbounds rowcol .+= (indxs[L - i] .- 1) .* prod(view(dims, 1:i)) + @inbounds rowcol .+= (indxs[L-i] .- 1) .* prod(view(dims, 1:i)) end end return rowcol @@ -653,7 +666,7 @@ function _Dnm_cnst_eta(ξ::Number, n::Int, m::Int) n -= 1 m -= 1 s = 1.0 - for i in (m + 1):n + for i in (m+1):n s *= i end ret = sqrt(1 / s) * ξ^(n - m) * exp(-abs2(ξ) / 2.0) * _alaguerre(abs2(ξ), m, n - m) diff --git a/src/ions.jl b/src/ions.jl index ee4eb52..0ef90ca 100644 --- a/src/ions.jl +++ b/src/ions.jl @@ -166,10 +166,7 @@ function sublevelalias!( sublevelalias!(I, sublevel, alias) end end -function sublevelalias!( - I::Ion, - aliasdict::Dict{String, Tuple{String, R}} where {R <: Real} -) +function sublevelalias!(I::Ion, aliasdict::Dict{String, Tuple{String, R}} where {R <: Real}) for (alias, sublevel) in aliasdict sublevelalias!(I, sublevel, alias) end @@ -335,7 +332,7 @@ Landé g-factor of fine structure energy level * `j`: electron total angular momentum quantum number * `s`: electronic spin angular momentum quantum number (defaults to 1/2) """ -landegj(l::Real, j::Real, s::Real = 1 // 2) = +landegj(l::Real, j::Real, s::Real=1 // 2) = 3 // 2 + (s * (s + 1) - l * (l + 1)) / (2j * (j + 1)) """ @@ -349,7 +346,7 @@ Landé g-factor of hyperfine energy level * `i`: nuclear spin angular momentum quantum number * `s`: electronic spin angular momentum quantum number (defaults to 1/2) """ -landegf(l::Real, j::Real, f::Real, i::Real, s::Real = 1 // 2) = +landegf(l::Real, j::Real, f::Real, i::Real, s::Real=1 // 2) = landegj(l, j, s) / 2 * (1 + ((j * (j + 1) - i * (i + 1)) / (f * (f + 1)))) landegf(qnums::NamedTuple) = landegf(qnums.l, qnums.j, qnums.f, qnums.i, qnums.s) """ @@ -395,7 +392,7 @@ function zeemanshift(I::Ion, sublevel::Tuple{String, Real}, B::Real) return zeemanshift(B, landegf(I, sublevel[1]), sublevel[2]) + nonlinear end zeemanshift(B::Real, g::Real, m::Real) = (μB / ħ) * g * B * m / 2π -zeemanshift(B::Real, l::Real, j::Real, f::Real, m::Real, i::Real, s::Real = 1 // 2) = +zeemanshift(B::Real, l::Real, j::Real, f::Real, m::Real, i::Real, s::Real=1 // 2) = zeemanshift(B, landegf(l, j, f, i, s), m) zeemanshift(B::Real, qnums::NamedTuple) = zeemanshift(B, qnums.l, qnums.j, qnums.f, qnums.m, qnums.i, qnums.s) @@ -407,7 +404,7 @@ zeemanshift(I::Ion, alias::String, B::Real) = zeemanshift(I, sublevel(I, alias), energy(I::Ion, sublevel; B=0, ignore_manualshift=false) Returns energy of `sublevel` of `I`. A Zeeman shift may be included by setting the value of the magnetic field `B`. The manual shift may be omitted by setting `ignore_manualshift=true`. """ -function energy(I::Ion, sublevel::Tuple{String, Real}; B = 0, ignore_manualshift = false) +function energy(I::Ion, sublevel::Tuple{String, Real}; B=0, ignore_manualshift=false) validatesublevel(I, sublevel) E0 = speciesproperties(I).full_level_structure[sublevel[1]].E zeeman = zeemanshift(I, sublevel, B) @@ -418,7 +415,7 @@ end energy(I::Ion, level::String) Returns the energy of `level` of `I`. """ -function energy(I::Ion, level_or_alias::String; B = 0, ignore_manualshift = false) +function energy(I::Ion, level_or_alias::String; B=0, ignore_manualshift=false) # If the second argument is a String, it could be either a level name or the alias of a sublevel if level_or_alias in levels(I) # Second argument is a level name. Return the bare energy of that level. @@ -428,8 +425,8 @@ function energy(I::Ion, level_or_alias::String; B = 0, ignore_manualshift = fals return energy( I, sublevel(I, level_or_alias), - B = B, - ignore_manualshift = ignore_manualshift + B=B, + ignore_manualshift=ignore_manualshift ) end end @@ -442,11 +439,11 @@ Computes the absolute values of the difference in energies between `transition[1 If between sublevels, then the Zeeman shift may be included by setting the value of the magnetic field `B`, and manual shifts may be omitted by setting `ignore_manualshift=true`. """ -function transitionfrequency(I::Ion, transition::Tuple; B = 0, ignore_manualshift = false) +function transitionfrequency(I::Ion, transition::Tuple; B=0, ignore_manualshift=false) # Multidispatch of the function energy should make this work regardless of whether the transition is between levels or sublevels, and regardless of whether or not aliases are used return abs( - energy(I, transition[1], B = B, ignore_manualshift = ignore_manualshift) - - energy(I, transition[2], B = B, ignore_manualshift = ignore_manualshift) + energy(I, transition[1], B=B, ignore_manualshift=ignore_manualshift) - + energy(I, transition[2], B=B, ignore_manualshift=ignore_manualshift) ) end @@ -454,9 +451,9 @@ end transitionwavelength(I::Ion, transition::Tuple; B=0, ignore_manualshift=false) Returns the wavelength corresponding to `transitionfrequency(I::Ion, transition::Tuple; B=0, ignore_manualshift=false)`. """ -function transitionwavelength(I::Ion, transition::Tuple; B = 0, ignore_manualshift = false) +function transitionwavelength(I::Ion, transition::Tuple; B=0, ignore_manualshift=false) return c / - transitionfrequency(I, transition, B = B, ignore_manualshift = ignore_manualshift) + transitionfrequency(I, transition, B=B, ignore_manualshift=ignore_manualshift) end """ @@ -503,7 +500,7 @@ end einsteinA(I::Ion, leveltransition::Tuple) Returns Einstein A coefficient corresponding to the transition `leveltransition[1] -> leveltransition[2]`. The first level must be the lower level and the second must be the upper level. """ -function einsteinA(I::Ion, leveltransition::Tuple{String,String}) +function einsteinA(I::Ion, leveltransition::Tuple{String, String}) (L1, L2) = leveltransition @assert (L1, L2) in leveltransitions(I) "invalid transition $L1 -> $L2" return speciesproperties(I).full_transitions[(L1, L2)].einsteinA @@ -513,7 +510,7 @@ end transitionmultipole(I::Ion, leveltransition::Tuple) Returns the transition multiple (`'E1'`, `'E2'`, etc.) corresponding to the transition `leveltransition[1] -> leveltransition[2]`. The first level must be the lower level and the second must be the upper level. """ -function transitionmultipole(I::Ion, leveltransition::Tuple{String,String}) +function transitionmultipole(I::Ion, leveltransition::Tuple{String, String}) (L1, L2) = leveltransition @assert (L1, L2) in leveltransitions(I) "invalid transition $L1 -> $L2" return speciesproperties(I).full_transitions[(L1, L2)].multipole @@ -565,7 +562,7 @@ function matrixelement( I::Real, ϵhat::NamedTuple, khat::NamedTuple, - Bhat::NamedTuple = (; z = 1) + Bhat::NamedTuple=(; z=1) ) # Level 1 *must* be the lower level and level 2 *must* be the upper level # Note that in this function, i is the nuclear spin @@ -599,7 +596,7 @@ function matrixelement( geometric_factor = abs( sqrt(2j2 + 1) * wigner3j(f2, 1, f1, -m2, q, m1) * - (transpose(c_rank1[q + 2, :]) * ϵhat_rotated) + (transpose(c_rank1[q+2, :]) * ϵhat_rotated) ) units_factor = abs(e * E / (2ħ) * sqrt(3 * A12 / (α * c * k^3))) return units_factor * hyperfine_factor * geometric_factor / 2π @@ -614,7 +611,7 @@ function matrixelement( geometric_factor = abs( sqrt(2j2 + 1) * wigner3j(f2, 2, f1, -m2, q, m1) * - (transpose(khat_rotated) * c_rank2[:, :, q + 3] * ϵhat_rotated) + (transpose(khat_rotated) * c_rank2[:, :, q+3] * ϵhat_rotated) ) units_factor = abs(e * E / (2ħ) * sqrt(15 * A12 / (α * c * k^3))) return units_factor * hyperfine_factor * geometric_factor / 2π @@ -629,7 +626,7 @@ function matrixelement( I::Real, ϵhat::NamedTuple, khat::NamedTuple, - Bhat::NamedTuple = (; z = 1) + Bhat::NamedTuple=(; z=1) ) SL1 = transition[1] SL2 = transition[2] @@ -787,9 +784,9 @@ IonProperties(; nuclearspin, full_level_structure, full_transitions, - default_sublevel_selection = missing, - gfactors = missing, - nonlinear_zeeman = missing + default_sublevel_selection=missing, + gfactors=missing, + nonlinear_zeeman=missing ) = IonProperties( shortname, mass, @@ -855,8 +852,8 @@ mutable struct IonInstance{Species <: Any} <: Ion # constructors (overrides default) function IonInstance{Species}( properties, - selected_sublevels = nothing, - manualshift = Dict() + selected_sublevels=nothing, + manualshift=Dict() ) where {Species <: Any} (sublevels, aliases) = _construct_sublevels(selected_sublevels, properties) shape = [length(sublevels)] diff --git a/src/iontraps.jl b/src/iontraps.jl index 1a3cc18..e8ef279 100644 --- a/src/iontraps.jl +++ b/src/iontraps.jl @@ -41,7 +41,7 @@ function linear_equilibrium_positions(N::Int) function f!(F, x, N) for i in 1:N F[i] = ( - x[i] - sum([1 / (x[i] - x[j])^2 for j in 1:(i - 1)]) + sum([1 / (x[i] - x[j])^2 for j in (i + 1):N]) + x[i] - sum([1 / (x[i] - x[j])^2 for j in 1:(i-1)]) + sum([1 / (x[i] - x[j])^2 for j in (i+1):N]) ) end end @@ -50,23 +50,23 @@ function linear_equilibrium_positions(N::Int) for i in 1:N, j in 1:N if i ≡ j J[i, j] = ( - 1 + 2 * sum([1 / (x[i] - x[j])^3 for j in 1:(i - 1)]) - 2 * sum([1 / (x[i] - x[j])^3 for j in (i + 1):N]) + 1 + 2 * sum([1 / (x[i] - x[j])^3 for j in 1:(i-1)]) - 2 * sum([1 / (x[i] - x[j])^3 for j in (i+1):N]) ) else J[i, j] = ( - -2 * sum([1 / (x[i] - x[j])^3 for j in 1:(i - 1)]) + 2 * sum([1 / (x[i] - x[j])^3 for j in (i + 1):N]) + -2 * sum([1 / (x[i] - x[j])^3 for j in 1:(i-1)]) + 2 * sum([1 / (x[i] - x[j])^3 for j in (i+1):N]) ) end end end # see eq.8 in the ref to see where (2.018/N^0.559) comes from if isodd(N) - initial_x = [(2.018 / N^0.559) * i for i in (-N ÷ 2):(N ÷ 2)] + initial_x = [(2.018 / N^0.559) * i for i in (-N÷2):(N÷2)] else initial_x = - [(2.018 / N^0.559) * i for i in filter(x -> x ≠ 0, collect((-N ÷ 2):(N ÷ 2)))] + [(2.018 / N^0.559) * i for i in filter(x -> x ≠ 0, collect((-N÷2):(N÷2)))] end - sol = nlsolve((F, x) -> f!(F, x, N), (J, x) -> j!(J, x, N), initial_x, method = :newton) + sol = nlsolve((F, x) -> f!(F, x, N), (J, x) -> j!(J, x, N), initial_x, method=:newton) return sol.zero end @@ -121,7 +121,7 @@ function Anm(N::Int, com::NamedTuple{(:x, :y, :z)}, axis::NamedTuple{(:x, :y, :z end end -_sparsify!(x, eps) = @. x[abs(x) < eps] = 0 +_sparsify!(x, eps) = @. x[abs(x) getfield.(p, x), fieldnames(eltype(p))) @assert length(ion_num) == length(unique(ion_num)) ( "a laser is pointing at the same ion twice" @@ -209,7 +212,7 @@ end Returns the electric field (in V/m) corresponding to a light intensity of `I` (in W/m²) `E = √(2I/(cϵ₀))` """ -efield(I::Real) = √(2I/(c*ϵ₀)) +efield(I::Real) = √(2I / (c * ϵ₀)) """ efield(laser::Laser)::Function @@ -237,4 +240,4 @@ function Base.print(L::Laser) println("k̂: ", "(x=$(L.k.x), y=$(L.k.y), z=$(L.k.z))") println("I(t=0): ", "$(L.I(0.0)) W/m²") return println("ϕ(t=0): ", "$(L.ϕ(0.0)/(2π)) ⋅ 2π") -end \ No newline at end of file +end diff --git a/src/operators.jl b/src/operators.jl index 61cff1f..7f99576 100644 --- a/src/operators.jl +++ b/src/operators.jl @@ -46,7 +46,7 @@ dimension `modecutoff(v)+1`. Otherwise if `method="analytic"`, the matrix elements are computed assuming an infinite-dimension Hilbert space. In general, this option will not return a unitary operator. """ -function displace(v::VibrationalMode, α::Number; method = "truncated") +function displace(v::VibrationalMode, α::Number; method="truncated") # @assert v.N ≥ abs(α) "`α` must be less than `v.N`" # Above line commented out to allow for Hamiltonian construction even if vibrational mode N = 0. # May want to think of a different way to perform this check in the future. @@ -56,8 +56,8 @@ function displace(v::VibrationalMode, α::Number; method = "truncated") return one(v) elseif method ≡ "analytic" @inbounds begin - @simd for n in 1:(modecutoff(v) + 1) - @simd for m in 1:(modecutoff(v) + 1) + @simd for n in 1:(modecutoff(v)+1) + @simd for m in 1:(modecutoff(v)+1) D[n, m] = _Dnm(α, n, m) end end @@ -79,7 +79,7 @@ the thermal density matrix is generated according to the formula: assuming an infinite-dimensional Hilbert space, is used: ``[ρ_{th}]_{ij} = δ_{ij} \\frac{nⁱ}{(n+1)^{i+1}}.`` """ -function thermalstate(v::VibrationalMode, n̄::Real; method = "truncated") +function thermalstate(v::VibrationalMode, n̄::Real; method="truncated") @assert modecutoff(v) ≥ n̄ "`n̄` must be less than `modecutoff(v)`" @assert method in ["truncated", "analytic"] "method ∉ [truncated, analytic]" if n̄ == 0 @@ -88,7 +88,10 @@ function thermalstate(v::VibrationalMode, n̄::Real; method = "truncated") d = [(n̄ / (n̄ + 1))^i for i in 0:(modecutoff(v))] return DenseOperator(v, diagm(0 => d) ./ sum(d)) elseif method ≡ "analytic" - return DenseOperator(v, diagm(0 => [(n̄ / (n̄ + 1))^i / (n̄ + 1) for i in 0:(modecutoff(v))])) + return DenseOperator( + v, + diagm(0 => [(n̄ / (n̄ + 1))^i / (n̄ + 1) for i in 0:(modecutoff(v))]) + ) end end @@ -103,7 +106,7 @@ function coherentstate(v::VibrationalMode, α::Number) k = zeros(ComplexF64, modecutoff(v) + 1) k[1] = exp(-abs2(α) / 2) @inbounds for n in 1:(modecutoff(v)) - k[n + 1] = k[n] * α / √n + k[n+1] = k[n] * α / √n end return Ket(v, k) end @@ -117,13 +120,13 @@ complex amplitude of the displacement. `method` can be either `"truncated"` or `"analytic"` and this argument determines how the displacement operator is computed (see: [`displace`](@ref)) . """ -function coherentthermalstate(v::VibrationalMode, n̄::Real, α::Number; method = "truncated") +function coherentthermalstate(v::VibrationalMode, n̄::Real, α::Number; method="truncated") @assert (modecutoff(v) ≥ n̄ && modecutoff(v) ≥ abs(α)) "`n̄`, `α` must be less than `modecutoff(v)`" @assert method in ["truncated", "analytic"] "method ∉ [truncated, analytic]" if method ≡ "truncated" d = displace(v, α) elseif method ≡ "analytic" - d = displace(v, α, method = "analytic") + d = displace(v, α, method="analytic") end return d * thermalstate(v, n̄) * d' end @@ -167,8 +170,7 @@ function ionstate(iontrap::IonTrap, states::Vector) @assert L ≡ length(states) "wrong number of states" return tensor([ionstate(allions[i], states[i]) for i in 1:L]) end -ionstate(chamber::Chamber, states::Vector) = - ionstate(iontrap(chamber), states) +ionstate(chamber::Chamber, states::Vector) = ionstate(iontrap(chamber), states) """ sigma(ion::Ion, ψ1::sublevel[, ψ2::sublevel]) @@ -195,7 +197,7 @@ If instead `obj<:Chamber`, then this is the same as `obj = Chamber.iontrap`. function ionprojector( IC::IonTrap, sublevels::Union{Tuple{String, Real}, String, Int}...; - only_ions = false + only_ions=false ) allions = ions(IC) L = length(allions) @@ -212,9 +214,9 @@ end function ionprojector( T::Chamber, sublevels::Union{Tuple{String, Real}, String, Int}...; - only_ions = false + only_ions=false ) - return ionprojector(iontrap(T), sublevels..., only_ions = only_ions) + return ionprojector(iontrap(T), sublevels..., only_ions=only_ions) end ############################################################################################# @@ -228,7 +230,7 @@ function _pf(s::Int, n::Int, m::Int) s -= 1 @assert n <= s && m <= s val = 1.0 / (s + 1) - for i in 0:(s - 2) + for i in 0:(s-2) if (m - i > 0) && (n - i > 0) val *= (s - i) / (√((m - i) * (n - i))) elseif m - i > 0 @@ -248,14 +250,14 @@ function _He(n::Int) a[1, 1] = 1 a[2, 1] = 0 a[2, 2] = 1 - for i in 2:(n + 1), j in 1:(n + 1) + for i in 2:(n+1), j in 1:(n+1) if j ≡ 1 - a[i + 1, j] = -(i - 1) * a[i - 1, j] + a[i+1, j] = -(i - 1) * a[i-1, j] else - a[i + 1, j] = a[i, j - 1] - (i - 1) * a[i - 1, j] + a[i+1, j] = a[i, j-1] - (i - 1) * a[i-1, j] end end - return [a[n + 1, k + 1] for k in 0:n] + return [a[n+1, k+1] for k in 0:n] end # computes He_n(x) (nth order Hermite polynomial) @@ -263,7 +265,7 @@ function _fHe(x::Real, n::Int) n -= 1 He = 1.0, x if n < 2 - return He[n + 1] + return He[n+1] end for i in 2:n He = He[2], x * He[2] - (i - 1) * He[1] @@ -297,7 +299,7 @@ end function _alaguerre(x::Real, n::Int, k::Int) L = 1.0, -x + k + 1 if n < 2 - return L[n + 1] + return L[n+1] end for i in 2:n L = L[2], ((k + 2i - 1 - x) * L[2] - (k + i - 1) * L[1]) / i @@ -314,7 +316,7 @@ function _Dnm(ξ::Number, n::Int, m::Int) n -= 1 m -= 1 s = 1.0 - for i in (m + 1):n + for i in (m+1):n s *= i end ret = sqrt(1 / s) * ξ^(n - m) * exp(-abs2(ξ) / 2.0) * _alaguerre(abs2(ξ), m, n - m) diff --git a/src/species/be9.jl b/src/species/be9.jl index ccde254..13e1cc9 100644 --- a/src/species/be9.jl +++ b/src/species/be9.jl @@ -1,35 +1,35 @@ export Be9 const properties_be9 = IonProperties( - shortname = "⁹Be", - mass = 1.496508080073e-26, - charge = 1, - nuclearspin = 3 // 2, - full_level_structure = OrderedDict( - "S1/2f=1" => (n = 2, l = 0, j = 1 // 2, f = 1, E = 0.78126104631e9), - "S1/2f=2" => (n = 2, l = 0, j = 1 // 2, f = 2, E = -0.468756627786e9), - "P1/2f=1" => (n = 2, l = 1, j = 1 // 2, f = 1, E = 957.2010729076436e12), - "P1/2f=2" => (n = 2, l = 1, j = 1 // 2, f = 2, E = 957.2008357076436e12), - "P3/2f=0" => (n = 2, l = 1, j = 3 // 2, f = 0, E = 957.3965669467407e12), - "P3/2f=1" => (n = 2, l = 1, j = 3 // 2, f = 1, E = 957.3965659267407e12), - "P3/2f=2" => (n = 2, l = 1, j = 3 // 2, f = 2, E = 957.3965638867406e12), - "P3/2f=3" => (n = 2, l = 1, j = 3 // 2, f = 3, E = 957.3965608267406e12), + shortname="⁹Be", + mass=1.496508080073e-26, + charge=1, + nuclearspin=3 // 2, + full_level_structure=OrderedDict( + "S1/2f=1" => (n=2, l=0, j=1 // 2, f=1, E=0.78126104631e9), + "S1/2f=2" => (n=2, l=0, j=1 // 2, f=2, E=-0.468756627786e9), + "P1/2f=1" => (n=2, l=1, j=1 // 2, f=1, E=957.2010729076436e12), + "P1/2f=2" => (n=2, l=1, j=1 // 2, f=2, E=957.2008357076436e12), + "P3/2f=0" => (n=2, l=1, j=3 // 2, f=0, E=957.3965669467407e12), + "P3/2f=1" => (n=2, l=1, j=3 // 2, f=1, E=957.3965659267407e12), + "P3/2f=2" => (n=2, l=1, j=3 // 2, f=2, E=957.3965638867406e12), + "P3/2f=3" => (n=2, l=1, j=3 // 2, f=3, E=957.3965608267406e12), ), - full_transitions = Dict( - ("S1/2f=1", "P1/2f=1") => (multipole = "E1", einsteinA = 19.4e6), - ("S1/2f=1", "P1/2f=2") => (multipole = "E1", einsteinA = 19.4e6), - ("S1/2f=2", "P1/2f=1") => (multipole = "E1", einsteinA = 19.4e6), - ("S1/2f=2", "P1/2f=2") => (multipole = "E1", einsteinA = 19.4e6), - ("S1/2f=1", "P3/2f=0") => (multipole = "E1", einsteinA = 19.4e6), - ("S1/2f=1", "P3/2f=1") => (multipole = "E1", einsteinA = 19.4e6), - ("S1/2f=1", "P3/2f=2") => (multipole = "E1", einsteinA = 19.4e6), - ("S1/2f=2", "P1/2f=1") => (multipole = "E1", einsteinA = 19.4e6), - ("S1/2f=2", "P1/2f=2") => (multipole = "E1", einsteinA = 19.4e6), - ("S1/2f=2", "P1/2f=3") => (multipole = "E1", einsteinA = 19.4e6), + full_transitions=Dict( + ("S1/2f=1", "P1/2f=1") => (multipole="E1", einsteinA=19.4e6), + ("S1/2f=1", "P1/2f=2") => (multipole="E1", einsteinA=19.4e6), + ("S1/2f=2", "P1/2f=1") => (multipole="E1", einsteinA=19.4e6), + ("S1/2f=2", "P1/2f=2") => (multipole="E1", einsteinA=19.4e6), + ("S1/2f=1", "P3/2f=0") => (multipole="E1", einsteinA=19.4e6), + ("S1/2f=1", "P3/2f=1") => (multipole="E1", einsteinA=19.4e6), + ("S1/2f=1", "P3/2f=2") => (multipole="E1", einsteinA=19.4e6), + ("S1/2f=2", "P1/2f=1") => (multipole="E1", einsteinA=19.4e6), + ("S1/2f=2", "P1/2f=2") => (multipole="E1", einsteinA=19.4e6), + ("S1/2f=2", "P1/2f=3") => (multipole="E1", einsteinA=19.4e6), ), # Optional fields - default_sublevel_selection = [ + default_sublevel_selection=[ ("S1/2f=1", "all"), ("S1/2f=2", "all"), ("P1/2f=1", "all"), @@ -42,9 +42,7 @@ const properties_be9 = IonProperties( ) # boilerplate code -IonInstance{:Be9}( - selected_sublevels = nothing, - manualshift = Dict() -) = IonInstance{:Be9}(properties_be9, selected_sublevels, manualshift) +IonInstance{:Be9}(selected_sublevels=nothing, manualshift=Dict()) = + IonInstance{:Be9}(properties_be9, selected_sublevels, manualshift) Be9 = IonInstance{:Be9} diff --git a/src/species/ca40.jl b/src/species/ca40.jl index 27c781c..c9a2c68 100644 --- a/src/species/ca40.jl +++ b/src/species/ca40.jl @@ -1,39 +1,37 @@ export Ca40 const properties_ca40 = IonProperties(; - shortname = "⁴⁰Ca", - mass = 6.635943757345042e-26, - charge = 1, - nuclearspin = 0, - full_level_structure = OrderedDict( - "S1/2" => (n = 4, l = 0, j = 1 // 2, f = 1 // 2, E = 0), - "D3/2" => (n = 3, l = 2, j = 3 // 2, f = 3 // 2, E = 4.09335071228e14), - "D5/2" => (n = 3, l = 2, j = 5 // 2, f = 5 // 2, E = 4.1115503183857306e14), - "P1/2" => (n = 4, l = 1, j = 1 // 2, f = 1 // 2, E = 7.554e14), - "P3/2" => (n = 4, l = 1, j = 3 // 2, f = 3 // 2, E = 7.621e14), + shortname="⁴⁰Ca", + mass=6.635943757345042e-26, + charge=1, + nuclearspin=0, + full_level_structure=OrderedDict( + "S1/2" => (n=4, l=0, j=1 // 2, f=1 // 2, E=0), + "D3/2" => (n=3, l=2, j=3 // 2, f=3 // 2, E=4.09335071228e14), + "D5/2" => (n=3, l=2, j=5 // 2, f=5 // 2, E=4.1115503183857306e14), + "P1/2" => (n=4, l=1, j=1 // 2, f=1 // 2, E=7.554e14), + "P3/2" => (n=4, l=1, j=3 // 2, f=3 // 2, E=7.621e14), ), - full_transitions = Dict( - ("S1/2", "D5/2") => (multipole = "E2", einsteinA = 8.562e-1), - ("S1/2", "P1/2") => (multipole = "E1", einsteinA = 1.299e8), - ("D3/2", "P1/2") => (multipole = "E1", einsteinA = 1.060e7), - ("S1/2", "D3/2") => (multipole = "E2", einsteinA = 9.259e-1), - ("S1/2", "P3/2") => (multipole = "E1", einsteinA = 1.351e8), - ("D3/2", "P3/2") => (multipole = "E1", einsteinA = 1.110e6), - ("D5/2", "P3/2") => (multipole = "E1", einsteinA = 9.901e6), + full_transitions=Dict( + ("S1/2", "D5/2") => (multipole="E2", einsteinA=8.562e-1), + ("S1/2", "P1/2") => (multipole="E1", einsteinA=1.299e8), + ("D3/2", "P1/2") => (multipole="E1", einsteinA=1.060e7), + ("S1/2", "D3/2") => (multipole="E2", einsteinA=9.259e-1), + ("S1/2", "P3/2") => (multipole="E1", einsteinA=1.351e8), + ("D3/2", "P3/2") => (multipole="E1", einsteinA=1.110e6), + ("D5/2", "P3/2") => (multipole="E1", einsteinA=9.901e6), ), # Optional fields - default_sublevel_selection = [("S1/2", "all"), ("D5/2", "all"),], - gfactors = Dict("S1/2" => 2.00225664, "D5/2" => 1.2003340) + default_sublevel_selection=[("S1/2", "all"), ("D5/2", "all"),], + gfactors=Dict("S1/2" => 2.00225664, "D5/2" => 1.2003340) #nonlinear_zeeman = Dict(("S1/2", -1//2) => B->1.3e-4*B^2, # ("D5/2", -5//2) => B->4.5e-4*B^2) # Syntax example, not numerically accurate ) # boilerplate code -IonInstance{:Ca40}( - selected_sublevels = nothing, - manualshift = Dict() -) = IonInstance{:Ca40}(properties_ca40, selected_sublevels, manualshift) +IonInstance{:Ca40}(selected_sublevels=nothing, manualshift=Dict()) = + IonInstance{:Ca40}(properties_ca40, selected_sublevels, manualshift) Ca40 = IonInstance{:Ca40} diff --git a/src/species/mg25.jl b/src/species/mg25.jl index cd99b62..b27ce16 100644 --- a/src/species/mg25.jl +++ b/src/species/mg25.jl @@ -1,34 +1,34 @@ export Mg25 const properties_mg25 = IonProperties( - shortname = "²⁵Mg", - mass = 4.1489954812e-26, - charge = 1, - nuclearspin = 5 // 2, - full_level_structure = OrderedDict( - "S1/2f=2" => (n = 3, l = 0, j = 1 // 2, f = 2, E = 1.043445158e9), - "S1/2f=3" => (n = 3, l = 0, j = 1 // 2, f = 3, E = -0.745317970e9), - "P1/2f=2" => (n = 3, l = 1, j = 1 // 2, f = 2, E = 1069.3408690519877e12), - "P1/2f=3" => (n = 3, l = 1, j = 1 // 2, f = 3, E = 1069.3405639519879e12), - "P3/2f=1" => (n = 3, l = 1, j = 3 // 2, f = 1, E = 1072.0853411194748e12), - "P3/2f=2" => (n = 3, l = 1, j = 3 // 2, f = 2, E = 1072.0853033394746e12), - "P3/2f=3" => (n = 3, l = 1, j = 3 // 2, f = 3, E = 1072.0852466694748e12), - "P3/2f=4" => (n = 3, l = 1, j = 3 // 2, f = 4, E = 1072.0851711094747e12), + shortname="²⁵Mg", + mass=4.1489954812e-26, + charge=1, + nuclearspin=5 // 2, + full_level_structure=OrderedDict( + "S1/2f=2" => (n=3, l=0, j=1 // 2, f=2, E=1.043445158e9), + "S1/2f=3" => (n=3, l=0, j=1 // 2, f=3, E=-0.745317970e9), + "P1/2f=2" => (n=3, l=1, j=1 // 2, f=2, E=1069.3408690519877e12), + "P1/2f=3" => (n=3, l=1, j=1 // 2, f=3, E=1069.3405639519879e12), + "P3/2f=1" => (n=3, l=1, j=3 // 2, f=1, E=1072.0853411194748e12), + "P3/2f=2" => (n=3, l=1, j=3 // 2, f=2, E=1072.0853033394746e12), + "P3/2f=3" => (n=3, l=1, j=3 // 2, f=3, E=1072.0852466694748e12), + "P3/2f=4" => (n=3, l=1, j=3 // 2, f=4, E=1072.0851711094747e12), ), - full_transitions = Dict( - ("S1/2f=2", "P1/2f=2") => (multipole = "E1", einsteinA = 41.3e6), - ("S1/2f=2", "P1/2f=3") => (multipole = "E1", einsteinA = 41.3e6), - ("S1/2f=3", "P1/2f=2") => (multipole = "E1", einsteinA = 41.3e6), - ("S1/2f=3", "P1/2f=3") => (multipole = "E1", einsteinA = 41.3e6), - ("S1/2f=2", "P3/2f=1") => (multipole = "E1", einsteinA = 41.3e6), - ("S1/2f=2", "P3/2f=2") => (multipole = "E1", einsteinA = 41.3e6), - ("S1/2f=2", "P3/2f=3") => (multipole = "E1", einsteinA = 41.3e6), - ("S1/2f=3", "P1/2f=2") => (multipole = "E1", einsteinA = 41.3e6), - ("S1/2f=3", "P1/2f=3") => (multipole = "E1", einsteinA = 41.3e6), - ("S1/2f=3", "P1/2f=4") => (multipole = "E1", einsteinA = 41.3e6), + full_transitions=Dict( + ("S1/2f=2", "P1/2f=2") => (multipole="E1", einsteinA=41.3e6), + ("S1/2f=2", "P1/2f=3") => (multipole="E1", einsteinA=41.3e6), + ("S1/2f=3", "P1/2f=2") => (multipole="E1", einsteinA=41.3e6), + ("S1/2f=3", "P1/2f=3") => (multipole="E1", einsteinA=41.3e6), + ("S1/2f=2", "P3/2f=1") => (multipole="E1", einsteinA=41.3e6), + ("S1/2f=2", "P3/2f=2") => (multipole="E1", einsteinA=41.3e6), + ("S1/2f=2", "P3/2f=3") => (multipole="E1", einsteinA=41.3e6), + ("S1/2f=3", "P1/2f=2") => (multipole="E1", einsteinA=41.3e6), + ("S1/2f=3", "P1/2f=3") => (multipole="E1", einsteinA=41.3e6), + ("S1/2f=3", "P1/2f=4") => (multipole="E1", einsteinA=41.3e6), ), # Optional fields - default_sublevel_selection = [ + default_sublevel_selection=[ ("S1/2f=2", "all"), ("S1/2f=3", "all"), ("P1/2f=2", "all"), @@ -41,9 +41,7 @@ const properties_mg25 = IonProperties( ) # boilerplate code -IonInstance{:Mg25}( - selected_sublevels = nothing, - manualshift = Dict() -) = IonInstance{:Mg25}(properties_mg25, selected_sublevels, manualshift) +IonInstance{:Mg25}(selected_sublevels=nothing, manualshift=Dict()) = + IonInstance{:Mg25}(properties_mg25, selected_sublevels, manualshift) Mg25 = IonInstance{:Mg25} diff --git a/src/species/yb171.jl b/src/species/yb171.jl index 55b619c..273ced9 100644 --- a/src/species/yb171.jl +++ b/src/species/yb171.jl @@ -1,78 +1,72 @@ export Yb171 const properties_yb171 = IonProperties( - shortname = "¹⁷¹Yb", - mass = 28.8384644689030595108e-26, - charge = 1, - nuclearspin = 1 // 2, - full_level_structure = OrderedDict( + shortname="¹⁷¹Yb", + mass=28.8384644689030595108e-26, + charge=1, + nuclearspin=1 // 2, + full_level_structure=OrderedDict( # should this be negative, or should I increase the energy of every state by 9GHz? - "S1/2f=0" => (n = 6, l = 0, j = 1 // 2, f = 0, E = -9.48210909315e9), - "S1/2f=1" => (n = 6, l = 0, j = 1 // 2, f = 1, E = 3.16070303105e9), - "D3/2f=0" => (n = 5, l = 2, j = 3 // 2, f = 1, E = 6.88342460964640e14), - "D3/2f=1" => (n = 5, l = 2, j = 3 // 2, f = 2, E = 6.88350470564640e14), - "D5/2f=2" => (n = 5, l = 2, j = 5 // 2, f = 2, E = 7.29477895985202e14), - "D5/2f=3" => (n = 5, l = 2, j = 5 // 2, f = 3, E = 7.29474121985202e14), - "P1/2f=0" => (n = 6, l = 1, j = 1 // 2, f = 0, E = 8.112913749003559e14), - "P1/2f=1" => (n = 6, l = 1, j = 1 // 2, f = 1, E = 8.112934798003559e14), - "[3/2]1/2f=0" => - (n = 6, l = nothing, j = 1 // 2, f = 0, E = 1.008917341058788e15), - "[3/2]1/2f=1" => - (n = 6, l = nothing, j = 1 // 2, f = 1, E = 1.008917341058788e15), - "[3/2]3/2f=1" => - (n = 6, l = nothing, j = 3 // 2, f = 1, E = 8.621425511314839e14), - "[3/2]3/2f=2" => - (n = 6, l = nothing, j = 3 // 2, f = 2, E = 8.621425511314839e14), - "[5/2]5/2f=2" => - (n = 6, l = nothing, j = 5 // 2, f = 2, E = 9.70461163716380e14), - "[5/2]5/2f=3" => - (n = 6, l = nothing, j = 5 // 2, f = 3, E = 9.70461163716380e14), - "F7/2f=3" => (n = 6, l = 3, j = 7 // 2, f = 3, E = 6.42115934728750e14), - "F7/2f=4" => (n = 6, l = 3, j = 7 // 2, f = 4, E = 6.42119554728750e14), + "S1/2f=0" => (n=6, l=0, j=1 // 2, f=0, E=-9.48210909315e9), + "S1/2f=1" => (n=6, l=0, j=1 // 2, f=1, E=3.16070303105e9), + "D3/2f=0" => (n=5, l=2, j=3 // 2, f=1, E=6.88342460964640e14), + "D3/2f=1" => (n=5, l=2, j=3 // 2, f=2, E=6.88350470564640e14), + "D5/2f=2" => (n=5, l=2, j=5 // 2, f=2, E=7.29477895985202e14), + "D5/2f=3" => (n=5, l=2, j=5 // 2, f=3, E=7.29474121985202e14), + "P1/2f=0" => (n=6, l=1, j=1 // 2, f=0, E=8.112913749003559e14), + "P1/2f=1" => (n=6, l=1, j=1 // 2, f=1, E=8.112934798003559e14), + "[3/2]1/2f=0" => (n=6, l=nothing, j=1 // 2, f=0, E=1.008917341058788e15), + "[3/2]1/2f=1" => (n=6, l=nothing, j=1 // 2, f=1, E=1.008917341058788e15), + "[3/2]3/2f=1" => (n=6, l=nothing, j=3 // 2, f=1, E=8.621425511314839e14), + "[3/2]3/2f=2" => (n=6, l=nothing, j=3 // 2, f=2, E=8.621425511314839e14), + "[5/2]5/2f=2" => (n=6, l=nothing, j=5 // 2, f=2, E=9.70461163716380e14), + "[5/2]5/2f=3" => (n=6, l=nothing, j=5 // 2, f=3, E=9.70461163716380e14), + "F7/2f=3" => (n=6, l=3, j=7 // 2, f=3, E=6.42115934728750e14), + "F7/2f=4" => (n=6, l=3, j=7 // 2, f=4, E=6.42119554728750e14), ), - full_transitions = Dict( + full_transitions=Dict( # Laser lines # 411nm - ("S1/2f=0", "D5/2f=2") => (multipole = "E2", einsteinA = 22), - ("S1/2f=1", "D5/2f=2") => (multipole = "E2", einsteinA = 22), - ("S1/2f=1", "D5/2f=3") => (multipole = "E2", einsteinA = 22), + ("S1/2f=0", "D5/2f=2") => (multipole="E2", einsteinA=22), + ("S1/2f=1", "D5/2f=2") => (multipole="E2", einsteinA=22), + ("S1/2f=1", "D5/2f=3") => (multipole="E2", einsteinA=22), # 369nm - ("S1/2f=0", "P1/2f=1") => (multipole = "E1", einsteinA = 1.155e8), - ("S1/2f=1", "P1/2f=0") => (multipole = "E1", einsteinA = 1.155e8), - ("S1/2f=1", "P1/2f=1") => (multipole = "E1", einsteinA = 1.155e8), + ("S1/2f=0", "P1/2f=1") => (multipole="E1", einsteinA=1.155e8), + ("S1/2f=1", "P1/2f=0") => (multipole="E1", einsteinA=1.155e8), + ("S1/2f=1", "P1/2f=1") => (multipole="E1", einsteinA=1.155e8), # 935nm - ("D3/2f=1", "[3/2]1/2f=0") => (multipole = "E1", einsteinA = 1.2e5), - ("D3/2f=1", "[3/2]1/2f=1") => (multipole = "E1", einsteinA = 1.2e5), - ("D3/2f=2", "[3/2]1/2f=1") => (multipole = "E1", einsteinA = 1.2e5), + ("D3/2f=1", "[3/2]1/2f=0") => (multipole="E1", einsteinA=1.2e5), + ("D3/2f=1", "[3/2]1/2f=1") => (multipole="E1", einsteinA=1.2e5), + ("D3/2f=2", "[3/2]1/2f=1") => (multipole="E1", einsteinA=1.2e5), # 760nm - ("F7/2f=3", "[3/2]3/2f=1") => (multipole = "E2", einsteinA = 5e4), - ("F7/2f=3", "[3/2]3/2f=2") => (multipole = "E2", einsteinA = 5e4), - ("F7/2f=4", "[3/2]3/2f=2") => (multipole = "E2", einsteinA = 5e4), + ("F7/2f=3", "[3/2]3/2f=1") => (multipole="E2", einsteinA=5e4), + ("F7/2f=3", "[3/2]3/2f=2") => (multipole="E2", einsteinA=5e4), + ("F7/2f=4", "[3/2]3/2f=2") => (multipole="E2", einsteinA=5e4), # 638nm # 976nm # 861nm # Decay lines # P1/2 -> D3/2 - ("P1/2f=0", "D3/2f=1") => (multipole = "E1", einsteinA = 5.7e5), - ("P1/2f=1", "D3/2f=2") => (multipole = "E1", einsteinA = 5.7e5), - ("P1/2f=1", "D3/2f=2") => (multipole = "E1", einsteinA = 5.7e5), + ("P1/2f=0", "D3/2f=1") => (multipole="E1", einsteinA=5.7e5), + ("P1/2f=1", "D3/2f=2") => (multipole="E1", einsteinA=5.7e5), + ("P1/2f=1", "D3/2f=2") => (multipole="E1", einsteinA=5.7e5), # [3/2]1/2 -> S1/2 - ("[3/2]1/2f=0", "S1/2f=1") => (multipole = "E1", einsteinA = 8.05e7), - ("[3/2]1/2f=1", "S1/2f=0") => (multipole = "E1", einsteinA = 8.05e7), - ("[3/2]1/2f=1", "S1/2f=1") => (multipole = "E1", einsteinA = 8.05e7), + ("[3/2]1/2f=0", "S1/2f=1") => (multipole="E1", einsteinA=8.05e7), + ("[3/2]1/2f=1", "S1/2f=0") => (multipole="E1", einsteinA=8.05e7), + ("[3/2]1/2f=1", "S1/2f=1") => (multipole="E1", einsteinA=8.05e7), # [3/2]3/2 -> S1/2 - ("[3/2]1/2f=0", "S1/2f=1") => (multipole = "E1", einsteinA = 5.125e7), - ("[3/2]1/2f=1", "S1/2f=0") => (multipole = "E1", einsteinA = 5.125e7), - ("[3/2]1/2f=1", "S1/2f=1") => (multipole = "E1", einsteinA = 5.125e7), + ("[3/2]1/2f=0", "S1/2f=1") => (multipole="E1", einsteinA=5.125e7), + ("[3/2]1/2f=1", "S1/2f=0") => (multipole="E1", einsteinA=5.125e7), + ("[3/2]1/2f=1", "S1/2f=1") => (multipole="E1", einsteinA=5.125e7), ), # Optional fields - default_sublevel_selection = [ + default_sublevel_selection=[ ("S1/2f=0", "all"), ("S1/2f=1", "all"), ("P1/2f=0", "all"), ("P1/2f=1", "all") ], - gfactors = Dict( + gfactors=Dict( "S1/2f=0" => 1.998, "S1/2f=1" => 1.998, "D5/2f=2" => 1.202, @@ -84,16 +78,14 @@ const properties_yb171 = IonProperties( "[3/2]3/2f=1" => 1.44, "[3/2]3/2f=2" => 1.44, ), - nonlinear_zeeman = Dict( + nonlinear_zeeman=Dict( ("S1/2f=0", 0) => B -> -155.305 * B^2, ("S1/2f=1", 0) => B -> 155.305 * B^2 ) ) # boilerplate code -IonInstance{:Yb171}( - selected_sublevels = nothing, - manualshift = Dict() -) = IonInstance{:Yb171}(properties_yb171, selected_sublevels, manualshift) +IonInstance{:Yb171}(selected_sublevels=nothing, manualshift=Dict()) = + IonInstance{:Yb171}(properties_yb171, selected_sublevels, manualshift) Yb171 = IonInstance{:Yb171} diff --git a/src/timeevolution.jl b/src/timeevolution.jl index 6a96776..035cbbc 100644 --- a/src/timeevolution.jl +++ b/src/timeevolution.jl @@ -58,7 +58,7 @@ function timeevolution.schroedinger_dynamic( tspan, rho0::T, f::Function; - fout::Union{Function, Nothing} = nothing, + fout::Union{Function, Nothing}=nothing, kwargs... ) where {B <: Basis, T <: Operator{B, B}} check_bases(f(0.0, 0.0).basis_l, rho0.basis_l) @@ -67,7 +67,7 @@ function timeevolution.schroedinger_dynamic( tspan, rho0, (t, rho) -> (f(t, rho), Jvec, Jvec); - fout = fout, + fout=fout, kwargs... ) end diff --git a/src/vibrationalmodes.jl b/src/vibrationalmodes.jl index 61d5f50..bfaf7e3 100644 --- a/src/vibrationalmodes.jl +++ b/src/vibrationalmodes.jl @@ -40,13 +40,7 @@ mutable struct VibrationalMode <: IonSimBasis shape::Vector{Int} axis::NamedTuple{(:x, :y, :z)} _cnst_δν::Bool - function VibrationalMode( - ν, - modestructure; - δν::Tδν = 0.0, - N = 10, - axis = ẑ - ) where {Tδν} + function VibrationalMode(ν, modestructure; δν::Tδν=0.0, N=10, axis=ẑ) where {Tδν} if Tδν <: Number _cnst_δν = true δνt(t) = δν @@ -134,7 +128,7 @@ Sets the upper bound of the Hilbert space of `mode` to be the Fock state `N`. function modecutoff!(mode::VibrationalMode, N::Int) @assert N >= 0 "N must be a nonnegative integer" mode.N = N - mode.shape = Int[N + 1] + mode.shape = Int[N+1] end diff --git a/test/test_analytic_functions.jl b/test/test_analytic_functions.jl index 1c3c02c..9ceda64 100644 --- a/test/test_analytic_functions.jl +++ b/test/test_analytic_functions.jl @@ -44,5 +44,5 @@ end # test bluesideband, Ω̃ ≈ ηΩ η = 0.1 - @test rabiflop(t, Ω, η, n̄, s = 1) ≈ @.(sin(2π * η * t / 2)^2) rtol = 1e-2 + @test rabiflop(t, Ω, η, n̄, s=1) ≈ @.(sin(2π * η * t / 2)^2) rtol = 1e-2 end diff --git a/test/test_chambers.jl b/test/test_chambers.jl index 036ef71..1bfb5e3 100644 --- a/test/test_chambers.jl +++ b/test/test_chambers.jl @@ -8,16 +8,16 @@ using Suppressor # set up system C = Ca40() λ = transitionwavelength(C, ("S1/2", "D5/2")) - L1 = Laser(λ = λ) - L2 = Laser(λ = λ) + L1 = Laser(λ=λ) + L2 = Laser(λ=λ) chain = LinearChain( - ions = [C, C], - comfrequencies = (x = 3e6, y = 3e6, z = 1e6), - selectedmodes = (; z = [1]) + ions=[C, C], + comfrequencies=(x=3e6, y=3e6, z=1e6), + selectedmodes=(; z=[1]) ) @testset "chambers -- Chamber" begin - T = Chamber(iontrap = chain, lasers = [L1, L1]) + T = Chamber(iontrap=chain, lasers=[L1, L1]) # test construction of CompositeBasis @test basis(T).bases[1] ≡ T.iontrap.ions[1] @@ -36,10 +36,10 @@ using Suppressor # test for warning when lasers=[L, L, L, ...] where L point to the same thing warn = "Some lasers point to the same thing. Making copies." - @test_logs (:warn, warn) Chamber(iontrap = chain, lasers = [L1, L1, L1]) + @test_logs (:warn, warn) Chamber(iontrap=chain, lasers=[L1, L1, L1]) # and test that, in this case, copies are made - T1 = Chamber(iontrap = chain, lasers = [L1, L1, L1, L1], δB = t -> t) - for i in 1:4, j in (i + 1):4 + T1 = Chamber(iontrap=chain, lasers=[L1, L1, L1, L1], δB=t -> t) + for i in 1:4, j in (i+1):4 @test !(T1.lasers[i] ≡ T1.lasers[j]) end @@ -53,9 +53,9 @@ using Suppressor # @test basis(T) == basis # changing the iontrap should also update the basis chain1 = LinearChain( - ions = [C, C, C], - comfrequencies = (x = 3e6, y = 3e6, z = 1e6), - selectedmodes = (; z = [1]) + ions=[C, C, C], + comfrequencies=(x=3e6, y=3e6, z=1e6), + selectedmodes=(; z=[1]) ) @test length(basis(T).bases) ≡ 3 iontrap!(T, chain1) @@ -69,10 +69,10 @@ using Suppressor @testset "chambers -- general functions" begin Bhat = (x̂ + ŷ + ẑ) / √3 - T = Chamber(iontrap = chain, lasers = [L1], Bhat = Bhat) + T = Chamber(iontrap=chain, lasers=[L1], Bhat=Bhat) # globalbeam! - L = Laser(λ = λ) + L = Laser(λ=λ) globalbeam!(L, T) @test L.pointing == [(1, 1.0), (2, 1.0)] @@ -92,7 +92,7 @@ using Suppressor intensity_from_pitime!(laser_index, 1e-6, ion_index, transition, T) @test L1.I(0) == I1 # shouldn't be able to have a laser argument where laser.pointing = [] - L = Laser(λ = λ) + L = Laser(λ=λ) @test_throws AssertionError intensity_from_pitime(L, 1e-6, ion, transition, Bhat) L.pointing = [(1, 1.0)] @test isinf(intensity_from_pitime(L, 1e-6, ion, transition, x̂)) @@ -113,7 +113,7 @@ using Suppressor intensity_from_rabifrequency!(laser_index, 5e5, ion_index, transition, T) @test L1.I(0) == I1 # shouldn't be able to have a laser argument where laser.pointing = [] - L = Laser(λ = λ) + L = Laser(λ=λ) @test_throws AssertionError intensity_from_rabifrequency( L, 5e5, @@ -126,7 +126,7 @@ using Suppressor T.B = 4e-4 f = transitionfrequency(C, transition, T) @test f ≈ 4.111550340833542e14 - @test transitionfrequency(C, transition, B = T.B) == f + @test transitionfrequency(C, transition, B=T.B) == f @test transitionfrequency(1, transition, T) == f # set_gradient @@ -137,11 +137,11 @@ using Suppressor # test :(==) chain1 = LinearChain( - ions = [C, C], - comfrequencies = (x = 3e6, y = 3e6, z = 1e6), - selectedmodes = (x = [1], z = [1]) + ions=[C, C], + comfrequencies=(x=3e6, y=3e6, z=1e6), + selectedmodes=(x=[1], z=[1]) ) - T = Chamber(iontrap = chain1, lasers = [L1, L2]) + T = Chamber(iontrap=chain1, lasers=[L1, L2]) xmode = chain1.selectedmodes.x[1] zmode = chain1.selectedmodes.z[1] cb = basis(T) diff --git a/test/test_constants.jl b/test/test_constants.jl index 74d511f..34f8df2 100644 --- a/test/test_constants.jl +++ b/test/test_constants.jl @@ -46,15 +46,15 @@ using Test, IonSim, IonSim.PhysicalConstants, Suppressor @test IonSim._print_axis(x̂) == "x̂" @test IonSim._print_axis(ŷ) == "ŷ" @test IonSim._print_axis(ẑ) == "ẑ" - @test IonSim._print_axis((x = 1, y = 1, z = 1)) == string((x = 1, y = 1, z = 1)) + @test IonSim._print_axis((x=1, y=1, z=1)) == string((x=1, y=1, z=1)) end @testset "constants -- other constants" begin - @test x̂ / 2 == (x = 1 / 2, y = 0, z = 0) - @test 2x̂ == (x = 2, y = 0, z = 0) - @test x̂ * 2 == (x = 2, y = 0, z = 0) - @test x̂ + x̂ == (x = 2, y = 0, z = 0) - @test x̂ - x̂ == (x = 0, y = 0, z = 0) + @test x̂ / 2 == (x=1 / 2, y=0, z=0) + @test 2x̂ == (x=2, y=0, z=0) + @test x̂ * 2 == (x=2, y=0, z=0) + @test x̂ + x̂ == (x=2, y=0, z=0) + @test x̂ - x̂ == (x=0, y=0, z=0) @test ndot(x̂, x̂) == 1 @test ndot(x̂, ŷ) == 0 end diff --git a/test/test_dynamics.jl b/test/test_dynamics.jl index 43adb21..023d459 100644 --- a/test/test_dynamics.jl +++ b/test/test_dynamics.jl @@ -9,11 +9,11 @@ using Suppressor C = Ca40([("S1/2", -1 / 2), ("D5/2", -1 / 2)]) L = Laser() chain = LinearChain( - ions = [C], - comfrequencies = (x = 3e6, y = 3e6, z = 1e6), - selectedmodes = (; z = [1]) + ions=[C], + comfrequencies=(x=3e6, y=3e6, z=1e6), + selectedmodes=(; z=[1]) ) - T = Chamber(iontrap = chain, B = 4e-4, Bhat = ẑ, δB = 0, lasers = [L]) + T = Chamber(iontrap=chain, B=4e-4, Bhat=ẑ, δB=0, lasers=[L]) L.λ = transitionwavelength(C, (("S1/2", -1 / 2), ("D5/2", -1 / 2)), T) mode = T.iontrap.selectedmodes.z[1] L.k = (x̂ + ẑ) / √2 @@ -31,7 +31,7 @@ using Suppressor ) ex_ionsim_c0 = real.(expect(ionprojector(T, ("D5/2", -1 / 2)), sol)) ex_analyt_c0 = @.(sin(2π * Ω00 / 2 * tout * 1e-6)^2) - @test isapprox(ex_ionsim_c0, ex_analyt_c0, rtol = 1e-5) + @test isapprox(ex_ionsim_c0, ex_analyt_c0, rtol=1e-5) # This test serves to check for the presence of a bug -- The norm of the output seems to deviate from 1, moreso as the simulation goes on # If this test fails then the problem has likely been solved. @@ -50,7 +50,7 @@ using Suppressor ex_ionsim_d1 = real.(expect(ionprojector(T, ("D5/2", -1 / 2)), sol)) ex_analyt_d1 = @.((Ω00^2 / (Ω00^2 + Δ^2)) * sin(2π * √(Ω00^2 + Δ^2) / 2 * tout * 1e-6)^2) - @test isapprox(ex_ionsim_d1, ex_analyt_d1, rtol = 1e-4) + @test isapprox(ex_ionsim_d1, ex_analyt_d1, rtol=1e-4) # add detuning using ion's manualshift L.Δ = 0 @@ -65,7 +65,7 @@ using Suppressor ex_ionsim_d2 = real.(expect(ionprojector(T, ("D5/2", -1 / 2)), sol)) ex_analyt_d2 = @.((Ω00^2 / (Ω00^2 + Δ^2)) * sin(2π * √(Ω00^2 + Δ^2) / 2 * tout * 1e-6)^2) - @test isapprox(ex_ionsim_d2, ex_analyt_d2, rtol = 1e-4) + @test isapprox(ex_ionsim_d2, ex_analyt_d2, rtol=1e-4) # hot carrier # For this test, numerical result deviates from analytical quickly as nbar grows. @@ -76,19 +76,19 @@ using Suppressor n̄ = rand(1:10) ψi_mode = thermalstate(mode, n̄) ψi = ψi_ion ⊗ ψi_mode - h = hamiltonian(T, timescale=1e-6, lamb_dicke_order = 0) + h = hamiltonian(T, timescale=1e-6, lamb_dicke_order=0) tout, sol = timeevolution.schroedinger_dynamic(tspan, ψi, h) η = lambdicke(mode, C, L) ex_ionsim_cn = real.(expect(ionprojector(T, ("D5/2", -1 / 2)), sol)) ex_analyt_cn = analytical.rabiflop(1e-6 * tout, Ω, η, n̄) - @test isapprox(ex_ionsim_cn, ex_analyt_cn, rtol = 1e-2) + @test isapprox(ex_ionsim_cn, ex_analyt_cn, rtol=1e-2) # sideband transitions ## under an RWA, RSB transitions, when starting in motional ground state, ## should be suppressed L.Δ = -mode.ν - h = hamiltonian(T, timescale=1e-6, rwa_cutoff = 1e3) + h = hamiltonian(T, timescale=1e-6, rwa_cutoff=1e3) tspan_sb = 0:1:2000 tout, sol = timeevolution.schroedinger_dynamic( tspan_sb, @@ -101,7 +101,7 @@ using Suppressor ## a BSB should have a coupling strength ηΩ√(n+1) ## a RSB should have a coupling strength ηΩ√n L.Δ = -mode.ν - h = hamiltonian(T, timescale=1e-6, rwa_cutoff = 1e3) + h = hamiltonian(T, timescale=1e-6, rwa_cutoff=1e3) tout, sol = timeevolution.schroedinger_dynamic( tspan_sb, ionstate(T, [("S1/2", -1 / 2)]) ⊗ mode[1], @@ -109,10 +109,10 @@ using Suppressor ) ex_ionsim_rsb1 = real.(expect(ionprojector(T, ("D5/2", -1 / 2)), sol)) ex_analyt_rsb1 = @.(sin(2π * η * Ω00 / 2 * tout * 1e-6)^2) - @test isapprox(ex_ionsim_rsb1, ex_analyt_rsb1, rtol = 1e-5) + @test isapprox(ex_ionsim_rsb1, ex_analyt_rsb1, rtol=1e-5) L.Δ = mode.ν - h = hamiltonian(T, timescale=1e-6, rwa_cutoff = 1e3) + h = hamiltonian(T, timescale=1e-6, rwa_cutoff=1e3) tout, sol = timeevolution.schroedinger_dynamic( tspan, ionstate(T, [("S1/2", -1 / 2)]) ⊗ mode[1], @@ -120,7 +120,7 @@ using Suppressor ) ex_ionsim_bsb1 = real.(expect(ionprojector(T, ("D5/2", -1 / 2)), sol)) ex_analyt_bsb1 = @.(sin(2π * sqrt(2) * η * Ω00 / 2 * tout * 1e-6)^2) - @test isapprox(ex_ionsim_bsb1, ex_analyt_bsb1, rtol = 2e-2) + @test isapprox(ex_ionsim_bsb1, ex_analyt_bsb1, rtol=2e-2) end @testset "Dynamics -- time-dependent-parameters" begin @@ -130,11 +130,11 @@ using Suppressor C = Ca40([("S1/2", -1 / 2), ("D5/2", -1 / 2)]) L = Laser() chain = LinearChain( - ions = [C], - comfrequencies = (x = 3e6, y = 3e6, z = 1e6), - selectedmodes = (; z = [1]) + ions=[C], + comfrequencies=(x=3e6, y=3e6, z=1e6), + selectedmodes=(; z=[1]) ) - T = Chamber(iontrap = chain, B = 4e-4, Bhat = ẑ, δB = 0, lasers = [L]) + T = Chamber(iontrap=chain, B=4e-4, Bhat=ẑ, δB=0, lasers=[L]) L.λ = transitionwavelength(C, ("S1/2", "D5/2"), T) mode = T.iontrap.selectedmodes.z[1] L.k = (x̂ + ẑ) / √2 @@ -147,7 +147,7 @@ using Suppressor Ω00 = Ω * exp(-lambdicke(mode, C, L)^2 / 2) intensity_from_rabifrequency!(1, Ω, 1, (("S1/2", -1 / 2), ("D5/2", -1 / 2)), T) - h = hamiltonian(T, timescale=1e-6, lamb_dicke_order = 0) + h = hamiltonian(T, timescale=1e-6, lamb_dicke_order=0) tspan = 0:1e-3:4 tout, sol = timeevolution.schroedinger_dynamic( tspan, @@ -156,14 +156,14 @@ using Suppressor ) ex_ionsim_tdphi = real.(expect(ionprojector(T, ("D5/2", -1 / 2)), sol)) ex_analyt_tdphi = @.(sin(2π * 1e-6 * Ω00 / 2 * tout)^2) - @test isapprox(ex_ionsim_tdphi, ex_analyt_tdphi, rtol = 1e-5) + @test isapprox(ex_ionsim_tdphi, ex_analyt_tdphi, rtol=1e-5) # set Ω(t) to a step function L.λ = transitionwavelength(1, (("S1/2", -1 / 2), ("D5/2", -1 / 2)), T) phase!(L, 0) I = intensity_from_rabifrequency(1, Ω, 1, (("S1/2", -1 / 2), ("D5/2", -1 / 2)), T) L.I = t -> t < 1 ? 0 : I - h = hamiltonian(T, timescale=1e-6, lamb_dicke_order = 0) + h = hamiltonian(T, timescale=1e-6, lamb_dicke_order=0) tspan = 0:1e-3:3 tout, sol = timeevolution.schroedinger_dynamic( tspan, @@ -175,14 +175,14 @@ using Suppressor delayed_sin2(t) = t < 1 ? 0 : sin(2π * 1e-6 * Ω00 / 2 * (t - 1))^2 ex_ionsim_tdE = real.(expect(ionprojector(T, ("D5/2", -1 / 2)), sol)) ex_analyt_tdE = @.delayed_sin2(tout) - @test isapprox(ex_ionsim_tdE, ex_analyt_tdE, rtol = 1e-5) + @test isapprox(ex_ionsim_tdE, ex_analyt_tdE, rtol=1e-5) # check that δν(t) shifts sideband frequencies appropriately intensity!(L, I) tspan = 0:1e-3:30 mode.δν = t -> 20e3 L.Δ = mode.ν + 20e3 - h = hamiltonian(T, timescale=1e-6, rwa_cutoff = 1e5) + h = hamiltonian(T, timescale=1e-6, rwa_cutoff=1e5) tout, sol = timeevolution.schroedinger_dynamic( tspan, ionstate(T, [("S1/2", -1 / 2)]) ⊗ mode[1], @@ -191,7 +191,7 @@ using Suppressor η = lambdicke(mode, C, L) ex_ionsim_δν = real.(expect(ionprojector(T, ("D5/2", -1 / 2)), sol)) ex_analyt_δν = @.(sin(2π * sqrt(2) * (η / sqrt(1.02)) * Ω00 / 2 * 1e-6 * tout)^2) - @test isapprox(ex_ionsim_δν, ex_analyt_δν, rtol = 1e-1) + @test isapprox(ex_ionsim_δν, ex_analyt_δν, rtol=1e-1) # δB(t) ## add test @@ -203,11 +203,11 @@ using Suppressor L1 = Laser() L2 = Laser() chain = LinearChain( - ions = [C, C], - comfrequencies = (x = 3e6, y = 3e6, z = 1e6), - selectedmodes = (; z = [1]) + ions=[C, C], + comfrequencies=(x=3e6, y=3e6, z=1e6), + selectedmodes=(; z=[1]) ) - T = Chamber(iontrap = chain, B = 4e-4, Bhat = (x̂ + ẑ) / √2, lasers = [L1, L2]) + T = Chamber(iontrap=chain, B=4e-4, Bhat=(x̂ + ẑ) / √2, lasers=[L1, L2]) L1.λ = transitionwavelength(C, (("S1/2", -1 / 2), ("D5/2", -1 / 2)), T) L2.λ = transitionwavelength(C, (("S1/2", -1 / 2), ("D5/2", -1 / 2)), T) mode = T.iontrap.selectedmodes.z[1] @@ -229,14 +229,21 @@ using Suppressor intensity_from_rabifrequency!(1, Ω, 1, (("S1/2", -1 / 2), ("D5/2", -1 / 2)), T) intensity_from_rabifrequency!(2, Ω, 1, (("S1/2", -1 / 2), ("D5/2", -1 / 2)), T) ψi = ionstate(T, [("S1/2", -1 / 2), ("S1/2", -1 / 2)]) ⊗ mode[0] # initial state - h = hamiltonian(T, timescale=1e-6, rwa_cutoff = 5e5) + h = hamiltonian(T, timescale=1e-6, rwa_cutoff=5e5) tspan = 0:0.25:1000 tout, sol = timeevolution.schroedinger_dynamic(tspan, ψi, h) SS = expect(ionprojector(T, ("S1/2", -1 / 2), ("S1/2", -1 / 2)), sol) DD = expect(ionprojector(T, ("D5/2", -1 / 2), ("D5/2", -1 / 2)), sol) - ex = analytical.molmersorensen2ion(tspan, 1e-6Ω, 1e-6mode.ν, 1e-6mode.ν + 1e-6ϵ, η, 0) - @test isapprox(ex[1], SS, rtol = 1e-2) - @test isapprox(ex[2], DD, rtol = 1e-2) + ex = analytical.molmersorensen2ion( + tspan, + 1e-6Ω, + 1e-6mode.ν, + 1e-6mode.ν + 1e-6ϵ, + η, + 0 + ) + @test isapprox(ex[1], SS, rtol=1e-2) + @test isapprox(ex[2], DD, rtol=1e-2) # add checks for fast and hot gates end diff --git a/test/test_hamiltonians.jl b/test/test_hamiltonians.jl index 21ac0eb..23bb27e 100644 --- a/test/test_hamiltonians.jl +++ b/test/test_hamiltonians.jl @@ -12,7 +12,7 @@ function get_indices( ion_participation, vibrational_mode_dimension, two_modes; - cutoff = [Inf, Inf] + cutoff=[Inf, Inf] ) N = vibrational_mode_dimension[1] σ1 = [0 0; "a" 0] @@ -40,8 +40,7 @@ function get_indices( if k[i, j] == 0 continue end - mn = - map(x -> [parse(Int, x[2]), parse(Int, x[3])], split(k[i, j], ",")[1:(end - 1)]) + mn = map(x -> [parse(Int, x[2]), parse(Int, x[3])], split(k[i, j], ",")[1:(end-1)]) if sum([abs(x[1] - x[2]) > cutoff[l] for (l, x) in enumerate(mn)]) > 0 continue end @@ -98,7 +97,7 @@ end ijk = [(i1, i2), (j1, j2), (k1, k2)] (I, J) = IonSim._get_kron_indxs([(i1, i2), (j1, j2), (k1, k2)], [3, 3, 3]) @test [ - (parse(Int, i[2]), parse(Int, i[3])) for i in split(abc[I, J], ",")[1:(end - 1)] + (parse(Int, i[2]), parse(Int, i[3])) for i in split(abc[I, J], ",")[1:(end-1)] ] == ijk # ^this tests that abc[I, J] = a[i1,i2] * a[j1,j2] * c[k1,k2] @@ -124,11 +123,11 @@ end intensity!(L2, 2) phase!(L2, 2) chain = LinearChain( - ions = [C, C1], - comfrequencies = (x = 3e6, y = 3e6, z = 1e6), - selectedmodes = (; z = [1]) + ions=[C, C1], + comfrequencies=(x=3e6, y=3e6, z=1e6), + selectedmodes=(; z=[1]) ) - T = Chamber(iontrap = chain, lasers = [L1, L2]) + T = Chamber(iontrap=chain, lasers=[L1, L2]) Ωnmkj = IonSim._Ωmatrix(T, 1) t = rand(0:1e-3:100) @@ -184,9 +183,9 @@ end # _ηmatrix chain = LinearChain( - ions = [C, C1], - comfrequencies = (x = 2e6, y = 2e6, z = 1e6), - selectedmodes = (x = [1], y = [2], z = [1]) + ions=[C, C1], + comfrequencies=(x=2e6, y=2e6, z=1e6), + selectedmodes=(x=[1], y=[2], z=[1]) ) L1.k = (x̂ + ẑ) / √2 L2.k = (ŷ + ẑ) / √2 @@ -194,27 +193,27 @@ end L3.pointing = [(1, 1.0), (2, 1.0)] L3.λ = transitionwavelength(C, ("S1/2", "D5/2")) L3.k = ẑ - T = Chamber(iontrap = chain, lasers = [L1, L2, L3]) + T = Chamber(iontrap=chain, lasers=[L1, L2, L3]) η = IonSim._ηmatrix(T) # η[1,1,1] corresponds laser1, ion1, mode1 (x̂) and η[1,1,3] corresponds laser1, ion1, # mode3 (ẑ). They both have the same projection on L1.k, but the frequency of mode1 is # twice as large as mode2, so it's L-D factor should be √2 times smaller - @test abs(η[1, 1, end](1.0)) ≈ abs(η[1, 1, end - 2](1.0) / √2) + @test abs(η[1, 1, end](1.0)) ≈ abs(η[1, 1, end-2](1.0) / √2) # With this setup, the L-D factor should be the same for ion1 and ion2 @test η[1, 1, end](1.0) ≈ η[2, 1, end](1.0) # η[1, 2, 2] is the 1st ion, 2nd laser and y-stretch-mode. The L-D factor should be # opposite in sign, equal in magnitude to the 2nd ion, 2nd laser and y-stretch-mode - @test η[1, 2, end - 1](1.0) ≈ -η[2, 2, end - 1](1.0) + @test η[1, 2, end-1](1.0) ≈ -η[2, 2, end-1](1.0) # L3, which is in the ẑ direction should only have projection on zmode (mode3) @test η[1, 3, end] ≡ 0 - @test η[1, 3, end - 1] ≡ 0 - @test η[1, 3, end - 2](0.0) != 0 + @test η[1, 3, end-1] ≡ 0 + @test η[1, 3, end-2](0.0) != 0 # test construction of time-dep δν. If δν = 1e6*t, then after 3e-6 seconds (and since # ν=1e6), √(ν+δν(t)) = √2 * √(ν) = √2 * √(ν+δν(0)) frequency_fluctuation!(zmodes(chain)[1], t -> 1e6t) η = IonSim._ηmatrix(T) - @test η[1, 1, end - 2](0.0) ≈ 2 * η[1, 1, end - 2](3) + @test η[1, 1, end-2](0.0) ≈ 2 * η[1, 1, end-2](3) end @testset "hamiltonians -- setup functions" begin @@ -225,11 +224,11 @@ end L = Laser() L.λ = transitionwavelength(C, ("S1/2", "D5/2")) chain = LinearChain( - ions = [C], - comfrequencies = (x = 3e6, y = 3e6, z = 1e6), - selectedmodes = (; z = [1]) + ions=[C], + comfrequencies=(x=3e6, y=3e6, z=1e6), + selectedmodes=(; z=[1]) ) - T = Chamber(iontrap = chain, lasers = [L], δB = 0) + T = Chamber(iontrap=chain, lasers=[L], δB=0) global_B_indices, global_B_scales, bfunc = IonSim._setup_global_B_hamiltonian(T, 1) # T.δB = 0 -> global_B_indices and global_B_scales should be empty arrays @@ -255,11 +254,11 @@ end # setup system chain = LinearChain( - ions = [C], - comfrequencies = (x = 3e6, y = 3e6, z = 1e6), - selectedmodes = (; z = [1]) + ions=[C], + comfrequencies=(x=3e6, y=3e6, z=1e6), + selectedmodes=(; z=[1]) ) - T = Chamber(iontrap = chain, lasers = [L], δB = 0) + T = Chamber(iontrap=chain, lasers=[L], δB=0) # should return empty arrays if δν=0 δν_indices, δν_functions = IonSim._setup_δν_hamiltonian(T, 1) @test length(δν_indices) == 0 && length(δν_functions) == 0 @@ -273,23 +272,23 @@ end # add another mode with δν=0 and test output chain = LinearChain( - ions = [C], - comfrequencies = (x = 3e6, y = 3e6, z = 1e6), - selectedmodes = (y = [1], z = [1]) + ions=[C], + comfrequencies=(x=3e6, y=3e6, z=1e6), + selectedmodes=(y=[1], z=[1]) ) - T = Chamber(iontrap = chain, lasers = [L], δB = 0) + T = Chamber(iontrap=chain, lasers=[L], δB=0) modecutoff!(ymodes(T)[1], 3) modecutoff!(zmodes(T)[1], 3) frequency_fluctuation!(zmodes(T)[1], 1) δν_indices, δν_functions = IonSim._setup_δν_hamiltonian(T, 1) - indxs1 = [[8 * 4 * j + i for i in 1:(8 * 4)] for j in 1:3] + indxs1 = [[8 * 4 * j + i for i in 1:(8*4)] for j in 1:3] @test δν_indices[1] == indxs1 # test output when both modes have nonzero δν frequency_fluctuation!(zmodes(T)[1], 1) frequency_fluctuation!(ymodes(T)[1], 1) δν_indices, δν_functions = IonSim._setup_δν_hamiltonian(T, 1) - indxs = [[8 * j + i for i in 1:(8 * 8)] for j in 1:3] + indxs = [[8 * j + i for i in 1:(8*8)] for j in 1:3] @test δν_indices[1][1] == [vcat([[8 * (4(j - 1) + 1) + i for i in 1:8] for j in 1:4]...);] @test δν_indices[2] == indxs1 @@ -305,7 +304,7 @@ end @test δν_functions[2].(t) == 2π .* sin.(t) # _setup_fluctuation_hamiltonian - T = Chamber(iontrap = chain, lasers = [L], δB = 1) + T = Chamber(iontrap=chain, lasers=[L], δB=1) frequency_fluctuation!(ymodes(T)[1], cos) frequency_fluctuation!(zmodes(T)[1], sin) δν_indices, δν_functions = IonSim._setup_δν_hamiltonian(T, 1) @@ -326,16 +325,11 @@ end L = Laser() L.λ = transitionwavelength(C, ("S1/2", "D5/2")) chain = LinearChain( - ions = [C, C], - comfrequencies = (x = 3e6, y = 3e6, z = 1e6), - selectedmodes = (; z = [1]) - ) - T = Chamber( - iontrap = chain, - B = 4e-4, - Bhat = (x̂ + ŷ + ẑ) / √3, - lasers = [L] + ions=[C, C], + comfrequencies=(x=3e6, y=3e6, z=1e6), + selectedmodes=(; z=[1]) ) + T = Chamber(iontrap=chain, B=4e-4, Bhat=(x̂ + ŷ + ẑ) / √3, lasers=[L]) mode = T.iontrap.selectedmodes.z[1] modecutoff!(mode, rand(1:8)) N = mode.N + 1 @@ -367,16 +361,11 @@ end L1 = Laser() L1.λ = transitionwavelength(C1, ("S1/2", "D5/2")) chain1 = LinearChain( - ions = [C1, C1], - comfrequencies = (x = 3e6, y = 3e6, z = 1e6), - selectedmodes = (; z = [1, 2]) - ) - T1 = Chamber( - iontrap = chain1, - B = 4e-4, - Bhat = (x̂ + ŷ + ẑ) / √3, - lasers = [L1] + ions=[C1, C1], + comfrequencies=(x=3e6, y=3e6, z=1e6), + selectedmodes=(; z=[1, 2]) ) + T1 = Chamber(iontrap=chain1, B=4e-4, Bhat=(x̂ + ŷ + ẑ) / √3, lasers=[L1]) mode1 = T1.iontrap.selectedmodes.z[1] mode2 = T1.iontrap.selectedmodes.z[2] modecutoff!(mode1, N - 1) @@ -392,7 +381,7 @@ end ## now set non-trivial lamb_dicke_order lamb_dicke_order = [rand(1:N), rand(1:M)] - ridxs, cidxs = get_indices(12, [N, M], true, cutoff = reverse(lamb_dicke_order)) + ridxs, cidxs = get_indices(12, [N, M], true, cutoff=reverse(lamb_dicke_order)) _, r, c = IonSim._setup_base_hamiltonian( T1, 1e-6, @@ -424,27 +413,28 @@ end L = Laser() L.pointing = [(1, 1.0), (2, 1.0)] chain = LinearChain( - ions = [C_a, C_b], - comfrequencies = (x = 3e6, y = 3e6, z = 1e6), - selectedmodes = (; z = [1, 2]) - ) - T = Chamber( - iontrap = chain, - B = 4e-4, - Bhat = (x̂ + ŷ + ẑ) / √3, - lasers = [L] + ions=[C_a, C_b], + comfrequencies=(x=3e6, y=3e6, z=1e6), + selectedmodes=(; z=[1, 2]) ) + T = Chamber(iontrap=chain, B=4e-4, Bhat=(x̂ + ŷ + ẑ) / √3, lasers=[L]) L.λ = transitionwavelength(C_a, (("S1/2", -1 / 2), ("D5/2", -1 / 2)), T) mode1 = T.iontrap.selectedmodes.z[1] mode2 = T.iontrap.selectedmodes.z[2] - Δ = round(randn(), digits = 5) * 1e5 # TODO: this begins to fail at below 1 Hz! + Δ = round(randn(), digits=5) * 1e5 # TODO: this begins to fail at below 1 Hz! L.Δ = Δ ϕ = randn() phase!(L, ϕ) modecutoff!(mode1, 10) modecutoff!(mode2, 9) Ω = abs(randn()) # absolute value ensures QO hamiltonian will match, since IonSim will always have positive Rabi frequency owing to using laser intensity rather than amplitude - intensity_from_rabifrequency!(1, Ω * 1e6, 1, (("S1/2", -1 / 2), ("D5/2", -1 / 2)), T) + intensity_from_rabifrequency!( + 1, + Ω * 1e6, + 1, + (("S1/2", -1 / 2), ("D5/2", -1 / 2)), + T + ) # Case 1a: full hamiltonian (w conj_repeated_indices); controlled by not specifying rwa_cutoff @@ -462,79 +452,75 @@ end ηb2 = lambdicke(mode2, C_b, L) # displacement operators for COM and stretch modes - mode_op1(t; η) = displace(mode1, im * η * exp(im * 2π * t), method = "truncated") - mode_op2(t; η) = - displace(mode2, im * η * exp(im * 2π * √3 * t), method = "truncated") + mode_op1(t; η) = displace(mode1, im * η * exp(im * 2π * t), method="truncated") + mode_op2(t; η) = displace(mode2, im * η * exp(im * 2π * √3 * t), method="truncated") Hp(t) = ( - ion_op(t) ⊗ one(C_b) ⊗ mode_op1(t, η = ηa1) ⊗ mode_op2(t, η = ηa2) + - one(C_a) ⊗ ion_op(t) ⊗ mode_op1(t, η = ηb1) ⊗ mode_op2(t, η = ηb2) + ion_op(t) ⊗ one(C_b) ⊗ mode_op1(t, η=ηa1) ⊗ mode_op2(t, η=ηa2) + + one(C_a) ⊗ ion_op(t) ⊗ mode_op1(t, η=ηb1) ⊗ mode_op2(t, η=ηb2) ) qoH(t) = Hp(t) + dagger(Hp(t)) tp = abs(51randn()) # define hamiltonians.jl - H = hamiltonian(T, timescale=1e-6, lamb_dicke_order = 101) - H1 = hamiltonian(T, timescale=1e-6, lamb_dicke_order = 101, time_dependent_eta = true) + H = hamiltonian(T, timescale=1e-6, lamb_dicke_order=101) + H1 = hamiltonian(T, timescale=1e-6, lamb_dicke_order=101, time_dependent_eta=true) # test similarity at a random time input @test norm(qoH(tp).data - H(tp, 0).data) < 1e-4 @test norm(qoH(tp).data - H1(tp, 0).data) < 1e-4 # Case 1b: analytic solution (as opposed to truncated solution) - mode_op1(t; η) = displace(mode1, im * η * exp(im * 2π * t), method = "analytic") - mode_op2(t; η) = - displace(mode2, im * η * exp(im * 2π * √3 * t), method = "analytic") + mode_op1(t; η) = displace(mode1, im * η * exp(im * 2π * t), method="analytic") + mode_op2(t; η) = displace(mode2, im * η * exp(im * 2π * √3 * t), method="analytic") H1 = hamiltonian( T, timescale=1e-6, - lamb_dicke_order = 101, - displacement = "analytic", - time_dependent_eta = false + lamb_dicke_order=101, + displacement="analytic", + time_dependent_eta=false ) H = hamiltonian( T, timescale=1e-6, - lamb_dicke_order = 101, - displacement = "analytic", - time_dependent_eta = true + lamb_dicke_order=101, + displacement="analytic", + time_dependent_eta=true ) @test norm(qoH(tp).data - H(tp, 0).data) < 1e-4 @test H1(tp, 0).data ≈ H(tp, 0).data # Case 2a: full hamiltonian (w/o conj_repeated_indices) - H = hamiltonian(T, timescale=1e-6, lamb_dicke_order = 101, rwa_cutoff = 1e10) + H = hamiltonian(T, timescale=1e-6, lamb_dicke_order=101, rwa_cutoff=1e10) H1 = hamiltonian( T, timescale=1e-6, - lamb_dicke_order = 101, - time_dependent_eta = true, - rwa_cutoff = 1e10 + lamb_dicke_order=101, + time_dependent_eta=true, + rwa_cutoff=1e10 ) - mode_op1(t; η) = displace(mode1, im * η * exp(im * 2π * t), method = "truncated") - mode_op2(t; η) = - displace(mode2, im * η * exp(im * 2π * √3 * t), method = "truncated") + mode_op1(t; η) = displace(mode1, im * η * exp(im * 2π * t), method="truncated") + mode_op2(t; η) = displace(mode2, im * η * exp(im * 2π * √3 * t), method="truncated") @test norm(qoH(tp).data - H(tp, 0).data) < 1e-4 @test norm(qoH(tp).data - H1(tp, 0).data) < 1e-4 # Case 2b: Like 2a, but with analytic solution - mode_op1(t; η) = displace(mode1, im * η * exp(im * 2π * t), method = "analytic") - mode_op2(t; η) = - displace(mode2, im * η * exp(im * 2π * √3 * t), method = "analytic") + mode_op1(t; η) = displace(mode1, im * η * exp(im * 2π * t), method="analytic") + mode_op2(t; η) = displace(mode2, im * η * exp(im * 2π * √3 * t), method="analytic") H1 = hamiltonian( T, timescale=1e-6, - lamb_dicke_order = 101, - displacement = "analytic", - time_dependent_eta = false, - rwa_cutoff = 1e10 + lamb_dicke_order=101, + displacement="analytic", + time_dependent_eta=false, + rwa_cutoff=1e10 ) H = hamiltonian( T, timescale=1e-6, - lamb_dicke_order = 101, - displacement = "analytic", - time_dependent_eta = true, - rwa_cutoff = 1e10 + lamb_dicke_order=101, + displacement="analytic", + time_dependent_eta=true, + rwa_cutoff=1e10 ) @test norm(qoH(tp).data - H(tp, 0).data) < 1e-4 @test H1(tp, 0).data ≈ H(tp, 0).data @@ -549,28 +535,28 @@ end Hp(t) = ion_op(t) ⊗ one(C_b) ⊗ mode_op_a + one(C_a) ⊗ ion_op(t) ⊗ mode_op_b qoH(t) = Hp(t) + dagger(Hp(t)) - H = hamiltonian(T, timescale=1e-6, lamb_dicke_order = 0, rwa_cutoff = Inf) + H = hamiltonian(T, timescale=1e-6, lamb_dicke_order=0, rwa_cutoff=Inf) H1 = hamiltonian( T, timescale=1e-6, - lamb_dicke_order = 0, - rwa_cutoff = Inf, - time_dependent_eta = true + lamb_dicke_order=0, + rwa_cutoff=Inf, + time_dependent_eta=true ) H2 = hamiltonian( T, timescale=1e-6, - lamb_dicke_order = 0, - rwa_cutoff = Inf, - displacement = "analytic" + lamb_dicke_order=0, + rwa_cutoff=Inf, + displacement="analytic" ) H3 = hamiltonian( T, timescale=1e-6, - lamb_dicke_order = 0, - rwa_cutoff = Inf, - displacement = "analytic", - time_dependent_eta = true + lamb_dicke_order=0, + rwa_cutoff=Inf, + displacement="analytic", + time_dependent_eta=true ) # only considering first order corrections to carrier (propto η^2) so this won't be perfect @test norm((qoH(tp) - H(tp, 0)).data) < 2.5 @@ -579,28 +565,28 @@ end @test norm((qoH(tp) - H3(tp, 0)).data) < 2.5 # Case 4: Try a normal example with a intermediate rwa_cutoff value - H = hamiltonian(T, timescale=1e-6, lamb_dicke_order = 30, rwa_cutoff = 3e5) + H = hamiltonian(T, timescale=1e-6, lamb_dicke_order=30, rwa_cutoff=3e5) H1 = hamiltonian( T, timescale=1e-6, - lamb_dicke_order = 30, - rwa_cutoff = 3e5, - time_dependent_eta = true + lamb_dicke_order=30, + rwa_cutoff=3e5, + time_dependent_eta=true ) H2 = hamiltonian( T, timescale=1e-6, - lamb_dicke_order = 30, - rwa_cutoff = 3e5, - displacement = "analytic" + lamb_dicke_order=30, + rwa_cutoff=3e5, + displacement="analytic" ) H3 = hamiltonian( T, timescale=1e-6, - lamb_dicke_order = 30, - rwa_cutoff = 3e5, - displacement = "analytic", - time_dependent_eta = true + lamb_dicke_order=30, + rwa_cutoff=3e5, + displacement="analytic", + time_dependent_eta=true ) # only considering first order corrections to carrier (propto η^2) so this won't be perfect @test norm((qoH(tp) - H(tp, 0)).data) < 2 @@ -612,8 +598,8 @@ end @test_throws AssertionError hamiltonian( T, timescale=1e-6, - lamb_dicke_order = [1, 2, 3], - rwa_cutoff = 3e5 + lamb_dicke_order=[1, 2, 3], + rwa_cutoff=3e5 ) end end # end suppress diff --git a/test/test_ion_traps.jl b/test/test_ion_traps.jl index 6ae2ed3..94feb8b 100644 --- a/test/test_ion_traps.jl +++ b/test/test_ion_traps.jl @@ -6,9 +6,9 @@ using Suppressor @testset "ion_traps -- LinearChain" begin C = Ca40() lc = LinearChain( - ions = [C, C, C, C], - comfrequencies = (x = 5, y = 5, z = 1), - selectedmodes = (y = [1], z = [4]) + ions=[C, C, C, C], + comfrequencies=(x=5, y=5, z=1), + selectedmodes=(y=[1], z=[4]) ) @test ions(lc) == lc.ions # test modes, which should return an array of the selected @@ -29,15 +29,15 @@ using Suppressor # should get warning if same ion is input multiple times to ion kwarg warning = "Some ions point to the same thing. Making copies." @test_logs (:warn, warning) LinearChain( - ions = [C, C], - comfrequencies = (x = 4, y = 4, z = 1), - selectedmodes = (x = [], y = [], z = [1, 2]) + ions=[C, C], + comfrequencies=(x=4, y=4, z=1), + selectedmodes=(x=[], y=[], z=[1, 2]) ) # and copies should be made of the repeated ions, so that they are no longer the same chain1 = LinearChain( - ions = [C, C], - comfrequencies = (x = 4, y = 4, z = 1), - selectedmodes = (x = [], y = [], z = [1, 2]) + ions=[C, C], + comfrequencies=(x=4, y=4, z=1), + selectedmodes=(x=[], y=[], z=[1, 2]) ) @test !(chain1.ions[1] ≡ chain1.ions[2]) @@ -51,21 +51,21 @@ using Suppressor # known values posL = [-2.8708, -2.10003, -1.4504, -0.85378, -0.2821] pos = [posL; -1 .* reverse(posL)] - @test any(isapprox.(linear_equilibrium_positions(10), pos, rtol = 1e-4)) + @test any(isapprox.(linear_equilibrium_positions(10), pos, rtol=1e-4)) # and test calculation of characterstic length scale for linear chain, equal mass C = Ca40() @test characteristic_length_scale(mass(C), 1e6) ≈ 4.449042804354206e-6 # and do the same for Anm, which computes the normal modes - @test_throws AssertionError Anm(2, (x = 0.5, y = 0.5, z = 1), (x = 1, y = 0, z = 0)) + @test_throws AssertionError Anm(2, (x=0.5, y=0.5, z=1), (x=1, y=0, z=0)) cst = [-0.2132, 0.6742, -0.6742, 0.2132] - freq, mode = Anm(4, (x = 2, y = 2, z = 1), (x = 0, y = 0, z = 1))[end] + freq, mode = Anm(4, (x=2, y=2, z=1), (x=0, y=0, z=1))[end] @test freq ≈ √9.308 rtol = 1e-4 if mode[1] > 0 cst = -1 .* cst end - @test any(isapprox.(mode, cst, rtol = 1e-4)) + @test any(isapprox.(mode, cst, rtol=1e-4)) # test '_sparsify', which should drop small values from an array x = [1e-6, 1e-5] diff --git a/test/test_ions.jl b/test/test_ions.jl index 0cc36d9..244b504 100644 --- a/test/test_ions.jl +++ b/test/test_ions.jl @@ -62,7 +62,7 @@ using Suppressor @test quantumnumbers(C, ("D5/2", 1 / 2)).m == 1 // 2 @test energy(C, "S1/2") == energy(C, "S") - @test energy(C, "S1/2") != energy(C, "S", B = 1e-4) + @test energy(C, "S1/2") != energy(C, "S", B=1e-4) @test transitionwavelength(C, ("S", "D")) ≈ 7.29147e-7 end diff --git a/test/test_lasers.jl b/test/test_lasers.jl index e449079..b59ff98 100644 --- a/test/test_lasers.jl +++ b/test/test_lasers.jl @@ -3,30 +3,30 @@ using Suppressor @testset "lasers -- Laser" begin # should not be able to set unnormalized k-vectors or polarizations - @test_throws AssertionError Laser(ϵ = (x = 2, y = 0, z = 0)) - @test_throws AssertionError Laser(k = (x = 2, y = 0, z = 0)) + @test_throws AssertionError Laser(ϵ=(x=2, y=0, z=0)) + @test_throws AssertionError Laser(k=(x=2, y=0, z=0)) # The below test is commented out because we have temporarily removed the non-orthogonal warning for convenience (see lasers.jl) # should not be able to initialize with non-orthogonal ϵ̂, k̂ # @test_throws AssertionError Laser(ϵ=(x=1, y=0, z=0), k=(x=1, y=0, z=0)) # pointing cannot be overspecified - @test_throws AssertionError Laser(pointing = [(1, 0.5), (1, 0.5)]) + @test_throws AssertionError Laser(pointing=[(1, 0.5), (1, 0.5)]) # pointing magnitude must be within [0,1] - @test_throws AssertionError Laser(pointing = [(1, 2.0)]) + @test_throws AssertionError Laser(pointing=[(1, 2.0)]) # should be able to set L.I, L.ϕ to a constant Real value t = 0:1:100 - L = Laser(I = 1, ϕ = 1) + L = Laser(I=1, ϕ=1) ones = [1 for _ in t] # should be able to set L.I, L.ϕ to a function of time @test L.I.(t) == ones @test L.ϕ.(t) == ones - L = Laser(I = sin, ϕ = sin) + L = Laser(I=sin, ϕ=sin) @test L.I.(t) == sin.(t) @test L.ϕ.(t) == sin.(t) # test that normalization is enforced for altered L.ϵ/L.k - @test_throws AssertionError polarization!(L, (x = 2, y = 0, z = 0)) - @test_throws AssertionError wavevector!(L, (x = 2, y = 0, z = 0)) + @test_throws AssertionError polarization!(L, (x=2, y=0, z=0)) + @test_throws AssertionError wavevector!(L, (x=2, y=0, z=0)) # test that pointing constraints are enforced for altered values @test_throws AssertionError pointing!(L, [(1, 0.5), (1, 0.5)]) @test_throws AssertionError pointing!(L, [(1, 2.0)]) diff --git a/test/test_operators.jl b/test/test_operators.jl index 5d359e1..69cc1cd 100644 --- a/test/test_operators.jl +++ b/test/test_operators.jl @@ -8,11 +8,11 @@ using Suppressor # setup system C = Ca40([("S1/2", -1 // 2), ("D5/2", -1 // 2)]) chain = LinearChain( - ions = [C, C], - comfrequencies = (x = 2, y = 2, z = 1), - selectedmodes = (x = [1], y = [], z = [1]) + ions=[C, C], + comfrequencies=(x=2, y=2, z=1), + selectedmodes=(x=[1], y=[], z=[1]) ) - T = Chamber(iontrap = chain) + T = Chamber(iontrap=chain) allmodes = modes(chain) @testset "operators -- VibrationalMode operators" begin @@ -31,14 +31,17 @@ using Suppressor modecutoff!(allmodes[1], 200) i = rand(1:5) j = rand(1:5) - @test displace(allmodes[1], α, method = "analytic").data[i, j] ≈ + @test displace(allmodes[1], α, method="analytic").data[i, j] ≈ qo.displace(fb2, α).data[i, j] atol = 1e-6 # test that mean excitation of thermalstate is as expected modecutoff!(allmodes[1], 500) n̄ = abs(2randn()) @test expect(number(allmodes[1]), thermalstate(allmodes[1], n̄)) ≈ n̄ - @test expect(number(allmodes[1]), thermalstate(allmodes[1], n̄, method = "analytic")) ≈ n̄ + @test expect( + number(allmodes[1]), + thermalstate(allmodes[1], n̄, method="analytic") + ) ≈ n̄ # test coherentstate matches QO results α = 10 * (randn() + im * randn()) @@ -48,11 +51,12 @@ using Suppressor N = 500 modecutoff!(allmodes[1], N) n̄ = rand(0:1e-6:10) - @test coherentthermalstate(allmodes[1], n̄, 0, method = "analytic").data ≈ + @test coherentthermalstate(allmodes[1], n̄, 0, method="analytic").data ≈ thermalstate(allmodes[1], n̄).data - @test coherentthermalstate(allmodes[1], 0, α, method = "analytic").data ≈ + @test coherentthermalstate(allmodes[1], 0, α, method="analytic").data ≈ dm(coherentstate(allmodes[1], α)).data rtol = 1e-3 * N^2 - @test coherentthermalstate(allmodes[1], n̄, 0).data ≈ thermalstate(allmodes[1], n̄).data + @test coherentthermalstate(allmodes[1], n̄, 0).data ≈ + thermalstate(allmodes[1], n̄).data @test coherentthermalstate(allmodes[1], 0, α).data ≈ dm(coherentstate(allmodes[1], α)).data rtol = 1e-3 * N^2 @@ -86,14 +90,14 @@ using Suppressor @test sigma(C, ("S1/2", -1 // 2)).data == sigma(C, 1).data == ComplexF64[1 0; 0 0] # test ionprojector for IonTrap input - ψ = ionprojector(chain, ("S1/2", -1 // 2), ("D5/2", -1 // 2), only_ions = true) + ψ = ionprojector(chain, ("S1/2", -1 // 2), ("D5/2", -1 // 2), only_ions=true) @test ψ.data == kron( ComplexF64[0; 1] * ComplexF64[0; 1]', ComplexF64[1; 0] * ComplexF64[1; 0]' ) @test ionprojector(chain, ("S1/2", -1 // 2), ("D5/2", -1 // 2)) == ψ ⊗ one(allmodes[1]) ⊗ one(allmodes[2]) - @test ψ == ionprojector(T, ("S1/2", -1 // 2), ("D5/2", -1 // 2), only_ions = true) + @test ψ == ionprojector(T, ("S1/2", -1 // 2), ("D5/2", -1 // 2), only_ions=true) end @testset "operators -- internal functions" begin diff --git a/test/test_vibrational_modes.jl b/test/test_vibrational_modes.jl index 4052fbe..2466772 100644 --- a/test/test_vibrational_modes.jl +++ b/test/test_vibrational_modes.jl @@ -3,7 +3,7 @@ using Suppressor @testset "vibrational_modes -- VibrationalMode" begin # setup system - vm = VibrationalMode(1, [1, 1], δν = 1) + vm = VibrationalMode(1, [1, 1], δν=1) # make sure δν can appropriately handle constant or function t = 0:1:100 @@ -13,11 +13,11 @@ using Suppressor frequency_fluctuation!(vm, sin) @test vm.δν.(t) == sin.(t) @test !vm._cnst_δν - vm = VibrationalMode(1, [1, 1], δν = t -> t) + vm = VibrationalMode(1, [1, 1], δν=t -> t) @test !vm._cnst_δν # test == - vm1 = VibrationalMode(1, [1, 1], δν = 1) + vm1 = VibrationalMode(1, [1, 1], δν=1) @test vm == vm1 # test that updating cutoff also appropriately updates shape field