temporary point on the inner loop wrapper

julia/src/inner_loop/find_closest.jl

@@ -0,0 +1,30 @@
+using NLsolve
+
+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
+  return ans/n
+end
+
+function nlp_solve(start::Vector{Float64},
+                   final::Vector{Float64},
+                   craft::Sc,
+                   μ::Float64,
+                   t0::Float64,
+                   tf::Float64,
+                   x0::Matrix{Float64};
+                   tol=1e-6)
+
+  function f!(F,x)
+    F .= 0.0
+    F[1:6, 1] .= prop_nlsolve(tanh.(x), start, craft, μ, tf-t0) .- final
+  end
+
+  return nlsolve(f!, atanh.(x0), ftol=tol, autodiff=:forward, iterations=1_000)
+
+end

julia/src/inner_loop/inner_loop.jl

@@ -0,0 +1,25 @@
+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,
+                    RLA::Float64,
+                    DLA::Float64,
+                    phases::Vector{Phase})
+   
+  # 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
+      @assert phases[i].from_planet == phases[i-1].to_planet
+      @assert phases[i].v∞_outgoing == phases[i-1].v∞_incoming
+    end
+  end
+      
+  return RLA + DLA + phases[1].v∞_incoming
+end

julia/src/inner_loop/laguerre-conway.jl

@@ -0,0 +1,57 @@
+function laguerre_conway(state::Vector{T}, μ::Float64, time::Float64) where T
+
+  n = 5                               # Choose LaGuerre-Conway "n"
+  i = 0
+
+  r0 = state[1:3]                     # Are we in elliptical or hyperbolic orbit?
+  r0_mag = √(state[1]^2 + state[2]^2 + state[3]^2)
+  v0 = state[4:6]
+  v0_mag = √(state[4]^2 + state[5]^2 + state[6]^2)
+  σ0 = (r0 ⋅ v0)/√(μ)
+  a = 1 / ( 2/r0_mag - v0_mag^2/μ )
+  coeff = 1 - r0_mag/a
+
+  if a > 0    # Elliptical
+    ΔM = ΔE_new = √(μ) / sqrt(a^3) * time
+    Δ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 - n*F / ( dF + sign * √(abs((n-1)^2*dF^2 - n*(n-1)*F*d2F )))
+      i += 1
+    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) * time
+    ΔH = 0
+    ΔH_new = time < 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 - n*F / ( dF + sign * √(abs((n-1)^2*dF^2 - n*(n-1)*F*d2F )))
+      i += 1
+    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
+
+  if i > 100
+    throw(ErrorException("LaGuerre-Conway did not converge!"))
+  end
+
+  return [ F*state[1:3] + G*state[4:6]; Ft*state[1:3] + Gt*state[4:6]]
+
+end

julia/src/inner_loop/monotonic_basin_hopping.jl

@@ -0,0 +1,74 @@
+export mbh
+
+function pareto(α::Float64, n::Int)
+   s = rand((-1,1), (n,3))
+   r = rand(Float64, (n,3))
+   ϵ = 1e-10
+   return (s./ϵ) * (α - 1.0) ./ (ϵ ./ (ϵ .+ r)).^(-α)
+end
+
+function perturb(x::AbstractMatrix, n::Int)
+  ans =  x + pareto(1.01, n)
+  map!(elem -> elem > 1.0 ? 1.0 : elem, ans, ans)
+  map!(elem -> elem < -1.0 ? -1.0 : elem, ans, ans)
+  return ans
+end
+
+
+function new_x(n::Int)
+  2.0 * rand(Float64, (n,3)) .- 1.
+end
+
+function mbh(start::AbstractVector,
+             final::AbstractVector,
+             craft::Sc,
+             μ::AbstractFloat,
+             t0::AbstractFloat,
+             tf::AbstractFloat,
+             n::Int;
+             num_iters=50,
+             patience_level::Int=400,
+             tol=1e-6,
+             verbose=false)
+
+  archive = []
+
+  i = 0
+  if verbose println("Current Iteration") end
+  while true
+    i += 1
+    if verbose print("\r",i) end
+    impatience = 0
+    x_star = nlp_solve(start, final, craft, μ, t0, tf, new_x(n), tol=tol)
+    while converged(x_star) == false
+      x_star = nlp_solve(start, final, craft, μ, t0, tf, new_x(n), tol=tol)
+    end
+    if converged(x_star)
+      x_current = x_star
+      while impatience < patience_level
+        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_current = x_star
+          impatience = 0
+        else
+          impatience += 1
+        end
+      end
+      push!(archive, x_current)
+    end
+    if i >= num_iters break end
+  end
+  if verbose println() end
+
+  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))
+      best = candidate
+    end
+  end
+
+  return best, archive
+
+end

julia/src/inner_loop/phase.jl

@@ -0,0 +1,9 @@
+export Phase
+
+struct Phase
+  from_planet::String
+  to_planet::String
+  time_of_flight::Float64   # seconds
+  v∞_outgoing::Float64      # Km/s
+  v∞_incoming::Float64      # Km/s
+end

julia/src/inner_loop/propagator.jl

@@ -0,0 +1,232 @@
+export prop
+
+"""
+Maximum ΔV that a spacecraft can impulse for a given single time step
+"""
+function max_ΔV(duty_cycle::T,
+                num_thrusters::Int,
+                max_thrust::T,
+                tf::T,
+                t0::T,
+                mass::S) where {T <: Real, S <: 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},
+                  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)
+end
+
+
+"""
+This propagates over a given time period, with a certain number of intermediate steps
+"""
+function prop(ΔVs::Matrix{T},
+              state::Vector{Float64},
+              duty_cycle::Float64,
+              num_thrusters::Int,
+              max_thrust::Float64,
+              mass::Float64,
+              mass_flow_rate::Float64,
+              μ::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]
+
+  for i in 1:n
+    state, craft = 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)
+  end
+
+  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]
+  for i in 1:n
+    state, craft = prop_one(ΔVs[i,:], state, craft, μ, time/n)
+  end
+
+  return state
+
+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