{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "93d0a4f9",
   "metadata": {},
   "source": [
    "Arenstorf orbit - a stable periodic 3-body orbit where two bodies are\n",
    "orbiting in a plane, along with a third body of negligible mass (e.g. the Earth, the Moon, and a satellite)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3bf0f05c",
   "metadata": {},
   "outputs": [],
   "source": [
    "using OrdinaryDiffEq\n",
    "using Printf\n",
    "using PyPlot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45fc58a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "    arenstorf!(dudt, u, z, p, t)\n",
    "\n",
    "The system of differential equations for the restricted 3-body problem\n",
    "\"\"\"\n",
    "function arenstorf!(dudt, u, p, t)\n",
    "    x, y, dxdt, dydt = u\n",
    "    mu = p\n",
    "    mu1 = 1.0 - mu\n",
    "    d(r) = ((x + r)^2 + y^2)^(3/2)\n",
    "    d1 = d(mu)\n",
    "    d2 = d(-mu1)\n",
    "    \n",
    "    dudt[1] = u[3]\n",
    "    dudt[2] = u[4]\n",
    "    dudt[3] = x + 2*dydt  - mu1*(x + mu)/d1 - mu*(x - mu1)/d2\n",
    "    dudt[4] = y - 2*dxdt  - mu1*y/d1 - mu*y/d2\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "761f8fff",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initial conditions and ode parameters\n",
    "initial_positions =  [0.994,  0.0]\n",
    "initial_velocities = [0.0,   -2.00158510637908252240537862224]\n",
    "tspan = (0., 4*17.0652165601579625588917206249)\n",
    "p = 0.012277471;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d5b7141",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solve ODE\n",
    "prob = ODEProblem(arenstorf!, \n",
    "                  vcat(initial_positions, initial_velocities), \n",
    "                  tspan, p)\n",
    "sol = solve(prob, Tsit5(), abstol=1e-15, reltol=1e-15);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56e612f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualization parameters\n",
    "n = 1001             # Number of points\n",
    "t = range(tspan[1], tspan[2], length=n)\n",
    "tstep = t[2] - t[1]\n",
    "f = 0.025             # Fraction of displayed points\n",
    "pts_disp = round(Int, f*n);  # Number of displayed points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d69f1e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "sol_aren = sol.(t)\n",
    "x = getindex.(sol_aren, 1)\n",
    "y = getindex.(sol_aren, 2);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87cc4535",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the plotting limits\n",
    "x0 = -1.5\n",
    "x1 = -x0\n",
    "y0 = -1.5\n",
    "y1 = -y0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3aeb3c21",
   "metadata": {},
   "outputs": [],
   "source": [
    "function snapshot(nt, pts_disp, x0, x1, y0, y1, tspan, tstep)\n",
    "    nt_start = max(1, nt-pts_disp) \n",
    "    axis(\"square\")\n",
    "    xlim(x0, x1)\n",
    "    ylim(y0, y1)\n",
    "    titl = @sprintf(\"Arensdorf orbit for t = %5.2f .. %5.2f\",\n",
    "                     (nt_start - 1)*tstep + tspan[1],\n",
    "                     (nt - 1)*tstep + tspan[1])\n",
    "    title(titl, family=\"monospace\")\n",
    "    xlabel(L\"$x$\")\n",
    "    ylabel(L\"$y$\")\n",
    "    grid(true)\n",
    "    plot(x[nt_start:nt], y[nt_start:nt], color=\"red\", linewidth=0.5,\n",
    "         label=L\"$\\mu=0.0122775$\")\n",
    "    plot(x[nt], y[nt], color=\"red\", marker=\"o\", markersize=3.0)\n",
    "    legend()\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1dae58f",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = figure()\n",
    "for nt in 1:n\n",
    "    snapshot(nt, pts_disp, x0, x1, y0, y1, tspan, tstep)\n",
    "    display(fig)\n",
    "    sleep(0.01)\n",
    "    IJulia.clear_output(true) \n",
    "    clf()\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "539ac7b0",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.6.2",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}