From 7296df65bc6b6042a867f0077ac17a4528dec8fa Mon Sep 17 00:00:00 2001 From: rconnorjohnstone Date: Mon, 17 May 2021 21:32:25 -0600 Subject: [PATCH] Began! --- julia/constants.jl | 100 +++++++++++++++++++++++++++++++++++++++ julia/conversions.jl | 76 +++++++++++++++++++++++++++++ julia/laguerre-conway.jl | 57 ++++++++++++++++++++++ julia/main.jl | 6 +++ julia/sft.jl | 12 +++++ 5 files changed, 251 insertions(+) create mode 100644 julia/constants.jl create mode 100644 julia/conversions.jl create mode 100644 julia/laguerre-conway.jl create mode 100644 julia/main.jl create mode 100644 julia/sft.jl diff --git a/julia/constants.jl b/julia/constants.jl new file mode 100644 index 0000000..d410086 --- /dev/null +++ b/julia/constants.jl @@ -0,0 +1,100 @@ +# ----------------------------------------------------------------------------- +# DEFINING CONSTANTS +# ----------------------------------------------------------------------------- + +export μs, G, GMs, μ, rs, as, es + +# 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) + + G = 6.67430e-20 + + function μ(m1::Float64, m2::Float64) + return m2/(m1+m2) + end + + 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 + + GMs = Dict( + "Sun" => 132712440041.93938, + "Earth" => 398600.435436, + "Moon" => 4902.800066) + +# Radii + 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 + as = Dict( + "Mercury" => 57909083., + "Venus" => 108208601., + "Earth" => 149598023., + "Moon" => 384400., + "Mars" => 227939186., + "Jupiter" => 778298361., + "Saturn" => 1429394133., + "Uranus" => 2875038615., + "Neptune" => 4504449769., + "Pluto" => 5915799000.) + +# Eccentricities + es = Dict( + "Earth" => 0.016708617, + "Moon" => 0.0549) + +# J2 for basic oblateness + 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 + 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) + + AU = 149597870.691 #km + init_STM = vec(Matrix{Float64}(I,6,6)) diff --git a/julia/conversions.jl b/julia/conversions.jl new file mode 100644 index 0000000..f18ca2f --- /dev/null +++ b/julia/conversions.jl @@ -0,0 +1,76 @@ +function oe_to_rθh(oe::Vector,μ::Real) :: Vector + + a,e,i,Ω,ω,ν = oe + return [a*(1-e^2)/(1+e*cos(ν)), + 0, + 0, + (μ/sqrt(μ*a*(1-e^2)))*e*sin(ν), + (μ/sqrt(μ*a*(1-e^2)))*(1+e*cos(ν)), + 0] + +end + +function rθh_to_xyz(rθh_vec::Vector,oe::Vector) + + a,e,i,Ω,ω,ν = oe + θ = ω+ν + 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] + DCM = kron(Matrix(I,2,2),DCM) + return DCM*rθh_vec + +end + +function oe_to_xyz(oe::Vector,μ::Real) + + return rθh_to_xyz(oe_to_rθh(oe,μ),oe) + +end + +function xyz_to_oe(cart_vec::Vector,μ::Real) + + r_xyz, v_xyz = cart_vec[1:3],cart_vec[4:6] + r = norm(r_xyz) + h_xyz = cross(r_xyz,v_xyz) #km^2/s + h = norm(h_xyz) #km^2/s + ξ = dot(v_xyz,v_xyz)/2 - μ/r #km^2 s^-2 + a = -μ/(2ξ) #km + e_xyz = cross(v_xyz,h_xyz)/μ - r_xyz/r + 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") + 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 Ω") + 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)") + 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 ν") + end + Ω = dot(n_xyz,[0,1,0]) > 0. ? Ω : -Ω + ω = dot(e_xyz,[0,0,1]) > 0. ? ω : -ω + ν = dot(r_xyz,v_xyz) > 0. ? ν : -ν + return [a,e,i,Ω,ω,ν] + +end diff --git a/julia/laguerre-conway.jl b/julia/laguerre-conway.jl new file mode 100644 index 0000000..99782d2 --- /dev/null +++ b/julia/laguerre-conway.jl @@ -0,0 +1,57 @@ +function laguerre_conway(state::Vector{Float64}, μ::Float64, time::Float64) + + n = 5 # Choose LaGuerre-Conway "n" + i = 0 + + r0 = state[1:3] # Are we in elliptical or hyperbolic orbit? + r0_mag = norm(r0) + v0 = state[4:6] + v0_mag = norm(v0) + σ0 = (r0 ⋅ v0)/√(μ) + a = 1 / ( 2/r0_mag - v0_mag^2/μ ) + coeff = 1 - r0_mag/a + + if a > 0 # Elliptical + ΔM = ΔE_new = √(μ) / sqrt(a^3) * time + ΔE = 1000 + while abs(ΔE - ΔE_new) > 1e-12 + ΔE = ΔE_new + F = ΔE - ΔM + σ0 / √(a) * (1-cos(ΔE)) - coeff * sin(ΔE) + dF = 1 + σ0 / √(a) * sin(ΔE) - coeff * cos(ΔE) + d2F = σ0 / √(a) * cos(ΔE) + coeff * sin(ΔE) + 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 + end + F = 1 - a/r0_mag * (1-cos(ΔE)) + G = a * σ0/ √(μ) * (1-cos(ΔE)) + r0_mag * √(a) / √(μ) * sin(ΔE) + r = a + (r0_mag - a) * cos(ΔE) + σ0 * √(a) * sin(ΔE) + Ft = -√(a)*√(μ) / (r*r0_mag) * sin(ΔE) + Gt = 1 - a/r * (1-cos(ΔE)) + else # Hyperbolic or Parabolic + ΔN = √(μ) / sqrt(-a^3) * time + ΔH = 0 + ΔH_new = time < 0 ? -1 : 1 + while abs(ΔH - ΔH_new) > 1e-12 + ΔH = ΔH_new + F = -ΔN - ΔH + σ0 / √(-a) * (cos(ΔH)-1) + coeff * sin(ΔH) + dF = -1 + σ0 / √(-a) * sin(ΔH) + coeff * cos(ΔH) + d2F = σ0 / √(-a) * cos(ΔH) + coeff * sin(ΔH) + 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 + end + F = 1 - a/r0_mag * (1-cos(ΔH)) + G = a * σ0/ √(μ) * (1-cos(ΔH)) + r0_mag * √(-a) / √(μ) * sin(ΔH) + r = a + (r0_mag - a) * cos(ΔH) + σ0 * √(-a) * sin(ΔH) + Ft = -√(-a)*√(μ) / (r*r0_mag) * sin(ΔH) + Gt = 1 - a/r * (1-cos(ΔH)) + end + + if i > 100 + throw(ErrorException("LaGuerre-Conway did not converge!")) + end + + return [ F*state[1:3] + G*state[4:6]; Ft*state[1:3] + Gt*state[4:6]] + +end diff --git a/julia/main.jl b/julia/main.jl new file mode 100644 index 0000000..6bcec68 --- /dev/null +++ b/julia/main.jl @@ -0,0 +1,6 @@ +using LinearAlgebra + +include("constants.jl") +include("conversions.jl") +include("sft.jl") +include("laguerre-conway.jl") diff --git a/julia/sft.jl b/julia/sft.jl new file mode 100644 index 0000000..6ada456 --- /dev/null +++ b/julia/sft.jl @@ -0,0 +1,12 @@ +""" +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 +