diff --git a/julia/Manifest.toml b/julia/Manifest.toml index 1dbf6e2..0caedc8 100644 --- a/julia/Manifest.toml +++ b/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"] diff --git a/julia/src/Thesis.jl b/julia/src/Thesis.jl index 468c411..69c5dcd 100644 --- a/julia/src/Thesis.jl +++ b/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") diff --git a/julia/src/conversions.jl b/julia/src/conversions.jl index 6bc40d2..cb66a64 100644 --- a/julia/src/conversions.jl +++ b/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 diff --git a/julia/src/find_closest.jl b/julia/src/find_closest.jl index 9d496c7..518c95b 100644 --- a/julia/src/find_closest.jl +++ b/julia/src/find_closest.jl @@ -1,10 +1,14 @@ using NLsolve -export nlp_solve +export nlp_solve, mass_est -function treat_inputs(x::AbstractVector) - n::Int = length(x)/3 - reshape(x,(3,n))' +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}, @@ -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[7:3n] .= 0. + F .= 0.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 \ No newline at end of file diff --git a/julia/src/laguerre-conway.jl b/julia/src/laguerre-conway.jl index cb2cca3..7ca4503 100644 --- a/julia/src/laguerre-conway.jl +++ b/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 diff --git a/julia/src/monotonic_basin_hopping.jl b/julia/src/monotonic_basin_hopping.jl index cca4c08..dbb3c41 100644 --- a/julia/src/monotonic_basin_hopping.jl +++ b/julia/src/monotonic_basin_hopping.jl @@ -1,19 +1,20 @@ -function perturb(x::AbstractVector, n::Int) - perturb_vector = 0.02 * rand(Float64, (3n)) .- 0.01 - return x + perturb_vector +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 mass_better(x_star::AbstractVector, - x_current::AbstractVector, - start::AbstractVector, - final::AbstractVector, - craft::Sc, - μ::AbstractFloat, - t0::AbstractFloat, - tf::AbstractFloat) - mass_star = prop(treat_inputs(x_star), start, craft, μ, tf-t0)[2] - mass_current = prop(treat_inputs(x_current), start, craft, μ, tf-t0)[2] - return mass_star > mass_current +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, @@ -22,35 +23,51 @@ function mbh(start::AbstractVector, μ::AbstractFloat, t0::AbstractFloat, tf::AbstractFloat, - n::Int, - num_iters::Int=10, - tol=1e-6) - i::Int = 0 + n::Int; + num_iters=50, + patience_level::Int=400, + 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 - push!(archive, x_current) - - 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 = 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 - 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 \ No newline at end of file diff --git a/julia/src/plotting.jl b/julia/src/plotting.jl index e9fccf3..45d5fac 100644 --- a/julia/src/plotting.jl +++ b/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), diff --git a/julia/src/propagator.jl b/julia/src/propagator.jl index e6ee9f6..fbf0c85 100644 --- a/julia/src/propagator.jl +++ b/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}, - state::Vector{S}, +function prop_one(thrust_unit::Vector{<:Real}, + 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 - - # if mag > 1 || mag < 0 - # 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 + Δ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 @@ -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}, - state::Vector{Float64}, +function prop(ΔVs::Matrix{T}, + 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' - masses = craft.mass + 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) - states = [states; state'] - masses = [masses, craft.mass] + 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 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 \ No newline at end of file diff --git a/julia/test/find_closest.jl b/julia/test/find_closest.jl index f7870a3..f053130 100644 --- a/julia/test/find_closest.jl +++ b/julia/test/find_closest.jl @@ -1,7 +1,6 @@ @testset "Find Closest" begin - using NLsolve - using Thesis: treat_inputs + using NLsolve, PlotlyJS # 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)) - final = prop(ΔVs, start, sc, μs["Earth"], prop_time)[1][end,:] +# T_craft = hcat(repeat([0.6], n), repeat([0.], n), repeat([0.], n)) + 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) - result = nlp_solve(start, final, sc, μs["Earth"], 0.0, prop_time, x0) + Tx, Ty, Tz = conv_T(repeat([0.59], n), repeat([0.01], n), repeat([0.], n), + 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) - path3 = prop(zeros((100,3)), path2[end,:], sc, μs["Earth"], new_T)[1] + path2, mass, calc_final = prop(tanh.(result.zero), start, sc, μs["Earth"], prop_time) + 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 diff --git a/julia/test/monotonic_basin_hopping.jl b/julia/test/monotonic_basin_hopping.jl index 45e8373..452053d 100644 --- a/julia/test/monotonic_basin_hopping.jl +++ b/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 - n = 25 + prop_time = 0.75T + 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)) - final = prop(ΔVs, start, sc, μs["Earth"], prop_time)[1][end,:] +# T_craft = hcat(repeat([0.6], n), repeat([0.], n), repeat([0.], n)) + 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 - best, archive = mbh(start, final, sc, μs["Earth"], 0.0, prop_time, n) + # Find the best solution + 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 - @test converged(path) + transit, best_masses, calc_final = prop(tanh.(best.zero), start, sc, μs["Earth"], prop_time) + 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 diff --git a/julia/test/runtests.jl b/julia/test/runtests.jl index dc3d99d..9d25283 100644 --- a/julia/test/runtests.jl +++ b/julia/test/runtests.jl @@ -6,10 +6,10 @@ using Thesis # Tests @testset "All Tests" begin - include("spacecraft.jl") - include("laguerre-conway.jl") - include("propagator.jl") - include("plotting.jl") +# include("spacecraft.jl") +# include("laguerre-conway.jl") +# include("propagator.jl") +# include("plotting.jl") include("find_closest.jl") include("monotonic_basin_hopping.jl") end