symplectic.nim

# Symplectic integrators for separable Hamiltonian systems
#
# (c) 2018-2021 [email protected] (Ian Smith), for licencing see the LICENCE file

import std/strutils
import math

type
  ArgumentError* = object of ValueError

func log_error*(error: float): float =
  10.0 * log10(if abs(error) >= 1e-36: abs(error) else: 1e-36)

proc solve*[T](args: seq[string], p: var T, u_q, u_p: proc(p: var T, cd: float), output: proc(p: T, t: float)) =
  let
    order: int64 = parseInt(args[1])
    h: float = parseFloat(args[2])
    steps: int64 = parseInt(args[3])
    z1: float = 1.0 / (4.0 - pow(4.0, (1.0 / 3.0)))
    y1: float = 1.0 / (4.0 - pow(4.0, (1.0 / 5.0)))
    x1: float = 1.0 / (4.0 - pow(4.0, (1.0 / 7.0)))
    z0: float = 1.0 - 4.0 * z1
    y0: float = 1.0 - 4.0 * y1
    x0: float = 1.0 - 4.0 * x1
    cd: seq[float] = case order
    of 2:  # Stormer-Verlet
      @[0.5 * h, h]
    of 4:  # m4r35n357_4
      @[0.5 * h * z1, h * z1, h * z1, h * z1, 0.5 * h * (z1 + z0), h * z0]
    of 6:  # m4r35n357_6
      @[0.5 * h * z1 * y1, h * z1 * y1, h * z1 * y1, h * z1 * y1, 0.5 * h * (z1 + z0) * y1, h * z0 * y1, 0.5 * h * (z0 + z1) * y1, h * z1 * y1, h * z1 * y1, h * z1 * y1, h * z1 * y1, h * z1 * y1, h * z1 * y1, h * z1 * y1, 0.5 * h * (z1 + z0) * y1, h * z0 * y1, 0.5 * h * (z0 + z1) * y1, h * z1 * y1, h * z1 * y1, h * z1 * y1, 0.5 * h * z1 * (y1 + y0), h * z1 * y0, h * z1 * y0, h * z1 * y0, 0.5 * h * (z1 + z0) * y0, h * z0 * y0]
    of 8:  # m4r35n357_8
      @[0.5 * h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, 0.5 * h * (z1 + z0) * y1 * x1, h * z0 * y1 * x1, 0.5 * h * (z0 + z1) * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, 0.5 * h * (z1 + z0) * y1 * x1, h * z0 * y1 * x1, 0.5 * h * (z0 + z1) * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, 0.5 * h * z1 * (y1 + y0) * x1, h * z1 * y0 * x1, h * z1 * y0 * x1, h * z1 * y0 * x1, 0.5 * h * (z1 + z0) * y0 * x1, h * z0 * y0 * x1, 0.5 * h * (z0 + z1) * y0 * x1, h * z1 * y0 * x1, h * z1 * y0 * x1, h * z1 * y0 * x1, 0.5 * h * z1 * (y1 + y0) * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, 0.5 * h * (z1 + z0) * y1 * x1, h * z0 * y1 * x1, 0.5 * h * (z0 + z1) * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, 0.5 * h * (z1 + z0) * y1 * x1, h * z0 * y1 * x1, 0.5 * h * (z0 + z1) * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, 0.5 * h * (z1 + z0) * y1 * x1, h * z0 * y1 * x1, 0.5 * h * (z0 + z1) * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, 0.5 * h * (z1 + z0) * y1 * x1, h * z0 * y1 * x1, 0.5 * h * (z0 + z1) * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, 0.5 * h * z1 * (y1 + y0) * x1, h * z1 * y0 * x1, h * z1 * y0 * x1, h * z1 * y0 * x1, 0.5 * h * (z1 + z0) * y0 * x1, h * z0 * y0 * x1, 0.5 * h * (z0 + z1) * y0 * x1, h * z1 * y0 * x1, h * z1 * y0 * x1, h * z1 * y0 * x1, 0.5 * h * z1 * (y1 + y0) * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, 0.5 * h * (z1 + z0) * y1 * x1, h * z0 * y1 * x1, 0.5 * h * (z0 + z1) * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, 0.5 * h * (z1 + z0) * y1 * x1, h * z0 * y1 * x1, 0.5 * h * (z0 + z1) * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, h * z1 * y1 * x1, 0.5 * h * z1 * y1 * (x1 + x0), h * z1 * y1 * x0, h * z1 * y1 * x0, h * z1 * y1 * x0, 0.5 * h * (z1 + z0) * y1 * x0, h * z0 * y1 * x0, 0.5 * h * (z0 + z1) * y1 * x0, h * z1 * y1 * x0, h * z1 * y1 * x0, h * z1 * y1 * x0, h * z1 * y1 * x0, h * z1 * y1 * x0, h * z1 * y1 * x0, h * z1 * y1 * x0, 0.5 * h * (z1 + z0) * y1 * x0, h * z0 * y1 * x0, 0.5 * h * (z0 + z1) * y1 * x0, h * z1 * y1 * x0, h * z1 * y1 * x0, h * z1 * y1 * x0, 0.5 * h * z1 * (y1 + y0) * x0, h * z1 * y0 * x0, h * z1 * y0 * x0, h * z1 * y0 * x0, 0.5 * h * (z1 + z0) * y0 * x0, h * z0 * y0 * x0]
    else:
      echo "invalid order value ", order, ", should be 2, 4, 6, or 8"
      raise newException(ArgumentError, "Invalid integrator")
  output(p, 0.0)
  for step in 1 .. steps:
    for i in countup(0, cd.len - 1):
      if (i mod 2) == 0: u_q(p, cd[i]) else: u_p(p, cd[i])
    for i in countdown(cd.len - 2, 0):
      if (i mod 2) == 0: u_q(p, cd[i]) else: u_p(p, cd[i])
    output(p, step.float * h)