let pi = acos(-1.);;
let omega = 2. *. pi *. 30.;;


let absfloat = function
  |  x when x > 0. -> x
  | x -> -1. *. x;;

let f t t_1 v_1 a = a*. (cos(omega *. t_1) -. cos(omega *. t)) +. v_1 *. (t -. t_1 ) -. 5. *. (t -. t_1 )**2. ;;

let dichot f a t_0 t_1 v_1 precision =
  let x = ref (t_0 +. 2. *. precision) and y = ref t_1 in
  let x_int = ref 0.  in
    while ( absfloat (!x-. !y) > precision) do
      x_int := (!x+. !y)/. 2.;
      if ((f !x t_0 v_1 a)*.(f !x_int t_0 v_1 a) > 0.) then x:=!x_int
      else
	y:=!x_int;
    done;
    !x;;


(* On part de v_0, et A est donné. t_0 = 0. *)

let main () =
  let v = ref 0. and t = ref 0. and t_init = ref 0. and a= ref 0.00040 and p = open_out "bifurc.don" and t_int = ref 0. in
    for j=1 to 4000 do
      let v_0 = 2. *. omega *. !a in
      let v_init = ref v_0 in
	t:= !t +. 0.0002;
(*recherche ci-dessous de la borne sup de l'intervalle pour la dichotomie *)
	while ((f !t !t_init !v_init !a) > 0.000000) do
	  t:= !t +. 0.0002;
	done;
	t_int := (dichot f !a !t_init !t  !v_init  0.000001);
	t:=!t_int;
	v := -1. *. omega *. !a *. (1. +. 0.25) *. sin(omega*. !t) -. 0.25*. !v_init +. 0.25*. 10. *. (!t -. !t_init);
	t_init := !t;
	v_init := !v;
	for i=2 to 300 do
	  t:= !t +. 0.0002;
	  while (f !t !t_init !v_init !a > 0.000000) do
	    t:= !t +. 0.0002;
	  done;
	  t_int := (dichot f !a !t_init !t  !v_init  0.000001);
	  t := !t_int;
	  v := -1. *. omega *. !a *. (1. +. 0.25) *. sin(omega*. !t) -. 0.25*. !v_init +. 0.25*. 10. *. (!t -. !t_init);
	  t_init := !t;
	  v_init := !v;
	  if (i>280) then
	    Printf.fprintf p "%f %f\n" !a (( mod_float !t (2. *. pi))/. (2. *. pi));
	done;
	a := !a +. 0.0000001;
	t := 0.;
	t_init := 0.;
	t_int := 0.;
	v := 0.;
    done;
    close_out p;;

main ();;