(** Zélus: Hybrid system programming language

    - Model discrete systems and their continuous environment
    - Research language, design space for hybrid system modelers
    - Compilation to OCaml, execution with an off-the-shelf ODE solver
    - Developed by M. Pouzet and T. Bourke in the Inria PARKAS team *)











(** Plan:

    - Zélus, the synchronous part
    - Physical simulation, ODEs, and ODE solvers
    - Zélus, the hybrid part
    - Mixing discrete and continuous time, zero-crossing detection
    - Informal operational semantics
    - Low-level representation and runtime *)









(** Synchronous language kernel _à la_ Lustre:
    - programs are Mealy machines (outputs on each transition)
    - variables represent streams of values in time *)

let node incr x = y where
  y = x + 1

(*        x │ 8  3  2  7  5  3  …
    ────────┼─────────────────────
     incr x │ 9  4  3  8  6  4  …   *)





(** - we can use values of the previous instants with [pre] and
      initialize streams with [->] *)

let node accumulate x = z where
  rec z = x -> (pre z) + x

(*              x │ 1  2  5  2  5  3  …
    ──────────────┼─────────────────────
            pre x │    1  2  5  2  5  …
           0 -> x │ 0  2  5  2  5  3  …
       0 -> pre x │ 0  1  2  5  2  5  …
     accumulate x │ 1  3  8 10 15 18  …   *)

let node fib () = n where
  rec n = 0 -> pre (1 -> pre(n) + n)

(** - causality loops are forbidden ([rec x = x]) *)

(** - we can reset streams at will *)

let node stay x = y where            (* output the first value forever        *)
  rec y = x -> pre y

let node from x = y where            (* count up from the first value         *)
  rec y = x -> pre y + 1

let node loop x = y where
  rec y = reset from 0 every z
  and z = false -> pre y >= w
  and w = stay x

(*        x │ 2  _  _  _  _  _  …
    ────────┼─────────────────────
     loop x │ 0  1  2  0  1  2  …   *)








(** Calculus reminder:

    - Ordinary differential equations (ODEs), initial value problems
    - A little background on solvers *)








(** Already able to model physical behaviours! *)

let dt = 0.01                         (* Integration step                     *)
let g  = 9.81                         (* Gravitational constant               *)

let node f_integr (x0, x') = x where  (* Forward Euler integrator             *)
  rec x = x0 -> pre (x +. x' *. dt)

let node b_integr (x0, x') = x where  (* Backward Euler integrator            *)
  rec x = x0 -> (pre x) +. x' *. dt

let node bouncing_ball (p0, v0) = p where
  rec p = reset f_integr (q, v)    every z
  and v = reset b_integr (w, -. g) every z
  and q = p0 -> 0.0 and w = v0 -> -0.8 *. (pre v)
  and z = false -> (pre p) < 0.0







(** Cumbersome, and error-prone! *)

let node sincos () = (sin, cos) where
  rec sin = f_integr (0.0, cos)        (** [(dsin/dt) t =  cos t, sin 0 = 0]  *)
  and cos = f_integr (1.0, -. sin)     (** [(dcos/dt) t = -sin t, cos 0 = 1]  *)









(** Enter continuous-time constructs:
    - express values with initial value problems with [der] and [init] *)

let hybrid integr (x0, x') = x where
  der x = x' init x0

let hybrid time () = t where
  der t = 1.0 init 0.0

let hybrid falling_ball (p0, v0) = p where
  rec der p = v    init p0
  and der v = -. g init v0

(** How do we mix discrete time and continuous time together? *)




(** - mix discrete and continuous code with [up], [present], [reset]
      and [last] *)

let hybrid bouncing_ball (p0, v0) = p where
  rec der p = v    init p0 reset z -> 0.0
  and der v = -. g init v0 reset z -> -0.8 *. last v
  and     z = up(-. p)

let hybrid time_bounces (p0, v0) = (p, b) where
  rec      p = bouncing_ball (p0, v0)
  and  der t = 1.0 init 0.0
  and init b = 0.0
  and present up(-. p) -> do
    b = t -. last b
  done