let delta = sqrt epsilon_float
type vec = { x: float; y: float; z: float }
let ( *| ) s r = {x = s *. r.x; y = s *. r.y; z = s *. r.z}
let ( +| ) a b = {x = a.x +. b.x; y = a.y +. b.y; z = a.z +. b.z}
let ( -| ) a b = {x = a.x -. b.x; y = a.y -. b.y; z = a.z -. b.z}
let dot a b = a.x *. b.x +. a.y *. b.y +. a.z *. b.z
let unitise r = (1. /. sqrt (dot r r)) *| r
type ray = { orig: vec; dir: vec }
type sphere = { center: vec; radius: float }
type scene = Sphere of sphere | Group of sphere * scene list
let ray_sphere ray sphere =
  let v = sphere.center -| ray.orig in
  let b = dot v ray.dir in
  let disc = b *. b -. dot v v +. sphere.radius *. sphere.radius in
  if disc < 0. then infinity else
    let disc = sqrt disc in
    (fun t2 -> if t2 < 0. then infinity else
       ((fun t1 -> if t1 > 0. then t1 else t2) (b -. disc))) (b +. disc)
let intersect ray scene =
  let rec aux ((lambda, _) as hit) = function
      Sphere sphere ->
	let lambda' = ray_sphere ray sphere in
	if lambda' >= lambda then hit else
	  lambda', unitise (ray.orig +| lambda' *| ray.dir -| sphere.center)
    | Group (sphere, scenes) ->
	if ray_sphere ray sphere >= lambda then hit else
	  List.fold_left aux hit scenes in
  aux (infinity, {x=0.; y=0.; z=0.}) scene
let rec ray_trace light ray scene =
  let lambda, normal = intersect ray scene in
  if lambda = infinity then 0. else
    let g = dot normal light in
    if g >= 0. then 0. else
      let p = ray.orig +| lambda *| ray.dir +| delta *| normal in
      (fun (lambda, _) -> if lambda < infinity then 0. else -.g)
	(intersect { orig = p; dir = -1. *| light } scene)
let rec create n c r =
  let obj = Sphere { center = c; radius = r } in
  if n = 1 then obj else
    let r' = 3. *. r /. sqrt 12. in
    let aux l (x, z) = create (n-1) (c +| {x=x; y=r'; z=z}) (r /. 2.) :: l in
    Group ({ center = c; radius = 3. *. r }, List.fold_left aux [obj]
	     [-.r', -.r'; r', -.r'; -.r', r'; r', r'])
let () =
  let level = match Sys.argv with [| _; l |] -> int_of_string l | _ -> 6 in
  let scene = create level { x = 0.; y = -1.; z = 0. } 1. in
  let n = 512 and ss = 4 and light = unitise {x= -1.; y= -3.; z=2.} in
  Printf.printf "P5\n%d %d\n255\n" n n;
  for y = n - 1 downto 0 do
    for x = 0 to n - 1 do
      let g = ref 0. in
      for dx = 0 to ss - 1 do
	for dy = 0 to ss - 1 do
	  let aux x d = float (x - n / 2) +. float d /. float ss in
	  let d = unitise { x = aux x dx; y = aux y dy; z = float n } in
	  g := !g +. ray_trace light {orig={x=0.; y=0.; z= -4.}; dir=d} scene
	done
      done;
      let g = 0.5 +. 255. *. !g /. float (ss*ss) in
      Printf.printf "%c" (char_of_int (int_of_float g));
    done
  done