diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index d6e513a..8bad5ca 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -12,12 +12,7 @@ unit-test-job: - julia --project=julia/Project.toml -E 'using Pkg; Pkg.test()' artifacts: paths: - - julia/plots/plot_test.html - - julia/plots/find_closest_test.html - - julia/plots/mbh_nominal.html - - julia/plots/mbh_best.html - - julia/plots/mbh_sun_initial.html - - julia/plots/mbh_sun_solved.html - - julia/plots/inner_loop_before.html - - julia/plots/inner_loop_after.html - expire_in: 1 week + - julia/plots/plot_test_earth.html + - julia/plots/plot_test_sun.html + - julia/plots/nlp_test.html + expire_in: 10 weeks diff --git a/julia/src/Thesis.jl b/julia/src/Thesis.jl index e85b4ca..a96a9ca 100644 --- a/julia/src/Thesis.jl +++ b/julia/src/Thesis.jl @@ -1,25 +1,30 @@ +using SPICE + +try + furnsh("../../spice_files/naif0012.tls") + furnsh("../../spice_files/de430.bsp") +catch + furnsh("spice_files/naif0012.tls") + furnsh("spice_files/de430.bsp") +end + module Thesis - using LinearAlgebra, ForwardDiff, PlotlyJS, SPICE, Distributed + using LinearAlgebra + using ForwardDiff + using PlotlyJS + using Distributed - try - furnsh("../../SPICE/naif0012.tls") - furnsh("../../SPICE/de430.bsp") - catch - furnsh("SPICE/naif0012.tls") - furnsh("SPICE/de430.bsp") - end + include("./errors.jl") include("./constants.jl") include("./spacecraft.jl") include("./conversions.jl") include("./plotting.jl") include("./inner_loop/laguerre-conway.jl") include("./inner_loop/propagator.jl") - include("./inner_loop/find_closest.jl") - include("./inner_loop/monotonic_basin_hopping.jl") include("./inner_loop/phase.jl") - include("./inner_loop/inner_loop.jl") - include("./outer_loop.jl") + # include("./inner_loop/monotonic_basin_hopping.jl") + # include("./outer_loop.jl") end diff --git a/julia/src/constants.jl b/julia/src/constants.jl index 1d54d7c..3156e6d 100644 --- a/julia/src/constants.jl +++ b/julia/src/constants.jl @@ -2,113 +2,36 @@ # DEFINING CONSTANTS # ----------------------------------------------------------------------------- -export μs, G, GMs, μ, rs, as, es, AU, ids +export Body, Sun, Mercury, Venus, Earth, Moon, Mars +export Jupiter, Saturn, Uranus, Neptune, Pluto +export G, AU, init_STM, hour, day, year, second +export Pathlist -# Gravitational Constants - μs = Dict( - "Sun" => 1.32712440018e11, - "Mercury" => 2.2032e4, - "Venus" => 3.257e5, - "Earth" => 3.986004415e5, - "Moon" => 4.902799e3, - "Mars" => 4.305e4, - "Jupiter" => 1.266865361e8, - "Saturn" => 3.794e7, - "Uranus" => 5.794e6, - "Neptune" => 6.809e6, - "Pluto" => 9e2) +struct Body + μ::Float64 + r::Float64 # radius + color::String + id::Int # SPICE id +end - G = 6.67430e-20 +const Sun = Body(1.32712440018e11, 696000., "Electric", 10) +const Mercury = Body(2.2032e4, 2439., "heat", 1) +const Venus = Body(3.257e5, 6052., "turbid", 2) +const Earth = Body(3.986004415e5, 6378.1363, "Blues", 399) +const Moon = Body(4.902799e3, 1738., "Greys", 301) +const Mars = Body(4.305e4, 3397.2, "Reds", 4) +const Jupiter = Body(1.266865361e8, 71492., "solar", 5) +const Saturn = Body(3.794e7, 60268., "turbid", 6) +const Uranus = Body(5.794e6, 25559., "haline", 7) +const Neptune = Body(6.809e6, 24764., "ice", 8) +const Pluto = Body(9e2, 1151., "matter", 9) - function μ(m1::Float64, m2::Float64) - return m2/(m1+m2) - end +const G = 6.67430e-20 #universal gravity parameter +const AU = 149597870.691 #km +const init_STM = vec(Matrix{Float64}(I,6,6)) +const second = 1. +const hour = 3600. +const day = 86400. +const year = 365 * day - function μ(GM1::Float64, GM2::Float64, Grav::Float64) - return μ(GM1/Grav, GM2/Grav) - end - - function μ(primary::String, secondary::String) - return μ(GMs[primary]/G, GMs[secondary]/G) - end - - const GMs = Dict( - "Sun" => 132712440041.93938, - "Earth" => 398600.435436, - "Moon" => 4902.800066) - -# Radii - const rs = Dict( - "Sun" => 696000., - "Mercury" => 2439., - "Venus" => 6052., - "Earth" => 6378.1363, - "Moon" => 1738., - "Mars" => 3397.2, - "Jupiter" => 71492., - "Saturn" => 60268., - "Uranus" => 25559., - "Neptune" => 24764., - "Pluto" => 1151.) - -# Semi Major Axes - const as = Dict( - "Mercury" => 57909083., - "Venus" => 108208601., - "Earth" => 149598023., - "Moon" => 384400., - "Mars" => 227939186., - "Jupiter" => 778298361., - "Saturn" => 1429394133., - "Uranus" => 2875038615., - "Neptune" => 4504449769., - "Pluto" => 5915799000.) - -# Eccentricities - const es = Dict( - "Earth" => 0.016708617, - "Moon" => 0.0549) - -# J2 for basic oblateness - const j2s = Dict( - "Mercury" => 0.00006, - "Venus" => 0.000027, - "Earth" => 0.0010826269, - "Moon" => 0.0002027, - "Mars" => 0.001964, - "Jupiter" => 0.01475, - "Saturn" => 0.01645, - "Uranus" => 0.012, - "Neptune" => 0.004, - "Pluto" => 0.) - -# These are just the colors for plots - const p_colors = Dict( - "Sun" => "Electric", - "Mercury" => "heat", - "Venus" => "turbid", - "Earth" => "Blues", - "Moon" => "Greys", - "Mars" => "Reds", - "Jupiter" => "solar", - "Saturn" => "turbid", - "Uranus" => "haline", - "Neptune" => "ice", - "Pluto" => "matter") - - const ids = Dict( - "Sun" => 10, - "Mercury" => 1, - "Venus" => 2, - "Earth" => 399, - "Moon" => 301, - "Mars" => 4, - "Jupiter" => 5, - "Saturn" => 6, - "Uranus" => 7, - "Neptune" => 8, - "Pluto" => 9, - ) - - const AU = 149597870.691 #km - const init_STM = vec(Matrix{Float64}(I,6,6)) +Pathlist = Vector{Vector{Vector{Float64}}} diff --git a/julia/src/conversions.jl b/julia/src/conversions.jl index 115e799..ee389db 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, conv_T +export oe_to_rθh, rθh_to_xyz, oe_to_xyz, xyz_to_oe, conv_T, spiral, gen_orbit function oe_to_rθh(oe::Vector,μ::Real) :: Vector @@ -69,6 +69,30 @@ function xyz_to_oe(cart_vec::Vector,μ::Real) end +""" +A convenience function for generating start conditions from orbital elements +Inputs: a body, a period, and a mass +Output: a random reasonable orbit +""" +function gen_orbit(T::Float64, mass::Float64, primary::Body=Sun) + μ = primary.μ + i = rand(0.0:0.01:0.4999π) + θ = rand(0.0:0.01:2π) + i = 0 + while true + i += 1 + e = rand(0.0:0.01:0.5) + a = ∛(μ * ( T/2π )^2 ) + a*(1-e) < 1.1primary.r || return [ oe_to_xyz([ a, e, i, 0., 0., θ ], μ); mass ] + i < 100 || throw(GenOrbit_Error) + end +end + +""" +A convenience function for generating spiral trajectories +""" +spiral(mag,n,init,sc,time,primary=Sun) = conv_T(fill(mag, n), zeros(n), zeros(n), init, sc, time, primary) + """ Converts a series of thrust vectors from R,Θ,H frame to cartesian """ @@ -78,7 +102,7 @@ function conv_T(Tm::Vector{Float64}, init_state::Vector{Float64}, craft::Sc, time::Float64, - μ::Float64) + primary::Body=Sun) Txs = Float64[] Tys = Float64[] @@ -100,7 +124,7 @@ function conv_T(Tm::Vector{Float64}, 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, μ) + _,_,i,Ω,ω,ν = xyz_to_oe(state, primary.μ) θ = ω+ν 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; @@ -108,10 +132,10 @@ function conv_T(Tm::Vector{Float64}, si*sθ si*cθ ci ] Tx, Ty, Tz = DCM*thrust_rθh - state = prop_one([Tx, Ty, Tz], state, copy(craft), μ, time/n) + state = prop_one([Tx, Ty, Tz], state, craft, time/n, primary) push!(Txs, Tx) push!(Tys, Ty) push!(Tzs, Tz) end - return Txs, Tys, Tzs + return hcat(Txs, Tys, Tzs) end diff --git a/julia/src/errors.jl b/julia/src/errors.jl new file mode 100644 index 0000000..fd41b51 --- /dev/null +++ b/julia/src/errors.jl @@ -0,0 +1,16 @@ +struct LaGuerreConway_Error <: Exception end + +struct ΔVsize_Error <: Exception end + +struct GenOrbit_Error <: Exception end +Base.showerror(io::IO, e::GenOrbit_Error) = print(io, "Infinite Loop trying to generate the init orbit") + +struct PropOne_Error <: Exception + ΔV_unit::AbstractVector +end +Base.showerror(io::IO, e::PropOne_Error) = print(io, "tried to prop a unit ΔV of: ", e.ΔV_unit) + +struct Mass_Error{T} <: Exception where T <: Real + mass::T +end +Base.showerror(io::IO, e::Mass_Error) = print(io, "Mass (", e.mass, ") got too low in propagation") diff --git a/julia/src/inner_loop/find_closest.jl b/julia/src/inner_loop/find_closest.jl deleted file mode 100644 index 06b447b..0000000 --- a/julia/src/inner_loop/find_closest.jl +++ /dev/null @@ -1,45 +0,0 @@ -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 - -struct Result - converged::Bool - zero::Matrix{Float64} -end - -function nlp_solve(start::Vector{Float64}, - final::Vector{Float64}, - craft::Sc, - μ::Float64, - t0::Float64, - tf::Float64, - x0::Matrix{Float64}; - tol=1e-6, - num_iters=1_000) - - function f!(F,x) - 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 - - 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 - - return result - -end diff --git a/julia/src/inner_loop/laguerre-conway.jl b/julia/src/inner_loop/laguerre-conway.jl index 0220bf8..b14fdce 100644 --- a/julia/src/inner_loop/laguerre-conway.jl +++ b/julia/src/inner_loop/laguerre-conway.jl @@ -1,5 +1,6 @@ -function laguerre_conway(state::Vector{<:Real}, μ::Float64, time::Float64) +function laguerre_conway(state::Vector{<:Real}, time::Float64, primary::Body=Sun) + μ = primary.μ n = 5 # Choose LaGuerre-Conway "n" i = 0 @@ -22,7 +23,7 @@ function laguerre_conway(state::Vector{<:Real}, μ::Float64, time::Float64) 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 - if i > 100 throw(ErrorException("LaGuerre-Conway did not converge!")) end + if i > 100 throw(LaGuerreConway_Error("LaGuerre-Conway did not converge!")) end end F = 1 - a/r0_mag * (1-cos(ΔE)) G = a * σ0/ √(μ) * (1-cos(ΔE)) + r0_mag * √(a) / √(μ) * sin(ΔE) @@ -41,7 +42,7 @@ function laguerre_conway(state::Vector{<:Real}, μ::Float64, time::Float64) 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 - if i > 100 throw(ErrorException("LaGuerre-Conway did not converge!")) end + if i > 100 throw(LaGuerreConway_Error("LaGuerre-Conway did not converge!")) end end F = 1 - a/r0_mag * (1-cos(ΔH)) G = a * σ0/ √(μ) * (1-cos(ΔH)) + r0_mag * √(-a) / √(μ) * sin(ΔH) diff --git a/julia/src/inner_loop/monotonic_basin_hopping.jl b/julia/src/inner_loop/monotonic_basin_hopping.jl index 926dfbe..e4d3cf6 100644 --- a/julia/src/inner_loop/monotonic_basin_hopping.jl +++ b/julia/src/inner_loop/monotonic_basin_hopping.jl @@ -1,5 +1,8 @@ export mbh +""" +Generates n pareto-distributed random numbers +""" function pareto(α::Float64, n::Int) s = rand((-1,1), (n,3)) r = rand(Float64, (n,3)) @@ -7,6 +10,10 @@ function pareto(α::Float64, n::Int) return (s./ϵ) * (α - 1.0) ./ (ϵ ./ (ϵ .+ r)).^(-α) end +""" +Perturbs the monotonic basin hopping decision vector +TODO: This needs to be updated +""" function perturb(x::AbstractMatrix, n::Int) ans = x + pareto(1.01, n) map!(elem -> elem > 1.0 ? 1.0 : elem, ans, ans) @@ -14,11 +21,13 @@ function perturb(x::AbstractMatrix, n::Int) return ans end - function new_x(n::Int) 2.0 * rand(Float64, (n,3)) .- 1. end +function mbh(flybys::Vector{Planet}) +end + function mbh(start::AbstractVector, final::AbstractVector, craft::Sc, diff --git a/julia/src/inner_loop/phase.jl b/julia/src/inner_loop/phase.jl index ae809a2..a8c364f 100644 --- a/julia/src/inner_loop/phase.jl +++ b/julia/src/inner_loop/phase.jl @@ -1,9 +1,38 @@ -export Phase +using NLsolve + +export solve_phase + +""" +This function will take a single phase (so an initial state, and a final state) and an initial guess +to the thrust profile and use an NLP solver to find the nearest thrust profile to that initial guess +that satisfies the final state condition +""" +function solve_phase( start::Vector{Float64}, + final::Vector{Float64}, + craft::Sc, + tof::Float64, + x0::Matrix{Float64}, + primary::Body=Sun; + tol=1e-6, + num_iters=1_000 ) + + function f!(F,x) + try + F .= 0.0 + F[1:6, 1] .= prop(tanh.(x), start, craft, tof, primary)[2][1:6] .- final[1:6] + catch e + # If the error is due to something natural, just imply a penalty + if isa(Mass_Error, e) || isa(PropOne_Error, e) + F .= 10000000.0 + else + rethrow() + end + end + end + + result = nlsolve(f!, atanh.(x0), ftol=tol, autodiff=:forward, iterations=num_iters) + result.zero = tanh.(result.zero) + + return result -struct Phase - from_planet::String - to_planet::String - time_of_flight::Float64 # seconds - v∞_outgoing::Vector{Float64} # Km/s - v∞_incoming::Vector{Float64} # Km/s end diff --git a/julia/src/inner_loop/propagator.jl b/julia/src/inner_loop/propagator.jl index 278ebe5..763a173 100644 --- a/julia/src/inner_loop/propagator.jl +++ b/julia/src/inner_loop/propagator.jl @@ -18,18 +18,15 @@ A convenience function for using spacecraft. Note that this function outputs a s function prop_one(ΔV_unit::Vector{<:Real}, state::Vector{<:Real}, craft::Sc, - μ::Float64, - time::Float64) + time::Float64, + primary::Body=Sun) for direction in ΔV_unit - if abs(direction) > 1.0 - println(direction) - error("ΔV is impossibly high") - end + abs(direction) <= 1.0 || throw(PropOne_Error(ΔV_unit)) 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) + halfway = laguerre_conway(state, time/2, primary) + [zeros(3); ΔV] + final = laguerre_conway(halfway, time/2, primary) return [final; state[7] - craft.mass_flow_rate*norm(ΔV_unit)*time] end @@ -40,10 +37,10 @@ The propagator function function prop(ΔVs::Matrix{T}, state::Vector{Float64}, craft::Sc, - μ::Float64, - time::Float64) where T <: Real + time::Float64, + primary::Body=Sun) where T <: Real - if size(ΔVs)[2] != 3 throw(ErrorException("ΔV input is wrong size")) end + size(ΔVs)[2] == 3 || throw(ΔVsize_Error()) n = size(ΔVs)[1] x_states = Vector{T}() @@ -63,7 +60,7 @@ function prop(ΔVs::Matrix{T}, push!(masses, state[7]) for i in 1:n - state = prop_one(ΔVs[i,:], state, craft, μ, time/n) + state = prop_one(ΔVs[i,:], state, craft, time/n, primary) push!(x_states, state[1]) push!(y_states, state[2]) push!(z_states, state[3]) @@ -71,12 +68,14 @@ function prop(ΔVs::Matrix{T}, push!(dy_states, state[5]) push!(dz_states, state[6]) push!(masses, state[7]) - if state[7] < craft.dry_mass - println(state[7]) - error("Mass is too low") - end + state[7] >= craft.dry_mass || throw(Mass_Error(state[7])) end return [x_states, y_states, z_states, dx_states, dy_states, dz_states, masses], state end + +""" +Convenience function for propagating a state with no thrust +""" +prop(x::Vector{Float64}, t::Float64, p::Body=Sun) = prop(zeros(1000,3), x, no_thrust, t, p)[1] diff --git a/julia/src/outer_loop.jl b/julia/src/outer_loop.jl index 37f2ea2..6a08d93 100644 --- a/julia/src/outer_loop.jl +++ b/julia/src/outer_loop.jl @@ -61,3 +61,9 @@ function gen_decision_vector(launch_range::Vector{DateTime}, end +""" +This is the binary crossover function, implemented as detailed in Englander. +It chooses the first n and last m phases from +""" +function crossover() +end diff --git a/julia/src/plotting.jl b/julia/src/plotting.jl index fea4c5c..18b5856 100644 --- a/julia/src/plotting.jl +++ b/julia/src/plotting.jl @@ -16,8 +16,8 @@ function get_true_max(mat::Vector{Array{Float64,2}}) return maximum(abs.(flatten(mat))) end -function plot_orbits(paths::Vector{Vector{Vector{Float64}}}; - primary::String="Earth", +function plot_orbits(paths::Vector{Vector{Vector{Float64}}}, + primary::Body=Sun; labels::Vector{String}=Vector{String}(), title::String="Spacecraft Position", colors::Vector{String}=Vector{String}()) @@ -28,7 +28,7 @@ function plot_orbits(paths::Vector{Vector{Vector{Float64}}}; x = cos.(θ) * sin.(ϕ)' y = sin.(θ) * sin.(ϕ)' z = repeat(cos.(ϕ)',outer=[N, 1]) - ps = rs[primary] .* (x,y,z) + ps = primary.r .* (x,y,z) x_p,y_p,z_p = ps t1 = [] @@ -52,7 +52,7 @@ function plot_orbits(paths::Vector{Vector{Vector{Float64}}}; y=(y_p), z=(z_p), showscale=false, - colorscale = p_colors[primary]) + colorscale = primary.color) layout = Layout(title=title, width=1000, diff --git a/julia/src/spacecraft.jl b/julia/src/spacecraft.jl index fbffbca..25620d3 100644 --- a/julia/src/spacecraft.jl +++ b/julia/src/spacecraft.jl @@ -1,4 +1,5 @@ -export Sc +export Sc, test_sc, bepi, no_thrust + mutable struct Sc dry_mass::Float64 mass_flow_rate::Float64 @@ -7,18 +8,6 @@ mutable struct Sc duty_cycle::Float64 end -function Sc(name::String) - # This has extra thrusters to make plots more visible (and most don't use fuel) - if name == "test" - return Sc(9000., 0.00025/(2000*0.00981), 0.00025, 50, 0.9) - # This is the normal one - elseif name == "bepi" - return Sc(9000., 2*0.00025/(2000*0.00981), 0.00025, 2, 0.9) - elseif name == "no_thrust" - return Sc(9000., 0.01, 0., 0, 0.) - else - throw(ErrorException("Bad sc name")) - end -end - -Base.copy(s::Sc) = Sc(s.dry_mass, s.mass_flow_rate, s.max_thrust, s.num_thrusters, s.duty_cycle) +const test_sc = Sc(8000., 0.00025/(2000*0.00981), 0.00025, 50, 0.9) +const bepi = Sc(8000., 2*0.00025/(2000*0.00981), 0.00025, 2, 0.9) +const no_thrust = Sc(8000., 0.01, 0., 0, 0.) diff --git a/julia/test/inner_loop/find_closest.jl b/julia/test/inner_loop/find_closest.jl deleted file mode 100644 index 199377e..0000000 --- a/julia/test/inner_loop/find_closest.jl +++ /dev/null @@ -1,50 +0,0 @@ -@testset "Find Closest" begin - - println("Testing NLP solver") - - using NLsolve, PlotlyJS - - # Initial Setup - sc = Sc("test") - fresh_sc = copy(sc) - a = rand(25000:1.:40000) - e = rand(0.01:0.01:0.05) - i = rand(0.01:0.01:π/6) - T = 2π*√(a^3/μs["Earth"]) - prop_time = 5T - n = 200 - - # A simple orbit raising - start_mass = 10_000. - start = [ oe_to_xyz([ a, e, i, 0., 0., 0. ], μs["Earth"]); start_mass ] - Tx, Ty, Tz = conv_T(repeat([0.9], n), repeat([0.], n), repeat([0.], n), - start, - sc, - prop_time, - μs["Earth"]) - final = prop(hcat(Tx, Ty, Tz), start, copy(sc), μs["Earth"], prop_time)[2] - new_T = 2π*√(xyz_to_oe(final, μs["Earth"])[1]^3/μs["Earth"]) - - # This should be close enough to converge - Tx, Ty, Tz = conv_T(repeat([0.89], n), repeat([0.], n), repeat([0.], n), - start, - 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 result.converged - path1 = prop(zeros((100,3)), start, sc, μs["Earth"], T)[1] - path2, calc_final = prop(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, fresh_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") - if result.converged - @test norm(calc_final[1:6] - final[1:6]) < 1e-4 - end - -end diff --git a/julia/test/inner_loop/laguerre-conway.jl b/julia/test/inner_loop/laguerre-conway.jl index 20a9b9d..7e1a03f 100644 --- a/julia/test/inner_loop/laguerre-conway.jl +++ b/julia/test/inner_loop/laguerre-conway.jl @@ -5,26 +5,29 @@ using Thesis: laguerre_conway # Test that the propagator produces good periodic orbits (forwards and backwards) - for T in rand(3600*1.5:3600*4, (5)) - start = oe_to_xyz([ (μs["Earth"]*(T/(2π))^2)^(1/3), rand(0.01:0.01:0.5), rand(0.01:0.01:0.45π), 0., 0., 1. ], μs["Earth"]) + μ = Earth.μ + for T in rand(3600*1.5:3600*4, 5) + e = rand(0.0:0.01:0.75) + i = rand(0.0:0.01:0.499π) + start = [ oe_to_xyz([ (μ*(T/(2π))^2)^(1/3), e, i, 0., 0., 1. ], μ); 12_000. ] orbit = start for _ in 1:5 i = 0. while i < T - orbit = laguerre_conway(orbit, μs["Earth"], 1.) + orbit = laguerre_conway(orbit, 1., Earth) i += 1 end @test i ≈ T - @test norm(orbit - start) < 1e-2 + @test norm(orbit[1:6] - start[1:6]) < 1e-2 end for _ in 1:5 i = 0. while i > -T - orbit = laguerre_conway(orbit, μs["Earth"], -1.) + orbit = laguerre_conway(orbit, -1., Earth) i -= 1 end @test i ≈ -T - @test norm(orbit - start) < 1e-2 + @test norm(orbit[1:6] - start[1:6]) < 1e-2 end end diff --git a/julia/test/inner_loop/phase.jl b/julia/test/inner_loop/phase.jl new file mode 100644 index 0000000..e422f9c --- /dev/null +++ b/julia/test/inner_loop/phase.jl @@ -0,0 +1,39 @@ +@testset "Phase" begin + + println("Testing NLP solver") + + using NLsolve, PlotlyJS + + # Initial Setup + T = rand( 2hour : second : 12hour) + revs = 5 + n = 10revs + start_mass = 12_000. + + # A simple orbit raising + start = gen_orbit(T, start_mass, Earth) + thrust = spiral(0.9, n, start, test_sc, revs*T, Earth) + final = prop(thrust, start, test_sc, revs*T, Earth)[2] + new_T = 2π * √( xyz_to_oe(final, Earth.μ)[1]^3 / Earth.μ ) + + # This should be close enough to converge + thrust_guess = spiral(0.88, n, start, test_sc, revs*T, Earth) + result = solve_phase(start, final, test_sc, revs*T, thrust_guess, Earth) + calc_path, calc_final = prop(result.zero, start, test_sc, revs*T, Earth) + + # Test + @test converged(result) + @test norm(calc_final[1:6] - final[1:6]) < 1e-5 + + # Plot + paths = Pathlist() + push!(paths, prop(start, T, Earth), + calc_path, + prop(calc_final, T, Earth), + prop(final, T, Earth) ) + fig = plot_orbits(paths, Earth, + labels=["init", "transit", "post-transit", "final"], + colors=["#FFF","#F44","#4F4","#44F"]) + savefig(fig, "../plots/nlp_test.html") + +end diff --git a/julia/test/inner_loop/propagator.jl b/julia/test/inner_loop/propagator.jl index 595ad68..68ef7c4 100644 --- a/julia/test/inner_loop/propagator.jl +++ b/julia/test/inner_loop/propagator.jl @@ -6,26 +6,19 @@ # Set up start_mass = 10_000. - start = [oe_to_xyz([ (μs["Earth"]*(rand(3600*1.5:0.01:3600*4)/(2π))^2)^(1/3), - rand(0.01:0.01:0.5), - rand(0.01:0.01:0.45π), - 0., - 0., - 1. ], μs["Earth"]); start_mass] - stepsize = rand(100.0:0.01:500.0) + start = gen_orbit(rand(.5year : hour : 2year), start_mass) + stepsize = rand(hour : second : 6hour) # Test that Laguerre-Conway is the default propagator for spacecrafts - craft = Sc("no_thrust") - state = prop_one([0., 0., 0.], start, craft, μs["Earth"], stepsize) - @test laguerre_conway(start, μs["Earth"], stepsize) ≈ state[1:6] + state = prop_one([0., 0., 0.], start, no_thrust, stepsize) + @test laguerre_conway(start, stepsize) ≈ state[1:6] @test state[7] == start_mass # Test that mass is reduced properly - craft = Sc("test") - state = prop_one([1., 0., 0.], start, craft, μs["Earth"], stepsize) - @test state[7] == start_mass - craft.mass_flow_rate*stepsize + state = prop_one([1., 0., 0.], start, bepi, stepsize) + @test state[7] == start_mass - bepi.mass_flow_rate*stepsize # Test that a bad ΔV throws an error - @test_throws ErrorException prop_one([1.5, 0., 0.], start, craft, μs["Earth"], stepsize) + @test_throws Thesis.PropOne_Error prop_one([1.5, 0., 0.], start, bepi, stepsize) end diff --git a/julia/test/plotting.jl b/julia/test/plotting.jl index 46b93f6..739a27e 100644 --- a/julia/test/plotting.jl +++ b/julia/test/plotting.jl @@ -4,21 +4,32 @@ using PlotlyJS - # First some setup - sc = Sc("test") - T = rand(3600*2:0.01:3600*4) - start = [oe_to_xyz([ (μs["Earth"]*(T/(2π))^2)^(1/3), - 0.1, - π/4, - 0., - 0., - 1. ], μs["Earth"]); 10_000.] + # Plot an earth plot + T = rand(2hour : 1 : 4hour) revs = 30 n = revs*100 - ΔVs = repeat([0.9, 0., 0.]', outer=(n,1)) - path = prop(ΔVs, start, copy(sc), μs["Earth"], revs*T)[1] - p = plot_orbits([path]) - savefig(p,"../plots/plot_test.html") + + start = gen_orbit(T, 12_000., Earth) + thrust = spiral(0.9, n, start, test_sc, revs*T, Earth) + path = prop(thrust, start, test_sc, revs*T, Earth)[1] + + p = plot_orbits([path], Earth) + savefig(p,"../plots/plot_test_earth.html") @test typeof(p) == PlotlyJS.SyncPlot + # Now change a little bit and plot around the Sun + # This also checks that the spacecraft are configured right: + # They really shouldn't run out of fuel in 4 years + T = rand(year : hour : 4year) + tof = 4year + start = gen_orbit(T, 12_000.) + thrust = spiral(0.9, n, start, bepi, tof) + sun_paths = Vector{Vector{Vector{Float64}}}() + push!(sun_paths, prop(zeros(100,3), start, bepi, tof)[1]) + push!(sun_paths, prop(thrust, start, bepi, tof)[1]) + p = plot_orbits(sun_paths) + savefig(p,"../plots/plot_test_sun.html") + @test typeof(p) == PlotlyJS.SyncPlot + + end diff --git a/julia/test/runtests.jl b/julia/test/runtests.jl index 11195ed..9dbd4d7 100644 --- a/julia/test/runtests.jl +++ b/julia/test/runtests.jl @@ -4,24 +4,13 @@ using LinearAlgebra using SPICE using Thesis -try - furnsh("../../SPICE/naif0012.tls") - furnsh("../../SPICE/de430.bsp") -catch - furnsh("SPICE/naif0012.tls") - furnsh("SPICE/de430.bsp") -end - # Tests @testset "All Tests" begin - include("spacecraft.jl") include("plotting.jl") include("inner_loop/laguerre-conway.jl") include("inner_loop/propagator.jl") - include("inner_loop/find_closest.jl") - include("inner_loop/monotonic_basin_hopping.jl") - include("inner_loop/inner_loop.jl") - include("outer_loop.jl") + include("inner_loop/phase.jl") + # include("inner_loop/monotonic_basin_hopping.jl") end print() diff --git a/julia/test/spacecraft.jl b/julia/test/spacecraft.jl deleted file mode 100644 index 274b764..0000000 --- a/julia/test/spacecraft.jl +++ /dev/null @@ -1,28 +0,0 @@ -@testset "Spacecraft Construction" begin - - println("Testing spacecraft") - - # Test that the standard spacecraft can be created - craft = Sc("test") - @test craft.dry_mass == 9000. - @test craft.mass_flow_rate == craft.max_thrust/(0.00981*2000) - @test craft.max_thrust == 0.00025 - @test craft.num_thrusters == 50 - @test craft.duty_cycle == 0.9 - - craft = Sc("no_thrust") - @test craft.dry_mass == 9000. - @test craft.mass_flow_rate == 0.01 - @test craft.max_thrust == 0. - @test craft.num_thrusters == 0 - @test craft.duty_cycle == 0. - - # Test that the standard spacecraft can be copied - new_craft = copy(craft) - @test new_craft.dry_mass == craft.dry_mass - @test new_craft.mass_flow_rate == craft.mass_flow_rate - @test new_craft.max_thrust == craft.max_thrust - @test new_craft.num_thrusters == craft.num_thrusters - @test new_craft.duty_cycle == craft.duty_cycle - -end diff --git a/notes.md b/notes.md index c1caa63..644f203 100644 --- a/notes.md +++ b/notes.md @@ -1,5 +1,30 @@ - - Prioritize the flowcharts - - Be sure to include the error checking into the flowcharts - - look up Jacob Englander's PhD thesis for further background (UIUC) - - can begin reaching out to committee member to get a sense of availability in November in the next couple weeks - - Upload outline to the google drive \ No newline at end of file +# Notes + +## Meeting with Bosanac (2021-09-15) + +- Prioritize the flowcharts +- Be sure to include the error checking into the flowcharts +- look up Jacob Englander's PhD thesis for further background (UIUC) +- can begin reaching out to committee member to get a sense of availability in November in the next couple weeks +- Upload outline to the google drive + +## Refactor Notes (2021-09-25) + +I need to determine the list of decision variables for each step: + +- Defined by mission + - Launch Window + - Max launch C3 + - Max incoming C3 (or ||v∞||) + - Latest arrival +- Decided by outer loop + - Flyby choices +- Decided by inner loop + - Launch date + - Launch v∞_out + - for each phase: + - tof + - v∞_in + - turning angle +- Decided by NLP + - thrust profile for each phase diff --git a/SPICE/de430.bsp b/spice_files/de430.bsp similarity index 100% rename from SPICE/de430.bsp rename to spice_files/de430.bsp diff --git a/SPICE/naif0012.tls b/spice_files/naif0012.tls similarity index 100% rename from SPICE/naif0012.tls rename to spice_files/naif0012.tls