module Deriv = struct 
type cnst = V of float | N of string
and  expr = Add of expr * expr
          | Mul of expr * expr
          | Sub of expr * expr
          | Div of expr * expr
          | Sin of expr | Cos of expr
          | Var of string
          | Const of cnst
let rec deriv e v = match e with
    Add (x, y) -> Add (deriv x v, deriv y v)
  | Sub (x, y) -> Sub (deriv x v, deriv y v)
  | Mul (x, y) -> Add (Mul (x, deriv y v), Mul (deriv x v, y))
  | Div (x, y) -> Sub (Div (deriv x v, y), Div (Mul (x, deriv y v), Mul (y, y)))
  | Sin x -> Mul (Cos x, (deriv x v))
  | Cos x -> Sub (Const (V 0.), Mul (Sin x, (deriv x v)))
  | Var string -> Const (V (if string = v then 1. else 0.))
  | Const c -> Const (V 0.)       
let rec ev e (env: string -> float) = match e with
    Add(x, y) -> (ev x env) +. (ev y env)
  | Sub(x, y) -> (ev x env) -. (ev y env)
  | Mul(x, y) -> (ev x env) *. (ev y env)
  | Div(x, y) -> (ev x env) /. (ev y env)
  | Sin x -> sin (ev x env)
  | Cos x -> cos (ev x env)
  | Var str -> env str
  | Const c -> (match c with V n -> n | N str -> env str)
let (+) a b = Add (a, b) and (-) a b = Sub (a, b)
and ( * ) a b = Mul (a, b) and (/) a b = Div (a, b)
and nc st = Const (N st) and lc f = Const (V f)
let rec pr e = let ps = print_string in
  match e with
    Add (a, b) -> ps "("; pr a; ps " + "; pr b; ps ")"
  | Sub (a, b) -> ps "("; pr a; ps " - "; pr b; ps ")" 
  | Mul (a, b) -> ps "("; pr a; ps " * "; pr b; ps ")"
  | Div (a, b) -> ps "("; pr a; ps " / "; pr b; ps ")"
  | Sin a -> print_string "sin ("; pr a; print_string ")"
  | Cos a -> print_string "cos ("; pr a; print_string ")"
  | Var s -> print_string s
  | Const c -> match c with V f -> print_float f | N s -> print_string s
end;;
(* open Deriv;;
let xy = (Mul ((Var "x"), (Var "y")))
and m str = if str = "x" then 2.3 else if str = "y" then 7. else 1. /. 0. in
Printf.printf "%e\n" (ev (deriv xy "x") m);; *)
(* http://en.wikipedia.org/wiki/Kerr-Newman_metric *)
open Deriv;;
let kn bM bQ bJ = let bM = lc bM and bQ = lc bQ and bJ = lc bJ 
and bSI = lc 2.99792458e8, (* Speed of Light *)
   lc 6.67384e-11, (* Gravitational Constant *)
   lc 8.9875517873681764e9 (* Coulomb's Constant *)
and gu = lc 1., lc 1., lc 1. in
let c, bG, ke = if 0<1 then gu else bSI (* choice of units *) in
let sq x = x * x and a = bJ / (bM * c) in
let rs = lc 2. * bG * bM / sq c
and rQ2 = sq bQ * bG * ke / sq (sq c) in
fun r th dt dr dth dp ->
  let r2 = sq r and a2 = sq a and sth2 = sq (Sin th) in
  let del = r2 - rs * r + a2 + rQ2
  and rh2 = r2 + a2 * (lc 1. - sth2) in (
  sq (c * dt - a * sth2 * dp) * (del / rh2)
  - (sq dr / del + sq dth) * rh2
  - sq ((r2 + a2) * dp - a * c * dt) * sth2 / rh2) / sq c
in Printf.printf "Boo\n";
let bh = kn 1. 0. 0. in
pr (deriv (bh (lc 2.2) (lc 0.2) (lc 0.01) (lc 0.014) (lc 0.009) (lc 0.008))
  "th");
print_newline ();
pr (deriv (Mul ((Var "s"), (Var "s"))) "s");
print_newline (); flush stdout