val delta = Math.sqrt (Real.nextAfter(1.0, 2.0) - 1.0)
val infinity = Real.posInf
type vec = real * real * real
infix 7 *| fun s *| (x, y, z) : vec = (s*x, s*y, s*z)
infix 6 +| fun (x1, y1, z1) +| (x2, y2, z2) : vec = (x1+x2, y1+y2, z1+z2)
infix 6 -| fun (x1, y1, z1) -| (x2, y2, z2) : vec = (x1-x2, y1-y2, z1-z2)
fun dot (x1, y1, z1) (x2, y2, z2) : real = x1*x2 + y1*y2 + z1*z2
fun unitise r = (1.0 / Real.Math.sqrt (dot r r)) *| r
type ray = { orig: vec, dir: vec }
type sphere = { center: vec, radius: real }
datatype scene = Sphere of sphere | Group of sphere * scene list
fun ray_sphere ray sphere =
    let val v = #center sphere -| #orig ray
	val b = dot v (#dir ray)
	val disc = b*b - dot v v + #radius sphere * #radius sphere in
	if disc < 0.0 then infinity else
	let val disc = Real.Math.sqrt disc
	    val t2 = b + disc in
	    if t2 < 0.0 then infinity else
	    (fn t1 => if t1 > 0.0 then t1 else t2) (b - disc)
	end
    end
fun intersect ray scene =
    let fun aux (Sphere sphere, (l, n)) =
	    let val l' = ray_sphere ray sphere in
		if l' >= l then (l, n) else
		(l', unitise (#orig ray +| l' *| #dir ray -| #center sphere))
	    end
	  | aux (Group (sphere, scenes), (hit as (l, _))) =
	    if ray_sphere ray sphere >= l then hit else foldl aux hit scenes
    in aux (scene, (infinity, (0.0, 0.0, 0.0))) end
fun ray_trace light ray scene =
    let val (lambda, normal) = intersect ray scene in
	if lambda >= infinity then 0.0 else
	let val g = dot normal light in
	    if g >= 0.0 then 0.0 else
	    let val p = #orig ray +| lambda *| #dir ray +| delta *| normal
		val (l, _) = intersect { orig=p, dir= ~1.0 *| light } scene in
		if l >= infinity then ~g else 0.0
	    end
	end
    end
fun create level r (c as (x, y, z)) =
    let val obj = Sphere { center = c, radius = r } in
	if level = 1 then obj else
	let val r' = 3.0 * r / Real.Math.sqrt 12.0
	    fun aux x' z' = create (level-1) (0.5 * r) (x-x', y+r', z+z')
	    val objs = [obj, aux (~r') (~r'), aux r' (~r'),
			aux (~r') r', aux r' r'] in
	    Group ({ center = c, radius = 3.0 * r }, objs)
	end
    end
val n = 512
val ss = 4
val level = case CommandLine.arguments () of
		[s] => (case Int.fromString s of
			    SOME n => n | _ => 6)
	      | _ => 6
val scene = create level 1.0 (0.0, ~1.0, 0.0)
val light = unitise (~1.0, ~3.0, 2.0)
fun loop f accu (x, y, n) =
    if y=n then accu else
    if x=n then loop f accu (0, y+1, n) else
    loop f (f accu (real x) (real y)) (x+1, y, n)
fun eye_ray n ss x y g dx dy =
    g + ray_trace light { orig = (0.0, 0.0, ~4.0),
			  dir = unitise (x + dx/ss, y + dy/ss, n) } scene
fun pixel n ss () x y =
    let val x = x - n / 2.0
	val y = (n - 1.0) / 2.0 - y
	val g = loop (eye_ray n (real ss) x y) 0.0 (0, 0, ss) in
	print (String.str(Char.chr(Real.trunc (255.0 * g / real (ss*ss)))))
    end
val () =
    (fn s => print ("P5\n"^s^" "^s^"\n255\n")) (Int.toString n);
    loop (pixel (real n) ss) () (0, 0, n);