# Parameter (E, L, Q) generation for Kerr spacetime
#
# Example:
#
# (c) 2018-2021 m4r35n357@gmail.com (Ian Smith), for licencing see the LICENCE file

import os
import math
#import std/strutils
import dual

type Parameters = tuple[r_min: float, r_max: float, theta_max: float, m: float, mu2: float, a: float, a2: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 =
    let
      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, r: float, E: Dual, L: Dual, Q: Dual): Dual =
  let
    r2: float = r * r
    delta: float = 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, r: float, E: Dual, L: Dual, Q: Dual): Dual =
  let
    r2: float = r * r
    delta: float = 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

func THETA(p: Parameters, float theta, E: Dual, L: Dual, Q: Dual): Dual =
  let
    cth2: float = cos(theta) * cos(theta)
    sth2: float = 1.0 - cth2
    Q - cth2 * ((p.mu2 - E * E) * p.a2 + L * L / sth2)

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)