(* hisui1.ml - Don Yang (uguu.org)
   Mandelbrot set renderer

   ocaml hisui1.ml

   10/23/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 render_point z0 =
   let rec recursive_render_point z i =
      if Complex.norm z > 2.0 || i > iteration_count
      then i
      else recursive_render_point (Complex.add (Complex.mul z z) z0) (i + 1) in
   recursive_render_point z0 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 interesting x0 x1 y0 y1 =
   let series =
      List.map (fun x -> float x)
         [render_point {Complex.re = x0; Complex.im = y0};
          render_point {Complex.re = x0; Complex.im = y1};
          render_point {Complex.re = x1; Complex.im = y0};
          render_point {Complex.re = x1; Complex.im = y1};
          render_point {Complex.re = (x0 +. x1) *. 0.5;
                        Complex.im = (y0 +. y1) *. 0.5}] in
   let average =
      (List.fold_left (fun sum x -> sum +. x) 0.0 series) /. 5.0 in
   ((List.fold_left (fun sum_sq x -> sum_sq +. x) 0.0
      (List.rev_map (fun d -> (d -. average) *. (d -. average))
                    series)) /. 5.0) >
   ((float iteration_count) *. 4.6);;

(* Render image *)
let render_image x0 x1 xsteps y0 y1 ysteps =
   for y = 1 to ysteps do
      for x = 1 to xsteps do
         print_char (iter_char (render_point
            {Complex.re = (x1 -. x0) *. (float x) /.
               (float (xsteps - 1)) +. x0;
             Complex.im = (y0 -. y1) *. (float y) /.
               (float (ysteps - 1)) +. y1}))
      done;
      print_newline ()
   done;;

(* 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
   if interesting x0 (x0 +. s *. aspect_ratio) y0 (y0 +. s)
   then render_image x0 (x0 +. s *. aspect_ratio) x_steps
                     y0 (y0 +. s) y_steps
   else render_random_image();;

render_random_image();;