import os
#import math
#import std/strutils
#import dual
#import symplectic

type Parameters = tuple[r_min: float, r_max: float, theta_max: float, E: float, L: float, Q: float]

type Vector3 = tuple[a: float, b: float, c: float]

type Matrix3x3 = tuple[a: float, b: float, c: float, d: float, e: float, f: float, g: float, h: float, i: float]

func invert(m: Matrix3x3): Matrix3x3 =
    c: Matrix3x3 = ((m.e * m.i - m.f * m.h), -(m.d * m.i - m.f * m.g),  (m.d * m.h - m.e * m.g),
                   -(m.b * m.i - m.c * m.h),  (m.a * m.i - m.c * m.g), -(m.a * m.h - m.b * m.g),
                    (m.b * m.f - m.c * m.e), -(m.a * m.f - m.c * m.d),  (m.a * m.e - m.b * m.d))
    d: float = m.a * c.a + m.b * c.b + m.c * c.c
    assert(d != 0.0);
    return (c.a / d, c.d / d, c.g / d,
            c.b / d, c.e / d, c.h / d,
            c.c / d, c.f / d, c.i / d)

func mv_mult(m: Matrix3x3, v: Vector3): Vector3 =
    return (m.a * v.a + m.b * v.b + m.c * v.c,
            m.d * v.a + m.e * v.b + m.f * v.c,
            m.g * v.a + m.h * v.b + m.i * v.c)

func R(p: Parameters, float r, E: Dual, L: Dual, Q: Dual): Dual =
  r2: float = r * r
  delta: Dual = r2 - 2.0 * p.m * r + p.a2
  P: Dual = E * (r2 + p.a2) - L * p.a)
  K: dual = (L - E * p.a) * (L - E * p.a) + Q
  P * P - delta * (p.mu2 * r2 + K)

func dR_dr(p: Parameters, float r, E: Dual, L: Dual, Q: Dual): Dual =
  r2: float = r * r
  delta: Dual = r2 - 2.0 * p.m * r + p.a2
  P: Dual = E * (r2 + p.a2) - L * p.a)
  K: dual = (L - E * p.a) * (L - E * p.a) + Q
  4.0 * E * r * P - (2.0 * r - 2.0 * p.m) * (p.mu2 * r2 + K) - 2.0 * p.mu2 * r * delta

proc plot(p: Parameters) =
  echo p.q_r * sin(p.q_phi), " ", p.q_r * cos(p.q_phi), " ", 0.0, " ", t, " ", 10.0 * log10(if e >= 1e-36 : e else: 1e-36)

proc generate(p: Parameters):
  echo "Generate"

let
  args: seq[string] = commandLineParams()
var
  p: Parameters = (parseFloat(args[1]), parseFloat(args[2]), parseFloat(args[3]), 1.0, 0.0, 0.0)
generate(p)