Currently working on refactor, much work to do

This commit is contained in:
Connor
2021-09-25 16:10:02 -06:00
parent 7d0037f38d
commit af545ba1a7
23 changed files with 286 additions and 353 deletions

View File

@@ -1,25 +1,30 @@
using SPICE
try
furnsh("../../spice_files/naif0012.tls")
furnsh("../../spice_files/de430.bsp")
catch
furnsh("spice_files/naif0012.tls")
furnsh("spice_files/de430.bsp")
end
module Thesis
using LinearAlgebra, ForwardDiff, PlotlyJS, SPICE, Distributed
using LinearAlgebra
using ForwardDiff
using PlotlyJS
using Distributed
try
furnsh("../../SPICE/naif0012.tls")
furnsh("../../SPICE/de430.bsp")
catch
furnsh("SPICE/naif0012.tls")
furnsh("SPICE/de430.bsp")
end
include("./errors.jl")
include("./constants.jl")
include("./spacecraft.jl")
include("./conversions.jl")
include("./plotting.jl")
include("./inner_loop/laguerre-conway.jl")
include("./inner_loop/propagator.jl")
include("./inner_loop/find_closest.jl")
include("./inner_loop/monotonic_basin_hopping.jl")
include("./inner_loop/phase.jl")
include("./inner_loop/inner_loop.jl")
include("./outer_loop.jl")
# include("./inner_loop/monotonic_basin_hopping.jl")
# include("./outer_loop.jl")
end

View File

@@ -2,113 +2,36 @@
# DEFINING CONSTANTS
# -----------------------------------------------------------------------------
export μs, G, GMs, μ, rs, as, es, AU, ids
export Body, Sun, Mercury, Venus, Earth, Moon, Mars
export Jupiter, Saturn, Uranus, Neptune, Pluto
export G, AU, init_STM, hour, day, year, second
export Pathlist
# Gravitational Constants
μs = Dict(
"Sun" => 1.32712440018e11,
"Mercury" => 2.2032e4,
"Venus" => 3.257e5,
"Earth" => 3.986004415e5,
"Moon" => 4.902799e3,
"Mars" => 4.305e4,
"Jupiter" => 1.266865361e8,
"Saturn" => 3.794e7,
"Uranus" => 5.794e6,
"Neptune" => 6.809e6,
"Pluto" => 9e2)
struct Body
μ::Float64
r::Float64 # radius
color::String
id::Int # SPICE id
end
G = 6.67430e-20
const Sun = Body(1.32712440018e11, 696000., "Electric", 10)
const Mercury = Body(2.2032e4, 2439., "heat", 1)
const Venus = Body(3.257e5, 6052., "turbid", 2)
const Earth = Body(3.986004415e5, 6378.1363, "Blues", 399)
const Moon = Body(4.902799e3, 1738., "Greys", 301)
const Mars = Body(4.305e4, 3397.2, "Reds", 4)
const Jupiter = Body(1.266865361e8, 71492., "solar", 5)
const Saturn = Body(3.794e7, 60268., "turbid", 6)
const Uranus = Body(5.794e6, 25559., "haline", 7)
const Neptune = Body(6.809e6, 24764., "ice", 8)
const Pluto = Body(9e2, 1151., "matter", 9)
function μ(m1::Float64, m2::Float64)
return m2/(m1+m2)
end
const G = 6.67430e-20 #universal gravity parameter
const AU = 149597870.691 #km
const init_STM = vec(Matrix{Float64}(I,6,6))
const second = 1.
const hour = 3600.
const day = 86400.
const year = 365 * day
function μ(GM1::Float64, GM2::Float64, Grav::Float64)
return μ(GM1/Grav, GM2/Grav)
end
function μ(primary::String, secondary::String)
return μ(GMs[primary]/G, GMs[secondary]/G)
end
const GMs = Dict(
"Sun" => 132712440041.93938,
"Earth" => 398600.435436,
"Moon" => 4902.800066)
# Radii
const rs = Dict(
"Sun" => 696000.,
"Mercury" => 2439.,
"Venus" => 6052.,
"Earth" => 6378.1363,
"Moon" => 1738.,
"Mars" => 3397.2,
"Jupiter" => 71492.,
"Saturn" => 60268.,
"Uranus" => 25559.,
"Neptune" => 24764.,
"Pluto" => 1151.)
# Semi Major Axes
const as = Dict(
"Mercury" => 57909083.,
"Venus" => 108208601.,
"Earth" => 149598023.,
"Moon" => 384400.,
"Mars" => 227939186.,
"Jupiter" => 778298361.,
"Saturn" => 1429394133.,
"Uranus" => 2875038615.,
"Neptune" => 4504449769.,
"Pluto" => 5915799000.)
# Eccentricities
const es = Dict(
"Earth" => 0.016708617,
"Moon" => 0.0549)
# J2 for basic oblateness
const j2s = Dict(
"Mercury" => 0.00006,
"Venus" => 0.000027,
"Earth" => 0.0010826269,
"Moon" => 0.0002027,
"Mars" => 0.001964,
"Jupiter" => 0.01475,
"Saturn" => 0.01645,
"Uranus" => 0.012,
"Neptune" => 0.004,
"Pluto" => 0.)
# These are just the colors for plots
const p_colors = Dict(
"Sun" => "Electric",
"Mercury" => "heat",
"Venus" => "turbid",
"Earth" => "Blues",
"Moon" => "Greys",
"Mars" => "Reds",
"Jupiter" => "solar",
"Saturn" => "turbid",
"Uranus" => "haline",
"Neptune" => "ice",
"Pluto" => "matter")
const ids = Dict(
"Sun" => 10,
"Mercury" => 1,
"Venus" => 2,
"Earth" => 399,
"Moon" => 301,
"Mars" => 4,
"Jupiter" => 5,
"Saturn" => 6,
"Uranus" => 7,
"Neptune" => 8,
"Pluto" => 9,
)
const AU = 149597870.691 #km
const init_STM = vec(Matrix{Float64}(I,6,6))
Pathlist = Vector{Vector{Vector{Float64}}}

View File

@@ -1,4 +1,4 @@
export oe_to_rθh, rθh_to_xyz, oe_to_xyz, xyz_to_oe, conv_T
export oe_to_rθh, rθh_to_xyz, oe_to_xyz, xyz_to_oe, conv_T, spiral, gen_orbit
function oe_to_rθh(oe::Vector,μ::Real) :: Vector
@@ -69,6 +69,30 @@ function xyz_to_oe(cart_vec::Vector,μ::Real)
end
"""
A convenience function for generating start conditions from orbital elements
Inputs: a body, a period, and a mass
Output: a random reasonable orbit
"""
function gen_orbit(T::Float64, mass::Float64, primary::Body=Sun)
μ = primary.μ
i = rand(0.0:0.01:0.4999π)
θ = rand(0.0:0.01:2π)
i = 0
while true
i += 1
e = rand(0.0:0.01:0.5)
a = (μ * ( T/2π )^2 )
a*(1-e) < 1.1primary.r || return [ oe_to_xyz([ a, e, i, 0., 0., θ ], μ); mass ]
i < 100 || throw(GenOrbit_Error)
end
end
"""
A convenience function for generating spiral trajectories
"""
spiral(mag,n,init,sc,time,primary=Sun) = conv_T(fill(mag, n), zeros(n), zeros(n), init, sc, time, primary)
"""
Converts a series of thrust vectors from R,Θ,H frame to cartesian
"""
@@ -78,7 +102,7 @@ function conv_T(Tm::Vector{Float64},
init_state::Vector{Float64},
craft::Sc,
time::Float64,
μ::Float64)
primary::Body=Sun)
Txs = Float64[]
Tys = Float64[]
@@ -100,7 +124,7 @@ function conv_T(Tm::Vector{Float64},
for i in 1:n
mag, α, β = Tm[i], Ta[i], Tb[i]
thrust_rθh = mag * [cos(β)*sin(α), cos(β)*cos(α), sin(β)]
_,_,i,Ω,ω,ν = xyz_to_oe(state, μ)
_,_,i,Ω,ω,ν = xyz_to_oe(state, primary.μ)
θ = ω+ν
,,ci,si,, = cos(Ω),sin(Ω),cos(i),sin(i),cos(θ),sin(θ)
DCM = [*-*ci* -*-*ci* *si;
@@ -108,10 +132,10 @@ function conv_T(Tm::Vector{Float64},
si* si* ci ]
Tx, Ty, Tz = DCM*thrust_rθh
state = prop_one([Tx, Ty, Tz], state, copy(craft), μ, time/n)
state = prop_one([Tx, Ty, Tz], state, craft, time/n, primary)
push!(Txs, Tx)
push!(Tys, Ty)
push!(Tzs, Tz)
end
return Txs, Tys, Tzs
return hcat(Txs, Tys, Tzs)
end

16
julia/src/errors.jl Normal file
View File

@@ -0,0 +1,16 @@
struct LaGuerreConway_Error <: Exception end
struct ΔVsize_Error <: Exception end
struct GenOrbit_Error <: Exception end
Base.showerror(io::IO, e::GenOrbit_Error) = print(io, "Infinite Loop trying to generate the init orbit")
struct PropOne_Error <: Exception
ΔV_unit::AbstractVector
end
Base.showerror(io::IO, e::PropOne_Error) = print(io, "tried to prop a unit ΔV of: ", e.ΔV_unit)
struct Mass_Error{T} <: Exception where T <: Real
mass::T
end
Base.showerror(io::IO, e::Mass_Error) = print(io, "Mass (", e.mass, ") got too low in propagation")

View File

@@ -1,45 +0,0 @@
using NLsolve
export nlp_solve, mass_est
function mass_est(T)
ans = 0
n = Int(length(T)/3)
for i in 1:n ans += norm(T[i,:]) end
return ans/n
end
struct Result
converged::Bool
zero::Matrix{Float64}
end
function nlp_solve(start::Vector{Float64},
final::Vector{Float64},
craft::Sc,
μ::Float64,
t0::Float64,
tf::Float64,
x0::Matrix{Float64};
tol=1e-6,
num_iters=1_000)
function f!(F,x)
try
F .= 0.0
F[1:6, 1] .= prop(tanh.(x), start, copy(craft), μ, tf-t0)[2][1:6] .- final[1:6]
catch e
F .= 10000000.0
end
end
result = Result(false, zeros(size(x0)))
try
nl_results = nlsolve(f!, atanh.(x0), ftol=tol, autodiff=:forward, iterations=num_iters)
result = Result(converged(nl_results), tanh.(nl_results.zero))
catch e
end
return result
end

View File

@@ -1,5 +1,6 @@
function laguerre_conway(state::Vector{<:Real}, μ::Float64, time::Float64)
function laguerre_conway(state::Vector{<:Real}, time::Float64, primary::Body=Sun)
μ = primary.μ
n = 5 # Choose LaGuerre-Conway "n"
i = 0
@@ -22,7 +23,7 @@ function laguerre_conway(state::Vector{<:Real}, μ::Float64, time::Float64)
sign = dF >= 0 ? 1 : -1
ΔE_new = ΔE - n*F / ( dF + sign * (abs((n-1)^2*dF^2 - n*(n-1)*F*d2F )))
i += 1
if i > 100 throw(ErrorException("LaGuerre-Conway did not converge!")) end
if i > 100 throw(LaGuerreConway_Error("LaGuerre-Conway did not converge!")) end
end
F = 1 - a/r0_mag * (1-cos(ΔE))
G = a * σ0/ (μ) * (1-cos(ΔE)) + r0_mag * (a) / (μ) * sin(ΔE)
@@ -41,7 +42,7 @@ function laguerre_conway(state::Vector{<:Real}, μ::Float64, time::Float64)
sign = dF >= 0 ? 1 : -1
ΔH_new = ΔH - n*F / ( dF + sign * (abs((n-1)^2*dF^2 - n*(n-1)*F*d2F )))
i += 1
if i > 100 throw(ErrorException("LaGuerre-Conway did not converge!")) end
if i > 100 throw(LaGuerreConway_Error("LaGuerre-Conway did not converge!")) end
end
F = 1 - a/r0_mag * (1-cos(ΔH))
G = a * σ0/ (μ) * (1-cos(ΔH)) + r0_mag * (-a) / (μ) * sin(ΔH)

View File

@@ -1,5 +1,8 @@
export mbh
"""
Generates n pareto-distributed random numbers
"""
function pareto(α::Float64, n::Int)
s = rand((-1,1), (n,3))
r = rand(Float64, (n,3))
@@ -7,6 +10,10 @@ function pareto(α::Float64, n::Int)
return (s./ϵ) * (α - 1.0) ./ (ϵ ./ (ϵ .+ r)).^(-α)
end
"""
Perturbs the monotonic basin hopping decision vector
TODO: This needs to be updated
"""
function perturb(x::AbstractMatrix, n::Int)
ans = x + pareto(1.01, n)
map!(elem -> elem > 1.0 ? 1.0 : elem, ans, ans)
@@ -14,11 +21,13 @@ function perturb(x::AbstractMatrix, n::Int)
return ans
end
function new_x(n::Int)
2.0 * rand(Float64, (n,3)) .- 1.
end
function mbh(flybys::Vector{Planet})
end
function mbh(start::AbstractVector,
final::AbstractVector,
craft::Sc,

View File

@@ -1,9 +1,38 @@
export Phase
using NLsolve
export solve_phase
"""
This function will take a single phase (so an initial state, and a final state) and an initial guess
to the thrust profile and use an NLP solver to find the nearest thrust profile to that initial guess
that satisfies the final state condition
"""
function solve_phase( start::Vector{Float64},
final::Vector{Float64},
craft::Sc,
tof::Float64,
x0::Matrix{Float64},
primary::Body=Sun;
tol=1e-6,
num_iters=1_000 )
function f!(F,x)
try
F .= 0.0
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)
F .= 10000000.0
else
rethrow()
end
end
end
result = nlsolve(f!, atanh.(x0), ftol=tol, autodiff=:forward, iterations=num_iters)
result.zero = tanh.(result.zero)
return result
struct Phase
from_planet::String
to_planet::String
time_of_flight::Float64 # seconds
v∞_outgoing::Vector{Float64} # Km/s
v∞_incoming::Vector{Float64} # Km/s
end

View File

@@ -18,18 +18,15 @@ A convenience function for using spacecraft. Note that this function outputs a s
function prop_one(ΔV_unit::Vector{<:Real},
state::Vector{<:Real},
craft::Sc,
μ::Float64,
time::Float64)
time::Float64,
primary::Body=Sun)
for direction in ΔV_unit
if abs(direction) > 1.0
println(direction)
error("ΔV is impossibly high")
end
abs(direction) <= 1.0 || throw(PropOne_Error(ΔV_unit))
end
ΔV = max_ΔV(craft.duty_cycle, craft.num_thrusters, craft.max_thrust, time, 0., state[7]) * ΔV_unit
halfway = laguerre_conway(state, μ, time/2) + [zeros(3); ΔV]
final = laguerre_conway(halfway, μ, time/2)
halfway = laguerre_conway(state, time/2, primary) + [zeros(3); ΔV]
final = laguerre_conway(halfway, time/2, primary)
return [final; state[7] - craft.mass_flow_rate*norm(ΔV_unit)*time]
end
@@ -40,10 +37,10 @@ The propagator function
function prop(ΔVs::Matrix{T},
state::Vector{Float64},
craft::Sc,
μ::Float64,
time::Float64) where T <: Real
time::Float64,
primary::Body=Sun) where T <: Real
if size(ΔVs)[2] != 3 throw(ErrorException("ΔV input is wrong size")) end
size(ΔVs)[2] == 3 || throw(ΔVsize_Error())
n = size(ΔVs)[1]
x_states = Vector{T}()
@@ -63,7 +60,7 @@ function prop(ΔVs::Matrix{T},
push!(masses, state[7])
for i in 1:n
state = prop_one(ΔVs[i,:], state, craft, μ, time/n)
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])
@@ -71,12 +68,14 @@ function prop(ΔVs::Matrix{T},
push!(dy_states, state[5])
push!(dz_states, state[6])
push!(masses, state[7])
if state[7] < craft.dry_mass
println(state[7])
error("Mass is too low")
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
end
"""
Convenience function for propagating a state with no thrust
"""
prop(x::Vector{Float64}, t::Float64, p::Body=Sun) = prop(zeros(1000,3), x, no_thrust, t, p)[1]

View File

@@ -61,3 +61,9 @@ function gen_decision_vector(launch_range::Vector{DateTime},
end
"""
This is the binary crossover function, implemented as detailed in Englander.
It chooses the first n and last m phases from
"""
function crossover()
end

View File

@@ -16,8 +16,8 @@ function get_true_max(mat::Vector{Array{Float64,2}})
return maximum(abs.(flatten(mat)))
end
function plot_orbits(paths::Vector{Vector{Vector{Float64}}};
primary::String="Earth",
function plot_orbits(paths::Vector{Vector{Vector{Float64}}},
primary::Body=Sun;
labels::Vector{String}=Vector{String}(),
title::String="Spacecraft Position",
colors::Vector{String}=Vector{String}())
@@ -28,7 +28,7 @@ function plot_orbits(paths::Vector{Vector{Vector{Float64}}};
x = cos.(θ) * sin.(ϕ)'
y = sin.(θ) * sin.(ϕ)'
z = repeat(cos.(ϕ)',outer=[N, 1])
ps = rs[primary] .* (x,y,z)
ps = primary.r .* (x,y,z)
x_p,y_p,z_p = ps
t1 = []
@@ -52,7 +52,7 @@ function plot_orbits(paths::Vector{Vector{Vector{Float64}}};
y=(y_p),
z=(z_p),
showscale=false,
colorscale = p_colors[primary])
colorscale = primary.color)
layout = Layout(title=title,
width=1000,

View File

@@ -1,4 +1,5 @@
export Sc
export Sc, test_sc, bepi, no_thrust
mutable struct Sc
dry_mass::Float64
mass_flow_rate::Float64
@@ -7,18 +8,6 @@ mutable struct Sc
duty_cycle::Float64
end
function Sc(name::String)
# This has extra thrusters to make plots more visible (and most don't use fuel)
if name == "test"
return Sc(9000., 0.00025/(2000*0.00981), 0.00025, 50, 0.9)
# This is the normal one
elseif name == "bepi"
return Sc(9000., 2*0.00025/(2000*0.00981), 0.00025, 2, 0.9)
elseif name == "no_thrust"
return Sc(9000., 0.01, 0., 0, 0.)
else
throw(ErrorException("Bad sc name"))
end
end
Base.copy(s::Sc) = Sc(s.dry_mass, s.mass_flow_rate, s.max_thrust, s.num_thrusters, s.duty_cycle)
const test_sc = Sc(8000., 0.00025/(2000*0.00981), 0.00025, 50, 0.9)
const bepi = Sc(8000., 2*0.00025/(2000*0.00981), 0.00025, 2, 0.9)
const no_thrust = Sc(8000., 0.01, 0., 0, 0.)