using Dates

export inner_loop

"""
This is it. The outer function call for the inner loop. After this is done,
there's only the outer loop left to do. And that's pretty easy.
"""
function inner_loop(launch_date::DateTime,
                    craft::Sc,
                    phases::Vector{Phase};
                    min_flyby::Float64=1000.,
                    mbh_specs=nothing,
                    verbose=false)
   
  # First we need to do some quick checks that the mission is well formed
  for i in 1:length(phases)
    if i == 1
      @assert phases[i].from_planet == "Earth"
    else
      # Check that the planet is valid
      if phases[i].from_planet ∉ keys(μs)
        error("Planet is not valid: ", phases[i].from_planet)
      # Check that there is only one flyby planet
      elseif phases[i].from_planet != phases[i-1].to_planet
        fromP, toP = phases[i].from_planet, phases[i-1].to_planet
        error("Planets don't match up: (phase $(i)) $(fromP) / (phase $(i-1)) $(toP)")
      # Check that v∞_in == v∞_out
      elseif !isapprox(norm(phases[i].v∞_outgoing), norm(phases[i-1].v∞_incoming))
        norm_incoming = norm(phases[i-1].v∞_incoming)
        norm_outgoing = norm(phases[i].v∞_outgoing)
        error("""Norms of vectors aren't equal:
                 Phase $(i-1): $(phases[i-1].v∞_incoming) / norm: $(norm_incoming)
                 Phase $(i): $(phases[i].v∞_outgoing) / norm: $(norm_outgoing)""")
      end
      # Check that the approach is not too low
      v∞ = norm(phases[i].v∞_outgoing)
      δ = acos((phases[i].v∞_outgoing ⋅ phases[i-1].v∞_incoming)/v∞^2)
      flyby = μs[phases[i].from_planet]/v∞^2 * (1/sin(δ/2) - 1)
      true_min = rs[phases[i].from_planet] + min_flyby
      if flyby <= true_min
        error("Flyby was too low between phase $(i-1) and $(i): $(flyby) / $(true_min)")
      end
    end
  end

  time = utc2et(Dates.format(launch_date,"yyyy-mm-ddTHH:MM:SS"))
  thrust_profiles = []
  try
    for phase in phases
      planet1_state = spkssb(ids[phase.from_planet], time, "J2000")
      time += phase.time_of_flight
      planet2_state = spkssb(ids[phase.to_planet], time, "J2000")
      # TODO: Come up with improved method of calculating "n"
      start = planet1_state + [0., 0., 0., phase.v∞_outgoing...]
      final = planet2_state + [0., 0., 0., phase.v∞_incoming...]
      if mbh_specs === nothing
        best = mbh(start, final, craft, μs["Sun"], 0.0, phase.time_of_flight, 10,
                  verbose=verbose)[1]
      else
        num_iters, sil, dil = mbh_specs
        best = mbh(start, final, craft, μs["Sun"], 0.0, phase.time_of_flight, 10,
                  verbose=verbose, num_iters=num_iters, search_patience_lim=sil,
                  drill_patience_lim=dil)
      end
      push!(thrust_profiles, best)
      craft.mass = prop(tanh.(best.zero), planet1_state, sc, μs["Sun"], prop_time)[2][end]
    end
    return craft.mass, thrust_profiles
  catch
    return "One path did not converge"
  end
end