Stars Simulation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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"
    },
    {
     "data": {
      "application/javascript": [
       "if (typeof Jupyter !== \"undefined\") {require.undef(\"nbextensions/vpython_libraries/glow.min\");}else{element.textContent = ' ';}"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "if (typeof Jupyter !== \"undefined\") {require.undef(\"nbextensions/vpython_libraries/glowcomm\");}else{element.textContent = ' ';}"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "if (typeof Jupyter !== \"undefined\") {require.undef(\"nbextensions/vpython_libraries/jquery-ui.custom.min\");}else{element.textContent = ' ';}"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "if (typeof Jupyter !== \"undefined\") {require([\"nbextensions/vpython_libraries/glow.min\"], function(){console.log(\"GLOW LOADED\");});}else{element.textContent = ' ';}"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "if (typeof Jupyter !== \"undefined\") {require([\"nbextensions/vpython_libraries/glowcomm\"], function(){console.log(\"GLOWCOMM LOADED\");});}else{element.textContent = ' ';}"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "if (typeof Jupyter !== \"undefined\") {require([\"nbextensions/vpython_libraries/jquery-ui.custom.min\"], function(){console.log(\"JQUERY LOADED\");});}else{element.textContent = ' ';}"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from vpython import *\n",
    "\n",
    "# GlowScript 2.9 VPython\n",
    "\n",
    "# Stars interacting gravitationally\n",
    "\n",
    "# Bruce Sherwood\n",
    "\n",
    "scene.width = scene.height = 600\n",
    "\n",
    "# Display text below the 3D graphics:\n",
    "scene.title = \"Stars interacting gravitationally\"\n",
    "\n",
    "scene.caption = \"\"\"To rotate \"camera\", drag with right button or Ctrl-drag.\n",
    "To zoom, drag with middle button or Alt/Option depressed, or use scroll wheel.\n",
    "  On a two-button mouse, middle is left + right.\n",
    "To pan left/right and up/down, Shift-drag.\n",
    "Touch screen: pinch/extend to zoom, swipe or two-finger rotate.\n",
    "\n",
    "Replot the web page to restart with different (random) initial conditions.\"\"\"\n",
    "\n",
    "Nstars = 20  # change this to have more or fewer stars\n",
    "\n",
    "G = 6.7e-11 # Universal gravitational constant\n",
    "\n",
    "# Typical values\n",
    "Msun = 2E30\n",
    "Rsun = 2E9\n",
    "L = 4e10\n",
    "vsun = 0.8*sqrt(G*Msun/Rsun)\n",
    "\n",
    "scene.range = 2*L\n",
    "scene.forward = vec(-1,-1,-1)\n",
    "\n",
    "xaxis = curve(color=color.gray(0.5), radius=3e8)\n",
    "xaxis.append(vec(0,0,0))\n",
    "xaxis.append(vec(L,0,0))\n",
    "yaxis = curve(color=color.gray(0.5), radius=3e8)\n",
    "yaxis.append(vec(0,0,0))\n",
    "yaxis.append(vec(0,L,0))\n",
    "zaxis = curve(color=color.gray(0.5), radius=3e8)\n",
    "zaxis.append(vec(0,0,0))\n",
    "zaxis.append(vec(0,0,L))\n",
    "\n",
    "Stars = []\n",
    "star_colors = [color.red, color.green, color.blue,\n",
    "              color.yellow, color.cyan, color.magenta]\n",
    "\n",
    "psum = vec(0,0,0)\n",
    "for i in range(Nstars):\n",
    "    star = sphere(pos=L*vec.random(), make_trail=True, retain=150, trail_radius=2e8)\n",
    "    R = Rsun/2+Rsun*random()\n",
    "    star.radius = R\n",
    "    star.mass = Msun*(R/Rsun)**3\n",
    "    star.momentum = vec.random()*vsun*star.mass\n",
    "    star.color = star.trail_color = star_colors[i % 6]\n",
    "    Stars.append( star )\n",
    "    psum = psum + star.momentum\n",
    "\n",
    "#make total initial momentum equal zero\n",
    "for i in range(Nstars):\n",
    "    Stars[i].momentum = Stars[i].momentum - psum/Nstars\n",
    "\n",
    "dt = 1000\n",
    "hitlist = []\n",
    "\n",
    "def computeForces():\n",
    "    global hitlist, Stars\n",
    "    hitlist = []\n",
    "    N = len(Stars)\n",
    "    for i in range(N):\n",
    "        si = Stars[i]\n",
    "        if si is None: continue\n",
    "        F = vec(0,0,0)\n",
    "        pos1 = si.pos\n",
    "        m1 = si.mass\n",
    "        radius = si.radius\n",
    "        for j in range(N):\n",
    "            if i == j: continue\n",
    "            sj = Stars[j]\n",
    "            if sj is None: continue\n",
    "            r = sj.pos - pos1\n",
    "            rmag2 = mag2(r)\n",
    "            if rmag2 <= (radius+sj.radius)**2: hitlist.append([i,j])\n",
    "            F = F + (G*m1*sj.mass/(rmag2**1.5))*r\n",
    "        si.momentum = si.momentum + F*dt\n",
    "\n",
    "while True:\n",
    "    rate(100)\n",
    "    \n",
    "    # Compute all forces on all stars\n",
    "    computeForces()\n",
    "\n",
    "    # Having updated all momenta, now update all positions\n",
    "    for star in Stars:\n",
    "        if star is None: continue\n",
    "        star.pos = star.pos + star.momentum*(dt/star.mass)\n",
    "\n",
    "    # If any collisions took place, merge those stars\n",
    "    hit = len(hitlist)-1\n",
    "    while hit > 0:\n",
    "        s1 = Stars[hitlist[hit][0]]\n",
    "        s2 = Stars[hitlist[hit][1]]\n",
    "        if (s1 is None) or (s2 is None): continue\n",
    "    \n",
    "        mass = s1.mass + s2.mass\n",
    "        momentum = s1.momentum + s2.momentum\n",
    "        pos = (s1.mass*s1.pos + s2.mass*s2.pos) / mass\n",
    "        s1.color = s1.trail_color = (s1.mass*s1.color + s2.mass*s2.color) / mass\n",
    "        R = Rsun*(mass / Msun)**(1/3)\n",
    "        \n",
    "        s1.clear_trail()\n",
    "        s2.clear_trail()\n",
    "        s2.visible = False\n",
    "        \n",
    "        s1.mass = mass\n",
    "        s1.momentum = momentum\n",
    "        s1.pos = pos\n",
    "        s1.radius = R\n",
    "        Stars[hitlist[hit][1]] = None\n",
    "        hit -= 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}