(* Discrete-time integration functions *)

type saturation = Upper | Between | Lower

(* [output(0) = x0(0)] and [output(n) = output(n−1) + (k * t) * u(n−1)] *)
let node forward_euler(t)(k, x0, u) = output where
  rec output = x0 fby (output +. (k *. t) *. u)
	
(* forall n in Nat.
 * [output(0) = x0(0)]
 * [output(n) = output(n-1) + i * t * u(n)] *)   
let node backward_euler(t)(k, x0, u) = output where
  rec output = x0 -> pre(output) +. k *. t *. u
		     
(* forall n in Nat.
 * [output(n) = y(n)]
 * [x(0) = x0(0)]
 * [x(n) = y(n-1) + (k * t / 2) * u(n-1)]
 * [y(n) = x(n) + (k * t / 2) * u(n)] *)
let node trapezoidal_fixed(t)(k, x0, u) = output where
   rec x = x0 fby (y +. gain *. u)
   and y = x +. gain *. u
   and gain = k *. t /. 2.0
   and output = y

(* forall n in Nat.
 * [output(n) = y(n)]
 * [x(0) = x0(0)]
 * [x(n) = y(n-1)]
 * [pu(0) = u(0)]
 * [pu(n) = u(n-1)]
 * [y(n) = x(n) + (k * t / 2) * (u(n) + pu(n))] *)
let node trapezoidal_variable_step(t)(k, x0, u) = output where
  rec x = x0 fby y
  and y = x +. gain *. (u +. u fby u)
  and gain = k *. t /. 2.0
  and output = y


let node forward_euler_complete(t)(upper, lower, res, k, x0, u) =
   (output, sport, saturation) where
   rec sport = x0 fby (output +. k *. t *. u)
   and v = if res then x0 else sport
   and (output, saturation) =
     if v < lower then lower, Lower
     else if v > upper then upper, Upper else v, Between

let node backward_euler_complete(t)(upper, lower, res, k, x0, u) =
   (output, state_port, saturation) where
   rec state_port = x0 -> pre(output) +. k *. t *. u
   and v = if res then x0 else state_port
   and output, saturation =
     if v < lower then lower, Lower else if v > upper then upper, Upper
     else v, Between
      
let node int(t)(k, x0, u) = forward_euler(t)(k, x0, u)

let node derivative(h)(k, u) = 0.0 -> k *. (u -. pre(u)) /. h

(* Integration with reset is not provided in this library *)
(* because it is simple to get. E.g., the forward Euler with reset is *)
let node forward_euler_reset(t)(k, x0, res, u) = o where
  rec reset
	o = forward_euler(t)(k, x0, u)
      every res