TP Transfert Orbital - Affichage
Vous pouvez remplacer la cellule contenant la fonction d'affichage par le code suivant (pour afficher l'animation il faut appeler la fonction plot_solution
avec l'argumetn animation=true
) :
# Corps du satellite
@userplot SatBody
@recipe function f(cp::SatBody)
x, y = cp.args
seriestype --> :shape
fillcolor --> :goldenrod1
linecolor --> :goldenrod1
x, y
end
@userplot SatLine
@recipe function f(cp::SatLine)
x, y = cp.args
linecolor --> :black
x, y
end
# Bras du satellite
@userplot SatBras
@recipe function f(cp::SatBras)
x, y = cp.args
seriestype --> :shape
fillcolor --> :goldenrod1
linecolor --> :white
x, y
end
# panneau
@userplot PanBody
@recipe function f(cp::PanBody)
x, y = cp.args
seriestype --> :shape
fillcolor --> :dodgerblue4
linecolor --> :black
x, y
end
# flamme
@userplot Flamme
@recipe function f(cp::Flamme)
x, y = cp.args
seriestype --> :shape
fillcolor --> :darkorange1
linecolor --> false
x, y
end
function satellite!(pl; position=[0;0], scale=1, rotate=0)
# Fonctions utiles
R(θ) = [ cos(θ) -sin(θ)
sin(θ) cos(θ)] # rotation
T(x, v) = x.+v # translation
H(λ, x) = λ.*x # homotéthie
SA(x, θ, c) = c .+ 2*[cos(θ);sin(θ)]'*(x.-c).*[cos(θ);sin(θ)]-(x.-c) # symétrie axiale
#
O = position
Rθ = R(rotate)
# Paramètres
α = π/10.0 # angle du tube, par rapport à l'horizontal
β = π/2-π/40 # angle des bras, par rapport à l'horizontal
Lp = scale.*1.4*cos(α) # longueur d'un panneau
lp = scale.*2.6*sin(α) # largueur d'un panneau
# Param bras du satellite
lb = scale.*3*cos(α) # longueur des bras
eb = scale.*cos(α)/30 # demi largeur des bras
xb = 0.0
yb = scale.*sin(α)
# Paramètres corps du satellite
t = range(-α, α, length = 50)
x = scale.*cos.(t);
y = scale.*sin.(t);
Δx = scale.*cos(α)
Δy = scale.*sin(α)
# Paramètres flamme
hF = yb # petite hauteur
HF = 2*hF # grande hauteur
LF = 5.5*Δx # longueur
# Dessin bras du satellite
M = T(Rθ*[ [xb-eb, xb+eb, xb+eb+lb*cos(β), xb-eb+lb*cos(β), xb-eb]';
[yb, yb, yb+lb*sin(β), yb+lb*sin(β), yb]'], O);
satbras!(pl, M[1,:], M[2,:])
M = T(Rθ*[ [xb-eb, xb+eb, xb+eb+lb*cos(β), xb-eb+lb*cos(β), xb-eb]';
-[yb, yb, yb+lb*sin(β), yb+lb*sin(β), yb']'], O);
satbras!(pl, M[1,:], M[2,:])
# Dessin flamme
M = [0 -LF -LF 0
-hF/2 -HF/2 HF/2 hF/2]
M = T(Rθ*T(M, [-Δx; 0]), O)
flamme!(pl, M[1, :], M[2, :])
# Dessin corps du satellite
M = T(Rθ*[ x'; y'], O); satbody!(pl, M[1,:], M[2,:]) # bord droit
M = T(Rθ*[-x'; y'], O); satbody!(pl, M[1,:], M[2,:]) # bord gauche
M = T(Rθ*[scale.*[-cos(α), cos(α), cos(α), -cos(α)]'
scale.*[-sin(α), -sin(α), sin(α), sin(α)]'], O);
satbody!(pl, M[1,:], M[2,:]) # interieur
M = T(Rθ*[ x'; y'], O); satline!(pl, M[1,:], M[2,:]) # bord droit
M = T(Rθ*[-x'; y'], O); satline!(pl, M[1,:], M[2,:]) # bord gauche
M = T(Rθ*[ x'.-2*Δx; y'], O); satline!(pl, M[1,:], M[2,:]) # bord gauche (droite)
M = T(Rθ*[ scale.*[-cos(α), cos(α)]';
scale.*[ sin(α), sin(α)]'], O);
satline!(pl, M[1,:], M[2,:]) # haut
M = T(Rθ*[ scale.*[-cos(α), cos(α)]';
-scale.*[ sin(α), sin(α)]'], O);
satline!(pl, M[1,:], M[2,:]) # bas
# Panneau
ep = (lb-3*lp)/6
panneau = [0 Lp Lp 0
0 0 lp lp]
ey = 3*eb # eloignement des panneaux au bras
vy = [cos(β-π/2); sin(β-π/2)] .* ey
v0 = [0; yb]
v1 = 2*ep*[cos(β); sin(β)]
v2 = (3*ep+lp)*[cos(β); sin(β)]
v3 = (4*ep+2*lp)*[cos(β); sin(β)]
pa1 = T(R(β-π/2)*panneau, v0+v1+vy); pa = T(Rθ*pa1, O); panbody!(pl, pa[1,:], pa[2,:])
pa2 = T(R(β-π/2)*panneau, v0+v2+vy); pa = T(Rθ*pa2, O); panbody!(pl, pa[1,:], pa[2,:])
pa3 = T(R(β-π/2)*panneau, v0+v3+vy); pa = T(Rθ*pa3, O); panbody!(pl, pa[1,:], pa[2,:])
pa4 = SA(pa1, β, [xb; yb]); pa = T(Rθ*pa4, O); panbody!(pl, pa[1,:], pa[2,:])
pa5 = SA(pa2, β, [xb; yb]); pa = T(Rθ*pa5, O); panbody!(pl, pa[1,:], pa[2,:])
pa6 = SA(pa3, β, [xb; yb]); pa = T(Rθ*pa6, O); panbody!(pl, pa[1,:], pa[2,:])
pa7 = SA(pa1, 0, [0; 0]); pa = T(Rθ*pa7, O); panbody!(pl, pa[1,:], pa[2,:])
pa8 = SA(pa2, 0, [0; 0]); pa = T(Rθ*pa8, O); panbody!(pl, pa[1,:], pa[2,:])
pa9 = SA(pa3, 0, [0; 0]); pa = T(Rθ*pa9, O); panbody!(pl, pa[1,:], pa[2,:])
pa10 = SA(pa7, -β, [xb; -yb]); pa = T(Rθ*pa10, O); panbody!(pl, pa[1,:], pa[2,:])
pa11 = SA(pa8, -β, [xb; -yb]); pa = T(Rθ*pa11, O); panbody!(pl, pa[1,:], pa[2,:])
pa12 = SA(pa9, -β, [xb; -yb]); pa = T(Rθ*pa12, O); panbody!(pl, pa[1,:], pa[2,:])
end;
@userplot TrajectoryPlot
@recipe function f(cp::TrajectoryPlot)
t, x, y, tf = cp.args
n = argmin(abs.(t.-tf))
inds = 1:n
seriescolor --> :white
linewidth --> range(0, 5, length = n)
seriesalpha --> range(0, 1, length = n)
aspect_ratio --> 1
label --> false
x[inds], y[inds]
end
# Fonction d'affichage d'une solution
function plot_solution(p0, tf; animation=false, fps=10, nFrame=200)
# On trace l'orbite de départ et d'arrivée
gr(dpi=300, size=(500,400), thickness_scaling=1)
r0 = norm(x0[1:2])
v0 = norm(x0[3:4])
a = 1.0/(2.0/r0-v0*v0/μ)
t1 = r0*v0*v0/μ - 1.0;
t2 = (x0[1:2]'*x0[3:4])/sqrt(a*μ);
e_ellipse = norm([t1 t2])
p_orb = a*(1-e_ellipse^2);
n_theta = 151
Theta = range(0, stop=2*pi, length=n_theta)
X1_orb_init = zeros(n_theta)
X2_orb_init = zeros(n_theta)
X1_orb_arr = zeros(n_theta)
X2_orb_arr = zeros(n_theta)
# Orbite initiale
for i in 1:n_theta
theta = Theta[i]
r_orb = p_orb/(1+e_ellipse*cos(theta));
X1_orb_init[i] = r_orb*cos(theta);
X2_orb_init[i] = r_orb*sin(theta);
end
# Orbite d'arrivée
for i in 1:n_theta
theta = Theta[i]
X1_orb_arr[i] = rf*cos(theta) ;
X2_orb_arr[i] = rf*sin(theta);
end;
# Calcul de la trajectoire
ode_sol = f((t0, tf), x0, p0)
t = ode_sol.t
n = size(t, 1)
x1 = [ode_sol[1, j] for j in 1:n ]
x2 = [ode_sol[2, j] for j in 1:n ]
v1 = [ode_sol[3, j] for j in 1:n ]
v2 = [ode_sol[4, j] for j in 1:n ]
u = zeros(2, length(t))
for j in 1:size(t, 1)
u[:,j] = control(ode_sol[5:8, j])
end
ee = 0.2
xmin = minimum([x1; X1_orb_init; X1_orb_arr]); xmin = xmin - ee * abs(xmin);
xmax = maximum([x1; X1_orb_init; X1_orb_arr]); xmax = xmax + ee * abs(xmax);
ymin = minimum([x2; X2_orb_init; X2_orb_arr]); ymin = ymin - ee * abs(ymin);
ymax = maximum([x2; X2_orb_init; X2_orb_arr]); ymax = ymax + ee * abs(ymax);
if animation
nFrame = min(nFrame, n);
anim = @animate for i ∈ 1:nFrame
px = plot(background_color=:gray26, xlims=(xmin, xmax), ylims=(ymin, ymax),
legend = false, framestyle = :none, aspect_ratio=:equal)
plot!(px, X1_orb_init, X2_orb_init, color = :olivedrab1, linewidth=2)
plot!(px, X1_orb_arr, X2_orb_arr, color = :turquoise1, linewidth=2)
plot!(px, [0.0], [0.0], color = :steelblue2, seriestype=:scatter, markersize = 20, markerstrokewidth=0)
trajectoryplot!(px, t, x1, x2, i*tf/nFrame)
tcur = i*tf/nFrame
indf = argmin(abs.(t.-tcur))
plot!(px, [x1[1]], [x2[1]], seriestype=:scatter, color = :white, markerstrokewidth=0)
satellite!(px, position=[x1[indf];x2[indf]], scale=2000, rotate=atan(u[2, indf], u[1, indf]))
end
gif(anim, "anim.gif", fps=fps);
else
pu1 = plot(t, u[1,:], color = :red, xlabel = "t", ylabel = "u1", legend = false)
pu2 = plot(t, u[2,:], color = :red, xlabel = "t", ylabel = "u2", legend = false)
px = plot(background_color=:gray26, xlims=(xmin, xmax), ylims=(ymin, ymax),
legend = false, framestyle = :none, aspect_ratio=:equal)
plot!(px, X1_orb_init, X2_orb_init, color = :olivedrab1, linewidth=2)
plot!(px, X1_orb_arr, X2_orb_arr, color = :turquoise1, linewidth=2)
plot!(px, [0.0], [0.0], color = :steelblue2, seriestype=:scatter, markersize = 20, markerstrokewidth=0)
plot!(px, [x1[1]], [x2[1]], seriestype=:scatter, color = :white, markerstrokewidth=0)
plot!(px, x1, x2, color = :white, size = (800, 800), linewidth = range(0, 5, length = n))
satellite!(px, position=[x1[end];x2[end]], scale=2000, rotate=atan(u[2, end], u[1, end]))
display(plot(pu1, pu2, layout = (1,2), size = (800,400)))
display(px)
end
end;