Source file odeiv.ml

(* gsl-ocaml - OCaml interface to GSL                       *)
(* Copyright (©) 2002-2012 - Olivier Andrieu                *)
(* Distributed under the terms of the GPL version 3         *)


let () = Error.init ()

type system 

external _alloc : 
  (float -> float array -> float array -> unit) ->
  ?jacobian:(float -> float array -> Matrix.matrix -> float array -> unit) ->
  int -> system = "ml_gsl_odeiv_alloc_system"

external _free : system -> unit = "ml_gsl_odeiv_free_system"

let make_system f ?jac dim =
  let s = _alloc f ?jacobian:jac dim in
  Gc.finalise _free s ;
  s

type step
type step_kind = 
  | RK2
  | RK4
  | RKF45 
  | RKCK
  | RK8PD 
  | RK2IMP
  | RK2SIMP
  | RK4IMP 
  | BSIMP
  | GEAR1 
  | GEAR2

external _step_alloc : step_kind -> dim:int -> step = "ml_gsl_odeiv_step_alloc"
external _step_free : step -> unit = "ml_gsl_odeiv_step_free"

let make_step kind ~dim = 
  let s = _step_alloc kind ~dim in
  Gc.finalise _step_free s ;
  s

external step_reset : step -> unit = "ml_gsl_odeiv_step_reset"
external step_name  : step -> string = "ml_gsl_odeiv_step_name"
external step_order : step -> int = "ml_gsl_odeiv_step_order"

external step_apply : step -> t:float -> h:float -> y:float array ->
  yerr:float array -> ?dydt_in:float array -> ?dydt_out:float array ->
    system -> unit
	= "ml_gsl_odeiv_step_apply_bc" "ml_gsl_odeiv_step_apply"


type control

external _control_free : control -> unit = "ml_gsl_odeiv_control_free"
external _control_standard_new : eps_abs:float -> eps_rel:float ->
  a_y:float -> a_dydt:float -> control = "ml_gsl_odeiv_control_standard_new"
external _control_y_new  : eps_abs:float -> eps_rel:float -> control = "ml_gsl_odeiv_control_y_new"
external _control_yp_new : eps_abs:float -> eps_rel:float -> control = "ml_gsl_odeiv_control_yp_new"
external _control_scaled_new : eps_abs:float -> eps_rel:float ->
  a_y:float -> a_dydt:float -> scale_abs:float array -> control = "ml_gsl_odeiv_control_scaled_new"

let make_control_standard_new ~eps_abs ~eps_rel ~a_y ~a_dydt =
  let c = _control_standard_new ~eps_abs ~eps_rel ~a_y ~a_dydt in
  Gc.finalise _control_free c ;
  c
let make_control_y_new ~eps_abs ~eps_rel =
  let c = _control_y_new ~eps_abs ~eps_rel in
  Gc.finalise _control_free c ;
  c
let make_control_yp_new ~eps_abs ~eps_rel =
  let c = _control_yp_new ~eps_abs ~eps_rel in
  Gc.finalise _control_free c ;
  c
let make_control_scaled_new ~eps_abs ~eps_rel ~a_y ~a_dydt ~scale_abs =
  let c = _control_scaled_new ~eps_abs ~eps_rel ~a_y ~a_dydt ~scale_abs in
  Gc.finalise _control_free c ;
  c

external control_name : control -> string = "ml_gsl_odeiv_control_name"

type hadjust =
  | HADJ_DEC
  | HADJ_NIL
  | HADJ_INC

external control_hadjust : control -> step -> y:float array ->
  yerr:float array -> dydt:float array -> h:float -> hadjust * float
      = "ml_gsl_odeiv_control_hadjust_bc" "ml_gsl_odeiv_control_hadjust"

type evolve

external _evolve_alloc : int -> evolve = "ml_gsl_odeiv_evolve_alloc"
external _evolve_free  : evolve -> unit = "ml_gsl_odeiv_evolve_free"
let make_evolve dim =
  let e = _evolve_alloc dim in
  Gc.finalise _evolve_free e ;
  e

external evolve_reset : evolve -> unit = "ml_gsl_odeiv_evolve_reset"

external evolve_apply : evolve -> control -> step -> system ->
  t:float -> t1:float -> h:float -> y:float array -> float * float
      = "ml_gsl_odeiv_evolve_apply_bc" "ml_gsl_odeiv_evolve_apply"