initial commit

Pipfile

@@ -0,0 +1,16 @@
+[[source]]
+url = "https://pypi.org/simple"
+verify_ssl = true
+name = "pypi"
+
+[packages]
+numpy = "*"
+ipython = "*"
+jupyter = "*"
+matplotlib = "*"
+scipy = "*"
+
+[dev-packages]
+
+[requires]
+python_version = "3.12"

integration/orbits.ipynb

@@ -0,0 +1,101 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "import scipy as sc\n",
+    "import scipy.integrate as integrate"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# constants\n",
+    "G = 1\n",
+    "M = 1"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [
+    {
+     "ename": "TypeError",
+     "evalue": "only length-1 arrays can be converted to Python scalars",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
+      "Cell \u001b[0;32mIn[12], line 26\u001b[0m\n\u001b[1;32m     24\u001b[0m x0 \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39marray([\u001b[38;5;241m1\u001b[39m, \u001b[38;5;241m0\u001b[39m, \u001b[38;5;241m0\u001b[39m, \u001b[38;5;241m0.1\u001b[39m])\n\u001b[1;32m     25\u001b[0m t \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mlinspace(\u001b[38;5;241m0\u001b[39m, \u001b[38;5;241m10\u001b[39m, \u001b[38;5;241m1000\u001b[39m)\n\u001b[0;32m---> 26\u001b[0m x \u001b[38;5;241m=\u001b[39m \u001b[43mf\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx0\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mt\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mintegrate\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mquad_vec\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     28\u001b[0m \u001b[38;5;66;03m# plot the solution\u001b[39;00m\n\u001b[1;32m     29\u001b[0m plt\u001b[38;5;241m.\u001b[39mplot(x[:, \u001b[38;5;241m0\u001b[39m], x[:, \u001b[38;5;241m1\u001b[39m])\n",
+      "Cell \u001b[0;32mIn[12], line 20\u001b[0m, in \u001b[0;36mf\u001b[0;34m(x, t, integrator)\u001b[0m\n\u001b[1;32m     19\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mf\u001b[39m(x: np\u001b[38;5;241m.\u001b[39mndarray, t: \u001b[38;5;28mfloat\u001b[39m, integrator: \u001b[38;5;28mcallable\u001b[39m) \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m>\u001b[39m np\u001b[38;5;241m.\u001b[39mndarray:\n\u001b[0;32m---> 20\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mintegrator\u001b[49m\u001b[43m(\u001b[49m\u001b[43mdf\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mt\u001b[49m\u001b[43m)\u001b[49m\n",
+      "File \u001b[0;32m~/.local/share/virtualenvs/projects-X-9bmgL6/lib/python3.12/site-packages/scipy/integrate/_quad_vec.py:225\u001b[0m, in \u001b[0;36mquad_vec\u001b[0;34m(f, a, b, epsabs, epsrel, norm, cache_size, limit, workers, points, quadrature, full_output, args)\u001b[0m\n\u001b[1;32m    105\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mquad_vec\u001b[39m(f, a, b, epsabs\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m1e-200\u001b[39m, epsrel\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m1e-8\u001b[39m, norm\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m2\u001b[39m\u001b[38;5;124m'\u001b[39m, cache_size\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m100e6\u001b[39m,\n\u001b[1;32m    106\u001b[0m              limit\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m10000\u001b[39m, workers\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m1\u001b[39m, points\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m, quadrature\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m, full_output\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m,\n\u001b[1;32m    107\u001b[0m              \u001b[38;5;241m*\u001b[39m, args\u001b[38;5;241m=\u001b[39m()):\n\u001b[1;32m    108\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124mr\u001b[39m\u001b[38;5;124;03m\"\"\"Adaptive integration of a vector-valued function.\u001b[39;00m\n\u001b[1;32m    109\u001b[0m \n\u001b[1;32m    110\u001b[0m \u001b[38;5;124;03m    Parameters\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    223\u001b[0m \n\u001b[1;32m    224\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m--> 225\u001b[0m     a \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mfloat\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43ma\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    226\u001b[0m     b \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mfloat\u001b[39m(b)\n\u001b[1;32m    228\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m args:\n",
+      "\u001b[0;31mTypeError\u001b[0m: only length-1 arrays can be converted to Python scalars"
+     ]
+    }
+   ],
+   "source": [
+    "### linearized ode\n",
+    "\n",
+    "def A(x: np.ndarray) -> np.ndarray:\n",
+    "    # assuming 2d data\n",
+    "    r = x[:x.size//2]\n",
+    "    v = x[x.size//2:]\n",
+    "    r_norm = np.linalg.norm(r)\n",
+    "    return np.array([\n",
+    "        [0, 0, 1, 0],\n",
+    "        [0, 0, 0, 1],\n",
+    "        [-G*M/r_norm**3 + 3*G*M*r[0]**2/r_norm**5, 3*G*M*r[0]*r[1]/r_norm**5, 0, 0],\n",
+    "        [3*G*M*r[0]*r[1]/r_norm**5, -G*M/r_norm**3 + 3*G*M*r[1]**2/r_norm**5, 0, 0]\n",
+    "    ])\n",
+    "\n",
+    "def df(x: np.ndarray, t) -> np.ndarray:\n",
+    "    return A(x) @ x\n",
+    "\n",
+    "\n",
+    "def f(x: np.ndarray, t: float, integrator: callable) -> np.ndarray:\n",
+    "    integrate.simpson()\n",
+    "    return integrator(df, x, t)\n",
+    "\n",
+    "\n",
+    "# solve the ode for some initial conditions\n",
+    "x0 = np.array([1, 0, 0, 0.1])\n",
+    "t = np.linspace(0, 10, 1000)\n",
+    "x = f(x0, t, integrate.quad_vec)\n",
+    "\n",
+    "# plot the solution\n",
+    "plt.plot(x[:, 0], x[:, 1])\n",
+    "plt.show()\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "projects-X-9bmgL6",
+   "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.12.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

nbody/notes.md

@@ -0,0 +1,23 @@
+###
+
+# Comparing with the analytical potential
+
+Plot a log-log plot of rho(r):
+- the discrete solution has poissonian error: 1/sqrt(N) where N is the number of particles in the bins
+
+=> ensure the discretization is consistent
+
+
+# Comparing forces
+- Use newtons law by considering only the inner shells:
+F_i = GM(R_i)/R_i^2
+
+
+=> this should all be done with the .ascii file (All particles, not yet meshed)
+
+
+
+
+## Task2
+a) start with only the monopole and compare with the direct summation
+only then build the tree code