Getting pretty close on the final review

julia/lamberts_fig.jl

@@ -0,0 +1,17 @@
+using PlotlyJS: savefig
+
+r1 = oe_to_xyz([Earth.a, 0.6, 0.001, 0., 0., deg2rad(10.)], Sun.μ)[1:3]
+r2 = oe_to_xyz([Earth.a, 0.6, 0.001, 0., 0., deg2rad(140.)], Sun.μ)[1:3]
+tof = 0.5year
+v1 = Thesis.lamberts1(r1,r2,tof)[1]
+v2 = Thesis.lamberts2(r1,r2,tof)[1]
+state1 = [ r1; v1; 1000. ]
+state2 = [ r1; v2; 1000. ]
+path1 = prop(state1, tof)
+path2 = prop(state2, tof)
+plot([path2, path1],
+     title="Lambert's Problem Solutions",
+     colors=["#F00", "#0FF"],
+     labels=["Type I", "Type II"],
+     mode="light")
+

julia/src/utilities/lamberts.jl

@@ -71,6 +71,120 @@ function lamberts(planet1::Body,planet2::Body,leave::DateTime,arrive::DateTime)
 
 end
 
+function lamberts1(r1::Vector{Float64},r2::Vector{Float64},tof_req::Float64)
+  μ = Sun.μ
+  r1mag = norm(r1)
+  r2mag = norm(r2)
+
+  cos_dθ = dot(r1,r2)/(r1mag*r2mag)
+  dθ = atan(r2[2],r2[1]) - atan(r1[2],r1[1])
+  dθ = dθ > 2π ? dθ-2π : dθ
+  dθ = dθ < 0.0 ? dθ+2π : dθ
+  DM = -1
+  A = DM * √(r1mag * r2mag * (1 + cos_dθ))
+  dθ == 0 || A == 0 && error("Can't solve Lambert's Problem")
+
+  ψ, c2, c3 = 0, 1//2, 1//6
+  ψ_down = -4π ; ψ_up = 4π^2
+  y = r1mag + r2mag + (A*(ψ*c3 - 1)) / √(c2) ; χ = √(y/c2)
+  tof = ( χ^3*c3 + A*√(y) ) / √(μ)
+
+  i = 0
+  while abs(tof-tof_req) > 1e-2
+    y = r1mag + r2mag + (A*(ψ*c3 - 1)) / √(c2)
+    while y/c2 <= 0
+      # println("You finally hit that weird issue... ")
+      ψ += 0.1
+      if ψ > 1e-6
+        c2 = (1 - cos(√(ψ))) / ψ ; c3 = (√(ψ) - sin(√(ψ))) / √(ψ^3)
+      elseif ψ < -1e-6
+        c2 = (1 - cosh(√(-ψ))) / ψ ; c3 = (-√(-ψ) + sinh(√(-ψ))) / √((-ψ)^3)
+      else
+        c2 = 1//2 ; c3 = 1//6
+      end
+      y = r1mag + r2mag + (A*(ψ*c3 - 1)) / √(c2)
+    end
+    χ = √(y/c2)
+
+    tof = ( c3*χ^3 + A*√(y) ) / √(μ)
+    tof < tof_req ? ψ_down = ψ : ψ_up = ψ
+    ψ = (ψ_up + ψ_down) / 2
+
+    if ψ > 1e-6
+      c2 = (1 - cos(√(ψ))) / ψ ; c3 = (√(ψ) - sin(√(ψ))) / √(ψ^3)
+    elseif ψ < -1e-6
+      c2 = (1 - cosh(√(-ψ))) / ψ ; c3 = (-√(-ψ) + sinh(√(-ψ))) / √((-ψ)^3)
+    else
+      c2 = 1//2 ; c3 = 1//6
+    end
+
+    i += 1
+    i > 500 && return [NaN,NaN,NaN],[NaN,NaN,NaN]
+  end
+
+  f = 1 - y/r1mag ; g_dot = 1 - y/r2mag ; g = A * √(y/μ)
+  v0t = (r2 - f*r1)/g ; vft = (g_dot*r2 - r1)/g
+  return v0t, vft, tof_req
+
+end
+
+function lamberts2(r1::Vector{Float64},r2::Vector{Float64},tof_req::Float64)
+  μ = Sun.μ
+  r1mag = norm(r1)
+  r2mag = norm(r2)
+
+  cos_dθ = dot(r1,r2)/(r1mag*r2mag)
+  dθ = atan(r2[2],r2[1]) - atan(r1[2],r1[1])
+  dθ = dθ > 2π ? dθ-2π : dθ
+  dθ = dθ < 0.0 ? dθ+2π : dθ
+  DM = 1
+  A = DM * √(r1mag * r2mag * (1 + cos_dθ))
+  dθ == 0 || A == 0 && error("Can't solve Lambert's Problem")
+
+  ψ, c2, c3 = 0, 1//2, 1//6
+  ψ_down = -4π ; ψ_up = 4π^2
+  y = r1mag + r2mag + (A*(ψ*c3 - 1)) / √(c2) ; χ = √(y/c2)
+  tof = ( χ^3*c3 + A*√(y) ) / √(μ)
+
+  i = 0
+  while abs(tof-tof_req) > 1e-2
+    y = r1mag + r2mag + (A*(ψ*c3 - 1)) / √(c2)
+    while y/c2 <= 0
+      # println("You finally hit that weird issue... ")
+      ψ += 0.1
+      if ψ > 1e-6
+        c2 = (1 - cos(√(ψ))) / ψ ; c3 = (√(ψ) - sin(√(ψ))) / √(ψ^3)
+      elseif ψ < -1e-6
+        c2 = (1 - cosh(√(-ψ))) / ψ ; c3 = (-√(-ψ) + sinh(√(-ψ))) / √((-ψ)^3)
+      else
+        c2 = 1//2 ; c3 = 1//6
+      end
+      y = r1mag + r2mag + (A*(ψ*c3 - 1)) / √(c2)
+    end
+    χ = √(y/c2)
+
+    tof = ( c3*χ^3 + A*√(y) ) / √(μ)
+    tof < tof_req ? ψ_down = ψ : ψ_up = ψ
+    ψ = (ψ_up + ψ_down) / 2
+
+    if ψ > 1e-6
+      c2 = (1 - cos(√(ψ))) / ψ ; c3 = (√(ψ) - sin(√(ψ))) / √(ψ^3)
+    elseif ψ < -1e-6
+      c2 = (1 - cosh(√(-ψ))) / ψ ; c3 = (-√(-ψ) + sinh(√(-ψ))) / √((-ψ)^3)
+    else
+      c2 = 1//2 ; c3 = 1//6
+    end
+
+    i += 1
+    i > 500 && return [NaN,NaN,NaN],[NaN,NaN,NaN]
+  end
+
+  f = 1 - y/r1mag ; g_dot = 1 - y/r2mag ; g = A * √(y/μ)
+  v0t = (r2 - f*r1)/g ; vft = (g_dot*r2 - r1)/g
+  return v0t, vft, tof_req
+
+end
+
 function porkchop(planet1::Body, planet2::Body, departures::Vector{DateTime}, arrivals::Vector{DateTime})
   v∞_in = [ norm(lamberts(planet1, planet2, depart, arrive)[2]) for depart in departures, arrive in arrivals ]
   v∞_out = [ norm(lamberts(planet1, planet2, depart, arrive)[1]) for depart in departures, arrive in arrivals ]

julia/src/utilities/plotting.jl

@@ -89,14 +89,22 @@ function standard_layout(limit::Float64, title::AbstractString; mode="dark")
       plot_bgcolor="rgba(255,255,255,0.0)",
       scene_aspectmode = "data",
       scene_xaxis = attr(
+        title = "X (km)",
+        exponentformat = "power",
         autorange = true,
         color="rgb(0,0,0)"
       ),
       scene_yaxis = attr(
+        title = "Y (km)",
+        exponentformat = "power",
         autorange = true,
         color="rgb(0,0,0)"
       ),
-      scene_zaxis_visible = false
+      scene_zaxis = attr(
+        title = "Z (km)",
+        exponentformat = "power",
+        visible = false,
+      ),
     )
   end
 end
@@ -287,7 +295,7 @@ function plot(paths::Vector{Matrix{Real}};
   traces = [ trace... ]
   for i = 2:length(paths)
     color = colors === nothing ? random_color() : colors[i]
-    trace, new_limit = gen_plot(paths[i],label=labels[i],color=colors[i],markers=markers)
+    trace, new_limit = gen_plot(paths[i],label=labels[i],color=color,markers=markers)
     push!(traces, trace...)
     limit = max(limit, new_limit)
   end