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. v0 = 99.9026463841554175 *)

let nonlinear a t_0 =
  let v_0 = 2. *. omega *. a and v = ref 0. and t = ref 0. and t_init = ref t_0 and t_int = ref 0. in
  let v_init = ref v_0 in
  let p = open_out "Fichier.don" and q = open_out "Fichier2.don" and zout = open_out "trace.don" in
    Printf.fprintf p "%d %f\n" 0 !t_init;
    Printf.fprintf q "%d %f\n" 0 !v_init;
    t:= !t +. 0.0002;
    while ((f !t !t_init !v_init a) > 0.000000) do
      print_float (f !t !t_init !v_init a);
      print_newline();
      t:= !t +. 0.0002;
      Printf.fprintf zout "%f %f\n" !t ((f !t !t_init !v_init a) +. a*. cos(omega*. !t));
    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;
    Printf.fprintf zout "%f %f\n" !t (a *. cos(omega*. !t));
    Printf.fprintf p "%d %f\n" 1 (( mod_float !t (2. *. pi))/. (2. *. pi));
    Printf.fprintf q "%d %f\n" 1 !v;
    for i=2 to 500 do
      t:= !t +. 0.0002;
      while (f !t !t_init !v_init a > 0.000000) do
	t:= !t +. 0.0002;
	Printf.fprintf zout "%f %f\n" !t ((f !t !t_init !v_init a) +. a*. cos(omega*. !t));
      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;
      Printf.fprintf zout "%f %f\n" !t (a*. cos(omega*. !t));
      Printf.fprintf p "%d %f\n" i (( mod_float !t (2. *. pi))/. (2. *. pi));
      Printf.fprintf q "%d %f\n" i !v;
    done;
    close_out q;
    close_out zout;
    close_out p;;


let main () =
  nonlinear 0.00053 0.;;

main ();;