diff --git a/day16.py b/day16.py
new file mode 100644
index 0000000..2db0cbd
--- /dev/null
+++ b/day16.py
@@ -0,0 +1,76 @@
+import sys
+
+def simulate(map, entry):
+    seen = set()
+    energised = set()
+    stack: list[tuple[tuple[int, int], tuple[int, int]]] = [entry]
+    while stack:
+        entry = stack.pop()
+        if entry in seen:
+            continue
+        seen.add(entry)
+
+        ((x, y), (dx, dy)) = entry
+        if not (0 <= y < len(map)) or not (0 <= x < len(map[y])):
+            continue
+
+        energised.add((x, y))
+
+        match map[y][x]:
+            case ".":
+                stack.append(((x + dx, y + dy), (dx, dy)))
+            case "-" if dy == 0:
+                stack.append(((x + dx, y + dy), (dx, dy)))
+            case "-":
+                stack.append(((x - 1, y), (-1, 0)))
+                stack.append(((x + 1, y), (1, 0)))
+            case "|" if dx == 0:
+                stack.append(((x + dx, y + dy), (dx, dy)))
+            case "|":
+                stack.append(((x, y - 1), (0, -1)))
+                stack.append(((x, y + 1), (0, 1)))
+            case "/" if (dx, dy) == (1, 0):
+                stack.append(((x, y - 1), (0, -1)))
+            case "/" if (dx, dy) == (-1, 0):
+                stack.append(((x, y + 1), (0, 1)))
+            case "/" if (dx, dy) == (0, 1):
+                stack.append(((x - 1, y), (-1, 0)))
+            case "/" if (dx, dy) == (0, -1):
+                stack.append(((x + 1, y), (1, 0)))
+            case "\\" if (dx, dy) == (1, 0):
+                stack.append(((x, y + 1), (0, 1)))
+            case "\\" if (dx, dy) == (-1, 0):
+                stack.append(((x, y - 1), (0, -1)))
+            case "\\" if (dx, dy) == (0, 1):
+                stack.append(((x + 1, y), (1, 0)))
+            case "\\" if (dx, dy) == (0, -1):
+                stack.append(((x - 1, y), (-1, 0)))
+                
+    return len(energised)
+
+def main():
+    with open(sys.argv[1]) as f:
+        map = [list(l) for l in f.read().strip().splitlines()]
+
+    part1 = simulate(map, ((0, 0), (1, 0)))
+
+    part2 = 0
+    for y in range(len(map)):
+        w = len(map[y])
+
+        xp = simulate(map, ((0, y), (1, 0)))
+        xm = simulate(map, ((w - 1, y), (-1, 0)))
+        part2 = max(part2, xp, xm)
+
+    for x in range(len(map)):
+        h = len(map)
+
+        yp = simulate(map, ((x, 0), (0, 1)))
+        ym = simulate(map, ((x, h - 1), (0, -1)))
+        part2 = max(part2, yp, ym)
+
+    print(part1, part2)
+
+
+if __name__ == "__main__":
+    main()