Began!

julia/constants.jl

@@ -0,0 +1,100 @@
+# -----------------------------------------------------------------------------
+# DEFINING CONSTANTS
+# -----------------------------------------------------------------------------
+
+export μs, G, GMs, μ, rs, as, es
+
+# 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/conversions.jl

@@ -0,0 +1,76 @@
+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
+  elseif h_xyz[3]/h < 1.000001
+    i = acos(1.)
+  else
+    error("bad i")
+  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))
+  elseif dot(n_xyz,[1,0,0])/norm(n_xyz) < 1.0001
+    Ω = acos(1.)
+  else
+    error("bad Ω")
+  end
+  if dot(n_xyz,e_xyz)/(norm(n_xyz)*e) < 1.
+    ω = acos(dot(n_xyz,e_xyz)/(norm(n_xyz)*e))
+  elseif dot(n_xyz,e_xyz)/(norm(n_xyz)*e) < 1.0001
+    ω = acos(1.)
+  else
+    error("bad ω: $(e_xyz)")
+  end
+  if abs((dot(r_xyz,e_xyz))/(r*norm(e_xyz))) < 1.
+    ν = acos((dot(r_xyz,e_xyz))/(r*norm(e_xyz)))
+  elseif (dot(r_xyz,e_xyz))/(r*norm(e_xyz)) < 1.0001
+    ν = acos(1.)
+  else
+    error("bad ν")
+  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/laguerre-conway.jl

@@ -0,0 +1,57 @@
+function laguerre_conway(state::Vector{Float64}, μ::Float64, time::Float64)
+
+  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-12
+      Δ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-12
+      Δ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/main.jl

@@ -0,0 +1,6 @@
+using LinearAlgebra
+
+include("constants.jl")
+include("conversions.jl")
+include("sft.jl")
+include("laguerre-conway.jl")

julia/sft.jl

@@ -0,0 +1,12 @@
+"""
+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::Float64)
+  return duty_cycle*num_thrusters*max_thrust*(tf-t0)/mass
+end
+