Slow but steady progress on the mbh

julia/src/lamberts.jl

@@ -0,0 +1,72 @@
+using Dates
+
+"""
+I didn't want to have to write this...I'm going to try to do this as quickly as possible from old,
+bad code
+
+Convenience function for solving lambert's problem
+"""
+function lamberts(planet1::Body,planet2::Body,leave::DateTime,arrive::DateTime)
+
+  time_leave = utc2et(Dates.format(leave,"yyyy-mm-ddTHH:MM:SS"))
+  time_arrive = utc2et(Dates.format(arrive,"yyyy-mm-ddTHH:MM:SS"))
+  tof_req = time_arrive - time_leave
+
+  state1 = [spkssb(planet1.id, time_leave, "ECLIPJ2000"); 0.0]
+  state2 = [spkssb(planet2.id, time_arrive, "ECLIPJ2000"); 0.0]
+
+  r1 = state1[1:3] ; r1mag = norm(r1)
+  r2 = state2[1:3] ; r2mag = norm(r2)
+  μ = Sun.μ
+
+  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 = abs(dθ) > π ? -1 : 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 - state1[4:6], vft - state2[4:6], tof_req
+
+end