diff --git a/julia/src/Thesis.jl b/julia/src/Thesis.jl index 9f863d8..b9cdedf 100644 --- a/julia/src/Thesis.jl +++ b/julia/src/Thesis.jl @@ -1,29 +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 using ForwardDiff using PlotlyJS using Distributed + 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 include("./errors.jl") include("./constants.jl") include("./spacecraft.jl") + include("./mission.jl") include("./conversions.jl") + include("./lamberts.jl") include("./plotting.jl") include("./inner_loop/laguerre-conway.jl") include("./inner_loop/propagator.jl") include("./inner_loop/phase.jl") - # include("./inner_loop/monotonic_basin_hopping.jl") + include("./inner_loop/monotonic_basin_hopping.jl") # include("./outer_loop.jl") end diff --git a/julia/src/inner_loop/monotonic_basin_hopping.jl b/julia/src/inner_loop/monotonic_basin_hopping.jl index e4d3cf6..7813886 100644 --- a/julia/src/inner_loop/monotonic_basin_hopping.jl +++ b/julia/src/inner_loop/monotonic_basin_hopping.jl @@ -1,86 +1,243 @@ +using Dates + export mbh """ -Generates n pareto-distributed random numbers +Generates pareto-distributed random numbers + +This is usually around one to two percent, but sometimes up to ten or so (fat tails) """ -function pareto(α::Float64, n::Int) - s = rand((-1,1), (n,3)) - r = rand(Float64, (n,3)) - ϵ = 1e-10 - return (s./ϵ) * (α - 1.0) ./ (ϵ ./ (ϵ .+ r)).^(-α) +function pareto(shape::Tuple{Int, Int}, α::Float64=1.01) + s = rand((-1,1), shape) + r = rand(Float64, shape) + ϵ = 1e-10 + return 1 .+ (s./ϵ) * (α - 1.0) ./ (ϵ ./ (ϵ .+ r)).^(-α) +end + +function pareto(shape::Int, α::Float64=1.01) + s = rand((-1,1), shape) + r = rand(Float64, shape) + ϵ = 1e-10 + return 1 .+ (s./ϵ) * (α - 1.0) ./ (ϵ ./ (ϵ .+ r)).^(-α) +end + +function pareto(α::Float64=1.01) + s = rand((-1,1)) + r = rand(Float64) + ϵ = 1e-10 + return 1 + (s./ϵ) * (α - 1.0) ./ (ϵ ./ (ϵ .+ r)).^(-α) +end + +function pareto_add(α::Float64=1.01) + s = rand((-1,1)) + r = rand(Float64) + ϵ = 1e-10 + return (s./ϵ) * (α - 1.0) ./ (ϵ ./ (ϵ .+ r)).^(-α) end """ -Perturbs the monotonic basin hopping decision vector -TODO: This needs to be updated +Returns a random date between two dates """ -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 +function gen_date(date_range::Tuple{DateTime, DateTime}) + l0, lf = date_range + l0 + Dates.Millisecond(floor(rand()*(lf-l0).value)) end -function new_x(n::Int) - 2.0 * rand(Float64, (n,3)) .- 1. +""" +Returns a random amount of time in a range +""" +function gen_period(date_range::Tuple{DateTime, DateTime}) + l0, lf = date_range + Dates.Millisecond(floor(rand()*(lf-l0).value)) end -function mbh(flybys::Vector{Planet}) -end - -function mbh(start::AbstractVector, - final::AbstractVector, - craft::Sc, - μ::AbstractFloat, - t0::AbstractFloat, - tf::AbstractFloat, - n::Int; - search_patience_lim::Int=2000, - drill_patience_lim::Int=40, - tol=1e-6, - verbose=false) - - archive = [] - - i = 0 - x_current = Result(false, 1e8*ones(n,3)) - while i < search_patience_lim - i += 1 - drill_impatience = 0 - 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(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 +""" +Perturbs a valid mission with pareto-distributed variables, generating a mission guess +""" +function perturb(mission::Mission) + new_launch_date = mission.launch_date + Dates.Second(floor(7day * pareto_add())) + new_launch_v∞ = mission.launch_v∞ .* pareto(3) + new_phases = Vector{Phase}() + for phase in mission.phases + new_v∞_in = phase.v∞_in .* pareto(3) + new_δ = phase.δ * pareto() + new_tof = phase.tof * pareto() + new_thrust_profile = phase.thrust_profile .* pareto(size(phase.thrust_profile)) + push!(new_phases, Phase(phase.planet, new_v∞_in, new_δ, new_tof, new_thrust_profile)) end - if verbose println() end + # TODO: Mission_Guess.validate() + Mission_Guess(mission.sc, mission.start_mass, new_launch_date, new_launch_v∞, new_phases) +end - if archive == [] error("there were no converged paths") end +""" +Generates a new randomized guess for the mission decision variables +""" +function mission_guess( flybys::Vector{Body}, + sc::Sc, + start_mass::Float64, + launch_window::Tuple{DateTime, DateTime}, + max_C3_out::Float64, + max_v∞_in_mag::Float64, + latest_arrival::DateTime, + primary::Body=Sun ) + # TODO: Eventually I can calculate n more intelligently + n = 20 - current_best_mass = 1e8 - best = archive[1] - for candidate in archive - if mass_est(candidate.zero) < current_best_mass - current_best_mass = mass_est(candidate.zero) - best = candidate + # Determine the launch conditions + launch_date = gen_date(launch_window) + launch_v∞_normalized = rand(-1:0.0001:1, 3) + launch_v∞ = rand(0:0.0001:√max_C3_out) * launch_v∞_normalized/norm(launch_v∞_normalized) + + # Determine the leg lengths + num_phases = length(flybys) - 1 + total_tof = 100year + max_tof = (latest_arrival - launch_date).value / 1000 + tofs = rand(30day : hour : 0.7max_tof, num_phases) + total_tof = sum(tofs) + while total_tof > max_tof + tofs = rand(30day : hour : 0.7max_tof, num_phases) + total_tof = sum(tofs) + end + + # For each phase, determine the v∞_in and δ + phases = Vector{Phase}() + for i in 1:num_phases + flyby = flybys[i+1] + v∞_normalized = rand(-1:0.0001:1, 3) + v∞ = rand(0:0.0001:10) * v∞_normalized/norm(v∞_normalized) + δ = rand(0:0.0001:2π) + periapsis = (flyby.μ/(v∞ ⋅ v∞)) * ( 1/sin(δ/2) - 1 ) + while periapsis < flyby.r + 100. + δ = rand(0:0.0001:2π) + periapsis = (flyby.μ/(v∞ ⋅ v∞)) * ( 1/sin(δ/2) - 1 ) + end + thrusts = rand(-1:0.0001:1,(n,3)) + push!(phases, Phase(flyby, v∞, δ, tofs[i], thrusts)) + end + + # Finally, determine the arrival v∞ + arrival_v∞_normalized = rand(-1:0.0001:1, 3) + arrival_v∞ = rand(0:0.0001:max_v∞_in_mag) * arrival_v∞_normalized/norm(arrival_v∞_normalized) + + # And we can construct a mission guess object with these values + Mission_Guess( sc, start_mass, launch_date, launch_v∞, phases ) +end + +""" +A convenience function for calculating mass usage given a certain thrust profile +""" +function mass_consumption(sc::Sc, phase::Phase) + weighted_thrusting_time = 0.0 + n = size(phase.thrust_profile)[1] + for i in 1:size(phase.thrust_profile,1) + weighted_thrusting_time += norm(phase.thrust_profile[i,:]) * phase.tof/n + end + return weighted_thrusting_time*sc.mass_flow_rate +end + +""" +This attempts to determine v∞_out from v∞_in and the turning angle, assuming we are staying in the +plane of the three planets +""" +function calc_turn(p1::Body, p2::Body, p3::Body, v∞_in::Vector{Float64}, δ::Float64) +end + +""" +Sequentially calls the NLP solver to attempt to solve based on the initial guess +""" +function inner_loop_solve(guess::Mission_Guess) + v∞_out = guess.launch_v∞ + time = utc2et(Dates.format(guess.launch_date,"yyyy-mm-ddTHH:MM:SS")) + current_planet = Earth + start = [spkssb(Earth.id, time, "ECLIPJ2000"); 0.0] + [ zeros(3); guess.launch_v∞; guess.start_mass ] + + corrected_phases = Vector{Phase}() + for i in 1:length(guess.phases) + phase = guess.phases[i] + time += phase.tof + goal = spkssb(phase.planet.id, time, "ECLIPJ2000") + [zeros(3); phase.v∞_in] + result = solve_phase( start, goal, guess.sc, phase.tof, phase.thrust_profile) + converged(result) || return Bad_Mission() # Drop if it's not working + corrected_phase = Phase(phase.planet, phase.v∞_in, phase.δ, phase.tof, result.zero) + push!(corrected_phases, corrected_phase) + mass_used = mass_consumption(guess.sc, corrected_phase) + if i != length(guess.phases) + v∞_out = calc_turn(current_planet, phase.planet, guess.phases[i+1].planet, phase.v∞_in, phase.δ) + current_planet = phase.planet + planet_state = [spkssb(current_planet.id, time, "ECLIPJ2000"); 0.0] + start = planet_state + [ zeros(3); v∞_out; start_mass - mass_used ] end end - return best, archive - + return Mission(guess.sc, guess.start_mass, guess.launch_date, guess.launch_v∞, corrected_phases) end + +""" +The cost function for the mission +TODO: This will probably move and eventually be passed as an argument +""" +function cost(mission::Mission, max_C3::Float64, max_v∞::Float64) + mass_used = 0.0 + for phase in mission.phases mass_used += mass_used(sc, phase) end + mass_percent = mass_used/mission.sc.dry_mass + C3_percent = ( mission.launch_v∞ ⋅ mission.launch_v∞ ) / max_C3 + v∞_percent = norm(mission.phases[end].v∞_in) / max_v∞ + return 3mass_percent + C3_percent + v∞_percent +end + + + +function mbh( flybys::Vector{Body}, + sc::Sc, + start_mass::Float64, + launch_window::Pair{Date}, + max_C3::Float64, + max_v∞::Float64, + latest_arrival::Date, + primary::Body=Sun; + search_patience::Int=1_000, + drill_patience::Int=50) + # Let's pseudo-code this bitch + # + # First, we need a function (mission_guess below) that randomly generates the: + # - Launch date + # - Launch v∞_out + # - for each phase: + # - tof + # - v∞_in + # - turning angle + # + # Also need a function (inner_loop_solve below) that takes the generated decision vector guess and + # attempts to correct it with the NLP solver. It will either converge or not. + # + # Also need a costing function (may be provided potentially) + # + # Also need a perturb function + # + guess = mission_guess(flybys, sc, start_mass, launch_window, max_C3, max_v∞, latest_arrival, primary) + inner_loop_solve(guess) + # search_count = 0 + # x_current = "Bad" + # while search_count < search_patience + # search_count += 1 + # drill_count = 0 + # x_star = inner_loop_solve(mission_guess()) + # if x_star.converged + # if cost(x_star) < cost(x_current) + # x_current = x_star + # while drill_count < drill_patience + # x_star = inner_loop_solve(perturb(x_current)) + # if x_star.converged and cost(x_star) < cost(x_current) + # x_current = x_star + # drill_count = 0 + # else + # drill_count += 1 + # end + # end + # push!(archive, x_current) + # end + # end + # + # +end + diff --git a/julia/src/inner_loop/phase.jl b/julia/src/inner_loop/phase.jl index 656c34d..5044a89 100644 --- a/julia/src/inner_loop/phase.jl +++ b/julia/src/inner_loop/phase.jl @@ -31,8 +31,7 @@ function solve_phase( start::Vector{Float64}, end result = nlsolve(f!, atanh.(x0), ftol=tol, autodiff=:forward, iterations=num_iters) - converged(result) || throw(Convergence_Error()) - result.zero = tanh.(result.zero) + if converged(result) result.zero = tanh.(result.zero) end return result diff --git a/julia/src/inner_loop/propagator.jl b/julia/src/inner_loop/propagator.jl index b0b2d9e..1acc502 100644 --- a/julia/src/inner_loop/propagator.jl +++ b/julia/src/inner_loop/propagator.jl @@ -68,4 +68,4 @@ 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] +prop(x::Vector{Float64}, t::Float64, p::Body=Sun) = prop(zeros(1000,3), [x;1e10], no_thrust, t, p)[1] diff --git a/julia/src/lamberts.jl b/julia/src/lamberts.jl new file mode 100644 index 0000000..1ae7bfd --- /dev/null +++ b/julia/src/lamberts.jl @@ -0,0 +1,72 @@ +using Dates + +""" +I didn't want to have to write this...I'm going to try to do this as quickly as possible from old, +bad code + +Convenience function for solving lambert's problem +""" +function lamberts(planet1::Body,planet2::Body,leave::DateTime,arrive::DateTime) + + time_leave = utc2et(Dates.format(leave,"yyyy-mm-ddTHH:MM:SS")) + time_arrive = utc2et(Dates.format(arrive,"yyyy-mm-ddTHH:MM:SS")) + tof_req = time_arrive - time_leave + + state1 = [spkssb(planet1.id, time_leave, "ECLIPJ2000"); 0.0] + state2 = [spkssb(planet2.id, time_arrive, "ECLIPJ2000"); 0.0] + + r1 = state1[1:3] ; r1mag = norm(r1) + r2 = state2[1:3] ; r2mag = norm(r2) + μ = Sun.μ + + cos_dθ = dot(r1,r2)/(r1mag*r2mag) + dθ = atan(r2[2],r2[1]) - atan(r1[2],r1[1]) + dθ = dθ > 2π ? dθ-2π : dθ + dθ = dθ < 0.0 ? dθ+2π : dθ + DM = abs(dθ) > π ? -1 : 1 + A = DM * √(r1mag * r2mag * (1 + cos_dθ)) + dθ == 0 || A == 0 && error("Can't solve Lambert's Problem") + + ψ, c2, c3 = 0, 1//2, 1//6 + ψ_down = -4π ; ψ_up = 4π^2 + y = r1mag + r2mag + (A*(ψ*c3 - 1)) / √(c2) ; χ = √(y/c2) + tof = ( χ^3*c3 + A*√(y) ) / √(μ) + + i = 0 + while abs(tof-tof_req) > 1e-2 + y = r1mag + r2mag + (A*(ψ*c3 - 1)) / √(c2) + while y/c2 <= 0 + # println("You finally hit that weird issue... ") + ψ += 0.1 + if ψ > 1e-6 + c2 = (1 - cos(√(ψ))) / ψ ; c3 = (√(ψ) - sin(√(ψ))) / √(ψ^3) + elseif ψ < -1e-6 + c2 = (1 - cosh(√(-ψ))) / ψ ; c3 = (-√(-ψ) + sinh(√(-ψ))) / √((-ψ)^3) + else + c2 = 1//2 ; c3 = 1//6 + end + y = r1mag + r2mag + (A*(ψ*c3 - 1)) / √(c2) + end + χ = √(y/c2) + + tof = ( c3*χ^3 + A*√(y) ) / √(μ) + tof < tof_req ? ψ_down = ψ : ψ_up = ψ + ψ = (ψ_up + ψ_down) / 2 + + if ψ > 1e-6 + c2 = (1 - cos(√(ψ))) / ψ ; c3 = (√(ψ) - sin(√(ψ))) / √(ψ^3) + elseif ψ < -1e-6 + c2 = (1 - cosh(√(-ψ))) / ψ ; c3 = (-√(-ψ) + sinh(√(-ψ))) / √((-ψ)^3) + else + c2 = 1//2 ; c3 = 1//6 + end + + i += 1 + i > 500 && return [NaN,NaN,NaN],[NaN,NaN,NaN] + end + + f = 1 - y/r1mag ; g_dot = 1 - y/r2mag ; g = A * √(y/μ) + v0t = (r2 - f*r1)/g ; vft = (g_dot*r2 - r1)/g + return v0t - state1[4:6], vft - state2[4:6], tof_req + +end diff --git a/julia/src/mission.jl b/julia/src/mission.jl new file mode 100644 index 0000000..6469f60 --- /dev/null +++ b/julia/src/mission.jl @@ -0,0 +1,61 @@ +using Dates + +export Phase, Mission_Guess, Mission, Bad_Mission +export test_phase1, test_phase2 +export test_mission_guess, test_mission_guess_simple +export test_mission, test_mission_simple + +struct Phase + planet::Body + v∞_in::Vector{Float64} + δ::Float64 + tof::Float64 + thrust_profile::Matrix{Float64} +end +const test_phase1 = Phase(Venus, [10.4321, -6.3015, -0.01978], 0.2, 1.30464e7, zeros(20,3)) +const test_phase2 = Phase(Jupiter, [0.3, 7.1, 0.2], 2π, 3.9year, zeros(20,3)) + +struct Mission_Guess + sc::Sc + start_mass::Float64 + launch_date::DateTime + launch_v∞::Vector{Float64} + phases::Vector{Phase} + converged::Bool +end +Mission_Guess(args...) = Mission_Guess(args..., false) +const test_mission_guess = Mission_Guess( bepi, + 12_000., + DateTime(1992, 11, 19), + [-3.4, 1.2, 0.1], + [test_phase1, test_phase2] ) +const test_mission_guess_simple = Mission_Guess(bepi, + 12_000., + DateTime(1992, 11, 19), + [-3.4, 1.2, 0.1], + [test_phase1]) + +struct Mission + sc::Sc + start_mass::Float64 + launch_date::DateTime + launch_v∞::Vector{Float64} + phases::Vector{Phase} + converged::Bool +end +Mission(args...) = Mission(args..., true) +const test_mission = Mission(bepi, + 12_000., + DateTime(1992, 11, 19), + [-3.4, 1.2, 0.1], + [test_phase1, test_phase2]) +const test_mission_simple = Mission(bepi, + 12_000., + DateTime(1992, 11, 19), + [4.2984, -4.3272668, 1.43752], + [test_phase1]) + +struct Bad_Mission + converged::Bool +end +Bad_Mission() = Bad_Mission(false) diff --git a/julia/test/inner_loop/monotonic_basin_hopping.jl b/julia/test/inner_loop/monotonic_basin_hopping.jl index 455efa4..70f25e3 100644 --- a/julia/test/inner_loop/monotonic_basin_hopping.jl +++ b/julia/test/inner_loop/monotonic_basin_hopping.jl @@ -4,111 +4,42 @@ println("Testing Monotonic Basin Hopper") - # Initial Setup - # sc = Sc("test") - # a = rand(50_000:1.:100_000) - # e = rand(0.01:0.01:0.5) - # i = rand(0.01:0.01:π/6) - # T = 2π*√(a^3/μs["Earth"]) - # prop_time = 0.5T - # n = 20 - # start_mass = 10_000. - - # # A simple orbit raising - # start = [oe_to_xyz([ a, e, i, 0., 0., 0. ], μs["Earth"]); start_mass] - # Tx, Ty, Tz = conv_T(repeat([0.8], n), repeat([0.], n), repeat([0.], n), - # start, - # sc, - # prop_time, - # μs["Earth"]) - # nominal_path, 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"]) - - # # Find the best solution - # best, archive = mbh(start, - # final, - # sc, - # μs["Earth"], - # 0.0, - # prop_time, - # n, - # search_patience_lim=25, - # drill_patience_lim=50, - # verbose=true) - - # # Test and plot - # @test best.converged - # transit, calc_final = prop(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 = calc_final[end] - # nominal_mass = final[end] - # masses = [] - # for candidate in archive - # @test candidate.converged - # path2, calc_final = prop(candidate.zero, start, sc, μs["Earth"], prop_time) - # push!(masses, calc_final[end]) - # @test norm(calc_final[1:6] - final[1:6]) < 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 (start_mass - best_mass) < 1.1 * (start_mass - nominal_mass) - - # Now let's test a sun case. This should be pretty close to begin with - start_mass = 10_000. - launch_date = DateTime(2016,3,28) - launch_j2000 = utc2et(Dates.format(launch_date,"yyyy-mm-ddTHH:MM:SS")) - earth_start = [spkssb(ids["Earth"], launch_j2000, "ECLIPJ2000"); start_mass] - earth_speed = earth_start[4:6] - v∞ = 3.0*earth_speed/norm(earth_speed) - start = earth_start + [zeros(3); v∞; 0.0] - tof = 3600*24*30*10.75 - mars_state = [spkssb(Thesis.ids["Mars"], launch_j2000+tof, "ECLIPJ2000"); start_mass] - final = mars_state + [ zeros(3); [-1.1, -3., -2.6]; 0.0 ] - a = xyz_to_oe(final, μs["Sun"])[1] - T = 2π*√(a^3/μs["Sun"]) - n = 20 + # First we test the random mission guess generator + flybys = [Earth, Venus, Jupiter] + launch_window = ( DateTime(2021,12,25), DateTime(2025,12,25) ) + max_C3 = 10. + max_v∞ = 8. + latest_arrival = DateTime(2035,12,25) + random_guess = Thesis.mission_guess(flybys, bepi, 12_000., launch_window, max_C3, max_v∞, latest_arrival) + @test typeof(random_guess) == Mission_Guess - # But we'll plot to see - beginning_path = prop(zeros(100,3), start, Sc("test"), μs["Sun"], tof)[1] - ending_path = prop(zeros(100,3), final, Sc("test"), μs["Sun"], T)[1] - savefig(plot_orbits([beginning_path, ending_path], - labels=["initial", "ending"], - colors=["#F2F", "#2F2"]), - "../plots/mbh_sun_initial.html") - - # Now we solve and plot the new case - best, archive = mbh(start, - final, - Sc("test"), - μs["Sun"], - 0.0, - tof, - n, - search_patience_lim=25, - drill_patience_lim=50, - verbose=true) - solved_path, solved_state = prop(best.zero, start, Sc("test"), μs["Sun"], tof) - ending_path = prop(zeros(100,3), final, Sc("test"), μs["Sun"], T)[1] - savefig(plot_orbits([solved_path, ending_path], - labels=["best", "ending"], - colors=["#C2F", "#2F2"]), - "../plots/mbh_sun_solved.html") - - # We'll just make sure that this at least converged correctly - @test norm(solved_state[1:6] - final[1:6]) < 1e-4 + # Then the perturb function + mission_guess = Thesis.perturb(test_mission) + @test mission_guess.launch_date != test_mission.launch_date + @test mission_guess.launch_v∞ != test_mission.launch_v∞ + for i in 1:2 + @test mission_guess.phases[i].v∞_in != test_mission.phases[i].v∞_in + @test mission_guess.phases[i].δ != test_mission.phases[i].δ + @test mission_guess.phases[i].tof != test_mission.phases[i].tof + end + @test !mission_guess.converged + + # # Then the inner loop builder function + mission = Thesis.inner_loop_solve(test_mission_guess) + @test !mission.converged + # For the valid case we need to use a lambert's solver + # TODO: This is probably not acceptable for how close I have to be + leave = DateTime(1992,11,19) + arrive = DateTime(1993,4,1) + time_leave = utc2et(Dates.format(leave,"yyyy-mm-ddTHH:MM:SS")) + time_arrive = utc2et(Dates.format(arrive,"yyyy-mm-ddTHH:MM:SS")) + earth_state = [spkssb(Earth.id, time_leave, "ECLIPJ2000"); 0.0] + venus_state = [spkssb(Venus.id, time_arrive, "ECLIPJ2000"); 0.0] + v∞_out, v∞_in, tof = Thesis.lamberts(Earth, Venus, leave, arrive) + phase = Phase(Venus, 1.0001v∞_in, 0.2, tof, 0.0*ones(20,3)) + mission_guess = Mission_Guess(bepi, 12_000., leave, v∞_out, [phase]) + mission = Thesis.inner_loop_solve(mission_guess) + @test mission.converged end diff --git a/julia/test/runtests.jl b/julia/test/runtests.jl index 7dcbc8d..97b7a2c 100644 --- a/julia/test/runtests.jl +++ b/julia/test/runtests.jl @@ -4,10 +4,18 @@ using LinearAlgebra using SPICE using Thesis -@testset "All Tests" begin - include("plotting.jl") - include("inner_loop/laguerre-conway.jl") - include("inner_loop/propagator.jl") - include("inner_loop/phase.jl") - # include("inner_loop/monotonic_basin_hopping.jl") +try + furnsh("../../spice_files/naif0012.tls") + furnsh("../../spice_files/de430.bsp") +catch + furnsh("spice_files/naif0012.tls") + furnsh("spice_files/de430.bsp") +end + +@testset "All Tests" begin + # include("plotting.jl") + # include("inner_loop/laguerre-conway.jl") + # include("inner_loop/propagator.jl") + # include("inner_loop/phase.jl") + include("inner_loop/monotonic_basin_hopping.jl") end