Chapter 17 Complex Model Building Tutorial: Landward Movement, Schooling & Selective Tidal Stream Transport

⬇ Download this chapter

This tutorial builds on Chapter 15 by adding Selective Tidal Stream Transport (STST) to the landward movement and schooling framework. The environment now includes a procedural generated tidal velocity signal and a least-cost routing field, but still uses no real-world GIS data, empirical time series, or site-specific coordinates. The integration of STST with movement is the central modeling lesson: two functions that both govern movement must hand off control cleanly without conflicting.

17.1 Module Integration

STST wraps around active movement rather than running alongside it. At any given tick, the question is: can this agent overcome the current? If yes, the agent follows the least-cost path toward home-patch (migrate-landward). If no — meaning the current exceeds the agent’s effective swimming speed and that speed is already at its minimum — the agent enters STST, drifts passively with the flow, and conserves energy. When the agent’s swimming capacity recovers or the current weakens, it exits STST and resumes directed movement.

This means selective-tidal-stream-transport-landward is the gatekeeper for movement: it either calls migrate-landward internally or handles passive drift itself. migrate-landward is never called directly from the main loop in this configuration — STST does it.

Temporal structure. Each tick represents one generic time step. The tidal velocity cycle uses a sawtooth wave that oscillates between −0.5 and +0.5 m/s over a repeating 200-tick period. Positive velocity values represent landward flow, negative values represent seaward flow.

Spatial scale. The environment is a generic water channel. Agents are assigned a home-patch (upstream destination) and a migration-patch (downstream entry point) at setup. The least-cost routing field (cost-to-home) is reset to 1e6 each tick for all water patches, then re-calculated outward from home-patch before movement executes.

Coupling logic. For each time step, a velocity reporter writes a flow value uniformly to every water patch. calculate-difficulty-landward reads that velocity and writes difficulty-factor to swim landward. calculate-swimming-speed-landward reads difficulty-factor and energy and calculates individual swimming speed. STST then reads both velocity (patch) and speed (agent) to decide whether to drift or swim. If swimming, STST calls migrate-landward, which reads the current cost-to-home field and speed to step along the least-cost path. Schooling runs after movement, adjusting heading and speed for group cohesion before the tick closes.

17.2 Function Dependencies

Function Required Inputs Output Variables Dependent On
update-tidal-velocity ticks velocity (all water patches) patch environment
calculate-cost-to-home home-patch, patch-terrain cost-to-home patch velocity
calculate-difficulty-landward velocity, weight, max-weight difficulty-factor update-tidal-velocity
calculate-swimming-speed-landward energy, difficulty-factor, max-speed, min-speed, velocity, prev-speed speed calculate-difficulty-landward
selective-tidal-stream-transport-landward velocity, speed, min-speed, selective-tidal-transport?, STST-start-tick xcor, ycor, selective-tidal-transport?, tidal-transport-in-patch calculate-swimming-speed-landward
migrate-landward speed, home-patch, cost-to-home xcor, ycor, trail called from within STST
school schoolmates, nearest-neighbor, schooling coefficients heading, speed movement complete
calculate-swim-energy difficulty-factor, swim-efficiency, swim-base E-swim, energy executes after STST/movement

Key dependency: selective-tidal-stream-transport-landward must be called after calculate-swimming-speed-landward because it compares the freshly updated speed against the current patch velocity to decide whether STST triggers. Calling it before speed is updated would use a stale value and produce incorrect gating decisions.

STST state persistence. STST-start-tick records when the agent entered STST. A 1-tick cooldown prevents immediate re-triggering after exiting, which would cause oscillation. The selective-tidal-transport? boolean (0,1) must be initialized to false (0) at setup and never reset, which will be updated dynamically each time step.

17.3 Implementation in NetLogo

This example uses a generic grid environment with procedurally generated tidal velocity. There are no GIS data, CSV inputs, or calendar — just agents, patches, a sawtooth velocity signal, and the behavioral functions called in the sequence used by swim. A complete working implementation is available at Complex Model NetLogo Implementation.

17.3.1 Variable declarations

17.3.2 Setup

17.3.3 The go procedure — time step sequence

The structure directly mirrors swim: cost field → difficulty → speed → STST (which calls migrate-landward with schooling nested inside) → swim energy.

17.3.4 What to observe

  • STST vs. active swimming: Color agents with ifelse selective-tidal-transport? [set color orange][set color cyan]. Pulses of orange should appear during the seaward half of each 200-tick cycle.
  • Schooling during active movement only: Plot mean [distance nearest-neighbor] of turtles with [any? schoolmates]. It should stabilize near minimum-separation (0.3) during active swimming ticks. Agents in STST will have empty schoolmates agentsets — this is correct behavior.
  • Routing: Agents should follow the low-cost corridor upstream. Watch trails diverge around high-velocity patches.
  • Net upstream progress: Plot mean [ycor] of turtles. Despite tidal oscillation, mean position should trend toward home-patch, producing a staircase pattern.

17.4 Implementation in R

This example builds directly on the NetLogo submodel files developed throughout the book — Landward-Migration.nls, Selective-Tidal-Stream-Transport.nls, and Schooling.nls and its component files. The tidal velocity signal and world setup are generic stand-ins for empirical inputs; everything else follows the variable conventions and function call sequence established in those files. A complete working implementation is available at Complex Model R Implementation.

17.4.1 Agent and patch data structures

library(dplyr)
library(ggplot2)

# 1-D channel: patch 1 = migration-patch (downstream), patch 60 = home-patch
n_patches <- 60
n_ticks   <- 600   # 3 full 200-tick velocity cycles

patches <- tibble(
  patch_id     = 1:n_patches,
  terrain      = "water",
  velocity     = 0,
  depth        = runif(n_patches, 0, 5),
  cost_to_home = 1e6
)

n_agents <- 20

agents <- tibble(
  id                   = 1:n_agents,
  x                    = 1,
  home_patch           = n_patches,
  start_migration      = TRUE,
  landward_migration   = TRUE,
  at_destination       = FALSE,
  energy               = 100,
  speed                = 0.1,
  prev_speed           = 0.1,
  min_speed            = 0.1,
  max_speed            = sample(100:200, n_agents, replace = TRUE),
  swim_efficiency      = 1,
  swim_base            = 0.02,
  difficulty_factor    = 1,
  E_swim               = 0,
  weight               = sample(50:100, n_agents, replace = TRUE),
  in_stst              = FALSE,
  stst_start_tick      = 0,
  heading              = 90,
  minimum_separation   = 0.3,
  cohere_coefficient   = 0.50,
  align_coefficient    = 0.50,
  separate_coefficient = 0.50
)

max_seaward_velocity  <-  1.5
max_landward_velocity <- -1.5

17.4.2 Submodel functions

# ── Tidal velocity — 200-tick sawtooth, ±0.5 m/s ─────────────────────────────
tidal_velocity <- function(tick) {
  cycle_tick <- tick %% 200
  if (cycle_tick < 100)
    0.5 - cycle_tick * (1.0 / 100)
  else
    -0.5 + (cycle_tick - 100) * (1.0 / 100)
}

# ── Least-cost routing field (mirrors calculate-cost-to-home) ─────────────────
# Cumulative difficulty-weighted cost from home outward.
# compute-travel-cost base = abs(velocity) + 0.5 * max(0, depth slope).
# Here velocity is uniform per tick so the slope term drives variation.
calc_cost_to_home <- function(n_patches, difficulty, depths, velocity) {
  cost <- rep(1e6, n_patches)
  cost[n_patches] <- 0
  for (i in (n_patches - 1):1) {
    v_resistance <- abs(velocity)
    slope        <- max(0, depths[i + 1] - depths[i])
    base_cost    <- v_resistance + 0.5 * slope
    cost[i]      <- cost[i + 1] + base_cost * difficulty
  }
  cost
}

# ── Hydrodynamic difficulty (mirrors calculate-difficulty-landward) ───────────
calc_difficulty <- function(velocity, weight, max_weight,
                            v_max = max_seaward_velocity,
                            v_min = max_landward_velocity,
                            k = 0.75) {
  norm_v <- (velocity - v_min) / (v_max - v_min)
  df_raw <- (norm_v / (weight / max_weight))^k
  pmin(pmax(1 + 9 * df_raw, 1), 10)
}

# ── Swimming speed (mirrors calculate-swimming-speed-landward) ────────────────
calc_swim_speed <- function(energy, difficulty, max_speed, min_speed,
                            velocity, prev_speed, swim_efficiency,
                            k = -0.75, max_change = 0.5) {
  energy_factor   <- energy / 100
  velocity_impact <- k * velocity * 300
  desired <- min(max_speed, max(min_speed,
                (max_speed * energy_factor / difficulty) + velocity_impact))
  max_change_eff <- max_change * swim_efficiency
  delta  <- desired - prev_speed
  speed  <- prev_speed + sign(delta) * min(abs(delta), max_change_eff)
  max(speed, min_speed)
}

# ── Swimming energy cost (mirrors calculate-swim-energy) ─────────────────────
calc_swim_energy <- function(difficulty, swim_efficiency, swim_base,
                             beta = 0.75) {
  energy_mult <- difficulty * (1 / swim_efficiency)
  swim_base * (energy_mult^beta)
}

# ── Schooling — mirrors the BOIDS block inside migrate-landward ───────────────
# Called from inside migrate_landward_lc, not from the main loop.
school_agent <- function(agent_id, agents) {
  focal <- agents[agents$id == agent_id, ]
  mates <- agents[agents$id != agent_id &
                    agents$start_migration &
                    agents$landward_migration == focal$landward_migration, ]

  if (nrow(mates) == 0)
    return(list(heading = focal$heading, speed = focal$speed))

  dists   <- abs(mates$x - focal$x)
  nn_dist <- min(dists)
  nn_x    <- mates$x[which.min(dists)]

  new_heading <- focal$heading
  new_speed   <- focal$speed

  # separate / cohere+align
  if (nn_dist < focal$minimum_separation) {
    if (nn_x > focal$x) {
      if (focal$speed > focal$min_speed)
        new_speed <- focal$speed - 0.1                          # separate: slow
    } else {
      new_heading <- focal$heading +                             # separate: turn
        (focal$separate_coefficient / 100) * sign(focal$x - nn_x) * 10
    }
  } else {
    centroid_x  <- mean(mates$x)
    new_heading <- focal$heading +                               # cohere
      (focal$cohere_coefficient / 100) * sign(centroid_x - focal$x) * 10
    mean_heading <- mean(mates$heading)
    new_heading  <- new_heading +                                # align
      (focal$align_coefficient / 100) * (mean_heading - new_heading)
  }

  # adjust-speed
  if (max(mates$speed) > focal$speed + 0.1) {
    new_speed <- focal$speed + focal$speed * 0.1
  } else if (focal$speed > focal$min_speed) {
    new_speed <- focal$speed - focal$speed * 0.1
  }
  new_speed <- max(new_speed, focal$min_speed)

  list(heading = new_heading %% 360, speed = new_speed)
}

# ── Landward movement with schooling nested inside ───────────────────────────
# Mirrors migrate-landward: BOIDS first, then move along least-cost path.
migrate_landward_lc <- function(agent_id, agents) {
  school_result <- school_agent(agent_id, agents)
  focal  <- agents[agents$id == agent_id, ]
  steps  <- max(1, ceiling(school_result$speed / 3))
  new_x  <- min(focal$x + steps, focal$home_patch)
  list(x = new_x, heading = school_result$heading, speed = school_result$speed)
}

# ── STST gatekeeper (mirrors selective-tidal-stream-transport-landward) ────────
stst_landward <- function(agent_id, agents, velocity, current_tick) {
  focal  <- agents[agents$id == agent_id, ]
  S_swim <- focal$speed / 300

  # Case 1: cooldown active — continue drifting
  if (focal$in_stst && (current_tick - focal$stst_start_tick < 1)) {
    drift_dist <- abs(velocity) * (1 - focal$swim_efficiency)
    new_x      <- min(focal$x + ceiling(drift_dist), focal$home_patch)
    return(list(x = new_x, heading = focal$heading, speed = focal$speed,
                in_stst = TRUE, stst_start_tick = focal$stst_start_tick,
                energy = focal$energy))
  }

  # Case 2: flow exceeds swim capacity — enter STST
  if (abs(velocity) > S_swim && S_swim <= focal$min_speed) {
    drift_dist <- abs(velocity) * (1 - focal$swim_efficiency)
    new_x      <- min(focal$x + ceiling(drift_dist), focal$home_patch)
    return(list(x = new_x, heading = focal$heading, speed = focal$speed,
                in_stst = TRUE, stst_start_tick = current_tick,
                energy = focal$energy))
  }

  # Case 3: swim capacity sufficient — exit STST, run migrate (with schooling)
  move   <- migrate_landward_lc(agent_id, agents)
  E_swim <- calc_swim_energy(focal$difficulty_factor,
                             focal$swim_efficiency, focal$swim_base)
  list(x = move$x, heading = move$heading, speed = move$speed,
       in_stst = FALSE, stst_start_tick = focal$stst_start_tick,
       energy = max(focal$energy - E_swim, 0))
}

17.4.3 Simulation loop

log_ticks  <- vector("list", n_ticks)
max_weight <- max(agents$weight)

for (tick in 1:n_ticks) {

  # 1. Tidal velocity
  v                <- tidal_velocity(tick)
  patches$velocity <- v

  # 2. Rebuild least-cost field
  mean_diff            <- calc_difficulty(v, mean(agents$weight), max_weight)
  patches$cost_to_home <- calc_cost_to_home(n_patches, mean_diff,
                                            patches$depth, v)

  # 3. Per-agent difficulty and swimming speed (vectorised)
  agents <- agents |>
    mutate(
      difficulty_factor = calc_difficulty(v, weight, max_weight),
      speed = calc_swim_speed(energy, difficulty_factor, max_speed,
                              min_speed, v, prev_speed, swim_efficiency),
      prev_speed = speed
    )

  # 4. STST gatekeeper — calls migrate_landward_lc (with schooling) for active agents
  for (i in seq_len(nrow(agents))) {
    if (!agents$start_migration[i] || !agents$landward_migration[i] ||
        agents$at_destination[i]) next

    result <- stst_landward(
      agent_id     = agents$id[i],
      agents       = agents,
      velocity     = v,
      current_tick = tick
    )

    agents$x[i]               <- result$x
    agents$heading[i]         <- result$heading
    agents$speed[i]           <- result$speed
    agents$in_stst[i]         <- result$in_stst
    agents$stst_start_tick[i] <- result$stst_start_tick
    agents$energy[i]          <- result$energy
  }

  # 5. Swimming energy cost (baseline cost for all ticks including drift)
  agents <- agents |>
    mutate(
      E_swim = calc_swim_energy(difficulty_factor, swim_efficiency, swim_base),
      energy = pmax(energy - E_swim, 0)
    )

  # 6. Arrival check
  agents <- agents |>
    mutate(
      at_destination     = x >= home_patch,
      landward_migration = !at_destination
    )

  log_ticks[[tick]] <- tibble(
    tick        = tick,
    velocity    = v,
    mean_x      = mean(agents$x),
    mean_energy = mean(agents$energy),
    n_stst      = sum(agents$in_stst),
    n_arrived   = sum(agents$at_destination)
  )
}

results <- bind_rows(log_ticks)

17.4.4 Visualizing outputs

# Migration progress with tidal velocity overlay
ggplot(results, aes(x = tick)) +
  geom_line(aes(y = mean_x, color = "Mean position"), linewidth = 1) +
  geom_line(aes(y = velocity * 20 + 30, color = "Tidal velocity (scaled)"),
            linetype = "dashed", linewidth = 0.6) +
  scale_color_manual(values = c("Mean position"           = "#2c7bb6",
                                "Tidal velocity (scaled)" = "#d7191c")) +
  labs(title = "Migration progress with tidal signal",
       x = "Tick", y = "Mean upstream position", color = NULL) +
  theme_minimal()

# Agents in STST per tick
ggplot(results, aes(x = tick, y = n_stst)) +
  geom_area(fill = "#fdae61", alpha = 0.6) +
  geom_line(color = "#f46d43", linewidth = 0.8) +
  labs(title = "Agents in STST per tick",
       x = "Tick", y = "Count in passive drift") +
  theme_minimal()

# Energy over time
ggplot(results, aes(x = tick, y = mean_energy)) +
  geom_line(color = "#1a9641", linewidth = 1) +
  labs(title = "Mean agent energy over time",
       x = "Tick", y = "Energy") +
  theme_minimal()

What to look for:

  • Position plot: The 200-tick sawtooth produces a regular alternation between positive and negative flow. During the seaward-flow half, mean position plateaus as agents drift in STST. Progress resumes during the landward-flow half, producing a staircase pattern.
  • STST count: Pulses should align with the negative-velocity half of each 200-tick cycle.
  • Energy: Depletion slows during STST pulses. Agents drifting passively pay only the baseline swim_base cost rather than the full difficulty-scaled cost of active swimming.

17.5 Implementation in Python

The Python version follows the same structure as the R implementation — a pandas DataFrame of agents updated in a loop, with schooling nested inside migrate_landward_lc, which is itself called from inside stst_landward. The function nesting mirrors the NetLogo structure exactly. A complete working implementation is available at Complex Model Python Implementation.

17.5.1 Agent and patch data structures

import numpy as np
import pandas as pd

rng = np.random.default_rng(seed=42)

# 1-D channel: index 0 = migration-patch (downstream),
#              index 59 = home-patch (upstream)
n_patches = 60
n_ticks   = 600   # 3 full 200-tick velocity cycles

patches = pd.DataFrame({
    "patch_id"    : np.arange(n_patches),
    "terrain"     : "water",
    "velocity"    : 0.0,
    "depth"       : rng.uniform(0, 5, n_patches),
    "cost_to_home": 1e6,
})

n_agents = 20

agents = pd.DataFrame({
    "id"                  : np.arange(n_agents),
    "x"                   : 0,                                        # start at patch index 0 (downstream)
    "home_patch"          : n_patches - 1,                            # upstream destination (patch index 59)
    "start_migration"     : True,
    "landward_migration"  : True,
    "at_destination"      : False,
    "energy"              : 100.0,
    "speed"               : 0.1,
    "prev_speed"          : 0.1,
    "min_speed"           : 0.1,
    "max_speed"           : rng.integers(100, 200, n_agents).astype(float),
    "swim_efficiency"     : 1.0,
    "swim_base"           : 0.02,
    "difficulty_factor"   : 1.0,
    "E_swim"              : 0.0,
    "weight"              : rng.integers(50, 100, n_agents).astype(float),
    "in_stst"             : False,
    "stst_start_tick"     : 0,
    "heading"             : 90.0,
    "minimum_separation"  : 0.3,
    "cohere_coefficient"  : 0.50,
    "align_coefficient"   : 0.50,
    "separate_coefficient": 0.50,
})

MAX_SEAWARD_VELOCITY  =  1.5
MAX_LANDWARD_VELOCITY = -1.5

17.5.2 Submodel functions

# ── Tidal velocity — 200-tick sawtooth, ±0.5 m/s ─────────────────────────────
# Positive = landward flow, negative = seaward flow.

def tidal_velocity(tick):
    cycle_tick = tick % 200
    if cycle_tick < 100:
        return 0.5 - cycle_tick * (1.0 / 100)
    else:
        return -0.5 + (cycle_tick - 100) * (1.0 / 100)


# ── Least-cost routing field (mirrors calculate-cost-to-home) ─────────────────
# Cumulative difficulty-weighted cost from home outward.
# base_cost = abs(velocity) + 0.5 * max(0, depth slope), mirroring
# compute-travel-cost. Returns a numpy array indexed by patch position.

def calc_cost_to_home(n_patches, difficulty, depths, velocity):
    cost = np.full(n_patches, 1e6)
    cost[n_patches - 1] = 0.0           # home patch has zero cost
    for i in range(n_patches - 2, -1, -1):
        v_resistance = abs(velocity)
        slope        = max(0.0, depths[i + 1] - depths[i])
        base_cost    = v_resistance + 0.5 * slope
        cost[i]      = cost[i + 1] + base_cost * difficulty
    return cost


# ── Hydrodynamic difficulty (mirrors calculate-difficulty-landward) ───────────

def calc_difficulty(velocity, weight, max_weight,
                    v_max=MAX_SEAWARD_VELOCITY,
                    v_min=MAX_LANDWARD_VELOCITY,
                    k=0.75):
    norm_v = (velocity - v_min) / (v_max - v_min)
    df_raw = (norm_v / (weight / max_weight)) ** k
    return np.clip(1 + 9 * df_raw, 1, 10)


# ── Swimming speed (mirrors calculate-swimming-speed-landward) ────────────────

def calc_swim_speed(energy, difficulty, max_speed, min_speed,
                    velocity, prev_speed, swim_efficiency,
                    k=-0.75, max_change=0.5):
    energy_factor   = energy / 100.0
    velocity_impact = k * velocity * 300
    desired = np.clip(
        (max_speed * energy_factor / difficulty) + velocity_impact,
        min_speed, max_speed
    )
    max_change_eff = max_change * swim_efficiency
    delta = desired - prev_speed
    speed = prev_speed + np.sign(delta) * min(abs(delta), max_change_eff)
    return max(speed, min_speed)


# ── Swimming energy cost (mirrors calculate-swim-energy) ─────────────────────

def calc_swim_energy(difficulty, swim_efficiency, swim_base, beta=0.75):
    energy_mult = difficulty * (1.0 / swim_efficiency)
    return swim_base * (energy_mult ** beta)


# ── Schooling — mirrors the BOIDS block inside migrate-landward ───────────────
# Called from inside migrate_landward_lc, not from the main loop.
# Agents in STST never reach this function.

def school_agent(agent_id, agents):
    focal = agents.loc[agents["id"] == agent_id].iloc[0]
    mates = agents[
        (agents["id"] != agent_id) &
        (agents["start_migration"]) &
        (agents["landward_migration"] == focal["landward_migration"])
    ]

    if mates.empty:
        return focal["heading"], focal["speed"]

    dists   = np.abs(mates["x"] - focal["x"])
    nn_dist = dists.min()
    nn_x    = mates.loc[dists.idxmin(), "x"]

    new_heading = focal["heading"]
    new_speed   = focal["speed"]

    # separate / cohere + align
    if nn_dist < focal["minimum_separation"]:
        if nn_x > focal["x"]:
            # separate: slow down if neighbor is ahead
            if focal["speed"] > focal["min_speed"]:
                new_speed = focal["speed"] - 0.1
        else:
            # separate: turn away
            new_heading = focal["heading"] + (
                focal["separate_coefficient"] / 100
            ) * np.sign(focal["x"] - nn_x) * 10
    else:
        # cohere: turn toward centroid
        centroid_x  = mates["x"].mean()
        new_heading = focal["heading"] + (
            focal["cohere_coefficient"] / 100
        ) * np.sign(centroid_x - focal["x"]) * 10
        # align: match mean heading
        mean_heading = mates["heading"].mean()
        new_heading  = new_heading + (
            focal["align_coefficient"] / 100
        ) * (mean_heading - new_heading)

    # adjust-speed
    if mates["speed"].max() > focal["speed"] + 0.1:
        new_speed = focal["speed"] + focal["speed"] * 0.1
    elif focal["speed"] > focal["min_speed"]:
        new_speed = focal["speed"] - focal["speed"] * 0.1
    new_speed = max(new_speed, focal["min_speed"])

    return new_heading % 360, new_speed


# ── Landward movement with schooling nested inside ───────────────────────────
# Mirrors migrate-landward: BOIDS first, then move along least-cost path.

def migrate_landward_lc(agent_id, agents):
    new_heading, new_speed = school_agent(agent_id, agents)
    focal  = agents.loc[agents["id"] == agent_id].iloc[0]
    steps  = max(1, int(np.ceil(new_speed / 3)))
    new_x  = min(focal["x"] + steps, focal["home_patch"])
    return new_x, new_heading, new_speed


# ── STST gatekeeper (mirrors selective-tidal-stream-transport-landward) ────────
# Case 1: cooldown active — continue drifting, no schooling
# Case 2: flow > swim capacity — enter STST, drift, no schooling
# Case 3: swim capacity OK — exit STST, call migrate_landward_lc (with schooling)

def stst_landward(agent_id, agents, velocity, current_tick):
    focal  = agents.loc[agents["id"] == agent_id].iloc[0]
    S_swim = focal["speed"] / 300.0

    # Case 1: cooldown active — continue drifting
    if focal["in_stst"] and (current_tick - focal["stst_start_tick"] < 1):
        drift_dist = abs(velocity) * (1 - focal["swim_efficiency"])
        new_x      = min(focal["x"] + int(np.ceil(drift_dist)), focal["home_patch"])
        return {
            "x": new_x, "heading": focal["heading"], "speed": focal["speed"],
            "in_stst": True, "stst_start_tick": focal["stst_start_tick"],
            "energy": focal["energy"]
        }

    # Case 2: flow exceeds swim capacity — enter STST
    if abs(velocity) > S_swim and S_swim <= focal["min_speed"]:
        drift_dist = abs(velocity) * (1 - focal["swim_efficiency"])
        new_x      = min(focal["x"] + int(np.ceil(drift_dist)), focal["home_patch"])
        return {
            "x": new_x, "heading": focal["heading"], "speed": focal["speed"],
            "in_stst": True, "stst_start_tick": current_tick,
            "energy": focal["energy"]
        }

    # Case 3: swim capacity sufficient — exit STST, migrate (with schooling)
    new_x, new_heading, new_speed = migrate_landward_lc(agent_id, agents)
    E_swim     = calc_swim_energy(focal["difficulty_factor"],
                                  focal["swim_efficiency"], focal["swim_base"])
    new_energy = max(focal["energy"] - E_swim, 0.0)
    return {
        "x": new_x, "heading": new_heading, "speed": new_speed,
        "in_stst": False, "stst_start_tick": focal["stst_start_tick"],
        "energy": new_energy
    }

17.5.3 Simulation loop

log_ticks  = []
max_weight = agents["weight"].max()

for tick in range(1, n_ticks + 1):

    # 1. Tidal velocity
    v                 = tidal_velocity(tick)
    patches["velocity"] = v

    # 2. Rebuild least-cost field
    mean_diff              = calc_difficulty(v, agents["weight"].mean(), max_weight)
    patches["cost_to_home"] = calc_cost_to_home(
        n_patches, mean_diff, patches["depth"].values, v
    )

    # 3. Per-agent difficulty and swimming speed (vectorised)
    agents["difficulty_factor"] = calc_difficulty(
        v, agents["weight"].values, max_weight
    )
    agents["speed"] = [
        calc_swim_speed(
            agents.at[i, "energy"],
            agents.at[i, "difficulty_factor"],
            agents.at[i, "max_speed"],
            agents.at[i, "min_speed"],
            v,
            agents.at[i, "prev_speed"],
            agents.at[i, "swim_efficiency"],
        )
        for i in agents.index
    ]
    agents["prev_speed"] = agents["speed"]

    # 4. STST gatekeeper — calls migrate_landward_lc (with schooling) for active agents
    active = (
        agents["start_migration"] &
        agents["landward_migration"] &
        ~agents["at_destination"]
    )
    for i in agents.index[active]:
        result = stst_landward(
            agent_id     = agents.at[i, "id"],
            agents       = agents,
            velocity     = v,
            current_tick = tick,
        )
        agents.at[i, "x"]              = result["x"]
        agents.at[i, "heading"]        = result["heading"]
        agents.at[i, "speed"]          = result["speed"]
        agents.at[i, "in_stst"]        = result["in_stst"]
        agents.at[i, "stst_start_tick"] = result["stst_start_tick"]
        agents.at[i, "energy"]         = result["energy"]

    # 5. Swimming energy cost (baseline cost for all ticks including drift)
    agents["E_swim"] = [
        calc_swim_energy(
            agents.at[i, "difficulty_factor"],
            agents.at[i, "swim_efficiency"],
            agents.at[i, "swim_base"],
        )
        for i in agents.index
    ]
    agents["energy"] = np.maximum(agents["energy"] - agents["E_swim"], 0)

    # 6. Arrival check
    agents["at_destination"]    = agents["x"] >= agents["home_patch"]
    agents["landward_migration"] = ~agents["at_destination"]

    log_ticks.append({
        "tick"       : tick,
        "velocity"   : v,
        "mean_x"     : agents["x"].mean(),
        "mean_energy": agents["energy"].mean(),
        "n_stst"     : agents["in_stst"].sum(),
        "n_arrived"  : agents["at_destination"].sum(),
    })

results = pd.DataFrame(log_ticks)

17.5.4 Visualizing outputs

import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 3, figsize=(14, 4))

# Migration progress with tidal velocity overlay
ax = axes[0]
ax.plot(results["tick"], results["mean_x"],
        color="#2c7bb6", linewidth=1.5, label="Mean position")
ax.plot(results["tick"], results["velocity"] * 20 + 30,
        color="#d7191c", linewidth=0.8, linestyle="--",
        label="Tidal velocity (scaled)")
ax.set_title("Migration progress with tidal signal")
ax.set_xlabel("Tick")
ax.set_ylabel("Mean upstream position")
ax.legend(fontsize=8)

# Agents in STST per tick
ax = axes[1]
ax.fill_between(results["tick"], results["n_stst"],
                color="#fdae61", alpha=0.6)
ax.plot(results["tick"], results["n_stst"],
        color="#f46d43", linewidth=0.8)
ax.set_title("Agents in STST per tick")
ax.set_xlabel("Tick")
ax.set_ylabel("Count in passive drift")

# Energy over time
ax = axes[2]
ax.plot(results["tick"], results["mean_energy"],
        color="#1a9641", linewidth=1.5)
ax.set_title("Mean agent energy over time")
ax.set_xlabel("Tick")
ax.set_ylabel("Energy")

plt.tight_layout()
plt.show()