Ok, now open loop is working, sc mass changed to state, and other updates

julia/src/inner_loop/find_closest.jl

@@ -5,16 +5,13 @@ export nlp_solve, mass_est
 function mass_est(T)
   ans = 0
   n = Int(length(T)/3)
-  for i in 1:n
-    ans += norm(T[i,:])
-  end
+  for i in 1:n ans += norm(T[i,:]) end
   return ans/n
 end
 
-converged(x) = NLsolve.converged(x)
-
-function converged(_::String)
-  return false
+struct Result
+  converged::Bool
+  zero::Matrix{Float64}
 end
 
 function nlp_solve(start::Vector{Float64},
@@ -28,10 +25,21 @@ function nlp_solve(start::Vector{Float64},
                    num_iters=1_000)
 
   function f!(F,x)
-    F .= 0.0
-    F[1:6, 1] .= prop_nlsolve(tanh.(x), start, craft, μ, tf-t0) .- final
+    try
+      F .= 0.0
+      F[1:6, 1] .= prop(tanh.(x), start, copy(craft), μ, tf-t0)[2][1:6] .- final[1:6]
+    catch e
+      F .= 10000000.0
+    end
   end
 
-  return nlsolve(f!, atanh.(x0), ftol=tol, autodiff=:forward, iterations=num_iters)
+  result = Result(false, zeros(size(x0)))
+  try
+    nl_results = nlsolve(f!, atanh.(x0), ftol=tol, autodiff=:forward, iterations=num_iters)
+    result = Result(converged(nl_results), tanh.(nl_results.zero))
+  catch e
+  end
 
-end
+  return result
+
+end

julia/src/inner_loop/inner_loop.jl

@@ -8,6 +8,7 @@ there's only the outer loop left to do. And that's pretty easy.
 """
 function inner_loop(launch_date::DateTime,
                     craft::Sc,
+                    start_mass::Float64,
                     phases::Vector{Phase};
                     min_flyby::Float64=1000.,
                     mbh_specs=nothing,
@@ -38,36 +39,32 @@ function inner_loop(launch_date::DateTime,
       δ = 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
+      flyby <= true_min || error("Flyby too low from phase $(i-1) to $(i): $(flyby) / $(true_min)")
     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]
+
+  for phase in phases
+    planet1_state = [spkssb(ids[phase.from_planet], time, "J2000"); 0.0]
+    time += phase.time_of_flight
+    planet2_state = [spkssb(ids[phase.to_planet], time, "J2000"); 0.0]
+    start = planet1_state + [0., 0., 0., phase.v∞_outgoing..., start_mass]
+    final = planet2_state + [0., 0., 0., phase.v∞_incoming..., start_mass]
+    println(start)
+    println(final)
+    # TODO: Come up with improved method of calculating "n"
+    if mbh_specs === nothing
+      best = mbh(start, final, craft, μs["Sun"], 0.0, phase.time_of_flight, 20,
+                 verbose=verbose)[1]
+    else
+      sil, dil = mbh_specs
+      best = mbh(start, final, craft, μs["Sun"], 0.0, phase.time_of_flight, 20,
+                 verbose=verbose, search_patience_lim=sil, drill_patience_lim=dil)[1]
     end
-    return craft.mass, thrust_profiles
-  catch
-    return "One path did not converge"
+    push!(thrust_profiles, best.zero)
   end
+  return thrust_profiles
+
 end

julia/src/inner_loop/laguerre-conway.jl

@@ -1,4 +1,4 @@
-function laguerre_conway(state::Vector{T}, μ::Float64, time::Float64) where T
+function laguerre_conway(state::Vector{<:Real}, μ::Float64, time::Float64)
 
   n = 5                               # Choose LaGuerre-Conway "n"
   i = 0

julia/src/inner_loop/monotonic_basin_hopping.jl

@@ -26,42 +26,38 @@ function mbh(start::AbstractVector,
              t0::AbstractFloat,
              tf::AbstractFloat,
              n::Int;
-             num_iters=25,
-             search_patience_lim::Int=200,
-             drill_patience_lim::Int=200,
+             search_patience_lim::Int=2000,
+             drill_patience_lim::Int=40,
              tol=1e-6,
              verbose=false)
 
   archive = []
 
   i = 0
-  if verbose println("Current Iteration") end
-  while true
+  x_current = Result(false, 1e8*ones(n,3))
+  while i < search_patience_lim
     i += 1
-    if verbose print("\r",i) end
-    search_impatience = 0
     drill_impatience = 0
-    x_star = nlp_solve(start, final, craft, μ, t0, tf, new_x(n), tol=tol, num_iters=100)
-    while converged(x_star) == false && search_impatience < search_patience_lim
-      search_impatience += 1
-      x_star = nlp_solve(start, final, craft, μ, t0, tf, new_x(n), tol=tol, num_iters=100)
-    end
-    if drill_impatience > drill_patience_lim break end
-    drill_impatience = 0
-    if converged(x_star)
+    if verbose print("\r\t", "search: ", i, " drill: ", drill_impatience, "    ") end
+    x_star = nlp_solve(start, final, craft, μ, t0, tf, new_x(n), tol=tol)
+    # If x_star is converged and better, set new x_current
+    if x_star.converged && mass_est(x_star.zero) < mass_est(x_current.zero)
       x_current = x_star
+    end
+    # If x_star is converged, drill down, otherwise, start over
+    if x_star.converged
       while drill_impatience < drill_patience_lim
-        x_star = nlp_solve(start, final, craft, μ, t0, tf, perturb(tanh.(x_current.zero),n), tol=tol)
-        if converged(x_star) && mass_est(tanh.(x_star.zero)) < mass_est(tanh.(x_current.zero))
+        x_star = nlp_solve(start, final, craft, μ, t0, tf, perturb(x_current.zero,n), tol=tol)
+        if x_star.converged && mass_est(x_star.zero) < mass_est(x_current.zero)
           x_current = x_star
           drill_impatience = 0
         else
+          if verbose print("\r\t", "search: ", i, " drill: ", drill_impatience, "    ") end
           drill_impatience += 1
         end
       end
       push!(archive, x_current)
     end
-    if i >= num_iters break end
   end
   if verbose println() end
 
@@ -70,8 +66,8 @@ function mbh(start::AbstractVector,
   current_best_mass = 1e8
   best = archive[1]
   for candidate in archive
-    if mass_est(tanh.(candidate.zero)) < current_best_mass
-      current_best_mass = mass_est(tanh.(candidate.zero))
+    if mass_est(candidate.zero) < current_best_mass
+      current_best_mass = mass_est(candidate.zero)
       best = candidate
     end
   end

julia/src/inner_loop/propagator.jl

@@ -3,234 +3,80 @@ export prop
 """
 Maximum ΔV that a spacecraft can impulse for a given single time step
 """
-function max_ΔV(duty_cycle::T,
+function max_ΔV(duty_cycle::Float64,
                 num_thrusters::Int,
-                max_thrust::T,
-                tf::T,
-                t0::T,
-                mass::S) where {T <: Real, S <: Real}
+                max_thrust::Float64,
+                tf::Float64,
+                t0::Float64,
+                mass::T) where T <: Real
   return duty_cycle*num_thrusters*max_thrust*(tf-t0)/mass
 end
 
-"""
-This function propagates the spacecraft forward in time 1 Sim-Flanagan step (of variable length of time),
-applying a thrust in the center.
-"""
-function prop_one(thrust_unit::Vector{<:Real},
-                  state::Vector{<:Real},
-                  duty_cycle::Float64,
-                  num_thrusters::Int,
-                  max_thrust::Float64,
-                  mass::T,
-                  mass_flow_rate::Float64,
-                  μ::Float64,
-                  time::Float64) where T<:Real
-
-  ΔV = max_ΔV(duty_cycle, num_thrusters, max_thrust, time, 0., mass) * thrust_unit
-  halfway = laguerre_conway(state, μ, time/2) + [0., 0., 0., ΔV...]
-  return laguerre_conway(halfway, μ, time/2), mass - mass_flow_rate*norm(thrust_unit)*time
-
-end
-
 """
 A convenience function for using spacecraft. Note that this function outputs a sc instead of a mass
 """
-function prop_one(ΔV_unit::Vector{T},
-                  state::Vector{S},
+function prop_one(ΔV_unit::Vector{<:Real},
+                  state::Vector{<:Real},
                   craft::Sc,
                   μ::Float64,
-                  time::Float64) where {T <: Real,S <: Real}
-  state, mass =  prop_one(ΔV_unit, state, craft.duty_cycle, craft.num_thrusters, craft.max_thrust,
-                          craft.mass, craft.mass_flow_rate, μ, time)
-  return state, Sc(mass, craft.mass_flow_rate, craft.max_thrust, craft.num_thrusters, craft.duty_cycle)
+                  time::Float64)
+
+  for direction in ΔV_unit
+    if abs(direction) > 1.0
+      println(direction)
+      error("ΔV is impossibly high")
+    end
+  end
+  ΔV = max_ΔV(craft.duty_cycle, craft.num_thrusters, craft.max_thrust, time, 0., state[7]) * ΔV_unit
+  halfway = laguerre_conway(state, μ, time/2) + [zeros(3); ΔV]
+  final = laguerre_conway(halfway, μ, time/2)
+  return [final; state[7] - craft.mass_flow_rate*norm(ΔV_unit)*time]
+
 end
 
-
 """
-This propagates over a given time period, with a certain number of intermediate steps
+The propagator function
 """
 function prop(ΔVs::Matrix{T},
               state::Vector{Float64},
-              duty_cycle::Float64,
-              num_thrusters::Int,
-              max_thrust::Float64,
-              mass::Float64,
-              mass_flow_rate::Float64,
+              craft::Sc,
               μ::Float64,
               time::Float64) where T <: Real
 
-  if size(ΔVs)[2] != 3 throw(ErrorException("ΔV input is wrong size")) end
-  n = size(ΔVs)[i]
-
-  for i in 1:n
-    state, mass = prop_one(ΔVs[i,:], state, duty_cycle, num_thrusters, max_thrust, mass,
-                           mass_flow_rate, μ, time/n)
-  end
-
-  return state, mass
-
-end
-
-"""
-The same function, using Scs
-"""
-function prop(ΔVs::Matrix{T},
-              state::Vector{S},
-              craft::Sc,
-              μ::Float64,
-              time::Float64) where {T <: Real, S <: Real}
-
   if size(ΔVs)[2] != 3 throw(ErrorException("ΔV input is wrong size")) end
   n = size(ΔVs)[1]
 
-  x_states = [state[1]]
-  y_states = [state[2]]
-  z_states = [state[3]]
-  dx_states = [state[4]]
-  dy_states = [state[5]]
-  dz_states = [state[6]]
-  masses = [craft.mass]
+  x_states = Vector{T}()
+  y_states = Vector{T}()
+  z_states = Vector{T}()
+  dx_states = Vector{T}()
+  dy_states = Vector{T}()
+  dz_states = Vector{T}()
+  masses = Vector{T}()
+
+  push!(x_states, state[1])
+  push!(y_states, state[2])
+  push!(z_states, state[3])
+  push!(dx_states, state[4])
+  push!(dy_states, state[5])
+  push!(dz_states, state[6])
+  push!(masses, state[7])
 
   for i in 1:n
-    state, craft = prop_one(ΔVs[i,:], state, craft, μ, time/n)
+    state = prop_one(ΔVs[i,:], state, craft, μ, time/n)
     push!(x_states, state[1])
     push!(y_states, state[2])
     push!(z_states, state[3])
     push!(dx_states, state[4])
     push!(dy_states, state[5])
     push!(dz_states, state[6])
-    push!(masses, craft.mass)
+    push!(masses, state[7])
+    if state[7] < craft.dry_mass
+      println(state[7])
+      error("Mass is too low")
+    end
   end
 
-  return [x_states, y_states, z_states, dx_states, dy_states, dz_states], masses, state
+  return [x_states, y_states, z_states, dx_states, dy_states, dz_states, masses], state
 
 end
-
-function prop_nlsolve(ΔVs::Matrix{T},
-                      state::Vector{S},
-                      craft::Sc,
-                      μ::Float64,
-                      time::Float64) where {T <: Real, S <: Real}
-
-  n = size(ΔVs)[1]
-  try
-    for i in 1:n
-      state, craft = prop_one(ΔVs[i,:], state, craft, μ, time/n)
-    end
-    return state
-  catch
-    return [0., 0., 0., 0., 0., 0.]
-  end
-
-
-end
-
-function prop_simple(ΔVs::AbstractMatrix,
-                     state::AbstractVector,
-                     craft::Sc,
-                     μ::Float64,
-                     time::Float64)
-
-  if size(ΔVs)[2] != 3 throw(ErrorException("ΔV input is wrong size")) end
-  n = size(ΔVs)[1]
-
-  for i in 1:n
-    state, craft = prop_one(ΔVs[i,:], state, craft, μ, time/n)
-  end
-
-  return state
-
-end
-
-function prop_one_simple(Tx, Ty, Tz, x, y, z, dx, dy, dz, t, μ)
-
-  # perform laguerre_conway once
-  r0_mag = √(x^2 + y^2 + z^2)
-  v0_mag = √(dx^2 + dy^2 + dz^2)
-  σ0 = ([x, y, z] ⋅ [dx, dy, dz])/√(μ)
-  a = 1 / ( 2/r0_mag - v0_mag^2/μ )
-  coeff = 1 - r0_mag/a
-
-  if a > 0    # Elliptical
-    ΔM = ΔE_new = √(μ) / sqrt(a^3) * t/2
-    ΔE = 1000
-    while abs(ΔE - ΔE_new) > 1e-10
-      ΔE = ΔE_new
-      F = ΔE - ΔM + σ0 / √(a) * (1-cos(ΔE)) - coeff * sin(ΔE)
-      dF = 1 + σ0 / √(a) * sin(ΔE) - coeff * cos(ΔE)
-      d2F = σ0 / √(a) * cos(ΔE) + coeff * sin(ΔE)
-      sign = dF >= 0 ? 1 : -1
-      ΔE_new = ΔE - 5*F / ( dF + sign * √(abs((5-1)^2*dF^2 - 5*(5-1)*F*d2F )))
-    end
-    F = 1 - a/r0_mag * (1-cos(ΔE))
-    G = a * σ0/ √(μ) * (1-cos(ΔE)) + r0_mag * √(a) / √(μ) * sin(ΔE)
-    r = a + (r0_mag - a) * cos(ΔE) + σ0 * √(a) * sin(ΔE)
-    Ft = -√(a)*√(μ) / (r*r0_mag) * sin(ΔE)
-    Gt = 1 - a/r * (1-cos(ΔE))
-  else      # Hyperbolic or Parabolic
-    ΔN = √(μ) / sqrt(-a^3) * t/2
-    ΔH = 0
-    ΔH_new = t/2 < 0 ? -1 : 1
-    while abs(ΔH - ΔH_new) > 1e-10
-      ΔH = ΔH_new
-      F = -ΔN - ΔH + σ0 / √(-a) * (cos(ΔH)-1) + coeff * sin(ΔH)
-      dF = -1 + σ0 / √(-a) * sin(ΔH) + coeff * cos(ΔH)
-      d2F = σ0 / √(-a) * cos(ΔH) + coeff * sin(ΔH)
-      sign = dF >= 0 ? 1 : -1
-      ΔH_new = ΔH - 5*F / ( dF + sign * √(abs((5-1)^2*dF^2 - 5*(5-1)*F*d2F )))
-    end
-    F = 1 - a/r0_mag * (1-cos(ΔH))
-    G = a * σ0/ √(μ) * (1-cos(ΔH)) + r0_mag * √(-a) / √(μ) * sin(ΔH)
-    r = a + (r0_mag - a) * cos(ΔH) + σ0 * √(-a) * sin(ΔH)
-    Ft = -√(-a)*√(μ) / (r*r0_mag) * sin(ΔH)
-    Gt = 1 - a/r * (1-cos(ΔH))
-  end
-
-  # add the thrust vector
-  x,y,z,dx,dy,dz = [F*[x,y,z] + G*[dx,dy,dz]; Ft*[x,y,z] + Gt*[dx,dy,dz] + [Tx, Ty, Tz]]
-
-  #perform again
-  r0_mag = √(x^2 + y^2 + z^2)
-  v0_mag = √(dx^2 + dy^2 + dz^2)
-  σ0 = ([x, y, z] ⋅ [dx, dy, dz])/√(μ)
-  a = 1 / ( 2/r0_mag - v0_mag^2/μ )
-  coeff = 1 - r0_mag/a
-
-  if a > 0    # Elliptical
-    ΔM = ΔE_new = √(μ) / sqrt(a^3) * t/2
-    ΔE = 1000
-    while abs(ΔE - ΔE_new) > 1e-10
-      ΔE = ΔE_new
-      F = ΔE - ΔM + σ0 / √(a) * (1-cos(ΔE)) - coeff * sin(ΔE)
-      dF = 1 + σ0 / √(a) * sin(ΔE) - coeff * cos(ΔE)
-      d2F = σ0 / √(a) * cos(ΔE) + coeff * sin(ΔE)
-      sign = dF >= 0 ? 1 : -1
-      ΔE_new = ΔE - 5*F / ( dF + sign * √(abs((5-1)^2*dF^2 - 5*(5-1)*F*d2F )))
-    end
-    F = 1 - a/r0_mag * (1-cos(ΔE))
-    G = a * σ0/ √(μ) * (1-cos(ΔE)) + r0_mag * √(a) / √(μ) * sin(ΔE)
-    r = a + (r0_mag - a) * cos(ΔE) + σ0 * √(a) * sin(ΔE)
-    Ft = -√(a)*√(μ) / (r*r0_mag) * sin(ΔE)
-    Gt = 1 - a/r * (1-cos(ΔE))
-  else      # Hyperbolic or Parabolic
-    ΔN = √(μ) / sqrt(-a^3) * t/2
-    ΔH = 0
-    ΔH_new = t/2 < 0 ? -1 : 1
-    while abs(ΔH - ΔH_new) > 1e-10
-      ΔH = ΔH_new
-      F = -ΔN - ΔH + σ0 / √(-a) * (cos(ΔH)-1) + coeff * sin(ΔH)
-      dF = -1 + σ0 / √(-a) * sin(ΔH) + coeff * cos(ΔH)
-      d2F = σ0 / √(-a) * cos(ΔH) + coeff * sin(ΔH)
-      sign = dF >= 0 ? 1 : -1
-      ΔH_new = ΔH - 5*F / ( dF + sign * √(abs((5-1)^2*dF^2 - 5*(5-1)*F*d2F )))
-    end
-    F = 1 - a/r0_mag * (1-cos(ΔH))
-    G = a * σ0/ √(μ) * (1-cos(ΔH)) + r0_mag * √(-a) / √(μ) * sin(ΔH)
-    r = a + (r0_mag - a) * cos(ΔH) + σ0 * √(-a) * sin(ΔH)
-    Ft = -√(-a)*√(μ) / (r*r0_mag) * sin(ΔH)
-    Gt = 1 - a/r * (1-cos(ΔH))
-  end
-
-  return [F*[x,y,z] + G*[dx,dy,dz]; Ft*[x,y,z] + Gt*[dx,dy,dz]]
-
-end