Skip to content
This repository has been archived by the owner on Feb 26, 2024. It is now read-only.

Commit

Permalink
Start the great cleanup.
Browse files Browse the repository at this point in the history
  • Loading branch information
helgee committed Feb 25, 2017
1 parent 16e8041 commit c8800f6
Show file tree
Hide file tree
Showing 4 changed files with 54 additions and 91 deletions.
2 changes: 0 additions & 2 deletions REQUIRE
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,5 @@ Dopri
ERFA
JPLEphemeris
NLopt
NLsolve
Parameters
PyPlot
Roots
33 changes: 0 additions & 33 deletions src/mission/solvers.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,5 @@
using NLopt

import BlackBoxOptim: bboptimize, best_candidate
import NLsolve: nlsolve
import NLopt: optimize

Expand All @@ -23,7 +22,6 @@ type NLoptSolver <: Solver
dx::Float64
end

type BBOptim <: Solver end
type NLsolveSolver <: Solver end

function NLoptSolver(;
Expand All @@ -34,37 +32,6 @@ function NLoptSolver(;
NLoptSolver(algorithm, differences, dx)
end

function minimize(mission, objective, sol::BBOptim)
output = deepcopy(mission)
cons = AbstractConstraint[objective]
#= cons = AbstractConstraint[] =#
for element in fieldnames(mission.stop)
if !isnull(getfield(mission.stop, element))
val = get(getfield(mission.stop, element))
con = KEPLERIAN_CONSTRAINTS[element]
push!(cons, con(val))
end
end
ranges = [(lower, upper) for (lower, upper) in zip(lowerbounds(output), upperbounds(output))]
res = bboptimize(x -> fitness(x, output, cons), SearchRange=ranges,
method=:separable_nes)
init = best_candidate(res)
setparameters!(output, init)
minimize(output, objective, NLoptSolver())
end

function fitness(x, mission, cons)
setparameters!(mission, x)
res = propagate(mission)
values = Float64[]
for con in cons
push!(values, evaluate(con, res))
#= push!(values, scale(evaluate(con, res))) =#
#= println(scale(evaluate(con, res))) =#
end
return sumabs2(values)
end

function minimize(mission, objective, sol::NLsolveSolver)
output = deepcopy(mission)
#= cons = AbstractConstraint[objective] =#
Expand Down
32 changes: 15 additions & 17 deletions src/trajectories.jl
Original file line number Diff line number Diff line change
@@ -1,10 +1,8 @@
using Dierckx
using PlotlyJS

import Base: getindex, endof, show
import PlotlyJS: plot

export Trajectory, plot
export Trajectory

abstract AbstractTrajectory

Expand Down Expand Up @@ -92,17 +90,17 @@ end

endof(tra::Trajectory) = tra.t[end]

function plot{F<:Frame, T<:Timescale, C<:CelestialBody}(tra::Trajectory{F,T,C})
re = equatorial_radius(constants(C))
rp = polar_radius(constants(C))
n = 100
θ = linspace(-π/2, π/2, n)
ϕ = linspace(0, 2π, n)
x = [re * cos(i) * cos(j) for i in θ, j in ϕ]
y = [re * cos(i) * sin(j) for i in θ, j in ϕ]
z = [rp * sin(i) for i in θ, j in ϕ];
s = surface(x=x, y=y, z=z, colorscale="Blues", showlegend=false)
times = linspace(tra.t[1], tra.t[end], length(tra.t)*100)
p = scatter3d(;x=evaluate(tra.xspl, times), y=evaluate(tra.yspl, times), z=evaluate(tra.zspl, times), mode="lines", line=attr(color="rgb(255,0,0)"))
plot([s, p])
end
#= function plot{F<:Frame, T<:Timescale, C<:CelestialBody}(tra::Trajectory{F,T,C}) =#
#= re = equatorial_radius(constants(C)) =#
#= rp = polar_radius(constants(C)) =#
#= n = 100 =#
#= θ = linspace(-π/2, π/2, n) =#
#= ϕ = linspace(0, 2π, n) =#
#= x = [re * cos(i) * cos(j) for i in θ, j in ϕ] =#
#= y = [re * cos(i) * sin(j) for i in θ, j in ϕ] =#
#= z = [rp * sin(i) for i in θ, j in ϕ]; =#
#= s = surface(x=x, y=y, z=z, colorscale="Blues", showlegend=false) =#
#= times = linspace(tra.t[1], tra.t[end], length(tra.t)*100) =#
#= p = scatter3d(;x=evaluate(tra.xspl, times), y=evaluate(tra.yspl, times), z=evaluate(tra.zspl, times), mode="lines", line=attr(color="rgb(255,0,0)")) =#
#= plot([s, p]) =#
#= end =#
78 changes: 39 additions & 39 deletions test/propagators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -24,43 +24,43 @@
tra = trajectory(s0, Kepler())
@test tra[end] State(ep+EpochDelta(seconds=period(s0)), r0, v0)
end
#= @testset "ODE Propagator" begin =#
#= r0 = [1131.340, -2282.343, 6672.423] =#
#= v0 = [-5.64305, 4.30333, 2.42879] =#
#= Δt = 40*60 =#
#= Δe = EpochDelta(seconds=Δt) =#
#= ep = Epoch(TT, now()) =#
#= s0 = State(ep, r0, v0) =#
#= r1 = [-4219.7527, 4363.0292, -3958.7666] =#
#= v1 = [3.689866, -1.916735, -6.112511] =#
#= ode = ODE(maxstep=10) =#
#= s2 = state(s0, Δt/4, ode) =#
#= s1 = state(s0, Δt, ode) =#
#= @test s1.rv ≈ [r1; v1] =#
#= s1 = state(s0, Δe, ode) =#
#= @test s1.rv ≈ [r1; v1] =#
#= tra = trajectory(s0, Δe, ode) =#
#= s1 = tra.trajectory[ep+Δe] =#
#= @test s1.rv ≈ [r1; v1] =#
#= @test tra.trajectory[Δt/4] ≈ s2 =#
#= ode = ODE(discontinuities=[Discontinuity(PericenterEvent(), Stop())], maxstep=100) =#
#= s1 = state(s0, period(s0), ode) =#
#= ano = trueanomaly(s1) =#
#= ano = ano > π ? abs(ano - 2π) : ano =#
#= @test round(ano, 10) ≈ 0.0 =#
#= ode = ODE(discontinuities=[Discontinuity(ApocenterEvent(), Stop())], maxstep=100) =#
#= s1 = state(s0, period(s0), ode) =#
#= ano = trueanomaly(s1) =#
#= @test ano ≈ π =#
#= ode = ODE(discontinuities=[Discontinuity(ApocenterEvent(), Abort())], maxstep=100) =#
#= @test_throws PropagatorAbort state(s0, period(s0), ode) =#
#= ode = ODE(discontinuities=[Discontinuity(StartEvent(), ImpulsiveManeuver(along=-1))]) =#
#= @test_throws PropagatorAbort state(s0, period(s0), ode) =#
#= ode = ODE(events=[ApocenterEvent()], maxstep=10) =#
#= tra = trajectory(s0, period(s0)*4, ode) =#
#= for t in tra.events[:index] =#
#= ano = trueanomaly(tra.trajectory[t]) =#
#= @test ano ≈ π =#
#= end =#
#= end =#
@testset "ODE Propagator" begin
r0 = [1131.340, -2282.343, 6672.423]
v0 = [-5.64305, 4.30333, 2.42879]
Δt = 40*60
Δe = EpochDelta(seconds=Δt)
ep = Epoch(TT, 2000, 1, 1)
s0 = State(ep, r0, v0)
r1 = [-4219.7527, 4363.0292, -3958.7666]
v1 = [3.689866, -1.916735, -6.112511]
ode = ODE(maxstep=10)
s2 = state(s0, Δt/4, ode)
s1 = state(s0, Δt, ode)
@test s1.rv [r1; v1]
s1 = state(s0, Δe, ode)
@test s1.rv [r1; v1]
tra = trajectory(s0, Δe, ode)
s1 = tra.trajectory[ep+Δe]
@test s1.rv [r1; v1]
@test tra.trajectory[Δt/4] s2
ode = ODE(discontinuities=[Discontinuity(PericenterEvent(), Stop())], maxstep=100)
s1 = state(s0, period(s0), ode)
ano = trueanomaly(s1)
ano = ano > π ? abs(ano - 2π) : ano
@test round(ano, 10) 0.0
ode = ODE(discontinuities=[Discontinuity(ApocenterEvent(), Stop())], maxstep=100)
s1 = state(s0, period(s0), ode)
ano = trueanomaly(s1)
@test ano π
ode = ODE(discontinuities=[Discontinuity(ApocenterEvent(), Abort())], maxstep=100)
@test_throws PropagatorAbort state(s0, period(s0), ode)
ode = ODE(discontinuities=[Discontinuity(StartEvent(), ImpulsiveManeuver(along=-1))])
@test_throws PropagatorAbort state(s0, period(s0), ode)
ode = ODE(events=[ApocenterEvent()], maxstep=10)
tra = trajectory(s0, period(s0)*4, ode)
for t in tra.events[:index]
ano = trueanomaly(tra.trajectory[t])
@test ano π
end
end
end

0 comments on commit c8800f6

Please sign in to comment.