let zero = 0., 0., 0. let ( *| ) s (rx, ry, rz) = s *. rx, s *. ry, s *. rz let ( +| ) (ax, ay, az) (bx, by, bz) = ax +. bx, ay +. by, az +. bz let ( -| ) (ax, ay, az) (bx, by, bz) = ax -. bx, ay -. by, az -. bz let dot (ax, ay, az) (bx, by, bz) = ax *. bx +. ay *. by +. az *. bz let unitise r = 1. /. sqrt (dot r r) *| r let rec intersect o d (l, _ as hit) (c, r, s) = let v = c -| o in let b = dot v d in let disc = sqrt(b *. b -. dot v v +. r *. r) in let t1 = b -. disc and t2 = b +. disc in let l' = if t2>0. then if t1>0. then t1 else t2 else infinity in if l' >= l then hit else match s with [] -> l', unitise (o +| l' *| d -| c) | ss -> List.fold_left (intersect o d) hit ss let light = unitise (1., 3., -2.) and ss = 4 let rec create level c r = let obj = c, r, [] in if level = 1 then obj else let a = 3. *. r /. sqrt 12. in let aux x' z' = create (level - 1) (c +| (x', a, z')) (0.5 *. r) in c, 3.*.r, [obj; aux (-.a) (-.a); aux a (-.a); aux (-.a) a; aux a a] let level, n = try int_of_string Sys.argv.(1), int_of_string Sys.argv.(2) with _ -> 9, 512 let scene = create level (0., -1., 4.) 1. let rec ray_trace dir = let l, n = intersect zero dir (infinity, zero) scene in let g = dot n light in if g <= 0. then 0. else let p = l *| dir +| sqrt epsilon_float *| n in if fst(intersect p light (infinity, zero) scene) < infinity then 0. else g let aux x d = float x -. float n /. 2. +. float d /. float ss;; 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 d = 0 to ss*ss - 1 do g := !g +. ray_trace(unitise(aux x (d mod ss), aux y (d/ss), float n)) done; let g = 0.5 +. 255. *. !g /. float(ss*ss) in print_char(char_of_int(int_of_float g)) done; done