I THINK that the single shooter is ok. I can always improve it later

julia/conversions.jl

@@ -41,32 +41,24 @@ function xyz_to_oe(cart_vec::Vector,μ::Real)
   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")
+    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))
-  elseif dot(n_xyz,[1,0,0])/norm(n_xyz) < 1.0001
-    Ω = acos(1.)
   else
-    error("bad Ω")
+    Ω = acos(1.)
   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)")
+    ω = 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)))
-  elseif (dot(r_xyz,e_xyz))/(r*norm(e_xyz)) < 1.0001
-    ν = acos(1.)
   else
-    error("bad ν")
+    ν = acos(1.)
   end
   Ω = dot(n_xyz,[0,1,0]) > 0. ? Ω : -Ω
   ω = dot(e_xyz,[0,0,1]) > 0. ? ω : -ω

julia/laguerre-conway.jl

@@ -1,4 +1,4 @@
-function laguerre_conway(state::Vector{Float64}, μ::Float64, time::Float64)
+function laguerre_conway(state::Vector{T}, μ::Float64, time::Float64) where T <: Real
 
   n = 5                               # Choose LaGuerre-Conway "n"
   i = 0
@@ -14,7 +14,7 @@ function laguerre_conway(state::Vector{Float64}, μ::Float64, time::Float64)
   if a > 0    # Elliptical
     ΔM = ΔE_new = √(μ) / sqrt(a^3) * time
     ΔE = 1000
-    while abs(ΔE - ΔE_new) > 1e-12
+    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)
@@ -32,7 +32,7 @@ function laguerre_conway(state::Vector{Float64}, μ::Float64, time::Float64)
     ΔN = √(μ) / sqrt(-a^3) * time
     ΔH = 0
     ΔH_new = time < 0 ? -1 : 1
-    while abs(ΔH - ΔH_new) > 1e-12
+    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)

julia/main.jl

@@ -1,4 +1,4 @@
-using LinearAlgebra, JuMP, Ipopt
+using LinearAlgebra, ForwardDiff, Blink, PlotlyJS
 
 include("constants.jl")
 include("conversions.jl")
@@ -6,3 +6,5 @@ include("spacecraft.jl")
 include("plotting.jl")
 include("laguerre-conway.jl")
 include("propagator.jl")
+include("result.jl")
+include("single_shoot.jl")

julia/propagator.jl

@@ -1,12 +1,12 @@
 """
 Maximum ΔV that a spacecraft can impulse for a given single time step
 """
-function max_ΔV(duty_cycle::Float64,
+function max_ΔV(duty_cycle::T,
                 num_thrusters::Int,
-                max_thrust::Float64,
-                tf::Float64,
-                t0::Float64,
-                mass::Float64)
+                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
 
@@ -14,33 +14,49 @@ 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_unit::Vector{Float64},
-                  state::Vector{Float64},
+function prop_one(ΔV::Vector{T},
+                  state::Vector{S},
                   duty_cycle::Float64,
                   num_thrusters::Int,
                   max_thrust::Float64,
-                  mass::Float64,
+                  mass::S,
                   mass_flow_rate::Float64,
                   μ::Float64,
-                  time::Float64)
+                  time::Float64) where {T <: Real, S <: Real}
 
-  if norm(ΔV_unit) > 1.
-    throw(ErrorException("ΔV input is too high"))
+  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
-  halfway = laguerre_conway(state, μ, time/2)
-  halfway[4:6] += ΔV_unit * max_ΔV(duty_cycle, num_thrusters, max_thrust, time, 0., mass)
-  return laguerre_conway(halfway, μ, time/2), mass - mass_flow_rate*norm(ΔV_unit)*time
+
+  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{Float64},
-                  state::Vector{Float64},
+function prop_one(ΔV_unit::Vector{T},
+                  state::Vector{S},
                   craft::Sc,
                   μ::Float64,
-                  time::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)
@@ -49,7 +65,7 @@ end
 """
 This propagates over a given time period, with a certain number of intermediate steps
 """
-function prop(ΔV_units::Vector{Vector{Float64}},
+function prop(ΔVs::Matrix{T},
               state::Vector{Float64},
               duty_cycle::Float64,
               num_thrusters::Int,
@@ -57,15 +73,13 @@ function prop(ΔV_units::Vector{Vector{Float64}},
               mass::Float64,
               mass_flow_rate::Float64,
               μ::Float64,
-              time::Float64,
-              n::Int)
+              time::Float64) where T <: Real
 
-  if length(ΔV_units) != n
-    throw(ExceptionError("Bad number of ΔV vectors"))
-  end
+  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(ΔV_units[i], state, duty_cycle, num_thrusters, max_thrust, mass,
+    state, mass = prop_one(ΔVs[i,:], state, duty_cycle, num_thrusters, max_thrust, mass,
                            mass_flow_rate, μ, time/n)
   end
 
@@ -76,22 +90,20 @@ end
 """
 The same function, using Scs
 """
-function prop(ΔV_units::Vector{Vector{Float64}},
+function prop(ΔVs::AbstractArray{T},
               state::Vector{Float64},
               craft::Sc,
               μ::Float64,
-              time::Float64,
-              n::Int)
+              time::Float64) where T <: Real
 
-  if length(ΔV_units) != n
-    throw(ExceptionError("Bad number of ΔV vectors"))
-  end
+  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(ΔV_units[i], state, craft, μ, time/n)
+    state, craft = prop_one(ΔVs[i,:], state, craft, μ, time/n)
     states = [states; state']
     masses = [masses, craft.mass]
   end

julia/result.jl

@@ -0,0 +1,16 @@
+struct Result
+  converged::Bool
+  ΔVs::Matrix{Float64}
+  start::Vector{Float64}
+  final::Vector{Float64}
+end
+
+function Result(name::String)
+  if name == "test_converged"
+    return Result(true, Matrix(undef,0,0), Vector{Float64}(), Vector{Float64}())
+  elseif name == "test_unconverged"
+    return Result(false, Matrix(undef,0,0), Vector{Float64}(), Vector{Float64}())
+  else
+    throw(ErrorException("Bad result name"))
+  end
+end

julia/single_shoot.jl

@@ -0,0 +1,33 @@
+using NLsolve
+
+function treat_inputs(x::AbstractVector, n::Int)
+  inputs = reshape(copy(x),(3,n))'
+  for i in 1:n
+    inputs[i,1] = 0.5*tanh(inputs[i,1]) + 0.5
+    inputs[i,2] = π*tanh(inputs[i,2])
+    inputs[i,3] = π*tanh(inputs[i,3])/2
+  end
+  return inputs
+end
+
+function single_shoot(start::Vector{Float64},
+                      final::Vector{Float64},
+                      craft::Sc,
+                      μ::Float64,
+                      t0::Float64,
+                      tf::Float64,
+                      n::Int,
+                      x0::AbstractVector,
+                      tol=1e-2)
+
+  function f!(F,x)
+    F[1:6] .= prop(treat_inputs(x,n), start, craft, μ, tf-t0)[1][end,:] - final
+    F[7:3n] .= 0.
+    # if typeof(F[1]) == Float64 println(F[1:6]) end
+    # if typeof(F[1]) == Float64 println(treat_inputs(x,n)[1:8,1]) end
+    
+  end
+
+  return nlsolve(f!, x0, ftol=tol, autodiff=:forward, iterations=10_000)
+
+end

julia/spacecraft.jl

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

julia/test/plotting.jl

@@ -2,15 +2,16 @@
 
   # First some setup
   sc = Sc("test")
-  T = rand(3600*1.5:0.01:3600*4)
+  T = rand(3600*2:0.01:3600*4)
   start = oe_to_xyz([ (μs["Earth"]*(T/(2π))^2)^(1/3),
-                       0.3,
+                       0.1,
                        π/4,
                        0.,
                        0.,
                        1. ], μs["Earth"])
-  ΔVs = [ [1., 1., 1.]/20 for _ in 1:40 ]
-  path = prop(ΔVs, start, sc, μs["Earth"], 0.9T, 40)[1]
+  n = 100
+  ΔVs = repeat([0.5, 0., 0.]', outer=(n,1))
+  path = prop(ΔVs, start, sc, μs["Earth"], 3T)[1]
   p = plot_orbits([path])
   savefig(p,"plot_test.html")
   @test typeof(p) == PlotlyJS.SyncPlot

julia/test/propagator.jl

@@ -15,19 +15,23 @@
 
   # Test that Laguerre-Conway is the default propagator for spacecrafts
   craft = Sc("no_thrust")
+  start_mass = craft.mass
   state, craft = prop_one([0., 0., 0.], start, craft, μs["Earth"], stepsize)
   @test laguerre_conway(start, μs["Earth"], stepsize) ≈ state
-  @test craft.mass == 1000.
+  @test craft.mass == start_mass
 
   # Test that mass is reduced properly
   craft = Sc("test")
   start_mass = craft.mass
-  state, craft = prop_one([1., 1., 1.]/√(3), start, craft, μs["Earth"], stepsize)
+  state, craft = prop_one([1., 0., 0.], start, craft, μs["Earth"], stepsize)
   @test craft.mass == start_mass - craft.mass_flow_rate*stepsize
   
   # Test that a bad ΔV throws an error
   craft = Sc("test")
   start_mass = craft.mass
-  @test_throws ErrorException prop_one([1., 1., -1.], start, craft, μs["Earth"], stepsize)
+  @test_throws ErrorException prop_one([1.5, 0., 0.], start, craft, μs["Earth"], stepsize)
+
+  # Test that a full propagation doesn't take too long
+
 
 end

julia/test/result.jl

@@ -0,0 +1,16 @@
+@testset "Result Construction" begin
+
+  # Test that the standard results can be created
+  result = Result("test_converged")
+  @test result.converged == true
+  @test length(result.ΔVs) == 0
+  @test length(result.start) == 0
+  @test length(result.final) == 0
+  
+  result = Result("test_unconverged")
+  @test result.converged == false
+  @test length(result.ΔVs) == 0
+  @test length(result.start) == 0
+  @test length(result.final) == 0
+  
+end

julia/test/single_shoot.jl

@@ -0,0 +1,31 @@
+@testset "Single Shooting" begin
+
+  # Initial Setup
+  sc = Sc("test")
+  a = rand(15000:1.:40000)
+  e = rand(0.01:0.01:0.5)
+  i = rand(0.01:0.01:π/6)
+  T = 2π*√(a^3/μs["Earth"])
+  n = 50
+
+  # A simple orbit raising
+  start = oe_to_xyz([ a, e, i, 0., 0., 0. ], μs["Earth"])
+  ΔVs = repeat([0.6, 0., 0.]', outer=(n,1))
+  final = prop(ΔVs, start, sc, μs["Earth"], T)[1][end,:]
+
+  # This should be close enough to 0.6
+  x0 = repeat([atanh((0.4-0.5)/0.5), 0., 0.], n)
+  result = single_shoot(start, final, sc, μs["Earth"], 0.0, T, n, x0)
+
+  # Test and plot
+  @test converged(result)
+  path1 = prop(zeros((100,3)), start, sc, μs["Earth"], T)[1]
+  path2, mass = prop(treat_inputs(result.zero, n), start, sc, μs["Earth"], T)
+  path3 = prop(zeros((100,3)), path2[end,:], sc, μs["Earth"], T)[1]
+  savefig(plot_orbits([path1, path2, path3]), "single_shoot_test.html")
+  if converged(result)
+    @test norm(path2[end,:] - final) < 2e-2
+    sc = Sc("test")
+  end
+
+end

julia/test/spacecraft.jl

@@ -2,14 +2,14 @@
 
   # Test that the standard spacecraft can be created
   craft = Sc("test")
-  @test craft.mass == 1000.
+  @test craft.mass == 10000.
   @test craft.mass_flow_rate == 0.01
-  @test craft.max_thrust == 0.1
+  @test craft.max_thrust == 0.05
   @test craft.num_thrusters == 2
   @test craft.duty_cycle == 1.
   
   craft = Sc("no_thrust")
-  @test craft.mass == 1000.
+  @test craft.mass == 10000.
   @test craft.mass_flow_rate == 0.01
   @test craft.max_thrust == 0.
   @test craft.num_thrusters == 0

julia/test/test.jl

@@ -10,6 +10,8 @@ include("../main.jl")
   include("laguerre-conway.jl")
   include("propagator.jl")
   include("plotting.jl")
+  include("result.jl")
+  include("single_shoot.jl")
 end
 
 print()