NLP seems to be doing pretty well

2021-10-10 20:53:14 -06:00
parent f83d73449f
commit b29afbce04
13 changed files with 524 additions and 198 deletions
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -12,8 +12,13 @@ unit-test-job:
    - julia --project=julia/Project.toml -E 'using Pkg; Pkg.test()'
  artifacts:
    paths:
-      - julia/plots/plot_test_earth.html
-      - julia/plots/plot_test_sun.html
-      - julia/plots/nlp_test_earth.html
-      - julia/plots/nlp_test_sun.html
+      - julia/plots/nlp_test_1_phase.html
+      - julia/plots/nlp_test_2_phase.html
+      - julia/plots/plot_test_mission.html
+      - julia/plots/plot_test_mission_guess.html
+      - julia/plots/plot_test_path.html
+      - julia/plots/plot_test_paths.html
+      - julia/plots/plot_test_planet_move.html
+      - julia/plots/plot_test_planet_no_move.html
+      - julia/plots/plot_test_planet_solar_system.html
    expire_in: 10 weeks
--- a/julia/src/Thesis.jl
+++ b/julia/src/Thesis.jl
@@ -19,14 +19,14 @@ module Thesis
  include("./bodies.jl")
  include("./constants.jl")
  include("./spacecraft.jl")
-  include("./mission.jl")
-  include("./conversions.jl")
-  include("./lamberts.jl")
-  include("./plotting.jl")
  include("./inner_loop/laguerre-conway.jl")
  include("./inner_loop/propagator.jl")
+  include("./mission.jl")
+  include("./conversions.jl")
+  include("./plotting.jl")
  include("./inner_loop/nlp_solver.jl")
  # include("./inner_loop/monotonic_basin_hopping.jl")
  # include("./outer_loop.jl")
+  include("./lamberts.jl")

 end
--- a/julia/src/bodies.jl
+++ b/julia/src/bodies.jl
@@ -1,17 +1,22 @@
-export state
+export state, period

 struct Body
  μ::Float64
  r::Float64 # radius
  color::String
  id::Int # SPICE id
+  a::Float64 # semi-major axis
+  line_color::AbstractString
+  name::AbstractString
 end

-function state(p::Body, t::DateTime, add_v∞::Vector{Float64}=[0., 0., 0.], add_mass::Float64=1e10)
+period(b::Body) = 2π * √(b.a^3 / Sun.μ)
+
+function state(p::Body, t::DateTime, add_v∞::Vector{T}=[0., 0., 0.], add_mass::Float64=1e10) where T <: Real
  time = utc2et(Dates.format(t,"yyyy-mm-ddTHH:MM:SS"))
  [ spkssb(p.id, time, "ECLIPJ2000"); 0.0 ] + [ zeros(3); add_v∞; add_mass ]
 end

-function state(p::Body, t::Float64, add_v∞::Vector{Float64}=[0., 0., 0.], add_mass::Float64=1e10)
+function state(p::Body, t::S, add_v∞::Vector{T}=[0., 0., 0.], add_mass::Float64=1e10) where {T,S <: Real}
  [ spkssb(p.id, t, "ECLIPJ2000"); 0.0 ] + [ zeros(3); add_v∞; add_mass ]
 end
--- a/julia/src/constants.jl
+++ b/julia/src/constants.jl
@@ -7,17 +7,16 @@ export Jupiter, Saturn, Uranus, Neptune, Pluto
 export G, AU, init_STM, hour, day, year, second
 export Pathlist

-const Sun = Body(1.32712440018e11, 696000., "Electric", 10)
-const Mercury = Body(2.2032e4, 2439., "heat", 1)
-const Venus = Body(3.257e5, 6052., "turbid", 2)
-const Earth = Body(3.986004415e5, 6378.1363, "Blues", 399)
-const Moon = Body(4.902799e3, 1738., "Greys", 301)
-const Mars = Body(4.305e4, 3397.2, "Reds", 4)
-const Jupiter = Body(1.266865361e8, 71492., "solar", 5)
-const Saturn = Body(3.794e7, 60268., "turbid", 6)
-const Uranus = Body(5.794e6, 25559., "haline", 7)
-const Neptune = Body(6.809e6, 24764., "ice", 8)
-const Pluto = Body(9e2, 1151., "matter", 9)
+const Sun = Body(1.32712440018e11, 696000., "YlOrRd", 10, 0., "", "Sun")
+const Mercury = Body(2.2032e4, 2439., "YlOrRd", 1, 57909083., "#FF2299", "Mercury")
+const Venus = Body(3.257e5, 6052., "Portland", 2, 108208601., "#FFCC00", "Venus")
+const Earth = Body(3.986004415e5, 6378.1363, "Earth", 399, 149598023., "#0055FF", "Earth")
+const Mars = Body(4.305e4, 3397.2, "Reds", 4, 227939186., "#FF0000", "Mars")
+const Jupiter = Body(1.266865361e8, 71492., "YlOrRd", 5, 778298361., "#FF7700", "Jupiter")
+const Saturn = Body(3.794e7, 60268., "Portland", 6, 1429394133., "#FFFF00", "Saturn")
+const Uranus = Body(5.794e6, 25559., "YlGnBu", 7, 2875038615., "#0000FF", "Uranus")
+const Neptune = Body(6.809e6, 24764., "Blues", 8, 4504449769., "#6666FF", "Neptune")
+const Pluto = Body(9e2, 1151., "Picnic", 9, 5915799000., "#FFCCCC", "Pluto")

 const G = 6.67430e-20 #universal gravity parameter
 const AU = 149597870.691 #km
--- a/julia/src/errors.jl
+++ b/julia/src/errors.jl
@@ -15,7 +15,10 @@ struct V∞_Error <: Exception
  v∞_in::AbstractVector
  v∞_out::AbstractVector
 end
-Base.showerror(io::IO, e::V∞_Error) = print(io, "v∞s weren't the same: ", e.v∞_in, " and ", e.v∞_out)
+Base.showerror(io::IO, e::V∞_Error) = print(io, """v∞s weren't the same:
+                                                   v1 = $(e.v∞_in) and
+                                                   v2 = $(e.v∞_out)
+                                                   difference: $(norm(e.v∞_in) - norm(e.v∞_out))""")

 struct PropOne_Error <: Exception
  ΔV_unit::AbstractVector
@@ -31,3 +34,13 @@ struct Mission_Creation_Error{T} <: Exception where T <: Real
  num_entries::T
 end
 Base.showerror(io::IO, e::Mission_Creation_Error) = print(io, "Tried to create a mission with $(e.num_entries) entries")
+
+struct Planet_Match_Error{T} <: Exception where T <: Real
+  state::Vector{T}
+  planet::Vector{T}
+end
+Base.showerror(io::IO, e::Planet_Match_Error) = print(io, """"
+    Spacecraft didn't make it to planet
+    spacecraft_pos = $(e.state)
+    planet_pos = $(e.planet)
+    difference = $(abs(norm(e.state - e.planet)))""")
--- a/julia/src/inner_loop/nlp_solver.jl
+++ b/julia/src/inner_loop/nlp_solver.jl
@@ -5,19 +5,36 @@ export solve_mission
 struct Result
  converged::Bool
  info::Symbol
-  sol::Matrix{Float64}
+  sol::Mission_Guess
 end

+function constraint_bounds(guess::Mission_Guess)
+  low_constraint = Vector{Float64}()
+  high_constraint = Vector{Float64}()
+
+  for phase in guess.phases
+    if phase != guess.phases[end]
+      push!(low_constraint, -100., -100., -100., -0.005, -0.005, -0.005, -30., 100.)
+      push!(high_constraint, 100., 100., 100., 0.005, 0.005, 0.005, 30., 100_000.)
+    else
+      push!(low_constraint, -100., -100., -100., 0.)
+      push!(high_constraint, 100., 100., 100., guess.start_mass - guess.sc.dry_mass)
+    end
+  end
+  return low_constraint, high_constraint
+
+end
 """
 This function will take an initial guess to a mission profile and leverage Ipopt to solve the mission
+
+NOTE: This function will output a proper mission if it succeeded and just re-return the original
+guess otherwise
 """
-function solve_mission( guess::Mission_Guess;
-                        tol=1e-8,
+function solve_mission( guess::Mission_Guess,
+                        launch_window::Vector{DateTime},
+                        latest_arrival::DateTime;
+                        tol=1e-10,
                        verbose::Bool=false )
-  # First we need to define a function that propagates the mission all the way through
-  # The function takes in a constraints vector (g) and an input vector (x) and spits out 
-  # the objective function. In my case, the objective will be constant.
-  #

  # First we define our starting point
  x0 = Vector(guess)
@@ -52,59 +69,39 @@ function solve_mission( guess::Mission_Guess;
      goal = state(current_planet, time, v∞_in)
      final = prop(reshape(thrusts, (n,3)), start, guess.sc, tof)[2]
      if flyby != flybys[end]
-        g[1+8i:6+8i] .= final[1:6] .- goal[1:6]
-        g[7+8i] = norm(v∞_out) - norm(v∞_in)
+        g[1+8i:6+8i] .= (final[1:6] .- goal[1:6]) .* [1., 1., 1., 10_000., 10_000., 10_000.]
+        g[7+8i] = (norm(v∞_out) - norm(v∞_in)) * 10_000.
        δ = acos( ( v∞_in ⋅ v∞_out ) / ( norm(v∞_in) * norm(v∞_out) ) )
-        g[8+8i] = (current_planet.μ/(v∞_in ⋅ v∞_in)) * ( 1/sin(δ/2) - 1 )
+        g[8+8i] = (current_planet.μ/(v∞_in ⋅ v∞_in)) * ( 1/sin(δ/2) - 1 ) - current_planet.r
      else
-        g[8length(flybys)+1:8length(flybys)+3] .= final[1:3] .- goal[1:3]
-        g[8length(flybys)+4] = final[7]
+        g[8*(length(flybys)-1)+1:8*(length(flybys)-1)+3] .= final[1:3] .- goal[1:3]
+        g[8*(length(flybys)-1)+4] = final[7] - guess.sc.dry_mass
      end
      start = state(current_planet, time, v∞_out, final[7])
      i += 1
-      return 1.0
    end
+    return 1.0
  end
-  g = zeros(8length(flybys)+4)
-  optimizer!(g,x0)
-  return g

+  max_time = Dates.datetime2unix(latest_arrival) - Dates.datetime2unix(launch_window[1])
+  lower_x = lowest_mission_vector(launch_window, length(guess.phases), n)
+  upper_x = highest_mission_vector(launch_window, max_time, length(guess.phases), n)
+  num_constraints = 8*(length(guess.phases)-1) + 4
+  g_low, g_high = constraint_bounds(guess)
+  ipopt_options = Dict("constr_viol_tol" => tol,
+                       "acceptable_constr_viol_tol" => 100tol,
+                       "max_iter" => 10_000,
+                       "max_cpu_time" => 300.,
+                       "print_level" => 5)
+  options = Options(solver=IPOPT(ipopt_options), derivatives=ForwardFD())

+  x, _, info = minimize(optimizer!, x0, num_constraints, lower_x, upper_x, g_low, g_high, options)
+  if verbose println(info) end

-  # n = size(x0)[1]
-
-  # function f!(g,x)
-    # try
-      # g[1:6] .= prop(reshape(x,(n,3)), start, craft, tof, primary)[2][1:6]
-      # return 1.
-    # catch e
-      # if isa(e,LaGuerreConway_Error) || isa(e, Mass_Error) || isa(e, PropOne_Error)
-        # return 10_000.
-      # else
-        # rethrow
-      # end
-    # end
-  # end
-
-  # lower_x = -1 * ones(3n)
-  # upper_x = ones(3n)
-  # num_constraints = 6
-  # filename = "ipopt_"*string(rand(UInt))
-  # ipopt_options = Dict("constr_viol_tol" => tol,
-                       # "acceptable_constr_viol_tol" => 1e-4,
-                       # "max_iter" => 10_000,
-                       # "max_cpu_time" => 30.,
-                       # "print_level" => 0,
-                       # "output_file" => filename)
-  # options = Options(solver=IPOPT(ipopt_options),
-                    # derivatives=ForwardAD())
-  # x, _, info = minimize(f!, vec(x0), num_constraints, lower_x, upper_x, final[1:6], final[1:6], options)
-  # rm(filename)
-  # if info in [:Solve_Succeeded, :Solved_To_Acceptable_Level]
-    # return Result(true, info, reshape(x,(n,3)))
-  # else
-    # if verbose println(info) end
-    # return Result(false, info, x0)
-  # end
+  if info in [:Solve_Succeeded, :Solved_To_Acceptable_Level]
+    return Mission(x, guess.sc, guess.start_mass, flybys)
+  else
+    return guess
+  end

 end
--- a/julia/src/inner_loop/propagator.jl
+++ b/julia/src/inner_loop/propagator.jl
@@ -56,7 +56,20 @@ function prop(ΔVs::Matrix{T},
        for k in 1:7 push!(states[k], interpolated_state[k]) end
      end
    end
-    state = prop_one(ΔVs[i,:], state, craft, time/n, primary)
+    try
+      state = prop_one(ΔVs[i,:], state, craft, time/n, primary)
+    catch e
+      if isa(e,PropOne_Error)
+        for val in e.ΔV_unit
+          # If this isn't true, then we just let it slide
+          if abs(val) > 1.0001
+            rethrow()
+          end
+        end
+      else
+        rethrow()
+      end
+    end
    for j in 1:7 push!(states[j], state[j]) end
    state[7] >= craft.dry_mass || throw(Mass_Error(state[7]))
  end
@@ -68,4 +81,4 @@ end
 """
 Convenience function for propagating a state with no thrust
 """
-prop(x::Vector{Float64}, t::Float64, p::Body=Sun) = prop(zeros(1000,3), [x;1e10], no_thrust, t, p)[1]
+prop(x::Vector{Float64}, t::Float64, p::Body=Sun) = prop(zeros(1000,3), [x;1.], no_thrust, t, p)[1]
--- a/julia/src/mission.jl
+++ b/julia/src/mission.jl
@@ -1,6 +1,6 @@
 export Phase, Mission_Guess, Mission, Bad_Mission
 export test_phase1, test_phase2
-export test_mg, test_mission
+export test_mg
 export Vector

 struct Phase
@@ -11,7 +11,7 @@ struct Phase
  thrust_profile::Matrix{Float64}
 end

-const test_phase1 = Phase(Venus, [9., 10., 0.], [10., 9, 0.], 1.30464e7, zeros(20,3))
+const test_phase1 = Phase(Venus, [√20., -√80., 0.], [√50., -√50., 0.], 1.30464e7, zeros(20,3))
 const test_phase2 = Phase(Jupiter, [0.3, 7.1, 0.2], [0.3, 7.1, 0.2], 3.9year, zeros(20,3))
 const test_phases = [test_phase1, test_phase2]

@@ -46,12 +46,6 @@ function Mission_Guess(sc::Sc, mass::Float64, date::DateTime, v∞::Vector{Float
    mass_used += mass_consumption(sc, phase)
    mass - mass_used > sc.dry_mass || throw(Mass_Error(mass - mass_used))
    v∞_in, v∞_out = phase.v∞_in, phase.v∞_out
-    if phase != phases[end]
-      norm(v∞_in) ≈ norm(v∞_out) || throw(V∞_Error(v∞_in, v∞_out))
-      δ = acos( ( v∞_in ⋅ v∞_out ) / ( norm(v∞_in) * norm(v∞_out) ) )
-      periapsis = (phase.planet.μ/(v∞_in ⋅ v∞_in)) * ( 1/sin(δ/2) - 1 )
-      periapsis > 1.1phase.planet.r || throw(HitPlanet_Error())
-    end
  end
  Mission_Guess(sc, mass, date, v∞, phases, false)
 end
@@ -101,6 +95,32 @@ function Base.Vector(g::Mission_Guess)
  return result
 end

+function lowest_mission_vector(launch_window::Vector{DateTime}, num_phases::Int, n::Int)
+  result = Vector{Float64}()
+  push!(result, Dates.datetime2unix(launch_window[1]))
+  push!(result, -10*ones(3)...)
+  for i in 1:num_phases
+    push!(result, -10*ones(3)...)
+    push!(result, -10*ones(3)...)
+    push!(result, 86_400.0)
+    push!(result, -1 * ones(n,3)...)
+  end
+  return result
+end
+
+function highest_mission_vector(launch_window::Vector{DateTime}, mission_length::Float64, num_phases::Int, n::Int)
+  result = Vector{Float64}()
+  push!(result, Dates.datetime2unix(launch_window[2]))
+  push!(result, 10*ones(3)...)
+  for i in 1:num_phases
+    push!(result, 10*ones(3)...)
+    push!(result, 10*ones(3)...)
+    push!(result, mission_length)
+    push!(result, ones(n,3)...)
+  end
+  return result
+end
+
 const test_mg = Mission_Guess(bepi, 12_000., DateTime(1992,11,19), [-3.4,1.2,0.1], test_phases)
            
 struct Mission
@@ -113,27 +133,47 @@ struct Mission
 end

 """
-Constructor for a mission. Generally mission guesses are converged
+Constructor for a mission. Generally mission guesses are converged, so we check everything
 """
-function Mission(sc::Sc, mass::Float64, date::DateTime, v∞::Vector{Float64}, phases::Vector{Phase})
+function Mission(sc::Sc, mass::Float64, date::DateTime, v∞::Vector{Float64}, phases::Vector{Phase};
+                 force=false)
  # First do some checks to make sure that it's valid
-  mass_used = 0
-  for phase in phases
-    mass_used += mass_consumption(sc, phase)
-    mass - mass_used > sc.dry_mass || throw(Mass_Error(mass - mass_used))
-    v∞_in, v∞_out = phase.v∞_in, phase.v∞_out
-    if phase != phases[end]
-      norm(v∞_in) ≈ norm(v∞_out) || throw(V∞_Error(v∞_in, v∞_out))
-      δ = acos( ( v∞_in ⋅ v∞_out ) / ( norm(v∞_in) * norm(v∞_out) ) )
-      periapsis = (phase.planet.μ/(v∞_in ⋅ v∞_in)) * ( 1/sin(δ/2) - 1 )
-      periapsis > 1.1phase.planet.r || throw(HitPlanet_Error())
+  if !force
+    time = date
+    current_planet = Earth
+    start = state(current_planet, time, v∞, mass)
+    for phase in phases
+      final = prop(phase.thrust_profile, start, sc, phase.tof)[2]
+      mass = final[7]
+      mass > sc.dry_mass || throw(Mass_Error(mass - mass_used))
+      current_planet = phase.planet
+      time += Dates.Second(floor(phase.tof))
+      p_pos = state(current_planet, time)[1:3]
+      abs(norm(final[1:3] - p_pos)) < √30_000. || throw(Planet_Match_Error(final[1:3], p_pos))
+      v∞_in, v∞_out = phase.v∞_in, phase.v∞_out
+      if phase != phases[end]
+        abs(norm(v∞_in) - norm(v∞_out)) < 0.005 || throw(V∞_Error(v∞_in, v∞_out))
+        δ = acos( ( v∞_in ⋅ v∞_out ) / ( norm(v∞_in) * norm(v∞_out) ) )
+        periapsis = (phase.planet.μ/(v∞_in ⋅ v∞_in)) * ( 1/sin(δ/2) - 1 )
+        periapsis > 1.1phase.planet.r || throw(HitPlanet_Error())
+      end
    end
  end
  Mission(sc, mass, date, v∞, phases, true)
 end

-const test_mission = Mission(bepi, 12_000., DateTime(1992,11,19), [-3.4,1.2,0.1], test_phases)
- 
+"""
+BE CAREFUL!! This just makes a guess converged, whether true or not
+"""
+function Mission(g::Mission_Guess)
+  return Mission(g.sc, g.start_mass, g.launch_date, g.launch_v∞, g.phases, force=true)
+end
+
+function Mission(x::Vector{Float64}, sc::Sc, mass::Float64, flybys::Vector{Body})
+  guess = Mission_Guess(x, sc, mass, flybys::Vector{Body})
+  return Mission(guess)
+end
+
 struct Bad_Mission
  message::String
  converged::Bool
--- a/julia/src/plotting.jl
+++ b/julia/src/plotting.jl
@@ -3,7 +3,7 @@

 using Random, PlotlyJS, Base.Iterators

-export plot_orbits
+export plot

 function random_color()
  num1 = rand(0:255)
@@ -12,60 +12,281 @@ function random_color()
  return "#"*string(num1, base=16, pad=2)*string(num2, base=16, pad=2)*string(num3, base=16, pad=2)
 end

-function get_true_max(mat::Vector{Array{Float64,2}})
-  return maximum(abs.(flatten(mat)))
+# Some convenience functions for maximums
+"""
+Gets the maximum absolute value of a path
+"""
+function Base.maximum(path::Vector{Vector{Float64}})
+  return maximum( abs.( reduce( vcat,path ) ) )
 end

-function plot_orbits(paths::Vector{Vector{Vector{Float64}}},
-                     primary::Body=Sun;
-                     labels::Vector{String}=Vector{String}(),
-                     title::String="Spacecraft Position",
-                     colors::Vector{String}=Vector{String}())
+"""
+Gets the maximum absolute value of a list of paths
+"""
+function Base.maximum(paths::Vector{Vector{Vector{Float64}}})
+  return maximum( abs.( reduce( vcat,reduce( vcat,paths ) ) ) )
+end

+"""
+Gets the maximum absolute value of a surface
+"""
+function Base.maximum(xs::Matrix{Float64},ys::Matrix{Float64},zs::Matrix{Float64})
+  return maximum([ vec(xs), vec(ys), vec(zs) ])
+end
+
+"""
+Basically rewriting Julia at this point
+"""
+Base.getindex(x::Nothing, i::Int) = nothing
+
+"""
+Generates a layout that works for most plots. Requires a title and a limit.
+"""
+function standard_layout(limit::Float64, title::AbstractString)
+  limit *= 1.1
+  Layout(title=title,
+         # width=1400,
+         # height=800,
+         paper_bgcolor="rgba(5,10,40,1.0)",
+         plot_bgcolor="rgba(100,100,100,0.01)",
+         scene = attr(xaxis = attr(autorange = false,range=[-limit,limit]),
+                      yaxis = attr(autorange = false,range=[-limit,limit]),
+                      zaxis = attr(autorange = false,range=[-limit,limit]),
+         aspectratio=attr(x=1,y=1,z=1),
+         aspectmode="manual"))
+end
+
+"""
+A fundamental function to generate a plot for a single orbit path
+"""
+function gen_plot(path::Vector{Vector{Float64}};
+                  label::Union{AbstractString, Nothing} = nothing,
+                  color::AbstractString = random_color(),
+                  markers=true)
+
+  # Plot the orbit
+  if label === nothing
+    path_trace = scatter3d(;x=(path[1]),y=(path[2]),z=(path[3]),
+                            mode="lines", showlegend=false, line_color=color, line_width=3)
+  else
+    path_trace = scatter3d(;x=(path[1]),y=(path[2]),z=(path[3]),
+                            mode="lines", name=label, line_color=color, line_width=3)
+  end
+  traces = [ path_trace ]
+
+  # Plot the markers
+  if markers
+    push!(traces, scatter3d(;x=([path[1][1]]),y=([path[2][1]]),z=([path[3][1]]),
+                             mode="markers", showlegend=false,
+                             marker=attr(color=color, size=3, symbol="circle")))
+    push!(traces, scatter3d(;x=([path[1][end]]),y=([path[2][end]]),z=([path[3][end]]),
+                             mode="markers", showlegend=false,
+                             marker=attr(color=color, size=3, symbol="diamond")))
+  end
+
+  return traces, maximum(path[1:6])
+
+end
+
+"""
+A fundamental function to plot a single orbit path
+"""
+function plot(path::Vector{Vector{Float64}};
+              color::AbstractString = random_color(),
+              title::AbstractString = "Single Orbit Path",
+              markers=true)
+
+  # Generate the plot
+  p, limit = gen_plot(path, color=color, markers=markers)
+
+  # Now use the standard layout
+  layout = standard_layout(limit, title)
+
+  # Generate the plot
+  PlotlyJS.plot( PlotlyJS.Plot( p, layout ) )
+
+end
+
+"""
+This will generate the plot of a planetary body in it's correct colors and position in a
+sun-centered frame at a certain time
+
+If the body is the sun, no time is necessary
+
+An optional prop time is allowed, which will also provide an orbit
+"""
+function gen_plot(b::Body,
+                  time::Union{DateTime, Nothing}=nothing,
+                  prop_time::Union{Float64, Nothing}=nothing;
+                  scale::Int=100)
+
+  # Gotta be a little weird with Sun Scale
+  if b == Sun && scale > 25 scale = 25 end
+
+  # This is how many points to poll to generate the sphere
  N = 32
  θ = collect(range(0,length=N,stop=2π))
  ϕ = collect(range(0,length=N,stop=π))
-  x = cos.(θ) * sin.(ϕ)'
-  y = sin.(θ) * sin.(ϕ)'
-  z = repeat(cos.(ϕ)',outer=[N, 1])
-  ps = primary.r .* (x,y,z)
-  x_p,y_p,z_p = ps

-  t1 = []
-  for i = 1:length(paths)
-    path = paths[i]
-    label = labels != [] ? labels[i] : "orbit"
-    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))
-    push!(t1, scatter3d(;x=([path[1][1]]),y=([path[2][1]]),z=([path[3][1]]),
-                        mode="markers", showlegend=false,
-                        marker=attr(color=color, size=3, symbol="circle")))
-    push!(t1, scatter3d(;x=([path[1][end]]),y=([path[2][end]]),z=([path[3][end]]),
-                        mode="markers", showlegend=false,
-                        marker=attr(color=color, size=3, symbol="diamond")))
+  # Convert polar -> cartesian
+  unscaled_xs = cos.(θ) * sin.(ϕ)'
+  unscaled_ys = sin.(θ) * sin.(ϕ)'
+  unscaled_zs = repeat(cos.(ϕ)',outer=[N, 1])
+
+  # Scale by the size of the body (and a scaling factor for visualization)
+  xs,ys,zs = scale * b.r .* (unscaled_xs,unscaled_ys,unscaled_zs)
+
+  # Add an offset if we're not plotting the sun
+  body_state = zeros(6)
+  if b != Sun body_state = state(b, time) end
+
+  max_value = maximum(xs .+ body_state[1], ys .+ body_state[2], zs .+ body_state[3])
+
+  # Generate the surface trace
+  body_surface = surface(;x = xs .+ body_state[1],
+                          y = ys .+ body_state[2],
+                          z = zs .+ body_state[3],
+                          showscale = false,
+                          colorscale = b.color)
+  outputs = [ body_surface ]
+
+  if prop_time !== nothing
+    p, maybe_max = gen_plot(prop(body_state, prop_time), label=b.name, color=b.line_color, markers=false)
+    push!(outputs, p...)
+    if maybe_max > max_value max_value = maybe_max end
  end
-  limit = max(maximum(abs.(flatten(flatten(paths)))),
-              maximum(abs.(flatten(ps)))) * 1.1

-  t2 = surface(;x=(x_p),
-                y=(y_p),
-                z=(z_p),
-                showscale=false,
-                colorscale = primary.color)
-
-  layout = Layout(title=title,
-                  width=1000,
-                  height=600,
-                  paper_bgcolor="rgba(5,10,40,1.0)",
-                  plot_bgcolor="rgba(100,100,100,0.01)",
-                  scene = attr(xaxis = attr(autorange = false,range=[-limit,limit]),
-                               yaxis = attr(autorange = false,range=[-limit,limit]),
-                               zaxis = attr(autorange = false,range=[-limit,limit]),
-                  aspectratio=attr(x=1,y=1,z=1),
-                  aspectmode="manual"))
-
-  p = Plot([t1...,t2],layout)
-  plot(p)
+  return outputs, max_value
+
+end
+
+"""
+This will plot a planetary body in it's correct colors and position in a
+sun-centered frame at a certain time
+
+If the body is the sun, no time is necessary
+
+An optional prop time is allowed, which will also provide an orbit
+"""
+function plot(b::Body,
+              time::Union{DateTime, Nothing}=nothing,
+              prop_time::Union{Float64, Nothing}=nothing;
+              title::AbstractString="Single Planet Plot",
+              scale::Int=100)
+
+  # Generate the plot
+  p, limit = gen_plot(b,time,prop_time,scale=scale)
+
+  # Now use the standard layout
+  layout = standard_layout(limit, title)
+
+  # Generate the plot
+  PlotlyJS.plot( PlotlyJS.Plot( p, layout ) )
+
+end
+
+"""
+This will plot a vector of planetary bodies at a certain time
+
+The Sun is plotted automatically
+
+An optional orbits bool is allowed, which shows the orbits
+"""
+function plot(bodies::Vector{Body},
+              time::Union{DateTime, Nothing};
+              orbits::Bool=true,
+              title::AbstractString="Multiple Planet Plot",
+              scale::Int=100)
+
+  # Generate the plots
+  trace, limit = gen_plot(Sun, scale=scale)
+  traces = [ trace... ]
+  for b in bodies
+    if orbits
+      trace, new_limit = gen_plot(b,time,period(b),scale=scale)
+    else
+      trace, new_limit = gen_plot(b,time,scale=scale)
+    end
+    push!(traces, trace...)
+    limit = max(limit, new_limit)
+  end
+
+  # Now use the standard layout
+  layout = standard_layout(limit, title)
+
+  # Generate the plot
+  PlotlyJS.plot( PlotlyJS.Plot( traces, layout ) )
+
+end
+
+"""
+Generic plotting function for a bunch of paths. Useful for debugging.
+"""
+function plot(paths::Vector{Vector{Vector{Float64}}};
+              labels::Union{Vector{String},Nothing}=nothing,
+              colors::Union{Vector{String},Nothing}=nothing,
+              markers::Bool=true,
+              title::String="Spacecraft Position")
+
+  # Generate the path traces
+  color = colors === nothing ? random_color() : colors[1]
+  trace, limit = gen_plot(paths[1],label=labels[1],color=color,markers=markers)
+  traces = [ trace... ]
+  for i = 2:length(paths)
+    color = colors === nothing ? random_color() : colors[i]
+    trace, new_limit = gen_plot(paths[i],label=labels[i],color=color,markers=markers)
+    push!(traces, trace...)
+    limit = max(limit, new_limit)
+  end
+
+  # Generate the sun trace
+  push!(traces, gen_plot(Sun)[1]...)
+
+  # Determine layout details
+  layout = standard_layout(limit, title)
+
+  # Plot
+  PlotlyJS.plot( PlotlyJS.Plot( traces, layout ) )
+
+end
+
+function plot(m::Union{Mission, Mission_Guess}; title::String="Mision Plot")
+  # First plot the earth
+  # Then plot phase plus planet phase
+  time = m.launch_date
+  mass = m.start_mass
+  current_planet = Earth
+  earth_trace, limit = gen_plot(current_planet, time, year)
+  start = state(current_planet, time, m.launch_v∞, mass)
+  traces = [ earth_trace... ]
+
+  i = 1
+  for phase in m.phases
+    # First do the path
+    path, final = prop(phase.thrust_profile, start, m.sc, phase.tof, interpolate=true)
+    mass = final[7]
+    time += Dates.Second(floor(phase.tof))
+    current_planet = phase.planet
+    start = state(current_planet, time, phase.v∞_out, mass)
+    path_trace, new_limit = gen_plot(path, label="Phase "*string(i))
+    push!(traces, path_trace...)
+    limit = max(limit, new_limit)
+
+    # Then the planet
+    planet_trace, new_limit = gen_plot(current_planet, time, -phase.tof)
+    push!(traces, planet_trace...)
+    limit = max(limit, new_limit)
+
+    i += 1
+  end
+  
+  # Generate the sun trace
+  push!(traces, gen_plot(Sun)[1]...)
+
+  # Determine layout details
+  layout = standard_layout(limit, title)
+
+  # Plot
+  PlotlyJS.plot( PlotlyJS.Plot( traces, layout ) )

 end
--- a/julia/src/spacecraft.jl
+++ b/julia/src/spacecraft.jl
@@ -10,4 +10,4 @@ end

 const test_sc = Sc(8000., 0.00025/(2000*0.00981), 0.00025, 50, 0.9)
 const bepi = Sc(8000., 2*0.00025/(2000*0.00981), 0.00025, 2, 0.9)
-const no_thrust = Sc(8000., 0.01, 0., 0, 0.)
+const no_thrust = Sc(0., 0.01, 0., 0, 0.)
--- a/julia/test/inner_loop/nlp_solver.jl
+++ b/julia/test/inner_loop/nlp_solver.jl
@@ -1,36 +1,48 @@
@testset "Phase" begin

+  using PlotlyJS: savefig
+
  println("Testing NLP solver")

-  # We'll start by testing the mission_guess -> vector function
-  vec = Vector(test_mg)
-  @test typeof(vec) == Vector{Float64}
-
-  # Now we go in the other direction
-  flybys = [ p.planet for p in test_mg.phases ]
-  guess = Mission_Guess(vec, test_mg.sc, test_mg.start_mass, flybys)
-  @test typeof(guess) == Mission_Guess
-  @test guess.sc == test_mg.sc
-  @test guess.start_mass == test_mg.start_mass
-  @test guess.launch_date == test_mg.launch_date
-  @test guess.launch_v∞ == test_mg.launch_v∞
-  @test guess.converged == test_mg.converged
-  for i in 1:length(guess.phases)
-    @test guess.phases[i].planet == test_mg.phases[i].planet
-    @test guess.phases[i].v∞_in == test_mg.phases[i].v∞_in
-    @test guess.phases[i].v∞_out == test_mg.phases[i].v∞_out
-    @test guess.phases[i].tof == test_mg.phases[i].tof
-    @test guess.phases[i].thrust_profile == test_mg.phases[i].thrust_profile
-  end
-
-  # Now we test an example run of the basic "inner function"
+  # 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)
-  phase = Phase(Venus, v∞_in, v∞_in, tof, zeros(20,3))
-  guess = Mission_Guess(bepi, 12_000., leave, v∞_out, [phase])
-  g = solve_mission(guess)
-  println(g)
+  # 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, 12_000., test_leave, 0.9*v∞_out, [phase])
+  m = solve_mission(guess, launch_window, latest_arrival, verbose=true)
+  @test typeof(m) == Mission
+  @test m.converged == true
+
+  # Now we can plot the results to check visually
+  p = plot(m, title="NLP Test Solution")
+  savefig(p,"../plots/nlp_test_1_phase.html")
+
+  # Now we can look at a 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])
+  println(launch_v∞)
+  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)))
+  # This isn't quite right. V∞ discrepancy is too big
+  guess = Mission_Guess(bepi, 12_000., dates[1], launch_v∞, phases)
+  m = solve_mission(guess, launch_window, latest_arrival, verbose=true)
+  p = plot(m, title="NLP Test Solution (2 Phases)")
+  savefig(p,"../plots/nlp_test_2_phase.html")

 end
--- a/julia/test/plotting.jl
+++ b/julia/test/plotting.jl
@@ -2,34 +2,55 @@

  println("Testing plotting features")

-  using PlotlyJS
+  using PlotlyJS: savefig, SyncPlot

-  # Plot an earth plot
-  T = rand(2hour : 1 : 4hour)
-  revs = 30
-  n = revs*100
-
-  start = gen_orbit(T, 12_000., Earth)
-  thrust = spiral(0.9, n, start, test_sc, revs*T, Earth)
-  path = prop(thrust, start, test_sc, revs*T, Earth)[1]
-
-  p = plot_orbits([path], Earth)
-  savefig(p,"../plots/plot_test_earth.html")
-  @test typeof(p) == PlotlyJS.SyncPlot
-
-  # Now change a little bit and plot around the Sun
-  # This also checks that the spacecraft are configured right:
-  #   They really shouldn't run out of fuel in 4 years
+  # Generate a couple of useful paths
  T = rand(year : hour : 4year)
  tof = 4year
+  n = 4_000
  start = gen_orbit(T, 12_000.)
  thrust = spiral(0.9, n, start, bepi, tof)
-  sun_paths = Vector{Vector{Vector{Float64}}}()
-  push!(sun_paths, prop(zeros(100,3), start, bepi, tof)[1])
-  push!(sun_paths, prop(thrust, start, bepi, tof)[1])
-  p = plot_orbits(sun_paths)
-  savefig(p,"../plots/plot_test_sun.html")
-  @test typeof(p) == PlotlyJS.SyncPlot
+  spiral_path = prop(thrust, start, bepi, tof)[1]
+  no_thrust_path = prop(zeros(n,3), start, bepi, tof)[1]
+  paths = [ spiral_path, no_thrust_path ]

+  # First let's plot a single basic path
+  p = plot(spiral_path)
+  savefig(p,"../plots/plot_test_path.html")
+  @test typeof(p) == SyncPlot
+
+  # Next I just want to plot a planet both not moving...
+  p = plot(Jupiter, Dates.now(), title="Jupiter")
+  savefig(p,"../plots/plot_test_planet_nomove.html")
+  @test typeof(p) == SyncPlot
+
+  # ...and moving
+  p = plot(Venus, Dates.now(), 0.95period(Venus), title="Venus")
+  savefig(p,"../plots/plot_test_planet_move.html")
+  @test typeof(p) == SyncPlot
+
+  # Next I want to plot the entire solar system
+  p = plot([Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, Neptune, Pluto],
+           Dates.now(),
+           title="The Solar System")
+  savefig(p,"../plots/plot_test_solar_system.html")
+  @test typeof(p) == SyncPlot
+
+  # Next I want to plot a vector of states
+  p = plot(paths,
+           labels = ["No Thrust Orbit","Thrust Spiral"],
+           title = "Vector of Paths Plot Test")
+  savefig(p,"../plots/plot_test_paths.html")
+  @test typeof(p) == SyncPlot
+
+  # And finally, I'd like to plot a mission guess...
+  p = plot(test_mg)
+  savefig(p,"../plots/plot_test_mission_guess.html")
+  @test typeof(p) == SyncPlot
+
+  # ...and a mission
+  p = plot(Mission(test_mg))
+  savefig(p,"../plots/plot_test_mission.html")
+  @test typeof(p) == SyncPlot

 end
--- a/julia/test/runtests.jl
+++ b/julia/test/runtests.jl
@@ -14,9 +14,9 @@ catch
 end

@testset "All Tests" begin
-  # include("plotting.jl")
-  # include("inner_loop/laguerre-conway.jl")
-  # include("inner_loop/propagator.jl")
+  include("plotting.jl")
+  include("inner_loop/laguerre-conway.jl")
+  include("inner_loop/propagator.jl")
  include("inner_loop/nlp_solver.jl")
  # include("inner_loop/monotonic_basin_hopping.jl")
 end