@testset "Phase" begin

  println("Testing NLP solver")

  using SNOW, PlotlyJS

  # First we'll test an Earth case
  # Initial Setup
  T = rand( 8hour : second : day )
  revs = 8
  n = 10 * revs
  start_mass = 10_000.

  # A simple orbit raising
  start = gen_orbit(T, start_mass, Earth)
  thrust = spiral(0.9, n, start, test_sc, revs*T, Earth)
  final = prop(thrust, start, test_sc, revs*T, Earth)[2]
  new_T = 2π * √( xyz_to_oe(final, Earth.μ)[1]^3 / Earth.μ )

  # This should be close enough to converge, but Earth orbits are picky
  thrust_guess = spiral(0.89, n, start, test_sc, revs*T, Earth)
  guess_final = prop(thrust_guess, start, test_sc, revs*T, Earth)[2]
  result = solve_phase(start, final, test_sc, revs*T, thrust_guess, Earth, verbose=true)
  calc_path, calc_final = prop(result.sol, start, test_sc, revs*T, Earth, interpolate=true)

  # Test
  @test norm(guess_final[1:6] - final[1:6]) > 1e-4
  @test result.converged
  @test norm(calc_final[1:6] - final[1:6]) < 1e-4

  # Plot
  paths = Pathlist()
  push!(paths, prop(start, T, Earth),
               calc_path,
               prop(calc_final, T, Earth),
               prop(final, T, Earth) )
  fig = plot_orbits(paths, Earth,
                    labels=["init", "transit", "post-transit", "final"],
                    colors=["#FFF","#F44","#4F4","#44F"])
  savefig(fig, "../plots/nlp_test_earth.html")

  # Now the sun case
  T = rand( 0.5year : hour : 1.5year )
  tof = 0.7T
  n = 20
  start_mass = 12_000.

  # A simple orbit raising again, to keep things simple
  start = gen_orbit(T, start_mass)
  thrust = spiral(0.6, n, start, bepi, tof)
  final = prop(thrust, start, bepi, tof)[2]
  new_T = 2π * √( xyz_to_oe(final, Sun.μ)[1]^3 / Sun.μ )

  # This should be close enough to converge
  thrust_guess = spiral(0.55, n, start, bepi, tof)
  guess_final = prop(thrust_guess, start, bepi, tof)[2]
  result = solve_phase(start, final, bepi, tof, thrust_guess, verbose=true)
  calc_path, calc_final = prop(result.sol, start, bepi, tof, interpolate=true)

  # Test
  @test norm(guess_final[1:6] - final[1:6]) > 1e-4
  @test result.converged
  @test norm(calc_final[1:6] - final[1:6]) < 1e-4

  # Plot
  paths = Pathlist()
  push!(paths, prop(start, T), calc_path, prop(calc_final, T), prop(final, T) )
  fig = plot_orbits(paths,
                    labels=["init", "transit", "post-transit", "final"],
                    colors=["#FFF","#F44","#4F4","#44F"])
  savefig(fig, "../plots/nlp_test_sun.html")


end