open Hsim.Types
open Solvers

let of_array a = Bigarray.Array1.of_array Bigarray.Float64 Bigarray.c_layout a

type state = { lx : Zls.carray; i : bool }

let mu = 5.0
let x0 = 1.0
let y0 = 1.0

let fder _ y yd =
  yd.{0} <- y.{1}; yd.{1} <- (mu *. (1.0 -. (y.{0} *. y.{0})) *. y.{1}) -. y.{0}
let fzero _ _ _ = ()
let fout _ y     = of_array [| y.{0}; y.{1} |]
let step { i; lx } _ =
  let lx = if i then of_array [| x0; y0 |] else lx in
  of_array [| lx.{0}; lx.{1} |], { lx; i = false }

let cget s     = s.lx
let cset s lx = { s with lx }
let zset s _   = s
let yd         = Zls.cmake 2
let zout       = Zls.cmake 1
let zsize      = 1
let s_init     = { lx = of_array [| x0; y0 |]; i = true }
let reset _ _  = s_init
let jump _     = true

let van_der_pol () =
  HNode { state = s_init;
          fder = (fun _ _ y -> fder 0.0 y yd; yd);
          fzer = (fun _ _ y -> fzero 0.0 y zout; zout);
          fout = (fun s _ y -> fout s y);
          step;
          horizon = (fun _ -> max_float);
          cset; cget; zset; zsize; reset; jump }

(* let van_der_pol_prop_record () = *)
(*   let s_init =  *)
(*       { lx = of_array [| x0; y0 |]; i = true } in *)
(*   { name = "van_der_pol_prop"; *)
(*     s = s_init; *)
(*     fder = (fun s a y -> fder 0.0 y yd; yd); *)
(*     fzero = (fun s a y -> fzero 0.0 y zout; zout); *)
(*     fout = (fun s { lx } y -> 1.0 /. (abs_float (s.lx.{0} -. lx.{0}))); *)
(*     fstep = (fun s { lx } -> *)
(*       let v, s = fstep s () in *)
(*       s.lx.{0} -. lx.{0}, s); *)
(*     horizon = (fun s -> max_float); *)
(*     cset; cget; zset; csize; zsize; reset; jump } *)
(* let van_der_pol_with_assert () = *)
(*   Fun_hybrid *)
(*     { body = van_der_pol_record (); *)
(*       assertions = *)
(*         [Fun_hybrid { body = *)
(*                         van_der_pol_prop_record (); *)
(*                       assertions = [] }] } *)