{
 "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 Roots\n",
    "using PyPlot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45fc58a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "    arenstorf!(dudt, u, 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",
    "    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": [
    "function deviation(u0, v0, tapprox)\n",
    "    initial_positions =  [u0,  0.0]\n",
    "    initial_velocities = [0.0,  v0]\n",
    "    tspan = (0., 1000.0)\n",
    "    p = 0.012277471\n",
    "    probl = ODEProblem(arenstorf!, \n",
    "        vcat(initial_positions, initial_velocities), tspan, p)\n",
    "    condition(u, t, integrator) = (t < tapprox) + u[2] + (u[4] > 0)\n",
    "    affect!(integrator) = terminate!(integrator)\n",
    "    cb = ContinuousCallback(condition, affect!)\n",
    "    sol = solve(probl, Tsit5(), callback=cb, \n",
    "        abstol=1e-15, reltol=1e-15)\n",
    "    ufin = sol(sol.t[end])[1]\n",
    "    return u0 - ufin, sol\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d5b7141",
   "metadata": {},
   "outputs": [],
   "source": [
    "tapprox = 5.0\n",
    "\n",
    "# u0 = 0.95\n",
    "# v0 = -2.0\n",
    "\n",
    "u0 = 1.15\n",
    "vinit = -1.0\n",
    "\n",
    "f(x) = deviation(u0, x, tapprox)[1];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e35adef",
   "metadata": {},
   "outputs": [],
   "source": [
    "v0 = find_zero(f, vinit, Order1())\n",
    "\n",
    "d, sol = deviation(u0, v0, tapprox);\n",
    "\n",
    "println(\"orbit deviation: $d\")\n",
    "println(\"orbit period: $(sol.t[end])\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d69f1e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = getindex.(sol.u, 1)\n",
    "y = getindex.(sol.u, 2);\n",
    "plot(x, y)\n",
    "axis(\"equal\")\n",
    "grid(true)"
   ]
  },
  {
   "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
}