C01_3Bodies.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b9d68d4c-ed1e-41b8-a4f2-bf2e1d6bf419",
   "metadata": {},
   "source": [
    "# Into our unknown 3 bodies jounery\n",
    "We use R to produce trajectories of three bodies according to Newton's gravitation. <br>\n",
    "We illustrates the three body problem where small changes to initial conditions cause large changes to later positions. <br>\n",
    "It can be **chaotic** 😂"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5194cbd-8483-4dc3-a35c-51acb9d1f2a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "#install.packages('plot3D')\n",
    "library(plot3D)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ca33b1b-32e7-4b89-8ebd-6071329229a6",
   "metadata": {},
   "source": [
    "### Set up our model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "80e70379-04d1-4cac-b39f-f2af95a23246",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "delta_t <- 0.001\n",
    "steps <- 100000\n",
    "\n",
    "# Masses of planets\n",
    "m_1 <- 10\n",
    "m_2 <- 20\n",
    "m_3 <- 30\n",
    "\n",
    "# Starting coordinates for planets\n",
    "p1_start <- c(-10, 10, -11)\n",
    "v1_start <- c(-3, 0, 0)\n",
    "\n",
    "p2_start <- c(0, 0, 0)\n",
    "v2_start <- c(0, 0, 0)\n",
    "\n",
    "p3_start <- c(10, 10, 12.00000)\n",
    "v3_start <- c(3, 0, 0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5500f6ef-cdb3-4146-b91b-1eebd3489fe1",
   "metadata": {},
   "source": [
    "### Set up the accelerations function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b499118f-97eb-424a-bc55-1fd6194e642a",
   "metadata": {},
   "outputs": [],
   "source": [
    "accelerations <- function(p1, p2, p3) {\n",
    "  planet_1_dv <- -9.8 * m_2 * (p1 - p2) / sum((p1 - p2)^2)^(3/2) - 9.8 * m_3 * (p1 - p3) / sum((p1 - p3)^2)^(3/2)\n",
    "  planet_2_dv <- -9.8 * m_3 * (p2 - p3) / sum((p2 - p3)^2)^(3/2) - 9.8 * m_1 * (p2 - p1) / sum((p2 - p1)^2)^(3/2)\n",
    "  planet_3_dv <- -9.8 * m_1 * (p3 - p1) / sum((p3 - p1)^2)^(3/2) - 9.8 * m_2 * (p3 - p2) / sum((p3 - p2)^2)^(3/2)\n",
    "  return(list(planet_1_dv, planet_2_dv, planet_3_dv))\n",
    "}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e49f7fd-047a-4a5d-9b63-8bdc88ff3e24",
   "metadata": {},
   "source": [
    "### Initialise solution arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4bcdf1ea-a10f-48d8-865f-aa9bfd6104e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "p1 <- matrix(0, nrow = steps, ncol = 3)\n",
    "v1 <- matrix(0, nrow = steps, ncol = 3)\n",
    "\n",
    "p2 <- matrix(0, nrow = steps, ncol = 3)\n",
    "v2 <- matrix(0, nrow = steps, ncol = 3)\n",
    "\n",
    "p3 <- matrix(0, nrow = steps, ncol = 3)\n",
    "v3 <- matrix(0, nrow = steps, ncol = 3)\n",
    "\n",
    "# Starting points\n",
    "p1[1,] <- p1_start\n",
    "v1[1,] <- v1_start\n",
    "\n",
    "p2[1,] <- p2_start\n",
    "v2[1,] <- v2_start\n",
    "\n",
    "p3[1,] <- p3_start\n",
    "v3[1,] <- v3_start\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40dbb10d-0630-4d8c-ac01-868d763d1958",
   "metadata": {},
   "source": [
    "### Compute evolution of the system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c94878e8-42ac-47fc-807b-ccf39c46947a",
   "metadata": {},
   "outputs": [],
   "source": [
    "for (i in 1:(steps - 1)) {\n",
    "  dv1 <- accelerations(p1[i,], p2[i,], p3[i,])[[1]]\n",
    "  dv2 <- accelerations(p1[i,], p2[i,], p3[i,])[[2]]\n",
    "  dv3 <- accelerations(p1[i,], p2[i,], p3[i,])[[3]]\n",
    "  \n",
    "  v1[i + 1,] <- v1[i,] + dv1 * delta_t\n",
    "  v2[i + 1,] <- v2[i,] + dv2 * delta_t\n",
    "  v3[i + 1,] <- v3[i,] + dv3 * delta_t\n",
    "  \n",
    "  p1[i + 1,] <- p1[i,] + v1[i,] * delta_t\n",
    "  p2[i + 1,] <- p2[i,] + v2[i,] * delta_t\n",
    "  p3[i + 1,] <- p3[i,] + v3[i,] * delta_t\n",
    "}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71aec04e-7f5d-4be9-8743-552005f86cd0",
   "metadata": {},
   "source": [
    "### Plot the trajectories (Here we go! Welcome to our geeky club)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70297c1b-5dca-4bad-92c8-3c03edee9d80",
   "metadata": {},
   "outputs": [],
   "source": [
    "scatter3D(p1[,1], p1[,2], p1[,3], colvar=NULL, col=\"dark green\", pch=19, cex=0.1, alpha=0.5, xlab=\"\", ylab=\"\", zlab=\"\")\n",
    "scatter3D(p2[,1], p2[,2], p2[,3], colvar=NULL, col=\"red\", pch=19, cex=0.1, alpha=0.5, xlab=\"\", ylab=\"\", zlab=\"\", add=TRUE)\n",
    "scatter3D(p3[,1], p3[,2], p3[,3], colvar=NULL, col=\"blue\", pch=19, cex=0.1, alpha=0.5, xlab=\"\", ylab=\"\", zlab=\"\", add=TRUE)\n",
    "title('3 Body Problem')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "RTrial",
   "language": "R",
   "name": "ir43"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.3.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}