diff --git a/website/index.html b/website/index.html
new file mode 100644
index 0000000000000000000000000000000000000000..c3b7bc21533059b89bb3209d9eb897d743d2f8f7
--- /dev/null
+++ b/website/index.html
@@ -0,0 +1,69 @@
+<!-- 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="angular_mom">
+		</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
diff --git a/website/site.css b/website/site.css
new file mode 100644
index 0000000000000000000000000000000000000000..cb1bd57af932968dd621d039f985c752464d631a
--- /dev/null
+++ b/website/site.css
@@ -0,0 +1,99 @@
+@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
diff --git a/website/threebody.py b/website/threebody.py
new file mode 100644
index 0000000000000000000000000000000000000000..272dc272154585e83e98d26ab89b19874f5bd506
--- /dev/null
+++ b/website/threebody.py
@@ -0,0 +1,172 @@
+"""
+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
+
+    COM = np.sum(globals['m']) * (globals['m'][0] * x[0, :] + globals['m'][1] * x[1, :] + globals['m'][2] * x[2, :])
+
+    angular_mom = 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}Ŷ"
+    angular_mom.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...")