Merge branch 'main' into jima

.gitlab-ci.yml

@@ -3,24 +3,39 @@ image: python:3.11-bullseye
 stages:
   - test-jobs
   - test-tasks
+  - build-page
 
-before_script:
+test-jobs:
+  stage: test-jobs
+  script:
     - pip install poetry
     - poetry --version
     - cd build/
     - poetry install --no-root
     - source $(poetry env info --path)/bin/activate
     - cd ..
-
-
-test-jobs:
-  stage: test-jobs
-  script:
     - export PYTHONPATH="$PYTHONPATH:$CI_PROJECT_DIR/jobs/src"
     - poetry run pytest jobs/tests/
 
 test-tasks:
   stage: test-tasks
   script:
+    - pip install poetry
+    - poetry --version
+    - cd build/
+    - poetry install --no-root
+    - source $(poetry env info --path)/bin/activate
+    - cd ..
     - export PYTHONPATH="$PYTHONPATH:$CI_PROJECT_DIR/tasks/src"
     - poetry run pytest tasks/tests/
+  
+pages:
+  stage: build-page
+  script:
+  - echo 'Nothing to do but build a website...'
+  artifacts:
+    paths:
+    - public
+    expire_in: 1 day
+  rules:
+    - if: $CI_COMMIT_BRANCH == "main"

README.md

@@ -51,4 +51,10 @@ For live experience (meaning you enter the container), run:
 
 If you change code in `src`, you need to rebuild the image with `bash build/build-image.sh`. The `src` folder can also be mounted in the image, but the assumption is that life is already difficult as it is...
 
+### interactive simulation
+
+- Run `python -m http.server`
+- Go to `http://localhost:8000/public`
+
+
 *Free and open-source software for all souls! Technical support is, unfortunately, only for group members.*
\ No newline at end of file

jobs/src/integrators.py

 """
-This module contains several integrators, which is at least 
-of second consistency order and approximately preserves 
+This module contains 2 integrators, which is at least 
+of second consistency order and approrimately preserves 
 the energy over long times.
 """
+from collections.abc import Callable
+
+import numpy as np
+
+
+def verlet(F: Callable, q0: np.ndarray, p0: np.ndarray, m: np.ndarray, dt: float):
+    """
+    Velocity-Verlet integrator for one time step
+    """
+    p = p0
+    q = q0
+
+    p = p + 1 / 2 * F(q) * dt
+    q = q + 1 / m * p * dt
+    p = p + 1 / 2 * F(q) * dt
+
+    return p, q

jobs/src/system.py

@@ -3,12 +3,55 @@ This module contains base class for a n-body gravitational system
 with 2 methods of simulations: direct and approximated force field.
 """
 
+from collections.abc import Callable
+from dataclasses import dataclass
 
-class gravitational_system():
-    pass
+import numpy as np
 
-    def direct_simulation():
-        pass
 
-    def force_field():
+@dataclass
+class GravitationalSystem:
+    r0: np.ndarray[float]  # shape = (N,d)
+    v0: np.ndarray[float]  # shape = (N,d)
+    m: np.ndarray[float]  # shape = N
+    t: np.ndarray[float]  # shape = n_time_steps
+    force: Callable
+    solver: Callable
+
+    def __post_init__(self):
+        "Checking dimensions of inputs"
+        assert self.r0.shape == self.v0.shape
+        assert self.m.shape[0] == self.r0.shape[0]
+
+        solvers = ['verlet']
+        assert self.solver.__name__ in solvers
+
+    def direct_simulation(self):
+        """
+        Using integrator to compute trajector in phase space
+        """
+        p = np.zeros((len(self.t), *self.v0.shape))
+        q = np.zeros((len(self.t), *self.r0.shape))
+        q[0] = self.r0
+        p[0] = self.m * self.v0
+
+        for i in range(1, len(self.t)):
+            dt = self.t[i] - self.t[i - 1]
+            p[i], q[i] = self.solver(self.force, q[i - 1], p[i - 1], self.m, dt)
+
+        return self.t, p, q
+
+    def force_field(self):
+        """
+        Using force field method to compute trajector in phase space
+        """
+        p = np.zeros((len(self.t), *self.v0.shape))
+        q = np.zeros((len(self.t), *self.r0.shape))
+        q[0] = self.r0
+        p[0] = self.m * self.v0
+        for i in range(1, len(self.t)):
+            # something with dt = self.t[i] - self.t[i-1]
+            # p = m*v
             pass
+
+        return self.t, p, q

jobs/tests/test_integrator.py

 """
 This unittest tests implementation of integrators
-for analytically well-understood configurations of
-gravitional n-body systems where n ∈ {2, 3}
+for the restricted three-body (sun-earth-moon) problem.
+The sun is fixed at the origin (center of mass). For 
+simplicity, the moon is assumed to be in the eliptic plane,
+so that vectors can be treated as two-dimensional. 
 """
+
+import logging
+import unittest
+
+import numpy as np
+
+from jobs.src.integrators import verlet
+from jobs.src.system import GravitationalSystem
+
+logger = logging.getLogger(__name__)
+logging.basicConfig(level=logging.INFO)
+
+# system settings
+R_SE = 1  # distance between Earth and the Sun, 1 astronomical unit (AU)
+R_L = 2.57e-3 * R_SE  # distance between Earth and the Moon
+
+M_S = 1  # solar mass
+M_E = 3.00e-6 * M_S
+M_L = 3.69e-8 * M_S
+
+T_E = 1  # earth yr
+T_L = 27.3 / 365.3 * T_E
+
+G = 4 * np.pi**2  # AU^3 / m / yr^2
+
+# simulation settings
+T = 0.3  # yr
+n_order = 6
+dt = 4**(-n_order)  # yr
+
+# vector of r0 and v0
+x0 = np.array([
+    R_SE,
+    0,
+    R_SE + R_L,
+    0,
+    0,
+    R_SE * 2 * np.pi / T_E,
+    0,
+    1 / M_E * M_E * R_SE * 2 * np.pi / T_E + 1 * R_L * 2 * np.pi / T_L,
+])
+
+
+# force computation
+def force(q: np.ndarray) -> np.ndarray:
+    r1 = q[0:2]
+    r2 = q[2:4]
+
+    return np.array([
+        -G * M_E * (M_S * r1 / np.linalg.norm(r1)**3 + M_L * (r1 - r2) / np.linalg.norm(r1 - r2)**3),
+        -G * M_L * (M_S * r2 / np.linalg.norm(r2)**3 + M_E * (r2 - r1) / np.linalg.norm(r2 - r1)**3),
+    ]).flatten()
+
+
+class IntegratorTest(unittest.TestCase):
+
+    def test_verlet(self):
+        """
+        Test functionalities of velocity-Verlet algorithm
+        """
+
+        system = GravitationalSystem(r0=x0[:4],
+                                     v0=x0[4:],
+                                     m=np.array([M_E, M_E, M_L, M_L]),
+                                     t=np.linspace(0, T, int(T // dt)),
+                                     force=force,
+                                     solver=verlet)
+
+        t, p, q = system.direct_simulation()
+
+        ## checking total energy conservation
+        H = np.linalg.norm(p[:,:2], axis=1)**2 / (2 * M_E) + np.linalg.norm(p[:,2:], axis=1)**2 / (2 * M_L) + \
+            -G * M_S * M_E / np.linalg.norm(q[:,:2], axis=1) - G * M_S * M_L / np.linalg.norm(q[:,2:], axis=1) + \
+            -G * M_E * M_L / np.linalg.norm(q[:,2:] - q[:,:2], axis=1)
+
+        self.assertTrue(np.greater(1e-10 + np.zeros(H.shape[0]), H - H[0]).all())
+
+        ## checking total linear momentum conservation
+        P = p[:, :2] + p[:, 2:]
+        self.assertTrue(np.greater(1e-10 + np.zeros(P[0].shape), P - P[0]).all())
+
+        ## checking total angular momentum conservation
+        L = np.cross(q[:, :2], p[:, :2]) + np.cross(q[:, 2:], p[:, 2:])
+        self.assertTrue(np.greater(1e-10 + np.zeros(L.shape[0]), L - L[0]).all())
+
+
+if __name__ == '__main__':
+    unittest.main()

poetry.lock

@@ -200,6 +200,17 @@ docs = ["furo", "jaraco.packaging (>=9.3)", "jaraco.tidelift (>=1.4)", "pygments
 testing = ["build[virtualenv]", "filelock (>=3.4.0)", "flake8-2020", "ini2toml[lite] (>=0.9)", "jaraco.develop (>=7.21)", "jaraco.envs (>=2.2)", "jaraco.path (>=3.2.0)", "pip (>=19.1)", "pytest (>=6)", "pytest-black (>=0.3.7)", "pytest-checkdocs (>=2.4)", "pytest-cov", "pytest-enabler (>=2.2)", "pytest-mypy (>=0.9.1)", "pytest-perf", "pytest-ruff", "pytest-timeout", "pytest-xdist", "tomli-w (>=1.0.0)", "virtualenv (>=13.0.0)", "wheel"]
 testing-integration = ["build[virtualenv] (>=1.0.3)", "filelock (>=3.4.0)", "jaraco.envs (>=2.2)", "jaraco.path (>=3.2.0)", "packaging (>=23.1)", "pytest", "pytest-enabler", "pytest-xdist", "tomli", "virtualenv (>=13.0.0)", "wheel"]
 
+[[package]]
+name = "toml"
+version = "0.10.2"
+description = "Python Library for Tom's Obvious, Minimal Language"
+optional = false
+python-versions = ">=2.6, !=3.0.*, !=3.1.*, !=3.2.*"
+files = [
+    {file = "toml-0.10.2-py2.py3-none-any.whl", hash = "sha256:806143ae5bfb6a3c6e736a764057db0e6a0e05e338b5630894a5f779cabb4f9b"},
+    {file = "toml-0.10.2.tar.gz", hash = "sha256:b3bda1d108d5dd99f4a20d24d9c348e91c4db7ab1b749200bded2f839ccbe68f"},
+]
+
 [[package]]
 name = "virtualenv"
 version = "20.25.0"
@@ -234,4 +245,4 @@ files = [
 [metadata]
 lock-version = "2.0"
 python-versions = "~3.11"
-content-hash = "79e1810a3524adfaaf1d9aa0a659f9dc3c9f35f46f7daa71963d06e011323054"
+content-hash = "8c0901f89e24d88032125c0e1c3297b9518cb1491f74361046bc8b771bae0fb9"

public/index.html

+<!-- HTML and Pyodide related sections are borrowed from Hendrik Weimer -->
+
+<!doctype html>
+<html>
+  <head>
+    <meta charset="utf-8">
+    <title>Three-body simulator</title>
+    <link rel="stylesheet" href="site.css">
+  </head>
+  <body>
+    <style>
+       body { font-size: 100%; font-family: monospace;}
+       .mtable table td tr { display: block; border-collapse: collapse;}
+       .mtable td { vertical-align: top; }
+       input { width: 180px; font-size: 100%; }
+       .error { color: red; font-size: 75%; }
+       .potential { width: 250px; }
+       .controls td { padding-left: 10px; padding-right: 10px;}
+       table.controls { border-radius: 8px; background: #c2c2c2;}
+       p.status { font-size: 100%; }
+       p.observables { font-size: 100%;}
+    </style>
+    <table class="mtable" style="width: 100%; max-width: 1024px;">
+      <tbody>
+	<tr>
+	  <td id="parent" style="width: 60%;">
+	    <p id="status" class="status">Loading Pyodide...</p>  
+	    <canvas id="canvas" style="width: 100%;"></canvas>
+	    <table class="controls">
+	      <tr>
+		<td>x<sub>0</sub></td>
+		<td>x<sub>1</sub></td>
+		<td>x<sub>2</sub></td>
+	      </tr>
+	      <tr>
+		<td><input id="x0" value="-1.5" ></td>
+		<td><input id="x1" value="0" ></td>
+		<td><input id="x2" value="2.0" ></td>
+	      </tr>
+	      <tr>
+		<td><p class="error" id="errorx0"></p></td>
+		<td><p class="error" id="errorx1"></p></td>
+		<td><p class="error" id="errorx2"></p></td>
+              </tr>
+	    </table>
+	  </td>
+	  <td style="width: 100%;">
+	    <table>
+	      <tr>
+		<td>
+                  <p class="observables" id="obs">
+		</td>
+	      </tr>
+	    </table>
+	  </td>
+	</tr>
+      </tbody>
+    </table>
+    <script type="text/javascript" src="https://cdn.jsdelivr.net/pyodide/v0.24.1/full/pyodide.js"></script>
+      <script type="text/javascript">
+      async function main() {
+        let pyodide = await loadPyodide();
+      await pyodide.loadPackage(["numpy", "scipy"]);
+      pyodide.runPython(await (await fetch("threebody.py")).text());
+      }
+      main();
+    </script>
+  </body>
+</html>
\ No newline at end of file

public/site.css

+@media (prefers-color-scheme: light) {
+    :root {
+	--background-color: #f2f2f2;
+	--softbox-color: #c2c2c2;
+	--foreground-color: #222222;
+	--link-color: blue;
+	--link-vis-color: rebeccapurple;
+	--emph-color: #116611;
+    }
+}
+
+@media (prefers-color-scheme: dark) {
+    :root {
+	--background-color: #222222;
+	--softbox-color: #333333;
+	--foreground-color: #f2f2f2;
+	--link-color: cyan;
+	--link-vis-color: rebeccapurple;
+	--emph-color: #99ff99;
+    }
+}
+
+body {
+    background: var(--background-color);
+    color: var(--foreground-color);
+    font-family: monospace;
+    max-width: 40em;
+}
+
+.box {
+    border-width: 1px;
+    border-style: solid;
+    border-color: var(--foreground-color);
+    border-radius: 8px;
+    margin: 1em;
+    padding: 1em;
+}
+
+.softbox {
+    background: var(--softbox-color);
+    border-radius: 8px;
+    margin: 1em;
+    padding: 1em;
+    width: fit-content;
+}
+
+.linkbox::before {
+    content: "=> ";
+}
+
+.shellbox > *::before {
+    content: "$ ";
+}
+
+.indent {
+    margin-left: 1em;
+}
+
+.list {
+    list-style-type: none;
+}
+
+.bold {
+    font-weight: bold;
+}
+
+a {
+    color: var(--link-color);
+}
+
+a:visited {
+    color: var(--link-vis-color);
+}
+
+.bmtitle {
+    font-weight: bold;
+}
+
+.bmtags {
+    font-style: italic;
+}
+
+.keepws {
+    white-space: pre;
+}
+
+em {
+    color: var(--emph-color);
+    font-style: normal;
+    font-weight: bold;
+}
+
+code {
+    color: var(--emph-color);
+}
+
+table {
+    border-spacing: 1em 0px;
+}
\ No newline at end of file

public/threebody.py

+"""
+HTML and Pyodide related sections are borrowed from Hendrik Weimer
+"""
+from js import window, document, console
+from pyodide.ffi import create_proxy
+
+import numpy as np
+from itertools import combinations
+from scipy.spatial.distance import pdist
+
+status = document.getElementById("status")
+
+proxy = dict()
+globals = dict()
+
+globals['N'] = 3
+globals['m'] = np.array([1, 2, 1])
+globals['x'] = np.zeros((globals['N'], 2))
+globals['p'] = np.zeros((globals['N'], 2))
+
+M_pixel = np.zeros((2, 2))
+globals['firstrun'] = True
+
+
+def E():
+    x = globals['x']
+    p = globals['p']
+    m = globals['m']
+
+    m_combinations = list(combinations(range(3), 2))
+    m_prod = np.array([m[i[0]] * m[i[1]] for i in m_combinations])
+    E_kin = np.sum(np.sum(p**2, axis=1) / (2 * m))
+    E_pot = np.sum(-m_prod / pdist(x))
+
+    return E_kin + E_pot
+
+
+def F(x):
+
+    x0 = x[0, :]
+    x1 = x[1, :]
+    x2 = x[2, :]
+
+    return np.array([
+        -globals['m'][0] * (globals['m'][1] * (x0 - x1) / np.linalg.norm(x0 - x1)**3 + globals['m'][2] *
+                            (x0 - x2) / np.linalg.norm(x0 - x2)**3),
+        -globals['m'][1] * (globals['m'][0] * (x1 - x0) / np.linalg.norm(x1 - x0)**3 + globals['m'][2] *
+                            (x1 - x2) / np.linalg.norm(x1 - x2)**3),
+        -globals['m'][2] * (globals['m'][0] * (x2 - x0) / np.linalg.norm(x2 - x0)**3 + globals['m'][1] *
+                            (x2 - x1) / np.linalg.norm(x2 - x1)**3),
+    ])
+
+
+def verlet(F, x0, p0, m, dt):
+    """
+    Verlet integrator for one time step
+    """
+    x = x0
+    p = p0
+
+    p = p + 1 / 2 * F(x) * dt
+    x = x + 1 / m[:, np.newaxis] * p * dt
+    p = p + 1 / 2 * F(x) * dt
+
+    return x, p
+
+
+def draw_circle(ctx, x, y, r=5):
+
+    pixelx = cwidth / 10 * (x + 5)
+    pixely = cheight / 10 * (y + 5)
+
+    ctx.beginPath()
+    ctx.arc(int(pixelx), int(pixely), r, 0, 2 * np.pi)
+    ctx.fillStyle = 'blue'
+    ctx.fill()
+    ctx.stroke()
+
+
+def replot_canvas():
+
+    dt = 1 / 25 / 10
+
+    canvas = document.getElementById("canvas")
+    ctx = canvas.getContext("2d")
+
+    ctx.clearRect(0, 0, cwidth, cheight)
+
+    x, p = verlet(F, globals['x'], globals['p'], globals['m'], dt)
+    globals['x'] = x
+    globals['p'] = p
+
+    for i in range(globals['N']):
+        draw_circle(ctx, x[i][0], x[i][1], 5 * globals['m'][i]**(1 / 3))
+
+    L0 = float(np.cross(x[0, :], p[0, :]))
+    L1 = float(np.cross(x[1, :], p[1, :]))
+    L2 = float(np.cross(x[2, :], p[2, :]))
+    L_tot = L0 + L1 + L2
+
+    P_tot = p[0, :] + p[1, :] + p[2, :]
+
+    COM = np.sum(globals['m']) * (globals['m'][0] * x[0, :] + globals['m'][1] * x[1, :] + globals['m'][2] * x[2, :])
+
+    obs = document.getElementById("angular_mom")
+    s = f"L0 = {L0:.6f}Ẑ<br>L1 = {L1:.6f}Ẑ<br>L2 = {L2:.6f}Ẑ<br>L_tot = {L_tot:.6f}Ẑ<br><br>E = {E():.6f}<br> COM = {COM[0]:.6f}X̂, {COM[1]:.6f}Ŷ<br> P_tot = {P_tot[0]:.6f}X̂, {P_tot[1]:.6f}Ŷ"
+    obs.innerHTML = s
+
+
+def push_queue(func, str):
+    status.innerHTML = str
+    proxy[func] = create_proxy(func)
+    window.setTimeout(proxy[func], 100)
+
+
+def set_main():
+
+    console.log("threebody: entering set_main")
+
+    global cwidth, cheight, pwidth
+
+    canvas = document.getElementById("canvas")
+
+    parentwidth = document.getElementById("parent").clientWidth
+    cwidth = parentwidth
+    cheight = int(cwidth * 0.75)
+
+    pwidth = cwidth
+
+    canvas.setAttribute("width", cwidth)
+    canvas.setAttribute("height", cheight)
+
+    console.log("threebody: exiting set_main")
+
+    push_queue(init_ode, "&nbsp;")
+
+
+def init_ode(foo=0):
+
+    console.log("threebody: entering init_ode")
+
+    if not globals['firstrun']:
+        window.clearInterval(globals['timer'])
+
+    exit_code = 0
+    for i in range(3):
+        error = document.getElementById("errorx" + str(i))
+        try:
+            globals['x'][i][0] = document.getElementById("x" + str(i)).value
+            error.innerHTML = ""
+        except ValueError:
+            error.innerHTML = "Input error"
+            status.innerHTML = ""
+            exit_code = 1
+
+    globals['p'][0][1] = 1
+    globals['p'][1][1] = 0
+    globals['p'][2][1] = -1
+
+    if exit_code == 0:
+        console.log("setting timer")
+        proxy[replot_canvas] = create_proxy(replot_canvas)
+
+        globals['timer'] = window.setInterval(proxy[replot_canvas], 40)
+
+    if globals['firstrun']:
+        for i in ["x0", "x1", "x2"]:
+            ii = document.getElementById(i)
+            ii.addEventListener("change", proxy[init_ode])
+
+        globals['firstrun'] = 0
+
+
+push_queue(set_main, "Setting up main window...")

pyproject.toml

@@ -7,6 +7,7 @@ readme = "README.md"
 
 [tool.poetry.dependencies]
 python = "~3.11"
+toml = "^0.10.2"
 
 [tool.poetry.group.dev.dependencies]
 ruff = "^0.0.267"