diff --git a/CHANGELOG b/CHANGELOG
new file mode 100644
index 0000000..f938118
--- /dev/null
+++ b/CHANGELOG
@@ -0,0 +1,3 @@
+- Created an abstract noise operation class `AbstractNoiseBellOp`
+- Implemented `PauliZOp` for phase damping errors and `MixedStateOp` for amplitude damping errors
+- Added tests for `PauliZOp` and `MixedStateOp` for comparison to exact density matrix computation
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index 95f251f..6f0e196 100644
--- a/Project.toml
+++ b/Project.toml
@@ -17,4 +17,4 @@ BPGatesQuantikzExt = "Quantikz"
 Quantikz = "1.3.1"
 QuantumClifford = "0.8, 0.9"
 Random = "1"
-julia = "1.9"
\ No newline at end of file
+julia = "1.9"
diff --git a/src/BPGates.jl b/src/BPGates.jl
index facaaec..9213c6a 100644
--- a/src/BPGates.jl
+++ b/src/BPGates.jl
@@ -9,7 +9,8 @@ export BellState,
     BellSinglePermutation, BellDoublePermutation, BellPauliPermutation,
     BellMeasure, bellmeasure!,
     BellGate, CNOTPerm, GoodSingleQubitPerm,
-    PauliNoiseOp, PauliNoiseBellGate, NoisyBellMeasure, NoisyBellMeasureNoisyReset
+    PauliNoiseOp, PauliZOp, PauliNoiseBellGate, NoisyBellMeasure, NoisyBellMeasureNoisyReset, AbstractNoiseBellOp,
+    MixedStateOp, AsymmDepOp
 
 function int_to_bit(int,digits)
     int = int - 1 # -1 so that we use julia indexing conventions
@@ -168,6 +169,7 @@ struct BellMeasure <: QuantumClifford.AbstractMeasurement
     sidx::Int
     function BellMeasure(p,s)
         1 <= p <= 3 || throw(ArgumentError("The basis measurement index needs to be between 1 and 3"))
+        # why is there no check here for the measured bell pair? ensuring that the index is >= 1?
         new(p,s)
     end
 end
@@ -214,9 +216,13 @@ end
 
 """The permutations realized by [`BellPauliPermutation`](@ref)."""
 const pauli_perm_tuple = (
+    # 3 is actually 2; 2 is actually 3
     (1, 2, 3, 4),
+    # Z gate?
     (3, 4, 1, 2), ## X flip
+    # X gate?
     (2, 1, 4, 3), ## Z flip
+    # bit and phase flip
     (4, 3, 2, 1)  ## Y flip
 )
 
@@ -226,9 +232,16 @@ const pauli_perm_qc = (sId1,sX,sZ,sY)
 function QuantumClifford.apply!(state::BellState, op::BellPauliPermutation) # TODO abstract away the permutation application as it is used by other gates too
     phase = state.phases
     @inbounds phase_idx = bit_to_int(phase[op.sidx*2-1],phase[op.sidx*2])
+    # println("QuantumClifford.apply!, BellPauliPermutation, op.sidx: $(op.sidx)")
+    # println("QuantumClifford.apply!, BellPauliPermutation, phase_idx: $phase_idx")
+    # println("QuantumClifford.apply!, BellPauliPermutation, op.pidx: $(op.pidx)")
     @inbounds perm = pauli_perm_tuple[op.pidx]
     @inbounds permuted_idx = perm[phase_idx]
+    # println("QuantumClifford.apply!, BellPauliPermutation, permuted_idx: $permuted_idx")
+    # this makes no sense. Make bit 2 go before bit1.
     bit1, bit2 = int_to_bit(permuted_idx,Val(2))
+    # println("QuantumClifford.apply!, BellPauliPermutation, bit1: $bit1")
+    # println("QuantumClifford.apply!, BellPauliPermutation, bit2: $bit2")
     @inbounds phase[op.sidx*2-1] = bit1
     @inbounds phase[op.sidx*2] = bit2
     return state
@@ -321,6 +334,7 @@ julia> bellmeasure!(BellState([0,1,1,1]), BellMeasure(2,1))
 function bellmeasure!(state::BellState, op::BellMeasure) # TODO document which index corresponds to which measurement
     phase = state.phases
     result = measure_tuple[bit_to_int(phase[op.sidx*2-1],phase[op.sidx*2])][op.midx]
+    # TODO: introduce latency here?
     phase[op.sidx*2-1:op.sidx*2] .= 0 # reset the measured pair to 00
     return state, result
 end
@@ -418,6 +432,8 @@ const good_perm_qc = ( # From the appendix of Optimized Entanglement Purificatio
     (h*ph,p*hp)
 )
 
+# Does this only modify the second bell pair? No phase kickback?-- interesting.
+# phase kickback here contributes to a global phase
 const cnot_perm = (1, 2, 11, 12, 6, 5, 16, 15, 9, 10, 3, 4, 14, 13, 8, 7)
 
 function QuantumClifford.apply!(state::BellState, g::CNOTPerm) # TODO abstract away the permutation application as it is used by other gates too
@@ -471,7 +487,10 @@ end
 
 """A wrapper for `BellGate` that implements Pauli noise in addition to the gate."""
 struct PauliNoiseBellGate{G} <: BellOp where {G<:BellOp} # TODO make it work with the QuantumClifford noise ops
+    # to do the above TODO, probably have to make it match some interface?
     g::G
+    # probability that the state changes in these bases after the noise is applied
+    # probability that the qubit at the specified index changes to an X gate.
     px::Float64
     py::Float64
     pz::Float64
@@ -518,6 +537,112 @@ function QuantumClifford.apply!(state::BellState, g::PauliNoiseOp)
     return state
 end
 
+# Type for abstracting away noisy operations on bell states.
+abstract type AbstractNoiseBellOp <: BellOp end
+
+"""
+PauliZOp(idx, pz) causes qubit at idx `idx` to have a Z flip with probabilty `pz`.
+"""
+struct PauliZOp <: AbstractNoiseBellOp
+    idx::Int
+    # assert that 0 <= pz <= 1? not strictly necessary. it wouldn't make sense, but it
+    # also wouldn't break the code. natural floor and ceiling here.
+    pz::Float64
+end
+
+function QuantumClifford.apply!(state::BellState, g::PauliZOp)
+    i = g.idx
+    r = rand()
+    if r<g.pz
+        # I believe this actually is a Z flip.
+        apply!(state, BellPauliPermutation(2, i))
+    end
+    return state
+end
+
+# TODO: continue from here
+# TODO: assert that the values add to 1
+# This is NOT at all 'extensible' (though, it is 'fast')
+# TODO: Honestly this doesn't make that much sense to me. Think about why
+# this works.
+# Defines the probabilities of transitioning from one bell state to another
+function get_mixed_transition_probs(λ, cur_state_idx)
+    # 0 <= λ <= 1 (λ = 1 - e ^(-t / T1))
+    mixed_state_tuple = (
+        (0.5 * λ^2 - λ + 1, 0.5 * λ * (1 - λ), 0.5 * λ^2, 0.5 * λ * (1 - λ)),
+        (0.5 * λ, 1 - λ, 0.5 * λ, 0.0),
+        (0.5 * λ^2, 0.5 * λ * (1 - λ), 0.5 * λ^2 - λ + 1, 0.5 * λ * (1 - λ)),
+        (0.5 * λ, 0.0, 0.5 * λ, 1 - λ)
+    )
+    return mixed_state_tuple[cur_state_idx]
+end
+
+"""
+MixedStateOp(idx, lambda) causes Bell state at idx `idx` to change to another
+Bell state determined by `mixed_state_tuple` and `lambda`.
+"""
+struct MixedStateOp <: AbstractNoiseBellOp
+    idx::Int
+    lambda::Float64
+end
+
+function select_rand_state(transition_probs)
+    # TODO: hardcoded for speed, BUT change this to a macro for readability.
+    rand_prob_val = rand()
+    new_phase_idx = 0
+    if rand_prob_val < transition_probs[1]
+        new_phase_idx = 1
+    elseif rand_prob_val < transition_probs[1] + transition_probs[2]
+        new_phase_idx = 2
+    elseif rand_prob_val < transition_probs[1] + transition_probs[2] + transition_probs[3]
+        new_phase_idx = 3
+    else
+        new_phase_idx = 4
+    end
+    return new_phase_idx
+end
+
+function QuantumClifford.apply!(state::BellState, g::MixedStateOp)
+    state_idx = g.idx
+    phase = state.phases
+    @inbounds phase_idx = bit_to_int(phase[state_idx * 2 - 1],phase[state_idx * 2])
+    transition_probs = get_mixed_transition_probs(g.lambda, phase_idx)
+    new_phase_idx = select_rand_state(transition_probs)
+    # this should never happen; can remove for speed
+    @assert new_phase_idx != 0
+    state_bit1, state_bit2 = int_to_bit(new_phase_idx, Val(2))
+    @inbounds phase[state_idx * 2 - 1] = state_bit1
+    @inbounds phase[state_idx * 2] = state_bit2
+    return state
+end
+
+struct AsymmDepOp <: AbstractNoiseBellOp
+    idx::Int
+    px::Float64
+    py::Float64
+    pz::Float64
+end
+
+function QuantumClifford.apply!(state::BellState, g::AsymmDepOp)
+    i = g.idx
+    # TODO repetition with ...NoisyReset and PauliNoise...
+    # ^ ?
+    r = rand()
+    # 2: Z, 3: X, 4: Y
+    # println("QuantumClifford.apply!, AsymmDepOp")
+    # print("QuantumClifford.apply!, AsymmDepOp, px_bell: $px_bell")
+    # print("QuantumClifford.apply!, AsymmDepOp, py_bell: $py_bell")
+    # print("QuantumClifford.apply!, AsymmDepOp, pz_bell: $pz_bell")
+    if r<g.px
+        apply!(state, BellPauliPermutation(3, i))
+    elseif r<g.px+g.pz
+        apply!(state, BellPauliPermutation(2, i))
+    elseif r<g.px+g.pz+g.py
+        apply!(state, BellPauliPermutation(4, i))
+    end
+    return state
+end
+
 """A wrapper for [`BellMeasure`](@ref) that implements measurement noise."""
 struct NoisyBellMeasure <: BellOp # TODO make it work with the QuantumClifford noise ops
     m::BellMeasure
@@ -527,7 +652,9 @@ end
 """A wrapper for [`BellMeasure`](@ref) that implements measurement noise and Pauli noise after the reset."""
 struct NoisyBellMeasureNoisyReset <: BellOp # TODO make it work with the QuantumClifford noise ops
     m::BellMeasure
+    # this p is the probability that the incorrect measurement is returned?
     p::Float64
+    # probability of resetting to a different Bell state after measurement
     px::Float64
     py::Float64
     pz::Float64
@@ -538,8 +665,10 @@ function QuantumClifford.applywstatus!(state::BellState, op::NoisyBellMeasure)
     state, result⊻(rand()<op.p) ? continue_stat : failure_stat
 end
 
+# This function is like a 'continuation' of the above function
 function QuantumClifford.applywstatus!(state::BellState, op::NoisyBellMeasureNoisyReset)
     state, result = bellmeasure!(state, op.m)
+    # pretty high probability of flipping to a correct measurement?
     cont = result⊻(rand()<op.p)
     cont && apply!(state, PauliNoiseOp(op.m.sidx,op.px,op.py,op.pz))
     state, cont ? continue_stat : failure_stat
diff --git a/test/Project.toml b/test/Project.toml
index ccf8fc1..c0e7bd2 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -6,6 +6,7 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 Quantikz = "b0d11df0-eea3-4d79-b4a5-421488cbf74b"
 QuantumClifford = "0525e862-1e90-11e9-3e4d-1b39d7109de1"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
diff --git a/test/runtests.jl b/test/runtests.jl
index 5eb00ca..80e8f20 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -16,17 +16,26 @@ end
 macro doset(descr)
     quote
         if doset($descr)
-            @safetestset $descr begin include("test_"*$descr*".jl") end
+            @safetestset $descr begin 
+            if $descr == "mixedstateop" || $descr == "paulizop"
+                include("test_kraus_mem_noise_helpers.jl")
+            end
+            include("test_"*$descr*".jl") end
         end
     end
 end
 
 println("Starting tests with $(Threads.nthreads()) threads out of `Sys.CPU_THREADS = $(Sys.CPU_THREADS)`...")
 
+# include("test_kraus_mem_noise_helpers.jl")
+
 @doset "quantumclifford"
 @doset "quantikz"
 get(ENV,"JET_TEST","")=="true" && @doset "jet"
 @doset "doctests"
+@doset "bpgates"
+@doset "mixedstateop"
+@doset "paulizop"
 
 using Aqua
 doset("aqua") && begin
diff --git a/test/test_bpgates.jl b/test/test_bpgates.jl
new file mode 100644
index 0000000..2b834dc
--- /dev/null
+++ b/test/test_bpgates.jl
@@ -0,0 +1,108 @@
+using Test
+using Logging
+
+using BPGates
+
+using QuantumClifford
+
+function test_pauliz_constructor()
+
+    test_pauliz = PauliZOp(1, 0.5)
+    
+    @test test_pauliz isa PauliZOp
+
+    @test test_pauliz.idx == 1
+
+    @test test_pauliz.pz == 0.5
+end
+
+function test_pauliz_application_guaranteed()
+    n_bellpairs = 3
+    changed_bp_idx = 2
+    test_state = BellState(n_bellpairs)
+    test_guaranteed_pauliz = PauliZOp(changed_bp_idx, 1)
+    
+    QuantumClifford.apply!(test_state, test_guaranteed_pauliz)
+
+    @test test_state.phases[2*changed_bp_idx - 1] == 0
+    @test test_state.phases[2*changed_bp_idx] == 1
+end
+
+function test_pauliz_application_guaranteed_none()
+    n_bellpairs = 3
+    changed_bp_idx = 2
+    test_state = BellState(n_bellpairs)
+    test_guaranteed_pauliz = PauliZOp(changed_bp_idx, 0)
+    
+    # TODO: do I have to import QuantumClifford to use it when testing?
+    QuantumClifford.apply!(test_state, test_guaranteed_pauliz)
+
+    @test test_state.phases[2*changed_bp_idx - 1] == 0
+    @test test_state.phases[2*changed_bp_idx] == 0
+end
+
+function test_mixed_state_op_constructor()
+    mixed_state_op = MixedStateOp(1, 0.5)
+    
+    @test mixed_state_op isa MixedStateOp
+
+    @test mixed_state_op.idx == 1
+
+    @test mixed_state_op.lambda == 0.5
+end
+
+function test_apply_mixed_state_op()
+    n_bellpairs = 1
+    changed_bp_idx = 1
+    test_state = BellState(n_bellpairs)
+    mixed_state_op = MixedStateOp(changed_bp_idx, 0.0)
+
+    QuantumClifford.apply!(test_state, mixed_state_op)
+
+    @test test_state.phases[2 * changed_bp_idx - 1] == 0
+    @test test_state.phases[2 * changed_bp_idx] == 0
+end
+
+function test_apply_mixed_state_op_diff_bellstate()
+    n_bellpairs = 1
+    changed_bp_idx = 1
+    test_state = BellState((0, 1))
+    mixed_state_op = MixedStateOp(changed_bp_idx, 0.0)
+    
+    QuantumClifford.apply!(test_state, mixed_state_op)
+    
+    @test test_state.phases[2 * changed_bp_idx - 1] == 0
+    @test test_state.phases[2 * changed_bp_idx] == 1
+end
+
+function test_apply_both_memory_errors()
+    n_bellpairs = 1
+    changed_bp_idx = 1
+    
+    test_state = BellState(n_bellpairs)
+
+    noise_ops = Vector{AbstractNoiseBellOp}()
+    
+    push!(noise_ops, PauliZOp(changed_bp_idx, 0))
+    push!(noise_ops, MixedStateOp(changed_bp_idx, 0.0))
+
+    for noise_op in noise_ops
+        QuantumClifford.apply!(test_state, noise_op)
+    end
+
+    @test test_state.phases[2 * changed_bp_idx - 1] == 0
+    @test test_state.phases[2 * changed_bp_idx] == 0
+end
+
+@testset "BPGates.jl, PauliZOp tests" begin
+    test_pauliz_constructor()
+    test_pauliz_application_guaranteed()
+    test_pauliz_application_guaranteed_none()
+end
+
+@testset "BPGates.jl, MixedStateOp tests" begin
+    test_mixed_state_op_constructor()
+    test_apply_mixed_state_op()
+    test_apply_mixed_state_op_diff_bellstate()
+    test_apply_both_memory_errors()
+end
diff --git a/test/test_kraus_mem_noise_helpers.jl b/test/test_kraus_mem_noise_helpers.jl
new file mode 100644
index 0000000..e0bc996
--- /dev/null
+++ b/test/test_kraus_mem_noise_helpers.jl
@@ -0,0 +1,259 @@
+module TestKrausNoiseHelpers
+
+using Test
+using Logging
+
+using BPGates: get_mixed_transition_probs, select_rand_state
+
+export test_bell_basis_probs, test_nonoise_bell_basis, test_almost_all_noise,
+test_random_lambdas, test_equal_mixed_probs_2, test_equal_mixed_probs_4,
+test_random_mixed_states, test_random_mixed_states_random_lambdas
+
+# using QuantumClifford: apply!
+
+# Get the bell states used for testing
+function get_bell_states()
+    return (1 / sqrt(2)) * [
+        [1.0; 0.0; 0.0; 1.0],
+        [0.0; 1.0; 1.0; 0.0],
+        [1.0; 0.0; 0.0; -1.0],
+        [0.0; 1.0; -1.0; 0.0]
+    ]
+end
+
+# Returns an index selected from the indices of a list, weighted by the probability at that index
+function get_weighted_index(prob_dist::Vector{Float64})
+    if !isapprox(sum(prob_dist), 1.0, atol=1e-9)
+        @error "get_weighted_index: prob_dist sums to $(sum(prob_dist)) and not 1"
+    end
+    rand_num = rand()
+    cur_sum = 0.0
+    for (idx, prob_val) in enumerate(prob_dist)
+        if rand_num < cur_sum + prob_val
+            return idx
+        end
+        cur_sum += prob_val
+    end
+end
+
+# Converts a bit to an int from left to right (LSB is the rightmost bit)
+function bit_to_int_testing(bits)
+    reduce(⊻,(bit<<(length(bits) - index) for (index,bit) in enumerate(bits))) + 1 # +1 so that we use julia indexing convenctions
+end
+
+# Returns the tensor product of a set of Kraus operators
+function get_tensored_kraus_ops(kraus_ops::Vector{Matrix{Complex}})
+    tensored_kraus_ops = []
+    for op1 in kraus_ops
+        for op2 in kraus_ops
+            push!(tensored_kraus_ops, kron(op1, op2))
+        end
+    end
+    return tensored_kraus_ops
+end
+
+# Returns the diagonal of a matrix
+function get_diagonal(matrix::Matrix{ComplexF64})::Vector{Float64}
+    diagonal_length = min(size(matrix)[1], size(matrix)[2])
+    if size(matrix)[1] != size(matrix)[2]
+        @error("get_diagonal: matrix is not square, rows: $(size(matrix)[1]) != cols: $(size(matrix)[2])")
+    end
+    diag_elements = []
+    for idx in 1:diagonal_length
+        push!(diag_elements, matrix[idx, idx])
+    end
+    return diag_elements
+end
+
+# Returns the bell state corresponding to the input phases
+function get_bell_state(phases::BitVector)
+    state_idx = bit_to_int_testing(phases)
+    bell_states = get_bell_states()
+    return bell_states[state_idx]
+end
+
+# Numerically computes the transition probabilities of one state to another,
+# given a mixed probability distribution of states and a set of kraus operators
+function get_numerical_trans_probs(mixed_state_probs::Vector{Float64}, kraus_ops::Vector{Matrix{Complex}})
+    bell_states = get_bell_states()
+    ρ_bell = zeros(ComplexF64, length(bell_states), length(bell_states))
+    for (bell_idx, bell_prob) in enumerate(mixed_state_probs)
+        ρ_bell += bell_prob * bell_states[bell_idx] * bell_states[bell_idx]'
+    end
+    @debug("get_numerical_trans_probs: ρ_bell: $ρ_bell")
+    tensored_ops = get_tensored_kraus_ops(kraus_ops)
+    # Apply the Kraus operators to the Bell state
+    ρ_new = zeros(ComplexF64, size(ρ_bell)[1], size(ρ_bell)[2])
+    for K in tensored_ops
+        ρ_new += K * ρ_bell * adjoint(K)
+    end
+    U_bell = hcat(get_bell_states()...)
+    @debug("get_numerical_trans_probs: ρ_new: $ρ_new")
+    @debug("get_numerical_trans_probs: U_bell: $U_bell")
+    ρ_Bell_basis = adjoint(U_bell) * ρ_new * U_bell
+    @debug("get_numerical_trans_probs: ρ_Bell_basis: $ρ_Bell_basis")
+    @debug("get_numerical_trans_probs: size(ρ_Bell_basis): $(size(ρ_Bell_basis))")
+    return get_diagonal(ρ_Bell_basis)
+end
+
+# Testing functions
+
+# Test our simulated noise for the pure states in the bell basis
+function test_bell_basis_probs(kraus_ops_factory::Function, get_simulated_trans_probs::Function)
+    statedist_list = [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]]
+    lambda = 1 - exp(-500e-9 / 260e-6)
+    num_samples = 1000000
+    kraus_ops = kraus_ops_factory(lambda)
+    for state_dist in statedist_list
+        @info "test_bell_basis_probs: state_dist: $state_dist"
+        numerical_trans_probs = get_numerical_trans_probs(state_dist, kraus_ops)
+        simulated_trans_probs = get_simulated_trans_probs(state_dist, lambda, num_samples)
+        # @info("test_A_state_probs: numerical_trans_probs: $numerical_trans_probs")
+        # @info("test_A_state_probs: simulated_trans_probs: $simulated_trans_probs")
+        @test all(isapprox.(numerical_trans_probs, simulated_trans_probs, atol=(2 / sqrt(num_samples))))
+    end
+end
+
+# Test our simulated noise when there should be no noise applied for states in the bell basis
+function test_nonoise_bell_basis(kraus_ops_factory::Function, get_simulated_trans_probs::Function)
+    statedist_list = [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]]
+    lambda = 0.0
+    num_samples = 1000000
+    kraus_ops = kraus_ops_factory(lambda)
+    for state_dist in statedist_list
+        @info "test_nonoise_bell_basis: state_dist: $state_dist"
+        numerical_trans_probs = get_numerical_trans_probs(state_dist, kraus_ops)
+        simulated_trans_probs = get_simulated_trans_probs(state_dist, lambda, num_samples)
+        # @info("test_A_state_probs: numerical_trans_probs: $numerical_trans_probs")
+        # @info("test_A_state_probs: simulated_trans_probs: $simulated_trans_probs")
+        @test all(isapprox.(numerical_trans_probs, simulated_trans_probs, atol=(2 / sqrt(num_samples))))
+    end
+end
+
+# Test our simulated noise with a noise factor close to 1 for states in the bell basis
+function test_almost_all_noise(kraus_ops_factory::Function, get_simulated_trans_probs::Function)
+    statedist_list = [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]]
+    lambda = 0.999999999
+    num_samples = 1000000
+    kraus_ops = kraus_ops_factory(lambda)
+    for state_dist in statedist_list
+        @info "test_almost_all_noise: state_dist: $state_dist"
+        numerical_trans_probs = get_numerical_trans_probs(state_dist, kraus_ops)
+        simulated_trans_probs = get_simulated_trans_probs(state_dist, lambda, num_samples)
+        @info("test_almost_all_noise: numerical_trans_probs: $numerical_trans_probs")
+        @info("test_almost_all_noise: simulated_trans_probs: $simulated_trans_probs")
+        @test all(isapprox.(numerical_trans_probs, simulated_trans_probs, atol=(2 / sqrt(num_samples))))
+    end
+end
+
+# Test random noise factors for states in the bell basis
+function test_random_lambdas(kraus_ops_factory::Function, get_simulated_trans_probs::Function)
+    statedist_list = [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]]
+    num_samples = 1000000
+    num_random_lambdas = 5
+    for _ in 1:num_random_lambdas
+        lambda = rand()
+        num_samples = 1000000
+        kraus_ops = kraus_ops_factory(lambda)
+        for state_dist in statedist_list
+            @info "test_random_lambdas: state_dist: $state_dist"
+            @info "test_random_lambdas: lambda: $lambda"
+            numerical_trans_probs = get_numerical_trans_probs(state_dist, kraus_ops)
+            simulated_trans_probs = get_simulated_trans_probs(state_dist, lambda, num_samples)
+            @info("test_random_lambdas: numerical_trans_probs: $numerical_trans_probs")
+            @info("test_random_lambdas: simulated_trans_probs: $simulated_trans_probs")
+            @test all(isapprox.(numerical_trans_probs, simulated_trans_probs, atol=(2 / sqrt(num_samples))))
+        end
+    end
+end
+
+# Test all mixed states composed of two states in the Bell basis with equal probability
+function test_equal_mixed_probs_2(kraus_ops_factory::Function, get_simulated_trans_probs::Function)
+    bell_states = get_bell_states()
+    lambda = 1 - exp(-500e-9 / 260e-6)
+    kraus_ops = kraus_ops_factory(lambda)
+    num_samples = 1000000
+    for first_phases_idx in 1:length(bell_states)
+        for second_phases_idx in 1:length(bell_states)
+            if first_phases_idx == second_phases_idx
+                continue
+            end
+            state_dist = [0.0 for _ in 1:length(bell_states)]
+            for basis_idx in 1:length(bell_states)
+                if basis_idx == first_phases_idx || basis_idx == second_phases_idx
+                    state_dist[basis_idx] = 0.5
+                end
+            end
+            numerical_trans_probs = get_numerical_trans_probs(state_dist, kraus_ops)
+            simulated_trans_probs = get_simulated_trans_probs(state_dist, lambda, num_samples)
+            @info("test_equal_mixed_probs_2: numerical_trans_probs: $numerical_trans_probs")
+            @info("test_equal_mixed_probs_2: simulated_trans_probs: $simulated_trans_probs")
+            @test all(isapprox.(numerical_trans_probs, simulated_trans_probs, atol=(2 / sqrt(num_samples))))
+        end
+    end
+end
+
+# Test a mixed state composed of all 4 bell basis states with equal probability
+function test_equal_mixed_probs_4(kraus_ops_factory::Function, get_simulated_trans_probs::Function)
+    state_dist = [0.25 for _ in 1:length(get_bell_states())]
+    lambda = 1 - exp(-500e-9 / 260e-6)
+    num_samples = 1000000
+    kraus_ops = kraus_ops_factory(lambda)
+    @info "test_equal_mixed_probs_4: state_dist: $state_dist"
+    numerical_trans_probs = get_numerical_trans_probs(state_dist, kraus_ops)
+    simulated_trans_probs = get_simulated_trans_probs(state_dist, lambda, num_samples)
+    @info("test_equal_mixed_probs_4: numerical_trans_probs: $numerical_trans_probs")
+    @info("test_equal_mixed_probs_4: simulated_trans_probs: $simulated_trans_probs")
+    @test all(isapprox.(numerical_trans_probs, simulated_trans_probs, atol=(2 / sqrt(num_samples))))
+end
+
+# Test random mixed states in the bell basis
+function test_random_mixed_states(kraus_ops_factory::Function, get_simulated_trans_probs::Function)
+    num_random_states = 20
+    statedist_list = []
+    for _ in 1:num_random_states
+        state_dist = [rand() for _ in 1:length(get_bell_states())]
+        state_dist = state_dist / sum(state_dist)
+        push!(statedist_list, state_dist)
+    end
+    @debug "length(statedist_list): $(length(statedist_list))"
+    lambda = 1 - exp(-500e-9 / 260e-6)
+    num_samples = 1000000
+    kraus_ops = kraus_ops_factory(lambda)
+    for state_dist in statedist_list
+        @info "test_random_mixed_states: state_dist: $state_dist"
+        numerical_trans_probs = get_numerical_trans_probs(state_dist, kraus_ops)
+        simulated_trans_probs = get_simulated_trans_probs(state_dist, lambda, num_samples)
+        @info("test_random_mixed_states: numerical_trans_probs: $numerical_trans_probs")
+        @info("test_random_mixed_states: simulated_trans_probs: $simulated_trans_probs")
+        @test all(isapprox.(numerical_trans_probs, simulated_trans_probs, atol=(2 / sqrt(num_samples))))
+    end
+end
+
+# Test random mixed states with random noise parameters in the Bell basis
+function test_random_mixed_states_random_lambdas(kraus_ops_factory::Function, get_simulated_trans_probs::Function)
+    num_random_states = 20
+    statedist_list = []
+    for _ in 1:num_random_states
+        state_dist = [rand() for _ in 1:length(get_bell_states())]
+        state_dist = state_dist / sum(state_dist)
+        push!(statedist_list, state_dist)
+    end
+    num_random_lambdas = 20
+    num_samples = 1000000
+    for _ in 1:num_random_lambdas
+        lambda = rand()
+        @info "test_random_mixed_states_random_lambdas: lambda: $lambda"
+        kraus_ops = kraus_ops_factory(lambda)
+        for state_dist in statedist_list
+            @info "test_random_mixed_states_random_lambdas: state_dist: $state_dist"
+            numerical_trans_probs = get_numerical_trans_probs(state_dist, kraus_ops)
+            simulated_trans_probs = get_simulated_trans_probs(state_dist, lambda, num_samples)
+            @info("test_random_mixed_states_random_lambdas: numerical_trans_probs: $numerical_trans_probs")
+            @info("test_random_mixed_states_random_lambdas: simulated_trans_probs: $simulated_trans_probs")
+            @test all(isapprox.(numerical_trans_probs, simulated_trans_probs, atol=(2 / sqrt(num_samples))))
+        end
+    end
+end
+
+end
\ No newline at end of file
diff --git a/test/test_mixedstateop.jl b/test/test_mixedstateop.jl
new file mode 100644
index 0000000..d8e13cc
--- /dev/null
+++ b/test/test_mixedstateop.jl
@@ -0,0 +1,44 @@
+using .TestKrausNoiseHelpers
+using Logging
+
+# TODO: have tolerance for numerical precision for noise calculations
+
+# Set up the logger to enable debug level
+function enable_debug_logging()
+    global_logger(ConsoleLogger(stderr, Logging.Info))
+end
+
+# Define kraus operators to test
+function get_ampdamp_kraus_ops(λ::Float64)::Vector{Matrix{Complex}}
+    # Define Kraus operators
+    E1 = [1 0; 0 sqrt(1 - λ)]
+    E2 = [0 sqrt(λ); 0 0]
+    return [E1, E2]
+end
+
+# Returns the simulated transition probabilities from the mixed state probabilities,
+# the noise parameter, simulated for num_samples times
+function get_simulated_trans_probs(mixed_state_probs::Vector{Float64}, noise_param::Float64, num_samples::Integer)
+    sampled_state_count = fill(0, length(TestKrausNoiseHelpers.get_bell_states()))
+    for _ in 1:num_samples
+        state_idx = TestKrausNoiseHelpers.get_weighted_index(mixed_state_probs)
+        transition_probs = TestKrausNoiseHelpers.get_mixed_transition_probs(noise_param, state_idx)
+        new_state_idx = TestKrausNoiseHelpers.select_rand_state(transition_probs)
+        # TODO: ensure that the index is inbounds for my list
+        sampled_state_count[new_state_idx] += 1
+    end
+    return sampled_state_count / num_samples
+end
+
+# Run all tests
+@testset "BPGates.jl, get_mixed_transition_probs tests" begin
+    enable_debug_logging()
+    test_bell_basis_probs(get_ampdamp_kraus_ops, get_simulated_trans_probs)
+    test_nonoise_bell_basis(get_ampdamp_kraus_ops, get_simulated_trans_probs)
+    test_almost_all_noise(get_ampdamp_kraus_ops, get_simulated_trans_probs)
+    test_random_lambdas(get_ampdamp_kraus_ops, get_simulated_trans_probs)
+    test_equal_mixed_probs_2(get_ampdamp_kraus_ops, get_simulated_trans_probs)
+    test_equal_mixed_probs_4(get_ampdamp_kraus_ops, get_simulated_trans_probs)
+    test_random_mixed_states(get_ampdamp_kraus_ops, get_simulated_trans_probs)
+    test_random_mixed_states_random_lambdas(get_ampdamp_kraus_ops, get_simulated_trans_probs)
+end
diff --git a/test/test_paulizop.jl b/test/test_paulizop.jl
new file mode 100644
index 0000000..0e49500
--- /dev/null
+++ b/test/test_paulizop.jl
@@ -0,0 +1,70 @@
+using Test
+using Logging
+
+using BPGates: get_mixed_transition_probs, select_rand_state
+
+using .TestKrausNoiseHelpers
+
+using QuantumClifford: apply!
+
+# TODO: have tolerance for numerical precision for noise calculations
+
+# Set up the logger to enable debug level
+function enable_debug_logging()
+    global_logger(ConsoleLogger(stderr, Logging.Info))
+end
+
+# Define kraus operators to test
+function get_phasedamp_kraus_ops(λ::Float64)::Vector{Matrix{Complex}}
+    # Define Kraus operators
+    E1 = sqrt(1 - (λ / 2)) * [1 0; 0 1]
+    E2 = sqrt(λ / 2) * [1 0; 0 -1]
+    return [E1, E2]
+end
+
+function do_z_flip(state_idx::Int)
+    if state_idx == 1
+        return 3
+    elseif state_idx == 2
+        return 4
+    elseif state_idx == 3
+        return 1
+    elseif state_idx == 4
+        return 2
+    else
+        @error "do_z_flip: state_idx is not between 1 and 4, state_idx: $state_idx"
+    end
+end
+
+# Returns the simulated transition probabilities from the mixed state probabilities,
+# the noise parameter, simulated for num_samples times
+function get_simulated_trans_probs(mixed_state_probs::Vector{Float64}, noise_param::Float64, num_samples::Integer)
+    sampled_state_count = fill(0, length(TestKrausNoiseHelpers.get_bell_states()))
+    for _ in 1:num_samples
+        state_idx = TestKrausNoiseHelpers.get_weighted_index(mixed_state_probs)
+        # transition_probs = TestKrausNoiseHelpers.get_mixed_transition_probs(noise_param, state_idx)
+        # new_state_idx = TestKrausNoiseHelpers.select_rand_state(transition_probs)
+        r = rand()
+        if r < 2 * (noise_param / 2) * (1 - (noise_param / 2))
+            new_state_idx = do_z_flip(state_idx)
+        else
+            new_state_idx = state_idx
+        end
+        # TODO: ensure that the index is inbounds for my list
+        sampled_state_count[new_state_idx] += 1
+    end
+    return sampled_state_count / num_samples
+end
+
+# Run all tests
+@testset "BPGates.jl, get_mixed_transition_probs tests" begin
+    enable_debug_logging()
+    test_bell_basis_probs(get_phasedamp_kraus_ops, get_simulated_trans_probs)
+    test_nonoise_bell_basis(get_phasedamp_kraus_ops, get_simulated_trans_probs)
+    test_almost_all_noise(get_phasedamp_kraus_ops, get_simulated_trans_probs)
+    test_random_lambdas(get_phasedamp_kraus_ops, get_simulated_trans_probs)
+    test_equal_mixed_probs_2(get_phasedamp_kraus_ops, get_simulated_trans_probs)
+    test_equal_mixed_probs_4(get_phasedamp_kraus_ops, get_simulated_trans_probs)
+    test_random_mixed_states(get_phasedamp_kraus_ops, get_simulated_trans_probs)
+    test_random_mixed_states_random_lambdas(get_phasedamp_kraus_ops, get_simulated_trans_probs)
+end