diff --git a/tp/space.jl b/tp/space.jl
new file mode 100644
index 0000000000000000000000000000000000000000..1f9c72aaf5b65ac8ff8f09316865a8daefccd9ee
--- /dev/null
+++ b/tp/space.jl
@@ -0,0 +1,849 @@
+module Space
+
+ using Plots
+ using Plots.PlotMeasures
+ using SparseArrays
+ using Printf
+ using DifferentialEquations
+ using LinearAlgebra # norm
+
+ export animation
+
+ # --------------------------------------------------------------------------------------------
+ x0 = [-42272.67, 0, 0, -5796.72] # état initial
+ μ = 5.1658620912*1e12
+ rf = 42165.0
+ t0 = 0
+ global γ_max = nothing
+ global F_max = nothing
+
+ # Maximizing control
+ function control(p)
+ u = zeros(eltype(p),2)
+ u[1] = p[3]*γ_max/norm(p[3:4])
+ u[2] = p[4]*γ_max/norm(p[3:4])
+ return u
+ end
+
+ # --------------------------------------------------------------------------------------------
+ # Données pour la trajectoire durant le transfert
+ mutable struct Transfert
+ initial_adjoint
+ duration
+ flow
+ end
+
+ # données pour la trajectoire pré-transfert
+ mutable struct PreTransfert
+ initial_point
+ duration
+ end
+
+ # données pour la trajectoire post-transfert
+ mutable struct PostTransfert
+ duration
+ end
+
+ # --------------------------------------------------------------------------------------------
+ # Default options for flows
+ # --------------------------------------------------------------------------------------------
+ function __abstol()
+ return 1e-10
+ end
+
+ function __reltol()
+ return 1e-10
+ end
+
+ function __saveat()
+ return []
+ end
+
+ # --------------------------------------------------------------------------------------------
+ #
+ # Vector Field
+ #
+ # --------------------------------------------------------------------------------------------
+ struct VectorField f::Function end
+
+ function (vf::VectorField)(x) # https://docs.julialang.org/en/v1/manual/methods/#Function-like-objects
+ return vf.f(x)
+ end
+
+ # Flow of a vector field
+ function Flow(vf::VectorField)
+
+ function rhs!(dx, x, dummy, t)
+ dx[:] = vf(x)
+ end
+
+ function f(tspan, x0; abstol=__abstol(), reltol=__reltol(), saveat=__saveat())
+ ode = ODEProblem(rhs!, x0, tspan)
+ sol = solve(ode, Tsit5(), abstol=abstol, reltol=reltol, saveat=saveat)
+ return sol
+ end
+
+ function f(t0, x0, t; abstol=__abstol(), reltol=__reltol(), saveat=__saveat())
+ sol = f((t0, t), x0; abstol=abstol, reltol=reltol, saveat=saveat)
+ n = size(x0, 1)
+ return sol[1:n, end]
+ end
+
+ return f
+
+ end
+
+ # Kepler's law
+ function kepler(x)
+ n = size(x, 1)
+ dx = zeros(eltype(x), n)
+ x1 = x[1]
+ x2 = x[2]
+ x3 = x[3]
+ x4 = x[4]
+ dsat = norm(x[1:2]); r3 = dsat^3;
+ dx[1] = x3
+ dx[2] = x4
+ dx[3] = -μ*x1/r3
+ dx[4] = -μ*x2/r3
+ return dx
+ end
+
+ KeplerFlow = Flow(VectorField(kepler));
+
+ function animation(p0, tf, f, γ_max; fps=10, nFrame=100)
+ #
+ t_pre_transfert = 2*5.302395297743802
+ x_pre_transfert = x0
+ pre_transfert_data = PreTransfert(x_pre_transfert, t_pre_transfert)
+ #
+ transfert_data = Transfert(p0, tf, f)
+ #
+ t_post_transfert = 20
+ post_transfert_data = PostTransfert(t_post_transfert)
+ #
+ __animation(pre_transfert_data, transfert_data, post_transfert_data, KeplerFlow, γ_max, fps=fps, nFrame=nFrame)
+ end
+
+ function __animation(pre_transfert_data, transfert_data, post_transfert_data, f0, _γ_max; fps=10, nFrame=200)
+
+ #
+ global F_max = _γ_max*2000*10^3/3600^2
+ global γ_max = _γ_max
+
+ # pré-transfert
+ t_pre_transfert = pre_transfert_data.duration
+ x_pre_transfert = pre_transfert_data.initial_point
+
+ # transfert
+ p0_transfert = transfert_data.initial_adjoint
+ tf_transfert = transfert_data.duration
+ f = transfert_data.flow
+
+ # post-trasfert
+ t_post_transfert = post_transfert_data.duration
+
+ # On calcule les orbites initiale et finale
+ 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.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)
+
+ for i in 1:n_theta
+ theta = Theta[i]
+ r_orb = p_orb/(1+e_ellipse*cos(theta));
+ # Orbite initiale
+ X1_orb_init[i] = r_orb*cos(theta);
+ X2_orb_init[i] = r_orb*sin(theta);
+ # Orbite finale
+ X1_orb_arr[i] = rf*cos(theta) ;
+ X2_orb_arr[i] = rf*sin(theta);
+ end
+
+ # Centre de la fenêtre
+ cx = 40000
+ cy = -7000
+
+ # Taille de la fenêtre
+ w = 1600
+ h = 900
+ ρ = h/w
+
+ # Limites de la fenêtre
+ ee = 0.5
+ xmin = minimum([X1_orb_init; X1_orb_arr]); xmin = xmin - ee * abs(xmin);
+ xmax = maximum([X1_orb_init; X1_orb_arr]); xmax = xmax + ee * abs(xmax);
+ ymin = minimum([X2_orb_init; X2_orb_arr]); ymin = ymin - ee * abs(ymin);
+ ymax = maximum([X2_orb_init; X2_orb_arr]); ymax = ymax + ee * abs(ymax);
+
+ Δy = ymax-ymin
+ Δx = Δy/ρ
+ Δn = Δx - (xmax-xmin)
+ xmin = xmin - Δn/2
+ xmax = xmax + Δn/2
+
+ # Trajectoire pré-transfert
+ traj_pre_transfert = f0((0.0, t_pre_transfert), x_pre_transfert)
+ times_pre_transfert = traj_pre_transfert.t
+ n_pre_transfert = size(times_pre_transfert, 1)
+ x1_pre_transfert = [traj_pre_transfert[1, j] for j in 1:n_pre_transfert ]
+ x2_pre_transfert = [traj_pre_transfert[2, j] for j in 1:n_pre_transfert ]
+ v1_pre_transfert = [traj_pre_transfert[3, j] for j in 1:n_pre_transfert ]
+ v2_pre_transfert = [traj_pre_transfert[4, j] for j in 1:n_pre_transfert ]
+
+ # Trajectoire pendant le transfert
+ traj_transfert = f((t0, tf_transfert), x0, p0_transfert)
+ times_transfert = traj_transfert.t
+ n_transfert = size(times_transfert, 1)
+ x1_transfert = [traj_transfert[1, j] for j in 1:n_transfert ]
+ x2_transfert = [traj_transfert[2, j] for j in 1:n_transfert ]
+ v1_transfert = [traj_transfert[3, j] for j in 1:n_transfert ]
+ v2_transfert = [traj_transfert[4, j] for j in 1:n_transfert ]
+ u_transfert = zeros(2, length(times_transfert))
+ for j in 1:n_transfert
+ u_transfert[:,j] = control(traj_transfert[5:8, j])
+ end
+ u_transfert = u_transfert./γ_max # ||u|| ≤ 1
+
+ # post-transfert
+ x_post_transfert = traj_transfert[:,end]
+
+ # Trajectoire post-transfert
+ traj_post_transfert = f0((0.0, t_post_transfert), x_post_transfert)
+ times_post_transfert = traj_post_transfert.t
+ n_post_transfert = size(times_post_transfert, 1)
+ x1_post_transfert = [traj_post_transfert[1, j] for j in 1:n_post_transfert ]
+ x2_post_transfert = [traj_post_transfert[2, j] for j in 1:n_post_transfert ]
+ v1_post_transfert = [traj_post_transfert[3, j] for j in 1:n_post_transfert ]
+ v2_post_transfert = [traj_post_transfert[4, j] for j in 1:n_post_transfert ]
+
+ # Angles de rotation du satellite pendant le pré-transfert
+ # Et poussée normalisée entre 0 et 1
+ θ0 = atan(u_transfert[2, 1], u_transfert[1, 1])
+ θ_pre_transfert = range(π/2, mod(θ0, 2*π), length=n_pre_transfert)
+ F_pre_transfert = zeros(1, n_pre_transfert)
+
+ # Angles de rotation du satellite pendant le transfert
+ θ_transfert = atan.(u_transfert[2, :], u_transfert[1, :])
+ F_transfert = zeros(1, n_transfert)
+ for j in 1:n_transfert
+ F_transfert[j] = norm(u_transfert[:,j])
+ end
+
+ # Angles de rotation du satellite pendant le post-transfert
+ θ1 = atan(u_transfert[2, end], u_transfert[1, end])
+ θ2 = atan(-x2_post_transfert[end], -x1_post_transfert[end])
+ θ_post_transfert = range(mod(θ1, 2*π), mod(θ2, 2*π), length=n_post_transfert)
+ F_post_transfert = zeros(1, n_post_transfert)
+
+ # Etoiles
+ Δx = xmax-xmin
+ Δy = ymax-ymin
+ ρ = Δy/Δx
+ S = stars(ρ)
+
+ # nombre total d'images
+ nFrame = min(nFrame, n_pre_transfert+n_transfert+n_post_transfert);
+
+ # Pour l'affichage de la trajectoire globale
+ times = [times_pre_transfert[1:end-1];
+ times_pre_transfert[end].+times_transfert[1:end-1];
+ times_pre_transfert[end].+times_transfert[end].+times_post_transfert[1:end]]
+ x1 = [x1_pre_transfert[1:end-1]; x1_transfert[1:end-1]; x1_post_transfert[:]]
+ x2 = [x2_pre_transfert[1:end-1]; x2_transfert[1:end-1]; x2_post_transfert[:]]
+ v1 = [v1_pre_transfert[1:end-1]; v1_transfert[1:end-1]; v1_post_transfert[:]]
+ v2 = [v2_pre_transfert[1:end-1]; v2_transfert[1:end-1]; v2_post_transfert[:]]
+ θ = [ θ_pre_transfert[1:end-1]; θ_transfert[1:end-1]; θ_post_transfert[:]]
+ F = [ F_pre_transfert[1:end-1]; F_transfert[1:end-1]; F_post_transfert[:]]
+
+ # plot thrust on/off
+ th = [BitArray(zeros(n_pre_transfert-1));
+ BitArray(ones(n_transfert-1));
+ BitArray(zeros(n_post_transfert))]
+
+ # plot trajectory
+ pt = [BitArray(ones(n_pre_transfert-1));
+ BitArray(ones(n_transfert-1));
+ BitArray(zeros(n_post_transfert))]
+
+ # Contrôle sur la trajectoire totale
+ u_total = hcat([zeros(2, n_pre_transfert-1),
+ u_transfert[:, 1:n_transfert-1],
+ zeros(2, n_post_transfert)]...)
+
+ # temps total
+ temps_transfert_global = times[end]
+
+ # pas de temps pour le transfert global
+ if nFrame>1
+ Δt = temps_transfert_global/(nFrame-1)
+ else
+ Δt = 0.0
+ end
+
+ # opacités des orbites initiale et finale
+ op_initi = [range(0.0, 1.0, length=n_pre_transfert-1);
+ range(1.0, 0.0, length=n_transfert-1);
+ zeros(n_post_transfert)]
+ op_final = [zeros(n_pre_transfert-1);
+ range(0.0, 1.0, length=n_transfert-1);
+ range(1.0, 0.0, length=int(n_post_transfert/4));
+ zeros(n_post_transfert-int(n_post_transfert/4))]
+
+ println("Fps = ", fps)
+ println("NFrame = ", nFrame)
+ println("Vitesse = ", nFrame/(temps_transfert_global*fps))
+ println("Durée totale de la mission = ", temps_transfert_global, " h")
+
+ # animation
+ anim = @animate for i ∈ 1:nFrame
+
+ # Δt : pas de temps
+ # time_current : temps courant de la mission totale à l'itération i
+ # i_current : indice tel que times[i_current] = time_current
+ # w, h : width, height de la fenêtre
+ # xmin, xmax, ymin, ymax : limites des axes du plot principal
+ # X1_orb_init, X2_orb_init : coordonnées de l'orbite initiale
+ # X1_orb_arr, X2_orb_arr : coordonnées de l'orbite finale
+ # cx, cy : coordonées du centre de l'affichage du tranfert
+ # S : data pour les étoiles
+ # Δx, Δy : xmax-xmin, ymax-ymin
+ # times : tous les temps de la mission complète, ie pre-transfert, transfert et post-transfert
+ # x1, x2 : vecteur de positions du satellite
+ # θ : vecteur d'orientations du satellite
+ # th : vecteur de booléens - thrust on/off
+ # u_total : vecteur de contrôles pour toute la mission
+ # F_max, γ_max : poussée max
+ # subplot_current : valeur du subplot courant
+ # cam_x, cam_y : position de la caméra
+ # cam_zoom : zoom de la caméra
+
+ cam_x = cx
+ cam_y = cy
+ cam_zoom = 1
+
+ time_current = (i-1)*Δt
+ i_current = argmin(abs.(times.-time_current))
+
+ px = background(w, h, xmin, xmax, ymin, ymax,
+ X1_orb_init, X2_orb_init, X1_orb_arr, X2_orb_arr,
+ cx, cy, S, Δx, Δy, cam_x, cam_y, cam_zoom,
+ op_initi[i_current], op_final[i_current], times, time_current)
+
+ trajectoire!(px, times, x1, x2, θ, F, th, time_current, cx, cy, pt)
+
+ subplot_current = 2
+ subplot_current = panneau_control!(px, time_current, times, u_total,
+ F_max, subplot_current)
+
+ subplot_current = panneau_information!(px, subplot_current, time_current,
+ i_current, x1, x2, v1, v2, θ, F_max, tf_transfert, X1_orb_init,
+ X2_orb_init, X1_orb_arr, X2_orb_arr)
+
+ end
+
+ # enregistrement
+ #gif(anim, "transfert-temps-min-original.mp4", fps=fps);
+ gif(anim, "transfert-temps-min.gif", fps=fps);
+
+ end;
+
+ # --------------------------------------------------------------------------------------------
+ # Preliminaries
+ int(x) = floor(Int, x);
+ rayon_Terre = 6371;
+ sc_sat = 1000; # scaling for the satellite
+
+ # --------------------------------------------------------------------------------------------
+ #
+ # Satellite
+ #
+ # --------------------------------------------------------------------------------------------
+ # 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; subplot=1, thrust=false, position=[0;0], scale=1, rotate=0, thrust_val=1.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,:], subplot=subplot)
+ 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,:], subplot=subplot)
+
+ # Dessin corps du satellite
+ M = T(Rθ*[ x'; y'], O); satbody!(pl, M[1,:], M[2,:], subplot=subplot) # bord droit
+ M = T(Rθ*[-x'; y'], O); satbody!(pl, M[1,:], M[2,:], subplot=subplot) # bord gauche
+ M = T(Rθ*[scale.*[-cos(α), cos(α), cos(α), -cos(α)]'
+ scale.*[-sin(α), -sin(α), sin(α), sin(α)]'], O);
+ satbody!(pl, M[1,:], M[2,:], subplot=subplot) # interieur
+
+ M = T(Rθ*[ x'; y'], O); satline!(pl, M[1,:], M[2,:], subplot=subplot) # bord droit
+ M = T(Rθ*[-x'; y'], O); satline!(pl, M[1,:], M[2,:], subplot=subplot) # bord gauche
+ M = T(Rθ*[ x'.-2*Δx; y'], O); satline!(pl, M[1,:], M[2,:], subplot=subplot) # bord gauche (droite)
+ M = T(Rθ*[ scale.*[-cos(α), cos(α)]';
+ scale.*[ sin(α), sin(α)]'], O);
+ satline!(pl, M[1,:], M[2,:], subplot=subplot) # haut
+ M = T(Rθ*[ scale.*[-cos(α), cos(α)]';
+ -scale.*[ sin(α), sin(α)]'], O);
+ satline!(pl, M[1,:], M[2,:], subplot=subplot) # 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,:], subplot=subplot)
+ pa2 = T(R(β-π/2)*panneau, v0+v2+vy); pa = T(Rθ*pa2, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
+ pa3 = T(R(β-π/2)*panneau, v0+v3+vy); pa = T(Rθ*pa3, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
+
+ pa4 = SA(pa1, β, [xb; yb]); pa = T(Rθ*pa4, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
+ pa5 = SA(pa2, β, [xb; yb]); pa = T(Rθ*pa5, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
+ pa6 = SA(pa3, β, [xb; yb]); pa = T(Rθ*pa6, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
+
+ pa7 = SA(pa1, 0, [0; 0]); pa = T(Rθ*pa7, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
+ pa8 = SA(pa2, 0, [0; 0]); pa = T(Rθ*pa8, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
+ pa9 = SA(pa3, 0, [0; 0]); pa = T(Rθ*pa9, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
+
+ pa10 = SA(pa7, -β, [xb; -yb]); pa = T(Rθ*pa10, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
+ pa11 = SA(pa8, -β, [xb; -yb]); pa = T(Rθ*pa11, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
+ pa12 = SA(pa9, -β, [xb; -yb]); pa = T(Rθ*pa12, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
+
+ # Dessin flamme
+ if thrust
+ lthrust_max = 8000.0
+ origin = T(Rθ*[Δx, 0.0], O)
+ thrust_ = thrust_val*Rθ*[lthrust_max, 0.0]
+ quiver!(pl, [origin[1]], [origin[2]], color=:red,
+ quiver=([thrust_[1]], [thrust_[2]]), linewidth=1, subplot=subplot)
+ end
+
+ end;
+
+ # --------------------------------------------------------------------------------------------
+ #
+ # Stars
+ #
+ # --------------------------------------------------------------------------------------------
+ @userplot StarsPlot
+ @recipe function f(cp::StarsPlot)
+ xo, yo, Δx, Δy, I, A, m, n = cp.args
+ seriestype --> :scatter
+ seriescolor --> :white
+ seriesalpha --> A
+ markerstrokewidth --> 0
+ markersize --> 2*I[3]
+ xo.+I[2].*Δx./n, yo.+I[1].*Δy./m
+ end
+
+ mutable struct Stars
+ s
+ i
+ m
+ n
+ a
+ end
+
+ # Etoiles
+ function stars(ρ)
+ n = 200 # nombre de colonnes
+ m = int(ρ*n) # nombre de lignes
+ d = 0.03
+ s = sprand(m, n, d)
+ i = findnz(s)
+ s = dropzeros(s)
+ # alpha
+ amin = 0.6
+ amax = 1.0
+ Δa = amax-amin
+ a = amin.+rand(length(i[3])).*Δa
+ return Stars(s, i, m, n, a)
+ end;
+
+ # --------------------------------------------------------------------------------------------
+ #
+ # Trajectory
+ #
+ # --------------------------------------------------------------------------------------------
+ @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.0, 3.0, length=n)
+ seriesalpha --> range(0.0, 1.0, length=n)
+ aspect_ratio --> 1
+ label --> false
+ x[inds], y[inds]
+ end
+
+ function trajectoire!(px, times, x1, x2, θ, F, thrust, t_current, cx, cy, pt)
+
+ i_current = argmin(abs.(times.-t_current))
+
+ # Trajectoire
+ if t_current>0 && pt[i_current]
+ trajectoryplot!(px, times, cx.+x1, cy.+x2, t_current)
+ end
+
+ # Satellite
+ satellite!(px, position=[cx+x1[i_current]; cy+x2[i_current]], scale=sc_sat,
+ rotate=θ[i_current], thrust=thrust[i_current], thrust_val=F[i_current])
+
+ end;
+
+ # --------------------------------------------------------------------------------------------
+ #
+ # Background
+ #
+ # --------------------------------------------------------------------------------------------
+ # Décor
+ function background(w, h, xmin, xmax, ymin, ymax,
+ X1_orb_init, X2_orb_init, X1_orb_arr, X2_orb_arr,
+ cx, cy, S, Δx, Δy, cam_x, cam_y, cam_zoom, opi, opf, times, time_current)
+
+ # Fond
+ px = plot(background_color=:gray8, legend = false, aspect_ratio=:equal,
+ size = (w, h), framestyle = :none,
+ left_margin=-25mm, bottom_margin=-10mm, top_margin=-10mm, right_margin=-10mm,
+ xlims=(xmin, xmax), ylims=(ymin, ymax) #,
+ #camera=(0, 90), xlabel = "x", ylabel= "y", zlabel = "z", right_margin=25mm
+ )
+
+ # Etoiles
+ starsplot!(px, xmin, ymin, Δx, Δy, S.i, S.a, S.m, S.n)
+
+ starsplot!(px, xmin-Δx, ymin-Δy, Δx, Δy, S.i, S.a, S.m, S.n)
+ starsplot!(px, xmin, ymin-Δy, Δx, Δy, S.i, S.a, S.m, S.n)
+ starsplot!(px, xmin, ymin+Δy, Δx, Δy, S.i, S.a, S.m, S.n)
+ starsplot!(px, xmin-Δx, ymin, Δx, Δy, S.i, S.a, S.m, S.n)
+ starsplot!(px, xmin+Δx, ymin, Δx, Δy, S.i, S.a, S.m, S.n)
+ starsplot!(px, xmin-Δx, ymin+Δy, Δx, Δy, S.i, S.a, S.m, S.n)
+ starsplot!(px, xmin, ymin+Δy, Δx, Δy, S.i, S.a, S.m, S.n)
+ starsplot!(px, xmin+Δx, ymin+Δy, Δx, Δy, S.i, S.a, S.m, S.n)
+
+ # Orbite initiale
+ nn = length(X2_orb_init)
+ plot!(px, cx.+X1_orb_init, cy.+X2_orb_init, color = :olivedrab1, linewidth=2, alpha=opi)
+
+ # Orbite finale
+ nn = length(X2_orb_init)
+ plot!(px, cx.+X1_orb_arr, cy.+X2_orb_arr, color = :turquoise1, linewidth=2, alpha=opf)
+
+ # Terre
+ θ = range(0.0, 2π, length=100)
+ rT = rayon_Terre #- 1000
+ xT = cx .+ rT .* cos.(θ)
+ yT = cy .+ rT .* sin.(θ)
+ nn = length(xT)
+ plot!(px, xT, yT, color = :dodgerblue1, seriestype=:shape, linewidth=0)
+
+ # Soleil
+ i_current = argmin(abs.(times.-time_current))
+ e = π/6
+ β = range(π/4+e, π/4-e, length=length(times))
+ ρ = sqrt((0.8*xmax-cx)^2+(0.8*ymax-cy)^2)
+ cxS = cx + ρ * cos(β[i_current])
+ cyS = cy + ρ * sin(β[i_current])
+ θ = range(0.0, 2π, length=100)
+ rS = 2000
+ xS = cxS .+ rS .* cos.(θ)
+ yS = cyS .+ rS .* sin.(θ)
+ plot!(px, xS, yS, color = :gold, seriestype=:shape, linewidth=0)
+
+ # Point de départ
+ #plot!(px, [cx+x0[1]], [cy+x0[2]], seriestype=:scatter, color = :white, markerstrokewidth=0)
+
+ # Cadre
+ #plot!(px, [xmin, xmax, xmax, xmin, xmin], [ymin, ymin, ymax, ymax, ymin], color=:white)
+
+ # Angle
+ #plot!(px, camera=(30,90))
+
+ return px
+
+ end;
+
+ function panneau_control!(px, time_current, times, u, F_max, subplot_current)
+
+ # We set the (optional) position relative to top-left of the 1st subplot.
+ # The call is `bbox(x, y, width, height, origin...)`,
+ # where numbers are treated as "percent of parent".
+ Δt_control = 2 #0.1*times[end]
+
+ xx = 0.08
+ yy = 0.1
+ plot!(px, times, F_max.*u[1,:],
+ inset = (1, bbox(xx, yy, 0.2, 0.15, :top, :left)),
+ subplot = subplot_current,
+ bg_inside = nothing,
+ color=:red, legend=:true, yguidefontrotation=-90, yguidefonthalign = :right,
+ xguidefontsize=14, legendfontsize=14,
+ ylims=(-F_max,F_max),
+ xlims=(time_current-Δt_control, time_current+1.0*Δt_control),
+ ylabel="u₁ [N]", label=""
+ )
+
+ plot!(px,
+ [time_current, time_current], [-F_max,F_max],
+ subplot = subplot_current, label="")
+
+ plot!(px, times, F_max.*u[2,:],
+ inset = (1, bbox(xx, yy+0.2, 0.2, 0.15, :top, :left)),
+ #ticks = nothing,
+ subplot = subplot_current+1,
+ bg_inside = nothing,
+ color=:red, legend=:true, yguidefontrotation=-90, yguidefonthalign = :right,
+ xguidefontsize=14, legendfontsize=14,
+ ylims=(-F_max,F_max),
+ xlims=(time_current-Δt_control, time_current+1.0*Δt_control),
+ ylabel="u₂ [N]", xlabel="temps [h]", label=""
+ )
+
+ plot!(px,
+ [time_current, time_current], [-F_max,F_max],
+ subplot = subplot_current+1, label="")
+
+ return subplot_current+2
+ end;
+
+ @enum PB tfmin=1 consomin=2
+
+ function panneau_information!(px, subplot_current, time_current, i_current,
+ x1, x2, v1, v2, θ, F_max, tf_transfert, X1_orb_init, X2_orb_init, X1_orb_arr, X2_orb_arr;
+ pb=tfmin, tf_min=0.0, conso=[])
+
+ # panneaux information
+ xx = 0.06
+ yy = 0.3
+ plot!(px,
+ inset = (1, bbox(xx, yy, 0.2, 0.15, :bottom, :left)),
+ subplot=subplot_current,
+ bg_inside = nothing, legend=:false,
+ framestyle=:none
+ )
+
+ a = 0.0
+ b = 0.0
+ Δtxt = 0.3
+
+ txt = @sprintf("Temps depuis départ [h] = %3.2f", time_current)
+ plot!(px, subplot = subplot_current,
+ annotation=((a, b+3*Δtxt, txt)),
+ annotationcolor=:white, annotationfontsize=14,
+ annotationhalign=:left)
+
+ txt = @sprintf("Vitesse satellite [km/h] = %5.2f", sqrt(v1[i_current]^2+v2[i_current]^2))
+ plot!(px, subplot = subplot_current,
+ annotation=((a, b+2*Δtxt, txt)),
+ annotationcolor=:white, annotationfontsize=14,
+ annotationhalign=:left)
+
+ txt = @sprintf("Distance à la Terre [km] = %5.2f", sqrt(x1[i_current]^2+x2[i_current]^2)-rayon_Terre)
+ plot!(px, subplot = subplot_current,
+ annotation=((a, b+Δtxt, txt)),
+ annotationcolor=:white, annotationfontsize=14,
+ annotationhalign=:left)
+
+ txt = @sprintf("Poussée maximale [N] = %5.2f", F_max)
+ plot!(px, subplot = subplot_current,
+ annotation=((a, b-0*Δtxt, txt)),
+ annotationcolor=:white, annotationfontsize=14,
+ annotationhalign=:left)
+
+ txt = @sprintf("Durée du transfert [h] = %5.2f", tf_transfert)
+ plot!(px, subplot = subplot_current,
+ annotation=((a, b-1*Δtxt, txt)),
+ annotationcolor=:white, annotationfontsize=14,
+ annotationhalign=:left)
+
+ if pb==consomin
+ txt = @sprintf("Durée et consommation comparées")
+ plot!(px, subplot = subplot_current,
+ annotation=((a, b-2*Δtxt, txt)),
+ annotationcolor=:green, annotationfontsize=16,
+ annotationhalign=:left)
+
+ txt = @sprintf("au problème en temps minimal :")
+ plot!(px, subplot = subplot_current,
+ annotation=((a, b-2.75*Δtxt, txt)),
+ annotationcolor=:green, annotationfontsize=16,
+ annotationhalign=:left)
+
+ txt = @sprintf("Durée du transfert [relatif] = %5.2f", tf_transfert/tf_min)
+ plot!(px, subplot = subplot_current,
+ annotation=((a, b-3.75*Δtxt, txt)),
+ annotationcolor=:white, annotationfontsize=14,
+ annotationhalign=:left)
+
+ txt = @sprintf("Consommation [relative] = %5.2f", conso[i_current]./tf_min)
+ plot!(px, subplot = subplot_current,
+ annotation=((a, b-4.75*Δtxt, txt)),
+ annotationcolor=:white, annotationfontsize=14,
+ annotationhalign=:left)
+ end
+
+ subplot_current = subplot_current+1
+ plot!(px,
+ inset = (1, bbox(0.3, 0.03, 0.5, 0.1, :top, :left)),
+ subplot = subplot_current,
+ bg_inside = nothing, legend=:false, aspect_ratio=:equal,
+ xlims=(-4, 4),
+ framestyle=:none
+ )
+
+ rmax = 3*maximum(sqrt.(X1_orb_init.^2+X2_orb_init.^2))
+ plot!(px, subplot = subplot_current,
+ 0.04.+X1_orb_init./rmax, X2_orb_init./rmax.-0.03,
+ color = :olivedrab1, linewidth=2)
+
+ rmax = 7*maximum(sqrt.(X1_orb_arr.^2+X2_orb_arr.^2))
+ plot!(px, subplot = subplot_current,
+ 3.3.+X1_orb_arr./rmax, X2_orb_arr./rmax.-0.03,
+ color = :turquoise1, linewidth=2)
+
+ satellite!(px, subplot=subplot_current,
+ position=[0.67; 0.28], scale=0.08,
+ rotate=π/2)
+
+ # Terre
+ θ = range(0.0, 2π, length=100)
+ rT = 0.1
+ xT = 3.55 .+ rT .* cos.(θ)
+ yT = 0.28 .+ rT .* sin.(θ)
+ plot!(px, subplot=subplot_current, xT, yT, color = :dodgerblue1, seriestype=:shape, linewidth=0)
+
+ s1 = "Transfert orbital du satellite autour de la Terre\n"
+ s2 = "de l'orbite elliptique vers l'orbite circulaire\n"
+ if pb==tfmin
+ s3 = "en temps minimal"
+ elseif pb==consomin
+ s3 = "en consommation minimale"
+ else
+ error("pb pb")
+ end
+ txt = s1 * s2 * s3
+ plot!(px, subplot = subplot_current,
+ annotation=((0, 0, txt)),
+ annotationcolor=:white, annotationfontsize=18,
+ annotationhalign=:center)
+
+ # Realisation
+ subplot_current = subplot_current+1
+ xx = 0.0
+ yy = 0.02
+ plot!(px,
+ inset = (1, bbox(xx, yy, 0.12, 0.05, :bottom, :right)),
+ subplot = subplot_current,
+ bg_inside = nothing, legend=:false,
+ framestyle=:none
+ )
+
+ s1 = "Réalisation : Olivier Cots (2022)"
+ txt = s1
+ plot!(px, subplot = subplot_current,
+ annotation=((0, 0, txt)),
+ annotationcolor=:gray, annotationfontsize=8, annotationhalign=:left)
+
+ return subplot_current+1
+
+ end;
+
+end # module
\ No newline at end of file