diff --git a/LaTeX/fig/EMJS_plot.png b/LaTeX/fig/EMJS_plot.png index e745295..d035262 100644 Binary files a/LaTeX/fig/EMJS_plot.png and b/LaTeX/fig/EMJS_plot.png differ diff --git a/LaTeX/fig/EMJS_plot_noplanets.png b/LaTeX/fig/EMJS_plot_noplanets.png index 50229e2..2aa8823 100644 Binary files a/LaTeX/fig/EMJS_plot_noplanets.png and b/LaTeX/fig/EMJS_plot_noplanets.png differ diff --git a/LaTeX/fig/EMJS_thrust_components.png b/LaTeX/fig/EMJS_thrust_components.png index 19f0cca..a39148c 100644 Binary files a/LaTeX/fig/EMJS_thrust_components.png and b/LaTeX/fig/EMJS_thrust_components.png differ diff --git a/LaTeX/fig/EMJS_thrust_mag.png b/LaTeX/fig/EMJS_thrust_mag.png index 342a067..4b8d4d4 100644 Binary files a/LaTeX/fig/EMJS_thrust_mag.png and b/LaTeX/fig/EMJS_thrust_mag.png differ diff --git a/LaTeX/fig/EMS_plot.png b/LaTeX/fig/EMS_plot.png index f314a5d..536519b 100644 Binary files a/LaTeX/fig/EMS_plot.png and b/LaTeX/fig/EMS_plot.png differ diff --git a/LaTeX/fig/EMS_plot_noplanets.png b/LaTeX/fig/EMS_plot_noplanets.png index d027f32..6509c23 100644 Binary files a/LaTeX/fig/EMS_plot_noplanets.png and b/LaTeX/fig/EMS_plot_noplanets.png differ diff --git a/LaTeX/fig/EMS_thrust_components.png b/LaTeX/fig/EMS_thrust_components.png index 3a22131..e5447e7 100644 Binary files a/LaTeX/fig/EMS_thrust_components.png and b/LaTeX/fig/EMS_thrust_components.png differ diff --git a/LaTeX/fig/EMS_thrust_mag.png b/LaTeX/fig/EMS_thrust_mag.png index d27910a..cfc62dd 100644 Binary files a/LaTeX/fig/EMS_thrust_mag.png and b/LaTeX/fig/EMS_thrust_mag.png differ diff --git a/LaTeX/fig/porkchop.png b/LaTeX/fig/porkchop.png index 347b776..a4a1299 100644 Binary files a/LaTeX/fig/porkchop.png and b/LaTeX/fig/porkchop.png differ diff --git a/LaTeX/thesis.tex b/LaTeX/thesis.tex index 3a3e37e..00ddace 100644 --- a/LaTeX/thesis.tex +++ b/LaTeX/thesis.tex @@ -1155,9 +1155,9 @@ 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. + numbers and normalized. This vector is then multiplied by a uniform random number + between 0 and the 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 @@ -1288,17 +1288,17 @@ a tuning parameter to determine the size of the tails and width of the distribution set to $1.01$, but easily tunable. - The perturbation function, then steps through each parameter of the mission, - generating a new mission guess with the parameters modified by the above Pareto - distribution. After this perturbation, the NLP solver is then called again to find - a valid solution in the vicinity of this new guess. If the solution is better than - the current basin solution, it replaces that value and the drill counter is reset to + The perturbation function then steps through each parameter of the mission, + generating a new guess with the parameters modified by the Pareto distribution. + After this perturbation, the NLP solver is then called again to find a valid + solution in the vicinity of this new guess. If the solution is better than the + current basin solution, it replaces that value and the drill counter is reset to zero. If it is better than the current total best, it replaces that value as well. Otherwise, the drill counter increments and the process is repeated. Therefore, the drill patience allows the mission designer to determine a maximum number of - iterations to perform without any improvements in a row before ending a given drill - loop. This process can be repeated essentially ''search patience`` number of times - in order to fully traverse all basins. + iterations to perform without improvement in a row before ending the drill loop. + This process can be repeated essentially ''search patience`` number of times in + order to fully traverse all basins. \chapter{Results Analysis} \label{results} @@ -1436,8 +1436,11 @@ with Mars after three and one half years to rendezvous in mid-December 2027. Unfortunately, the launch energy required to effectively used the gravity assist with Mars at this time is quite high. The $C_3$ value was found to be $60.4102$ kilometers - per second squared. However, for this phase, the thrusters are almost entirely turned - off, allowing for a nearly-natural trajectory to Mars rendezvous. + per second squared. However, for this phase, the thrust magnitudes are quite low, + raising slowly only as the spacecraft approaches Mars, allowing for a nearly-natural + trajectory to Mars rendezvous. Note also that the in-plane thrust direction was neither + zero nor $\pi$, implying that these thrusts were steering thrusts rather than + momentum-increasing thrusts. \begin{figure}[H] \centering @@ -1518,11 +1521,11 @@ previous trajectory. However, this time the launch energy is considerably lower, with a $C_3$ value of only $40.4386$ kilometer per second squared. Rather than employ an almost entirely natural coasting arc to Mars, however, this trajectory performs some thrusting - at about the apoapsis point of its orbit in order to raise the periapsis enough to - rendezvous at roughly the same incidence angle in Mars' orbit, but one revolution ahead. - In this case, the launch was a bit earlier, occurring in November of 2023, with the Mars - flyby occurring in mid-April of 2026. This will prove to be helpful in comparison with - the other result, as this mission profile is much longer. + almost entirely in the velocity direction, increasing its orbital energy in order to + achieve the same Mars rendezvous. In this case, the launch was a bit earlier, occurring + in November of 2023, with the Mars flyby occurring in mid-April of 2026. This will prove + to be helpful in comparison with the other result, as this mission profile is much + longer. \begin{figure}[H] \centering @@ -1565,9 +1568,11 @@ Finally, this mission also has a third phase. The Jupiter flyby provides quite a strong $\Delta V$ for the spacecraft, allowing the following phase to largely be a coasting arc - to Saturn almost one revolution later. Because of this long coasting period, the mission - length increases considerably during this leg, arriving at Saturn in December of 2037, - over 8 years after the Jupiter flyby. + to Saturn almost one revolution later. During the most efficient part of the journey, + some thrust in the velocity direction accounts for a little bit of orbit-raising, but + the phase is largely natural. Because of this long coasting period, the mission length + increases considerably during this leg, arriving at Saturn in December of 2037, over 8 + years after the Jupiter flyby. However, there are many advantages to this approach relative to the other trajectory. While the fuel use is also slightly higher at $530.668$ kilograms, plenty of payload diff --git a/julia/ingest_mission.jl b/julia/ingest_mission.jl index 4f9ee56..e71f392 100644 --- a/julia/ingest_mission.jl +++ b/julia/ingest_mission.jl @@ -2,6 +2,7 @@ using Thesis using SPICE using LinearAlgebra using Dates +using PlotlyJS try furnsh("../../spice_files/naif0012.tls") @@ -48,16 +49,59 @@ end sorted_missions1 = sort(missions, by=cost_1) sorted_missions2 = sort(missions, by=cost_2) -println(get_id.(sorted_missions1[1:10])) -println(get_id.(sorted_missions2[1:10])) +println(get_id.(sorted_missions1[1:5])) +println(get_id.(sorted_missions2[1:5])) -display_mission = improve(sorted_missions2[1], n=7) -# display(plot(display_mission, +# long, short = improve.(sorted_missions2[1:2], n=7)#, use_cost=true) +display_mission = long + +# Display Orbit Plot +# display(Thesis.plot(long, # title="Sample Algorithm Result Mission", # mode="light", # planet_colors=false, # phase_colors=["#F00","#F0F","#6AF"], -# markers=false)); -display(plot_thrust(display_mission, title="Thrust Magnitude vs Time", mode="light")) -display(plot_thrust_components(display_mission, title="Thrust Components vs Time", mode="light")) +# markers=false, +# camera=((-0.2, -0.3, -0.2),(1.75,-1.75,1.5)))); + +# Save Orbit Plot +savefig(Thesis.plot(short, title="EMS Mission Profile", mode="light", planet_colors=false, + phase_colors=["#F00","#F0F","#6AF"], markers=false, + camera=((-0.2, 0.1, -0.3),(-2.,-2.,1.))), + "/home/connor/projects/thesis/LaTeX/fig/EMS_plot.png", + width=850, height=400) +savefig(Thesis.plot(short, title="EMS Mission Profile", mode="light", planet_colors=false, + phase_colors=["#F00","#F0F","#6AF"], markers=false, planets=false, + camera=((-0.2, 0.1, -0.3),(-2.,-2.,1.))), + "/home/connor/projects/thesis/LaTeX/fig/EMS_plot_noplanets.png", + width=850, height=400) +savefig(Thesis.plot(long, title="EMS Mission Profile", mode="light", planet_colors=false, + phase_colors=["#F00","#F0F","#6AF"], markers=false, + camera=((-0.2, -0.3, -0.2),(1.75,-1.75,1.5))), + "/home/connor/projects/thesis/LaTeX/fig/EMJS_plot.png", + width=850, height=400) +savefig(Thesis.plot(long, title="EMS Mission Profile", mode="light", planet_colors=false, + phase_colors=["#F00","#F0F","#6AF"], markers=false, planets=false, + camera=((-0.2, -0.3, -0.2),(2.,-2.,1.5))), + "/home/connor/projects/thesis/LaTeX/fig/EMJS_plot_noplanets.png", + width=850, height=400) +# Display Thrust Plots +# display(plot_thrust(display_mission, title="Thrust Magnitude vs Time", mode="light")) +# display(plot_thrust_components(display_mission, title="Thrust Components vs Time", mode="light")) + +# Save Thrust Plots +savefig(plot_thrust(short, title="Thrust Magnitude vs Time", mode="light"), + "/home/connor/projects/thesis/LaTeX/fig/EMS_thrust_mag.png", + width=850, height=400) +savefig(plot_thrust_components(short, title="Thrust Components vs Time", mode="light"), + "/home/connor/projects/thesis/LaTeX/fig/EMS_thrust_components.png", + width=850, height=400) +savefig(plot_thrust(long, title="Thrust Magnitude vs Time", mode="light"), + "/home/connor/projects/thesis/LaTeX/fig/EMJS_thrust_mag.png", + width=850, height=400) +savefig(plot_thrust_components(long, title="Thrust Components vs Time", mode="light"), + "/home/connor/projects/thesis/LaTeX/fig/EMJS_thrust_components.png", + width=850, height=400) + +# Save mission # store(display_mission, "/home/connor/projects/thesis/archive/best/long_mission") diff --git a/julia/porkchop.jl b/julia/porkchop.jl new file mode 100644 index 0000000..e9aafe2 --- /dev/null +++ b/julia/porkchop.jl @@ -0,0 +1,7 @@ +using Dates +using Thesis + +departures = collect(DateTime(2023,7,1):Dates.Day(1):DateTime(2023,12,1)) +arrivals = collect(DateTime(2025,6,1):Dates.Day(1):DateTime(2026,1,1)) +p = porkchop(Earth, Mars, departures, arrivals) +display(p); diff --git a/julia/src/types/mission.jl b/julia/src/types/mission.jl index 76dcd14..6c14b24 100644 --- a/julia/src/types/mission.jl +++ b/julia/src/types/mission.jl @@ -344,7 +344,7 @@ end """ Increases the fidelity of a mission by using it as a base and resolving with increased number of points """ -function improve(mission::Mission; n::Int=5) +function improve(mission::Mission; n::Int=5, use_cost=false) new_phases = Vector{Phase}() for phase in mission.phases new_thrusts = repeat(phase.thrust_profile, inner=(n,1)) @@ -366,7 +366,8 @@ function improve(mission::Mission; n::Int=5) arrival + Dates.Day(7), 200., 20., - CPU_multiplier = 30. + CPU_multiplier = 60., + use_cost=use_cost, ) new_mission.converged || error("Mission improvement failed") return new_mission diff --git a/julia/src/utilities/conversions.jl b/julia/src/utilities/conversions.jl index 6e7a3a9..170d3d2 100644 --- a/julia/src/utilities/conversions.jl +++ b/julia/src/utilities/conversions.jl @@ -1,4 +1,4 @@ -export oe_to_rθh, rθh_to_xyz, oe_to_xyz, xyz_to_oe, conv_T, spiral, gen_orbit +export oe_to_rθh, rθh_to_xyz, oe_to_xyz, xyz_to_oe, xyz_to_rθh, conv_T, spiral, gen_orbit function oe_to_rθh(oe::Vector,μ::Real) :: Vector @@ -25,6 +25,18 @@ function rθh_to_xyz(rθh_vec::Vector,oe::Vector) end +function xyz_to_rθh(xyz_vec::Vector,oe::Vector) + + a,e,i,Ω,ω,ν = oe + θ = ω+ν + cΩ,sΩ,ci,si,cθ,sθ = cos(Ω),sin(Ω),cos(i),sin(i),cos(θ),sin(θ) + DCM = [cΩ*cθ-sΩ*ci*sθ -cΩ*sθ-sΩ*ci*cθ sΩ*si; + sΩ*cθ+cΩ*ci*sθ -sΩ*sθ+cΩ*ci*cθ -cΩ*si; + si*sθ si*cθ ci] + return inv(DCM)*xyz_vec + +end + function oe_to_xyz(oe::Vector,μ::Real) return rθh_to_xyz(oe_to_rθh(oe,μ),oe) diff --git a/julia/src/utilities/lamberts.jl b/julia/src/utilities/lamberts.jl index dc2e2a9..67bb787 100644 --- a/julia/src/utilities/lamberts.jl +++ b/julia/src/utilities/lamberts.jl @@ -1,3 +1,5 @@ +export lamberts, porkchop + """ 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 @@ -68,3 +70,47 @@ function lamberts(planet1::Body,planet2::Body,leave::DateTime,arrive::DateTime) return v0t - state1[4:6], vft - state2[4:6], tof_req end + +function porkchop(planet1::Body, planet2::Body, departures::Vector{DateTime}, arrivals::Vector{DateTime}) + v∞_in = [ norm(lamberts(planet1, planet2, depart, arrive)[2]) for depart in departures, arrive in arrivals ] + v∞_out = [ norm(lamberts(planet1, planet2, depart, arrive)[1]) for depart in departures, arrive in arrivals ] + c3 = [ norm(lamberts(planet1, planet2, depart, arrive)[1])^2 for depart in departures, arrive in arrivals ] + traces = Vector{PlotlyJS.AbstractTrace}() + outputs = (v∞_in, v∞_out) + names = ("Earth Departure v∞", "Mars Arrival v∞") + colors = ("#0A0", "#A0A") + for (output, name, color) in zip(outputs, names, colors) + push!(traces, contour(x=departures, + y=arrivals, + z=output, + name=name, + contours_coloring="none", + contours_showlabels=true, + contours_start=12, + contours_end=20, + contours_size=2, + line_color=color, + line_width=2, + showscale=false, + )) + push!(traces, contour(x=departures, + y=arrivals, + z=output, + showlegend=false, + contours_coloring="none", + contours_showlabels=true, + contours_start=0, + contours_end=10, + contours_size=1, + line_color=color, + line_width=2, + showscale=false, + )) + end + layout = Layout(title_text="Sample Porkchop Earth->Mars Plot", + xaxis_title="Departure Date", + yaxis_title="Arrival Date", + margin_l=120, + ) + return PlotlyJS.plot(traces, layout) +end diff --git a/julia/src/utilities/plotting.jl b/julia/src/utilities/plotting.jl index 3bd20f2..3333c9d 100644 --- a/julia/src/utilities/plotting.jl +++ b/julia/src/utilities/plotting.jl @@ -309,7 +309,10 @@ function plot(m::Union{Mission, Mission_Guess}; mode::String="dark", planet_colors=true, phase_colors=nothing, - markers=true) + markers=true, + planets=true, + camera=nothing, + ) # First plot the earth # Then plot phase plus planet phase time = datetime2julian(m.launch_date) @@ -325,7 +328,11 @@ function plot(m::Union{Mission, Mission_Guess}; end end start = state(current_planet, time, m.launch_v∞, mass) - traces = [ earth_trace... ] + if planets + traces = [ earth_trace... ] + else + traces = Vector{PlotlyJS.AbstractTrace}() + end i = 1 for phase in m.phases @@ -365,8 +372,10 @@ function plot(m::Union{Mission, Mission_Guess}; color="#000") end end - push!(traces, planet_trace...) - limit = max(limit, new_limit) + if planets + push!(traces, planet_trace...) + limit = max(limit, new_limit) + end i += 1 end @@ -376,6 +385,14 @@ function plot(m::Union{Mission, Mission_Guess}; # Determine layout details layout = standard_layout(limit, title, mode=mode) + if camera !== nothing + layout["scene_camera_center_x"] = camera[1][1] + layout["scene_camera_center_y"] = camera[1][2] + layout["scene_camera_center_z"] = camera[1][3] + layout["scene_camera_eye_x"] = camera[2][1] + layout["scene_camera_eye_y"] = camera[2][2] + layout["scene_camera_eye_z"] = camera[2][3] + end # Plot PlotlyJS.plot( PlotlyJS.Plot( traces, layout ) ) @@ -410,30 +427,40 @@ function plot_thrust(m::Union{Mission, Mission_Guess}; title::String="Mission Plot", mode::String="dark", ) - times = Vector{Float64}() + times = Vector{String}() vals = Vector{Float64}() time = datetime2julian(m.launch_date) start = state(Earth, time, m.launch_v∞, m.start_mass) + traces = Vector{PlotlyJS.AbstractTrace}() + limits = Vector{Float64}() for phase in m.phases - # First get the path + n = size(phase.thrust_profile)[1] + for i in 1:n + push!(times, format(julian2datetime(time + i*phase.tof/(86400n)),"yyyy-mm-dd HH:MM:SS")) + push!(vals, norm(phase.thrust_profile[i,:])) + end path = prop(phase.thrust_profile, start, m.sc, phase.tof) time += phase.tof/86400. + time_string = format(julian2datetime(time),"yyyy-mm-dd HH:MM:SS") + push!(traces, scatter(;x=[time_string, time_string], + y=[-1.0,1.0], + mode="lines", + line_color="#000", + name=phase.planet.name*" Flyby")) mass = path[end,end] start = state(phase.planet, time, phase.v∞_out, mass) - n = size(phase.thrust_profile)[1] - println(n) - for i in 1:n - push!(times, time + i*phase.tof/(86400n)) - push!(vals, norm(phase.thrust_profile[i])) - end end - trace = scatter(;x=times, y=vals) + push!(traces, scatter(;x=times, y=vals, name="Thrust Mag")) # Determine layout details - layout = Layout(;title=title, xaxis_title="Time (JD)", yaxis_title="Thrust (% of Max)") + layout = Layout(;title=title, + xaxis_title="Year", + yaxis_title="Thrust (% of Max)", + yaxis_range=[0.0,maximum(vals)*1.02], + ) # Plot - PlotlyJS.plot( PlotlyJS.Plot( [trace], layout ) ) + PlotlyJS.plot( PlotlyJS.Plot( traces, layout ) ) end @@ -441,37 +468,50 @@ function plot_thrust_components(m::Union{Mission, Mission_Guess}; title::String="Mission Plot", mode::String="dark", ) - times = Vector{Float64}() - vals1 = Vector{Float64}() - vals2 = Vector{Float64}() - vals3 = Vector{Float64}() + traces = Vector{PlotlyJS.AbstractTrace}() + times = Vector{String}() + xs = Vector{Float64}() + ys = Vector{Float64}() + zs = Vector{Float64}() + αs = Vector{Float64}() + βs = Vector{Float64}() time = datetime2julian(m.launch_date) start = state(Earth, time, m.launch_v∞, m.start_mass) for phase in m.phases - # First get the path + n = size(phase.thrust_profile)[1] path = prop(phase.thrust_profile, start, m.sc, phase.tof) + for i in 1:n + push!(times, format(julian2datetime(time + i*phase.tof/(86400n)),"yyyy-mm-dd HH:MM:SS")) + push!(xs, phase.thrust_profile[i,1]) + push!(ys, phase.thrust_profile[i,2]) + push!(zs, phase.thrust_profile[i,3]) + r, θ, h = xyz_to_rθh(phase.thrust_profile[i,1:3], xyz_to_oe(path[i+1,1:6],Sun.μ)) + β = asin(h) + α = atan(r/cos(β),θ/cos(β)) + push!(αs, α) + push!(βs, β) + end time += phase.tof/86400. + time_string = format(julian2datetime(time),"yyyy-mm-dd HH:MM:SS") + push!(traces, scatter(;x=[time_string, time_string], + y=[-10.0,10.0], + mode="lines", + line_color="#000", + name=phase.planet.name*" Flyby")) mass = path[end,end] start = state(phase.planet, time, phase.v∞_out, mass) - n = size(phase.thrust_profile)[1] - println(n) - for i in 1:n - push!(times, time + i*phase.tof/(86400n)) - push!(vals1, phase.thrust_profile[i,1]) - push!(vals2, phase.thrust_profile[i,2]) - push!(vals3, phase.thrust_profile[i,3]) - end end - trace1 = scatter(;x=times, y=vals1, name="x component") - trace2 = scatter(;x=times, y=vals2, name="y component") - trace3 = scatter(;x=times, y=vals3, name="z component") + push!(traces, scatter(;x=times, y=αs, name="in-plane angle")) + push!(traces, scatter(;x=times, y=βs, name="out-of-plane angle")) # Determine layout details layout = Layout(;title=title, - xaxis_title="Time (JD)", - yaxis_title="Thrust (% of Max)") + xaxis_title="Year", + yaxis_title="Thrust Angle (Radians)", + yaxis_range=[min(minimum(αs),minimum(βs))*1.02,max(maximum(αs),maximum(βs))*1.02], + ) # Plot - PlotlyJS.plot( PlotlyJS.Plot( [trace1, trace2, trace3], layout ) ) + PlotlyJS.plot( PlotlyJS.Plot( traces, layout ) ) end diff --git a/notes.md b/notes.md index 85c6bb9..9c3d630 100644 --- a/notes.md +++ b/notes.md @@ -64,8 +64,9 @@ Todo: - [x] Include how mass is being decremented in the SFT propagator section - [ ] Think about how to include plots that show the MBH algorithm jumping through the space and looking at various basins. +- [ ] Porkchop plot - [x] Include some plots of the thrust vectors rather than the raw data in the appendix - - [ ] Format dates better and include lines where the flybys are - - [ ] Reconsider wording on analysis - - [ ] Reconsider plotting the angles instead + - [x] Format dates better and include lines where the flybys are + - [x] Reconsider wording on analysis + - [x] Reconsider plotting the angles instead - [ ] Submit thesis defense form