diff --git a/julia/plotting.jl b/julia/plotting.jl
index 055c383..3662953 100644
--- a/julia/plotting.jl
+++ b/julia/plotting.jl
@@ -18,7 +18,8 @@ function plot_orbits(paths::Vector{Array{Float64,2}};
                      primary::String="Earth",
                      plot_theme::Symbol=:juno,
                      labels::Vector{String}=Vector{String}(),
-                     title::String="Spacecraft Position")
+                     title::String="Spacecraft Position",
+                     colors::Vector{String}=Vector{String}())
 
   N = 32
   θ = collect(range(0,length=N,stop=2π))
@@ -33,7 +34,7 @@ function plot_orbits(paths::Vector{Array{Float64,2}};
   for i = 1:length(paths)
     path = [ x for x in paths[i] ]
     label = labels != [] ? labels[i] : "orbit"
-    color = random_color()
+    color = colors != [] ? colors[i] : random_color()
     push!(t1, scatter3d(;x=(path[:,1]),y=(path[:,2]),z=(path[:,3]),
                          mode="lines", name=label, line_color=color, line_width=3))
   end
diff --git a/julia/propagator.jl b/julia/propagator.jl
index 161d562..b6e21de 100644
--- a/julia/propagator.jl
+++ b/julia/propagator.jl
@@ -26,13 +26,13 @@ function prop_one(ΔV::Vector{T},
 
   mag, α, β = ΔV
 
-  if mag > 1 || mag < 0
-    throw(ErrorException("ΔV input is too high: $mag"))
-  elseif α > π || α < -π
-    throw(ErrorException("α angle is incorrect: $α"))
-  elseif β > π/2 || β < -π/2
-    throw(ErrorException("β angle is incorrect: $β"))
-  end
+  # if mag > 1 || mag < 0
+    # throw(ErrorException("ΔV input is too high: $mag"))
+  # elseif α > π || α < -π
+    # throw(ErrorException("α angle is incorrect: $α"))
+  # elseif β > π/2 || β < -π/2
+    # throw(ErrorException("β angle is incorrect: $β"))
+  # end
 
   thrust_rθh = mag * [cos(β)*sin(α), cos(β)*cos(α), sin(β)]
   a,e,i,Ω,ω,ν = xyz_to_oe(state, μ)
diff --git a/julia/single_shoot.jl b/julia/single_shoot.jl
index 983f732..5413188 100644
--- a/julia/single_shoot.jl
+++ b/julia/single_shoot.jl
@@ -1,13 +1,8 @@
-using NLsolve
+using NLsolve, NLopt
 
-function treat_inputs(x::AbstractVector, n::Int)
-  inputs = reshape(copy(x),(3,n))'
-  for i in 1:n
-    inputs[i,1] = 0.5*tanh(inputs[i,1]) + 0.5
-    inputs[i,2] = π*tanh(inputs[i,2])
-    inputs[i,3] = π*tanh(inputs[i,3])/2
-  end
-  return inputs
+function treat_inputs(x::AbstractVector)
+  n::Int = length(x)/3
+  reshape(x,(3,n))'
 end
 
 function single_shoot(start::Vector{Float64},
@@ -21,13 +16,72 @@ function single_shoot(start::Vector{Float64},
                       tol=1e-2)
 
   function f!(F,x)
-    F[1:6] .= prop(treat_inputs(x,n), start, craft, μ, tf-t0)[1][end,:] - final
+    F[1:6] .= prop(treat_inputs(x), start, craft, μ, tf-t0)[1][end,:] - final
     F[7:3n] .= 0.
-    # if typeof(F[1]) == Float64 println(F[1:6]) end
-    # if typeof(F[1]) == Float64 println(treat_inputs(x,n)[1:8,1]) end
-    
   end
 
   return nlsolve(f!, x0, ftol=tol, autodiff=:forward, iterations=10_000)
 
 end
+
+function single_shoot2(start::Vector,
+                       final::Vector,
+                       craft::Sc,
+                       μ::AbstractFloat,
+                       t0::AbstractFloat,
+                       tf::AbstractFloat,
+                       x0::Vector,
+                       tol=1e-8)
+
+  n::Int = length(x0)/3
+  m0 = craft.mass
+
+  f(x::Vector) = m0 - prop(treat_inputs(x), start, craft, μ, tf-t0)[2][end]
+  f_constraint(x::Vector) = norm(prop(treat_inputs(x), start, craft, μ, tf-t0)[1][end,:] - final)
+
+  function nlfunc(x::Vector, grad::Vector)
+    try
+      if length(grad) != 0
+        ForwardDiff.gradient!(grad, f, x)
+      end
+      f(x)
+    catch e
+      println("Error was $e")
+      throw(e)
+    end
+  end
+
+  function nlconstraint(x::Vector, grad::Vector)
+    if length(grad) != 0
+      ForwardDiff.gradient!(grad, f_constraint, x)
+    end
+    f_constraint(x)
+  end
+
+  opt = Opt(:LD_MMA, 3n)
+  lower_bounds = Vector{Float64}()
+  upper_bounds = Vector{Float64}()
+  for i in 1:3n
+    if i%3 == 1
+      push!(lower_bounds, 0.)
+      push!(upper_bounds, 1.)
+    elseif i%3 == 2
+      push!(lower_bounds, -π)
+      push!(upper_bounds, π)
+    elseif i%3 == 0
+      push!(lower_bounds, -π/2)
+      push!(upper_bounds, π/2)
+    end
+  end
+  opt.lower_bounds = lower_bounds
+  opt.upper_bounds = upper_bounds
+  opt.xtol_rel = 1e-4
+  opt.min_objective = nlfunc
+  inequality_constraint!(opt, nlconstraint, 1e-8)
+
+  (minf, minx, ret) = optimize(opt, x0)
+  numevals = opt.numevals
+
+  return minf, minx, ret, numevals
+
+end
diff --git a/julia/test/propagator.jl b/julia/test/propagator.jl
index 51b7e34..13b8ab1 100644
--- a/julia/test/propagator.jl
+++ b/julia/test/propagator.jl
@@ -27,9 +27,9 @@
   @test craft.mass == start_mass - craft.mass_flow_rate*stepsize
   
   # Test that a bad ΔV throws an error
-  craft = Sc("test")
-  start_mass = craft.mass
-  @test_throws ErrorException prop_one([1.5, 0., 0.], start, craft, μs["Earth"], stepsize)
+  # craft = Sc("test")
+  # start_mass = craft.mass
+  # @test_throws ErrorException prop_one([1.5, 0., 0.], start, craft, μs["Earth"], stepsize)
 
   # Test that a full propagation doesn't take too long
 
diff --git a/julia/test/single_shoot.jl b/julia/test/single_shoot.jl
index 7194635..96c36dc 100644
--- a/julia/test/single_shoot.jl
+++ b/julia/test/single_shoot.jl
@@ -6,26 +6,31 @@
   e = rand(0.01:0.01:0.5)
   i = rand(0.01:0.01:π/6)
   T = 2π*√(a^3/μs["Earth"])
+  prop_time = 2T
   n = 50
 
   # A simple orbit raising
   start = oe_to_xyz([ a, e, i, 0., 0., 0. ], μs["Earth"])
   ΔVs = repeat([0.6, 0., 0.]', outer=(n,1))
-  final = prop(ΔVs, start, sc, μs["Earth"], T)[1][end,:]
+  final = prop(ΔVs, start, sc, μs["Earth"], prop_time)[1][end,:]
+  new_T = 2π*√(xyz_to_oe(final, μs["Earth"])[1]^3/μs["Earth"])
 
   # This should be close enough to 0.6
-  x0 = repeat([atanh((0.4-0.5)/0.5), 0., 0.], n)
-  result = single_shoot(start, final, sc, μs["Earth"], 0.0, T, n, x0)
+  x0 = repeat([0.59, 0., 0.], n)
+  result = single_shoot2(start, final, sc, μs["Earth"], 0.0, prop_time, x0)
 
   # Test and plot
-  @test converged(result)
+  @test result[3] == :XTOL_REACHED
   path1 = prop(zeros((100,3)), start, sc, μs["Earth"], T)[1]
-  path2, mass = prop(treat_inputs(result.zero, n), start, sc, μs["Earth"], T)
-  path3 = prop(zeros((100,3)), path2[end,:], sc, μs["Earth"], T)[1]
-  savefig(plot_orbits([path1, path2, path3]), "single_shoot_test.html")
-  if converged(result)
-    @test norm(path2[end,:] - final) < 2e-2
-    sc = Sc("test")
+  path2, mass = prop(treat_inputs(result[2]), start, sc, μs["Earth"], prop_time)
+  path3 = prop(zeros((100,3)), path2[end,:], sc, μs["Earth"], new_T)[1]
+  path4 = prop(zeros((100,3)), final, sc, μs["Earth"], new_T)[1]
+  savefig(plot_orbits([path1, path2, path3, path4],
+                      labels=["inital", "transit", "after transit", "final"],
+                      colors=["#FFFFFF","#FF4444","#44FF44","#4444FF"]),
+          "single_shoot_test.html")
+  if result[3] == :XTOL_REACHED
+    @test norm(path2[end,:] - final) < 1e-6
   end
 
 end