Re-packaged

julia/src/Thesis.jl

@@ -0,0 +1,12 @@
+module Thesis
+
+  using LinearAlgebra, ForwardDiff, Blink, PlotlyJS
+
+  include("./constants.jl")
+  include("./conversions.jl")
+  include("./spacecraft.jl")
+  include("./plotting.jl")
+  include("./laguerre-conway.jl")
+  include("./propagator.jl")
+  include("./find_closest.jl")
+end

julia/src/constants.jl

@@ -0,0 +1,100 @@
+# -----------------------------------------------------------------------------
+# DEFINING CONSTANTS
+# -----------------------------------------------------------------------------
+
+export μs, G, GMs, μ, rs, as, es, AU
+
+# Gravitational Constants
+  μs = Dict(
+  "Sun" => 1.32712440018e11,
+  "Mercury" => 2.2032e4,
+  "Venus" => 3.257e5,
+  "Earth" => 3.986004415e5,
+  "Moon" => 4.902799e3,
+  "Mars" => 4.305e4,
+  "Jupiter" => 1.266865361e8,
+  "Saturn" => 3.794e7,
+  "Uranus" => 5.794e6,
+  "Neptune" => 6.809e6,
+  "Pluto" => 9e2)
+
+  G = 6.67430e-20
+
+  function μ(m1::Float64, m2::Float64)
+    return m2/(m1+m2)
+  end
+
+  function μ(GM1::Float64, GM2::Float64, Grav::Float64)
+    return μ(GM1/Grav, GM2/Grav)
+  end
+
+  function μ(primary::String, secondary::String)
+    return μ(GMs[primary]/G, GMs[secondary]/G)
+  end
+
+  GMs = Dict(
+  "Sun" => 132712440041.93938,
+  "Earth" => 398600.435436,
+  "Moon" => 4902.800066)
+
+# Radii
+  rs = Dict(
+  "Sun" => 696000.,
+  "Mercury" => 2439.,
+  "Venus" => 6052.,
+  "Earth" => 6378.1363,
+  "Moon" => 1738.,
+  "Mars" => 3397.2,
+  "Jupiter" => 71492.,
+  "Saturn" => 60268.,
+  "Uranus" => 25559.,
+  "Neptune" => 24764.,
+  "Pluto" => 1151.)
+
+# Semi Major Axes
+  as = Dict(
+  "Mercury" => 57909083.,
+  "Venus" => 108208601.,
+  "Earth" => 149598023.,
+  "Moon" => 384400.,
+  "Mars" => 227939186.,
+  "Jupiter" => 778298361.,
+  "Saturn" => 1429394133.,
+  "Uranus" => 2875038615.,
+  "Neptune" => 4504449769.,
+  "Pluto" => 5915799000.)
+
+# Eccentricities
+  es = Dict(
+  "Earth" => 0.016708617,
+  "Moon" => 0.0549)
+
+# J2 for basic oblateness
+  j2s = Dict(
+  "Mercury" => 0.00006,
+  "Venus" => 0.000027,
+  "Earth" => 0.0010826269,
+  "Moon" => 0.0002027,
+  "Mars" => 0.001964,
+  "Jupiter" => 0.01475,
+  "Saturn" => 0.01645,
+  "Uranus" => 0.012,
+  "Neptune" => 0.004,
+  "Pluto" => 0.)
+
+# These are just the colors for plots
+  p_colors = Dict(
+  "Sun" => :Electric,
+  "Mercury" => :heat,
+  "Venus" => :turbid,
+  "Earth" => :Blues,
+  "Moon" => :Greys,
+  "Mars" => :Reds,
+  "Jupiter" => :solar,
+  "Saturn" => :turbid,
+  "Uranus" => :haline,
+  "Neptune" => :ice,
+  "Pluto" => :matter)
+
+  AU = 149597870.691 #km
+  init_STM = vec(Matrix{Float64}(I,6,6))

julia/src/conversions.jl

@@ -0,0 +1,70 @@
+export oe_to_rθh, rθh_to_xyz, oe_to_xyz, xyz_to_oe
+
+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

julia/src/find_closest.jl

@@ -0,0 +1,24 @@
+using NLsolve
+
+function treat_inputs(x::AbstractVector)
+  n::Int = length(x)/3
+  reshape(x,(3,n))'
+end
+
+function single_shoot(start::Vector{Float64},
+                      final::Vector{Float64},
+                      craft::Sc,
+                      μ::Float64,
+                      t0::Float64,
+                      tf::Float64,
+                      x0::AbstractVector,
+                      tol=1e-2)
+
+  function f!(F,x)
+    F[1:6] .= prop(treat_inputs(x), start, craft, μ, tf-t0)[1][end,:] - final
+    F[7:3n] .= 0.
+  end
+
+  return nlsolve(f!, x0, ftol=tol, autodiff=:forward, iterations=10_000)
+
+end

julia/src/laguerre-conway.jl

@@ -0,0 +1,57 @@
+function laguerre_conway(state::Vector{T}, μ::Float64, time::Float64) where T <: Real
+
+  n = 5                               # Choose LaGuerre-Conway "n"
+  i = 0
+
+  r0 = state[1:3]                     # Are we in elliptical or hyperbolic orbit?
+  r0_mag = norm(r0)
+  v0 = state[4:6]
+  v0_mag = norm(v0)
+  σ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
+    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
+    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
+
+  if i > 100
+    throw(ErrorException("LaGuerre-Conway did not converge!"))
+  end
+
+  return [ F*state[1:3] + G*state[4:6]; Ft*state[1:3] + Gt*state[4:6]]
+
+end

julia/src/plotting.jl

@@ -0,0 +1,64 @@
+# ------------------------------------------------------------------------------
+# PLOTTING FUNCTIONS
+
+using Random, PlotlyJS, Base.Iterators
+
+export plot_orbits
+
+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
+
+function get_true_max(mat::Vector{Array{Float64,2}})
+  return maximum(abs.(flatten(mat)))
+end
+
+function plot_orbits(paths::Vector{Array{Float64,2}};
+                     primary::String="Earth",
+                     plot_theme::Symbol=:juno,
+                     labels::Vector{String}=Vector{String}(),
+                     title::String="Spacecraft Position",
+                     colors::Vector{String}=Vector{String}())
+
+  N = 32
+  θ = collect(range(0,length=N,stop=2π))
+  ϕ = collect(range(0,length=N,stop=π))
+  x = cos.(θ) * sin.(ϕ)'
+  y = sin.(θ) * sin.(ϕ)'
+  z = repeat(cos.(ϕ)',outer=[N, 1])
+  ps = rs[primary] .* (x,y,z)
+  x_p,y_p,z_p = ps
+
+  t1 = []
+  for i = 1:length(paths)
+    path = [ x for x in paths[i] ]
+    label = labels != [] ? labels[i] : "orbit"
+    color = colors != [] ? colors[i] : random_color()
+    push!(t1, scatter3d(;x=(path[:,1]),y=(path[:,2]),z=(path[:,3]),
+                         mode="lines", name=label, line_color=color, line_width=3))
+  end
+  limit = max(maximum(abs.(flatten(paths))),maximum(abs.(flatten(ps)))) * 1.1
+
+  t2 = surface(;x=(x_p),
+                y=(y_p),
+                z=(z_p),
+                showscale=false,
+                colorscale = p_colors[primary])
+
+  layout = Layout(;title="Orbit Plot",
+                   width=1000,
+                   height=600,
+                   paper_bgcolor="#222529",
+                   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"))
+
+  p = Plot([t1...,t2],layout)
+  plot(p)
+
+end

julia/src/propagator.jl

@@ -0,0 +1,115 @@
+export prop
+
+"""
+Maximum ΔV that a spacecraft can impulse for a given single time step
+"""
+function max_ΔV(duty_cycle::T,
+                num_thrusters::Int,
+                max_thrust::T,
+                tf::T,
+                t0::T,
+                mass::S) where {T <: Real, S <: Real}
+  return duty_cycle*num_thrusters*max_thrust*(tf-t0)/mass
+end
+
+"""
+This function propagates the spacecraft forward in time 1 Sim-Flanagan step (of variable length of time),
+applying a thrust in the center.
+"""
+function prop_one(ΔV::Vector{T},
+                  state::Vector{S},
+                  duty_cycle::Float64,
+                  num_thrusters::Int,
+                  max_thrust::Float64,
+                  mass::S,
+                  mass_flow_rate::Float64,
+                  μ::Float64,
+                  time::Float64) where {T <: Real, S <: Real}
+
+  mag, α, β = ΔV
+
+  # if mag > 1 || mag < 0
+    # throw(ErrorException("ΔV input is too high: $mag"))
+  # elseif α > π || α < -π
+    # throw(ErrorException("α angle is incorrect: $α"))
+  # elseif β > π/2 || β < -π/2
+    # throw(ErrorException("β angle is incorrect: $β"))
+  # end
+
+  thrust_rθh = mag * [cos(β)*sin(α), cos(β)*cos(α), sin(β)]
+  a,e,i,Ω,ω,ν = xyz_to_oe(state, μ)
+  θ = ω+ν
+  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]
+  ΔV = DCM*thrust_rθh
+
+  thrust = max_ΔV(duty_cycle, num_thrusters, max_thrust, time, 0., mass) * ΔV
+  halfway = laguerre_conway(state, μ, time/2) + [0., 0., 0., thrust[1], thrust[2], thrust[3]]
+  return laguerre_conway(halfway, μ, time/2), mass - mass_flow_rate*norm(ΔV)*time
+
+end
+
+"""
+A convenience function for using spacecraft. Note that this function outputs a sc instead of a mass
+"""
+function prop_one(ΔV_unit::Vector{T},
+                  state::Vector{S},
+                  craft::Sc,
+                  μ::Float64,
+                  time::Float64) where {T <: Real,S <: Real}
+  state, mass =  prop_one(ΔV_unit, state, craft.duty_cycle, craft.num_thrusters, craft.max_thrust,
+                          craft.mass, craft.mass_flow_rate, μ, time)
+  return state, Sc(mass, craft.mass_flow_rate, craft.max_thrust, craft.num_thrusters, craft.duty_cycle)
+end
+
+"""
+This propagates over a given time period, with a certain number of intermediate steps
+"""
+function prop(ΔVs::Matrix{T},
+              state::Vector{Float64},
+              duty_cycle::Float64,
+              num_thrusters::Int,
+              max_thrust::Float64,
+              mass::Float64,
+              mass_flow_rate::Float64,
+              μ::Float64,
+              time::Float64) where T <: Real
+
+  if size(ΔVs)[2] != 3 throw(ErrorException("ΔV input is wrong size")) end
+  n = size(ΔVs)[i]
+
+  for i in 1:n
+    state, mass = prop_one(ΔVs[i,:], state, duty_cycle, num_thrusters, max_thrust, mass,
+                           mass_flow_rate, μ, time/n)
+  end
+
+  return state, mass
+
+end
+
+"""
+The same function, using Scs
+"""
+function prop(ΔVs::AbstractArray{T},
+              state::Vector{Float64},
+              craft::Sc,
+              μ::Float64,
+              time::Float64) where T <: Real
+
+  if size(ΔVs)[2] != 3 throw(ErrorException("ΔV input is wrong size")) end
+  n = size(ΔVs)[1]
+
+  states = state'
+  masses = craft.mass
+
+  for i in 1:n
+    state, craft = prop_one(ΔVs[i,:], state, craft, μ, time/n)
+    states = [states; state']
+    masses = [masses, craft.mass]
+  end
+
+  return states, masses
+
+end

julia/src/spacecraft.jl

@@ -0,0 +1,19 @@
+export Sc
+
+struct Sc{T <: Real}
+  mass::T
+  mass_flow_rate::Float64
+  max_thrust::Float64
+  num_thrusters::Int
+  duty_cycle::Float64
+end
+
+function Sc(name::String)
+  if name == "test"
+    return Sc(10000., 0.01, 0.05, 2, 1.)
+  elseif name == "no_thrust"
+    return Sc(10000., 0.01, 0., 0, 0.)
+  else
+    throw(ErrorException("Bad sc name"))
+  end
+end