diff --git a/LaTeX/flowcharts/mbh.png b/LaTeX/flowcharts/mbh.png
index 2e272c3..40dbadb 100644
Binary files a/LaTeX/flowcharts/mbh.png and b/LaTeX/flowcharts/mbh.png differ
diff --git a/LaTeX/flowcharts/mbh.svg b/LaTeX/flowcharts/mbh.svg
index 7fe629f..42c4390 100644
--- a/LaTeX/flowcharts/mbh.svg
+++ b/LaTeX/flowcharts/mbh.svg
@@ -7,7 +7,7 @@
viewBox="0 0 279.4 215.9"
version="1.1"
id="svg5"
- inkscape:version="1.1.1 (3bf5ae0d25, 2021-09-20)"
+ inkscape:version="1.1.1 (3bf5ae0d25, 2021-09-20, custom)"
sodipodi:docname="mbh.svg"
inkscape:export-filename="/home/connor/projects/thesis/LaTeX/flowcharts/mbh.png"
inkscape:export-xdpi="150"
@@ -27,15 +27,18 @@
inkscape:document-units="mm"
showgrid="true"
units="in"
- inkscape:zoom="1.1148803"
- inkscape:cx="548.48938"
- inkscape:cy="287.02633"
- inkscape:window-width="1916"
- inkscape:window-height="1045"
+ inkscape:zoom="0.78833942"
+ inkscape:cx="525.15451"
+ inkscape:cy="378.64401"
+ inkscape:window-width="1912"
+ inkscape:window-height="1040"
inkscape:window-x="1920"
- inkscape:window-y="31"
+ inkscape:window-y="32"
inkscape:window-maximized="1"
- inkscape:current-layer="layer1">
+ inkscape:current-layer="layer1"
+ inkscape:snap-smooth-nodes="false"
+ inkscape:snap-midpoints="false"
+ inkscape:snap-object-midpoints="false">
+
+
+
+
+
Generate Random Generate Random Mission Guess
+ id="tspan10731">Mission Guess
Correct the Guess
+ y="45.095099"
+ id="tspan10733">Correct the Guess
+
+ Correct the Guess
Correct the Guess
+ id="tspan10737">Correct the Guess
Store as Current Store as Current Best, Set Drill=0
+ id="tspan10741">Best, Set Drill=0
Perturb the Best Perturb the Best Guess (in Basin)
+ id="tspan10745">Guess (in Basin)
Hit SearchHit Search Limit?
+ id="tspan10751">Limit?
Hit DrillHit Drill Limit?
+ id="tspan10757">Limit?
BetterBetter thanthan currentcurrent best?
+ id="tspan10771">best?
BetterBetter thanthan currentcurrent best?
+ id="tspan10785">best?
Output Output Best Best Mission
+ id="tspan10791">Mission
+ inkscape:export-ydpi="150"
+ sodipodi:nodetypes="cc" />
+
+ inkscape:export-ydpi="150"
+ sodipodi:nodetypes="cc" />
+ inkscape:export-ydpi="150"
+ sodipodi:nodetypes="cc" />
True
+ id="tspan10793">True
False
+ id="tspan10795">False
False
+ id="tspan10797">False
False
+ id="tspan10799">False
False
+ id="tspan10801">False
False
+ id="tspan10803">False
False
+ id="tspan10805">False
True
+ id="tspan10807">True
True
+ id="tspan10809">True
True
+ id="tspan10811">True
True
+ id="tspan10813">True
True
+ id="tspan10815">True
Search += 1
+ id="tspan10817">Search += 1
+ Lambert's
+ Standard
Drill += 1
+ id="tspan10823">Drill += 1
Drill += 1
+ id="tspan10825">Drill += 1
+
+
+
+
+
diff --git a/LaTeX/thesis.tex b/LaTeX/thesis.tex
index 87dc440..297a44f 100644
--- a/LaTeX/thesis.tex
+++ b/LaTeX/thesis.tex
@@ -12,8 +12,8 @@
\degree{Master of Science}{M.S., Aerospace Engineering}
\dept{Department of}{Aerospace Engineering}
\advisor{Prof.}{Natasha Bosanac}
-\reader{TBD: Kathryn Davis}
-\readerThree{TBD: Daniel Scheeres}
+\reader{Kathryn Davis}
+\readerThree{Daniel Scheeres}
\abstract{ \OnePageChapter
There are a variety of approaches to finding and optimizing low-thrust trajectories in
@@ -573,7 +573,7 @@
Englander in his Hybrid Optimal Control Problem paper, but does not allow for missions
with multiple targets, simplifying the optimization.
- \section{Inner Loop Implementation}
+ \section{Inner Loop Implementation}\label{inner_loop_section}
The optimization routine can be reasonable separated into two separate ``loops'' wherein
the first loop is used, given an initial guess, to find valid trajectories within the
@@ -761,13 +761,117 @@
in order to achieve maximum precision.
\section{Outer Loop Implementation}
- Overview the outer loop. This may require a final flowchart, but might potentially be too
- simple to lend itself to one.
- \subsection{Inner Loop Calling Function}
- The primary reason for including this section is to discuss the error handling.
+ Now we have the tools in place for, given a potential ''mission guess`` in the
+ vicinity of a valid guess, attempting to find a valid and optimal solution in that
+ vicinity. Now what remains is to develop a routine for efficiently generating these
+ random mission guesses in such a way that thoroughly searches the entirety of the
+ solution space with enough granularity that all spaces are considered by the inner loop
+ solver.
- \subsection{Monotonic Basin Hopping}
+ Once that has been accomplished, all that remains is an ''outer loop`` that can generate
+ new guesses or perturb existing valid missions as needed in order to drill down into a
+ specific local minimum. In this thesis, that is accomplished through the use of a
+ Monotonic Basin Hopping algorithm. This will be described more thoroughly in
+ Section~\ref{mbh_subsection}, but Figure~\ref{mbh_flow} outlines the process steps of
+ the algorithm.
+
+ \begin{figure}
+ \centering
+ \includegraphics[width=\textwidth]{flowcharts/mbh}
+ \caption{A flowchart visualizing the steps in the monotonic basin hopping
+ algorithm}
+ \label{mbh_flow}
+ \end{figure}
+
+ \subsection{Random Mission Generation}
+
+ At a basic level, the algorithm needs to produce a mission guess (represented by all
+ of the values described in Section~\ref{inner_loop_section}) that contains random
+ values within reasonable bounds in the space. This leaves a number of variables open
+ to for implementation. For instance, it remains to be determined which distribution
+ function to use for the random values over each of those variables, which bounds to
+ use, as well as the possibilities for any improvements to a purely random search.
+
+ Currently, the first value set for the mission guess is that of $n$, which is the
+ number of sub-trajectories that each arc will be broken into for the Sims-Flanagan
+ based propagator. For this implementation, that was chosen to be 20, based upon a
+ number of tests in which the calculation time for the propagation was compared
+ against the accuracy of a much higher $n$ value for some known thrust controls, such
+ as a simple spiral orbit trajectory. This value of 20 tends to perform well and
+ provide reasonable accuracy, without producing too many variables for the NLP
+ optimizer to control for (since the impulsive thrust at the center of each of the
+ sub-trajectories is a control variable). This leaves some room for future
+ improvements, as will be discussed in Section~\ref{improvement_section}.
+
+ The bounds for the launch date are provided by the user in the form of a launch
+ window, so the initial launch date is just chosen as a standard random value from a
+ uniform distribution within those bounds.
+
+ A unit launch direction is then also chosen as a 3-length vector of uniform random
+ numbers, then normalized. This unit vector is then multiplied by a uniform random
+ number between 0 and the square root of the maximum launch $C_3$ specified by the
+ user to generate an initial $\vec{v_\infty}$ vector at launch.
+
+ Next, the times of flight of each phase of the mission is then decided. Since launch
+ date has already been selected, the maximum time of flight can be calculated by
+ subtracting the launch date from the latest arrival date provided by the mission
+ designer. Then, each leg is chosen from a uniform distribution with bounds currently
+ set to a minimum flight time of 30 days and a maximum of 70\% of the maximum time of
+ flight. These leg flight times are then iteratively re-generated until the total
+ time of flight (represented by the sum of the leg flight times) is less than the
+ maximum time of flight. This allows for a lot of flexibility in the leg flight
+ times, but does tend toward more often producing longer missions, particularly for
+ missions with more flybys.
+
+ Then, the internal components for each phase are generated. It is at this step, that
+ the mission guess generator splits the outputs into two separate outputs. The first
+ is meant to be truly random, as is generally used as input for a monotonic basin
+ hopping algorithm. The second utilizes a Lambert's solver to determine the
+ appropriate hyperbolic velocities (both in and out) at each flyby to generate a
+ natural trajectory arc. For this Lambert's case, the mission guess is simply seeded
+ with zero thrust controls and outputted to the monotonic basin hopper. The intention
+ here is that if the time of flights are randomly chosen so as to produce a
+ trajectory that is possible with a control in the vicinity of a natural trajectory,
+ we want to be sure to find that trajectory. More detail on how this is handled is
+ available in Section~\ref{mbh_subsection}.
+
+ However, for the truly random mission guess, there are still the $v_\infty$ values
+ and the initial thrust guesses to generate. For each of the phases, the incoming
+ excess hyperbolic velocity is calculated in much the same way that the launch
+ velocity was calculated. However, instead of multiplying the randomly generate unit
+ direction vector by a random number between 0 and the square root of the maximum
+ launch $C_3$, bounds of 0 and 10 kilometers per second are used instead, to provide
+ realistic flyby values.
+
+ The outgoing excess hyperbolic velocity at infinity is then calculated by again
+ choosing a uniform random unit direction vector, then by multiplying this value by
+ the magnitude of the incoming $v_{\infty}$ since this is a constraint of a
+ non-powered flyby.
+
+ From these two velocity vectors the turning angle, and thus the periapsis of the
+ flyby, can then be calculated by the following equations:
+
+ \begin{align}
+ \delta &= \arccos \left( \frac{\vec{v}_{\infty,in} \cdot
+ \vec{v}_{\infty,out}}{|v_{\infty,in}| \cdot {|v_{\infty,out}}|} \right) \\
+ r_p &= \frac{\mu}{\vec{v}_{\infty,in} \cdot \vec{v}_{\infty,out}} \cdot \left(
+ \frac{1}{\sin(\delta/2)} - 1 \right)
+ \end{align}
+
+ If this radius of periapse is then found to be less than the minimum safe radius
+ (currently set to the radius of the planet plus 100 kilometers), then the process is
+ repeated with new random flyby velocities until a valid seed flyby is found. These
+ checks are also performed each time a mission is perturbed or generated by the nlp
+ solver.
+
+ The final requirement then, is the thrust controls, which are actually quite simple.
+ Since the thrust is defined as a 3-vector of values between -1 and 1 representing
+ some percentage of the full thrust producible by the spacecraft thrusters in that
+ direction, the initial thrust controls can then be generated as a $20 \times 3$
+ matrix of uniform random numbers within that bound.
+
+ \subsection{Monotonic Basin Hopping}\label{mbh_subsection}
Outline the MBH algorithm, going into detail at each step. Mention the long-tailed PDF being
used and go into quite a bit of detail. Englander's paper on the MBH algorithm specifically
should be a good guide. Mention validation.
@@ -798,7 +902,7 @@
Talk a bit about why this work is valuable. Missions that could have benefited, missions that
this enables, etc.
- \section{Recommendations for Future Work}
+ \section{Recommendations for Future Work}\label{improvement_section}
Recommend future work, obviously. There are a \emph{ton} of opportunities for improvement
including parallelization, cluster computing, etc.
diff --git a/julia/Manifest.toml b/julia/Manifest.toml
index 5a9fc7d..0d4a718 100644
--- a/julia/Manifest.toml
+++ b/julia/Manifest.toml
@@ -99,6 +99,12 @@ git-tree-sha1 = "bdc0937269321858ab2a4f288486cb258b9a0af7"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.3.0"
+[[CodeTracking]]
+deps = ["InteractiveUtils", "UUIDs"]
+git-tree-sha1 = "9aa8a5ebb6b5bf469a7e0e2b5202cf6f8c291104"
+uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
+version = "1.0.6"
+
[[CodecBzip2]]
deps = ["Bzip2_jll", "Libdl", "TranscodingStreams"]
git-tree-sha1 = "2e62a725210ce3c3c2e1a3080190e7ca491f18d7"
@@ -329,6 +335,12 @@ git-tree-sha1 = "2f49f7f86762a0fbbeef84912265a1ae61c4ef80"
uuid = "7d188eb4-7ad8-530c-ae41-71a32a6d4692"
version = "0.3.4"
+[[JuliaInterpreter]]
+deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"]
+git-tree-sha1 = "b55aae9a2bf436fc797d9c253a900913e0e90178"
+uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a"
+version = "0.9.3"
+
[[Kaleido_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "2ef87eeaa28713cb010f9fb0be288b6c1a4ecd53"
@@ -390,6 +402,12 @@ version = "0.3.0"
[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+[[LoweredCodeUtils]]
+deps = ["JuliaInterpreter"]
+git-tree-sha1 = "6b0440822974cab904c8b14d79743565140567f6"
+uuid = "6f1432cf-f94c-5a45-995e-cdbf5db27b0b"
+version = "2.2.1"
+
[[METIS_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "2dc1a9fc87e57e32b1fc186db78811157b30c118"
@@ -559,6 +577,12 @@ git-tree-sha1 = "63ee24ea0689157a1113dbdab10c6cb011d519c4"
uuid = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
version = "1.9.0"
+[[Revise]]
+deps = ["CodeTracking", "Distributed", "FileWatching", "JuliaInterpreter", "LibGit2", "LoweredCodeUtils", "OrderedCollections", "Pkg", "REPL", "Requires", "UUIDs", "Unicode"]
+git-tree-sha1 = "2f9d4d6679b5f0394c52731db3794166f49d5131"
+uuid = "295af30f-e4ad-537b-8983-00126c2a3abe"
+version = "3.3.1"
+
[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
diff --git a/julia/Project.toml b/julia/Project.toml
index 8bfc3ca..c66204b 100644
--- a/julia/Project.toml
+++ b/julia/Project.toml
@@ -12,6 +12,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
PlotlyBase = "a03496cd-edff-5a9b-9e67-9cda94a718b5"
PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
SNOW = "105f6ee8-0889-46b1-9d4b-65e03cbc8f13"
SPICE = "5bab7191-041a-5c2e-a744-024b9c3a5062"
SymbolServer = "cf896787-08d5-524d-9de7-132aaa0cb996"
diff --git a/julia/src/Thesis.jl b/julia/src/Thesis.jl
index 0287548..5885b1b 100644
--- a/julia/src/Thesis.jl
+++ b/julia/src/Thesis.jl
@@ -7,7 +7,7 @@ module Thesis
using PlotlyJS
using Distributed
using SPICE
- using Dates: DateTime, Millisecond, Dates, Second, format, datetime2unix, unix2datetime
+ using Dates: DateTime, Millisecond, Dates, Second, format, datetime2julian, julian2datetime
try
furnsh("../../spice_files/naif0012.tls")
diff --git a/julia/src/nlp_solver.jl b/julia/src/nlp_solver.jl
index 0c869a9..0e8290f 100644
--- a/julia/src/nlp_solver.jl
+++ b/julia/src/nlp_solver.jl
@@ -2,30 +2,27 @@ using SNOW
export solve_mission
-const pos_scale = ones(3)
-# const vel_scale = 100. * ones(3)
-const vel_scale = ones(3)
-# const v∞_scale = norm(vel_scale)
-const v∞_scale = 1.
-const state_scale = [pos_scale; vel_scale ]
-
function constraint_bounds(guess::Mission_Guess)
low_constraint = Vector{Float64}()
high_constraint = Vector{Float64}()
for phase in guess.phases
- state_constraint = [100., 100., 100., 0.005, 0.005, 0.005] .* state_scale
+ # Constraints on the final state
+ state_constraint = [10_000., 10_000., 10_000., 0.05, 0.05, 0.05]
push!(low_constraint, (-1 * state_constraint)... )
push!(high_constraint, state_constraint... )
if phase != guess.phases[end]
- push!(low_constraint, -0.005*v∞_scale, 100.)
- push!(high_constraint, 0.005*v∞_scale, 1e7)
+ # Constraints on the v∞ and turning angle
+ push!(low_constraint, -0.05, 100.)
+ push!(high_constraint, 0.05, 1e7)
else
+ # Constraints on mass
push!(low_constraint, 0.)
push!(high_constraint, guess.start_mass - guess.sc.dry_mass)
end
end
-
+
+ # Constraints on C3 and V∞
push!(low_constraint, 0., 0.)
push!(high_constraint, 1e3, 1e3)
@@ -57,8 +54,8 @@ function solve_mission( guess::Mission_Guess,
# Establish initial conditions
v∞_out = x[2:4]
current_planet = Earth
- launch_date = unix2datetime(x[1])
- time = utc2et(format(launch_date,"yyyy-mm-ddTHH:MM:SS"))
+ launch_date = julian2datetime(x[1])
+ time = Dates.datetime2julian(launch_date)
start = state(current_planet, time, v∞_out, guess.start_mass)
final = zeros(7)
@@ -80,16 +77,16 @@ function solve_mission( guess::Mission_Guess,
thrusts = reshape(phase_params[8:end], (n,3))
# Propagate
- final = prop(thrusts, start, guess.sc, tof)[2]
+ final = prop(thrusts, start, guess.sc, tof)[end,:]
current_planet = flyby
- time += tof
+ time += tof/86400
goal = state(current_planet, time, v∞_in)[1:6]
start = state(current_planet, time, v∞_out, final[7])
# Do Checks
- g[1+8i:6+8i] .= (final[1:6] .- goal[1:6]) .* state_scale
+ g[1+8i:6+8i] .= (final[1:6] .- goal[1:6])
if flyby != flybys[end]
- g[7+8i] = (norm(v∞_out) - norm(v∞_in)) * v∞_scale
+ g[7+8i] = (norm(v∞_out) - norm(v∞_in))
δ = acos( ( v∞_in ⋅ v∞_out ) / ( norm(v∞_in) * norm(v∞_out) ) )
g[8+8i] = (current_planet.μ/(v∞_in ⋅ v∞_in)) * ( 1/sin(δ/2) - 1 ) - current_planet.r
else
@@ -100,7 +97,7 @@ function solve_mission( guess::Mission_Guess,
i += 1
end
- return 1.0
+ return -final[end]
catch e
@@ -114,7 +111,7 @@ function solve_mission( guess::Mission_Guess,
end
end
- max_time = datetime2unix(latest_arrival) - datetime2unix(launch_window[1])
+ max_time = (datetime2julian(latest_arrival) - datetime2julian(launch_window[1]))*86400
lower_x = lowest_mission_vector(launch_window, length(guess.phases), n)
upper_x = highest_mission_vector(launch_window, max_time, length(guess.phases), n)
num_constraints = 8*(length(guess.phases)-1) + 9
@@ -123,7 +120,7 @@ function solve_mission( guess::Mission_Guess,
"acceptable_constr_viol_tol" => 100tol,
"bound_relax_factor" => 0.,
"max_iter" => 100_000,
- "max_cpu_time" => 5. * length(guess.phases),
+ "max_cpu_time" => 10. * length(guess.phases),
"print_level" => print_level)
options = Options(solver=IPOPT(ipopt_options), derivatives=ForwardFD())
diff --git a/julia/src/types/bodies.jl b/julia/src/types/bodies.jl
index 4da940b..1cd6034 100644
--- a/julia/src/types/bodies.jl
+++ b/julia/src/types/bodies.jl
@@ -1,22 +1,106 @@
export state, period
struct Body
- μ::Float64
- r::Float64 # radius
+ μ::Real
+ r::Real # radius
color::String
id::Int # SPICE id
- a::Float64 # semi-major axis
+ a::Real # semi-major axis
line_color::AbstractString
name::AbstractString
end
period(b::Body) = 2π * √(b.a^3 / Sun.μ)
-function state(p::Body, t::DateTime, add_v∞::Vector{T}=[0., 0., 0.], add_mass::Float64=1e10) where T <: Real
- time = utc2et(format(t,"yyyy-mm-ddTHH:MM:SS"))
- [ spkssb(p.id, time, "ECLIPJ2000"); 0.0 ] + [ zeros(3); add_v∞; add_mass ]
+function state(p::Body, t::DateTime, add_v∞::AbstractVector=[0., 0., 0.], add_mass::Float64=1e10)
+ time = Dates.datetime2julian(t)
+ [ basic_ephemeris(p,time); 0.0 ] + [ zeros(3); add_v∞; add_mass ]
end
-function state(p::Body, t::S, add_v∞::Vector{T}=[0., 0., 0.], add_mass::Float64=1e10) where {T,S <: Real}
- [ spkssb(p.id, t, "ECLIPJ2000"); 0.0 ] + [ zeros(3); add_v∞; add_mass ]
+function state(p::Body, t::Real, add_v∞::AbstractVector=[0., 0., 0.], add_mass::Real=1e10)
+ [ basic_ephemeris(p,t); 0.0 ] + [ zeros(3); add_v∞; add_mass ]
+end
+
+"""
+Provides a really basic, but fast, ephemeris for most planets.
+"""
+function basic_ephemeris(planet::Body,date::Real)
+ century = (date-2451545)/36525
+ if planet.name == "Mercury"
+ table = [ 252.250906 149472.6746358 0.00000535 0.000000002 ;
+ 0.387098310 0 0 0 ;
+ 0.20563175 0.000020406 0.0000000284 0.00000000017 ;
+ 7.004986 -0.0059516 0.00000081 0.000000041 ;
+ 48.330893 -0.1254229 -0.00008833 -0.000000196 ;
+ 77.456119 0.1588643 -0.00001343 0.000000039]
+ elseif planet.name == "Venus"
+ table = [ 181.979801 58517.8156760 0.00000165 -0.000000002 ;
+ 0.72332982 0 0 0 ;
+ 0.00677188 -0.000047766 0.0000000975 0.00000000044 ;
+ 3.394662 -0.0008568 -0.00003244 0.000000010 ;
+ 76.679920 -0.2780080 -0.00014256 -0.000000198 ;
+ 131.563707 0.0048646 -0.00138232 -0.000005332]
+ elseif planet.name == "Earth"
+ table = [ 100.466449 35999.3728519 -0.00000568 0 ;
+ 1.000001018 0 0 0 ;
+ 0.01670862 -0.000042037 -0.0000001236 0.00000000004 ;
+ 0 0.0130546 -0.00000931 -0.000000034 ;
+ 174.873174 -0.2410908 0.00004067 -0.000001327 ;
+ 102.937348 0.3225557 0.00015026 0.000000478]
+ elseif planet.name == "Mars"
+ table = [ 355.433275 19140.2993313 0.00000261 -0.000000003 ;
+ 1.523679342 0 0 0 ;
+ 0.09340062 0.000090483 -0.0000000806 -0.00000000035 ;
+ 1.849726 -0.0081479 -0.00002255 -0.000000027 ;
+ 49.558093 -0.2949846 -0.00063993 -0.000002143 ;
+ 336.060234 0.4438898 -0.00017321 0.000000300]
+ elseif planet.name == "Jupiter"
+ table = [ 34.351484 3034.9056746 -0.00008501 0.000000004 ;
+ 5.202603191 0.0000001913 0 0 ;
+ 0.04849485 0.000163244 -0.0000004719 -0.00000000197 ;
+ 1.303270 -0.0019872 0.00003318 0.000000092 ;
+ 100.464441 0.1766828 0.00090387 -0.000007032 ;
+ 14.331309 0.2155525 0.00072252 -0.000004590]
+ elseif planet.name == "Saturn"
+ table = [ 50.077471 1222.1137943 0.00021004 -0.000000019 ;
+ 9.554909596 -0.0000021389 0 0 ;
+ 0.05550862 -0.000346818 -0.0000006456 0.00000000338 ;
+ 2.488878 0.0025515 -0.00004903 0.000000018 ;
+ 113.665524 -0.2566649 -0.00018345 0.000000357 ;
+ 93.056787 0.5665496 0.00052809 0.000004882]
+ elseif planet.name == "Uranus"
+ table = [ 314.055005 429.8640561 0.00030434 0.000000026 ;
+ 19.218446062 -0.0000000372 0.00000000098 0 ;
+ 0.04629590 -0.000027337 0.0000000790 0.00000000025 ;
+ 0.773196 0.0007744 0.00003749 -0.000000092 ;
+ 74.005947 0.5211258 0.00133982 0.000018516 ;
+ 173.005159 1.4863784 0.0021450 0.000000433]
+ elseif planet.name == "Neptune"
+ table = [ 304.348665 219.8833092 0.00030926 0.000000018 ;
+ 30.110386869 -0.0000001663 0.00000000069 0 ;
+ 0.00898809 0.000006408 -0.0000000008 -0.00000000005 ;
+ 1.769952 -0.0093082 -0.00000708 0.000000028 ;
+ 131.784057 1.1022057 0.00026006 -0.000000636 ;
+ 48.123691 1.4262677 0.00037918 -0.000000003]
+ elseif planet.name == "Pluto"
+ table = [ 238.92903833 145.20780515 0 0 ;
+ 39.48211675 -0.00031596 0 0 ;
+ 0.24882730 0.00005170 0 0 ;
+ 17.14001206 0.00004818 0 0 ;
+ 110.30393684 -0.01183482 0 0 ;
+ 224.06891629 -0.04062942 0 0]
+ else
+ return "No such planet"
+ end
+ poly = [1,century,century^2,century^3]
+ L,a,e,i,Ω,P = table * poly
+ L,i,Ω,P = [L,i,Ω,P] * (π/180)
+ ω = P-Ω
+ M = L-P
+ θ = M + (2e - 0.25e^3 + (5/96)*e^5)*sin(M) +
+ (1.25e^2 - (11/24)*e^4)*sin(2M) + ((13/12)*e^3 - (43/64)*e^5)*sin(3M) +
+ (103/96)*e^4*sin(4M) + (1097/960)*e^5*sin(5M)
+ oe = [a*AU,e,i,Ω,ω,θ]
+ return oe_to_xyz(oe, Sun.μ)
+
end
diff --git a/julia/src/types/mission.jl b/julia/src/types/mission.jl
index 8476fb1..02af410 100644
--- a/julia/src/types/mission.jl
+++ b/julia/src/types/mission.jl
@@ -49,7 +49,7 @@ This is the opposite of the function below. It takes a vector and any other nece
"""
function Mission_Guess(x::Vector{Float64}, sc::Sc, mass::Float64, flybys::Vector{Body})
# Variable mission params
- launch_date = unix2datetime(x[1])
+ launch_date = julian2datetime(x[1])
launch_v∞ = x[2:4]
# Try to intelligently determine n
@@ -78,7 +78,7 @@ information from the guess though.
"""
function Base.Vector(g::Mission_Guess)
result = Vector{Float64}()
- push!(result, datetime2unix(g.launch_date))
+ push!(result, datetime2julian(g.launch_date))
push!(result, g.launch_v∞...)
for phase in g.phases
push!(result,phase.v∞_in...)
@@ -91,7 +91,7 @@ end
function lowest_mission_vector(launch_window::Tuple{DateTime,DateTime}, num_phases::Int, n::Int)
result = Vector{Float64}()
- push!(result, datetime2unix(launch_window[1]))
+ push!(result, datetime2julian(launch_window[1]))
push!(result, -10*ones(3)...)
for i in 1:num_phases
push!(result, -10*ones(3)...)
@@ -104,7 +104,7 @@ end
function highest_mission_vector(launch_window::Tuple{DateTime,DateTime}, mission_length::Float64, num_phases::Int, n::Int)
result = Vector{Float64}()
- push!(result, datetime2unix(launch_window[2]))
+ push!(result, datetime2julian(launch_window[2]))
push!(result, 10*ones(3)...)
for i in 1:num_phases
push!(result, 10*ones(3)...)
@@ -133,14 +133,14 @@ function Mission(sc::Sc, mass::Float64, date::DateTime, v∞::Vector{Float64}, p
force=false)
# First do some checks to make sure that it's valid
if !force
- time = utc2et(format(date,"yyyy-mm-ddTHH:MM:SS"))
+ time = datetime2julian(date)
current_planet = Earth
start = state(current_planet, time, v∞, mass)
for phase in phases
#Propagate
- final = prop(phase.thrust_profile, start, sc, phase.tof)[2]
+ final = prop(phase.thrust_profile, start, sc, phase.tof)[end,:]
current_planet = phase.planet
- time += phase.tof
+ time += phase.tof/86400
goal = state(current_planet, time, phase.v∞_in)
start = state(current_planet, time, phase.v∞_out, final[7])
@@ -178,7 +178,7 @@ end
function Mission(x::Vector{Float64}, sc::Sc, mass::Float64, flybys::Vector{Body})
# Variable mission params
- launch_date = unix2datetime(x[1])
+ launch_date = julian2datetime(x[1])
launch_v∞ = x[2:4]
# Try to intelligently determine n
diff --git a/julia/src/utilities/laguerre-conway.jl b/julia/src/utilities/laguerre-conway.jl
index 12472a6..8652abf 100644
--- a/julia/src/utilities/laguerre-conway.jl
+++ b/julia/src/utilities/laguerre-conway.jl
@@ -1,4 +1,4 @@
-function laguerre_conway(state::Vector{<:Real}, time::Float64, primary::Body=Sun)
+function laguerre_conway(state::AbstractArray, time::Real, primary::Body=Sun)
μ = primary.μ
n = 5 # Choose LaGuerre-Conway "n"
diff --git a/julia/src/utilities/plotting.jl b/julia/src/utilities/plotting.jl
index ea2f6e2..b5f01de 100644
--- a/julia/src/utilities/plotting.jl
+++ b/julia/src/utilities/plotting.jl
@@ -59,27 +59,27 @@ end
"""
A fundamental function to generate a plot for a single orbit path
"""
-function gen_plot(path::Vector{Vector{Float64}};
+function gen_plot(path::AbstractArray;
label::Union{AbstractString, Nothing} = nothing,
color::AbstractString = random_color(),
markers=true)
# Plot the orbit
if label === nothing
- path_trace = scatter3d(;x=(path[1]),y=(path[2]),z=(path[3]),
+ path_trace = scatter3d(;x=(path[:,1]),y=(path[:,2]),z=(path[:,3]),
mode="lines", showlegend=false, line_color=color, line_width=3)
else
- path_trace = scatter3d(;x=(path[1]),y=(path[2]),z=(path[3]),
+ path_trace = scatter3d(;x=(path[:,1]),y=(path[:,2]),z=(path[:,3]),
mode="lines", name=label, line_color=color, line_width=3)
end
traces = [ path_trace ]
# Plot the markers
if markers
- push!(traces, scatter3d(;x=([path[1][1]]),y=([path[2][1]]),z=([path[3][1]]),
+ push!(traces, scatter3d(;x=([path[1,1]]),y=([path[1,2]]),z=([path[1,3]]),
mode="markers", showlegend=false,
marker=attr(color=color, size=3, symbol="circle")))
- push!(traces, scatter3d(;x=([path[1][end]]),y=([path[2][end]]),z=([path[3][end]]),
+ push!(traces, scatter3d(;x=([path[end,1]]),y=([path[end,2]]),z=([path[end,3]]),
mode="markers", showlegend=false,
marker=attr(color=color, size=3, symbol="diamond")))
end
@@ -222,7 +222,7 @@ end
"""
Generic plotting function for a bunch of paths. Useful for debugging.
"""
-function plot(paths::Vector{Vector{Vector{Float64}}};
+function plot(paths::Vector{Matrix{Real}};
labels::Union{Vector{String},Nothing}=nothing,
colors::Union{Vector{String},Nothing}=nothing,
markers::Bool=true,
@@ -243,7 +243,8 @@ function plot(paths::Vector{Vector{Vector{Float64}}};
push!(traces, gen_plot(Sun)[1]...)
# Determine layout details
- layout = standard_layout(limit, title)
+ # layout = standard_layout(limit, title)
+ layout = standard_layout(1e9, title)
# Plot
PlotlyJS.plot( PlotlyJS.Plot( traces, layout ) )
@@ -253,19 +254,20 @@ end
function plot(m::Union{Mission, Mission_Guess}; title::String="Mision Plot")
# First plot the earth
# Then plot phase plus planet phase
- time = m.launch_date
+ time = datetime2julian(m.launch_date)
mass = m.start_mass
current_planet = Earth
- earth_trace, limit = gen_plot(current_planet, time, year)
+ earth_trace, limit = gen_plot(current_planet, julian2datetime(time), year)
start = state(current_planet, time, m.launch_v∞, mass)
traces = [ earth_trace... ]
i = 1
for phase in m.phases
# First do the path
- path, final = prop(phase.thrust_profile, start, m.sc, phase.tof, interpolate=true)
+ path = prop(phase.thrust_profile, start, m.sc, phase.tof)
+ final = path[end,:]
mass = final[7]
- time += Second(floor(phase.tof))
+ time += phase.tof/86400
current_planet = phase.planet
start = state(current_planet, time, phase.v∞_out, mass)
path_trace, new_limit = gen_plot(path, label="Phase "*string(i))
@@ -273,7 +275,7 @@ function plot(m::Union{Mission, Mission_Guess}; title::String="Mision Plot")
limit = max(limit, new_limit)
# Then the planet
- planet_trace, new_limit = gen_plot(current_planet, time, -phase.tof)
+ planet_trace, new_limit = gen_plot(current_planet, julian2datetime(time), -phase.tof)
push!(traces, planet_trace...)
limit = max(limit, new_limit)
diff --git a/julia/src/utilities/propagator.jl b/julia/src/utilities/propagator.jl
index f74c8b0..29a636b 100644
--- a/julia/src/utilities/propagator.jl
+++ b/julia/src/utilities/propagator.jl
@@ -3,83 +3,52 @@ export prop
"""
Maximum ΔV that a spacecraft can impulse for a given single time step
"""
-function max_ΔV(duty_cycle::Float64,
+function max_ΔV(duty_cycle::Real,
num_thrusters::Int,
- max_thrust::Float64,
- tf::Float64,
- t0::Float64,
- mass::T) where T <: Real
- return duty_cycle*num_thrusters*max_thrust*(tf-t0)/mass
-end
-
-"""
-A convenience function for using spacecraft. Note that this function outputs a sc instead of a mass
-"""
-function prop_one(ΔV_unit::Vector{<:Real},
- state::Vector{<:Real},
- craft::Sc,
- time::Float64,
- primary::Body=Sun)
-
- ΔV = max_ΔV(craft.duty_cycle, craft.num_thrusters, craft.max_thrust, time, 0., state[7]) * ΔV_unit
- halfway = laguerre_conway(state, time/2, primary) + [zeros(3); ΔV]
- final = laguerre_conway(halfway, time/2, primary)
- return [final; state[7] - mfr(craft)*norm(ΔV_unit)*time]
-
+ max_thrust::Real,
+ tf::Real,
+ t0::Real)
+ return duty_cycle*num_thrusters*max_thrust*(tf-t0)
end
"""
The propagator function
"""
-function prop(ΔVs::Matrix{T},
- state::Vector{Float64},
+function prop(ΔVs::AbstractArray,
+ state::AbstractArray,
craft::Sc,
- time::Float64,
- primary::Body=Sun;
- interpolate::Bool=false) where T <: Real
+ time::Real;
+ primary::Body=Sun)
- size(ΔVs)[2] == 3 || throw(ΔVsize_Error())
+ size(ΔVs)[2] == 3 || throw(ΔVsize_Error(size(ΔVs)))
n = size(ΔVs)[1]
-
- states = [ Vector{T}(),Vector{T}(),Vector{T}(),Vector{T}(),Vector{T}(),Vector{T}(),Vector{T}() ]
-
- for i in 1:7 push!(states[i], state[i]) end
+ output = Array{Real,2}(undef,n,7)
+ current = copy(state)
+ full_thrust = max_ΔV(craft.duty_cycle, craft.num_thrusters, craft.max_thrust, time/n, 0.)
for i in 1:n
- if interpolate
- interpolated_state = copy(state)
- for j in 1:49
- interpolated_state = [laguerre_conway(interpolated_state, time/50n, primary); 0.0]
- for k in 1:7 push!(states[k], interpolated_state[k]) end
- end
+ halfway = laguerre_conway(current[1:6], time/(2n), primary)
+ if time > 0
+ halfway += [zeros(3); full_thrust * ΔVs[i,:]] / current[7]
+ current[7] -= mfr(craft) * norm(ΔVs[i,:]) * time/n
+ else
+ current[7] -= mfr(craft) * norm(ΔVs[i,:]) * time/n
+ halfway += [zeros(3); full_thrust * ΔVs[i,:]] / current[7]
end
- try
- state = prop_one(ΔVs[i,:], state, craft, time/n, primary)
- catch e
- if isa(e,PropOne_Error)
- for val in e.ΔV_unit
- # If this isn't true, then we just let it slide
- if abs(val) > 1.0001
- rethrow()
- end
- end
- else
- rethrow()
- end
- end
- for j in 1:7 push!(states[j], state[j]) end
+ current[1:6] = laguerre_conway(halfway, time/(2n), primary)
+ output[i,:] = current
end
- state[7] > craft.dry_mass || throw(Mass_Error(state[7]))
+ current[7] > craft.dry_mass || throw(Mass_Error(current[7]))
- return states, state
+ return output
end
"""
Convenience function for propagating a state with no thrust
"""
-prop(x::Vector{Float64}, t::Float64, p::Body=Sun) = prop(zeros(1000,3), [x;1.], no_thrust, t, p)[1]
+prop(x::AbstractArray, t::Real, p::Body=Sun) = prop(zeros(1000,3), x, no_thrust, t, p)
"""
This is solely for the purposes of getting the final state of a mission or guess
@@ -91,7 +60,7 @@ function prop(m::Union{Mission, Mission_Guess})
final = zeros(7)
for phase in m.phases
- final = prop(phase.thrust_profile, start, m.sc, phase.tof)[2]
+ final = prop(phase.thrust_profile, start, m.sc, phase.tof)[end,:]
mass = final[7]
current_planet = phase.planet
time += Second(floor(phase.tof))
diff --git a/julia/test/missions/nlp_2_phase b/julia/test/missions/nlp_2_phase
index 86788f8..fd1dde5 100644
--- a/julia/test/missions/nlp_2_phase
+++ b/julia/test/missions/nlp_2_phase
@@ -7,23 +7,23 @@ Spacecraft: bepi
num_thrusters: 2
duty_cycle: 0.9
Launch Mass: 3600.0 kg
-Launch Date: 2021-10-31T23:59:59.769
-Launch V∞: [1.0462564431708832, -2.520754077506101, -2.301043507881933] km/s
+Launch Date: 2021-10-31T22:00:16.400
+Launch V∞: [1.0097399011270518, -2.6447886914618515, -2.5394649406128433] km/s
Phase 1:
Planet: Venus
- V∞_in: [0.4205842845952983, -3.5241183574614556, 4.879873784266574] km/s
- V∞_out: [2.657685712690303, -4.920354075306563, 2.266369032728228] km/s
- time of flight: 1.2614399691127092e7 seconds
- arrival date: 2022-03-26T23:59:58.769
- thrust profile: [0.08791530712557315 -0.009101317973001748 0.010254083007822207; 0.07287521745776969 -0.011381949725375708 0.051332448819202944; 0.04543116258472458 -0.01531732267940295 0.08273458014460809; 0.011034188560216448 -0.020310595889307323 0.10100599446320643; -0.009607377273254756 -0.02382633371411493 0.1081982710061048; -0.028510355701058584 -0.030611944159923403 0.11101434996597628; -0.039048295910837236 -0.03820244320795566 0.11279184708362161; -0.048252651406879885 -0.04651039104792323 0.11301481830933677; -0.05288674350674192 -0.0561882917884983 0.11293454239992265; -0.05458628096513579 -0.06643801004024381 0.11223764822140121; -0.053564283975236535 -0.07751051107973227 0.11061646758637526; -0.049500823222769094 -0.0891415249668692 0.10769319474620971; -0.04185476746663511 -0.10066361474827486 0.10355236756997019; -0.03087924291992901 -0.1112479139000727 0.09781871133696314; -0.015564164113106642 -0.11996802006040436 0.09064682025194883; 0.001781758545597473 -0.12613079909578606 0.0809801955196938; 0.021345925406049154 -0.12876486287710442 0.06930029766487243; 0.041502826765638504 -0.12763804654948588 0.054677487677216745; 0.06024411885991524 -0.12236183464053872 0.038231302757203295; 0.07560839256024353 -0.114157581119008 0.018886579096999003] %
+ V∞_in: [0.6971321781569763, -3.000415336232889, 5.257930893162622] km/s
+ V∞_out: [3.155471614002537, -4.846243818499916, 1.9261368837603878] km/s
+ time of flight: 1.261440000905646e7 seconds
+ arrival date: 2022-03-26T22:00:16.400
+ thrust profile: [0.034282174989117516 0.030630284530489906 0.0516767016545959; 0.04055472535721235 0.027376197088858283 0.04809039664320572; 0.043660143604191994 0.03047637966598728 0.04358760046457853; 0.04125099543115894 0.03016179866323979 0.04499932191472562; 0.04495441021193453 0.03348083036100891 0.03861183062100992; 0.04710965203705168 0.02635602095857307 0.04054560175736578; 0.05033947480940243 0.03062684297509654 0.032922298467773295; 0.04663650635412875 0.030369216537824293 0.03703971691255565; 0.048565701285880614 0.029118981813448557 0.034508738912923845; 0.05056448644255328 0.028722504694661054 0.031034420379757226; 0.04932798248085456 0.03007397062288681 0.030824783232637425; 0.047996060872376366 0.028952951818463326 0.032857369566278626; 0.04951821298908887 0.028758116977346378 0.027958339954661; 0.0489313400699588 0.03046347487498191 0.027859841906091212; 0.04520420340375004 0.031670198439770064 0.030514442243397236; 0.04566034232248259 0.034381957552538395 0.024900845759794794; 0.04082535412813151 0.03357236769501879 0.03170011587741194; 0.04621217362755022 0.03182636410505487 0.02410773528455805; 0.045645717216652795 0.032921517920520164 0.02079135127615303; 0.042597542211632096 0.03558968396172082 0.021282994074205266; 0.03731044485129105 0.03994289623933655 0.02130495307758324; 0.04247833588674327 0.03457796735712127 0.01899569141913135; 0.04312850336659755 0.03399248745544899 0.013149159955936978; 0.03754604612989521 0.03763011819702751 0.018318693954420588; 0.03535094668642273 0.036858320028499134 0.02199684830247042; 0.028891076662331875 0.04324155785176078 0.014413619206886816; 0.03267849036426984 0.0401634806637448 0.011840587979407397; 0.024706839495172658 0.043042163807245896 0.01826522977160073; 0.028601237657836467 0.041588848578667854 0.014411946517623132; 0.028228090581187873 0.04015749743865513 0.01508529741887872; 0.024783781293558266 0.04189396738290361 0.014616887877185162; 0.024764235231400287 0.0412235929856152 0.014206631004299218; 0.022700981741250053 0.04285522888573239 0.008652897464076466; 0.023149708367887163 0.04204784290805528 0.011482287191821549; 0.030060665695868752 0.03682561335933804 0.01406740384301527; 0.020327198179371554 0.03967929127546447 0.01907556971774192; 0.01948576963473395 0.04150171458022886 0.01440277138047399; 0.014609183554173198 0.04193091541988154 0.0190142349738259; 0.017338241200120032 0.0405327064144073 0.019397828837718817; 0.013519064235294736 0.043194207588669445 0.012618271521276186; 0.01416527483071525 0.042360112796518024 0.018565046588499192; 0.016192504354863486 0.04321790744316709 0.014857969976184007; 0.005926517746692747 0.04385012579593121 0.01848227965737381; 0.007881472501745706 0.045008039792868576 0.01571801630208125; 0.006544092936553959 0.04561895672702495 0.016342515520064185; 0.0018483989093118227 0.046533973409849064 0.014988164886431584; -0.000568247332158013 0.047073416094890685 0.01286929122897256; -0.008004594524198529 0.046599777614930456 0.013010850227684606; -0.00696671229499173 0.048339673337854945 0.010999881717700421; -0.0070218437672999665 0.04911375570173229 0.010996777473604348] %
Phase 2:
Planet: Mars
- V∞_in: [3.866657617879469, -3.8561104436565112, -1.1528388074396114] km/s
- V∞_out: [-0.04156826938628659, 0.04537147320396694, 0.010962101613372877] km/s
- time of flight: 1.3305600070108933e7 seconds
- arrival date: 2022-08-27T23:59:58.769
- thrust profile: [0.12294683695617006 -0.14210501591303906 -0.015953927996763824; 0.1346036754636979 -0.11703322342276389 -0.02118233265872474; 0.14035562458965456 -0.09253436643679969 -0.026037109750586127; 0.14159974389839877 -0.07034252588185486 -0.02965559917950358; 0.13983758574258884 -0.0508280123903436 -0.03226072974284784; 0.13628355001371623 -0.03381714643542805 -0.033706183284138075; 0.13174907912062978 -0.01944710462973476 -0.03378159713957984; 0.12690191547644916 -0.006666552590772735 -0.032314321099804075; 0.12190989012800546 0.005350412631176013 -0.02896415432946389; 0.11640547480173373 0.01766009561198736 -0.024295908236754368; 0.11004779022898785 0.029624240108957346 -0.018629744776725948; 0.10255142565747595 0.041391197147168296 -0.011260773047897207; 0.09257939116538691 0.054055648527399226 -0.0018819573006649628; 0.07862129598678515 0.06762185940978795 0.006827208192964654; 0.06184044671572194 0.07873199678299247 0.015671300335402687; 0.04256729407158805 0.08711502955932304 0.023547740065060118; 0.022644089962366266 0.0922573542104121 0.029300773045221138; 0.0012772594401787602 0.09453861142625342 0.0335188440659034; -0.017213511457067864 0.09403439603239908 0.036105272081528195; -0.03466637624188391 0.0910940423628138 0.03867124120601387] %
+ V∞_in: [3.934071361708631, -4.0556866988998195, -0.9569218994864468] km/s
+ V∞_out: [2.835079207012675, -3.038506396757051, -0.8239067509398897] km/s
+ time of flight: 1.3305600002135515e7 seconds
+ arrival date: 2022-08-27T22:00:16.400
+ thrust profile: [0.07875872819429838 -0.07602928888291517 -0.002461076427925589; 0.08368497242841927 -0.06670986115534407 0.001665869512245407; 0.08577981435177857 -0.05793119355590426 0.011165188393882227; 0.08758545440549516 -0.0339816665457332 0.026794441364176573; 0.08293133779344392 -0.02340157704539268 0.034551607127364624; 0.07757395448326561 -0.006354368804274271 0.03802615353405709; 0.0644275187408338 0.00817115132820835 0.04477293215146323; 0.052643603336670175 0.021208603088669686 0.04641877730052604; 0.0422696791291819 0.03249686419539815 0.04293455453504036; 0.0299947244754212 0.0342146829206795 0.04610828643246764; 0.028375921792895057 0.031929570942726804 0.04880758017415033; 0.021457075341902082 0.03512409953676916 0.047848639128187716; 0.014050820576587325 0.03598753239044044 0.048208234639668036; 0.01239374458074181 0.03519088938079009 0.04991578465780201; 0.005554038390401058 0.03590306275800209 0.04943448284625346; 0.004594710645977571 0.03473960951290201 0.05085953544892566; 0.004962419167461687 0.03677569856463667 0.05029354676804025; 0.0007318114245317572 0.03710676170960615 0.050422609329575686; -0.0008906817624837742 0.03770220933797602 0.05041455076147805; -0.006778334873847526 0.03779829345387088 0.04962125740173904; -0.007046493202932243 0.03568043355312452 0.051960594399799505; -0.008049628911291552 0.03913802154188755 0.05028104660474625; -0.010537903266724359 0.03672325453597655 0.051884804336631986; -0.00818351261490794 0.038741730202093635 0.05210658416215128; -0.010666656787731375 0.03849972475725258 0.05243476797310733; -0.012555727410020852 0.0382137309957785 0.05268364578100002; -0.009820005328087349 0.03916769798664862 0.05359353372069135; -0.012387176669361927 0.03970296584000698 0.05331339060382668; -0.01244834483921497 0.041886262658702045 0.05253893772334037; -0.013444072160084112 0.04262744032016184 0.05253898626998434; -0.012523862311046558 0.04142301768047668 0.05454027471095643; -0.012895225834978194 0.042511978835118515 0.0544877324498131; -0.012679086383971097 0.04485961538708221 0.05361843203762159; -0.01233563279407355 0.04412582031595734 0.05489553351056654; -0.011291838361870397 0.04737557615034447 0.05363956736166156; -0.011230502857283789 0.04704594850205801 0.05467457973688944; -0.008797763456055051 0.04928744387907065 0.05415511568853677; -0.008048802844087589 0.05075747730843418 0.0538825677403911; -0.0056970113123549744 0.05224183281724393 0.05369505177198884; -0.006070834938122434 0.052936925948144475 0.05379731537648658; -0.0024653651357053295 0.05567255160163994 0.05213940324181525; -0.000733453786943943 0.057067189219935754 0.05145199179114529; 0.0012456777227569263 0.059943358064031824 0.04907287600209174; 0.004373513720293685 0.06169575905345249 0.047206183474067936; 0.009873790970193171 0.06486640320654036 0.04255101782330187; 0.013319573203495716 0.06692522008245975 0.03887179390451873; 0.018201455983236366 0.06835502892414509 0.034298233864195185; 0.023662275014910684 0.06904022589813187 0.029150000528514425; 0.02774195908046571 0.07005784655018363 0.02165291283778826; 0.031997497603024734 0.06946489515884077 0.015502921841459603] %
-Mass Used: 121.86921799846004 kg
-Launch C3: 12.743654889305812 km²/s²
-||V∞_in||: 5.58118860132149 km/s
+Mass Used: 61.61570260672124 kg
+Launch C3: 14.463364075014354 km²/s²
+||V∞_in||: 5.730725224643533 km/s
diff --git a/julia/test/nlp_solver.jl b/julia/test/nlp_solver.jl
index a6f78ad..642fcc4 100644
--- a/julia/test/nlp_solver.jl
+++ b/julia/test/nlp_solver.jl
@@ -4,72 +4,37 @@
println("Testing NLP solver")
- # Test the optimizer for a one-phase mission
- # The lambert's solver said this should be pretty valid
- launch_window = DateTime(1992,11,1), DateTime(1992,12,1)
- latest_arrival = DateTime(1993,6,1)
- leave, arrive = DateTime(1992,11,19), DateTime(1993,4,1)
- test_leave = DateTime(1992,11,12)
- earth_state = state(Earth, leave)
- venus_state = state(Venus, arrive)
- v∞_out, v∞_in, tof = Thesis.lamberts(Earth, Venus, leave, arrive)
- # We can get the thrust profile and tof pretty wrong and still be ok
- phase = Phase(Venus, 1.1v∞_in, v∞_in, 0.9*tof, 0.1*ones(20,3))
- guess = Mission_Guess(bepi, 3_600., test_leave, 0.9*v∞_out, [phase])
- m = solve_mission(guess, launch_window, latest_arrival, 200., 20., verbose=true)
- @test typeof(m) == Mission
-
- # Now we can plot the results to check visually
- p = plot(m, title="NLP Test Solution")
- savefig(p,"../plots/nlp_test_1_phase.html")
- store(m, "missions/nlp_1_phase")
-
- # Now we can look at a slightly more complicated trajectory
- flybys = [Earth, Venus, Mars]
+ # This will be the only test to verify rather than testing different complexities
+ n = 50
+ flybys = [Earth, Venus, Mars]#, Jupiter]
launch_window = DateTime(2021,10,1), DateTime(2021,12,1)
- latest_arrival = DateTime(2023,1,1)
- dates = [DateTime(2021,11,1), DateTime(2022,3,27), DateTime(2022,8,28)]
+ latest_arrival = DateTime(2028,1,1)
+ dates = [DateTime(2021,11,1),
+ DateTime(2022,3,27),
+ DateTime(2022,8,28)]
+ # DateTime(2028,3,1)]
phases = Vector{Phase}()
launch_v∞, _, tof1 = Thesis.lamberts(flybys[1], flybys[2], dates[1], dates[2])
for i in 1:length(dates)-2
- v∞_out1, v∞_in1, tof1 = Thesis.lamberts(flybys[i], flybys[i+1], dates[i], dates[i+1])
- v∞_out2, v∞_in2, tof2 = Thesis.lamberts(flybys[i+1], flybys[i+2], dates[i+1], dates[i+2])
- push!(phases, Phase(flybys[i+1], v∞_in1, v∞_out2, tof1, 0.01*ones(20,3)))
+ v∞_in1, tof1 = Thesis.lamberts(flybys[i], flybys[i+1], dates[i], dates[i+1])[2:3]
+ v∞_out2 = Thesis.lamberts(flybys[i+1], flybys[i+2], dates[i+1], dates[i+2])[1]
+ push!(phases, Phase(flybys[i+1], v∞_in1, v∞_out2, tof1, 0.01*ones(50,3)))
end
- v∞_out, v∞_in, tof = Thesis.lamberts(flybys[end-1], flybys[end], dates[end-1], dates[end])
- push!(phases, Phase(flybys[end], v∞_in, v∞_in, tof, 0.01*ones(20,3)))
+ v∞_in, tof = Thesis.lamberts(flybys[end-1], flybys[end], dates[end-1], dates[end])[2:3]
+ push!(phases, Phase(flybys[end], v∞_in, v∞_in, tof, 0.01*ones(50,3)))
guess = Mission_Guess(bepi, 3_600., dates[1], launch_v∞, phases)
- m = solve_mission(guess, launch_window, latest_arrival, 200., 20., verbose=true)
+ m = solve_mission(guess,
+ launch_window,
+ latest_arrival,
+ 200.,
+ 20.,
+ verbose=true,
+ print_level=4)
@test typeof(m) == Mission
- p = plot(m, title="NLP Test Solution (2 Phases)")
- savefig(p,"../plots/nlp_test_2_phase.html")
- store(m, "missions/nlp_2_phase")
-
- # Here is the final, most complicated, trajectory to test
- # Ignoring for now as the initial guess makes the test take too long to converge with mbh settings
- # flybys = [Earth, Venus, Earth, Mars, Earth, Jupiter]
- # launch_window = DateTime(2023,1,1), DateTime(2024,1,1)
- # latest_arrival = DateTime(2031,1,1)
- # dates = [DateTime(2023,5,23),
- # DateTime(2023,10,21),
- # DateTime(2024,8,24),
- # DateTime(2025,2,13),
- # DateTime(2026,11,22),
- # DateTime(2032,1,1)]
- # phases = Vector{Phase}()
- # launch_v∞, _, tof1 = Thesis.lamberts(flybys[1], flybys[2], dates[1], dates[2])
- # for i in 1:length(dates)-2
- # v∞_out1, v∞_in1, tof1 = Thesis.lamberts(flybys[i], flybys[i+1], dates[i], dates[i+1])
- # v∞_out2, v∞_in2, tof2 = Thesis.lamberts(flybys[i+1], flybys[i+2], dates[i+1], dates[i+2])
- # push!(phases, Phase(flybys[i+1], 1.02v∞_in1, 0.98v∞_out2, 1.02tof1, 0.02*ones(20,3)))
- # end
- # v∞_out, v∞_in, tof = Thesis.lamberts(flybys[end-1], flybys[end], dates[end-1], dates[end])
- # push!(phases, Phase(flybys[end], v∞_in, v∞_in, tof, 0.01*ones(20,3)))
- # guess = Mission_Guess(bepi, 3_600., dates[1], launch_v∞, phases)
- # m = solve_mission(guess, launch_window, latest_arrival, 200., 20., verbose=true)
- # @test typeof(m) == Mission
- # p = plot(m, title="NLP Test Solution (5 Phases)")
- # savefig(p,"../plots/nlp_test_5_phase.html")
- # store(m, "missions/nlp_5_phase")
+ if typeof(m) == Mission
+ p = plot(m, title="NLP Test Solution (2 Phases)")
+ savefig(p,"../plots/nlp_test_2_phase.html")
+ store(m, "missions/nlp_2_phase")
+ end
end
diff --git a/julia/test/propagator.jl b/julia/test/propagator.jl
index bcc6210..1ea57c1 100644
--- a/julia/test/propagator.jl
+++ b/julia/test/propagator.jl
@@ -1,21 +1,27 @@
@testset "Propagator" begin
- println("Testing propagator")
+ using PlotlyJS: savefig
- using Thesis: prop_one
+ println("Testing Propagator")
# Set up
start_mass = 10_000.
start = gen_orbit(rand(.5year : hour : 2year), start_mass)
- stepsize = rand(hour : second : 6hour)
+ time = rand(0.2year : second : 0.5year)
+ n = 1_000
- # Test that Laguerre-Conway is the default propagator for spacecrafts
- 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 Propagation works
+ states = prop(zeros(n,3), start, no_thrust, time)
+ @test states[1,:] != states[end,:]
# Test that mass is reduced properly
- state = prop_one([1., 0., 0.], start, bepi, stepsize)
- @test state[7] == start_mass - mfr(bepi)*stepsize
-
+ states = prop(ones(n,3)/√3, start, bepi, time)
+ @test states[end,7] ≈ start_mass - mfr(bepi)*time
+
+ # Test that the propagator works backwards
+ backwards = prop(ones(n,3)/√3, states[end,:], bepi, -time)
+ p = plot([states, backwards])
+ savefig(p,"../plots/prop_back.html")
+ @test backwards[end,:] ≈ start
+
end
diff --git a/julia/test/runtests.jl b/julia/test/runtests.jl
index 5cc4122..b80aea9 100644
--- a/julia/test/runtests.jl
+++ b/julia/test/runtests.jl
@@ -3,7 +3,7 @@ using Random
using LinearAlgebra
using SPICE
using Thesis
-using Dates: DateTime, Millisecond, Dates, Second, format, datetime2unix, unix2datetime
+using Dates: DateTime, Millisecond, Dates, Second, format, datetime2julian, julian2datetime
try
furnsh("../../spice_files/naif0012.tls")
@@ -14,10 +14,10 @@ catch
end
@testset "All Tests" begin
- include("plotting.jl")
- include("laguerre-conway.jl")
+ # include("plotting.jl")
+ # include("laguerre-conway.jl")
include("propagator.jl")
- include("nlp_solver.jl")
- include("mbh.jl")
- include("genetic_algorithm.jl")
+ # include("nlp_solver.jl")
+ # include("mbh.jl")
+ # include("genetic_algorithm.jl")
end