diff --git a/.gitignore b/.gitignore
index c4ea72c..282da38 100644
--- a/.gitignore
+++ b/.gitignore
@@ -30,3 +30,6 @@ Manifest.toml
 **/Manifest.toml
 # ...except the one inside docs/src/assets
 !docs/src/assets/Manifest.toml
+
+# Ignore references 
+references/
\ No newline at end of file
diff --git a/save/3body_L1.jl b/save/3body_L1.jl
new file mode 100644
index 0000000..7829aa8
--- /dev/null
+++ b/save/3body_L1.jl
@@ -0,0 +1,79 @@
+using OptimalControl
+using NLPModelsIpopt
+using OrdinaryDiffEq
+using Plots
+#using MINPACK
+using ForwardDiff
+using LinearAlgebra
+
+ε = 0.01 
+
+## Problem definition
+
+G = 6.61e10-11
+m₁ = 5.972e24          # Mass of the Earth
+m₂ = 7.348e22          # Mass of the Moon
+μ = m₂ / (m₁ + m₂)     # Ratio of primary masses
+M = 1000               # Mass of the satellite (kg)
+# !!!!!!!!!!!!
+γ_max = 1              # Poussée max
+T_max = 1              # Maximal thrust
+# !!!!!!!!!!!!
+
+# 1st body is located at (-μ, 0)
+# 2nd body is located at (1-μ, 0)
+
+# Kepler law speed in circular movement: v² = G * M / d
+d_earth = 0.1
+d_moon = 0.1
+
+v_earth_circular = sqrt(G * m₁ / d_earth)
+v_moon_circular = sqrt(G * m₂ / d_moon)
+
+q0 = [-μ - d_earth, 0.0]   # Initial position (km) (e.g. GEO Earth)
+v0 = [0, v_earth_circular]       # Initial velocity (km/s)
+qf = [1 - μ + d_moon, 0.0]    # Target position (km) (e.g. GEO Moon)
+vf = [0.0, -v_moon_circular]       # Target velocity (km/s)
+x0 = vcat(q0, v0)
+xf = vcat(qf, vf)
+
+tf_guess = 20000.0          # Time guess in seconds (5.5 hours)
+
+# =============== INITIAL GUESS ===============
+x(t) = x0 + (xf - x0) * t / tf_guess
+u(t) = [0.01, 0.01] # initial control guess
+
+
+nlp_init = (state = x, control = u, variable = tf_guess)
+
+# =============== OCP DEFINITION ===============
+@def ocp begin
+    tf ∈ R, variable
+    t ∈ [0, tf], time
+    x ∈ R⁴, state
+    u ∈ R², control
+    x(0) == x0
+    x(tf) == xf
+    r₁₃ = sqrt((x₁(t) + μ)^2 + x₂(t)^2) + 1e-9  # Avoid division by zero
+    r₂₃ = sqrt((x₁(t) - 1 + μ)^2 + x₂(t)^2) + 1e-9  # Avoid division by zero
+    F0 = [  x₃(t), 
+            x₄(t), 
+            2*x₄(t) + x₁(t) - ((1 - μ) / r₁₃^3)*(x₁(t) + μ) - (μ / r₂₃^3)*(x₁(t) - 1 + μ),
+            -2*x₃(t) + x₂(t) - ((1 - μ) / r₁₃^3)*x₂(t) - (μ / r₂₃^3)*x₂(t) ]
+    F1 = [0, 0, 1, 0]
+    F2 = [0, 0, 0, 1]
+    ẋ(t) == F0 + (T_max/M)*F1*u1(t) + (T_max/M)*F2*u₂(t)
+    u₁(t)^2 + u₂(t)^2 ≤ 1
+    # Fonction objectif avec régularisation pour éviter les problèmes de domaine
+    u_norm = sqrt(u₁(t)^2 + u₂(t)^2) + 1e-12  # Évite log(0)
+    u_complement = max(1e-12, 1 - u_norm)  # Évite log(≤0)
+    ∫(u_norm - (1-ε)*log(u_norm) - (1-ε)*log(u_complement)) → min
+end
+
+
+# =============== SOLVE ===============
+
+nlp_sol = OptimalControl.solve(ocp; init=nlp_init, grid_size=50, print_level=3)
+plot(nlp_sol, vars=:x, title="3-body Controlled Transfer")
+
+
diff --git a/save/three_body_time.jl b/save/three_body_time.jl
new file mode 100644
index 0000000..91d27c2
--- /dev/null
+++ b/save/three_body_time.jl
@@ -0,0 +1,53 @@
+using OptimalControl
+using CTBase
+using NLPModelsIpopt
+using OrdinaryDiffEq
+using Plots
+using ForwardDiff
+using LinearAlgebra
+
+## Problem definition (CR3BP)
+
+μ = 0.0121505856              # Earth-Moon reduced mass ratio
+ϵ = 1.0                       # Maximum thrust magnitude (normalized)
+q0 = [1.0 - μ - 0.01, 0.0]    # Initial position (near Earth)
+v0 = [0.0, 0.0]               # Initial velocity
+qf = [μ + 0.05, 0.0]          # Final position (near Moon)
+vf = [0.0, 0.0]               # Final velocity
+x0 = vcat(q0, v0)
+xf = vcat(qf, vf)
+tf_guess = 5.0               # Initial guess for final time
+
+# Gradient of the effective potential
+function grad_omega(q)
+    x, y = q
+    r1 = sqrt((x + μ)^2 + y^2) + 1e-6
+    r2 = sqrt((x - 1 + μ)^2 + y^2) + 1e-6
+    dΩdx = x - (1 - μ)*(x + μ)/r1^3 - μ*(x - 1 + μ)/r2^3
+    dΩdy = y - (1 - μ)*y/r1^3 - μ*y/r2^3
+    return [dΩdx, dΩdy]
+end
+
+
+x(t) = x0 + (xf - x0) * t / tf_guess
+u(t) = [0.01, 0.01] 
+nlp_init = (state = x, control = u, variable = tf_guess)
+
+@def ocp begin
+    tf ∈ R, variable
+    t ∈ [0, tf], time
+    x = (q1, q2, v1, v2) ∈ R⁴, state
+    u ∈ R², control
+    x(0) == x0
+    x(tf) == xf
+    r1 = sqrt((q1(t) + μ)^2 + q2(t)^2) + 1e-6
+    r2 = sqrt((q1(t) - 1 + μ)^2 + q2(t)^2) + 1e-6
+    dΩ1 = q1(t) - (1 - μ)*(q1(t) + μ)/r1^3 - μ*(q1(t) - 1 + μ)/r2^3
+    dΩ2 = q2(t) - (1 - μ)*q2(t)/r1^3 - μ*q2(t)/r2^3
+    ẋ(t) == [v1(t), v2(t), 2v2(t) + dΩ1 + ϵ*u₁(t), -2v1(t) + dΩ2 + ϵ*u₂(t)]
+    u₁(t)^2 + u₂(t)^2 ≤ 1
+    tf → min
+end
+
+nlp_sol = OptimalControl.solve(ocp; init=nlp_init, grid_size=100, print_level=3)
+plot(nlp_sol, vars=:x, title="CR3BP - Minimum Time Transfer")
diff --git a/save/two_body.jl b/save/two_body.jl
new file mode 100644
index 0000000..27412ea
--- /dev/null
+++ b/save/two_body.jl
@@ -0,0 +1,138 @@
+using OptimalControl
+using NLPModelsIpopt
+using OrdinaryDiffEq
+using Plots
+using MINPACK
+using ForwardDiff
+using LinearAlgebra
+
+Tmax = 60                                  # Maximum thrust in Newtons
+cTmax = 3600^2 / 1e6; T = Tmax * cTmax     # Conversion from Newtons to kg x Mm / h²
+mass0 = 1500                               # Initial mass of the spacecraft
+β = 1.42e-02                               # Engine specific impulsion
+μ = 5165.8620912                           # Earth gravitation constant
+P0 = 11.625                                # Initial semilatus rectum
+ex0, ey0 = 0.75, 0                         # Initial eccentricity
+hx0, hy0 = 6.12e-2, 0                      # Initial ascending node and inclination
+L0 = π                                     # Initial longitude
+Pf = 42.165                                # Final semilatus rectum
+exf, eyf = 0, 0                            # Final eccentricity
+hxf, hyf = 0, 0                            # Final ascending node and inclination
+ε = 1e-1                                   # Regularization parameter for logarithmic barrier
+
+asqrt(x; ε=1e-9) = sqrt(sqrt(x^2 + ε^2))
+function F0(x)
+    P, ex, ey, hx, hy, L = x
+    pdm = asqrt(P / μ)
+    cl = cos(L)
+    sl = sin(L)
+    w = 1 + ex * cl + ey * sl
+    F = zeros(eltype(x), 6)
+    F[6] = w^2 / (P * pdm)
+    return F
+end
+
+function F1(x)
+    P, ex, ey, hx, hy, L = x
+    pdm = asqrt(P / μ)
+    cl = cos(L)
+    sl = sin(L)
+    F = zeros(eltype(x), 6)
+    F[2] = pdm *   sl
+    F[3] = pdm * (-cl)
+    return F
+end
+
+function F2(x)
+    P, ex, ey, hx, hy, L = x
+    pdm = asqrt(P / μ)
+    cl = cos(L)
+    sl = sin(L)
+    w = 1 + ex * cl + ey * sl
+    F = zeros(eltype(x), 6)
+    F[1] = pdm * 2 * P / w
+    F[2] = pdm * (cl + (ex + cl) / w)
+    F[3] = pdm * (sl + (ey + sl) / w)
+    return F
+end
+
+function F3(x)
+    P, ex, ey, hx, hy, L = x
+    pdm = asqrt(P / μ)
+    cl = cos(L)
+    sl = sin(L)
+    w = 1 + ex * cl + ey * sl
+    pdmw = pdm / w
+    zz = hx * sl - hy * cl
+    uh = (1 + hx^2 + hy^2) / 2
+    F = zeros(eltype(x), 6)
+    F[2] = pdmw * (-zz * ey)
+    F[3] = pdmw *   zz * ex
+    F[4] = pdmw *   uh * cl
+    F[5] = pdmw *   uh * sl
+    F[6] = pdmw *   zz
+    return F
+end
+
+tf = 15
+Lf = 3π
+x0 = [P0, ex0, ey0, hx0, hy0, L0]
+xf = [Pf, exf, eyf, hxf, hyf, Lf]
+x(t) = x0 + (xf - x0) * t / tf
+u = [0.1, 0.5, 0.]
+nlp_init = (state=x, control=u, variable=tf)
+
+function min_tf()
+    @def ocp begin
+        tf ∈ R, variable
+        t ∈ [0, tf], time
+        x = (P, ex, ey, hx, hy, L) ∈ R⁶, state
+        u ∈ R³, control
+        x(0) == x0
+        x[1:5](tf) == xf[1:5]
+        mass = mass0 - β * T * t
+        ẋ(t) == F0(x(t)) + T / mass * (u₁(t) * F1(x(t)) + u₂(t) * F2(x(t)) + u₃(t) * F3(x(t)))
+        u₁(t)^2 + u₂(t)^2 + u₃(t)^2 ≤ 1
+        tf → min
+    end
+    return ocp
+end
+
+function min_conso(e)
+    @def ocp begin
+        t ∈ [0, tf], time
+        x = (P, ex, ey, hx, hy, L) ∈ R⁶, state
+        u ∈ R³, control
+        x(0) == x0
+        x[1:5](tf) == xf[1:5]
+        mass = mass0 - β * T * t
+        u_norm = sqrt(u₁(t)^2 + u₂(t)^2 + u₃(t)^2)
+        ẋ(t) == F0(x(t)) + T / mass * (u₁(t) * F1(x(t)) + u₂(t) * F2(x(t)) + u₃(t) * F3(x(t)))
+        1e-3 ≤ u_norm^2 ≤ 1
+        u_g = max(1e-10, 1 - u_norm)
+        ∫(u_norm - e * (log(u_norm) + log(u_g))) → min
+    end
+    return ocp
+end
+
+ocp1 = min_tf()
+sol = solve(ocp1; init=nlp_init, grid_size=100)
+Tmax = 100.0
+T = Tmax * cTmax
+tf = variable(sol)
+ocp = min_conso(ε)
+nlp_init = (state = state(sol), control = control(sol))
+
+eps = 1e-1
+i = nlp_init
+s = sol
+ocp = min_conso(eps)
+global s = solve(ocp; init=i, grid_size=500)
+global i = s
+
+
+plot(s; label="ε = $eps")
+
+t = time_grid(s)
+u = control(s)
+plot(t, norm∘u; label="‖u‖", xlabel="t")
\ No newline at end of file
diff --git a/save/two_body_JuMP.jl b/save/two_body_JuMP.jl
new file mode 100644
index 0000000..e35a036
--- /dev/null
+++ b/save/two_body_JuMP.jl
@@ -0,0 +1,227 @@
+using JuMP
+using Ipopt
+using Plots
+using LinearAlgebra
+
+Tmax = 60                                  # Maximum thrust in Newtons
+cTmax = 3600^2 / 1e6; T = Tmax * cTmax     # Conversion from Newtons to kg x Mm / h²
+mass0 = 1500                               # Initial mass of the spacecraft
+β = 1.42e-02                               # Engine specific impulse
+μ = 5165.8620912                           # Earth gravitational constant
+P0 = 11.625                                # Initial semilatus rectum
+ex0, ey0 = 0.75, 0                         # Initial eccentricity
+hx0, hy0 = 6.12e-2, 0                      # Initial ascending node and inclination
+L0 = π                                     # Initial longitude
+Pf = 42.165                                # Final semilatus rectum
+exf, eyf = 0, 0                            # Final eccentricity
+hxf, hyf = 0, 0                            # Final ascending node and inclination
+ε = 1e-1                                   # Regularization parameter for logarithmic barrier
+tf = 20
+
+#1 reduce epsilon 
+#2 increase tf u
+#reduce Tmax
+
+N = 100  # Number of discretization points
+dt = tf / (N-1)
+
+# ===================== Initial and Final States =====================
+Lf = 3π                                      # Estimate of final longitude
+x0 = [P0, ex0, ey0, hx0, hy0, L0]            # Initial state
+xf = [Pf, exf, eyf, hxf, hyf, Lf]            # Final state
+
+# Time grid
+t_grid = range(0, tf, length=N)
+
+# ===================== JuMP Model =====================
+model = JuMP.Model(Ipopt.Optimizer)
+JuMP.set_silent(model)
+
+# State variables [P, ex, ey, hx, hy, L] at N time points
+JuMP.@variable(model, P[1:N])
+JuMP.@variable(model, ex[1:N])
+JuMP.@variable(model, ey[1:N])
+JuMP.@variable(model, hx[1:N])
+JuMP.@variable(model, hy[1:N])
+JuMP.@variable(model, L[1:N])
+
+# Control variables [u1, u2, u3] at N time points
+JuMP.@variable(model, u1[1:N])
+JuMP.@variable(model, u2[1:N])
+JuMP.@variable(model, u3[1:N])
+
+# Auxiliary variable for control norm
+JuMP.@variable(model, u_norm[1:N] >= 1e-6)
+
+# ===================== Constraints =====================
+
+# Initial conditions
+JuMP.@constraint(model, P[1] == x0[1])
+JuMP.@constraint(model, ex[1] == x0[2])
+JuMP.@constraint(model, ey[1] == x0[3])
+JuMP.@constraint(model, hx[1] == x0[4])
+JuMP.@constraint(model, hy[1] == x0[5])
+JuMP.@constraint(model, L[1] == x0[6])
+
+# Final conditions (only the first 5 states)
+JuMP.@constraint(model, P[N] == xf[1])
+JuMP.@constraint(model, ex[N] == xf[2])
+JuMP.@constraint(model, ey[N] == xf[3])
+JuMP.@constraint(model, hx[N] == xf[4])
+JuMP.@constraint(model, hy[N] == xf[5])
+
+# Constraints on control norm
+for i in 1:N
+    JuMP.@constraint(model, u1[i]^2 + u2[i]^2 + u3[i]^2 <= u_norm[i]^2)
+    JuMP.@constraint(model, 1e-3 <= u_norm[i]^2 <= 1)
+    JuMP.@constraint(model, u_norm[i] >= 1e-6)  # Avoid u_norm = 0 for log
+    JuMP.@constraint(model, u_norm[i] <= 1.0 - 1e-6)  # Avoid u_norm = 1 for log
+end
+
+# Dynamics (Euler discretization) - CORRECTION: direct calculation in JuMP
+for i in 1:N-1
+    # Mass at point i
+    mass_i = mass0 - β * T * t_grid[i]
+    
+    # State variables at point i
+    P_i = P[i]
+    ex_i = ex[i]
+    ey_i = ey[i]
+    hx_i = hx[i]
+    hy_i = hy[i]
+    L_i = L[i]
+    
+    # Control variables at point i
+    u1_i = u1[i]
+    u2_i = u2[i]
+    u3_i = u3[i]
+    
+    # Direct calculation of Gauss equations in JuMP
+    # Smoothed sqrt to avoid differentiation problems
+    asqrt_P = sqrt(sqrt(P_i^2 + 1e-18))  # P/μ with smoothed sqrt
+    pdm = asqrt_P / sqrt(μ)
+    
+    cl = cos(L_i)
+    sl = sin(L_i)
+    w = 1 + ex_i * cl + ey_i * sl
+    
+    # F0 contribution
+    dL_dt_F0 = w^2 / (P_i * pdm)
+    
+    # F1 contribution
+    dex_dt_F1 = pdm * sl
+    dey_dt_F1 = pdm * (-cl)
+    
+    # F2 contribution
+    dP_dt_F2 = pdm * 2 * P_i / w
+    dex_dt_F2 = pdm * (cl + (ex_i + cl) / w)
+    dey_dt_F2 = pdm * (sl + (ey_i + sl) / w)
+    
+    # F3 contribution
+    pdmw = pdm / w
+    zz = hx_i * sl - hy_i * cl
+    uh = (1 + hx_i^2 + hy_i^2) / 2
+    
+    dex_dt_F3 = pdmw * (-zz * ey_i)
+    dey_dt_F3 = pdmw * zz * ex_i
+    dhx_dt_F3 = pdmw * uh * cl
+    dhy_dt_F3 = pdmw * uh * sl
+    dL_dt_F3 = pdmw * zz
+    
+    # Total dynamics: ẋ = F0 + T/mass * (u1*F1 + u2*F2 + u3*F3)
+    thrust_factor = T / mass_i
+    
+    dP_dt = thrust_factor * u2_i * dP_dt_F2
+    dex_dt = thrust_factor * (u1_i * dex_dt_F1 + u2_i * dex_dt_F2 + u3_i * dex_dt_F3)
+    dey_dt = thrust_factor * (u1_i * dey_dt_F1 + u2_i * dey_dt_F2 + u3_i * dey_dt_F3)
+    dhx_dt = thrust_factor * u3_i * dhx_dt_F3
+    dhy_dt = thrust_factor * u3_i * dhy_dt_F3
+    dL_dt = dL_dt_F0 + thrust_factor * u3_i * dL_dt_F3
+    
+    # Euler discretization
+    JuMP.@constraint(model, P[i+1] == P[i] + dt * dP_dt)
+    JuMP.@constraint(model, ex[i+1] == ex[i] + dt * dex_dt)
+    JuMP.@constraint(model, ey[i+1] == ey[i] + dt * dey_dt)
+    JuMP.@constraint(model, hx[i+1] == hx[i] + dt * dhx_dt)
+    JuMP.@constraint(model, hy[i+1] == hy[i] + dt * dhy_dt)
+    JuMP.@constraint(model, L[i+1] == L[i] + dt * dL_dt)
+end
+
+# ===================== Objective function with logarithmic barrier =====================
+# Objective: ∫(u_norm - ε * (log(u_norm) + log(1 - u_norm))) dt
+JuMP.@objective(model, Min, 
+    sum(u_norm[i] - ε * (log(u_norm[i]) + log(max(1e-10, 1 - u_norm[i]))) 
+        for i in 1:N) * dt)
+
+# ===================== Initialization =====================
+for i in 1:N
+    t_frac = (i-1) / (N-1)
+    
+    # State initialization (linear interpolation)
+    JuMP.set_start_value(P[i], x0[1] + t_frac * (xf[1] - x0[1]))
+    JuMP.set_start_value(ex[i], x0[2] + t_frac * (xf[2] - x0[2]))
+    JuMP.set_start_value(ey[i], x0[3] + t_frac * (xf[3] - x0[3]))
+    JuMP.set_start_value(hx[i], x0[4] + t_frac * (xf[4] - x0[4]))
+    JuMP.set_start_value(hy[i], x0[5] + t_frac * (xf[5] - x0[5]))
+    JuMP.set_start_value(L[i], x0[6] + t_frac * (xf[6] - x0[6]))
+    
+    # Control initialization
+    JuMP.set_start_value(u1[i], 0.1)
+    JuMP.set_start_value(u2[i], 0.5)
+    JuMP.set_start_value(u3[i], 0.0)
+    JuMP.set_start_value(u_norm[i], 0.5)
+end
+
+# ===================== Solve =====================
+println("Solving JuMP problem...")
+JuMP.optimize!(model)
+
+if JuMP.termination_status(model) == JuMP.MOI.LOCALLY_SOLVED
+    println("Solution found!")
+else
+    println("Convergence issue: ", JuMP.termination_status(model))
+end
+
+# ===================== Extract results =====================
+P_opt = JuMP.value.(P)
+ex_opt = JuMP.value.(ex)
+ey_opt = JuMP.value.(ey)
+hx_opt = JuMP.value.(hx)
+hy_opt = JuMP.value.(hy)
+L_opt = JuMP.value.(L)
+u1_opt = JuMP.value.(u1)
+u2_opt = JuMP.value.(u2)
+u3_opt = JuMP.value.(u3)
+u_norm_opt = JuMP.value.(u_norm)
+
+# ===================== Plots =====================
+
+# Plot states
+plt_states = plot(layout=(2,3), size=(1200, 600))
+plot!(plt_states[1], t_grid, P_opt, label="P", title="P", linewidth=2)
+plot!(plt_states[2], t_grid, ex_opt, label="ex", title="ex", linewidth=2)
+plot!(plt_states[3], t_grid, ey_opt, label="ey", title="ey", linewidth=2)
+plot!(plt_states[4], t_grid, hx_opt, label="hx", title="hx", linewidth=2)
+plot!(plt_states[5], t_grid, hy_opt, label="hy", title="hy", linewidth=2)
+plot!(plt_states[6], t_grid, L_opt, label="L", title="L", linewidth=2)
+display(plt_states)
+
+# Plot controls
+plt_control = plot(layout=(2,2), size=(800, 600))
+plot!(plt_control[1], t_grid, u1_opt, label="u1", title="u1", linewidth=2)
+plot!(plt_control[2], t_grid, u2_opt, label="u2", title="u2", linewidth=2)
+plot!(plt_control[3], t_grid, u3_opt, label="u3", title="u3", linewidth=2)
+plot!(plt_control[4], t_grid, u_norm_opt, label="||u||", title="||u||", linewidth=2)
+display(plt_control)
+
+# ===================== Results =====================
+println("\n========== JUMP RESULTS ==========")
+println("Status: ", JuMP.termination_status(model))
+println("Transfer time: ", tf, " h")
+println("Objective value: ", JuMP.objective_value(model))
+println("Number of points: ", N)
+println("Max control norm: ", maximum(u_norm_opt))
+println("Min control norm: ", minimum(u_norm_opt))
+println("Initial states: ", x0)
+println("Final states: ", [P_opt[end], ex_opt[end], ey_opt[end], hx_opt[end], hy_opt[end], L_opt[end]])
+println("Final error: ", norm([P_opt[end], ex_opt[end], ey_opt[end], hx_opt[end], hy_opt[end]] - xf[1:5]))