#set page(width: 120mm, height: 137mm, margin: 8mm, fill: rgb("#F5F3ED"))
#set text(font: "Libertinus Serif", fill: rgb("#183343"))
#set par(spacing: 0pt)

#let damping = float(sys.inputs.at("damping", default: "0.35"))
#assert(damping >= 0 and damping <= 2, message: "damping must be 0..2")
#let field(x, y) = (y, -x - damping * y)
#let point(x, y) = (50mm + x * 13mm, 50mm - y * 13mm)

#let arrow(x, y, dx, dy) = {
  let norm = calc.sqrt(dx * dx + dy * dy)
  if norm > 0 {
    let (px, py) = point(x, y)
    let (u, v) = (dx / norm, -dy / norm)
    let shaft = 3.5mm
    let (ex, ey) = (px + u * shaft, py + v * shaft)
    place(top + left, dx: px, dy: py,
      line(end: (u * shaft, v * shaft), stroke: 0.6pt + rgb("#087F8C")))
    place(top + left, polygon(
      (ex, ey),
      (ex - u * 1.2mm + v * 0.7mm, ey - v * 1.2mm - u * 0.7mm),
      (ex - u * 1.2mm - v * 0.7mm, ey - v * 1.2mm + u * 0.7mm),
      fill: rgb("#087F8C"), stroke: none,
    ))
  }
}

#text(size: 17pt, weight: "bold")[A computed phase portrait]
#v(2mm)
$ dot(x) = y, quad dot(y) = -x - #damping y $
#v(3mm)

#block(width: 100mm, height: 100mm, {
  place(top + left, dx: 4mm, dy: 50mm, line(end: (92mm, 0mm), stroke: 0.35pt + rgb("#B5C3C6")))
  place(top + left, dx: 50mm, dy: 4mm, line(end: (0mm, 92mm), stroke: 0.35pt + rgb("#B5C3C6")))
  for x in range(-3, 4) {
    for y in range(-3, 4) {
      let (dx, dy) = field(x, y)
      arrow(x, y, dx, dy)
    }
  }
  // Explicit Euler samples: an illustrative path, not a high-accuracy solver.
  let (x, y) = (2.7, 0.0)
  let trajectory = (point(x, y),)
  for _ in range(450) {
    let (dx, dy) = field(x, y)
    (x, y) = (x + 0.025 * dx, y + 0.025 * dy)
    trajectory.push(point(x, y))
  }
  place(top + left, curve(stroke: 1pt + rgb("#B65C32"),
    curve.move(trajectory.first()),
    ..trajectory.slice(1).map(p => curve.line(p))))
  place(top + left, dx: 96mm, dy: 48mm, text(size: 9pt)[$x$])
  place(top + left, dx: 52mm, dy: 0mm, text(size: 9pt)[$y$])
})