diff --git a/packages/engine/stdlib/src/py/dev_requirements.txt b/packages/engine/stdlib/src/py/dev_requirements.txt
index 6ba5375bfe6..17156bdee73 100644
--- a/packages/engine/stdlib/src/py/dev_requirements.txt
+++ b/packages/engine/stdlib/src/py/dev_requirements.txt
@@ -1,4 +1,4 @@
-pytest == "6.2.2"
-black == "20.8b1"
-pylint == "2.7.2"
+pytest == 6.2.2
+black == 20.8b1
+pylint == 2.7.2
 mypy==0.812
\ No newline at end of file
diff --git a/packages/engine/stdlib/src/py/hstd/agent.py b/packages/engine/stdlib/src/py/hstd/agent.py
index 42505476c6d..d77c8c7ed54 100644
--- a/packages/engine/stdlib/src/py/hstd/agent.py
+++ b/packages/engine/stdlib/src/py/hstd/agent.py
@@ -13,9 +13,16 @@ def generate_agent_id():
 @dataclass
 class AgentState:
     agent_id: str = field(default_factory=generate_agent_id)
+    agent_name: Optional[str] = None
     position: Optional[List[float]] = None
     direction: Optional[List[float]] = None
 
+    def __setitem__(self, key, value):
+        setattr(self, key, value)
+
+    def __getitem__(self, key):
+        return getattr(self, key)
+
 
 class AgentFieldError(Exception):
     def __init__(self, agent_id: str, field: str, msg: str = ""):
diff --git a/packages/engine/stdlib/src/py/hstd/context.py b/packages/engine/stdlib/src/py/hstd/context.py
new file mode 100644
index 00000000000..fcb5d154c19
--- /dev/null
+++ b/packages/engine/stdlib/src/py/hstd/context.py
@@ -0,0 +1,9 @@
+from dataclasses import dataclass, field
+from typing import Optional, List
+
+
+@dataclass
+class Topology:
+    x_bounds: List[float]
+    y_bounds: List[float]
+    z_bounds: Optional[List[float]] = field(default_factory=lambda: [0.0, 0.0])
diff --git a/packages/engine/stdlib/src/py/hstd/init.py b/packages/engine/stdlib/src/py/hstd/init.py
new file mode 100644
index 00000000000..55f73b0e631
--- /dev/null
+++ b/packages/engine/stdlib/src/py/hstd/init.py
@@ -0,0 +1,130 @@
+"""
+Initialization utility functions.
+"""
+import math
+import random
+from copy import deepcopy
+from typing import Dict, List, Union, Callable, Mapping
+
+from .agent import AgentState
+from .context import Topology
+
+# AgentTemplate can be an AgentState, or function which returns an AgentState
+AgentFunction = Callable[[], AgentState]
+AgentTemplate = Union[AgentState, AgentFunction]
+
+
+def create_agent(template: AgentTemplate) -> AgentState:
+    if callable(template):
+        return template()
+    else:
+        return deepcopy(template)
+
+
+def scatter(count: int, topology: Topology, template: AgentTemplate) -> List[AgentState]:
+    """
+    Generate `count` agents using the `template`, assigning them random positions within
+    the `topology` bounds.
+
+    Args:
+        count: the number of agents to generate.
+        topology: the `context.globals()["topology"]` value.
+        template: an agent definition, or a function which returns an agent definition.
+    """
+    x_bounds = topology.x_bounds
+    y_bounds = topology.y_bounds
+
+    width = x_bounds[1] - x_bounds[0]
+    height = y_bounds[1] - y_bounds[0]
+
+    def assign_random_position() -> AgentState:
+        x = random.uniform(0, width) + x_bounds[0]
+        y = random.uniform(0, height) + y_bounds[0]
+
+        agent = create_agent(template)
+        agent["position"] = [x, y]
+
+        return agent
+
+    agents = [assign_random_position() for i in range(count)]
+
+    return agents
+
+
+def stack(count: int, template: AgentTemplate) -> List[AgentState]:
+    """
+    Generate `count` agents using the `template`.
+
+    Args:
+        count: the number of agents to generate.
+        template: an agent definition, or a function which returns an agent definition.
+    """
+    agents = [create_agent(template) for i in range(count)]
+
+    return agents
+
+
+def grid(topology: Topology, template: AgentTemplate) -> List[AgentState]:
+    """
+    Generate agents on every integer location within the `topology` bounds.
+
+    Args:
+        topology: the `context.globals()["topology"]` value.
+        template: an agent definition, or a function which returns an agent definition.
+    """
+    x_bounds = topology.x_bounds
+    y_bounds = topology.y_bounds
+
+    width = x_bounds[1] - x_bounds[0]
+    height = y_bounds[1] - y_bounds[0]
+    count = width * height
+
+    def assign_grid_position(ind: int) -> AgentState:
+        x = (ind % width) + x_bounds[0]
+        y = math.floor(ind / width) + y_bounds[0]
+
+        agent = create_agent(template)
+        agent["position"] = [x, y]
+
+        return agent
+
+    agents = [assign_grid_position(i) for i in range(int(count))]
+
+    return agents
+
+
+def create_layout(
+    layout: List[List[str]], templates: Mapping[str, AgentState], offset: List[float] = [0, 0, 0]
+) -> List[AgentState]:
+    """
+    Generate agents with positions based on a `layout`, and definitions
+    based on the `templates`.
+
+    Args:
+        layout: the locations of agents, typically uploaded as a csv dataset
+        templates: the definitions for each type of agent refernced in the layout
+        offset: optional offset specifying the position of the bottom right corner of the `layout`
+    """
+
+    height = len(layout)
+    agents: Dict[str, List[AgentState]] = {}
+
+    for pos_y, row in enumerate(layout):
+        for pos_x, template_type in enumerate(row):
+            if template_type in templates:
+                if template_type not in agents:
+                    agents[template_type] = []
+
+                agent_name = (templates[template_type].agent_name or template_type) + str(
+                    len(agents[template_type])
+                )
+
+                agent = templates[template_type]
+                agent["agent_name"] = agent_name
+                agent["position"] = [pos_x + offset[0], height - pos_y + offset[1], offset[2]]
+
+                agents[template_type].append(agent)
+
+    agent_list = [agent for sublist in agents.values() for agent in sublist]
+
+    return agent_list
diff --git a/packages/engine/stdlib/src/py/hstd/neighbor.py b/packages/engine/stdlib/src/py/hstd/neighbor.py
new file mode 100644
index 00000000000..db64979e60b
--- /dev/null
+++ b/packages/engine/stdlib/src/py/hstd/neighbor.py
@@ -0,0 +1,131 @@
+"""
+Neighbor utility functions.
+"""
+from typing import List
+
+from .spatial import distance_between
+from .agent import AgentFieldError, AgentState
+
+
+def neighbors_on_position(agent: AgentState, neighbors: List[AgentState]) -> List[AgentState]:
+    """
+    Returns all `neighbors` whose position is identical to the `agent`.
+    """
+    if agent.position is None:
+        raise AgentFieldError(agent.agent_id, "position", "cannot be None")
+
+    return [n for n in neighbors if n.position == agent.position]
+
+
+def neighbors_in_radius(
+    agent: AgentState,
+    neighbors: List[AgentState],
+    max_radius: float = 1,
+    min_radius: float = 0,
+    distance_function: str = "euclidean",
+    z_axis: bool = False,
+) -> List[AgentState]:
+    """
+    Returns all neighbors within a certain vision radius of an agent.
+    Default is 2D (`z_axis` set to false). Set `z_axis` to true for 3D positions.
+
+    Args:
+        agent: central agent
+        neighbors: context.neighbors() array, or an array of agents
+        max_radius: minimum radius for valid neighbors
+        min_radius: maximum radius for valid neighbors
+        distance_function: type of distance function to use
+        z_axis: include z-axis in distance calculations
+    """
+
+    if agent.position is None:
+        raise AgentFieldError(agent.agent_id, "position", "cannot be None")
+
+    def in_radius(neighbor: AgentState) -> bool:
+        if neighbor.position is None:
+            return False
+
+        d = distance_between(neighbor, agent, distance_function, z_axis)
+        if d is None:
+            return False
+
+        return (d <= max_radius) and (d >= min_radius)
+
+    return [n for n in neighbors if in_radius(n)]
+
+
+def difference_vector(vec1: List[float], vec2: List[float]):
+    """
+    Calculate the difference vector `vec2` - `vec1`.
+    """
+    return [vec2[ind] - vec1[ind] for ind in range(len(vec1))]
+
+
+def in_front_planar(agent: AgentState, neighbor: AgentState) -> bool:
+    """
+    Return True if a neighbor is anywhere in front of the agent.
+    """
+
+    a_dir = agent["direction"]
+
+    [dx, dy, dz] = difference_vector(agent["position"], neighbor["position"])
+    D = a_dir[0] * dx + a_dir[1] * dy + a_dir[2] * dz
+
+    return D > 0
+
+
+def is_linear(agent: AgentState, neighbor: AgentState, front: bool) -> bool:
+    """
+    Check if a neighbor lies along the direction vector of the agent, and is
+    in front of or behind the agent, based on `front`.
+    """
+    [dx, dy, dz] = difference_vector(agent["position"], neighbor["position"])
+    [ax, ay, az] = agent["direction"]
+
+    cross_product = [dy * az - dz * ay, dx * az - dz * ax, dx * ay - dy * ax]
+
+    if cross_product != [0, 0, 0]:
+        return False
+
+    # check if same direction
+    same_dir = (ax * dx > 0) or (ay * dy > 0) or (az * dz > 0)
+
+    return same_dir is front
+
+
+def neighbors_in_front(
+    agent: AgentState, neighbors: List[AgentState], colinear: bool = False
+) -> List[AgentState]:
+    """
+    Return all `neighbors` in front of the `agent`. If `colinear` is True
+    check that the neighbor lies along the agent's direction vector.
+    """
+
+    if agent.position is None:
+        raise AgentFieldError(agent.agent_id, "position", "cannot be None")
+    if agent.direction is None:
+        raise AgentFieldError(agent.agent_id, "direction", "cannot be None")
+
+    if colinear:
+        return [n for n in neighbors if is_linear(agent, n, True)]
+    else:
+        return [n for n in neighbors if in_front_planar(agent, n)]
+
+
+def neighbors_behind(
+    agent: AgentState, neighbors: List[AgentState], colinear: bool = False
+) -> List[AgentState]:
+    """
+    Return all `neighbors` behind the `agent`. If `colinear` is True
+    check that the neighbor lies along the agent's direction vector.
+    """
+
+    if agent.position is None:
+        raise AgentFieldError(agent.agent_id, "position", "cannot be None")
+    if agent.direction is None:
+        raise AgentFieldError(agent.agent_id, "direction", "cannot be None")
+
+    if colinear:
+        return [n for n in neighbors if is_linear(agent, n, False)]
+    else:
+        return [n for n in neighbors if not in_front_planar(agent, n)]
diff --git a/packages/engine/stdlib/src/py/hstd/rand.py b/packages/engine/stdlib/src/py/hstd/rand.py
new file mode 100644
index 00000000000..ac388d18ac9
--- /dev/null
+++ b/packages/engine/stdlib/src/py/hstd/rand.py
@@ -0,0 +1,15 @@
+"""
+Random utility functions.
+"""
+
+import random as rand
+
+
+def set_seed(s: str):
+    """ Set the random seed for Python's random library """
+    rand.seed(s)
+
+
+def random():
+    """ Returns a random number between 0 and 1 """
+    return rand.random()
diff --git a/packages/engine/stdlib/src/py/hstd/spatial.py b/packages/engine/stdlib/src/py/hstd/spatial.py
index 7e035bb5903..dc40110c7ea 100644
--- a/packages/engine/stdlib/src/py/hstd/spatial.py
+++ b/packages/engine/stdlib/src/py/hstd/spatial.py
@@ -3,45 +3,40 @@
 """
 import math
 import random
-from dataclasses import dataclass
 from typing import Optional, List
 
 from .agent import AgentState, AgentFieldError
+from .context import Topology
 
 
-@dataclass
-class Topology:
-    x_bounds: List[float]
-    y_bounds: List[float]
-    z_bounds: Optional[List[float]]
-
-
-def manhattan_distance(p1: List[float], p2: List[float]) -> float:
+def manhattan_distance(p1: List[float], p2: List[float], z_axis: bool = True) -> float:
     dx = abs(p1[0] - p2[0])
     dy = abs(p1[1] - p2[1])
     dz = abs(p1[2] - p2[2])
-    return dx + dy + dz
+    return dx + dy + (dz if z_axis else 0)
 
 
-def euclidean_squared_distance(p1: List[float], p2: List[float]) -> float:
+def euclidean_squared_distance(p1: List[float], p2: List[float], z_axis: bool = True) -> float:
     dx = p1[0] - p2[0]
     dy = p1[1] - p2[1]
     dz = p1[2] - p2[2]
-    return dx * dx + dy * dy + dz * dz
+    return dx * dx + dy * dy + (dz * dz if z_axis else 0)
 
 
-def euclidean_distance(p1: List[float], p2: List[float]) -> float:
-    return math.sqrt(euclidean_squared_distance(p1, p2))
+def euclidean_distance(p1: List[float], p2: List[float], z_axis: bool = True) -> float:
+    return math.sqrt(euclidean_squared_distance(p1, p2, z_axis))
 
 
-def chebyshev_distance(p1: List[float], p2: List[float]) -> float:
+def chebyshev_distance(p1: List[float], p2: List[float], z_axis: bool = True) -> float:
     dx = abs(p1[0] - p2[0])
     dy = abs(p1[1] - p2[1])
     dz = abs(p1[2] - p2[2])
-    return max(dx, dy, dz)
+    return max(dx, dy, (dz if z_axis else 0))
 
 
-def distance_between(a: AgentState, b: AgentState, distance="euclidean") -> Optional[float]:
+def distance_between(
+    a: AgentState, b: AgentState, distance="euclidean", z_axis=True
+) -> Optional[float]:
     """
     Returns the specified distance between two agents. The parameter `distance` must be one
     of 'euclidean', 'euclidean_sq', 'manhattan' or 'chebyshev'.
@@ -52,13 +47,13 @@ def distance_between(a: AgentState, b: AgentState, distance="euclidean") -> Opti
         raise AgentFieldError(b.agent_id, "position", "cannot be None")
 
     if distance == "euclidean":
-        return euclidean_distance(a.position, b.position)
+        return euclidean_distance(a.position, b.position, z_axis)
     elif distance == "euclidean_sq":
-        return euclidean_squared_distance(a.position, b.position)
+        return euclidean_squared_distance(a.position, b.position, z_axis)
     elif distance == "manhattan":
-        return manhattan_distance(a.position, b.position)
+        return manhattan_distance(a.position, b.position, z_axis)
     elif distance == "chebyshev":
-        return chebyshev_distance(a.position, b.position)
+        return chebyshev_distance(a.position, b.position, z_axis)
 
     raise ValueError(
         "distance must be one of 'euclidean', 'euclidean_sq', 'manhattan' or 'chebyshev'"
diff --git a/packages/engine/stdlib/src/py/hstd/test_init.py b/packages/engine/stdlib/src/py/hstd/test_init.py
new file mode 100644
index 00000000000..888ca2f35f8
--- /dev/null
+++ b/packages/engine/stdlib/src/py/hstd/test_init.py
@@ -0,0 +1,65 @@
+from .agent import AgentState
+from .spatial import Topology
+from .init import scatter, grid, stack, create_layout
+
+init_topology = Topology([0, 2], [0, 2], [])
+
+agent = AgentState(agent_name="test")
+agent_function = lambda: AgentState(agent_name="test")
+
+num_agents = 4
+
+
+def test_scatter():
+    scatter_agents = scatter(num_agents, init_topology, agent)
+    scatter_agents_function = scatter(num_agents, init_topology, agent_function)
+
+    assert len(scatter_agents) == num_agents
+    assert len(scatter_agents_function) == num_agents
+
+    def subtest(a):
+        assert a["position"][0] >= init_topology.x_bounds[0]
+        assert a["position"][0] <= init_topology.x_bounds[1]
+        assert a["position"][1] >= init_topology.y_bounds[0]
+        assert a["position"][1] <= init_topology.y_bounds[1]
+
+        assert a["agent_name"] == "test"
+
+    [subtest(agent) for agent in scatter_agents]
+    [subtest(agent) for agent in scatter_agents_function]
+
+
+def test_stack():
+    stack_agents = stack(num_agents, agent)
+    stack_agents_function = stack(num_agents, agent_function)
+
+    assert len(stack_agents) == num_agents
+    assert len(stack_agents_function) == num_agents
+
+    def subtest(a):
+        assert a["agent_name"] == "test"
+
+    [subtest(agent) for agent in stack_agents]
+    [subtest(agent) for agent in stack_agents_function]
+
+
+def test_grid():
+    grid_agents = grid(init_topology, agent)
+    grid_agents_function = grid(init_topology, agent)
+
+    assert len(grid_agents) == num_agents
+    assert len(grid_agents_function) == num_agents
+
+    def subtest(a):
+        assert a["position"][0] >= init_topology.x_bounds[0]
+        assert a["position"][0] <= init_topology.x_bounds[1]
+        assert a["position"][1] >= init_topology.y_bounds[0]
+        assert a["position"][1] <= init_topology.y_bounds[1]
+
+        assert int(a["position"][0]) == a["position"][0]
+        assert int(a["position"][1]) == a["position"][1]
+
+        assert a["agent_name"] == "test"
+
+    [subtest(agent) for agent in grid_agents]
+    [subtest(agent) for agent in grid_agents_function]
diff --git a/packages/engine/stdlib/src/py/hstd/test_neighbor.py b/packages/engine/stdlib/src/py/hstd/test_neighbor.py
new file mode 100644
index 00000000000..a3d48197115
--- /dev/null
+++ b/packages/engine/stdlib/src/py/hstd/test_neighbor.py
@@ -0,0 +1,72 @@
+from .agent import AgentState
+from .neighbor import (
+    neighbors_on_position,
+    neighbors_in_radius,
+    neighbors_in_front,
+    neighbors_behind,
+)
+
+na = AgentState(position=[1, 1, 0], direction=[1, 0, 0])
+nb = AgentState(position=[1, 2, 0], direction=[1, 1, 0])
+nc = AgentState(position=[-1, 1, 0])
+nd = AgentState(position=[1, 1, 0])
+ne = AgentState(position=[2, 3, 0])
+nf = AgentState(position=[3, 2, 0])
+ng = AgentState(position=[6, 6, -1], direction=[1, 0, 0])
+nh = AgentState(position=[6, 9, 0])
+ni = AgentState(position=[4, 9, 0])
+nj = AgentState(position=[3, 2, 2])
+nk = AgentState(position=[3, 1, 0])
+nl = AgentState(position=[1, 0, 0], direction=[1, 1, 1])
+nm = AgentState(position=[0, 1, 0])
+nn = AgentState(position=[0, -1, -1])
+
+
+def same_position_test():
+    same_pos = neighbors_on_position(na, [nb, nc, nd, ne, nf])
+    assert same_pos == [nd]
+
+
+def test_max_radius():
+    in_radius = neighbors_in_radius(
+        ng, [na, nb, nc, nd, ne, nf, nh, ni], 3, distance_function="chebyshev"
+    )
+    assert in_radius == [nh, ni]
+
+
+def test_max_min_radius():
+    in_radius_1 = neighbors_in_radius(ng, [na, nb, nc, nd, nf, nh, ni, nj], 4, 3.5)
+    assert in_radius_1 == [ni]
+
+    in_radius_2 = neighbors_in_radius(
+        ng, [na, nb, nc, nd, nf, nh, ni, nj], 7, 3.5, "euclidean", True
+    )
+    assert in_radius_2 == [nb, nf, ni, nj]
+
+
+def test_front():
+    in_front_1 = neighbors_in_front(nb, [na, nc, ne, nf, ng, nh, nj])
+    assert in_front_1 == [ne, nf, ng, nh, nj]
+
+    in_front_2 = neighbors_in_front(na, [nb, nc, ne, nf, ng, nh, nj, nk, nl], True)
+    assert in_front_2 == [nk]
+
+    in_front_3 = neighbors_in_front(nl, [na, nb, nc, ne, nf, ng, nh, nj, nk], True)
+    assert in_front_3 == [nj]
+
+    in_front_4 = neighbors_in_front(nb, [na, nc, ne, nf, ng, nh, nj, nk], True)
+    assert in_front_4 == [ne]
+
+
+def test_behind():
+    behind_1 = neighbors_behind(nb, [na, ne, ng, nh, nj, nm])
+    assert behind_1 == [na, nm]
+
+    behind_2 = neighbors_behind(nb, [na, ne, ng, nh, nj, nm], True)
+    assert behind_2 == [nm]
+
+    behind_3 = neighbors_behind(na, [nb, nc, ne, ng, nh, nj, nm], True)
+    assert behind_3 == [nc, nm]
+
+    behind_4 = neighbors_behind(nl, [na, nb, nc, ne, ng, nh, nj, nm, nn], True)
+    assert behind_4 == [nn]
diff --git a/packages/engine/stdlib/src/py/hstd/test_rand.py b/packages/engine/stdlib/src/py/hstd/test_rand.py
new file mode 100644
index 00000000000..2e7dedc5e97
--- /dev/null
+++ b/packages/engine/stdlib/src/py/hstd/test_rand.py
@@ -0,0 +1,16 @@
+# type: ignore
+import pytest
+
+from .rand import set_seed, random
+
+
+def test_random():
+    n = random()
+    nn = random()
+    assert n != nn
+
+    set_seed("test")
+    n = random()
+    set_seed("test")
+    nn = random()
+    assert n == pytest.approx(nn)
diff --git a/packages/engine/stdlib/src/py/hstd/test_spatial.py b/packages/engine/stdlib/src/py/hstd/test_spatial.py
index 33faf481b1e..39dc85b0ebc 100644
--- a/packages/engine/stdlib/src/py/hstd/test_spatial.py
+++ b/packages/engine/stdlib/src/py/hstd/test_spatial.py
@@ -1,3 +1,4 @@
+# type: ignore
 import pytest
 
 from .agent import AgentState
@@ -36,8 +37,10 @@ def test_chebyshev_distance_between_tests():
 
 
 def test_normalize_direction():
-    assert normalize_vector(a.direction) == [0.7071067811865475, 0.7071067811865475]
-    assert normalize_vector(b.direction) == [0.31622776601683794, 0.9486832980505138]
+    if a.direction:
+        assert normalize_vector(a.direction) == [0.7071067811865475, 0.7071067811865475]
+    if b.direction:
+        assert normalize_vector(b.direction) == [0.31622776601683794, 0.9486832980505138]
 
 
 def test_random_position():
diff --git a/packages/engine/stdlib/src/ts/neighbor.spec.ts b/packages/engine/stdlib/src/ts/neighbor.spec.ts
index e122ce5c649..bb0d84053cb 100644
--- a/packages/engine/stdlib/src/ts/neighbor.spec.ts
+++ b/packages/engine/stdlib/src/ts/neighbor.spec.ts
@@ -25,7 +25,7 @@ test("find neighbors with same position", () => {
 });
 
 test("find neighbors within a radius of 3", () => {
-  expect(neighborsInRadius(ng, [na, nb, nc, nd, ne, nf, nh, ni], 3)).toEqual([
+  expect(neighborsInRadius(ng, [na, nb, nc, nd, ne, nf, nh, ni], 3, 0, "chebyshev")).toEqual([
     { position: [6, 9, 0] },
     { position: [4, 9, 0] },
   ]);
@@ -33,12 +33,12 @@ test("find neighbors within a radius of 3", () => {
 
 test("find neighbors within a max radius of 4 and min radius of 3", () => {
   expect(
-    neighborsInRadius(ng, [na, nb, nc, nd, nf, nh, ni, nj], 4, 3)
-  ).toEqual([{ position: [3, 2, 0] }, { position: [3, 2, 2] }]);
+    neighborsInRadius(ng, [na, nb, nc, nd, nf, nh, ni, nj], 4, 3.5)
+  ).toEqual([ni]);
 
   expect(
-    neighborsInRadius(ng, [na, nb, nc, nd, nf, nh, ni, nj], 4, 3, true)
-  ).toEqual([{ position: [3, 2, 2] }]);
+    neighborsInRadius(ng, [na, nb, nc, nd, nf, nh, ni, nj], 7, 3.5, "euclidean", true)
+  ).toEqual([nb, nf, ni, nj]);
 });
 
 test("find neighbors in front of agent tests", () => {
diff --git a/packages/engine/stdlib/src/ts/neighbor.ts b/packages/engine/stdlib/src/ts/neighbor.ts
index 793414e24d2..f4b0b5a8c4c 100644
--- a/packages/engine/stdlib/src/ts/neighbor.ts
+++ b/packages/engine/stdlib/src/ts/neighbor.ts
@@ -1,5 +1,6 @@
 /** Neighbor Functions */
 import { PotentialAgent } from "./agent";
+import { Distance, distanceBetween } from "./spatial";
 
 const posError = new Error("agent must have a position");
 const dirError = new Error("agent must have a direction");
@@ -38,39 +39,26 @@ export function neighborsOnPosition(
  * @param neighbors - context.neighbors() array, or an array of agents
  * @param max_radius - defaults to 1
  * @param min_radius - defaults to 0
+ * @param distanceFunction - defaults to "euclidean"
  * @param z_axis - defaults to false
  */
 export function neighborsInRadius(
   agent: PotentialAgent,
   neighbors: PotentialAgent[],
-  max_radius = 1,
-  min_radius = 0,
-  z_axis = false
+  max_radius: number = 1,
+  min_radius: number = 0,
+  distanceFunction: Distance = "euclidean",
+  z_axis: boolean = false
 ) {
+  const aPos = agent.position;
+  if (!aPos) { throw posError; }
+
   return neighbors.filter((neighbor) => {
-    const aPos = agent.position;
     const nPos = neighbor.position;
+    if (!nPos) { return false; }
 
-    if (!aPos || !nPos) {
-      throw posError;
-    }
-
-    const notZ: number = z_axis ? 0 : 1;
-
-    for (let i = 0; i < aPos.length - 1 * notZ; i++) {
-      const max = [aPos[i] + max_radius, aPos[i] - max_radius];
-      const min = [aPos[i] + min_radius, aPos[i] - min_radius];
-      if (
-        !(
-          (nPos[i] <= max[0] && nPos[i] >= min[0]) ||
-          (nPos[i] >= max[1] && nPos[i] <= min[1])
-        )
-      ) {
-        return false;
-      }
-    }
-
-    return true;
+    const d = distanceBetween(neighbor, agent, distanceFunction, z_axis);
+    return (d <= max_radius) && (d >= min_radius);
   });
 }
 
@@ -85,7 +73,7 @@ export function neighborsInRadius(
 export function neighborsInFront(
   agent: PotentialAgent,
   neighbors: PotentialAgent[],
-  colinear = false
+  colinear: boolean = false
 ) {
   return neighbors.filter((neighbor) => {
     const aPos = agent.position;
diff --git a/packages/engine/stdlib/src/ts/spatial.ts b/packages/engine/stdlib/src/ts/spatial.ts
index 88e45139ea8..0dc45b1d7a6 100644
--- a/packages/engine/stdlib/src/ts/spatial.ts
+++ b/packages/engine/stdlib/src/ts/spatial.ts
@@ -9,6 +9,36 @@ export type Distance =
   "euclidean_sq" |
   "chebyshev";
 
+const { abs, pow, max, sqrt } = Math;
+
+export function manhattan_distance(a_pos: number[], b_pos: number[], z_axis: boolean = true) {
+  const dx = abs(a_pos[0] - b_pos[0]);
+  const dy = abs(a_pos[1] - b_pos[1]);
+  const dz = abs(a_pos[2] - b_pos[2]);
+
+  return dx + dy + (z_axis ? dz : 0);
+}
+
+export function euclidean_squared_distance(a_pos: number[], b_pos: number[], z_axis: boolean = true) {
+  const dx = pow(a_pos[0] - b_pos[0], 2);
+  const dy = pow(a_pos[1] - b_pos[1], 2);
+  const dz = pow(a_pos[2] - b_pos[2], 2);
+
+  return dx + dy + (z_axis ? dz : 0);
+}
+
+export function euclidean_distance(a_pos: number[], b_pos: number[], z_axis: boolean = true) {
+  return sqrt(euclidean_squared_distance(a_pos, b_pos, z_axis));
+}
+
+export function chebyshev_distance(a_pos: number[], b_pos: number[], z_axis: boolean = true) {
+  const dx = abs(a_pos[0] - b_pos[0]);
+  const dy = abs(a_pos[1] - b_pos[1]);
+  const dz = abs(a_pos[2] - b_pos[2]);
+
+  return max(dx, dy, (z_axis ? dz : 0));
+}
+
 /**
  * Returns the specified distance between two agents.
  *  distance is one of the four distance functions supported by HASH,
@@ -20,49 +50,28 @@ export type Distance =
 export function distanceBetween(
   agentA: PotentialAgent,
   agentB: PotentialAgent,
-  distance: Distance = "euclidean"
+  distance: Distance = "euclidean",
+  zAxis: boolean = true
 ) {
-  type IdFuncs = {
-    // eslint-disable-next-line no-unused-vars
-    [index in Distance]: (a_pos: number[], b_pos: number[]) => number;
-  };
-
-  const { abs, pow, max, sqrt } = Math;
-
-  const dFuncs: IdFuncs = {
-    manhattan: (a_pos: number[], b_pos: number[]) =>
-      abs(a_pos[0] - b_pos[0]) +
-      abs(a_pos[1] - b_pos[1]) +
-      abs(a_pos[2] - b_pos[2]),
-    euclidean: (a_pos: number[], b_pos: number[]) =>
-      sqrt(
-        pow(a_pos[0] - b_pos[0], 2) +
-          pow(a_pos[1] - b_pos[1], 2) +
-          pow(a_pos[2] - b_pos[2], 2)
-      ),
-    euclidean_sq: (a_pos: number[], b_pos: number[]) =>
-      pow(a_pos[0] - b_pos[0], 2) +
-      pow(a_pos[1] - b_pos[1], 2) +
-      pow(a_pos[2] - b_pos[2], 2),
-    chebyshev: (a_pos: number[], b_pos: number[]) =>
-      max(
-        abs(a_pos[0] - b_pos[0]),
-        abs(a_pos[1] - b_pos[1]),
-        abs(a_pos[2] - b_pos[2])
-      ),
-  };
-
   const aPos = agentA.position;
   const bPos = agentB.position;
 
   if (!aPos || !bPos) {
     throw posError;
   }
-  if (!dFuncs[distance]) {
-    throw new Error("distance must be one of 'euclidean', 'manhattan', 'euclidean_sq' or 'chebyshev'");
-  }
 
-  return dFuncs[distance](aPos, bPos);
+  switch (distance) {
+    case "manhattan":
+      return manhattan_distance(aPos, bPos, zAxis);
+    case "euclidean":
+      return euclidean_distance(aPos, bPos, zAxis);
+    case "euclidean_sq":
+      return euclidean_squared_distance(aPos, bPos, zAxis);
+    case "chebyshev":
+      return chebyshev_distance(aPos, bPos, zAxis);
+    default:
+      throw new Error("distance must be one of 'euclidean', 'manhattan', 'euclidean_sq' or 'chebyshev'");
+  }
 }
 
 /**