Update for ky

julia/src/utilities/constants.jl

@@ -0,0 +1,29 @@
+# -----------------------------------------------------------------------------
+# DEFINING CONSTANTS
+# -----------------------------------------------------------------------------
+
+export Body, Sun, Mercury, Venus, Earth, Moon, Mars
+export Jupiter, Saturn, Uranus, Neptune, Pluto
+export G, AU, init_STM, hour, day, year, second
+export Pathlist
+
+const Sun = Body(1.32712440018e11, 696000., "YlOrRd", 10, 0., "", "Sun")
+const Mercury = Body(2.2032e4, 2439., "YlOrRd", 1, 57909083., "#FF2299", "Mercury")
+const Venus = Body(3.257e5, 6052., "Portland", 2, 108208601., "#FFCC00", "Venus")
+const Earth = Body(3.986004415e5, 6378.1363, "Earth", 399, 149598023., "#0055FF", "Earth")
+const Mars = Body(4.305e4, 3397.2, "Reds", 4, 227939186., "#FF0000", "Mars")
+const Jupiter = Body(1.266865361e8, 71492., "YlOrRd", 5, 778298361., "#FF7700", "Jupiter")
+const Saturn = Body(3.794e7, 60268., "Portland", 6, 1429394133., "#FFFF00", "Saturn")
+const Uranus = Body(5.794e6, 25559., "YlGnBu", 7, 2875038615., "#0000FF", "Uranus")
+const Neptune = Body(6.809e6, 24764., "Blues", 8, 4504449769., "#6666FF", "Neptune")
+const Pluto = Body(9e2, 1151., "Picnic", 9, 5915799000., "#FFCCCC", "Pluto")
+
+const G = 6.67430e-20 #universal gravity parameter
+const AU = 149597870.691 #km
+const init_STM = vec(Matrix{Float64}(I,6,6))
+const second = 1.
+const hour = 3600.
+const day = 86400.
+const year = 365 * day
+
+Pathlist = Vector{Vector{Vector{Float64}}}

julia/src/utilities/conversions.jl

@@ -0,0 +1,141 @@
+export oe_to_rθh, rθh_to_xyz, oe_to_xyz, xyz_to_oe, conv_T, spiral, gen_orbit
+
+function oe_to_rθh(oe::Vector,μ::Real) :: Vector
+
+  a,e,i,Ω,ω,ν = oe
+  return [a*(1-e^2)/(1+e*cos(ν)),
+  0,
+  0,
+  (μ/sqrt(μ*a*(1-e^2)))*e*sin(ν),
+  (μ/sqrt(μ*a*(1-e^2)))*(1+e*cos(ν)),
+  0]
+
+end
+
+function rθh_to_xyz(rθh_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]
+  DCM = kron(Matrix(I,2,2),DCM)
+  return DCM*rθh_vec
+
+end
+
+function oe_to_xyz(oe::Vector,μ::Real)
+
+  return rθh_to_xyz(oe_to_rθh(oe,μ),oe)
+
+end
+
+function xyz_to_oe(cart_vec::Vector,μ::Real)
+
+  r_xyz, v_xyz = cart_vec[1:3],cart_vec[4:6]
+  r = norm(r_xyz)
+  h_xyz = cross(r_xyz,v_xyz) #km^2/s
+  h = norm(h_xyz) #km^2/s
+  ξ = dot(v_xyz,v_xyz)/2 - μ/r #km^2 s^-2
+  a = -μ/(2ξ) #km
+  e_xyz = cross(v_xyz,h_xyz)/μ - r_xyz/r
+  e = norm(e_xyz)
+  if h_xyz[3]/h < 1.
+    i = acos(h_xyz[3]/h) #rad
+  else
+    i = acos(1.)
+  end
+  n_xyz = cross([0,0,1],h_xyz)
+  if dot(n_xyz,[1,0,0])/norm(n_xyz) < 1.
+    Ω = acos(dot(n_xyz,[1,0,0])/norm(n_xyz))
+  else
+    Ω = acos(1.)
+  end
+  if dot(n_xyz,e_xyz)/(norm(n_xyz)*e) < 1.
+    ω = acos(dot(n_xyz,e_xyz)/(norm(n_xyz)*e))
+  else
+    ω = acos(1.)
+  end
+  if abs((dot(r_xyz,e_xyz))/(r*norm(e_xyz))) < 1.
+    ν = acos((dot(r_xyz,e_xyz))/(r*norm(e_xyz)))
+  else
+    ν = acos(1.)
+  end
+  Ω = dot(n_xyz,[0,1,0]) > 0. ? Ω : -Ω
+  ω = dot(e_xyz,[0,0,1]) > 0. ? ω : -ω
+  ν = dot(r_xyz,v_xyz) > 0. ? ν : -ν
+  return [a,e,i,Ω,ω,ν]
+
+end
+
+"""
+A convenience function for generating start conditions from orbital elements
+Inputs: a body, a period, and a mass
+Output: a random reasonable orbit
+"""
+function gen_orbit(T::Float64, mass::Float64, primary::Body=Sun)
+  μ = primary.μ
+  i = rand(0.0:0.01:0.4999π)
+  θ = rand(0.0:0.01:2π)
+  i = 0
+  while true
+    i += 1
+    e = rand(0.0:0.01:0.5)
+    a = ∛(μ * ( T/2π )^2 )
+    a*(1-e) < 1.1primary.r || return [ oe_to_xyz([ a, e, i, 0., 0., θ ], μ); mass ]
+    i < 100 || throw(GenOrbit_Error)
+  end
+end
+
+"""
+A convenience function for generating spiral trajectories
+"""
+spiral(mag,n,init,sc,time,primary=Sun) = conv_T(fill(mag, n), zeros(n), zeros(n), init, sc, time, primary)
+
+"""
+Converts a series of thrust vectors from R,Θ,H frame to cartesian
+"""
+function conv_T(Tm::Vector{Float64},
+                Ta::Vector{Float64},
+                Tb::Vector{Float64},
+                init_state::Vector{Float64},
+                craft::Sc,
+                time::Float64,
+                primary::Body=Sun)
+
+  Txs = Float64[]
+  Tys = Float64[]
+  Tzs = Float64[]
+  n::Int = length(Tm)
+
+  for i in 1:n
+    mag, α, β = Tm[i], Ta[i], Tb[i]
+    if mag > 1 || mag < 0
+      @error "ΔV input is too high: $mag"
+    elseif α > π || α < -π
+      @error "α angle is incorrect: $α"
+    elseif β > π/2 || β < -π/2
+      @error "β angle is incorrect: $β"
+    end
+  end
+
+  state = init_state
+  for i in 1:n
+    mag, α, β = Tm[i], Ta[i], Tb[i]
+    thrust_rθh = mag * [cos(β)*sin(α), cos(β)*cos(α), sin(β)]
+    _,_,i,Ω,ω,ν = xyz_to_oe(state, primary.μ)
+    θ = ω+ν
+    cΩ,sΩ,ci,si,cθ,sθ = cos(Ω),sin(Ω),cos(i),sin(i),cos(θ),sin(θ)
+    DCM = [cΩ*cθ-sΩ*ci*sθ -cΩ*sθ-sΩ*ci*cθ sΩ*si;
+           sΩ*cθ+cΩ*ci*sθ -sΩ*sθ+cΩ*ci*cθ -cΩ*si;
+           si*sθ          si*cθ           ci      ]
+    Tx, Ty, Tz = DCM*thrust_rθh
+
+    state = prop_one([Tx, Ty, Tz], state, craft, time/n, primary)
+    push!(Txs, Tx)
+    push!(Tys, Ty)
+    push!(Tzs, Tz)
+  end
+  return hcat(Txs, Tys, Tzs)
+end

julia/src/utilities/laguerre-conway.jl

@@ -0,0 +1,57 @@
+function laguerre_conway(state::Vector{<:Real}, time::Float64, primary::Body=Sun)
+
+  μ = primary.μ
+  n = 5                               # Choose LaGuerre-Conway "n"
+  i = 0
+
+  r0 = state[1:3]                     # Are we in elliptical or hyperbolic orbit?
+  r0_mag = √(state[1]^2 + state[2]^2 + state[3]^2)
+  v0 = state[4:6]
+  v0_mag = √(state[4]^2 + state[5]^2 + state[6]^2)
+  σ0 = (r0 ⋅ v0)/√(μ)
+  a = 1 / ( 2/r0_mag - v0_mag^2/μ )
+  coeff = 1 - r0_mag/a
+
+  if a > 0    # Elliptical
+    ΔM = ΔE_new = √(μ) / sqrt(a^3) * time
+    ΔE = 1000
+    while abs(ΔE - ΔE_new) > 1e-10
+      ΔE = ΔE_new
+      F = ΔE - ΔM + σ0 / √(a) * (1-cos(ΔE)) - coeff * sin(ΔE)
+      dF = 1 + σ0 / √(a) * sin(ΔE) - coeff * cos(ΔE)
+      d2F = σ0 / √(a) * cos(ΔE) + coeff * sin(ΔE)
+      sign = dF >= 0 ? 1 : -1
+      ΔE_new = ΔE - n*F / ( dF + sign * √(abs((n-1)^2*dF^2 - n*(n-1)*F*d2F )))
+      i += 1
+      if i > 100 throw(LaGuerreConway_Error()) end
+    end
+    F = 1 - a/r0_mag * (1-cos(ΔE))
+    G = a * σ0/ √(μ) * (1-cos(ΔE)) + r0_mag * √(a) / √(μ) * sin(ΔE)
+    r = a + (r0_mag - a) * cos(ΔE) + σ0 * √(a) * sin(ΔE)
+    Ft = -√(a)*√(μ) / (r*r0_mag) * sin(ΔE)
+    Gt = 1 - a/r * (1-cos(ΔE))
+  else      # Hyperbolic or Parabolic
+    ΔN = √(μ) / sqrt(-a^3) * time
+    ΔH = 0
+    ΔH_new = time < 0 ? -1 : 1
+    while abs(ΔH - ΔH_new) > 1e-10
+      ΔH = ΔH_new
+      F = -ΔN - ΔH + σ0 / √(-a) * (cos(ΔH)-1) + coeff * sin(ΔH)
+      dF = -1 + σ0 / √(-a) * sin(ΔH) + coeff * cos(ΔH)
+      d2F = σ0 / √(-a) * cos(ΔH) + coeff * sin(ΔH)
+      sign = dF >= 0 ? 1 : -1
+      ΔH_new = ΔH - n*F / ( dF + sign * √(abs((n-1)^2*dF^2 - n*(n-1)*F*d2F )))
+      i += 1
+      if i > 100 throw(LaGuerreConway_Error()) end
+    end
+    F = 1 - a/r0_mag * (1-cos(ΔH))
+    G = a * σ0/ √(μ) * (1-cos(ΔH)) + r0_mag * √(-a) / √(μ) * sin(ΔH)
+    r = a + (r0_mag - a) * cos(ΔH) + σ0 * √(-a) * sin(ΔH)
+    Ft = -√(-a)*√(μ) / (r*r0_mag) * sin(ΔH)
+    Gt = 1 - a/r * (1-cos(ΔH))
+  end
+
+
+  return [ F*state[1:3] + G*state[4:6]; Ft*state[1:3] + Gt*state[4:6]]
+
+end

julia/src/utilities/lamberts.jl

@@ -0,0 +1,70 @@
+"""
+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(format(leave,"yyyy-mm-ddTHH:MM:SS"))
+  time_arrive = utc2et(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

julia/src/utilities/plotting.jl

@@ -0,0 +1,292 @@
+# ------------------------------------------------------------------------------
+# PLOTTING FUNCTIONS
+
+using Random, PlotlyJS, Base.Iterators
+
+export plot
+
+function random_color()
+  num1 = rand(0:255)
+  num2 = rand(0:255)
+  num3 = rand(0:255)
+  return "#"*string(num1, base=16, pad=2)*string(num2, base=16, pad=2)*string(num3, base=16, pad=2)
+end
+
+# Some convenience functions for maximums
+"""
+Gets the maximum absolute value of a path
+"""
+function Base.maximum(path::Vector{Vector{Float64}})
+  return maximum( abs.( reduce( vcat,path ) ) )
+end
+
+"""
+Gets the maximum absolute value of a list of paths
+"""
+function Base.maximum(paths::Vector{Vector{Vector{Float64}}})
+  return maximum( abs.( reduce( vcat,reduce( vcat,paths ) ) ) )
+end
+
+"""
+Gets the maximum absolute value of a surface
+"""
+function Base.maximum(xs::Matrix{Float64},ys::Matrix{Float64},zs::Matrix{Float64})
+  return maximum([ vec(xs), vec(ys), vec(zs) ])
+end
+
+"""
+Basically rewriting Julia at this point
+"""
+Base.getindex(x::Nothing, i::Int) = nothing
+
+"""
+Generates a layout that works for most plots. Requires a title and a limit.
+"""
+function standard_layout(limit::Float64, title::AbstractString)
+  limit *= 1.1
+  Layout(title=title,
+         # width=1400,
+         # height=800,
+         paper_bgcolor="rgba(5,10,40,1.0)",
+         plot_bgcolor="rgba(100,100,100,0.01)",
+         scene = attr(xaxis = attr(autorange = false,range=[-limit,limit]),
+                      yaxis = attr(autorange = false,range=[-limit,limit]),
+                      zaxis = attr(autorange = false,range=[-limit,limit]),
+         aspectratio=attr(x=1,y=1,z=1),
+         aspectmode="manual"))
+end
+
+"""
+A fundamental function to generate a plot for a single orbit path
+"""
+function gen_plot(path::Vector{Vector{Float64}};
+                  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]),
+                            mode="lines", showlegend=false, line_color=color, line_width=3)
+  else
+    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]]),
+                             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]]),
+                             mode="markers", showlegend=false,
+                             marker=attr(color=color, size=3, symbol="diamond")))
+  end
+
+  return traces, maximum(path[1:6])
+
+end
+
+"""
+A fundamental function to plot a single orbit path
+"""
+function plot(path::Vector{Vector{Float64}};
+              color::AbstractString = random_color(),
+              title::AbstractString = "Single Orbit Path",
+              markers=true)
+
+  # Generate the plot
+  p, limit = gen_plot(path, color=color, markers=markers)
+
+  # Now use the standard layout
+  layout = standard_layout(limit, title)
+
+  # Generate the plot
+  PlotlyJS.plot( PlotlyJS.Plot( p, layout ) )
+
+end
+
+"""
+This will generate the plot of a planetary body in it's correct colors and position in a
+sun-centered frame at a certain time
+
+If the body is the sun, no time is necessary
+
+An optional prop time is allowed, which will also provide an orbit
+"""
+function gen_plot(b::Body,
+                  time::Union{DateTime, Nothing}=nothing,
+                  prop_time::Union{Float64, Nothing}=nothing;
+                  scale::Int=100)
+
+  # Gotta be a little weird with Sun Scale
+  if b == Sun && scale > 25 scale = 25 end
+
+  # This is how many points to poll to generate the sphere
+  N = 32
+  θ = collect(range(0,length=N,stop=2π))
+  ϕ = collect(range(0,length=N,stop=π))
+
+  # Convert polar -> cartesian
+  unscaled_xs = cos.(θ) * sin.(ϕ)'
+  unscaled_ys = sin.(θ) * sin.(ϕ)'
+  unscaled_zs = repeat(cos.(ϕ)',outer=[N, 1])
+
+  # Scale by the size of the body (and a scaling factor for visualization)
+  xs,ys,zs = scale * b.r .* (unscaled_xs,unscaled_ys,unscaled_zs)
+
+  # Add an offset if we're not plotting the sun
+  body_state = zeros(6)
+  if b != Sun body_state = state(b, time) end
+
+  max_value = maximum(xs .+ body_state[1], ys .+ body_state[2], zs .+ body_state[3])
+
+  # Generate the surface trace
+  body_surface = surface(;x = xs .+ body_state[1],
+                          y = ys .+ body_state[2],
+                          z = zs .+ body_state[3],
+                          showscale = false,
+                          colorscale = b.color)
+  outputs = [ body_surface ]
+
+  if prop_time !== nothing
+    p, maybe_max = gen_plot(prop(body_state, prop_time), label=b.name, color=b.line_color, markers=false)
+    push!(outputs, p...)
+    if maybe_max > max_value max_value = maybe_max end
+  end
+
+  return outputs, max_value
+
+end
+
+"""
+This will plot a planetary body in it's correct colors and position in a
+sun-centered frame at a certain time
+
+If the body is the sun, no time is necessary
+
+An optional prop time is allowed, which will also provide an orbit
+"""
+function plot(b::Body,
+              time::Union{DateTime, Nothing}=nothing,
+              prop_time::Union{Float64, Nothing}=nothing;
+              title::AbstractString="Single Planet Plot",
+              scale::Int=100)
+
+  # Generate the plot
+  p, limit = gen_plot(b,time,prop_time,scale=scale)
+
+  # Now use the standard layout
+  layout = standard_layout(limit, title)
+
+  # Generate the plot
+  PlotlyJS.plot( PlotlyJS.Plot( p, layout ) )
+
+end
+
+"""
+This will plot a vector of planetary bodies at a certain time
+
+The Sun is plotted automatically
+
+An optional orbits bool is allowed, which shows the orbits
+"""
+function plot(bodies::Vector{Body},
+              time::Union{DateTime, Nothing};
+              orbits::Bool=true,
+              title::AbstractString="Multiple Planet Plot",
+              scale::Int=100)
+
+  # Generate the plots
+  trace, limit = gen_plot(Sun, scale=scale)
+  traces = [ trace... ]
+  for b in bodies
+    if orbits
+      trace, new_limit = gen_plot(b,time,period(b),scale=scale)
+    else
+      trace, new_limit = gen_plot(b,time,scale=scale)
+    end
+    push!(traces, trace...)
+    limit = max(limit, new_limit)
+  end
+
+  # Now use the standard layout
+  layout = standard_layout(limit, title)
+
+  # Generate the plot
+  PlotlyJS.plot( PlotlyJS.Plot( traces, layout ) )
+
+end
+
+"""
+Generic plotting function for a bunch of paths. Useful for debugging.
+"""
+function plot(paths::Vector{Vector{Vector{Float64}}};
+              labels::Union{Vector{String},Nothing}=nothing,
+              colors::Union{Vector{String},Nothing}=nothing,
+              markers::Bool=true,
+              title::String="Spacecraft Position")
+
+  # Generate the path traces
+  color = colors === nothing ? random_color() : colors[1]
+  trace, limit = gen_plot(paths[1],label=labels[1],color=color,markers=markers)
+  traces = [ trace... ]
+  for i = 2:length(paths)
+    color = colors === nothing ? random_color() : colors[i]
+    trace, new_limit = gen_plot(paths[i],label=labels[i],color=color,markers=markers)
+    push!(traces, trace...)
+    limit = max(limit, new_limit)
+  end
+
+  # Generate the sun trace
+  push!(traces, gen_plot(Sun)[1]...)
+
+  # Determine layout details
+  layout = standard_layout(limit, title)
+
+  # Plot
+  PlotlyJS.plot( PlotlyJS.Plot( traces, layout ) )
+
+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
+  mass = m.start_mass
+  current_planet = Earth
+  earth_trace, limit = gen_plot(current_planet, 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)
+    mass = final[7]
+    time += Second(floor(phase.tof))
+    current_planet = phase.planet
+    start = state(current_planet, time, phase.v∞_out, mass)
+    path_trace, new_limit = gen_plot(path, label="Phase "*string(i))
+    push!(traces, path_trace...)
+    limit = max(limit, new_limit)
+
+    # Then the planet
+    planet_trace, new_limit = gen_plot(current_planet, time, -phase.tof)
+    push!(traces, planet_trace...)
+    limit = max(limit, new_limit)
+
+    i += 1
+  end
+  
+  # Generate the sun trace
+  push!(traces, gen_plot(Sun)[1]...)
+
+  # Determine layout details
+  layout = standard_layout(limit, title)
+
+  # Plot
+  PlotlyJS.plot( PlotlyJS.Plot( traces, layout ) )
+
+end

julia/src/utilities/propagator.jl

@@ -0,0 +1,103 @@
+export prop
+
+"""
+Maximum ΔV that a spacecraft can impulse for a given single time step
+"""
+function max_ΔV(duty_cycle::Float64,
+                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]
+
+end
+
+"""
+The propagator function
+"""
+function prop(ΔVs::Matrix{T},
+              state::Vector{Float64},
+              craft::Sc,
+              time::Float64,
+              primary::Body=Sun;
+              interpolate::Bool=false) where T <: Real
+
+  size(ΔVs)[2] == 3 || throw(ΔVsize_Error())
+  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
+
+  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
+    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
+  end
+
+  state[7] > craft.dry_mass || throw(Mass_Error(state[7]))
+
+  return states, state
+
+end
+
+"""
+Convenience function for propagating a state with no thrust
+"""
+prop(x::Vector{Float64}, t::Float64, p::Body=Sun) = prop(zeros(1000,3), [x;1.], no_thrust, t, p)[1]
+
+"""
+This is solely for the purposes of getting the final state of a mission or guess
+"""
+function prop(m::Union{Mission, Mission_Guess})
+    time = m.launch_date
+    current_planet = Earth
+    start = state(current_planet, time, m.launch_v∞, m.start_mass)
+
+    final = zeros(7)
+    for phase in m.phases
+      final = prop(phase.thrust_profile, start, m.sc, phase.tof)[2]
+      mass = final[7]
+      current_planet = phase.planet
+      time += Second(floor(phase.tof))
+      start = state(current_planet, time, phase.v∞_out, mass)
+    end
+
+    return final
+
+end