N-body_Simulation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div id=\"glowscript\" class=\"glowscript\"></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "if (typeof Jupyter !== \"undefined\") { window.__context = { glowscript_container: $(\"#glowscript\").removeAttr(\"id\")};}else{ element.textContent = ' ';}"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import vpython as vp\n",
    "\n",
    "# Gravitational constant\n",
    "G = 1\n",
    "# Small number added to avoid singularities and to improve the numerical stability\n",
    "SOFTENING = 0.0001\n",
    "dt = 10e-4\n",
    "\n",
    "# Initial conditions\n",
    "position = np.array(\n",
    "    [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 3.0, 0.0], [0.0, -4.0, 0.0]]\n",
    ")\n",
    "\n",
    "velocity = np.array(\n",
    "    [[0.0, 0.0, 0.0], [0.0, 100.0, 0.0], [-57.0, 0.0, 0.0], [50.0, 0.0, 0.0]]\n",
    ")\n",
    "\n",
    "mass = np.array([[10000, 1, 5, 10]])\n",
    "\n",
    "# Create a canvas for the 3D scene\n",
    "scene = vp.canvas()\n",
    "\n",
    "\n",
    "def get_vp_vector(np_array):\n",
    "    return vp.vector(np_array[0], np_array[1], np_array[2])\n",
    "\n",
    "\n",
    "# Create bodies\n",
    "star = vp.sphere(\n",
    "    pos=get_vp_vector(position[0]), radius=0.4, color=vp.color.yellow, make_trail=True\n",
    ")\n",
    "planet1 = vp.sphere(\n",
    "    pos=get_vp_vector(position[1]), radius=0.25, color=vp.color.magenta, make_trail=True\n",
    ")\n",
    "planet2 = vp.sphere(\n",
    "    pos=get_vp_vector(position[2]), radius=0.15, color=vp.color.blue, make_trail=True\n",
    ")\n",
    "planet3 = vp.sphere(\n",
    "    pos=get_vp_vector(position[3]),\n",
    "    radius=0.20,\n",
    "    color=vp.color.green,\n",
    "    make_trail=True,\n",
    ")\n",
    "\n",
    "\n",
    "def get_accelaration(position, mass):\n",
    "    # Calculates the accelation\n",
    "    x = position[:, 0:1]\n",
    "    y = position[:, 1:2]\n",
    "    z = position[:, 2:3]\n",
    "\n",
    "    dx = x - x.T\n",
    "    dy = y - y.T\n",
    "    dz = z - z.T\n",
    "\n",
    "    r_norm = np.sqrt(dx**2 + dy**2 + dz**2 + SOFTENING**2)\n",
    "    r_norm_cube = r_norm**3\n",
    "\n",
    "    ax = (G * dx * mass.T) / r_norm_cube\n",
    "    ay = (G * dy * mass.T) / r_norm_cube\n",
    "    az = (G * dz * mass.T) / r_norm_cube\n",
    "    accelaration = np.vstack(\n",
    "        [np.sum(ax, axis=0), np.sum(ay, axis=0), np.sum(az, axis=0)]\n",
    "    ).T\n",
    "    return accelaration\n",
    "\n",
    "\n",
    "# Simulation's main loop\n",
    "for i in range(500):\n",
    "    vp.rate(30)\n",
    "\n",
    "    accelaration = get_accelaration(position, mass)\n",
    "\n",
    "    # Update the position and velocity\n",
    "    velocity += accelaration * dt\n",
    "    position += velocity * dt\n",
    "\n",
    "    star.pos = get_vp_vector(position[0])\n",
    "    planet1.pos = get_vp_vector(position[1])\n",
    "    planet2.pos = get_vp_vector(position[2])\n",
    "    planet3.pos = get_vp_vector(position[3])\n",
    "\n",
    "\n",
    "scene.capture(\"orbits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}