Merge branch 'nlsolve' into 'main'

Got all the way to MBH with NLsolve

See merge request school/thesis!1

julia/Manifest.toml

@@ -5,9 +5,9 @@ uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
 
 [[ArrayInterface]]
 deps = ["IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
-git-tree-sha1 = "baf4ef9082070477046bd98306952292bfcb0af9"
+git-tree-sha1 = "85d03b60274807181bae7549bb22b2204b6e5a0e"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "3.1.25"
+version = "3.1.30"
 
 [[Artifacts]]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
@@ -416,9 +416,9 @@ version = "1.6.1"
 
 [[Static]]
 deps = ["IfElse"]
-git-tree-sha1 = "62701892d172a2fa41a1f829f66d2b0db94a9a63"
+git-tree-sha1 = "854b024a4a81b05c0792a4b45293b85db228bd27"
 uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
-version = "0.3.0"
+version = "0.3.1"
 
 [[StaticArrays]]
 deps = ["LinearAlgebra", "Random", "Statistics"]

julia/src/Thesis.jl

@@ -3,8 +3,8 @@ module Thesis
   using LinearAlgebra, ForwardDiff, PlotlyJS
 
   include("./constants.jl")
-  include("./conversions.jl")
   include("./spacecraft.jl")
+  include("./conversions.jl")
   include("./plotting.jl")
   include("./laguerre-conway.jl")
   include("./propagator.jl")

julia/src/conversions.jl

@@ -1,4 +1,4 @@
-export oe_to_rθh, rθh_to_xyz, oe_to_xyz, xyz_to_oe
+export oe_to_rθh, rθh_to_xyz, oe_to_xyz, xyz_to_oe, conv_T
 
 function oe_to_rθh(oe::Vector,μ::Real) :: Vector
 
@@ -68,3 +68,51 @@ function xyz_to_oe(cart_vec::Vector,μ::Real)
   return [a,e,i,Ω,ω,ν]
 
 end
+
+"""
+Converts a series of thrust vectors from R,Θ,H frame to cartesian
+"""
+function conv_T(Tm::Vector{Float64},
+                Ta::Vector{Float64},
+                Tb::Vector{Float64},
+                init_state::Vector{Float64},
+                m::Float64,
+                craft::Sc,
+                time::Float64,
+                μ::Float64)
+
+  Txs = Float64[]
+  Tys = Float64[]
+  Tzs = Float64[]
+  n::Int = length(Tm)
+
+  for i in 1:n
+    mag, α, β = Tm[i], Ta[i], Tb[i]
+    if mag > 1 || mag < 0
+      @error "ΔV input is too high: $mag"
+    elseif α > π || α < -π
+      @error "α angle is incorrect: $α"
+    elseif β > π/2 || β < -π/2
+      @error "β angle is incorrect: $β"
+    end
+  end
+
+  state = init_state
+  for i in 1:n
+    mag, α, β = Tm[i], Ta[i], Tb[i]
+    thrust_rθh = mag * [cos(β)*sin(α), cos(β)*cos(α), sin(β)]
+    _,_,i,Ω,ω,ν = xyz_to_oe(state, μ)
+    θ = ω+ν
+    cΩ,sΩ,ci,si,cθ,sθ = cos(Ω),sin(Ω),cos(i),sin(i),cos(θ),sin(θ)
+    DCM = [cΩ*cθ-sΩ*ci*sθ -cΩ*sθ-sΩ*ci*cθ sΩ*si;
+           sΩ*cθ+cΩ*ci*sθ -sΩ*sθ+cΩ*ci*cθ -cΩ*si;
+           si*sθ          si*cθ           ci      ]
+    Tx, Ty, Tz = DCM*thrust_rθh
+
+    state = prop_one([Tx, Ty, Tz], state, craft, μ, time/n)[1]
+    push!(Txs, Tx)
+    push!(Tys, Ty)
+    push!(Tzs, Tz)
+  end
+  return Txs, Tys, Tzs
+end

julia/src/find_closest.jl

@@ -1,10 +1,14 @@
 using NLsolve
 
-export nlp_solve
+export nlp_solve, mass_est
 
-function treat_inputs(x::AbstractVector)
+function mass_est(T)
-  n::Int = length(x)/3
+  ans = 0
-  reshape(x,(3,n))'
+  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},
@@ -13,15 +17,14 @@ function nlp_solve(start::Vector{Float64},
                    μ::Float64,
                    t0::Float64,
                    tf::Float64,
-                   x0::AbstractVector;
+                   x0::Matrix{Float64};
                    tol=1e-6)
 
-  n::Int = length(x0)/3
   function f!(F,x)
-    F[1:6] .= prop(treat_inputs(x), start, craft, μ, tf-t0)[1][end,:] - final
+    F .= 0.0
-    F[7:3n] .= 0.
+    F[1:6, 1] .= prop_nlsolve(tanh.(x), start, craft, μ, tf-t0) .- final
   end
 
-  return nlsolve(f!, x0, ftol=tol, autodiff=:forward, iterations=1_000)
+  return nlsolve(f!, atanh.(x0), ftol=tol, autodiff=:forward, iterations=1_000)
 
 end

julia/src/laguerre-conway.jl

@@ -1,12 +1,12 @@
-function laguerre_conway(state::Vector{T}, μ::Float64, time::Float64) where T <: Real
+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 = norm(r0)
+  r0_mag = √(state[1]^2 + state[2]^2 + state[3]^2)
   v0 = state[4:6]
-  v0_mag = norm(v0)
+  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

julia/src/monotonic_basin_hopping.jl

@@ -1,19 +1,20 @@
-function perturb(x::AbstractVector, n::Int)
+function pareto(α::Float64, n::Int)
-  perturb_vector = 0.02 * rand(Float64, (3n)) .- 0.01
+   s = rand((-1,1), (n,3))
-  return x + perturb_vector
+   r = rand(Float64, (n,3))
+   ϵ = 1e-10
+   return (s./ϵ) * (α - 1.0) ./ (ϵ ./ (ϵ .+ r)).^(-α)
 end
 
-function mass_better(x_star::AbstractVector,
+function perturb(x::AbstractMatrix, n::Int)
-                     x_current::AbstractVector,
+  ans =  x + pareto(1.01, n)
-                     start::AbstractVector,
+  map!(elem -> elem > 1.0 ? 1.0 : elem, ans, ans)
-                     final::AbstractVector,
+  map!(elem -> elem < -1.0 ? -1.0 : elem, ans, ans)
-                     craft::Sc,
+  return ans
-                     μ::AbstractFloat,
+end
-                     t0::AbstractFloat,
+
-                     tf::AbstractFloat)
+
-  mass_star = prop(treat_inputs(x_star), start, craft, μ, tf-t0)[2]
+function new_x(n::Int)
-  mass_current = prop(treat_inputs(x_current), start, craft, μ, tf-t0)[2]
+  2.0 * rand(Float64, (n,3)) .- 1.
-  return mass_star > mass_current
 end
 
 function mbh(start::AbstractVector,
@@ -22,35 +23,51 @@ function mbh(start::AbstractVector,
              μ::AbstractFloat,
              t0::AbstractFloat,
              tf::AbstractFloat,
-             n::Int,
+             n::Int;
-             num_iters::Int=10,
+             num_iters=50,
-             tol=1e-6)
+             patience_level::Int=400,
-  i::Int = 0
+             tol=1e-6,
+             verbose=false)
+
   archive = []
-  x_star = nlp_solve(start, final, craft, μ, t0, tf, rand(Float64,(3n)), tol=tol)
-  while converged(x_star) == false
-    x_star = nlp_solve(start, final, craft, μ, t0, tf, rand(Float64,(3n)), tol=tol)
-  end
 
-  x_current = x_star
+  i = 0
-  push!(archive, x_current)
+  if verbose println("Current Iteration") end
-
+  while true
-  while i < num_iters
-    x_star = nlp_solve(start, final, craft, μ, t0, tf, perturb(x_current.zero,n), tol=tol)
-    if converged(x_star) && mass_better(x_star.zero, x_current.zero, start, final, craft, μ, t0, tf)
-      x_current = x_star
-      push!(archive, x_star)
-    else
-      while converged(x_star) == false
-        x_star = nlp_solve(start, final, craft, μ, t0, tf, rand(Float64,(3n)), tol=tol)
-      end
-      if mass_better(x_star.zero, x_current.zero, start, final, craft, μ, t0, tf)
-        x_current = x_star
-        push!(archive, x_star)
-      end
-    end
     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
 
-  return x_current, archive
+  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/plotting.jl

@@ -16,7 +16,7 @@ function get_true_max(mat::Vector{Array{Float64,2}})
   return maximum(abs.(flatten(mat)))
 end
 
-function plot_orbits(paths::Vector{Array{Float64,2}};
+function plot_orbits(paths::Vector{Vector{Vector{Float64}}};
                      primary::String="Earth",
                      plot_theme::Symbol=:juno,
                      labels::Vector{String}=Vector{String}(),
@@ -34,13 +34,14 @@ function plot_orbits(paths::Vector{Array{Float64,2}};
 
   t1 = []
   for i = 1:length(paths)
-    path = [ x for x in paths[i] ]
+    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]),
+    push!(t1, scatter3d(;x=(path[1]),y=(path[2]),z=(path[3]),
                          mode="lines", name=label, line_color=color, line_width=3))
   end
-  limit = max(maximum(abs.(flatten(paths))),maximum(abs.(flatten(ps)))) * 1.1
+  limit = max(maximum(abs.(flatten(flatten(paths)))),
+              maximum(abs.(flatten(ps)))) * 1.1
 
   t2 = surface(;x=(x_p),
                 y=(y_p),

julia/src/propagator.jl

@@ -16,38 +16,19 @@ 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(ΔV::Vector{T},
+function prop_one(thrust_unit::Vector{<:Real},
-                  state::Vector{S},
+                  state::Vector{<:Real},
                   duty_cycle::Float64,
                   num_thrusters::Int,
                   max_thrust::Float64,
-                  mass::S,
+                  mass::T,
                   mass_flow_rate::Float64,
                   μ::Float64,
-                  time::Float64) where {T <: Real, S <: Real}
+                  time::Float64) where T<:Real
 
-  mag, α, β = ΔV
+  ΔV = max_ΔV(duty_cycle, num_thrusters, max_thrust, time, 0., mass) * thrust_unit
-
+  halfway = laguerre_conway(state, μ, time/2) + [0., 0., 0., ΔV...]
-  # if mag > 1 || mag < 0
+  return laguerre_conway(halfway, μ, time/2), mass - mass_flow_rate*norm(thrust_unit)*time
-    # throw(ErrorException("ΔV input is too high: $mag"))
-  # elseif α > π || α < -π
-    # throw(ErrorException("α angle is incorrect: $α"))
-  # elseif β > π/2 || β < -π/2
-    # throw(ErrorException("β angle is incorrect: $β"))
-  # end
-
-  thrust_rθh = mag * [cos(β)*sin(α), cos(β)*cos(α), sin(β)]
-  a,e,i,Ω,ω,ν = xyz_to_oe(state, μ)
-  θ = ω+ν
-  cΩ,sΩ,ci,si,cθ,sθ = cos(Ω),sin(Ω),cos(i),sin(i),cos(θ),sin(θ)
-  DCM = [cΩ*cθ-sΩ*ci*sθ -cΩ*sθ-sΩ*ci*cθ sΩ*si;
-  sΩ*cθ+cΩ*ci*sθ -sΩ*sθ+cΩ*ci*cθ -cΩ*si;
-  si*sθ si*cθ ci]
-  ΔV = DCM*thrust_rθh
-
-  thrust = max_ΔV(duty_cycle, num_thrusters, max_thrust, time, 0., mass) * ΔV
-  halfway = laguerre_conway(state, μ, time/2) + [0., 0., 0., thrust[1], thrust[2], thrust[3]]
-  return laguerre_conway(halfway, μ, time/2), mass - mass_flow_rate*norm(ΔV)*time
 
 end
 
@@ -64,6 +45,7 @@ function prop_one(ΔV_unit::Vector{T},
   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
 """
@@ -92,24 +74,159 @@ end
 """
 The same function, using Scs
 """
-function prop(ΔVs::AbstractArray{T},
+function prop(ΔVs::Matrix{T},
-              state::Vector{Float64},
+              state::Vector{S},
               craft::Sc,
               μ::Float64,
-              time::Float64) where T <: Real
+              time::Float64) where {T <: Real, S <: Real}
 
   if size(ΔVs)[2] != 3 throw(ErrorException("ΔV input is wrong size")) end
   n = size(ΔVs)[1]
 
-  states = state'
+  x_states = [state[1]]
-  masses = craft.mass
+  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)
-    states = [states; state']
+    push!(x_states, state[1])
-    masses = [masses, craft.mass]
+    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 states, masses
+  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

julia/test/find_closest.jl

@@ -1,7 +1,6 @@
 @testset "Find Closest" begin
 
-   using NLsolve
+  using NLsolve, PlotlyJS
-   using Thesis: treat_inputs
 
   # Initial Setup
   sc = Sc("test")
@@ -10,30 +9,41 @@
   i = rand(0.01:0.01:π/6)
   T = 2π*√(a^3/μs["Earth"])
   prop_time = 2T
-  n = 30
+  n = 20
 
   # A simple orbit raising
   start = oe_to_xyz([ a, e, i, 0., 0., 0. ], μs["Earth"])
-  ΔVs = repeat([0.6, 0., 0.]', outer=(n,1))
+#   T_craft = hcat(repeat([0.6], n), repeat([0.], n), repeat([0.], n))
-  final = prop(ΔVs, start, sc, μs["Earth"], prop_time)[1][end,:]
+  Tx, Ty, Tz = conv_T(repeat([0.6], n), repeat([0.], n), repeat([0.], n),
+                      start,
+                      sc.mass,
+                      sc,
+                      prop_time,
+                      μs["Earth"])
+  final = prop(hcat(Tx, Ty, Tz), start, sc, μs["Earth"], prop_time)[3]
   new_T = 2π*√(xyz_to_oe(final, μs["Earth"])[1]^3/μs["Earth"])
 
   # This should be close enough to 0.6
-  x0 = repeat([0.55, 0., 0.], n)
+  Tx, Ty, Tz = conv_T(repeat([0.59], n), repeat([0.01], n), repeat([0.], n),
-  result = nlp_solve(start, final, sc, μs["Earth"], 0.0, prop_time, x0)
+                      start,
+                      sc.mass,
+                      sc,
+                      prop_time,
+                      μs["Earth"])
+  result = nlp_solve(start, final, sc, μs["Earth"], 0.0, prop_time, hcat(Tx, Ty, Tz))
 
   # Test and plot
   @test converged(result)
   path1 = prop(zeros((100,3)), start, sc, μs["Earth"], T)[1]
-  path2, mass = prop(treat_inputs(result.zero), start, sc, μs["Earth"], prop_time)
+  path2, mass, calc_final = prop(tanh.(result.zero), start, sc, μs["Earth"], prop_time)
-  path3 = prop(zeros((100,3)), path2[end,:], sc, μs["Earth"], new_T)[1]
+  path3 = prop(zeros((100,3)), calc_final, sc, μs["Earth"], new_T)[1]
   path4 = prop(zeros((100,3)), final, sc, μs["Earth"], new_T)[1]
   savefig(plot_orbits([path1, path2, path3, path4],
                       labels=["initial", "transit", "after transit", "final"],
                       colors=["#FFFFFF","#FF4444","#44FF44","#4444FF"]),
-          "../plots/find_closest_test.html")
+                      "../plots/find_closest_test.html")
   if converged(result)
-    @test norm(path2[end,:] - final) < 1e-4
+    @test norm(calc_final - final) < 1e-4
   end
 
 end

julia/test/monotonic_basin_hopping.jl

@@ -8,22 +8,63 @@
   e = rand(0.01:0.01:0.5)
   i = rand(0.01:0.01:π/6)
   T = 2π*√(a^3/μs["Earth"])
-  prop_time = 2T
+  prop_time = 0.75T
-  n = 25
+  n = 10
 
   # A simple orbit raising
   start = oe_to_xyz([ a, e, i, 0., 0., 0. ], μs["Earth"])
-  ΔVs = repeat([0.6, 0., 0.]', outer=(n,1))
+#   T_craft = hcat(repeat([0.6], n), repeat([0.], n), repeat([0.], n))
-  final = prop(ΔVs, start, sc, μs["Earth"], prop_time)[1][end,:]
+  Tx, Ty, Tz = conv_T(repeat([0.8], n), repeat([0.], n), repeat([0.], n),
+                      start,
+                      sc.mass,
+                      sc,
+                      prop_time,
+                      μs["Earth"])
+  nominal_path, normal_mass, final = prop(hcat(Tx, Ty, Tz), start, sc, μs["Earth"], prop_time)
   new_T = 2π*√(xyz_to_oe(final, μs["Earth"])[1]^3/μs["Earth"])
 
-  # This should be close enough to 0.6
+  # Find the best solution
-  best, archive = mbh(start, final, sc, μs["Earth"], 0.0, prop_time, n)
+  best, archive = mbh(start,
+                      final,
+                      sc,
+                      μs["Earth"],
+                      0.0,
+                      prop_time,
+                      n,
+                      num_iters=5,
+                      patience_level=50,
+                      verbose=true)
 
   # Test and plot
   @test converged(best)
-  for path in archive
+  transit, best_masses, calc_final = prop(tanh.(best.zero), start, sc, μs["Earth"], prop_time)
-   @test converged(path)
+  initial_path = prop(zeros((100,3)), start, sc, μs["Earth"], T)[1]
+  after_transit = prop(zeros((100,3)), calc_final, sc, μs["Earth"], new_T)[1]
+  final_path = prop(zeros((100,3)), final, sc, μs["Earth"], new_T)[1]
+  savefig(plot_orbits([initial_path, nominal_path, final_path],
+                      labels=["initial", "nominal transit", "final"],
+                      colors=["#FF4444","#44FF44","#4444FF"]),
+                      "../plots/mbh_nominal.html")
+  savefig(plot_orbits([initial_path, transit, after_transit, final_path],
+                      labels=["initial", "transit", "after transit", "final"],
+                      colors=["#FFFFFF", "#FF4444","#44FF44","#4444FF"]),
+                      "../plots/mbh_best.html")
+  i = 0
+  best_mass = best_masses[end]
+  nominal_mass = normal_mass[end]
+  masses = []
+  for candidate in archive
+    @test converged(candidate)
+    path2, cand_ms, calc_final = prop(tanh.(candidate.zero), start, sc, μs["Earth"], prop_time)
+    push!(masses, cand_ms[end])
+    @test norm(calc_final - final) < 1e-4
   end
+  @test best_mass == maximum(masses)
+
+  # This won't always work since the test is reduced in fidelity,
+  # but hopefully will usually work:
+  @test (sc.mass - best_mass) < 1.1 * (sc.mass - nominal_mass)
+  @show best_mass
+  @show nominal_mass
 
 end

julia/test/runtests.jl

@@ -6,10 +6,10 @@ using Thesis
 
 # Tests
 @testset "All Tests" begin
-  include("spacecraft.jl")
+#   include("spacecraft.jl")
-  include("laguerre-conway.jl")
+#   include("laguerre-conway.jl")
-  include("propagator.jl")
+#   include("propagator.jl")
-  include("plotting.jl")
+#   include("plotting.jl")
   include("find_closest.jl")
   include("monotonic_basin_hopping.jl")
 end