Merge branch 'master' into gh-pages

ode/cr_3-body_problem/julia/3-body_orbit.jl

@@ -4,20 +4,25 @@ using Plots
 
 include("config.jl")
 
+r_1 = u -> hypot(u[1] + mu, u[2], u[3])
+r_2 = u -> hypot(u[1] - 1 + mu, u[2], u[3])
+
 # Numerical Methods for Ordinary Differential Equations
 # J. C. Butcher
 # ISBN: 9781119121503
-function r3body(du, u, p, t)
+function cr3body(du, u, p, t)
     du[1] = u[4]
     du[2] = u[5]
     du[3] = u[6]
-    du[4] = 2*u[5] + u[1] - mu*(u[1] + mu - 1)/hypot(u[2], u[3], u[1] + mu - 1)^3 - (1 - mu)*(u[1] + mu)/hypot(u[2], u[3], u[1] + mu)^3
-    du[5] = -2*u[4] + u[2] - mu*u[2]/hypot(u[2], u[3], u[1] + mu - 1)^3 - (1 - mu)*u[2]/hypot(u[2], u[3], u[1] + mu)^3
-    du[6] = - mu*u[3]/hypot(u[2], u[3], u[1] + mu - 1)^3 - (1 - mu)*u[3]/hypot(u[2], u[3], u[1] + mu)^3
+    du[4] = 2*u[5] + u[1] - (1 - mu)*(u[1] + mu)/r_1(u)^3 - mu*(u[1] + mu - 1)/r_2(u)^3
+    du[5] = -2*u[4] + u[2] - (1 - mu)*u[2]/r_1(u)^3 - mu*u[2]/r_2(u)^3
+    du[6] =  -(1 - mu)*u[3]/r_1(u)^3 - mu*u[3]/r_2(u)^3
 end
 
+jacobi(u) = 0.5*(u[4]^2 + u[5]^2 + u[6]^2) - (1 - mu)/r_1(u) - mu/r_2(u) - 0.5*((1 - mu)*r_1(u)^2 + mu*r_2(u)^2)
+
 t_step = t_period/N
-prob = ODEProblem(r3body, u0, (0.0, t_period))
+prob = ODEProblem(cr3body, u0, (0.0, t_period))
 sol = solve(prob, Vern7(), adaptive=false, dt=t_step)
 
 println("mu = ", mu)
@@ -63,5 +68,24 @@ plot!(plotOrbit, x_foreground_color_axis=:lightgrey, y_foreground_color_axis=:li
 plot!(plotOrbit, x_foreground_color_text=:lightgrey, y_foreground_color_text=:lightgrey)
 plot!(plotOrbit, x_foreground_color_border=:lightgrey, y_foreground_color_border=:lightgrey)
 
-gui(), readline()
-savefig(plotOrbit, "orbit.png")
+#gui(), readline()
+savefig(plotOrbit, "orbit.png")
+
+inv0 = jacobi(u0)
+println("jacobi constant = ", inv0)
+inv = map(jacobi, sol.u)
+
+plotError = scatter(
+    sol.t, (inv .- inv0)./inv0,
+    xlabel = "time",
+    ylabel = "error",
+    plot_title = "Relative invariant error in restricted 3-body problem",
+    markersize = 0.8,
+    alpha = 0.15,
+    legend = false,
+    dpi = 300,
+    annotation = (t_period - 3.0, 0.0, text("\$\\mu= $(mu)\$\n\$inv = $(inv0)\$", :crimson, :bottom, 4))
+)
+
+# gui(), readline()
+savefig(plotError, "error.png")

ode/cr_3-body_problem/julia/config.jl

@@ -1,13 +1,13 @@
 # For mu = 1/81.45 (Earth-Moon system) there are some quasi-periodic orbits
 
 # t_period = 17.06521656015796
-# u0 = [0.994, 0.0, 0.0, 0.0, −2.0015851063790825224, 0.0]
+# u0 = [0.994, 0.0, 0.0, 0.0, -2.0015851063790825224, 0.0]
 
 # t_period = 19.14045706162071
-# u0 = [0.87978, 0.0, 0.0, 0.0, −0.3797, 0.0]
+# u0 = [0.87978, 0.0, 0.0, 0.0, -0.3797, 0.0]
 
-mu = 1/200
-N = 32000
+mu = 1/81.45
+N = 6000
 
-t_period = 960
+t_period = 480
 u0 = [0.5, 0.85, 0.0, 0.0, 0.0, 0.0]

ode/cr_3-body_problem/julia/energy_contour.jl

@@ -0,0 +1,33 @@
+using Plots; pythonplot(size = (300, 300))
+
+include("config.jl")
+
+function energy(x, y)
+    r_1 = sqrt((x + mu)^2 + y^2)
+    r_2 = sqrt((x + mu - 1)^2 + y^2)
+    e = -(1 - mu)/r_1 - mu/r_2 - 0.5*((1 - mu)*r_1^2 + mu*r_2^2)
+    if e < -1.76
+        return -1.76
+    end
+    return e
+end
+
+limit = 1.5
+x = range(-limit, limit, length = 720)
+y = range(-limit, limit, length = 720)
+
+z = @. energy(x', y)
+println("mu = ", mu)
+
+plotContour = contour(x, y, z,
+    levels = range(-1.75, -1.5, length = 120),
+    lw = 0.4,
+    alpha = 0.95,
+    aspect_ratio = 1,
+    dpi = 300,
+    annotation = (0.15, -0.25, text("\$\\mu= $(mu)\$", :crimson, :bottom, 4))
+)
+
+gui(), readline()
+# ERROR: unable to save file
+# savefig(plotContour, "pseudo-potential.png")

ode/cr_3-body_problem/octave/config.m

@@ -0,0 +1,7 @@
+function cfg = config
+    cfg.mu = 1/81.45;
+    cfg.x_init = [0.87978; 0; 0; 0; -0.3797; 0];
+    cfg.t_max = 19.14045706162071;
+    cfg.N = 6000;
+end
+

ode/cr_3-body_problem/octave/cr3bp_orbit.m

@@ -0,0 +1,90 @@
+% Circular restricted 3-body problem
+function cr3bp_orbit
+
+    % Read config params
+    cfg = config;
+
+    % Equations of motion in corotating frame
+    function xdot = f(x, t)
+        mu = cfg.mu;
+        r1 = hypot((x(1) + mu), x(2), x(3));
+        r2 = hypot((x(1) - 1 + mu), x(2), x(3));
+        xdot(1) = x(4);
+        xdot(2) = x(5);
+        xdot(3) = x(6);
+        xdot(4) = 2*x(5) + x(1) - (1 - mu)*(x(1) + mu)/r1^3 - mu*(x(1) - 1 + mu)/r2^3;
+        xdot(5) = -2*x(4) + x(2) - (1 - mu)*x(2)/r1^3 - mu*x(2)/r2^3;
+        xdot(6) = -(1 - mu)*x(3)/r1^3 - mu*x(3)/r2^3;
+    endfunction
+
+    % Constant of motion
+    function e = energy(x)
+        mu = cfg.mu;
+        r1 = hypot((x(1) + mu), x(2), x(3));
+        r2 = hypot((x(1) - 1 + mu), x(2), x(3));
+        e = 0.5*(x(4)^2 + x(5)^2 + x(6)^2) - (1 - mu)/r1 - mu/r2 - 0.5*((1 - mu)*r1^2 + mu*r2^2);
+    endfunction
+
+    % Plot orbit
+    function configuration_plot(x)
+        gap = 0.1;
+        marker_size = 0.5;
+
+        hf = figure();
+        hold on
+        scatter(x(:, [1]), x(:, [2]), marker_size, 'r');
+        axis equal
+        axis([-1 - gap, 1 + gap, -1 - gap, 1 + gap]);
+        set(gca, 'fontsize', 11);
+        grid on
+        title "Circular restricted 3-body problem"
+        xlabel "x"
+        ylabel "y"
+        plot(-cfg.mu, 0, 'bo', 'MarkerFaceColor', 'b', 'MarkerSize', 8)
+        plot(1 - cfg.mu, 0, 'bo', 'MarkerFaceColor', '#FF8C00', 'MarkerSize', 4)
+
+        glim = axis();
+        text(
+            glim(1) + 0.05*(glim(2) - glim(1)),
+            glim(4) - 0.05*(glim(4) - glim(3)),
+            sprintf("\\mu = %d\nt = [0, %d]", cfg.mu, cfg.t_max)
+        );
+        print(hf, "orbit.jpg", "-S960,720");
+        print(hf, "orbit.eps", "-S960,720");
+    endfunction
+
+    % Plot relative error in Jacobi constant
+    function invariant_plot(x, t)
+        marker_size = 1.5;
+
+        e0 = energy(cfg.x_init);
+        m = 1;
+        for el = x'
+            sample(m).e = energy(el) - e0;
+            m++;
+        endfor
+
+        hf = figure();
+        scatter(t, [sample.e], marker_size);
+        set(gca, 'fontsize', 11);
+        grid on
+        title "Relative error in Jacobi constant"
+        xlabel "Time"
+        ylabel "Relative error"
+        legend("Jacobi constant");
+
+        glim = axis();
+        text(
+            glim(1) + 0.05*(glim(2) - glim(1)),
+            glim(4) - 0.05*(glim(4) - glim(3)),
+            sprintf("\\mu = %d\nJ = %d", cfg.mu, e0)
+        );
+        print(hf, "error.jpg", "-S960,720");
+        print(hf, "error.eps", "-S960,720");
+    endfunction
+
+    x = lsode("f", cfg.x_init, (t = linspace (0, cfg.t_max, cfg.N)'));
+    configuration_plot(x);
+    invariant_plot(x, t);
+end
+

ode/cr_3-body_problem/python/energy_contour.py

@@ -0,0 +1,34 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+from matplotlib.patches import Circle
+
+import config as cfg
+
+r_1 = lambda x, y: np.sqrt((x + cfg.mu)**2 + y**2)
+r_2 = lambda x, y: np.sqrt((x - 1 + cfg.mu)**2 + y**2)
+
+limit = 1.5
+delta = 0.005
+x = np.arange(-limit, limit, delta)
+y = np.arange(-limit, limit, delta)
+X, Y = np.meshgrid(x, y)
+
+Z = - (1 - cfg.mu)/r_1(X, Y) - cfg.mu/r_2(X, Y) - 0.5*((1 - cfg.mu)*r_1(X, Y)**2 + cfg.mu*r_2(X, Y)**2)
+
+fig = plt.figure(figsize=(8, 8))
+plt.title("Potential energy in circular restricted 3-body problem")
+plt.axis("equal")
+plt.xlabel("x")
+plt.ylabel("y")
+plt.grid(True, linestyle = '--', linewidth = 0.2)
+
+plt.contour(X, Y, Z, levels=np.arange(-1.75, -1.5, 0.005), linewidths=0.1)
+body1 = Circle((-cfg.mu, 0), radius=0.03, color='darkslateblue', alpha=0.5)
+body2 = Circle((1 - cfg.mu, 0), radius=0.01, color='crimson', alpha=0.5)
+plt.gca().add_artist(body1)
+plt.gca().add_artist(body2)
+xmin, xmax, ymin, ymax = plt.axis()
+plt.text(xmin + 0.02*(xmax - xmin), ymax - 0.06*(xmax - xmin), "mu = {:.6f}".format(cfg.mu))
+
+fig.savefig('pseudo-potential.pdf', dpi=150)