diff --git a/Project.toml b/Project.toml
index 42aae1d..b49002d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,3 +5,5 @@ version = "0.1.0"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
+QEDbase = "10e22c08-3ccb-4172-bfcf-7d7aa3d04d93"
+QEDprocesses = "46de9c38-1bb3-4547-a1ec-da24d767fdad"
diff --git a/src/FeynmanDiagramGenerator.jl b/src/FeynmanDiagramGenerator.jl
index 38de335..95a127d 100644
--- a/src/FeynmanDiagramGenerator.jl
+++ b/src/FeynmanDiagramGenerator.jl
@@ -10,9 +10,12 @@ export external_particles, virtual_particles, process, vertices
 
 export Forest
 
+export plane_trees, labelled_plane_trees
+
 export can_interact
 export QED
 
+include("trees/labelled_plane_trees.jl")
 include("trees/trees.jl")
 include("trees/iterator.jl")
 include("trees/print.jl")
@@ -23,5 +26,6 @@ include("diagrams/interface.jl")
 include("diagrams/vertex.jl")
 include("diagrams/diagrams.jl")
 include("diagrams/iterator.jl")
+include("diagrams/QED.jl")
 
 end # module FeynmanDiagramGenerator
diff --git a/src/diagrams/QED.jl b/src/diagrams/QED.jl
index bc3f110..5a46acc 100644
--- a/src/diagrams/QED.jl
+++ b/src/diagrams/QED.jl
@@ -16,7 +16,9 @@ function plane_trees(n::Int64)
 end
 
 function iterate(lpt::LabelledPlaneTreeIterator)
-    return ()
+    state = LabelledPlaneTreeIteratorState()
+
+    
 end
 
 function iterate(lpt::LabelledPlaneTreeIterator, state)
@@ -26,5 +28,6 @@ function iterate(lpt::LabelledPlaneTreeIterator, state)
             i += 1
         end
     end
-end
 
+    return (1, )
+end
diff --git a/src/qft/qft.jl b/src/qft/qft.jl
index fe6e341..cd6875c 100644
--- a/src/qft/qft.jl
+++ b/src/qft/qft.jl
@@ -32,7 +32,5 @@ end
 can_interact(theory::AbstractQFT, p1::AbstractParticle, p2::AbstractParticle, p3::AbstractParticle) = can_interact(theory, Vertex(p1, p2, p3))
 
 QED = QFT([
-    Vertex(Electron(), Positron(), Photon()),
-    Vertex(Muon(), AntiMuon(), Photon()),
-    Vertex(Tauon(), AntiTauon(), Photon())
+    Vertex(Electron(), Positron(), Photon())
 ])
diff --git a/src/trees/labelled_plane_trees.jl b/src/trees/labelled_plane_trees.jl
new file mode 100644
index 0000000..b10970d
--- /dev/null
+++ b/src/trees/labelled_plane_trees.jl
@@ -0,0 +1,134 @@
+using Combinatorics
+
+function _gen_trees(vec::Vector{Tuple{Int64,Vector{Int64}}})
+    # input is one tree structure
+    # node neighbours are partially filled from the front
+
+    n = length(vec)
+
+    # recursion ends if only two nodes are left with 1 
+    sum_open_edges = sum(getindex.(vec, 1)) # total number of open edges
+    sum_open_edges -= sum(length.(getindex.(vec, 2))) # minus total number of assigned edges
+
+    @assert sum_open_edges != 0 "Cannot have empty tree"
+    @assert sum_open_edges % 2 == 0 "Cannot have uneven number of open edges"
+
+    # find first entry that has only one open edge
+    single_open_edge_nodes = Vector{Int64}()
+    i = 0
+    for node in vec
+        i += 1
+        if node[1] - length(node[2]) == 1
+            push!(single_open_edge_nodes, i)
+        end
+    end
+
+    @assert length(single_open_edge_nodes) >= 2 "There were less than two nodes with only one open edge, which cannot happen for valid trees"
+
+    results = Vector{Vector{Vector{Int64}}}()
+
+    if sum_open_edges == 2
+        # assign open edges
+        n1 = single_open_edge_nodes[1]
+        n2 = single_open_edge_nodes[2]
+        push!(vec[n1][2], n2)
+        push!(vec[n2][2], n1)
+
+        push!(results, getindex.(vec, 2))
+
+        return results
+    end
+
+    # choose first node
+    n1 = single_open_edge_nodes[1]
+
+    # iterate through second nodes: all nodes that have at least 
+    for n2 in 1:n
+        if vec[n2][1] - length(vec[n2][2]) <= 1
+            # must have at least 2 open slots to be a partner for the edge
+            continue
+        end
+
+        # make new vec and connect edges
+        recurse_vec = deepcopy(vec)
+        push!(recurse_vec[n1][2], n2)
+        push!(recurse_vec[n2][2], n1)
+
+        generated_trees = _gen_trees(recurse_vec)
+        append!(results, generated_trees)
+    end
+
+    return results
+end
+
+function _labelled_plane_trees_unpermuted(N::Int64)
+    all_trees = Vector{Vector{Vector{Int64}}}()
+
+    parts = partitions(N - 2)
+    if N == 1
+        push!(all_trees, Vector{Vector{Vector{Int64}}}())
+        return all_trees
+    end
+
+    if N == 2
+        # fix empty partition when n = 2
+        parts = [Vector{Int64}()]
+    end
+
+    for origp in parts
+        p = copy(origp) # can't change original partition because of the iteration
+
+        n = length(p)
+        resize!(p, N)
+        for i in 1:N
+            if i > n
+                p[i] = 0
+            end
+            p[i] += 1
+        end
+
+        for mp in multiset_permutations(p, N)
+            # representation as vector of vectors
+            # each vector is a node
+            # each vector has the indices of the connected nodes
+            tree = [(l, Vector{Int64}()) for l in mp]
+
+            generated_trees = _gen_trees(tree)
+            append!(all_trees, generated_trees)
+        end
+    end
+
+    return all_trees
+end
+
+function _plane_tree_permutations_helper(t::Vector{Vector{Int64}}, depth::Int64)
+    if depth > length(t)
+        return [t]
+    end
+
+    all_perms = Vector{Vector{Vector{Int64}}}()
+    for permutation in permutations(t[depth])
+        permuted_t = deepcopy(t)
+        permuted_t[depth] = permutation
+        append!(all_perms, _plane_tree_permutations_helper(permuted_t, depth + 1))
+    end
+
+    return all_perms
+end
+
+# return all the possible edge permutations of a given labelled plane tree
+function _plane_tree_permutations(t::Vector{Vector{Int64}})
+    return _plane_tree_permutations_helper(t, 1)
+end
+
+function labelled_plane_trees(N::Int64)
+    all_trees = Vector{Vector{Vector{Int64}}}()
+
+    sizehint!(all_trees, factorial(3N - 3, 2N - 1))
+
+    for t in _labelled_plane_trees_unpermuted(N)
+        append!(all_trees, _plane_tree_permutations(t))
+    end
+
+    return all_trees
+end