diff --git a/flows.jl b/flows.jl
new file mode 100644
index 0000000000000000000000000000000000000000..d785a113fc15371e6161dbf6ba06c7ded1db7203
--- /dev/null
+++ b/flows.jl
@@ -0,0 +1,562 @@
+using ForwardDiff
+using OrdinaryDiffEq
+#using DifferentialEquations
+
+# Alias for gradient and jacobian
+grad(f, x) = ForwardDiff.gradient(f, x)
+jac(f, x) = ForwardDiff.jacobian(f, x);
+
+# types of controls
+@enum CONTROL umin=1 umax=2 usingular=3 uboundary=4 #uregular=5
+
+# --------------------------------------------------------------------------------------------
+# Default options for flows
+# --------------------------------------------------------------------------------------------
+function __abstol()
+ return 1e-10
+end
+
+function __reltol()
+ return 1e-10
+end
+
+function __saveat()
+ return []
+end
+
+# --------------------------------------------------------------------------------------------
+# Hamiltonian
+# --------------------------------------------------------------------------------------------
+struct Hamiltonian f::Function end
+
+function (h::Hamiltonian)(x, p) # https://docs.julialang.org/en/v1/manual/methods/#Function-like-objects
+ return h.f(x,p)
+end
+
+# Flow from a Hamiltonian
+function Flow(h::Hamiltonian)
+
+ function rhs!(dz, z, dummy, t)
+ n = size(z, 1)÷2
+ foo = z -> h(z[1:n], z[n+1:2*n])
+ dh = grad(foo, z)
+ dz[1:n] = dh[n+1:2n]
+ dz[n+1:2n] = -dh[1:n]
+ end
+
+ function f(tspan, x0, p0; abstol=__abstol(), reltol=__reltol(), saveat=__saveat())
+ z0 = [ x0 ; p0 ]
+ ode = ODEProblem(rhs!, z0, tspan)
+ sol = solve(ode, Tsit5(), abstol=abstol, reltol=reltol, saveat=saveat)
+ return sol
+ end
+
+ function f(t0, x0, p0, tf; abstol=__abstol(), reltol=__reltol(), saveat=__saveat())
+ sol = f((t0, tf), x0, p0, abstol=abstol, reltol=reltol, saveat=saveat)
+ n = size(x0, 1)
+ return sol[1:n, end], sol[n+1:2*n, end]
+ end
+
+ return f
+
+end;
+
+# --------------------------------------------------------------------------------------------
+#
+# Single input and affine Mayer system
+#
+# --------------------------------------------------------------------------------------------
+struct SIMayer
+ f₀::Function
+ f₁::Function
+ control_bounds::Tuple{Number, Number}
+ constraint::Union{Function, Nothing}
+end
+
+function __vector_fields(ocp::SIMayer)
+ f₀ = ocp.f₀
+ f₁ = ocp.f₁
+ return f₀, f₁
+end
+
+function __hamiltonian_lifts(ocp::SIMayer)
+ f₀ = ocp.f₀
+ f₁ = ocp.f₁
+ h₀(x, p) = p'*f₀(x)
+ h₁(x, p) = p'*f₁(x)
+ return h₀, h₁
+end
+
+function __min_control(ocp::SIMayer)
+ u(x, p) = ocp.control_bounds[1]
+end
+
+function __max_control(ocp::SIMayer)
+ u(x, p) = ocp.control_bounds[2]
+end
+
+function __singular_control(ocp::SIMayer)
+
+ h₀, h₁ = __hamiltonian_lifts(ocp)
+ H₀₁ = Poisson(h₀, h₁)
+ H₀₀₁ = Poisson(h₀, H₀₁)
+ H₁₀₁ = Poisson(h₁, H₀₁)
+ us(x, p) = -H₀₀₁(x, p)/H₁₀₁(x, p)
+
+ return us
+
+end
+
+function __boundary_control(ocp::SIMayer)
+
+ f₀, f₁ = __vector_fields(ocp)
+ g = ocp.constraint
+ ub(x) = -Lie(f₀, g)(x) / Lie(f₁, g)(x)
+
+ return ub
+
+end
+
+function __boundary_multiplier(ocp::SIMayer)
+
+ f₀ = ocp.f₀
+ f₁ = ocp.f₁
+ h₀, h₁ = __hamiltonian_lifts(ocp)
+ H₀₁ = Poisson(h₀, h₁)
+ g = ocp.constraint
+ μb(x, p) = H₀₁(x, p) / Lie(f₁, g)(x)
+
+ return μb
+
+end
+
+#function __gradient_constrain(ocp::SIMayer)
+# g = ocp.constraint
+# ∇g(x) = grad(g, x)
+# return ∇g
+#end
+
+function __hamiltonian_min(ocp::SIMayer)
+ h₀, h₁ = __hamiltonian_lifts(ocp)
+ u = __min_control(ocp)
+ h(x, p) = h₀(x, p) + u(x, p) * h₁(x, p)
+ return h
+end
+
+function __hamiltonian_max(ocp::SIMayer)
+ h₀, h₁ = __hamiltonian_lifts(ocp)
+ u = __max_control(ocp)
+ h(x, p) = h₀(x, p) + u(x, p) * h₁(x, p)
+ return h
+end
+
+function __hamiltonian_sing(ocp::SIMayer)
+ h₀, h₁ = __hamiltonian_lifts(ocp)
+ us = __singular_control(ocp)
+ h(x, p) = h₀(x, p) + us(x, p) * h₁(x, p)
+ return h
+end
+
+function __hamiltonian_bound(ocp::SIMayer)
+ h₀, h₁ = __hamiltonian_lifts(ocp)
+ ub = __boundary_control(ocp)
+ μb = __boundary_multiplier(ocp)
+ g = ocp.constraint
+ h(x, p) = h₀(x, p) + ub(x) * h₁(x, p) + μb(x,p)*g(x)
+ return h
+end
+
+function Hamiltonian(ocp::SIMayer, control::CONTROL)
+ if control==usingular
+ return Hamiltonian(__hamiltonian_sing(ocp))
+ elseif control==umin
+ return Hamiltonian(__hamiltonian_min(ocp))
+ elseif control==umax
+ return Hamiltonian(__hamiltonian_max(ocp))
+ elseif control==uboundary
+ return Hamiltonian(__hamiltonian_bound(ocp))
+ else
+ nothing
+ end
+end;
+
+function Control(ocp::SIMayer, control::CONTROL)
+ if control==usingular
+ return __singular_control(ocp)
+ elseif control==umin
+ return __min_control(ocp)
+ elseif control==umax
+ return __max_control(ocp)
+ elseif control==uboundary
+ return __boundary_control(ocp)
+ else
+ nothing
+ end
+end;
+
+function Multiplier(ocp::SIMayer)
+ return __boundary_multiplier(ocp)
+end;
+
+# --------------------------------------------------------------------------------------------
+#
+# Hamiltonian Vector Field
+#
+# --------------------------------------------------------------------------------------------
+struct HamiltonianVectorField f::Function end
+
+function (hv::HamiltonianVectorField)(x, p) # https://docs.julialang.org/en/v1/manual/methods/#Function-like-objects
+ return hv.f(x,p)
+end
+
+# Fonction permettant de calculer le flot d'un système hamiltonien
+function Flow(hv::HamiltonianVectorField)
+
+ function rhs!(dz, z, dummy, t)
+ n = size(z, 1)÷2
+ dz[:] = hv(z[1:n], z[n+1:2*n])
+ end
+
+ function f(tspan, x0, p0; abstol=__abstol(), reltol=__reltol(), saveat=__saveat())
+ z0 = [ x0 ; p0 ]
+ ode = ODEProblem(rhs!, z0, tspan)
+ sol = solve(ode, Tsit5(), abstol=abstol, reltol=reltol, saveat=saveat)
+ return sol
+ end
+
+ function f(t0, x0, p0, t; abstol=__abstol(), reltol=__reltol(), saveat=__saveat())
+ sol = f((t0, t), x0, p0; abstol=abstol, reltol=reltol, saveat=saveat)
+ n = size(x0, 1)
+ return sol[1:n, end], sol[n+1:2*n, end]
+ end
+
+ return f
+
+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;
+
+# --------------------------------------------------------------------------------------------
+#
+# Pseudo Hamiltonian
+#
+# --------------------------------------------------------------------------------------------
+struct PseudoHamiltonian f::Function end
+
+function (h::PseudoHamiltonian)(x, p, u) # https://docs.julialang.org/en/v1/manual/methods/#Function-like-objects
+ return h.f(x, p, u)
+end
+
+function (h::PseudoHamiltonian)(x, p, p⁰, u)
+ return h.f(x, p, p⁰, u)
+end
+
+# Flow from a pseudo-Hamiltonian
+function Flow(H::PseudoHamiltonian; p⁰::Number=Inf)
+
+ if p⁰==Inf
+ h = H.f
+ else
+ h(x, p, u) = H.f(x, p, p⁰, u)
+ end
+
+ function ∂h∂z(z, u)
+ n = size(z, 1)÷2
+ foo(z) = h(z[1:n], z[n+1:2*n], u)
+ return grad(foo, z)
+ end
+
+ function ∂h∂u(z, u)
+ n = size(z, 1)÷2
+ foo(u) = h(z[1:n], z[n+1:2*n], u)
+ return grad(foo, u)
+ end
+
+ ∂h∂u(x, p, u) = ∂h∂u([x; p], u)
+
+ ∂²h∂u²(z, u) = jac(u->∂h∂u(z, u), u)
+ ∂²h∂z∂u(z, u) = jac(u->∂h∂z(z, u), u)
+
+ function rhs!(dw, w, λ, t)
+ # w = (z, u) = (x, p, u)
+ n, m = λ
+ z = w[1:2n]
+ u = w[2n+1:2n+m]
+ dh = ∂h∂z(z, u)
+ hv = [dh[n+1:2n]; -dh[1:n]]
+ dw[1:2n] = hv
+ dw[2n+1:2n+m] = -∂²h∂u²(z, u)\(∂²h∂z∂u(z, u)'*hv)
+ end
+
+ function f(tspan, x0, p0, u0; abstol=__abstol(), reltol=__reltol(), saveat=__saveat())
+ w0 = [ x0 ; p0; u0 ]
+ λ = [size(x0, 1), size(u0, 1)]
+ ode = ODEProblem(rhs!, w0, tspan, λ)
+ sol = solve(ode, Tsit5(), abstol=abstol, reltol=reltol, saveat=saveat)
+ return sol
+ end
+
+ function f(t0, x0, p0, u0, tf; abstol=__abstol(), reltol=__reltol(), saveat=__saveat())
+ sol = f((t0, tf), x0, p0, u0, abstol=abstol, reltol=reltol, saveat=saveat)
+ n = size(x0, 1)
+ return sol[1:n, end], sol[n+1:2n, end], sol[2n+1:end, end]
+ end
+
+ return f, ∂h∂u
+
+end;
+
+# --------------------------------------------------------------------------------------------
+#
+# Mayer
+#
+# --------------------------------------------------------------------------------------------
+struct Mayer f::Function end
+
+# Flow from Mayer system
+Flow(Σu::Mayer) = Flow(PseudoHamiltonian((x, p, u) -> p'*Σu.f(x,u)));
+
+# --------------------------------------------------------------------------------------------
+#
+# Lagrange
+#
+# --------------------------------------------------------------------------------------------
+struct Lagrange
+ f::Function
+ f⁰::Function
+end
+
+# Flow from Lagrange system
+Flow(Σu::Lagrange, p⁰::Number=-1.0) = Flow(PseudoHamiltonian((x, p, p⁰, u) -> p⁰*Σu.f⁰(x,u)+p'*Σu.f(x,u)), p⁰=p⁰);
+
+# --------------------------------------------------------------------------------------------
+#
+# SIMayer
+#
+# --------------------------------------------------------------------------------------------
+function __flow(rhs!, tspan, x0, p0; abstol, reltol, saveat)
+ z0 = [ x0 ; p0 ]
+ n = size(x0, 1)
+ ode = ODEProblem(rhs!, z0, tspan, n)
+ sol = solve(ode, Tsit5(), abstol=abstol, reltol=reltol, saveat=saveat)
+ return sol
+end
+
+function __flow(rhs!, t0, x0, p0, tf; abstol, reltol, saveat)
+ sol = __flow(rhs!, (t0, tf), x0, p0, abstol=abstol, reltol=reltol, saveat=saveat)
+ n = size(x0, 1)
+ return sol[1:n, end], sol[n+1:2n, end]
+end
+
+function __dh(ocp::SIMayer)
+
+ f₀ = ocp.f₀
+ f₁ = ocp.f₁
+ h₀(x, p) = p'*f₀(x)
+ h₁(x, p) = p'*f₁(x)
+
+ function dh₀(x, p)
+ n = size(x, 1)
+ foo(z) = h₀(z[1:n], z[n+1:2n])
+ return grad(foo, [x; p])
+ end
+
+ function dh₁(x, p)
+ n = size(x, 1)
+ foo(z) = h₁(z[1:n], z[n+1:2n])
+ return grad(foo, [x; p])
+ end
+
+ return dh₀, dh₁
+
+end
+
+function __FlowMIN(ocp::SIMayer)
+
+ u_bounds = ocp.control_bounds
+ dh₀, dh₁ = __dh(ocp)
+
+ function rhs!(dz, z, n, t)
+ x = z[1:n]
+ p = z[n+1:2n]
+ u = u_bounds[1]
+ dh0 = dh₀(x, p)
+ dh1 = dh₁(x, p)
+ hv0 = [dh0[n+1:2n]; -dh0[1:n]]
+ hv1 = [dh1[n+1:2n]; -dh1[1:n]]
+ dz[1:2n] = hv0 + u*hv1
+ end
+
+ function f(tspan, x0, p0; abstol=__abstol(), reltol=__reltol(), saveat=__saveat())
+ return __flow(rhs!, tspan, x0, p0, abstol=abstol, reltol=reltol, saveat=saveat)
+ end
+
+ function f(t0, x0, p0, tf; abstol=__abstol(), reltol=__reltol(), saveat=__saveat())
+ return __flow(rhs!, t0, x0, p0, tf, abstol=abstol, reltol=reltol, saveat=saveat)
+ end
+
+ return f
+
+end;
+
+function __FlowMAX(ocp::SIMayer)
+
+ u_bounds = ocp.control_bounds
+ dh₀, dh₁ = __dh(ocp)
+
+ function rhs!(dz, z, n, t)
+ x = z[1:n]
+ p = z[n+1:2n]
+ u = u_bounds[2]
+ dh0 = dh₀(x, p)
+ dh1 = dh₁(x, p)
+ hv0 = [dh0[n+1:2n]; -dh0[1:n]]
+ hv1 = [dh1[n+1:2n]; -dh1[1:n]]
+ dz[1:2n] = hv0 + u*hv1
+ end
+
+ function f(tspan, x0, p0; abstol=__abstol(), reltol=__reltol(), saveat=__saveat())
+ return __flow(rhs!, tspan, x0, p0, abstol=abstol, reltol=reltol, saveat=saveat)
+ end
+
+ function f(t0, x0, p0, tf; abstol=__abstol(), reltol=__reltol(), saveat=__saveat())
+ return __flow(rhs!, t0, x0, p0, tf, abstol=abstol, reltol=reltol, saveat=saveat)
+ end
+
+ return f
+
+end;
+
+function __us(ocp::SIMayer)
+
+ f₀ = ocp.f₀
+ f₁ = ocp.f₁
+ h₀(x, p) = p'*f₀(x)
+ h₁(x, p) = p'*f₁(x)
+
+ # singular control
+ H₀₁ = Poisson(h₀, h₁)
+ H₀₀₁ = Poisson(h₀, H₀₁)
+ H₁₀₁ = Poisson(h₁, H₀₁)
+ us(x, p) = -H₀₀₁(x, p)/H₁₀₁(x, p)
+
+ return us
+
+end
+
+function __FlowSING(ocp::SIMayer)
+
+ dh₀, dh₁ = __dh(ocp)
+ us = __us(ocp)
+
+ function rhs!(dz, z, n, t)
+ x = z[1:n]
+ p = z[n+1:2n]
+ u = us(x, p)
+ dh0 = dh₀(x, p)
+ dh1 = dh₁(x, p)
+ hv0 = [dh0[n+1:2n]; -dh0[1:n]]
+ hv1 = [dh1[n+1:2n]; -dh1[1:n]]
+ dz[1:2n] = hv0 + u*hv1
+ end
+
+ function f(tspan, x0, p0; abstol=__abstol(), reltol=__reltol(), saveat=__saveat())
+ return __flow(rhs!, tspan, x0, p0, abstol=abstol, reltol=reltol, saveat=saveat)
+ end
+
+ function f(t0, x0, p0, tf; abstol=__abstol(), reltol=__reltol(), saveat=__saveat())
+ return __flow(rhs!, t0, x0, p0, tf, abstol=abstol, reltol=reltol, saveat=saveat)
+ end
+
+ return f, us
+
+end;
+
+function __u_μ_boundary(ocp::SIMayer)
+ f₀ = ocp.f₀
+ f₁ = ocp.f₁
+ h₀(x, p) = p'*f₀(x)
+ h₁(x, p) = p'*f₁(x)
+ H₀₁ = Poisson(h₀, h₁)
+ g = ocp.constraint
+ ub(x) = -Lie(f₀, g)(x) / Lie(f₁, g)(x)
+ μb(x, p) = H₀₁(x, p) / Lie(f₁, g)(x)
+ ∇g(x) = grad(g, x)
+ return ub, μb, ∇g
+end
+
+function __FlowBOUND(ocp::SIMayer)
+
+ dh₀, dh₁ = __dh(ocp)
+ ub, μb, ∇g = __u_μ_boundary(ocp)
+
+ function rhs!(dz, z, n, t)
+ x = z[1:n]
+ p = z[n+1:2n]
+ u = ub(x)
+ μ = μb(x, p)
+ dh0 = dh₀(x, p)
+ dh1 = dh₁(x, p)
+ hv0 = [dh0[n+1:2n]; -dh0[1:n]]
+ hv1 = [dh1[n+1:2n]; -dh1[1:n]]
+ dz[1:2n] = hv0 + u*hv1
+ dz[n+1:2n] = dz[n+1:2n] - μ*∇g(x)
+ end
+
+ function f(tspan, x0, p0; abstol=__abstol(), reltol=__reltol(), saveat=__saveat())
+ return __flow(rhs!, tspan, x0, p0, abstol=abstol, reltol=reltol, saveat=saveat)
+ end
+
+ function f(t0, x0, p0, tf; abstol=__abstol(), reltol=__reltol(), saveat=__saveat())
+ return __flow(rhs!, t0, x0, p0, tf, abstol=abstol, reltol=reltol, saveat=saveat)
+ end
+
+ return f, ub, μb
+
+end;
+
+function Flow(ocp::SIMayer, control::CONTROL)
+ if control==usingular
+ return __FlowSING(ocp)
+ elseif control==umin
+ return __FlowMIN(ocp)
+ elseif control==umax
+ return __FlowMAX(ocp)
+ elseif control==uboundary
+ return __FlowBOUND(ocp)
+ else
+ nothing
+ end
+end;
diff --git a/minimum-time-orbital-transfer.ipynb b/minimum-time-orbital-transfer.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..dee2087226353838fc0850817928ba04999a8611
--- /dev/null
+++ b/minimum-time-orbital-transfer.ipynb
@@ -0,0 +1,803 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "baa07c4f-8a11-49d2-b052-e64d34594294",
+ "metadata": {},
+ "source": [
+ "# Orbital transfer in minimum time\n",
+ "\n",
+ "* Author: Olivier Cots\n",
+ "* Date: 05/06/2022\n",
+ "* For the \"Maison de Fermat\""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "b32c3c7c-dcc2-4dd9-9143-b4daf05f1c88",
+ "metadata": {},
+ "source": [
+ "We consider the orbital transfer in minimum time of a satellite around the Earth, \n",
+ "from an initial elliptic orbit to a circular target orbit.\n",
+ "This optimal control problem can be written in the form\n",
+ "\n",
+ "$$\n",
+ "\\left\\lbrace\n",
+ "\\begin{array}{l}\n",
+ " \\min J(x, u, t_f) = t_f \\\\[1.0em]\n",
+ " \\ \\ \\dot{x}_{1}(t) = ~ x_{3}(t) \\\\\n",
+ " \\ \\ \\dot{x}_{2}(t) = ~ x_{4}(t) \\\\\n",
+ " \\ \\ \\dot{x}_{3}(t) = -\\dfrac{\\mu\\, x_{1}(t)}{r^{3}(t)} + u_{1}(t) \\\\\n",
+ " \\ \\ \\dot{x}_{4}(t) = -\\dfrac{\\mu\\, x_{2}(t)}{r^{3}(t)}+u_{2}(t), ~~ ||u(t)|| \\leq \\gamma_\\mathrm{max}, ~~ t \\in [0,t_f] ~~ \\text{a. e}, ~~ u(t) = (u_1(t),u_2(t)), \\\\[1.0em]\n",
+ " \\ \\ x_{1}(0)=x_{0,1}, ~~ x_{2}(0)=x_{0,2}, ~~ x_{3}(0)=x_{0,3}, ~~ x_{4}(0)=x_{0,4}, \\\\\n",
+ " \\ \\ r^2(t_f) = r_{f}^2, ~~ x_{3}(t_f)=-\\sqrt{\\dfrac{\\mu}{r_{f}^{3}}}x_{2}(t_f), ~~ x_{4} (t_f)= \\sqrt{\\dfrac{\\mu}{r_{f}^{3}}}x_{1}(t_f), \\\\\n",
+ "\\end{array}\n",
+ "\\right.\n",
+ "$$\n",
+ "\n",
+ "where $r(t)=\\sqrt{x_{1}^{2}(t)+x_{2}^{2}(t)}$.\n",
+ "The chosen units are the kilometer for distances et the hour for the times. We set the following parameters:\n",
+ " \n",
+ "$$\n",
+ "\\mu=5.1658620912 \\times 10^{12} \\ \\mathrm{km}^{3}.\\mathrm{h}^{-2}, \\quad r_{f} = 42165 \\ \\mathrm{km}.\n",
+ "$$\n",
+ "\n",
+ "The parameter $\\gamma_\\mathrm{max}$ depends on the maximal thrust $F_\\mathrm{max}$ according to\n",
+ "\n",
+ "$$\n",
+ "\\gamma_\\mathrm{max} = \\frac{F_\\mathrm{max}\\times 3600^2}{m} \n",
+ "$$\n",
+ "\n",
+ "where $m$ is the mass of the satellite that we fix to $m=2000\\ \\mathrm{kg}$."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "a552de56-944a-412b-97b4-a71c52123786",
+ "metadata": {},
+ "source": [
+ "## Packages"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "104fc457-beaa-4783-8ce2-bb32ff0ab9e3",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "using MINPACK # fsolve\n",
+ "using LinearAlgebra # norm"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "0995cf52-227b-4c1f-88cd-ec3f6bcde259",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "include(\"flows.jl\");\n",
+ "include(\"plottings.jl\");"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "9b5d4e41-f03e-4e62-9970-8c0c24642119",
+ "metadata": {},
+ "source": [
+ "## Constants of the problems"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "f39bc943-a1d2-47a7-81ea-3d766c9d11a6",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "x0 = [-42272.67, 0, 0, -5796.72] # état initial\n",
+ "μ = 5.1658620912*1e12\n",
+ "rf = 42165.0 ;\n",
+ "F_max = 20.0 #100.0\n",
+ "γ_max = F_max*3600.0^2/(2000.0*10^3)\n",
+ "t0 = 0.0\n",
+ "rf3 = rf^3 ;\n",
+ "α = sqrt(μ/rf3);"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "80679ed8-77de-43aa-b542-6a5f7733e86b",
+ "metadata": {},
+ "source": [
+ "## The Hamiltonian system"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "2139105f-a494-4692-8225-46c165e1fae0",
+ "metadata": {},
+ "source": [
+ "The pseudo-Hamiltonian is\n",
+ "\n",
+ "$$\n",
+ "H(x, p, u) = p_{1}x_3 + p_{2}x_4 + p_{3}\\left(-\\dfrac{\\mu\\, x_{1}}{r^{3}} + u_{1}\\right) + p_{4}\\left(-\\dfrac{\\mu\\, x_{2}}{r^{3}} + u_{2}\\right)\n",
+ "$$\n",
+ "\n",
+ "and the pseudo-Hamiltonian system is given by\n",
+ "\n",
+ "$$\n",
+ " \\vec{H}(x, p, u) = \\left(\\frac{\\partial H}{\\partial p}(x, p, u), \n",
+ " -\\frac{\\partial H}{\\partial x}(x, p, u) \\right)\n",
+ "$$\n",
+ "\n",
+ "where \n",
+ "\n",
+ "$$\n",
+ "\\begin{aligned}\n",
+ "\\frac{\\partial H}{\\partial p}(x, p, u) &= \\left( x_3, x_4, -\\dfrac{\\mu\\, x_{1}}{r^{3}} + u_{1}, -\\dfrac{\\mu\\, x_{2}}{r^{3}} + u_{2}\\right), \\\\\n",
+ "\\frac{\\partial H}{\\partial x}(x, p, u) &= \\left(\\left(\\frac{\\mu}{r^3} - \\frac{3\\mu x_1^2}{r^5}\\right)p_3 + \\frac{3\\mu x_1x_2}{r^5}p_4, \\frac{3\\mu x_1x_2}{r^5}p_3 + \\left(\\frac{\\mu}{r^3} + \\frac{3\\mu x_2^2}{r^5}\\right)p_4, p_1, p_2\\right).\n",
+ "\\end{aligned}\n",
+ "$$\n",
+ "\n",
+ "The maximization condition from the Pontryagin Maximum Principle gives the maximizing control\n",
+ "(we assume that $(p_3,p_4)$ never vanishes) :\n",
+ "\n",
+ "$$\n",
+ " u(x(t),p(t)) = \\frac{\\gamma_\\mathrm{max}}{\\sqrt{p_3^2(t)+p_4^2(t)}}(p_3(t),p_4(t)).\n",
+ "$$\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "650273cc-a4fc-46f9-ba02-a65f1a57ed90",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Maximizing control\n",
+ "function control(p)\n",
+ " u = zeros(eltype(p),2)\n",
+ " u[1] = p[3]*γ_max/norm(p[3:4])\n",
+ " u[2] = p[4]*γ_max/norm(p[3:4])\n",
+ " return u\n",
+ "end;\n",
+ "\n",
+ "# Maximized Hamiltonian\n",
+ "function hfun(x, p)\n",
+ " u = control(p)\n",
+ " nor = norm(x[1:2]); nor3 = nor^3\n",
+ " h = p[1]*x[3] + p[2]*x[4] + p[3]*(-μ*x[1]/nor3 + u[1]) + p[4]*(-μ*x[2]/nor^3 + u[2])\n",
+ " return h\n",
+ "end;\n",
+ "\n",
+ "# Hamiltonian system\n",
+ "function hv(x, p)\n",
+ " n = size(x, 1)\n",
+ " hv = zeros(eltype(x), 2*n)\n",
+ " x_1 = x[1]; x_2 = x[2]; x_3 = x[3]; x_4 = x[4]\n",
+ " p_1 = p[1]; p_2 = p[2]; p_3 = p[3]; p_4 = p[4]\n",
+ " norme = norm(x[1:2]); r3 = norme^3; r5 = norme^5\n",
+ " u = control(p)\n",
+ " hv[1] = x_3;\n",
+ " hv[2] = x_4\n",
+ " hv[3] = -μ*x_1/r3 + u[1]\n",
+ " hv[4] = -μ*x_2/r3 + u[2]\n",
+ " hv[5] = μ*((1/r3 - (3*x_1^2)/r5)*p_3 - ((3*x_1*x_2)/r5)*p_4)\n",
+ " hv[6] = μ*(-((3*x_1*x_2)/r5)*p_3 + (1/r3 - (3*x_2^2)/r5)*p_4)\n",
+ " hv[7] = -p_1\n",
+ " hv[8] = -p_2\n",
+ " return hv\n",
+ "end;\n",
+ "\n",
+ "f = Flow(HamiltonianVectorField(hv));"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "6bda152c-5c5f-4de3-aad9-6e7b0f7975f4",
+ "metadata": {},
+ "source": [
+ "## The shooting function"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "7a35da8c-771e-4e57-8413-fc615d140474",
+ "metadata": {},
+ "source": [
+ "We introduce\n",
+ "$\n",
+ " \\alpha := \\sqrt{\\frac{\\mu}{r_f^3}}.\n",
+ "$\n",
+ "The final condition can be written $c(x(t_f)) = 0$, with $c \\colon \\mathbb{R}^4 \\to \\mathbb{R}^3$\n",
+ "given by\n",
+ "\n",
+ "$$\n",
+ " c(x) := (r^2(x) - r_f^2, ~~ x_{3} + \\alpha\\, x_{2}, ~~ x_{4} - \\alpha\\, x_{1}).\n",
+ "$$\n",
+ "\n",
+ "The final time being free, we get the condition at the final time\n",
+ "\n",
+ "$$\n",
+ " H(x(t_f), p(t_f), u(t_f)) = -p^0 = 1. \\quad \\text{(cas normal)}\n",
+ "$$\n",
+ "\n",
+ "Besides, the transversality condition\n",
+ "\n",
+ "$$\n",
+ "p(t_f) = c'(x(t_f))^T \\lambda, ~~ \\lambda \\in \\mathbb{R}^3,\n",
+ "$$\n",
+ "\n",
+ "leads to\n",
+ "\n",
+ "$$\n",
+ "\\Phi(x(t_f), p(t_f)) := x_2(t_f) \\Big( p_1(t_f) + \\alpha\\, p_4(t_f) \\Big) - x_1(t_f) \\Big( p_2(t_f) - \\alpha\\, p_3(t_f) \\Big) = 0.\n",
+ "$$\n",
+ "\n",
+ "Considering the limit conditions, the final condition on the pseudo-Hamiltonian and the transversality condition,\n",
+ "the shooting function is given by\n",
+ "\n",
+ "\\begin{equation*}\n",
+ " \\begin{array}{rlll}\n",
+ " S \\colon & \\mathbb{R}^5 & \\longrightarrow & \\mathbb{R}^5 \\\\\n",
+ " & (p_0, t_f) & \\longmapsto &\n",
+ " S(p_0, t_f) := \\begin{pmatrix}\n",
+ " c(x(t_f, x_0, p_0)) \\\\[0.5em]\n",
+ " \\Phi(z(t_f, x_0, p_0)) \\\\[0.5em]\n",
+ " H(z(t_f, x_0, p_0), u(z(t_f, x_0, p_0))) - 1\n",
+ " \\end{pmatrix}\n",
+ " \\end{array}\n",
+ "\\end{equation*}\n",
+ "\n",
+ "where $z(t_f, x_0, p_0)$ is the solution at time $t_f$ of the feedback pseudo-Hamiltonian system (that is the \n",
+ "pseudo-Hamiltonian system with the maximizing control), starting at time $t_0$ from $(x_0, p_0)$. We recall that\n",
+ "we note $z=(x, p)$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "4deeb191-35a8-4f5b-9dca-efa1e7b1ac8b",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Iter f(x) inf-norm Step 2-norm Step time\n",
+ "------ -------------- -------------- --------------\n",
+ " 1 7.026435e+00 0.000000e+00 5.059817\n",
+ " 2 7.334182e-03 4.372294e-11 16.066451\n",
+ " 3 3.692279e-03 1.109605e-11 0.018698\n",
+ " 4 9.853102e-08 1.134859e-11 0.022441\n",
+ " 5 1.590524e-08 8.354949e-22 0.019856\n",
+ "Results of Nonlinear Solver Algorithm\n",
+ " * Algorithm: Modified Powell (User Jac, Expert)\n",
+ " * Starting Point: [-0.0013615, -7.34989e-6, -5.359923e-5, -0.00858271, 59.8551668]\n",
+ " * Zero: [-0.0013615711107251762, -7.349896988606744e-6, -5.35992371862138e-5, -0.008582713089817779, 59.85516688789379]\n",
+ " * Inf-norm of residuals: 0.000000\n",
+ " * Convergence: true\n",
+ " * Message: algorithm estimates that the relative error between x and the solution is at most tol\n",
+ " * Total time: 21.195272 seconds\n",
+ " * Function Calls: 5\n",
+ " * Jacobian Calls (df/dx): 1\n"
+ ]
+ },
+ {
+ "data": {
+ "text/plain": [
+ "Transfert(59.85516688789379, [-0.0013615711107251762, -7.349896988606744e-6, -5.35992371862138e-5, -0.008582713089817779])"
+ ]
+ },
+ "execution_count": 5,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "# Données pour la trajectoire durant le transfert\n",
+ "mutable struct Transfert\n",
+ " duration\n",
+ " initial_adjoint\n",
+ "end\n",
+ "\n",
+ "# Fonction de tir\n",
+ "function shoot(p0, tf)\n",
+ " \n",
+ " # Integration\n",
+ " xf, pf = f(t0, x0, p0, tf)\n",
+ " \n",
+ " # Conditions\n",
+ " s = zeros(eltype(p0), 5)\n",
+ " s[1] = norm(xf[1:2]) - rf\n",
+ " s[2] = xf[3] + α*xf[2]\n",
+ " s[3] = xf[4] - α*xf[1]\n",
+ " s[4] = xf[2]*(pf[1] + α*pf[4]) - xf[1]*(pf[2] - α*pf[3])\n",
+ " s[5] = hfun(xf, pf) - 1\n",
+ " \n",
+ " return s\n",
+ "\n",
+ "end;\n",
+ "\n",
+ "# Initial guess\n",
+ "#y_guess = [1.0323e-4, 4.915e-5, 3.568e-4, -1.554e-4, 13.4] # for F_max = 100N\n",
+ "ξ_guess = [-0.0013615, -7.34989e-6, -5.359923e-5, -0.00858271, 59.8551668] # for F_max = 20N\n",
+ "\n",
+ "# Solve\n",
+ "foo(ξ) = shoot(ξ[1:4], ξ[5])\n",
+ "jfoo(ξ) = ForwardDiff.jacobian(foo, ξ)\n",
+ "foo!(s, ξ) = ( s[:] = foo(ξ); nothing )\n",
+ "jfoo!(js, ξ) = ( js[:] = jfoo(ξ); nothing )\n",
+ "\n",
+ "nl_sol = fsolve(foo!, jfoo!, ξ_guess, show_trace=true); println(nl_sol)\n",
+ "\n",
+ "# Retrieves solution\n",
+ "if nl_sol.converged\n",
+ " p0_sol = nl_sol.x[1:4]\n",
+ " tf_sol = nl_sol.x[5]\n",
+ " transfert_data = Transfert(tf_sol, p0_sol);\n",
+ "else\n",
+ " error(\"Not converged\")\n",
+ "end"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ae9cd7d1-97b2-465f-8f05-f2eeea9b1def",
+ "metadata": {},
+ "source": [
+ "## Pre and post-transfer"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "id": "73dcf5ff-4ff2-4894-9407-720bed9b4983",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Kepler's law\n",
+ "function kepler(x)\n",
+ " n = size(x, 1)\n",
+ " dx = zeros(eltype(x), n)\n",
+ " x1 = x[1]\n",
+ " x2 = x[2]\n",
+ " x3 = x[3]\n",
+ " x4 = x[4]\n",
+ " dsat = norm(x[1:2]); r3 = dsat^3;\n",
+ " dx[1] = x3\n",
+ " dx[2] = x4\n",
+ " dx[3] = -μ*x1/r3\n",
+ " dx[4] = -μ*x2/r3\n",
+ " return dx\n",
+ "end\n",
+ "\n",
+ "f0 = Flow(VectorField(kepler));"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "id": "04f4e32f-ce8a-4e1a-90a1-b3dc0af2a956",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Iter f(x) inf-norm Step 2-norm Step time\n",
+ "------ -------------- -------------- --------------\n",
+ " 1 3.021414e-07 0.000000e+00 3.061420\n",
+ " 2 6.839400e-10 1.399191e-13 4.821815\n",
+ "Results of Nonlinear Solver Algorithm\n",
+ " * Algorithm: Modified Powell (User Jac, Expert)\n",
+ " * Starting Point: [6738.200652584018, 1.3618799608695503e-21, 2.8839710428042565e-7, 36366.2117341436, 5.302395297743802]\n",
+ " * Zero: [6738.200652696158, 2.2925434185945747e-21, 1.4403503883997278e-7, 36366.211733817254, 5.302395297772681]\n",
+ " * Inf-norm of residuals: 0.000000\n",
+ " * Convergence: true\n",
+ " * Message: algorithm estimates that the relative error between x and the solution is at most tol\n",
+ " * Total time: 7.883518 seconds\n",
+ " * Function Calls: 2\n",
+ " * Jacobian Calls (df/dx): 1\n"
+ ]
+ },
+ {
+ "data": {
+ "text/plain": [
+ "PreTransfert(10.604790595545362, [-42272.67, 0.0, 0.0, -5796.72])"
+ ]
+ },
+ "execution_count": 7,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "# On cherche le point le plus proche de la Terre sur l'orbite initiale\n",
+ "# Ce sera le point de départ de la pré-mission\n",
+ "# On cherche aussi le temps\n",
+ "\n",
+ "# données pour la trajectoire pré-transfert\n",
+ "mutable struct PreTransfert\n",
+ " duration\n",
+ " initial_point\n",
+ "end\n",
+ "\n",
+ "# Fonction de tir\n",
+ "function depart(xp, tp)\n",
+ " \n",
+ " # Integration\n",
+ " x0_ = f0(0.0, xp, tp)\n",
+ " \n",
+ " # Conditions\n",
+ " s = zeros(eltype(xp), 5)\n",
+ " s[1] = xp[2]\n",
+ " s[2] = x0_[1] - x0[1]\n",
+ " s[3] = x0_[2] - x0[2]\n",
+ " s[4] = x0_[3] - x0[3]\n",
+ " s[5] = x0_[4] - x0[4]\n",
+ " \n",
+ " return s\n",
+ "\n",
+ "end;\n",
+ "\n",
+ "# Itéré initial\n",
+ "ξ_guess = [6738.200652584018, 1.3618799608695503e-21, 2.8839710428042565e-7, 36366.2117341436, 5.302395297743802]\n",
+ "\n",
+ "# Solve\n",
+ "foo(ξ) = depart(ξ[1:4], ξ[5])\n",
+ "jfoo(ξ) = ForwardDiff.jacobian(foo, ξ)\n",
+ "foo!(s, ξ) = ( s[:] = foo(ξ); nothing )\n",
+ "jfoo!(js, ξ) = ( js[:] = jfoo(ξ); nothing )\n",
+ "\n",
+ "nl_sol = fsolve(foo!, jfoo!, ξ_guess, show_trace=true); println(nl_sol)\n",
+ "\n",
+ "# Retrieves solution\n",
+ "if nl_sol.converged\n",
+ " x_pre_transfert = x0\n",
+ " t_pre_transfert = 2.0*nl_sol.x[5]\n",
+ " pre_transfert_data = PreTransfert(t_pre_transfert, x_pre_transfert);\n",
+ "else\n",
+ " error(\"Not converged\")\n",
+ "end"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "id": "00a00a7d-a161-4327-9ea3-f8e4a34bfb77",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# données pour la trajectoire post-transfert\n",
+ "mutable struct PostTransfert\n",
+ " duration\n",
+ "end\n",
+ "\n",
+ "t_post_transfert = 20.0\n",
+ "\n",
+ "post_transfert_data = PostTransfert(t_post_transfert);"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "0686c07f-aea9-405a-95bd-af42f09370b6",
+ "metadata": {},
+ "source": [
+ "## Animation"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "id": "e61c1bee-88b4-4b57-8a09-1249f31978e0",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "function animation(pre_transfert_data, transfert_data, post_transfert_data; fps=10, nFrame=200)\n",
+ " \n",
+ " # pré-transfert\n",
+ " t_pre_transfert = pre_transfert_data.duration\n",
+ " x_pre_transfert = pre_transfert_data.initial_point\n",
+ " \n",
+ " # transfert\n",
+ " p0_transfert = transfert_data.initial_adjoint\n",
+ " tf_transfert = transfert_data.duration\n",
+ " \n",
+ " # post-trasfert\n",
+ " t_post_transfert = post_transfert_data.duration\n",
+ " \n",
+ " # On calcule les orbites initiale et finale\n",
+ " r0 = norm(x0[1:2])\n",
+ " v0 = norm(x0[3:4])\n",
+ " a = 1.0/(2.0/r0-v0*v0/μ)\n",
+ " t1 = r0*v0*v0/μ - 1.0;\n",
+ " t2 = (x0[1:2]'*x0[3:4])/sqrt(a*μ);\n",
+ " e_ellipse = norm([t1 t2])\n",
+ " p_orb = a*(1-e_ellipse^2);\n",
+ " n_theta = 151\n",
+ " Theta = range(0.0, stop=2*pi, length=n_theta)\n",
+ " X1_orb_init = zeros(n_theta)\n",
+ " X2_orb_init = zeros(n_theta)\n",
+ " X1_orb_arr = zeros(n_theta)\n",
+ " X2_orb_arr = zeros(n_theta)\n",
+ "\n",
+ " for i in 1:n_theta\n",
+ " theta = Theta[i]\n",
+ " r_orb = p_orb/(1+e_ellipse*cos(theta));\n",
+ " # Orbite initiale\n",
+ " X1_orb_init[i] = r_orb*cos(theta);\n",
+ " X2_orb_init[i] = r_orb*sin(theta);\n",
+ " # Orbite finale\n",
+ " X1_orb_arr[i] = rf*cos(theta) ;\n",
+ " X2_orb_arr[i] = rf*sin(theta);\n",
+ " end\n",
+ " \n",
+ " # Centre de la fenêtre\n",
+ " cx = 40000\n",
+ " cy = -7000\n",
+ "\n",
+ " # Taille de la fenêtre\n",
+ " w = 1600\n",
+ " h = 900\n",
+ " ρ = h/w\n",
+ " \n",
+ " # Limites de la fenêtre\n",
+ " ee = 0.5\n",
+ " xmin = minimum([X1_orb_init; X1_orb_arr]); xmin = xmin - ee * abs(xmin);\n",
+ " xmax = maximum([X1_orb_init; X1_orb_arr]); xmax = xmax + ee * abs(xmax);\n",
+ " ymin = minimum([X2_orb_init; X2_orb_arr]); ymin = ymin - ee * abs(ymin);\n",
+ " ymax = maximum([X2_orb_init; X2_orb_arr]); ymax = ymax + ee * abs(ymax);\n",
+ " \n",
+ " Δy = ymax-ymin\n",
+ " Δx = Δy/ρ\n",
+ " Δn = Δx - (xmax-xmin)\n",
+ " xmin = xmin - Δn/2\n",
+ " xmax = xmax + Δn/2\n",
+ " \n",
+ " # Trajectoire pré-transfert\n",
+ " traj_pre_transfert = f0((0.0, t_pre_transfert), x_pre_transfert)\n",
+ " times_pre_transfert = traj_pre_transfert.t\n",
+ " n_pre_transfert = size(times_pre_transfert, 1)\n",
+ " x1_pre_transfert = [traj_pre_transfert[1, j] for j in 1:n_pre_transfert ]\n",
+ " x2_pre_transfert = [traj_pre_transfert[2, j] for j in 1:n_pre_transfert ]\n",
+ " v1_pre_transfert = [traj_pre_transfert[3, j] for j in 1:n_pre_transfert ]\n",
+ " v2_pre_transfert = [traj_pre_transfert[4, j] for j in 1:n_pre_transfert ] \n",
+ " \n",
+ " # Trajectoire pendant le transfert\n",
+ " traj_transfert = f((t0, tf_transfert), x0, p0_transfert)\n",
+ " times_transfert = traj_transfert.t\n",
+ " n_transfert = size(times_transfert, 1)\n",
+ " x1_transfert = [traj_transfert[1, j] for j in 1:n_transfert ]\n",
+ " x2_transfert = [traj_transfert[2, j] for j in 1:n_transfert ]\n",
+ " v1_transfert = [traj_transfert[3, j] for j in 1:n_transfert ]\n",
+ " v2_transfert = [traj_transfert[4, j] for j in 1:n_transfert ]\n",
+ " u_transfert = zeros(2, length(times_transfert))\n",
+ " for j in 1:n_transfert\n",
+ " u_transfert[:,j] = control(traj_transfert[5:8, j])\n",
+ " end\n",
+ "\n",
+ " # post-transfert\n",
+ " x_post_transfert = traj_transfert[:,end]\n",
+ "\n",
+ " # Trajectoire post-transfert\n",
+ " traj_post_transfert = f0((0.0, t_post_transfert), x_post_transfert)\n",
+ " times_post_transfert = traj_post_transfert.t\n",
+ " n_post_transfert = size(times_post_transfert, 1)\n",
+ " x1_post_transfert = [traj_post_transfert[1, j] for j in 1:n_post_transfert ]\n",
+ " x2_post_transfert = [traj_post_transfert[2, j] for j in 1:n_post_transfert ]\n",
+ " v1_post_transfert = [traj_post_transfert[3, j] for j in 1:n_post_transfert ]\n",
+ " v2_post_transfert = [traj_post_transfert[4, j] for j in 1:n_post_transfert ] \n",
+ "\n",
+ " # Angles de rotation du satellite pendant le pré-transfert\n",
+ " θ0 = atan(u_transfert[2, 1], u_transfert[1, 1])\n",
+ " θ_pre_transfert = range(π/2, mod(θ0, 2*π), length=n_pre_transfert)\n",
+ "\n",
+ " # Angles de rotation du satellite pendant le transfert\n",
+ " θ_transfert = atan.(u_transfert[2, :], u_transfert[1, :])\n",
+ "\n",
+ " # Angles de rotation du satellite pendant le post-transfert\n",
+ " θ1 = atan(u_transfert[2, end], u_transfert[1, end])\n",
+ " θ2 = atan(-x2_post_transfert[end], -x1_post_transfert[end])\n",
+ " θ_post_transfert = range(mod(θ1, 2*π), mod(θ2, 2*π), length=n_post_transfert)\n",
+ "\n",
+ " # Etoiles\n",
+ " Δx = xmax-xmin\n",
+ " Δy = ymax-ymin \n",
+ " ρ = Δy/Δx\n",
+ " S = stars(ρ)\n",
+ "\n",
+ " # nombre total d'images\n",
+ " nFrame = min(nFrame, n_pre_transfert+n_transfert+n_post_transfert);\n",
+ " \n",
+ " # Pour l'affichage de la trajectoire globale\n",
+ " times = [times_pre_transfert[1:end-1];\n",
+ " times_pre_transfert[end].+times_transfert[1:end-1];\n",
+ " times_pre_transfert[end].+times_transfert[end].+times_post_transfert[1:end]]\n",
+ " x1 = [x1_pre_transfert[1:end-1]; x1_transfert[1:end-1]; x1_post_transfert[:]]\n",
+ " x2 = [x2_pre_transfert[1:end-1]; x2_transfert[1:end-1]; x2_post_transfert[:]]\n",
+ " v1 = [v1_pre_transfert[1:end-1]; v1_transfert[1:end-1]; v1_post_transfert[:]]\n",
+ " v2 = [v2_pre_transfert[1:end-1]; v2_transfert[1:end-1]; v2_post_transfert[:]]\n",
+ " θ = [ θ_pre_transfert[1:end-1]; θ_transfert[1:end-1]; θ_post_transfert[:]]\n",
+ "\n",
+ " # plot thrust on/off\n",
+ " th = [BitArray(zeros(n_pre_transfert-1)); \n",
+ " BitArray(ones(n_transfert-1));\n",
+ " BitArray(zeros(n_post_transfert))]\n",
+ " \n",
+ " # plot trajectory\n",
+ " pt = [BitArray(ones(n_pre_transfert-1)); \n",
+ " BitArray(ones(n_transfert-1));\n",
+ " BitArray(zeros(n_post_transfert))]\n",
+ " \n",
+ " # Contrôle sur la trajectoire totale\n",
+ " u_total = hcat([zeros(2, n_pre_transfert-1), \n",
+ " u_transfert[:, 1:n_transfert-1],\n",
+ " zeros(2, n_post_transfert)]...)\n",
+ "\n",
+ " # temps total\n",
+ " temps_transfert_global = times[end]\n",
+ "\n",
+ " # pas de temps pour le transfert global\n",
+ " if nFrame>1\n",
+ " Δt = temps_transfert_global/(nFrame-1)\n",
+ " else\n",
+ " Δt = 0.0\n",
+ " end\n",
+ " \n",
+ " # opacités des orbites initiale et finale\n",
+ " op_initi = [range(0.0, 1.0, length=n_pre_transfert-1);\n",
+ " range(1.0, 0.0, length=n_transfert-1);\n",
+ " zeros(n_post_transfert)]\n",
+ " op_final = [zeros(n_pre_transfert-1);\n",
+ " range(0.0, 1.0, length=n_transfert-1);\n",
+ " range(1.0, 0.0, length=int(n_post_transfert/4));\n",
+ " zeros(n_post_transfert-int(n_post_transfert/4))]\n",
+ " \n",
+ " # animation\n",
+ " anim = @animate for i ∈ 1:nFrame\n",
+ " \n",
+ " # Δt : pas de temps\n",
+ " # time_current : temps courant de la mission totale à l'itération i\n",
+ " # i_current : indice tel que times[i_current] = time_current\n",
+ " # w, h : width, height de la fenêtre\n",
+ " # xmin, xmax, ymin, ymax : limites des axes du plot principal\n",
+ " # X1_orb_init, X2_orb_init : coordonnées de l'orbite initiale\n",
+ " # X1_orb_arr, X2_orb_arr : coordonnées de l'orbite finale\n",
+ " # cx, cy : coordonées du centre de l'affichage du tranfert\n",
+ " # S : data pour les étoiles\n",
+ " # Δx, Δy : xmax-xmin, ymax-ymin\n",
+ " # times : tous les temps de la mission complète, ie pre-transfert, transfert et post-transfert\n",
+ " # x1, x2 : vecteur de positions du satellite\n",
+ " # θ : vecteur d'orientations du satellite\n",
+ " # th : vecteur de booléens - thrust on/off\n",
+ " # u_total : vecteur de contrôles pour toute la mission\n",
+ " # F_max, γ_max : poussée max\n",
+ " # subplot_current : valeur du subplot courant\n",
+ " # cam_x, cam_y : position de la caméra\n",
+ " # cam_zoom : zoom de la caméra\n",
+ " \n",
+ " cam_x = cx\n",
+ " cam_y = cy\n",
+ " cam_zoom = 1\n",
+ " \n",
+ " time_current = (i-1)*Δt\n",
+ " i_current = argmin(abs.(times.-time_current))\n",
+ " \n",
+ " px = background(w, h, xmin, xmax, ymin, ymax, \n",
+ " X1_orb_init, X2_orb_init, X1_orb_arr, X2_orb_arr,\n",
+ " cx, cy, S, Δx, Δy, cam_x, cam_y, cam_zoom, \n",
+ " op_initi[i_current], op_final[i_current], times, time_current)\n",
+ "\n",
+ " trajectoire!(px, times, x1, x2, θ, th, time_current, cx, cy, pt) \n",
+ " \n",
+ " subplot_current = 2\n",
+ " subplot_current = panneau_control!(px, time_current, times, u_total, \n",
+ " F_max, γ_max, subplot_current)\n",
+ " \n",
+ " #subplot_current = panneau_information!(px, subplot_current, time_current, \n",
+ " # i_current, x1, x2, v1, v2, θ, F_max, tf_transfert, X1_orb_init, \n",
+ " # X2_orb_init, X1_orb_arr, X2_orb_arr)\n",
+ " \n",
+ " end\n",
+ "\n",
+ " display(plot(px))\n",
+ " # enregistrement\n",
+ " #gif(anim, \"transfert-temps-min-original.mp4\", fps=fps);\n",
+ " #gif(anim, \"transfert-temps-min.gif\", fps=fps);\n",
+ " \n",
+ "end;"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "id": "a23d5c91-72e5-4a93-b689-8227887f0551",
+ "metadata": {},
+ "outputs": [
+ {
+ "ename": "LoadError",
+ "evalue": "Cannot convert Measures.AbsoluteLength to series data for plotting",
+ "output_type": "error",
+ "traceback": [
+ "Cannot convert Measures.AbsoluteLength to series data for plotting",
+ "",
+ "Stacktrace:",
+ " [1] error(s::String)",
+ " @ Base ./error.jl:33",
+ " [2] _prepare_series_data(x::Measures.AbsoluteLength)",
+ " @ RecipesPipeline ~/opt/miniconda3/envs/julia/share/julia/packages/RecipesPipeline/kJPXU/src/series.jl:8",
+ " [3] _series_data_vector(x::Measures.AbsoluteLength, plotattributes::Dict{Symbol, Any})",
+ " @ RecipesPipeline ~/opt/miniconda3/envs/julia/share/julia/packages/RecipesPipeline/kJPXU/src/series.jl:35",
+ " [4] macro expansion",
+ " @ ~/opt/miniconda3/envs/julia/share/julia/packages/RecipesPipeline/kJPXU/src/series.jl:135 [inlined]",
+ " [5] apply_recipe(plotattributes::AbstractDict{Symbol, Any}, #unused#::Type{RecipesPipeline.SliceIt}, x::Any, y::Any, z::Any)",
+ " @ RecipesPipeline ~/opt/miniconda3/envs/julia/share/julia/packages/RecipesBase/qpxEX/src/RecipesBase.jl:289",
+ " [6] _process_userrecipes!(plt::Any, plotattributes::Any, args::Any)",
+ " @ RecipesPipeline ~/opt/miniconda3/envs/julia/share/julia/packages/RecipesPipeline/kJPXU/src/user_recipe.jl:36",
+ " [7] recipe_pipeline!(plt::Any, plotattributes::Any, args::Any)",
+ " @ RecipesPipeline ~/opt/miniconda3/envs/julia/share/julia/packages/RecipesPipeline/kJPXU/src/RecipesPipeline.jl:70",
+ " [8] _plot!(plt::Plots.Plot, plotattributes::Any, args::Any)",
+ " @ Plots ~/opt/miniconda3/envs/julia/share/julia/packages/Plots/GGa6i/src/plot.jl:208",
+ " [9] plot(args::Any; kw::Base.Pairs{Symbol, V, Tuple{Vararg{Symbol, N}}, NamedTuple{names, T}} where {V, N, names, T<:Tuple{Vararg{Any, N}}})",
+ " @ Plots ~/opt/miniconda3/envs/julia/share/julia/packages/Plots/GGa6i/src/plot.jl:91",
+ " [10] plot",
+ " @ ~/opt/miniconda3/envs/julia/share/julia/packages/Plots/GGa6i/src/plot.jl:85 [inlined]",
+ " [11] animation(pre_transfert_data::PreTransfert, transfert_data::Transfert, post_transfert_data::PostTransfert; fps::Int64, nFrame::Int64)",
+ " @ Main ./In[9]:208",
+ " [12] top-level scope",
+ " @ In[10]:4",
+ " [13] eval",
+ " @ ./boot.jl:373 [inlined]",
+ " [14] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)",
+ " @ Base ./loading.jl:1196"
+ ]
+ }
+ ],
+ "source": [
+ "nFrame = 1; fps = 1\n",
+ "#nFrame = 10; fps = 1\n",
+ "#nFrame = 2000; fps = 20\n",
+ "animation(pre_transfert_data, transfert_data, post_transfert_data, nFrame=nFrame, fps=fps)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "id": "acc4155b-8c6a-4e06-a823-8c22fd9ff72a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "if nFrame == 2000\n",
+ " convert = `ffmpeg -y -i transfert-temps-min-original.mp4 -vf format=yuv420p transfert-temps-min.mp4`\n",
+ " run(convert)\n",
+ "end"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "c663a64a-1e53-43ba-8fbf-86fb073cc526",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Julia 1.7.3",
+ "language": "julia",
+ "name": "julia-1.7"
+ },
+ "language_info": {
+ "file_extension": ".jl",
+ "mimetype": "application/julia",
+ "name": "julia",
+ "version": "1.7.3"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/plottings.jl b/plottings.jl
new file mode 100644
index 0000000000000000000000000000000000000000..6c475c7d9c828d553821cd622858278288f79478
--- /dev/null
+++ b/plottings.jl
@@ -0,0 +1,449 @@
+using Plots
+using Plots.PlotMeasures
+using SparseArrays
+using Printf
+
+# 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)
+
+ # 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
+ origin = T(Rθ*[Δx, 0.0], O)
+ thrust = Rθ*[8000, 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, θ, 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])
+
+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, γ_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,:]./γ_max,
+ 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], [-γ_max,γ_max],
+ subplot = subplot_current, label="")
+
+ plot!(px, times, F_max.*u[2,:]./γ_max,
+ 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], [-γ_max,γ_max],
+ subplot = subplot_current+1, label="")
+
+ return subplot_current+2
+end;
+
+
+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)
+
+ # 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)
+ annotate!(px, subplot = subplot_current, a, b+2*Δ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))
+ annotate!(px, subplot = subplot_current, a, b+Δ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)
+ annotate!(px, subplot = subplot_current, a, b, txt,
+ annotationcolor=:white, annotationfontsize=14,
+ annotationhalign=:left)
+
+ txt = @sprintf("Poussée maximale [N] = %5.2f", F_max)
+ annotate!(px, subplot = subplot_current, a, b-2*Δtxt, txt,
+ annotationcolor=:white, annotationfontsize=14,
+ annotationhalign=:left)
+
+ txt = @sprintf("Durée du transfert [h] = %5.2f", tf_transfert)
+ annotate!(px, subplot = subplot_current, a, b-3*Δtxt, txt,
+ annotationcolor=:white, annotationfontsize=14,
+ annotationhalign=:left)
+
+ 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
+ )
+
+ s1 = "Transfert orbital du satellite autour de la Terre\n"
+ s2 = "de l'orbite elliptique vers l'orbite circulaire\n"
+ s3 = "en temps minimal"
+ txt = s1 * s2 * s3
+ annotate!(px, subplot = subplot_current, 0, 0, txt,
+ annotationcolor=:white, annotationfontsize=18,
+ annotationhalign=:center)
+
+ 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)
+
+ 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
+ annotate!(px, subplot = subplot_current, 0, 0, txt,
+ annotationcolor=:gray, annotationfontsize=8, annotationhalign=:left)
+
+ return subplot_current+1
+
+end;
\ No newline at end of file
diff --git a/transfert-temps-min.ipynb b/transfert-temps-min.ipynb
index 9f622f93ea9902b2bb45f060c63e1f907bc2abe1..cbf9553ed4f8c5f0a0e657642950b6b08b233a6b 100644
--- a/transfert-temps-min.ipynb
+++ b/transfert-temps-min.ipynb
@@ -61,7 +61,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
@@ -81,7 +81,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
@@ -97,7 +97,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
@@ -114,7 +114,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
@@ -193,7 +193,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
@@ -354,7 +354,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
@@ -403,7 +403,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
@@ -430,7 +430,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
@@ -504,7 +504,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
@@ -661,7 +661,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
@@ -690,7 +690,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 29,
"metadata": {},
"outputs": [],
"source": [
@@ -927,7 +927,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 30,
"metadata": {},
"outputs": [],
"source": [
@@ -952,9 +952,26 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 31,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Iter f(x) inf-norm Step 2-norm \n",
+ "------ -------------- --------------\n",
+ " 0 8.710000e+01 NaN\n",
+ " 1 4.664840e+01 2.860952e+02\n",
+ " 2 3.948627e-02 1.733579e+00\n",
+ " 3 1.210757e-05 1.418640e-01\n",
+ " 4 3.710738e-10 3.830941e-07\n",
+ "\n",
+ "Final solution:\n",
+ "[6738.200652737065, 1.5644718641515079e-24, 1.1947020276724434e-9, 36366.21173381801, 5.302395297764675]\n"
+ ]
+ }
+ ],
"source": [
"# On cherche le point le plus proche de la Terre sur l'orbite initiale\n",
"# Ce sera le point de départ de la pré-mission\n",
@@ -1048,7 +1065,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
@@ -1144,9 +1161,25 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 33,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Iter f(x) inf-norm Step 2-norm \n",
+ "------ -------------- --------------\n",
+ " 0 7.026523e+00 NaN\n",
+ " 1 7.334383e-03 6.609950e-06\n",
+ " 2 1.068202e-08 6.699837e-06\n",
+ " 3 6.846676e-09 1.949020e-11\n",
+ "\n",
+ "Final solution:\n",
+ "[-0.0013615711124372242, -7.349896979829465e-6, -5.359923708254862e-5, -0.008582713096395626, 59.85516689032128]\n"
+ ]
+ }
+ ],
"source": [
"# Données pour la trajectoire durant le transfert\n",
"mutable struct Transfert\n",
@@ -1203,7 +1236,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 34,
"metadata": {},
"outputs": [],
"source": [
@@ -1226,9 +1259,50 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 35,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "ename": "LoadError",
+ "evalue": "MethodError: no method matching process_annotation(::Plots.Subplot{Plots.GRBackend}, ::Float64)\n\u001b[0mClosest candidates are:\n\u001b[0m process_annotation(::Plots.Subplot, ::Any, \u001b[91m::Any\u001b[39m, \u001b[91m::Any\u001b[39m) at ~/opt/miniconda3/envs/julia/share/julia/packages/Plots/GGa6i/src/components.jl:618\n\u001b[0m process_annotation(::Plots.Subplot, ::Any, \u001b[91m::Any\u001b[39m, \u001b[91m::Any\u001b[39m, \u001b[91m::Any\u001b[39m) at ~/opt/miniconda3/envs/julia/share/julia/packages/Plots/GGa6i/src/components.jl:618\n\u001b[0m process_annotation(::Plots.Subplot, \u001b[91m::Union{Symbol, Tuple, AbstractVector{Symbol}}\u001b[39m, \u001b[91m::Any\u001b[39m) at ~/opt/miniconda3/envs/julia/share/julia/packages/Plots/GGa6i/src/components.jl:632\n\u001b[0m ...",
+ "output_type": "error",
+ "traceback": [
+ "MethodError: no method matching process_annotation(::Plots.Subplot{Plots.GRBackend}, ::Float64)\n\u001b[0mClosest candidates are:\n\u001b[0m process_annotation(::Plots.Subplot, ::Any, \u001b[91m::Any\u001b[39m, \u001b[91m::Any\u001b[39m) at ~/opt/miniconda3/envs/julia/share/julia/packages/Plots/GGa6i/src/components.jl:618\n\u001b[0m process_annotation(::Plots.Subplot, ::Any, \u001b[91m::Any\u001b[39m, \u001b[91m::Any\u001b[39m, \u001b[91m::Any\u001b[39m) at ~/opt/miniconda3/envs/julia/share/julia/packages/Plots/GGa6i/src/components.jl:618\n\u001b[0m process_annotation(::Plots.Subplot, \u001b[91m::Union{Symbol, Tuple, AbstractVector{Symbol}}\u001b[39m, \u001b[91m::Any\u001b[39m) at ~/opt/miniconda3/envs/julia/share/julia/packages/Plots/GGa6i/src/components.jl:632\n\u001b[0m ...",
+ "",
+ "Stacktrace:",
+ " [1] _update_subplot_periphery(sp::Plots.Subplot{Plots.GRBackend}, anns::Vector{Any})",
+ " @ Plots ~/opt/miniconda3/envs/julia/share/julia/packages/Plots/GGa6i/src/args.jl:1928",
+ " [2] _update_subplot_args(plt::Plots.Plot{Plots.GRBackend}, sp::Plots.Subplot{Plots.GRBackend}, plotattributes_in::Dict{Symbol, Any}, subplot_index::Int64, remove_pair::Bool)",
+ " @ Plots ~/opt/miniconda3/envs/julia/share/julia/packages/Plots/GGa6i/src/args.jl:2061",
+ " [3] _subplot_setup(plt::Plots.Plot{Plots.GRBackend}, plotattributes::Dict{Symbol, Any}, kw_list::Vector{Dict{Symbol, Any}})",
+ " @ Plots ~/opt/miniconda3/envs/julia/share/julia/packages/Plots/GGa6i/src/pipeline.jl:290",
+ " [4] plot_setup!(plt::Plots.Plot{Plots.GRBackend}, plotattributes::Dict{Symbol, Any}, kw_list::Vector{Dict{Symbol, Any}})",
+ " @ Plots ~/opt/miniconda3/envs/julia/share/julia/packages/Plots/GGa6i/src/pipeline.jl:151",
+ " [5] recipe_pipeline!(plt::Any, plotattributes::Any, args::Any)",
+ " @ RecipesPipeline ~/opt/miniconda3/envs/julia/share/julia/packages/RecipesPipeline/kJPXU/src/RecipesPipeline.jl:87",
+ " [6] _plot!(plt::Plots.Plot, plotattributes::Any, args::Any)",
+ " @ Plots ~/opt/miniconda3/envs/julia/share/julia/packages/Plots/GGa6i/src/plot.jl:208",
+ " [7] #plot!#152",
+ " @ ~/opt/miniconda3/envs/julia/share/julia/packages/Plots/GGa6i/src/plot.jl:198 [inlined]",
+ " [8] #annotate!#525",
+ " @ ~/opt/miniconda3/envs/julia/share/julia/packages/Plots/GGa6i/src/Plots.jl:249 [inlined]",
+ " [9] panneau_information!(px::Plots.Plot{Plots.GRBackend}, subplot_current::Int64, time_current::Float64, i_current::Int64, x1::Vector{Float64}, x2::Vector{Float64}, v1::Vector{Float64}, v2::Vector{Float64}, θ::Vector{Float64}, F_max::Float64, tf_transfert::Float64, X1_orb_init::Vector{Float64}, X2_orb_init::Vector{Float64}, X1_orb_arr::Vector{Float64}, X2_orb_arr::Vector{Float64})",
+ " @ Main ./In[27]:63",
+ " [10] macro expansion",
+ " @ ./In[29]:202 [inlined]",
+ " [11] macro expansion",
+ " @ ~/opt/miniconda3/envs/julia/share/julia/packages/Plots/GGa6i/src/animation.jl:208 [inlined]",
+ " [12] animation(pre_transfert_data::PreTransfert, transfert_data::Transfert, post_transfert_data::PostTransfert; fps::Int64, nFrame::Int64)",
+ " @ Main ./In[29]:162",
+ " [13] top-level scope",
+ " @ In[35]:3",
+ " [14] eval",
+ " @ ./boot.jl:373 [inlined]",
+ " [15] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)",
+ " @ Base ./loading.jl:1196"
+ ]
+ }
+ ],
"source": [
"#nFrame = 10; fps = 1\n",
"nFrame = 2000; fps = 20\n",
@@ -1237,9 +1311,87 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 36,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "ffmpeg version 5.0.1 Copyright (c) 2000-2022 the FFmpeg developers\n",
+ " built with Apple clang version 13.0.0 (clang-1300.0.29.30)\n",
+ " configuration: --prefix=/usr/local/Cellar/ffmpeg/5.0.1_3 --enable-shared --enable-pthreads --enable-version3 --cc=clang --host-cflags= --host-ldflags= --enable-ffplay --enable-gnutls --enable-gpl --enable-libaom --enable-libbluray --enable-libdav1d --enable-libmp3lame --enable-libopus --enable-librav1e --enable-librist --enable-librubberband --enable-libsnappy --enable-libsrt --enable-libtesseract --enable-libtheora --enable-libvidstab --enable-libvmaf --enable-libvorbis --enable-libvpx --enable-libwebp --enable-libx264 --enable-libx265 --enable-libxml2 --enable-libxvid --enable-lzma --enable-libfontconfig --enable-libfreetype --enable-frei0r --enable-libass --enable-libopencore-amrnb --enable-libopencore-amrwb --enable-libopenjpeg --enable-libspeex --enable-libsoxr --enable-libzmq --enable-libzimg --disable-libjack --disable-indev=jack --enable-videotoolbox\n",
+ " libavutil 57. 17.100 / 57. 17.100\n",
+ " libavcodec 59. 18.100 / 59. 18.100\n",
+ " libavformat 59. 16.100 / 59. 16.100\n",
+ " libavdevice 59. 4.100 / 59. 4.100\n",
+ " libavfilter 8. 24.100 / 8. 24.100\n",
+ " libswscale 6. 4.100 / 6. 4.100\n",
+ " libswresample 4. 3.100 / 4. 3.100\n",
+ " libpostproc 56. 3.100 / 56. 3.100\n",
+ "Input #0, mov,mp4,m4a,3gp,3g2,mj2, from 'transfert-temps-min-original.mp4':\n",
+ " Metadata:\n",
+ " major_brand : isom\n",
+ " minor_version : 512\n",
+ " compatible_brands: isomiso2avc1mp41\n",
+ " encoder : Lavf58.76.100\n",
+ " Duration: 00:01:30.70, start: 0.000000, bitrate: 284 kb/s\n",
+ " Stream #0:0[0x1](und): Video: h264 (High 4:4:4 Predictive) (avc1 / 0x31637661), yuv444p(progressive), 1600x900, 283 kb/s, 20 fps, 20 tbr, 10240 tbn (default)\n",
+ " Metadata:\n",
+ " handler_name : VideoHandler\n",
+ " vendor_id : [0][0][0][0]\n",
+ "Stream mapping:\n",
+ " Stream #0:0 -> #0:0 (h264 (native) -> h264 (libx264))\n",
+ "Press [q] to stop, [?] for help\n",
+ "[libx264 @ 0x7fde1b529080] using cpu capabilities: MMX2 SSE2Fast SSSE3 SSE4.2 AVX FMA3 BMI2 AVX2\n",
+ "[libx264 @ 0x7fde1b529080] profile High, level 4.0, 4:2:0, 8-bit\n",
+ "[libx264 @ 0x7fde1b529080] 264 - core 164 r3095 baee400 - H.264/MPEG-4 AVC codec - Copyleft 2003-2022 - http://www.videolan.org/x264.html - options: cabac=1 ref=3 deblock=1:0:0 analyse=0x3:0x113 me=hex subme=7 psy=1 psy_rd=1.00:0.00 mixed_ref=1 me_range=16 chroma_me=1 trellis=1 8x8dct=1 cqm=0 deadzone=21,11 fast_pskip=1 chroma_qp_offset=-2 threads=6 lookahead_threads=1 sliced_threads=0 nr=0 decimate=1 interlaced=0 bluray_compat=0 constrained_intra=0 bframes=3 b_pyramid=2 b_adapt=1 b_bias=0 direct=1 weightb=1 open_gop=0 weightp=2 keyint=250 keyint_min=20 scenecut=40 intra_refresh=0 rc_lookahead=40 rc=crf mbtree=1 crf=23.0 qcomp=0.60 qpmin=0 qpmax=69 qpstep=4 ip_ratio=1.40 aq=1:1.00\n",
+ "Output #0, mp4, to 'transfert-temps-min.mp4':\n",
+ " Metadata:\n",
+ " major_brand : isom\n",
+ " minor_version : 512\n",
+ " compatible_brands: isomiso2avc1mp41\n",
+ " encoder : Lavf59.16.100\n",
+ " Stream #0:0(und): Video: h264 (avc1 / 0x31637661), yuv420p(tv, progressive), 1600x900, q=2-31, 20 fps, 10240 tbn (default)\n",
+ " Metadata:\n",
+ " handler_name : VideoHandler\n",
+ " vendor_id : [0][0][0][0]\n",
+ " encoder : Lavc59.18.100 libx264\n",
+ " Side data:\n",
+ " cpb: bitrate max/min/avg: 0/0/0 buffer size: 0 vbv_delay: N/A\n",
+ "frame= 1814 fps= 66 q=-1.0 Lsize= 3013kB time=00:01:30.55 bitrate= 272.6kbits/s speed=3.31x \n",
+ "video:2991kB audio:0kB subtitle:0kB other streams:0kB global headers:0kB muxing overhead: 0.723257%\n",
+ "[libx264 @ 0x7fde1b529080] frame I:8 Avg QP:15.52 size: 64104\n",
+ "[libx264 @ 0x7fde1b529080] frame P:544 Avg QP:20.28 size: 2445\n",
+ "[libx264 @ 0x7fde1b529080] frame B:1262 Avg QP:26.77 size: 966\n",
+ "[libx264 @ 0x7fde1b529080] consecutive B-frames: 4.1% 6.6% 8.1% 81.1%\n",
+ "[libx264 @ 0x7fde1b529080] mb I I16..4: 7.9% 64.6% 27.5%\n",
+ "[libx264 @ 0x7fde1b529080] mb P I16..4: 0.3% 0.4% 0.2% P16..4: 3.5% 0.7% 0.5% 0.0% 0.0% skip:94.4%\n",
+ "[libx264 @ 0x7fde1b529080] mb B I16..4: 0.0% 0.1% 0.0% B16..8: 4.0% 0.5% 0.1% direct: 0.1% skip:95.1% L0:46.8% L1:51.5% BI: 1.7%\n",
+ "[libx264 @ 0x7fde1b529080] 8x8 transform intra:56.6% inter:26.6%\n",
+ "[libx264 @ 0x7fde1b529080] coded y,uvDC,uvAC intra: 15.0% 5.8% 4.9% inter: 0.5% 0.3% 0.2%\n",
+ "[libx264 @ 0x7fde1b529080] i16 v,h,dc,p: 58% 39% 3% 0%\n",
+ "[libx264 @ 0x7fde1b529080] i8 v,h,dc,ddl,ddr,vr,hd,vl,hu: 56% 5% 39% 0% 0% 0% 0% 0% 0%\n",
+ "[libx264 @ 0x7fde1b529080] i4 v,h,dc,ddl,ddr,vr,hd,vl,hu: 48% 19% 16% 3% 3% 3% 3% 3% 3%\n",
+ "[libx264 @ 0x7fde1b529080] i8c dc,h,v,p: 91% 6% 3% 0%\n",
+ "[libx264 @ 0x7fde1b529080] Weighted P-Frames: Y:0.0% UV:0.0%\n",
+ "[libx264 @ 0x7fde1b529080] ref P L0: 67.7% 2.4% 19.2% 10.8%\n",
+ "[libx264 @ 0x7fde1b529080] ref B L0: 78.1% 17.9% 4.0%\n",
+ "[libx264 @ 0x7fde1b529080] ref B L1: 94.6% 5.4%\n",
+ "[libx264 @ 0x7fde1b529080] kb/s:270.10\n"
+ ]
+ },
+ {
+ "data": {
+ "text/plain": [
+ "Process(`\u001b[4mffmpeg\u001b[24m \u001b[4m-y\u001b[24m \u001b[4m-i\u001b[24m \u001b[4mtransfert-temps-min-original.mp4\u001b[24m \u001b[4m-vf\u001b[24m \u001b[4mformat=yuv420p\u001b[24m \u001b[4mtransfert-temps-min.mp4\u001b[24m`, ProcessExited(0))"
+ ]
+ },
+ "execution_count": 36,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
"source": [
"convert = `ffmpeg -y -i transfert-temps-min-original.mp4 -vf format=yuv420p transfert-temps-min.mp4`\n",
"if nFrame > 1\n",