Might change the structure again...

2021-10-05 20:14:39 -06:00
parent 8e8a23c52a
commit cd1806058d
10 changed files with 280 additions and 176 deletions
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,4 @@
 pdf/
 **/*.pdf
 julia/**/*.html
+julia/test/ipopt.out
--- a/julia/src/Thesis.jl
+++ b/julia/src/Thesis.jl
@@ -16,6 +16,7 @@ module Thesis
  end

  include("./errors.jl")
+  include("./bodies.jl")
  include("./constants.jl")
  include("./spacecraft.jl")
  include("./mission.jl")
--- a/julia/src/bodies.jl
+++ b/julia/src/bodies.jl
@@ -0,0 +1,17 @@
+export state
+
+struct Body
+  μ::Float64
+  r::Float64 # radius
+  color::String
+  id::Int # SPICE id
+end
+
+function state(p::Body, t::DateTime, add_v∞::Vector{Float64}=[0., 0., 0.], add_mass::Float64=1e10)
+  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)
+  [ 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,13 +7,6 @@ export Jupiter, Saturn, Uranus, Neptune, Pluto
 export G, AU, init_STM, hour, day, year, second
 export Pathlist

-struct Body
-  μ::Float64
-  r::Float64 # radius
-  color::String
-  id::Int # SPICE id
-end
-
 const Sun = Body(1.32712440018e11, 696000., "Electric", 10)
 const Mercury = Body(2.2032e4, 2439., "heat", 1)
 const Venus = Body(3.257e5, 6052., "turbid", 2)
--- a/julia/src/errors.jl
+++ b/julia/src/errors.jl
@@ -2,12 +2,21 @@ struct LaGuerreConway_Error <: Exception end

 struct ΔVsize_Error <: Exception end

+struct HitPlanet_Error <: Exception end
+Base.showerror(io::IO, e::HitPlanet_Error) = print(io, "spacecraft hit the planet...")
+
 struct Convergence_Error <: Exception end
 Base.showerror(io::IO, e::Convergence_Error) = print(io, "NLP solver didn't converge...")

 struct GenOrbit_Error <: Exception end
 Base.showerror(io::IO, e::GenOrbit_Error) = print(io, "infinite Loop trying to generate the init orbit")

+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)
+
 struct PropOne_Error <: Exception
  ΔV_unit::AbstractVector
 end
--- a/julia/src/inner_loop/monotonic_basin_hopping.jl
+++ b/julia/src/inner_loop/monotonic_basin_hopping.jl
@@ -41,30 +41,35 @@ function gen_date(date_range::Tuple{DateTime, DateTime})
  l0 + Millisecond(floor(rand()*(lf-l0).value))
 end

-"""
-Returns a random amount of time in a range
-"""
-function gen_period(date_range::Tuple{DateTime, DateTime})
-  l0, lf = date_range
-  Millisecond(floor(rand()*(lf-l0).value))
-end
-
 """
 Perturbs a valid mission with pareto-distributed variables, generating a mission guess
 """
 function perturb(mission::Mission)
-  new_launch_date = mission.launch_date + Dates.Second(floor(7day * pareto_add()))
-  new_launch_v∞ = mission.launch_v∞ .* pareto(3)
-  new_phases = Vector{Phase}()
-  for phase in mission.phases
-    new_v∞_in = phase.v∞_in .* pareto(3)
-    new_δ = phase.δ * pareto()
-    new_tof = phase.tof * pareto()
-    new_thrust_profile = phase.thrust_profile .* pareto(size(phase.thrust_profile))
-    push!(new_phases, Phase(phase.planet, new_v∞_in, new_δ, new_tof, new_thrust_profile))
+  mission_guess = Bad_Mission("Starting point")
+  while typeof(mission_guess) == Bad_Mission
+    new_launch_date = mission.launch_date + Dates.Second(floor(7day * pareto_add()))
+    new_launch_v∞ = mission.launch_v∞ .* pareto(3)
+    new_phases = Vector{Phase}()
+    for phase in mission.phases
+      new_v∞_in = phase.v∞_in .* pareto(3)
+      temp_v∞_out = phase.v∞_out .* pareto(3)
+      new_v∞_out = norm(new_v∞_in) * temp_v∞_out/norm(temp_v∞_out)
+      new_tof = phase.tof * pareto()
+      new_thrust_profile = phase.thrust_profile .* pareto(size(phase.thrust_profile))
+      push!(new_phases, Phase(phase.planet, new_v∞_in, new_v∞_out, new_tof, new_thrust_profile))
+    end
+    # TODO: Mission_Guess.validate()
+    try
+      mission_guess = Mission_Guess(mission.sc, mission.start_mass, new_launch_date, new_launch_v∞, new_phases)
+    catch e
+      if isa(e, Mass_Error) ||isa(e,V∞_Error) || isa(e,HitPlanet_Error)
+        continue
+      else
+        rethrow()
+      end
+    end
  end
-  # TODO: Mission_Guess.validate()
-  Mission_Guess(mission.sc, mission.start_mass, new_launch_date, new_launch_v∞, new_phases)
+  return mission_guess
 end

 """
@@ -78,66 +83,63 @@ function mission_guess( flybys::Vector{Body},
                        max_v∞_in_mag::Float64,
                        latest_arrival::DateTime,
                        primary::Body=Sun  )
-  # TODO: Eventually I can calculate n more intelligently
-  n = 20
+  mission_guess = Bad_Mission("Keep trying to generate a guess")
+  while mission_guess == Bad_Mission("Keep trying to generate a guess")
+    # TODO: Eventually I can calculate n more intelligently
+    n = 20

-  # Determine the launch conditions
-  launch_date = gen_date(launch_window)
-  launch_v∞_normalized = rand(-1:0.0001:1, 3)
-  launch_v∞ = rand(0:0.0001:√max_C3_out) * launch_v∞_normalized/norm(launch_v∞_normalized)
+    # Determine the launch conditions
+    launch_date = gen_date(launch_window)
+    launch_v∞_normalized = rand(-1:0.0001:1, 3)
+    launch_v∞ = rand(0:0.0001:√max_C3_out) * launch_v∞_normalized/norm(launch_v∞_normalized)

-  # Determine the leg lengths
-  num_phases = length(flybys) - 1
-  total_tof = 100year
-  max_tof = (latest_arrival - launch_date).value / 1000
-  tofs = rand(30day : hour : 0.7max_tof, num_phases)
-  total_tof = sum(tofs)
-  while total_tof > max_tof
+    # Determine the leg lengths
+    num_phases = length(flybys) - 1
+    total_tof = 100year
+    max_tof = (latest_arrival - launch_date).value / 1000
    tofs = rand(30day : hour : 0.7max_tof, num_phases)
    total_tof = sum(tofs)
-  end
-
-  # For each phase, determine the v∞_in and δ
-  phases = Vector{Phase}()
-  for i in 1:num_phases
-    flyby = flybys[i+1]
-    v∞_normalized = rand(-1:0.0001:1, 3)
-    v∞ = rand(0:0.0001:10) * v∞_normalized/norm(v∞_normalized)
-    δ = rand(0:0.0001:2π)
-    periapsis = (flyby.μ/(v∞ ⋅ v∞)) * ( 1/sin(δ/2) - 1 )
-    while periapsis < flyby.r + 100.
-      δ = rand(0:0.0001:2π)
-      periapsis = (flyby.μ/(v∞ ⋅ v∞)) * ( 1/sin(δ/2) - 1 )
+    while total_tof > max_tof
+      tofs = rand(30day : hour : 0.7max_tof, num_phases)
+      total_tof = sum(tofs)
+    end
+
+    # For each phase, determine the v∞_in and δ
+    phases = Vector{Phase}()
+    for i in 1:num_phases
+      flyby = flybys[i+1]
+      v∞_in_normalized = rand(-1:0.0001:1, 3)
+      v∞_in = rand(0:0.0001:10) * v∞_in_normalized/norm(v∞_in_normalized)
+      v∞_out_normalized = rand(-1:0.0001:1, 3)
+      v∞_out = norm(v∞_in) * v∞_out_normalized/norm(v∞_out_normalized)
+      δ = acos( ( v∞_in ⋅ v∞_out ) / ( norm(v∞_in) * norm(v∞_out) ) )
+      periapsis = (flyby.μ/(v∞_in ⋅ v∞_in)) * ( 1/sin(δ/2) - 1 )
+      while periapsis < flyby.r + 100.
+        v∞_out_normalized = rand(-1:0.0001:1, 3)
+        v∞_out = norm(v∞_in) * v∞_out_normalized/norm(v∞_out_normalized)
+        δ = acos( ( v∞_in ⋅ v∞_out ) / ( norm(v∞_in) * norm(v∞_out) ) )
+        periapsis = (flyby.μ/(v∞_in ⋅ v∞_in)) * ( 1/sin(δ/2) - 1 )
+      end
+      thrusts = rand(-1:0.0001:1,(n,3))
+      push!(phases, Phase(flyby, v∞_in, v∞_out, tofs[i], thrusts))
+    end
+
+    # Finally, determine the arrival v∞
+    arrival_v∞_normalized = rand(-1:0.0001:1, 3)
+    arrival_v∞ = rand(0:0.0001:max_v∞_in_mag) * arrival_v∞_normalized/norm(arrival_v∞_normalized)
+
+    # And we can construct a mission guess object with these values
+    try
+      mission_guess = Mission_Guess( sc, start_mass, launch_date, launch_v∞, phases )
+    catch e
+      if isa(e, Mass_Error) ||isa(e,V∞_Error) || isa(e,HitPlanet_Error)
+        continue
+      else
+        rethrow()
+      end
    end
-    thrusts = rand(-1:0.0001:1,(n,3))
-    push!(phases, Phase(flyby, v∞, δ, tofs[i], thrusts))
  end
-
-  # Finally, determine the arrival v∞
-  arrival_v∞_normalized = rand(-1:0.0001:1, 3)
-  arrival_v∞ = rand(0:0.0001:max_v∞_in_mag) * arrival_v∞_normalized/norm(arrival_v∞_normalized)
-
-  # And we can construct a mission guess object with these values
-  Mission_Guess( sc, start_mass, launch_date, launch_v∞, phases )
-end
-
-"""
-A convenience function for calculating mass usage given a certain thrust profile
-"""
-function mass_consumption(sc::Sc, phase::Phase)
-  weighted_thrusting_time = 0.0
-  n = size(phase.thrust_profile)[1]
-  for i in 1:size(phase.thrust_profile,1)
-    weighted_thrusting_time += norm(phase.thrust_profile[i,:]) * phase.tof/n
-  end
-  return weighted_thrusting_time*sc.mass_flow_rate
-end
-
-"""
-This attempts to determine v∞_out from v∞_in and the turning angle, assuming we are staying in the
-plane of the three planets
-"""
-function calc_turn(p1::Body, p2::Body, p3::Body, v∞_in::Vector{Float64}, δ::Float64)
+  return mission_guess
 end

 """
@@ -145,25 +147,26 @@ Sequentially calls the NLP solver to attempt to solve based on the initial guess
 """
 function inner_loop_solve(guess::Mission_Guess)
  v∞_out = guess.launch_v∞
-  time = utc2et(Dates.format(guess.launch_date,"yyyy-mm-ddTHH:MM:SS"))
+  mass = guess.start_mass
  current_planet = Earth
-  start = [spkssb(Earth.id, time, "ECLIPJ2000"); 0.0] + [ zeros(3); guess.launch_v∞; guess.start_mass ]
+  time = utc2et(Dates.format(guess.launch_date,"yyyy-mm-ddTHH:MM:SS"))
+  start = state(current_planet, time, v∞_out, mass)

  corrected_phases = Vector{Phase}()
  for i in 1:length(guess.phases)
    phase = guess.phases[i]
+    current_planet = phase.planet
    time += phase.tof
-    goal = spkssb(phase.planet.id, time, "ECLIPJ2000") + [zeros(3); phase.v∞_in] 
+    goal = state(current_planet, time, phase.v∞_in)
    result = solve_phase( start, goal, guess.sc, phase.tof, phase.thrust_profile)
    result.converged || return Bad_Mission(result.info) # Drop if it's not working
-    corrected_phase = Phase(phase.planet, phase.v∞_in, phase.δ, phase.tof, result.sol)
+    corrected_phase = Phase(phase.planet, phase.v∞_in, phase.v∞_out, phase.tof, result.sol)
    push!(corrected_phases, corrected_phase)
    mass_used = mass_consumption(guess.sc, corrected_phase)
+    mass -= mass_used
    if i != length(guess.phases)
-      v∞_out = calc_turn(current_planet, phase.planet, guess.phases[i+1].planet, phase.v∞_in, phase.δ)
-      current_planet = phase.planet
-      planet_state = [spkssb(current_planet.id, time, "ECLIPJ2000"); 0.0]
-      start = planet_state + [ zeros(3); v∞_out; start_mass - mass_used ]
+      v∞_out = phase.v∞_out
+      start = state(current_planet, time, v∞_out, mass)
    end
  end

@@ -183,8 +186,13 @@ function cost(mission::Mission, max_C3::Float64, max_v∞::Float64)
  return 3mass_percent + C3_percent + v∞_percent
 end

+cost(_::Bad_Mission) = Inf

-
+"""
+This is the main monotonic basin hopping function. There's a lot going on here, but the general idea
+is that hopefully you can provide mission parameters and a list of flybys and get the optimal
+mission that uses those flybys.
+"""
 function mbh( flybys::Vector{Body},
              sc::Sc,
              start_mass::Float64,
@@ -193,49 +201,64 @@ function mbh( flybys::Vector{Body},
              max_v∞::Float64,
              latest_arrival::DateTime,
              primary::Body=Sun;
-              search_patience::Int=1_000,
-              drill_patience::Int=50)
-  # Let's pseudo-code this bitch
-  #
-  # First, we need a function (mission_guess below) that randomly generates the:
-  #   - Launch date
-  #   - Launch v∞_out
-  #   - for each phase:
-  #       - tof
-  #       - v∞_in
-  #       - turning angle
-  #
-  # Also need a function (inner_loop_solve below) that takes the generated decision vector guess and
-  # attempts to correct it with the NLP solver. It will either converge or not.
-  #
-  # Also need a costing function (may be provided potentially)
-  #
-  # Also need a perturb function
-  #
-  guess = mission_guess(flybys, sc, start_mass, launch_window, max_C3, max_v∞, latest_arrival, primary)
-  inner_loop_solve(guess)
-  # search_count = 0
-  # x_current = "Bad"
-  # while search_count < search_patience
-    # search_count += 1
-    # drill_count = 0
-    # x_star = inner_loop_solve(mission_guess())
-    # if x_star.converged
-      # if cost(x_star) < cost(x_current)
-        # x_current = x_star
-      # while drill_count < drill_patience
-        # x_star = inner_loop_solve(perturb(x_current))
-        # if x_star.converged and cost(x_star) < cost(x_current)
-          # x_current = x_star
-          # drill_count = 0
-        # else
-          # drill_count += 1
-        # end
-      # end
-      # push!(archive, x_current)
-    # end
-  # end
-  #  
-  #
+              search_patience::Int=10_000,
+              drill_patience::Int=50,
+              verbose::Bool=false)
+
+  # A convenience function
+  random_guess() = mission_guess(flybys,sc,start_mass,launch_window,max_C3,max_v∞,latest_arrival,primary)
+  cost(m::Mission) = cost(m, max_C3, max_v∞)
+  status(s,d,l) = print("\r\t", "search: ", s, " drill: ", d, " archive length: ", l, "     ")
+
+  # Initialize stuff
+  search_count = 0
+  x_current = Bad_Mission("Starting point")
+  archive = Vector{Mission}()
+
+  # The main loop
+  if verbose println("Starting Monotonic Basin Hopper") end
+  while search_count < search_patience
+
+    # Intialize an x_star, if it doesn't converge, hop on to the next basin
+    search_count += 1
+    drill_count = 0
+    if verbose status(search_count, drill_count, length(archive)) end
+    x_star = inner_loop_solve(random_guess())
+    x_star.converged || continue
+    x_basin = x_star
+
+
+    # If it does, though, we check to see if it's better than the current
+    if cost(x_star) < cost(x_current) x_current = x_star end
+
+    # Either way, we need to drill into this particular basin, since it's valid
+    while drill_count < drill_patience
+
+      if verbose status(search_count, drill_count, length(archive)) end
+
+      #Perturb to generate a slightly different x_star
+      x_star = inner_loop_solve(perturb(x_basin))
+
+      # If better than the best, then keep it as current
+      if x_star.converged && cost(x_star) < cost(x_current)
+        x_current = x_star
+        x_basin = x_star
+        drill_count = 0
+      # If better than the best in this particular basin, keep it as basin
+      elseif x_star.converged && cost(x_star) < cost(x_basin)
+        x_basin = x_star
+        drill_count = 0
+      # Otherwise, keep drilling
+      else
+        drill_count += 1
+      end
+    end
+
+    x_current in archive || push!(archive, x_current)
+
+  end
+
+  return x_current, archive
+
 end

--- a/julia/src/inner_loop/phase.jl
+++ b/julia/src/inner_loop/phase.jl
@@ -39,15 +39,17 @@ function solve_phase( start::Vector{Float64},
  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)
-  # if verbose ipopt_options["print_level"] = 5 end
+                       "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
--- a/julia/src/mission.jl
+++ b/julia/src/mission.jl
@@ -1,17 +1,30 @@
 export Phase, Mission_Guess, Mission, Bad_Mission
 export test_phase1, test_phase2
-export test_mission_guess, test_mission_guess_simple
-export test_mission, test_mission_simple
+export test_mg, test_mission

 struct Phase
  planet::Body
  v∞_in::Vector{Float64}
-  δ::Float64
+  v∞_out::Vector{Float64}
  tof::Float64
  thrust_profile::Matrix{Float64}
 end
-const test_phase1 = Phase(Venus, [10.4321, -6.3015, -0.01978], 0.2, 1.30464e7, zeros(20,3))
-const test_phase2 = Phase(Jupiter, [0.3, 7.1, 0.2], 2π, 3.9year, zeros(20,3))
+
+const test_phase1 = Phase(Venus, [9., 10., 0.], [10., 9, 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]
+
+"""
+A convenience function for calculating mass usage given a certain thrust profile
+"""
+function mass_consumption(sc::Sc, phase::Phase)
+  weighted_thrusting_time = 0.0
+  n = size(phase.thrust_profile)[1]
+  for i in 1:size(phase.thrust_profile,1)
+    weighted_thrusting_time += norm(phase.thrust_profile[i,:]) * phase.tof/n
+  end
+  return weighted_thrusting_time*sc.mass_flow_rate
+end

 struct Mission_Guess
  sc::Sc
@@ -21,17 +34,28 @@ struct Mission_Guess
  phases::Vector{Phase}
  converged::Bool
 end
-Mission_Guess(args...) = Mission_Guess(args..., false)
-const test_mission_guess = Mission_Guess( bepi,
-                                          12_000.,
-                                          DateTime(1992, 11, 19),
-                                          [-3.4, 1.2, 0.1],
-                                          [test_phase1, test_phase2] )
-const test_mission_guess_simple = Mission_Guess(bepi,
-                                                12_000.,
-                                                DateTime(1992, 11, 19),
-                                                [-3.4, 1.2, 0.1],
-                                                [test_phase1])
+
+"""
+Constructor for a mission guess. Generally mission guesses are not converged
+"""
+function Mission_Guess(sc::Sc, mass::Float64, date::DateTime, v∞::Vector{Float64}, phases::Vector{Phase})
+  # 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())
+    end
+  end
+  Mission_Guess(sc, mass, date, v∞, phases, false)
+end
+
+const test_mg = Mission_Guess(bepi, 12_000., DateTime(1992,11,19), [-3.4,1.2,0.1], test_phases)
            
 struct Mission
  sc::Sc
@@ -41,20 +65,33 @@ struct Mission
  phases::Vector{Phase}
  converged::Bool
 end
-Mission(args...) = Mission(args..., true)
-const test_mission = Mission(bepi,
-                             12_000.,
-                             DateTime(1992, 11, 19),
-                             [-3.4, 1.2, 0.1],
-                             [test_phase1, test_phase2])
-const test_mission_simple = Mission(bepi,
-                                    12_000.,
-                                    DateTime(1992, 11, 19),
-                                    [4.2984, -4.3272668, 1.43752],
-                                    [test_phase1])
+
+"""
+Constructor for a mission. Generally mission guesses are converged
+"""
+function Mission(sc::Sc, mass::Float64, date::DateTime, v∞::Vector{Float64}, phases::Vector{Phase})
+  # 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())
+    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)
 
 struct Bad_Mission
-  message::Symbol
+  message::String
  converged::Bool
 end
-Bad_Mission(s) = Bad_Mission(s,false)
+Bad_Mission(s::String) = Bad_Mission(s,false)
+Bad_Mission(s::Symbol) = Bad_Mission(String(s),false)
+
--- a/julia/test/inner_loop/monotonic_basin_hopping.jl
+++ b/julia/test/inner_loop/monotonic_basin_hopping.jl
@@ -5,42 +5,54 @@
  println("Testing Monotonic Basin Hopper")

  # First we test the random mission guess generator
+  println("Testing guess generator function")
  flybys = [Earth, Venus, Jupiter]
  launch_window = ( DateTime(2021,12,25), DateTime(2025,12,25) )
  max_C3 = 10.
  max_v∞ = 8.
-  latest_arrival = DateTime(2035,12,25)
+  latest_arrival = DateTime(2030,12,25)
  random_guess = Thesis.mission_guess(flybys, bepi, 12_000., launch_window, max_C3, max_v∞, latest_arrival)
  @test typeof(random_guess) == Mission_Guess
  
  # Then the perturb function
+  println("Testing perturb function")
  mission_guess = Thesis.perturb(test_mission)
  @test mission_guess.launch_date != test_mission.launch_date
  @test mission_guess.launch_v∞ != test_mission.launch_v∞
  for i in 1:2
    @test mission_guess.phases[i].v∞_in != test_mission.phases[i].v∞_in
-    @test mission_guess.phases[i].δ != test_mission.phases[i].δ
+    @test mission_guess.phases[i].v∞_out != test_mission.phases[i].v∞_out
    @test mission_guess.phases[i].tof != test_mission.phases[i].tof
  end
  @test !mission_guess.converged
  
  # # Then the inner loop builder function
-  mission = Thesis.inner_loop_solve(test_mission_guess)
+  println("Testing inner loop solver function")
+  mission = Thesis.inner_loop_solve(test_mg)
  @test !mission.converged

  # For the valid case we need to use a lambert's solver
-  # TODO: This is probably not acceptable for how close I have to be
-  leave = DateTime(1992,11,19)
-  arrive = DateTime(1993,4,1)
-  time_leave = utc2et(Dates.format(leave,"yyyy-mm-ddTHH:MM:SS"))
-  time_arrive = utc2et(Dates.format(arrive,"yyyy-mm-ddTHH:MM:SS"))
-  earth_state = [spkssb(Earth.id, time_leave, "ECLIPJ2000"); 0.0]
-  venus_state = [spkssb(Venus.id, time_arrive, "ECLIPJ2000"); 0.0]
+  leave, arrive = DateTime(1992,11,19), DateTime(1993,4,1)
+  earth_state = state(Earth, leave)
+  venus_state = state(Venus, arrive)
  v∞_out, v∞_in, tof = Thesis.lamberts(Earth, Venus, leave, arrive)
-  phase = Phase(Venus, 1.01v∞_in, 0.2, tof, 0.01*ones(20,3))
+  phase = Phase(Venus, 1.01v∞_in, -1.01v∞_in, tof, 0.01*ones(20,3))
  mission_guess = Mission_Guess(bepi, 12_000., leave, v∞_out, [phase])
  mission = Thesis.inner_loop_solve(mission_guess)
  @test mission.converged
  if !mission.converged println(mission.message) end

+  # Now we're ready to test the whole thing!
+  println("Testing whole thing")
+  flybys = [ Earth, Venus ]
+  start_mass = 12_000.
+  launch_window = DateTime(1997, 10, 1), DateTime(1997, 10, 31)
+  max_C3 = 15.
+  max_v∞ = 20.
+  latest_arrival = DateTime(2001, 1, 1)
+  best, archive = mbh(flybys, bepi, start_mass, launch_window, max_C3, max_v∞, latest_arrival,
+                      verbose=true)
+  @test typeof(best) == Mission
+
+
 end
--- a/notes.md
+++ b/notes.md
@@ -33,4 +33,13 @@ Todo:

 1. ~~Add interpolate function~~
 2. ~~Add NLP solver sun test~~
-3. Start on MBH
+3. ~~Start on MBH~~
+
+## Refactor Notes (2021-09-26)
+
+Todo:
+
+1. Add mission validate function
+2. Add turning angle function
+3. Finish mbh
+4. Test it