diff --git a/AutFN_orbit.jl b/AutFN_orbit.jl
new file mode 100644
index 0000000..b820d31
--- /dev/null
+++ b/AutFN_orbit.jl
@@ -0,0 +1,108 @@
+include("Orb_AutFN.jl")
+
+###############################################################################
+#
+#  Action of WreathProductElems on AutGroupElem
+#
+###############################################################################
+
+function AutFG_emb(A::AutGroup, g::WreathProductElem)
+    isa(A.objectGroup, FreeGroup) || throw("Not an Aut(Fₙ)")
+    parent(g).P.n == length(A.objectGroup.gens) || throw("No natural embedding of $(parent(g)) into $A")
+    powers = [(elt == parent(elt)() ? 0: 1) for elt in g.n.elts]
+    elt = reduce(*, [A(Groups.flip_autsymbol(i))^pow for (i,pow) in enumerate(powers)])
+    Groups.r_multiply!(elt, [Groups.perm_autsymbol(g.p)])
+    return elt
+end
+
+function (g::WreathProductElem)(a::AutGroupElem)
+    g = AutFG_emb(parent(a),g)
+    return g*a*g^-1
+end
+
+###############################################################################
+#
+#  Generating set
+#
+###############################################################################
+
+function SOutFN_generating_set(N::Int)
+
+   SOutFN = AutGroup(FreeGroup(N), special=true, outer=true)
+   S = generators(SOutFN);
+   S = [S; [inv(s) for s in S]]
+
+   return SOutFN, unique(S)
+end
+
+###############################################################################
+#
+#  Parsing command line
+#
+###############################################################################
+
+function parse_commandline()
+    settings = ArgParseSettings()
+
+    @add_arg_table settings begin
+        "--tol"
+            help = "set numerical tolerance for the SDP solver (default: 1e-5)"
+            arg_type = Float64
+            default = 1e-5
+        "--iterations"
+            help = "set maximal number of iterations for the SDP solver (default: 20000)"
+            arg_type = Int
+            default = 20000
+        "--upper-bound"
+            help = "Set an upper bound for the spectral gap (default: Inf)"
+            arg_type = Float64
+            default = Inf
+        "--cpus"
+            help = "Set number of cpus used by solver (default: auto)"
+            arg_type = Int
+            required = false
+        "-N"
+            help = "Consider automorphisms of free group on N generators (default: N=2)"
+            arg_type = Int
+            default = 2
+        "--radius"
+            help = "Find the decomposition over B_r(e,S)"
+            arg_type = Int
+            default = 2
+    end
+
+    return parse_args(settings)
+end
+
+###############################################################################
+#
+#  main
+#
+###############################################################################
+
+function main()
+   parsed_args = parse_commandline()
+
+   N = parsed_args["N"]
+   dirname = "SOutF$(N)_E_r=$radius"
+   radius = parsed_args["radius"]
+
+   isdir(dirname) || mkdir(dirname)
+   logger = PropertyT.setup_logging(dirname)
+
+   G, S = SOutFN_generating_set(logger, dirname, N)
+   info(logger, G)
+   info(logger, "Symmetric generating set of size $(length(S))")
+   info(logger, S)
+   AutS = WreathProduct(Nemo.FiniteField(2,1, "a")[1], PermutationGroup(N))
+
+   tol = parsed_args["tol"]
+   iterations = parsed_args["iterations"]
+   upper_bound = parsed_args["upper-bound"]
+
+   solver = SCS.SCSSolver(eps=tol, max_iters=iterations, verbose=true, linearsolver=SCS.Indirect)
+
+   orbit_check_propertyT(logger, dirname, G, S, AutS, solver, upper_bound)
+end
+
+main()