Added a sun test case and interpolation

julia/src/errors.jl

@@ -2,6 +2,9 @@ struct LaGuerreConway_Error <: Exception end
 
 struct ΔVsize_Error <: Exception end
 
+struct Convergence_Error <: Exception end
+Base.showerror(io::IO, e::Convergence_Error) = print(io, "NLP solver didn't converge...")
+
 struct GenOrbit_Error <: Exception end
 Base.showerror(io::IO, e::GenOrbit_Error) = print(io, "Infinite Loop trying to generate the init orbit")
 

julia/src/inner_loop/phase.jl

@@ -22,7 +22,7 @@ function solve_phase( start::Vector{Float64},
       F[1:6, 1] .= prop(tanh.(x), start, craft, tof, primary)[2][1:6] .- final[1:6]
     catch e
       # If the error is due to something natural, just imply a penalty
-      if isa(Mass_Error, e) || isa(PropOne_Error, e)
+      if isa(e, Mass_Error) || isa(e, PropOne_Error)
         F .= 10000000.0
       else
         rethrow()
@@ -31,6 +31,7 @@ function solve_phase( start::Vector{Float64},
   end
 
   result = nlsolve(f!, atanh.(x0), ftol=tol, autodiff=:forward, iterations=num_iters)
+  converged(result) || throw(Convergence_Error())
   result.zero = tanh.(result.zero)
 
   return result

julia/src/inner_loop/propagator.jl

@@ -38,40 +38,30 @@ function prop(ΔVs::Matrix{T},
               state::Vector{Float64},
               craft::Sc,
               time::Float64,
-              primary::Body=Sun) where T <: Real
+              primary::Body=Sun;
+              interpolate::Bool=false) where T <: Real
 
   size(ΔVs)[2] == 3 || throw(ΔVsize_Error())
   n = size(ΔVs)[1]
 
-  x_states = Vector{T}()
-  y_states = Vector{T}()
-  z_states = Vector{T}()
-  dx_states = Vector{T}()
-  dy_states = Vector{T}()
-  dz_states = Vector{T}()
-  masses = Vector{T}()
+  states = [ Vector{T}(),Vector{T}(),Vector{T}(),Vector{T}(),Vector{T}(),Vector{T}(),Vector{T}() ]
 
-  push!(x_states, state[1])
-  push!(y_states, state[2])
-  push!(z_states, state[3])
-  push!(dx_states, state[4])
-  push!(dy_states, state[5])
-  push!(dz_states, state[6])
-  push!(masses, state[7])
+  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
     state = prop_one(ΔVs[i,:], state, craft, time/n, primary)
-    push!(x_states, state[1])
-    push!(y_states, state[2])
-    push!(z_states, state[3])
-    push!(dx_states, state[4])
-    push!(dy_states, state[5])
-    push!(dz_states, state[6])
-    push!(masses, state[7])
+    for j in 1:7 push!(states[j], state[j]) end
     state[7] >= craft.dry_mass || throw(Mass_Error(state[7]))
   end
 
-  return [x_states, y_states, z_states, dx_states, dy_states, dz_states, masses], state
+  return states, state
 
 end
 

julia/test/inner_loop/phase.jl

@@ -4,11 +4,12 @@
 
   using NLsolve, PlotlyJS
 
+  # First we'll test an Earth case
   # Initial Setup
-  T = rand( 2hour : second : 12hour)
-  revs = 5
-  n = 10revs
-  start_mass = 12_000.
+  T = rand( 8hour : second : day )
+  revs = 3
+  n = 10 * revs
+  start_mass = 10_000.
 
   # A simple orbit raising
   start = gen_orbit(T, start_mass, Earth)
@@ -19,7 +20,7 @@
   # This should be close enough to converge
   thrust_guess = spiral(0.88, n, start, test_sc, revs*T, Earth)
   result = solve_phase(start, final, test_sc, revs*T, thrust_guess, Earth)
-  calc_path, calc_final = prop(result.zero, start, test_sc, revs*T, Earth)
+  calc_path, calc_final = prop(result.zero, start, test_sc, revs*T, Earth, interpolate=true)
 
   # Test
   @test converged(result)
@@ -34,6 +35,35 @@
   fig = plot_orbits(paths, Earth,
                     labels=["init", "transit", "post-transit", "final"],
                     colors=["#FFF","#F44","#4F4","#44F"])
-  savefig(fig, "../plots/nlp_test.html")
+  savefig(fig, "../plots/nlp_test_earth.html")
+
+  # Now the sun case
+  T = rand( 0.5year : hour : 1.5year )
+  n = 20
+  start_mass = 12_000.
+
+  # A simple orbit raising again, to keep things simple
+  start = gen_orbit(T, start_mass)
+  thrust = spiral(0.6, n, start, bepi, revs*T)
+  final = prop(thrust, start, bepi, revs*T)[2]
+  new_T = 2π * √( xyz_to_oe(final, Sun.μ)[1]^3 / Sun.μ )
+
+  # This should be close enough to converge
+  thrust_guess = spiral(0.55, n, start, bepi, revs*T)
+  result = solve_phase(start, final, bepi, revs*T, thrust_guess)
+  calc_path, calc_final = prop(result.zero, start, bepi, revs*T, interpolate=true)
+
+  # Test
+  @test converged(result)
+  @test norm(calc_final[1:6] - final[1:6]) < 1e-5
+
+  # Plot
+  paths = Pathlist()
+  push!(paths, prop(start, T), calc_path, prop(calc_final, T), prop(final, T) )
+  fig = plot_orbits(paths,
+                    labels=["init", "transit", "post-transit", "final"],
+                    colors=["#FFF","#F44","#4F4","#44F"])
+  savefig(fig, "../plots/nlp_test_sun.html")
+
 
 end