@testset "NLP Solver" begin

  using PlotlyJS: savefig

  println("Testing NLP solver")

  # Test the optimizer for a one-phase mission
  # The lambert's solver said this should be pretty valid
  launch_window = DateTime(1992,11,1), DateTime(1992,12,1)
  latest_arrival = DateTime(1993,6,1)
  leave, arrive = DateTime(1992,11,19), DateTime(1993,4,1)
  test_leave = DateTime(1992,11,12)
  earth_state = state(Earth, leave)
  venus_state = state(Venus, arrive)
  v∞_out, v∞_in, tof = Thesis.lamberts(Earth, Venus, leave, arrive)
  # We can get the thrust profile and tof pretty wrong and still be ok
  phase = Phase(Venus, 1.1v∞_in, v∞_in, 0.9*tof, 0.1*ones(20,3))
  guess = Mission_Guess(bepi, 3_600., test_leave, 0.9*v∞_out, [phase])
  m = solve_mission(guess, launch_window, latest_arrival, verbose=true)
  @test typeof(m) == Mission

  # Now we can plot the results to check visually
  p = plot(m, title="NLP Test Solution")
  savefig(p,"../plots/nlp_test_1_phase.html")
  store(m, "missions/nlp_1_phase")

  # Now we can look at a slightly more complicated trajectory
  flybys = [Earth, Venus, Mars]
  launch_window = DateTime(2021,10,1), DateTime(2021,12,1)
  latest_arrival = DateTime(2023,1,1)
  dates = [DateTime(2021,11,1), DateTime(2022,3,27), DateTime(2022,8,28)]
  phases = Vector{Phase}()
  launch_v∞, _, tof1 = Thesis.lamberts(flybys[1], flybys[2], dates[1], dates[2])
  for i in 1:length(dates)-2
    v∞_out1, v∞_in1, tof1 = Thesis.lamberts(flybys[i], flybys[i+1], dates[i], dates[i+1])
    v∞_out2, v∞_in2, tof2 = Thesis.lamberts(flybys[i+1], flybys[i+2], dates[i+1], dates[i+2])
    push!(phases, Phase(flybys[i+1], v∞_in1, v∞_out2, tof1, 0.01*ones(20,3)))
  end
  v∞_out, v∞_in, tof = Thesis.lamberts(flybys[end-1], flybys[end], dates[end-1], dates[end])
  push!(phases, Phase(flybys[end], v∞_in, v∞_in, tof, 0.01*ones(20,3)))
  guess = Mission_Guess(bepi, 3_600., dates[1], launch_v∞, phases)
  m = solve_mission(guess, launch_window, latest_arrival, verbose=true)
  @test typeof(m) == Mission
  p = plot(m, title="NLP Test Solution (2 Phases)")
  savefig(p,"../plots/nlp_test_2_phase.html")
  store(m, "missions/nlp_2_phase")

  # Here is the final, most complicated, trajectory to test
  flybys = [Earth, Venus, Earth, Mars, Earth, Jupiter]
  launch_window = DateTime(2023,1,1), DateTime(2024,1,1)
  latest_arrival = DateTime(2031,1,1)
  dates = [DateTime(2023,5,23),
           DateTime(2023,10,21),
           DateTime(2024,8,24),
           DateTime(2025,2,13),
           DateTime(2026,11,22),
           DateTime(2032,1,1)]
  phases = Vector{Phase}()
  launch_v∞, _, tof1 = Thesis.lamberts(flybys[1], flybys[2], dates[1], dates[2])
  for i in 1:length(dates)-2
    v∞_out1, v∞_in1, tof1 = Thesis.lamberts(flybys[i], flybys[i+1], dates[i], dates[i+1])
    v∞_out2, v∞_in2, tof2 = Thesis.lamberts(flybys[i+1], flybys[i+2], dates[i+1], dates[i+2])
    push!(phases, Phase(flybys[i+1], 1.02v∞_in1, 0.98v∞_out2, 1.02tof1, 0.02*ones(20,3)))
  end
  v∞_out, v∞_in, tof = Thesis.lamberts(flybys[end-1], flybys[end], dates[end-1], dates[end])
  push!(phases, Phase(flybys[end], v∞_in, v∞_in, tof, 0.01*ones(20,3)))
  guess = Mission_Guess(bepi, 3_600., dates[1], launch_v∞, phases)
  m = solve_mission(guess, launch_window, latest_arrival, verbose=true)
  @test typeof(m) == Mission
  p = plot(m, title="NLP Test Solution (5 Phases)")
  savefig(p,"../plots/nlp_test_5_phase.html")
  store(m, "missions/nlp_5_phase")

end