diff --git a/src/QuantumClifford.jl b/src/QuantumClifford.jl
index d54ec03b6..9ad6488c1 100644
--- a/src/QuantumClifford.jl
+++ b/src/QuantumClifford.jl
@@ -439,6 +439,114 @@ See also: [`stabilizerview`](@ref), [`destabilizerview`](@ref), [`logicalxview`]
 tab(s::Stabilizer{T}) where {T} = s.tab
 tab(s::AbstractStabilizer) = s.tab
 
+"""
+Embed a Tableau in a larger Tableau.
+
+```jldoctest
+julia> embed(10, 5, T"Z X")
++ ____Z_____
++ ____X_____
+
+julia> embed(10, (1, 5), T"-XZ IZ YI YX")
+- X___Z_____
++ ____Z_____
++ Y_________
++ Y___X_____
+```
+"""
+function embed(n::Int, i::Int, t::Tableau)
+     if nqubits(t) == 1
+        tout = zero(typeof(t), size(t)[1], n)
+        for j in 1:size(t)[1]
+            tout[j, i] = t[j][1]
+            tout.phases[j] = t.phases[j]
+        end
+        return tout
+    else
+        throw(ArgumentError("""
+        You are trying to embed a small Tableau into a larger Tableau.
+        However, you have not given all the positions at which the Tableau needs to be embedded.
+        If you are directly calling `embed`, use the form `embed(nqubits, indices::Tuple, t::Tableau)`.
+        If you are not using `embed` directly, then `embed` must have been incorrectly called
+        by one of the functions you have called.
+        """))
+    end
+end
+
+function embed(n::Int, indices, t::Tableau)
+    if nqubits(t) == length(indices)
+        tout = zero(typeof(t), size(t)[1], n)
+        @inbounds @simd for j in 1:size(t)[1]
+            @inbounds @simd for i in eachindex(indices)
+                tout[j, indices[i]] = t[j][i]
+                tout.phases[j] = t.phases[j]
+            end
+        end
+        return tout
+    else
+        throw(ArgumentError(lazy"""
+        You are trying to embed a small Tableau into a larger Tableau.
+        However, you have not given all the positions at which the Tableau needs to be embedded.
+        If you are directly calling `embed`, use the form `embed(nqubits, indices::Tuple, t::Tableau)`.
+        If you are not using `embed` directly, then `embed` must have been incorrectly called
+        by one of the functions you have called.
+        """))
+    end
+end
+
+"""
+Embed a Stabilizer in a larger Stabilizer.
+
+```jldoctest
+julia> embed(10, 5, S"-Y Z")
+- ____Y_____
++ ____Z_____
+
+julia> embed(10, (1,9), S"XX -YZ")
++ X_______X_
+- Y_______Z_
+```
+"""
+function embed(n::Int, i::Int, s::AbstractStabilizer)
+     if nqubits(s) == 1
+        sout = zero(typeof(s), size(s)[1], n)
+        for j in 1:size(s)[1]
+            tab(sout)[j, i] = tab(s)[j][1]
+            tab(sout).phases[j] = tab(s).phases[j]
+        end
+        return sout
+    else
+        throw(ArgumentError("""
+        You are trying to embed a small Stabilizer into a larger Stabilizer.
+        However, you have not given all the positions at which the Stabilizer needs to be embedded.
+        If you are directly calling `embed`, use the form `embed(nqubits, indices::Tuple, s::AbstractStabilizer)`.
+        If you are not using `embed` directly, then `embed` must have been incorrectly called
+        by one of the functions you have called.
+        """))
+    end
+end
+
+function embed(n::Int, indices, s::AbstractStabilizer)
+    if nqubits(s) == length(indices)
+        sout = zero(typeof(s), size(s)[1], n)
+        @inbounds @simd for j in 1:size(s)[1]
+            @inbounds @simd for i in eachindex(indices)
+                tab(sout)[j, indices[i]] = tab(s)[j][i]
+                tab(sout).phases[j] = tab(s).phases[j]
+            end
+        end
+        return sout
+    else
+        throw(ArgumentError(lazy"""
+        You are trying to embed a small Stabilizer into a larger Stabilizer.
+        However, you have not given all the positions at which the Stabilizer needs to be embedded.
+        If you are directly calling `embed`, use the form `embed(nqubits, indices::Tuple, s::AbstractStabilizer)`.
+        If you are not using `embed` directly, then `embed` must have been incorrectly called
+        by one of the functions you have called.
+        """))
+    end
+end
+
 ##############################
 # Destabilizer formalism
 ##############################
diff --git a/test/test_embed.jl b/test/test_embed.jl
index 37d5c0df7..083631892 100644
--- a/test/test_embed.jl
+++ b/test/test_embed.jl
@@ -1,7 +1,30 @@
-@testitem "embed PauliOperator" begin
-    @test embed(5,3,P"-Y") == P"-__Y__"
-    @test embed(5,(3,5),P"-YZ") == P"-__Y_Z"
-    @test embed(5,[3,5],P"-YZ") == P"-__Y_Z"
-    @test_throws ArgumentError embed(5,[3,5],P"-YZX")
-    @test_throws ArgumentError embed(5,3,P"-ZX")
+@testitem "embed methods" begin
+
+    @testset "embed PauliOperator" begin
+        @test embed(5,3,P"-Y") == P"-__Y__"
+        @test embed(5,(3,5),P"-YZ") == P"-__Y_Z"
+        @test embed(5,[3,5],P"-YZ") == P"-__Y_Z"
+        @test_throws ArgumentError embed(5,[3,5],P"-YZX")
+        @test_throws ArgumentError embed(5,3,P"-ZX")
+    end
+
+    @testset "embed Stabilizer" begin
+        @test embed(5,3,S"-Y") == S"-__Y__"
+        @test embed(5,(3,5),S"-YZ") == S"-__Y_Z"
+        @test embed(5,[3,5],S"-YZ") == S"-__Y_Z"
+        @test embed(5, 5, S"-Y Z") == S"-____Y ____Z"
+        @test embed(5, (2,4), S"XX -YZ") == S"_X_X_ -_Y_Z_"
+        @test_throws ArgumentError embed(5,[3,5],S"-YZX")
+        @test_throws ArgumentError embed(5,3,S"-ZX")
+    end
+
+    @testset "embed Tableau" begin
+        @test embed(5,3,T"-Y") == T"-__Y__"
+        @test embed(5,(3,5),T"-YZ") == T"-__Y_Z"
+        @test embed(5,[3,5],T"-YZ") == T"-__Y_Z"
+        @test embed(5, 5, T"-Y Z") == T"-____Y ____Z"
+        @test embed(5, (2,4), T"XX -YZ") == T"_X_X_ -_Y_Z_"
+        @test_throws ArgumentError embed(5,[3,5],T"-YZX")
+        @test_throws ArgumentError embed(5,3,T"-ZX")
+    end
 end
diff --git a/test/test_throws.jl b/test/test_throws.jl
index 2c67f4f4a..e49c08dcd 100644
--- a/test/test_throws.jl
+++ b/test/test_throws.jl
@@ -1,66 +1,72 @@
 @testitem "throws" begin
-  using QuantumClifford: rank, mul_left!, mul_right!
-  using InteractiveUtils: subtypes
+    using QuantumClifford: rank, mul_left!, mul_right!
+    using InteractiveUtils: subtypes
 
-  @test_throws DimensionMismatch CliffordOperator(T"XXX ZZ_")
+    @test_throws DimensionMismatch CliffordOperator(T"XXX ZZ_")
 
-  @test_throws DimensionMismatch tCNOT*S"X"
+    @test_throws DimensionMismatch tCNOT*S"X"
 
-  #@test_throws DomainError bigram(random_stabilizer(50), clip=false)
+    #@test_throws DomainError bigram(random_stabilizer(50), clip=false)
 
-  @test_throws DomainError logdot(S"XX", S"XX ZZ")
-  @test_throws DimensionMismatch logdot(S"X", S"XX ZZ")
+    @test_throws DomainError logdot(S"XX", S"XX ZZ")
+    @test_throws DimensionMismatch logdot(S"X", S"XX ZZ")
 
-  @test_throws BadDataStructure rank(S"X")
-  @test_throws BadDataStructure rank(Destabilizer(S"X"))
+    @test_throws BadDataStructure rank(S"X")
+    @test_throws BadDataStructure rank(Destabilizer(S"X"))
 
-  @test_throws DimensionMismatch mul_left!(P"X", P"XX")
-  @test_throws DimensionMismatch mul_right!(P"X", P"XX")
+    @test_throws DimensionMismatch mul_left!(P"X", P"XX")
+    @test_throws DimensionMismatch mul_right!(P"X", P"XX")
 
-  @test_throws ArgumentError StabMixture(S"XX")
+    @test_throws ArgumentError StabMixture(S"XX")
 
-  @test_throws ArgumentError UnitaryPauliChannel([P"X"], [1,2])
-  @test_throws ArgumentError UnitaryPauliChannel([P"X",P"XX"], [1,2])
+    @test_throws ArgumentError UnitaryPauliChannel([P"X"], [1,2])
+    @test_throws ArgumentError UnitaryPauliChannel([P"X",P"XX"], [1,2])
 
-  @test_throws ArgumentError embed(10,2,P"XX")
-  @test_throws ArgumentError embed(10,[2,3],P"X")
+    @test_throws ArgumentError embed(10,2,P"XX")
+    @test_throws ArgumentError embed(10,[2,3],P"X")
 
-  struct A <: QuantumClifford.AbstractOperation end
-  @test_throws ArgumentError applybranches(S"X",A())
+    @test_throws ArgumentError embed(10,2,S"XX")
+    @test_throws ArgumentError embed(10,[2,3],S"X")
 
-  @test_throws BadDataStructure project!(Destabilizer(S"XX"), P"ZZ")
+    @test_throws ArgumentError embed(10,2,T"XX")
+    @test_throws ArgumentError embed(10,[2,3],T"X")
 
-  @test_throws DimensionMismatch reset_qubits!(ghz(4), ghz(3), [1,2])
-  @test_throws DimensionMismatch reset_qubits!(ghz(3), ghz(4), [1,2,3,4])
-  @test_throws DimensionMismatch reset_qubits!(MixedStabilizer(ghz(4)), MixedStabilizer(ghz(3)), [1,2])
-  @test_throws DimensionMismatch reset_qubits!(MixedStabilizer(ghz(3)), MixedStabilizer(ghz(4)), [1,2,3,4])
-  @test_throws DimensionMismatch reset_qubits!(MixedDestabilizer(ghz(4)), MixedDestabilizer(ghz(3)), [1,2])
-  @test_throws DimensionMismatch reset_qubits!(MixedDestabilizer(ghz(3)), MixedDestabilizer(ghz(4)), [1,2,3,4])
+    struct A <: QuantumClifford.AbstractOperation end
+    @test_throws ArgumentError applybranches(S"X",A())
 
-  #TODO broken in other ways @test_throws DomainError MixedDestabilizer(Destabilizer(S"XX"))
+    @test_throws BadDataStructure project!(Destabilizer(S"XX"), P"ZZ")
 
-  @test_throws DomainError 2*P"X"
+    @test_throws DimensionMismatch reset_qubits!(ghz(4), ghz(3), [1,2])
+    @test_throws DimensionMismatch reset_qubits!(ghz(3), ghz(4), [1,2,3,4])
+    @test_throws DimensionMismatch reset_qubits!(MixedStabilizer(ghz(4)), MixedStabilizer(ghz(3)), [1,2])
+    @test_throws DimensionMismatch reset_qubits!(MixedStabilizer(ghz(3)), MixedStabilizer(ghz(4)), [1,2,3,4])
+    @test_throws DimensionMismatch reset_qubits!(MixedDestabilizer(ghz(4)), MixedDestabilizer(ghz(3)), [1,2])
+    @test_throws DimensionMismatch reset_qubits!(MixedDestabilizer(ghz(3)), MixedDestabilizer(ghz(4)), [1,2,3,4])
 
-  @test_throws DimensionMismatch P"X" * S"XX"
+    #TODO broken in other ways @test_throws DomainError MixedDestabilizer(Destabilizer(S"XX"))
 
-  @test_throws ArgumentError one(typeof(T"X"), 1, basis=:U)
+    @test_throws DomainError 2*P"X"
 
-  for gt in subtypes(QuantumClifford.AbstractSingleQubitOperator)
-      gt == SingleQubitOperator && continue
-      @test_throws ArgumentError gt(0)
-      @test_throws ArgumentError gt(-1)
-  end
+    @test_throws DimensionMismatch P"X" * S"XX"
 
-  for gt in subtypes(QuantumClifford.AbstractTwoQubitOperator)
-      @test_throws ArgumentError gt(0,1)
-      @test_throws ArgumentError gt(-1,1)
-      @test_throws ArgumentError gt(2,2)
-  end
+    @test_throws ArgumentError one(typeof(T"X"), 1, basis=:U)
 
-  for m in [sMX,sMZ,sMY,sMRX,sMRZ,sMRY]
-      @test_throws ArgumentError m(0)
-      @test_throws ArgumentError m(-1)
-      @test_throws ArgumentError m(0,1)
-      @test_throws ArgumentError m(-1,0)
-  end
+    for gt in subtypes(QuantumClifford.AbstractSingleQubitOperator)
+        gt == SingleQubitOperator && continue
+        @test_throws ArgumentError gt(0)
+        @test_throws ArgumentError gt(-1)
+    end
+
+    for gt in subtypes(QuantumClifford.AbstractTwoQubitOperator)
+        @test_throws ArgumentError gt(0,1)
+        @test_throws ArgumentError gt(-1,1)
+        @test_throws ArgumentError gt(2,2)
+    end
+
+    for m in [sMX,sMZ,sMY,sMRX,sMRZ,sMRY]
+        @test_throws ArgumentError m(0)
+        @test_throws ArgumentError m(-1)
+        @test_throws ArgumentError m(0,1)
+        @test_throws ArgumentError m(-1,0)
+    end
 end