Got plotting working

julia/main.jl

@@ -3,6 +3,6 @@ using LinearAlgebra, JuMP, Ipopt
 include("constants.jl")
 include("conversions.jl")
 include("spacecraft.jl")
-include("sft.jl")
+include("plotting.jl")
 include("laguerre-conway.jl")
 include("propagator.jl")

julia/plotting.jl

@@ -0,0 +1,61 @@
+# ------------------------------------------------------------------------------
+# PLOTTING FUNCTIONS
+
+using Random, PlotlyJS, Base.Iterators
+
+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")
+
+  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 = 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/propagator.jl

@@ -1,3 +1,15 @@
+"""
+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
+
 """
 This function propagates the spacecraft forward in time 1 Sim-Flanagan step (of variable length of time),
 applying a thrust in the center.
@@ -26,11 +38,64 @@ A convenience function for using spacecraft. Note that this function outputs a s
 """
 function prop_one(ΔV_unit::Vector{Float64},
                   state::Vector{Float64},
-                  craft::sc,
+                  craft::Sc,
                   μ::Float64,
                   time::Float64)
   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)
+  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(ΔV_units::Vector{Vector{Float64}},
+              state::Vector{Float64},
+              duty_cycle::Float64,
+              num_thrusters::Int,
+              max_thrust::Float64,
+              mass::Float64,
+              mass_flow_rate::Float64,
+              μ::Float64,
+              time::Float64,
+              n::Int)
+
+  if length(ΔV_units) != n
+    throw(ExceptionError("Bad number of ΔV vectors"))
+  end
+
+  for i in 1:n
+    state, mass = prop_one(ΔV_units[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(ΔV_units::Vector{Vector{Float64}},
+              state::Vector{Float64},
+              craft::Sc,
+              μ::Float64,
+              time::Float64,
+              n::Int)
+
+  if length(ΔV_units) != n
+    throw(ExceptionError("Bad number of ΔV vectors"))
+  end
+
+  states = state'
+  masses = craft.mass
+
+  for i in 1:n
+    state, craft = prop_one(ΔV_units[i], state, craft, μ, time/n)
+    states = [states; state']
+    masses = [masses, craft.mass]
+  end
+
+  return states, masses
+
+end

julia/sft.jl

@@ -1,12 +0,0 @@
-"""
-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
-

julia/spacecraft.jl

@@ -1,4 +1,4 @@
-struct sc
+struct Sc
   mass::Float64
   mass_flow_rate::Float64
   max_thrust::Float64
@@ -6,11 +6,11 @@ struct sc
   duty_cycle::Float64
 end
 
-function sc(name::String)
+function Sc(name::String)
   if name == "test"
-    return sc(1000., 0.01, 0.1, 2, 1.)
+    return Sc(1000., 0.01, 0.1, 2, 1.)
   elseif name == "no_thrust"
-    return sc(1000., 0.01, 0., 0, 0.)
+    return Sc(1000., 0.01, 0., 0, 0.)
   else
     throw(ErrorException("Bad sc name"))
   end

julia/test/laguerre-conway.jl

@@ -11,7 +11,7 @@
         i += 1
       end
       @test i ≈ T
-      @test orbit ≈ start
+      @test norm(orbit - start) < 1e-2
     end
     for _ in 1:5
       i = 0.
@@ -20,7 +20,7 @@
         i -= 1
       end
       @test i ≈ -T
-      @test orbit ≈ start
+      @test norm(orbit - start) < 1e-2
     end
   end
 

julia/test/plotting.jl

@@ -0,0 +1,18 @@
+@testset "Plotting" begin
+
+  # First some setup
+  sc = Sc("test")
+  T = rand(3600*1.5:0.01:3600*4)
+  start = oe_to_xyz([ (μs["Earth"]*(T/(2π))^2)^(1/3),
+                       0.3,
+                       π/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]
+  p = plot_orbits([path])
+  savefig(p,"plot_test.html")
+  @test typeof(p) == PlotlyJS.SyncPlot
+
+end

julia/test/propagator.jl

@@ -14,19 +14,19 @@
   @test laguerre_conway(start, μs["Earth"], stepsize) ≈ propped[1]
 
   # Test that Laguerre-Conway is the default propagator for spacecrafts
-  craft = sc("no_thrust")
+  craft = Sc("no_thrust")
   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 that mass is reduced properly
-  craft = sc("test")
+  craft = Sc("test")
   start_mass = craft.mass
   state, craft = prop_one([1., 1., 1.]/√(3), 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")
+  craft = Sc("test")
   start_mass = craft.mass
   @test_throws ErrorException prop_one([1., 1., -1.], start, craft, μs["Earth"], stepsize)
 

julia/test/spacecraft.jl

@@ -1,19 +1,18 @@
 @testset "Spacecraft Construction" begin
 
-  # Test that the standard spacecrafts can be created
-  craft = sc("test")
+  # Test that the standard spacecraft can be created
+  craft = Sc("test")
   @test craft.mass == 1000.
   @test craft.mass_flow_rate == 0.01
   @test craft.max_thrust == 0.1
   @test craft.num_thrusters == 2
   @test craft.duty_cycle == 1.
   
-  craft = sc("no_thrust")
+  craft = Sc("no_thrust")
   @test craft.mass == 1000.
   @test craft.mass_flow_rate == 0.01
   @test craft.max_thrust == 0.
   @test craft.num_thrusters == 0
   @test craft.duty_cycle == 0.
 
-
 end

julia/test/test.jl

@@ -2,19 +2,14 @@ using Test
 using Random
 using LinearAlgebra
 
-# Includes 
-include("../constants.jl")
-include("../conversions.jl")
-include("../spacecraft.jl")
-include("../sft.jl")
-include("../laguerre-conway.jl")
-include("../propagator.jl")
+include("../main.jl")
 
 # Tests
 @testset "All Tests" begin
   include("spacecraft.jl")
   include("laguerre-conway.jl")
   include("propagator.jl")
+  include("plotting.jl")
 end
 
 print()