Both phase networks in VPython

Local Jupyter is required for the native VPython display. The checked default cells inspect routes without opening a scene. Set open_viewer=True and run the remaining cells interactively; interrupt the final event loop when finished.

I display connected witness graphs for both phases and their independent wrapping loops. All positions are in voxel coordinates. Opacity, lighting and camera controls affect display only. The witness graphs are subsets of the measured full phase networks.

I develop this calculation in small steps. I run a cell, inspect its result and test it before adding the next part. Run the cells in order in Jupyter. A code cell defines a part of the calculation or carries it out; defining a function does not run it.

I permit free noncommercial teaching, learning and demonstrations with acknowledgment. Research and commercial use require my written permission. Earlier licensed releases retain their permissions. Read the code-use terms.

Copyright (c) 2026 Saswata Bhattacharyya.

Step 1: Prepare the Python environment

I install the packages listed in the next cell once in the same environment as the notebook. In Jupyter or Colab, remove the comment marker from the installation line below if a package is missing. NumPy stores arrays and evaluates numerical operations. Matplotlib displays plots. Pillow reads and writes images. This installation does not change the model parameters.

# %pip install numpy vpython "setuptools<82"
import numpy as np
assert np.isfinite(np.array([0.0, 1.0])).all()
print("NumPy is ready:", np.__version__)

NumPy is ready: 2.1.3

Step 2: Load data and decide when to open the scene

I download the public viewer data, then check the route definitions without opening VPython. I use local Jupyter for the native viewer; Colab does not provide the same interactive VPython environment. I leave the window disabled initially so Run All does not enter an endless display loop.

from pathlib import Path
from urllib.request import urlretrieve

work = Path("lesson-data")
work.mkdir(exist_ok=True)

def get_file(relative_path):
    # I download a public example only when its local copy is missing.
    destination = work / Path(relative_path).name
    if not destination.exists():
        urlretrieve("https://saswataiith.github.io/" + relative_path, destination)
    return destination

viewer_path = get_file("assets/examples/bicontinuity-3d/viewer.json")

viewer_data = __import__("json").loads(viewer_path.read_text())
assert viewer_data["n"] == 256
for phase in ("red", "blue"):
    for axis in "xyz":
        route = np.asarray(viewer_data["routes"][phase][axis])
        steps = (np.diff(route, axis=0) + 128) % 256 - 128
        assert np.all(np.abs(steps).sum(axis=1) == 1)
        assert np.array_equal(steps.sum(axis=0), np.eye(3, dtype=int)["xyz".index(axis)] * 256)
print("All six saved routes have face-sharing steps and the stated winding.")

All six saved routes have face-sharing steps and the stated winding.

Step 3: The optional interactive viewer

I first run the data checks with this setting left False. To use the viewer in local Jupyter, change it to True, run this cell, then run the cells below it in order. The last cell runs until I interrupt the kernel. I must restart the kernel before rebuilding a second scene.

open_viewer = False

Step 4: Prepare the next part of the calculation

I define d, n, scene, scene.forward, scene.lights, scene.ambient. These are the same operations as in the original source, written on separate lines so their order is clear. Inspect their values before running the next cell.

if open_viewer:
    import json, sys
    from pathlib import Path
    from vpython import canvas, vec, vertex, triangle, compound, curve, distant_light, menu, checkbox, slider, button, rate
    d = json.loads(Path(viewer_path).read_text())
    n = d['n']
    scene = canvas(title='Both percolating phase networks', width=820, height=540, background=vec(0.98, 0.98, 0.96), center=vec(n / 2, n / 2, n / 2), range=n * 0.75)
    scene.forward = vec(-1, -0.7, -1)
    scene.lights = []
    scene.ambient = vec(0.22, 0.22, 0.22)

Step 5: Repeat the operation and inspect each case

I repeat the operation over the specified cases. Read the loop variable first, then follow what changes on each pass. The assertions check the stated numerical or topological condition.

if open_viewer:
    for direction, intensity in [((-0.6, 0.8, -1), 0.45), ((1, -0.2, 0.4), 0.18), ((0.3, 0.7, 1), 0.12)]:
        distant_light(direction=vec(*direction), color=vec(intensity, intensity, intensity))

Step 6: Prepare the next part of the calculation

I define surfaces, lines, loops, colors, m. These are the same operations as in the original source, written on separate lines so their order is clear. Inspect their values before running the next cell.

if open_viewer:
    surfaces = []
    lines = {'red': [], 'blue': []}
    loops = {'red': [], 'blue': []}
    colors = {'red': vec(0.8, 0.08, 0.06), 'blue': vec(0.03, 0.28, 0.85)}
    m = d['meshes']['red']

Step 7: Repeat the operation and inspect each case

I repeat the operation over the specified cases. Read the loop variable first, then follow what changes on each pass. The assertions check the stated numerical or topological condition.

if open_viewer:
    for start in range(0, len(m['triangles']), 1500):
        vertices = {}

        def point(i):
            if i not in vertices:
                vertices[i] = vertex(pos=vec(*m['vertices'][i]), normal=vec(*m['normals'][i]), color=vec(0.5, 0.55, 0.58))
            return vertices[i]
        faces = [triangle(v0=point(t[0]), v1=point(t[1]), v2=point(t[2])) for t in m['triangles'][start:start + 1500]]
        obj = compound(faces)
        obj.opacity = 0.04
        obj.shininess = 0.25
        surfaces.append(obj)

Step 8: Draw periodic paths without a false shortcut

I split a displayed line at a periodic face and resume it at the paired opposite face. I retain turns in the voxel route. A line across the middle of the box would incorrectly suggest passage through another phase.

if open_viewer:
    def draw_path(path, color, radius, out):
        points = [[x + 0.5 for x in path[0]]]

        def flush():
            if len(points) < 2:
                return
            simple = [points[0]]
            for i in range(1, len(points) - 1):
                a, b, c = points[i - 1:i + 2]
                u = vec(*(b[j] - a[j] for j in range(3)))
                v = vec(*(c[j] - b[j] for j in range(3)))
                if u.cross(v).mag > 1e-08:
                    simple.append(b)
            simple.append(points[-1])
            out.append(curve(pos=[vec(*p) for p in simple], radius=radius, color=color, emissive=True))
        for a, b in zip(path, path[1:]):
            a = [x + 0.5 for x in a]
            b = [x + 0.5 for x in b]
            axis = next((j for j in range(3) if abs(b[j] - a[j]) > n / 2), None)
            if axis is not None:
                face = list(a)
                face[axis] = n if b[axis] < a[axis] else 0
                points.append(face)
                flush()
                partner = list(face)
                partner[axis] = 0 if face[axis] == n else n
                points = [partner, b]
            else:
                points.append(b)
        flush()

Step 9: Repeat the operation and inspect each case

I repeat the operation over the specified cases. Read the loop variable first, then follow what changes on each pass. The assertions check the stated numerical or topological condition.

if open_viewer:
    for phase in ('red', 'blue'):
        net = d['networks'][phase]
        adj = [[] for _ in net['nodes']]
        seen = set()
        for u, v in net['edges']:
            adj[u].append(v)
            adj[v].append(u)
        for root in range(len(adj)):
            if len(adj[root]) == 2:
                continue
            for first in adj[root]:
                edge = tuple(sorted((root, first)))
                if edge in seen:
                    continue
                prev = root
                current = first
                path = [net['nodes'][root]]
                seen.add(edge)
                while True:
                    path.append(net['nodes'][current])
                    if len(adj[current]) != 2:
                        break
                    nxt = next((v for v in adj[current] if v != prev))
                    edge = tuple(sorted((current, nxt)))
                    if edge in seen:
                        break
                    seen.add(edge)
                    prev, current = (current, nxt)
                draw_path(path, colors[phase], n / 650, lines[phase])
        assert len(seen) == len(net['edges'])
        for route in d['routes'][phase].values():
            draw_path(route, colors[phase], n / 330, loops[phase])

Step 10: Apply the display controls

I apply phase selection, interface visibility and loop highlighting to the scene. These are display operations on saved data.

if open_viewer:
    def update(_=None):
        for obj in surfaces:
            obj.visible = interface.checked
            obj.opacity = opacity.value
        for phase in ('red', 'blue'):
            visible = selection.selected in ('both', phase)
            for obj in lines[phase]:
                obj.visible = visible
            for obj in loops[phase]:
                obj.visible = visible and emphasized.checked

Step 11: Restore the initial camera

I reset the camera center, range and viewing direction. Right-drag rotates the view, and the wheel changes the zoom.

if open_viewer:
    def reset_camera(_):
        scene.center = vec(n / 2, n / 2, n / 2)
        scene.range = n * 0.75
        scene.forward = vec(-1, -0.7, -1)

Step 12: Prepare the next part of the calculation

I define selection, interface, emphasized, opacity. These are the same operations as in the original source, written on separate lines so their order is clear. Inspect their values before running the next cell.

if open_viewer:
    scene.append_to_caption('\nRed and blue connected witness networks. Thick lines emphasize independent wrapping loops.\n')
    selection = menu(choices=['both', 'red', 'blue'], selected='both', bind=update)
    interface = checkbox(text='Show c=0.5 interface', checked=True, bind=update)
    emphasized = checkbox(text='Emphasize wrapping loops', checked=True, bind=update)
    opacity = slider(min=0, max=1, step=0.02, value=0.04, bind=update)
    button(text='Reset camera', bind=reset_camera)
    scene.append_to_caption('\nRotate: right-drag / Ctrl-drag. Zoom: wheel. Periodic edges resume at the paired face.\nEach displayed graph connects 125 distributed sites and has winding rank 3.\nFull-volume connectivity: each entire phase is one connected component.\n')
    update()

Step 13: The viewer event loop

I keep the VPython clock running so mouse and button events can update the scene. Run this cell interactively and interrupt the kernel when finished. This is an event loop, not a simulation of phase separation.

if open_viewer:
    while True:
        rate(30)

Step 14: Review and extend the calculation

I first restart the kernel and run all the required cells. Then I change one input, predict its effect and compare the result with that prediction. A successful check establishes the tested property; it does not establish convergence or equilibrium. I retain the complete script for command-line use. The notebook embeds the calculation directly and does not need that script to run.

Read the notebook execution checks → · More Python, Julia and C lessons →