(* kohaku0.ml - Don Yang (uguu.org)
   Julia set renderer

   ocaml kohaku0.ml

   10/13/05
                                                                           *)

(* Constants *)
let iteration_count = 26;;
let x_steps = 79;;
let y_steps = 24;;
let aspect_ratio = 1.5;;


(* Render a point by exponentiating it repeatedly.  Returns number of
   iteratations it takes for the number to diverge, clamped by
   iteration_count.                                                        *)
let iterative_render_point z0 c =
   let mutable_z = z0 in let z = ref mutable_z in
   let mutable_i = 0 in let i = ref mutable_i in
   while Complex.norm !z <= 2.0 && !i <= iteration_count do
      z := Complex.add (Complex.mul !z !z) c;
      i := !i + 1
   done;
   !i;;

let rec recursive_render_point z c i =
   if Complex.norm z > 2.0 || i > iteration_count
   then i
   else recursive_render_point (Complex.add (Complex.mul z z) c) c (i + 1);;

let render_point z0 c =
   recursive_render_point z0 c 0;;


(* Convert iteration value to character *)
let iter_char i =
   if i >= 26 || i < 0
   then ' '
   else (Char.chr ((Char.code 'A') + i));;


(* Check if a particular set of coordinates is interesting by checking
   for high variety of values in sampled points.                           *)
let average series =
   (List.fold_left (fun sum x -> sum +. x) 0.0 series) /.
   (float (List.length series));;

let stddev series =
   let a = average series in
   (List.fold_left (fun sum_sq x -> sum_sq +. x) 0.0
      (List.rev_map (fun d -> (d -. a) *. (d -. a)) series)) /.
   (float (List.length series));;

let interesting x0 x1 y0 y1 c =
   (stddev
      (List.map (fun x -> float x)
         [render_point {Complex.re = x0; Complex.im = y0} c;
          render_point {Complex.re = x0; Complex.im = y1} c;
          render_point {Complex.re = x1; Complex.im = y0} c;
          render_point {Complex.re = x1; Complex.im = y1} c;
          render_point {Complex.re = (x0 +. x1) *. 0.5;
                        Complex.im = (y0 +. y1) *. 0.5} c])) >
   ((float iteration_count) *. 4.6);;


(* Generate a list of evenly spaced values *)
let iterative_generate_series x0 x1 steps =
   Array.to_list (Array.init steps
      (fun i -> x0 +. (x1 -. x0) *. (float i) /. (float (steps - 1))));;

let rec recursive_generate_series x0 x1 steps i =
   if i >= steps
   then []
   else (x0 +. (x1 -. x0) *. (float i) /. (float (steps - 1))) ::
        (recursive_generate_series x0 x1 steps (i + 1));;

let generate_series x0 x1 steps =
   recursive_generate_series x0 x1 steps 0;;


(* Render image *)
let render_scanline x0 x1 xsteps y c =
   List.iter
      (fun x -> print_char (iter_char
         (render_point {Complex.re = x; Complex.im = y} c)))
      (generate_series x0 x1 xsteps);
   print_newline();;

let render_image x0 x1 xsteps y0 y1 ysteps c =
   List.iter
      (fun y -> render_scanline x0 x1 xsteps y c)
      (generate_series y1 y0 ysteps);;


(* Random numbers *)
Random.self_init();;

let random_range x0 x1 =
   (Random.float (x1 -. x0)) +. x0;;

let random_complex () =
   {Complex.re = random_range (-1.0) 1.0;
    Complex.im = random_range (-1.0) 1.0};;


(* Generate coordinates repeatedly until an interesting one is found,
   then render image.                                                      *)
let rec render_random_image () =
   let x0 = random_range (-1.0) 1.0 in
   let y0 = random_range (-1.0) 1.0 in
   let s = random_range 0.1 1.0 in
   let c = random_complex() in
   if interesting x0 (x0 +. s *. aspect_ratio) y0 (y0 +. s) c
   then render_image x0 (x0 +. s *. aspect_ratio) x_steps
                     y0 (y0 +. s) y_steps c
   else render_random_image();;

render_random_image();;