diff --git a/SPICE/de430.bsp b/SPICE/de430.bsp new file mode 100644 index 0000000..fac78c6 Binary files /dev/null and b/SPICE/de430.bsp differ diff --git a/SPICE/naif0012.tls b/SPICE/naif0012.tls new file mode 100644 index 0000000..e1afdee --- /dev/null +++ b/SPICE/naif0012.tls @@ -0,0 +1,152 @@ +KPL/LSK + + +LEAPSECONDS KERNEL FILE +=========================================================================== + +Modifications: +-------------- + +2016, Jul. 14 NJB Modified file to account for the leapsecond that + will occur on December 31, 2016. + +2015, Jan. 5 NJB Modified file to account for the leapsecond that + will occur on June 30, 2015. + +2012, Jan. 5 NJB Modified file to account for the leapsecond that + will occur on June 30, 2012. + +2008, Jul. 7 NJB Modified file to account for the leapsecond that + will occur on December 31, 2008. + +2005, Aug. 3 NJB Modified file to account for the leapsecond that + will occur on December 31, 2005. + +1998, Jul 17 WLT Modified file to account for the leapsecond that + will occur on December 31, 1998. + +1997, Feb 22 WLT Modified file to account for the leapsecond that + will occur on June 30, 1997. + +1995, Dec 14 KSZ Corrected date of last leapsecond from 1-1-95 + to 1-1-96. + +1995, Oct 25 WLT Modified file to account for the leapsecond that + will occur on Dec 31, 1995. + +1994, Jun 16 WLT Modified file to account for the leapsecond on + June 30, 1994. + +1993, Feb. 22 CHA Modified file to account for the leapsecond on + June 30, 1993. + +1992, Mar. 6 HAN Modified file to account for the leapsecond on + June 30, 1992. + +1990, Oct. 8 HAN Modified file to account for the leapsecond on + Dec. 31, 1990. + + +Explanation: +------------ + +The contents of this file are used by the routine DELTET to compute the +time difference + +[1] DELTA_ET = ET - UTC + +the increment to be applied to UTC to give ET. + +The difference between UTC and TAI, + +[2] DELTA_AT = TAI - UTC + +is always an integral number of seconds. The value of DELTA_AT was 10 +seconds in January 1972, and increases by one each time a leap second +is declared. Combining [1] and [2] gives + +[3] DELTA_ET = ET - (TAI - DELTA_AT) + + = (ET - TAI) + DELTA_AT + +The difference (ET - TAI) is periodic, and is given by + +[4] ET - TAI = DELTA_T_A + K sin E + +where DELTA_T_A and K are constant, and E is the eccentric anomaly of the +heliocentric orbit of the Earth-Moon barycenter. Equation [4], which ignores +small-period fluctuations, is accurate to about 0.000030 seconds. + +The eccentric anomaly E is given by + +[5] E = M + EB sin M + +where M is the mean anomaly, which in turn is given by + +[6] M = M + M t + 0 1 + +where t is the number of ephemeris seconds past J2000. + +Thus, in order to compute DELTA_ET, the following items are necessary. + + DELTA_TA + K + EB + M0 + M1 + DELTA_AT after each leap second. + +The numbers, and the formulation, are taken from the following sources. + + 1) Moyer, T.D., Transformation from Proper Time on Earth to + Coordinate Time in Solar System Barycentric Space-Time Frame + of Reference, Parts 1 and 2, Celestial Mechanics 23 (1981), + 33-56 and 57-68. + + 2) Moyer, T.D., Effects of Conversion to the J2000 Astronomical + Reference System on Algorithms for Computing Time Differences + and Clock Rates, JPL IOM 314.5--942, 1 October 1985. + +The variable names used above are consistent with those used in the +Astronomical Almanac. + +\begindata + +DELTET/DELTA_T_A = 32.184 +DELTET/K = 1.657D-3 +DELTET/EB = 1.671D-2 +DELTET/M = ( 6.239996D0 1.99096871D-7 ) + +DELTET/DELTA_AT = ( 10, @1972-JAN-1 + 11, @1972-JUL-1 + 12, @1973-JAN-1 + 13, @1974-JAN-1 + 14, @1975-JAN-1 + 15, @1976-JAN-1 + 16, @1977-JAN-1 + 17, @1978-JAN-1 + 18, @1979-JAN-1 + 19, @1980-JAN-1 + 20, @1981-JUL-1 + 21, @1982-JUL-1 + 22, @1983-JUL-1 + 23, @1985-JUL-1 + 24, @1988-JAN-1 + 25, @1990-JAN-1 + 26, @1991-JAN-1 + 27, @1992-JUL-1 + 28, @1993-JUL-1 + 29, @1994-JUL-1 + 30, @1996-JAN-1 + 31, @1997-JUL-1 + 32, @1999-JAN-1 + 33, @2006-JAN-1 + 34, @2009-JAN-1 + 35, @2012-JUL-1 + 36, @2015-JUL-1 + 37, @2017-JAN-1 ) + +\begintext + + diff --git a/julia/Manifest.toml b/julia/Manifest.toml index 0caedc8..86b56f6 100644 --- a/julia/Manifest.toml +++ b/julia/Manifest.toml @@ -33,6 +33,12 @@ git-tree-sha1 = "08d0b679fd7caa49e2bca9214b131289e19808c0" uuid = "ad839575-38b3-5650-b840-f874b8c74a25" version = "0.12.5" +[[CSPICE_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "56fbc7d843c24bb3a3689e81b064aa730f72db8c" +uuid = "07f52509-e9d9-513c-a20d-3b911885bf96" +version = "66.1.0+0" + [[ChainRulesCore]] deps = ["Compat", "LinearAlgebra", "SparseArrays"] git-tree-sha1 = "bdc0937269321858ab2a4f288486cb258b9a0af7" @@ -394,6 +400,12 @@ version = "1.1.3" [[SHA]] uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" +[[SPICE]] +deps = ["CSPICE_jll", "LinearAlgebra"] +git-tree-sha1 = "dcbc85c97ea75057e96547a189d819b5e80af86a" +uuid = "5bab7191-041a-5c2e-a744-024b9c3a5062" +version = "0.2.2" + [[Serialization]] uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" diff --git a/julia/Project.toml b/julia/Project.toml index 3e6ef18..611de22 100644 --- a/julia/Project.toml +++ b/julia/Project.toml @@ -11,4 +11,5 @@ NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" PlotlyBase = "a03496cd-edff-5a9b-9e67-9cda94a718b5" PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SPICE = "5bab7191-041a-5c2e-a744-024b9c3a5062" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/julia/src/Thesis.jl b/julia/src/Thesis.jl index 88a2fd4..f510c67 100644 --- a/julia/src/Thesis.jl +++ b/julia/src/Thesis.jl @@ -1,6 +1,9 @@ module Thesis - using LinearAlgebra, ForwardDiff, PlotlyJS + using LinearAlgebra, ForwardDiff, PlotlyJS, SPICE + + furnsh("../../SPICE/naif0012.tls") + furnsh("../../SPICE/de430.bsp") include("./constants.jl") include("./spacecraft.jl") @@ -12,4 +15,5 @@ module Thesis include("./inner_loop/monotonic_basin_hopping.jl") include("./inner_loop/phase.jl") include("./inner_loop/inner_loop.jl") + end \ No newline at end of file diff --git a/julia/src/constants.jl b/julia/src/constants.jl index cf459a3..a973363 100644 --- a/julia/src/constants.jl +++ b/julia/src/constants.jl @@ -32,13 +32,13 @@ export μs, G, GMs, μ, rs, as, es, AU return μ(GMs[primary]/G, GMs[secondary]/G) end - GMs = Dict( + const GMs = Dict( "Sun" => 132712440041.93938, "Earth" => 398600.435436, "Moon" => 4902.800066) # Radii - rs = Dict( + const rs = Dict( "Sun" => 696000., "Mercury" => 2439., "Venus" => 6052., @@ -52,7 +52,7 @@ export μs, G, GMs, μ, rs, as, es, AU "Pluto" => 1151.) # Semi Major Axes - as = Dict( + const as = Dict( "Mercury" => 57909083., "Venus" => 108208601., "Earth" => 149598023., @@ -65,12 +65,12 @@ export μs, G, GMs, μ, rs, as, es, AU "Pluto" => 5915799000.) # Eccentricities - es = Dict( + const es = Dict( "Earth" => 0.016708617, "Moon" => 0.0549) # J2 for basic oblateness - j2s = Dict( + const j2s = Dict( "Mercury" => 0.00006, "Venus" => 0.000027, "Earth" => 0.0010826269, @@ -83,7 +83,7 @@ export μs, G, GMs, μ, rs, as, es, AU "Pluto" => 0.) # These are just the colors for plots - p_colors = Dict( + const p_colors = Dict( "Sun" => :Electric, "Mercury" => :heat, "Venus" => :turbid, @@ -96,5 +96,19 @@ export μs, G, GMs, μ, rs, as, es, AU "Neptune" => :ice, "Pluto" => :matter) - AU = 149597870.691 #km - init_STM = vec(Matrix{Float64}(I,6,6)) + 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)) diff --git a/julia/src/inner_loop/inner_loop.jl b/julia/src/inner_loop/inner_loop.jl index d4b7012..7297a7c 100644 --- a/julia/src/inner_loop/inner_loop.jl +++ b/julia/src/inner_loop/inner_loop.jl @@ -7,19 +7,50 @@ This is it. The outer function call for the inner loop. After this is done, there's only the outer loop left to do. And that's pretty easy. """ function inner_loop(launch_date::DateTime, - RLA::Float64, - DLA::Float64, - phases::Vector{Phase}) + phases::Vector{Phase}; + min_flyby::Float64=1000.) # First we need to do some quick checks that the mission is well formed for i in 1:length(phases) if i == 1 @assert phases[i].from_planet == "Earth" else - @assert phases[i].from_planet == phases[i-1].to_planet - @assert phases[i].v∞_outgoing == phases[i-1].v∞_incoming + # Check that the planet is valid + if phases[i].from_planet ∉ keys(μs) + error("Planet is not valid: ", phases[i].from_planet) + # Check that there is only one flyby planet + elseif phases[i].from_planet != phases[i-1].to_planet + fromP, toP = phases[i].from_planet, phases[i-1].to_planet + error("Planets don't match up: (phase $(i)) $(fromP) / (phase $(i-1)) $(toP)") + # Check that v∞_in == v∞_out + elseif !isapprox(norm(phases[i].v∞_outgoing), norm(phases[i-1].v∞_incoming)) + norm_incoming = norm(phases[i-1].v∞_incoming) + norm_outgoing = norm(phases[i].v∞_outgoing) + error("""Norms of vectors aren't equal: + Phase $(i-1): $(phases[i-1].v∞_incoming) / norm: $(norm_incoming) + Phase $(i): $(phases[i].v∞_outgoing) / norm: $(norm_outgoing)""") + end + # Check that the approach is not too low + v∞ = norm(phases[i].v∞_outgoing) + δ = acos((phases[i].v∞_outgoing ⋅ phases[i-1].v∞_incoming)/v∞^2) + flyby = μs[phases[i].from_planet]/v∞^2 * (1/sin(δ/2) - 1) + true_min = rs[phases[i].from_planet] + min_flyby + if flyby <= true_min + error("Flyby was too low between phase $(i-1) and $(i): $(flyby) / $(true_min)") + end + end end + + time = utc2et(Dates.format(launch_date,"yyyy-mm-ddTHH:MM:SS")) + for phase in phases + planet1_state = spkssb(ids[phase.from_planet], time, "J2000") + println(planet1_state) + time += phase.time_of_flight + planet2_state = spkssb(ids[phase.to_planet], time, "J2000") + println(planet2_state) + end + - return RLA + DLA + phases[1].v∞_incoming + return phases[1].v∞_incoming[1] end \ No newline at end of file diff --git a/julia/src/inner_loop/phase.jl b/julia/src/inner_loop/phase.jl index 96a380d..e246a18 100644 --- a/julia/src/inner_loop/phase.jl +++ b/julia/src/inner_loop/phase.jl @@ -3,7 +3,7 @@ export Phase struct Phase from_planet::String to_planet::String - time_of_flight::Float64 # seconds - v∞_outgoing::Float64 # Km/s - v∞_incoming::Float64 # Km/s + time_of_flight::Float64 # seconds + v∞_outgoing::Vector{Float64} # Km/s + v∞_incoming::Vector{Float64} # Km/s end \ No newline at end of file diff --git a/julia/test/inner_loop/inner_loop.jl b/julia/test/inner_loop/inner_loop.jl index 3646f0d..07f976e 100644 --- a/julia/test/inner_loop/inner_loop.jl +++ b/julia/test/inner_loop/inner_loop.jl @@ -4,9 +4,18 @@ using Dates - phase1 = Phase("Earth", "Mars", 3600*24*365*1.5, 5., 2.) - phase2 = Phase("Mars", "Jupiter", 3600*24*365*3.5, 2., 0.1) - inner_loop(DateTime(2024,3,5), 0.3, 0.4, [phase1, phase2]) + phase1 = Phase("Earth", + "Mars", + 3600*24*365*1.5, + [1., 0.3, 0.], + [3., 3., 0.]) + + phase2 = Phase("Mars", + "Jupiter", + 3600*24*365*3.5, + [2., 3.7416573867739413, 0.], + [0.3, 1., 0.]) + inner_loop(DateTime(2024,3,5), [phase1, phase2]) @test true diff --git a/julia/test/runtests.jl b/julia/test/runtests.jl index 2eecafa..5a4eaa0 100644 --- a/julia/test/runtests.jl +++ b/julia/test/runtests.jl @@ -1,17 +1,20 @@ using Test using Random using LinearAlgebra +using SPICE using Thesis +furnsh("../../SPICE/naif0012.tls") +furnsh("../../SPICE/de430.bsp") # Tests @testset "All Tests" begin - include("spacecraft.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("spacecraft.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/inner_loop.jl") end diff --git a/readme.md b/readme.md index ab9102d..7a86851 100644 --- a/readme.md +++ b/readme.md @@ -23,9 +23,9 @@ To generate the PDFs (currently just the notes, but soon to include the thesis p In order to run the Julia code, you'll need Julia v1.6 or higher. This project is packaged as a package, so installing the dependencies (and the project itself) is simple. Just run: ```julia -julia> using Pkg -(Thesis) pkg> activate julia -(Thesis) pkg> build +using Pkg +Pkg.activate("julia") +Pkg.build() ``` ## Installation