Switching to my linux box

julia/src/Thesis.jl

@@ -1,6 +1,9 @@
 module Thesis
 
-  using LinearAlgebra, ForwardDiff, PlotlyJS
+  using LinearAlgebra, ForwardDiff, PlotlyJS, SPICE
+
+  furnsh("../../SPICE/naif0012.tls")
+  furnsh("../../SPICE/de430.bsp")
 
   include("./constants.jl")
   include("./spacecraft.jl")
@@ -12,4 +15,5 @@ module Thesis
   include("./inner_loop/monotonic_basin_hopping.jl")
   include("./inner_loop/phase.jl")
   include("./inner_loop/inner_loop.jl")
+
 end

julia/src/constants.jl

@@ -32,13 +32,13 @@ export μs, G, GMs, μ, rs, as, es, AU
     return μ(GMs[primary]/G, GMs[secondary]/G)
   end
 
-  GMs = Dict(
+  const GMs = Dict(
   "Sun" => 132712440041.93938,
   "Earth" => 398600.435436,
   "Moon" => 4902.800066)
 
 # Radii
-  rs = Dict(
+  const rs = Dict(
   "Sun" => 696000.,
   "Mercury" => 2439.,
   "Venus" => 6052.,
@@ -52,7 +52,7 @@ export μs, G, GMs, μ, rs, as, es, AU
   "Pluto" => 1151.)
 
 # Semi Major Axes
-  as = Dict(
+  const as = Dict(
   "Mercury" => 57909083.,
   "Venus" => 108208601.,
   "Earth" => 149598023.,
@@ -65,12 +65,12 @@ export μs, G, GMs, μ, rs, as, es, AU
   "Pluto" => 5915799000.)
 
 # Eccentricities
-  es = Dict(
+  const es = Dict(
   "Earth" => 0.016708617,
   "Moon" => 0.0549)
 
 # J2 for basic oblateness
-  j2s = Dict(
+  const j2s = Dict(
   "Mercury" => 0.00006,
   "Venus" => 0.000027,
   "Earth" => 0.0010826269,
@@ -83,7 +83,7 @@ export μs, G, GMs, μ, rs, as, es, AU
   "Pluto" => 0.)
 
 # These are just the colors for plots
-  p_colors = Dict(
+  const p_colors = Dict(
   "Sun" => :Electric,
   "Mercury" => :heat,
   "Venus" => :turbid,
@@ -96,5 +96,19 @@ export μs, G, GMs, μ, rs, as, es, AU
   "Neptune" => :ice,
   "Pluto" => :matter)
 
-  AU = 149597870.691 #km
-  init_STM = vec(Matrix{Float64}(I,6,6))
+  const ids = Dict(
+    "Sun" => 10,
+    "Mercury" => 1,
+    "Venus" => 2,
+    "Earth" => 399,
+    "Moon" => 301,
+    "Mars" => 4,
+    "Jupiter" => 5,
+    "Saturn" => 6,
+    "Uranus" => 7,
+    "Neptune" => 8,
+    "Pluto" => 9,
+  )
+
+  const AU = 149597870.691 #km
+  const init_STM = vec(Matrix{Float64}(I,6,6))

julia/src/inner_loop/inner_loop.jl

@@ -7,19 +7,50 @@ This is it. The outer function call for the inner loop. After this is done,
 there's only the outer loop left to do. And that's pretty easy.
 """
 function inner_loop(launch_date::DateTime,
-                    RLA::Float64,
-                    DLA::Float64,
-                    phases::Vector{Phase})
+                    phases::Vector{Phase};
+                    min_flyby::Float64=1000.)
    
   # First we need to do some quick checks that the mission is well formed
   for i in 1:length(phases)
     if i == 1
       @assert phases[i].from_planet == "Earth"
     else
-      @assert phases[i].from_planet == phases[i-1].to_planet
-      @assert phases[i].v∞_outgoing == phases[i-1].v∞_incoming
+      # Check that the planet is valid
+      if phases[i].from_planet ∉ keys(μs)
+        error("Planet is not valid: ", phases[i].from_planet)
+      # Check that there is only one flyby planet
+      elseif phases[i].from_planet != phases[i-1].to_planet
+        fromP, toP = phases[i].from_planet, phases[i-1].to_planet
+        error("Planets don't match up: (phase $(i)) $(fromP) / (phase $(i-1)) $(toP)")
+      # Check that v∞_in == v∞_out
+      elseif !isapprox(norm(phases[i].v∞_outgoing), norm(phases[i-1].v∞_incoming))
+        norm_incoming = norm(phases[i-1].v∞_incoming)
+        norm_outgoing = norm(phases[i].v∞_outgoing)
+        error("""Norms of vectors aren't equal:
+                 Phase $(i-1): $(phases[i-1].v∞_incoming) / norm: $(norm_incoming)
+                 Phase $(i): $(phases[i].v∞_outgoing) / norm: $(norm_outgoing)""")
+      end
+      # Check that the approach is not too low
+      v∞ = norm(phases[i].v∞_outgoing)
+      δ = acos((phases[i].v∞_outgoing ⋅ phases[i-1].v∞_incoming)/v∞^2)
+      flyby = μs[phases[i].from_planet]/v∞^2 * (1/sin(δ/2) - 1)
+      true_min = rs[phases[i].from_planet] + min_flyby
+      if flyby <= true_min
+        error("Flyby was too low between phase $(i-1) and $(i): $(flyby) / $(true_min)")
+      end
+
     end
   end
+
+  time = utc2et(Dates.format(launch_date,"yyyy-mm-ddTHH:MM:SS"))
+  for phase in phases
+    planet1_state = spkssb(ids[phase.from_planet], time, "J2000")
+    println(planet1_state)
+    time += phase.time_of_flight
+    planet2_state = spkssb(ids[phase.to_planet], time, "J2000")
+    println(planet2_state)
+  end
+
       
-  return RLA + DLA + phases[1].v∞_incoming
+  return phases[1].v∞_incoming[1]
 end

julia/src/inner_loop/phase.jl

@@ -3,7 +3,7 @@ export Phase
 struct Phase
   from_planet::String
   to_planet::String
-  time_of_flight::Float64   # seconds
-  v∞_outgoing::Float64      # Km/s
-  v∞_incoming::Float64      # Km/s
+  time_of_flight::Float64           # seconds
+  v∞_outgoing::Vector{Float64}      # Km/s
+  v∞_incoming::Vector{Float64}      # Km/s
 end