// based on http://www.ffconsultancy.com/free/ray_tracer/code/1/ray.ml
#pragma indent
using System.Math

[Record] \
struct Vec
  x : double
  y : double
  z : double

  public static @* (s : double, r : Vec) : Vec
    Vec (s * r.x, s * r.y, s * r.z)
  public static @+ (s : Vec, r : Vec) : Vec
    Vec (s.x + r.x, s.y + r.y, s.z + r.z)
  public static @- (s : Vec, r : Vec) : Vec
    Vec (s.x - r.x, s.y - r.y, s.z - r.z)
  public static @** (s : Vec, r : Vec) : double
    s.x * r.x + s.y * r.y + s.z * r.z

  public Unitise () : Vec
    (1 / Sqrt (this ** this)) * this


variant Scene
  | Sphere { m : Vec; r : double; }
  | Group { m : Vec; r : double; l : list [Scene]; }

def delta = 1.49011611938476562e-08
def inf = double.PositiveInfinity

def intersect (orig, dir : Vec, scene)
  def ray_sphere (center, radius : double)
    def v = center - orig
    def b = v ** dir
    def disc = b * b - (v ** v) + radius * radius
    if (disc < 0) inf
    else
      def disc = Sqrt (disc)
      def t2 = b + disc
      def t1 = b - disc
      if (t2 < 0) inf
      else if (t1 > 0) t1
      else t2

  def aux (elt, first)
    def l = first [0]
    match (elt)
      | Scene.Sphere (center, radius) =>
        def l' = ray_sphere (center, radius)
        if (l' >= l) first
        else (l', (orig + l' * dir - center).Unitise ())
      | Scene.Group (center, radius, scenes) =>
        def l' = ray_sphere (center, radius)
        if (l' >= l) first
        else scenes.FoldLeft (first, aux)
  aux (scene, (inf, Vec ()))


def ray_trace (light, orig, dir, scene)
  def (lambda, normal) = intersect (orig, dir, scene)
  if (lambda == inf) 0.0
  else
    def g = normal ** light
    if (g >= 0) 0.0
    else
      def p = orig + lambda * dir + delta * normal
      if ((intersect (p, -1 * light, scene)) [0] < inf) 0.0
      else -g

def create (level, c, r)
  def obj = Scene.Sphere (c, r)
  if (level == 1) obj
  else
    def a = 3 * r / Sqrt (12)
    def aux (x', z')
      create (level - 1, c + Vec (x', a, z'), 0.5 * r)
    Scene.Group (c, 3 * r , [obj, aux (-a, -a), aux (a, -a), aux (-a, a), aux (a, a)])

def main (level, n)
  def scene = create (level, Vec (0, -1, 0), 1)
  def light = Vec (-1, -3, 2).Unitise ()
  def ss = 4
  System.Console.Write ($ "P5\n$n $n\n255\n")
  def s = System.Console.OpenStandardOutput ()
  for (mutable y = n - 1; y >= 0; y--)
    for (mutable x = 0; x < n; x++)
      mutable g = 0.0
      for (mutable dx = 0; dx < ss; dx++)
        for (mutable dy = 0; dy < ss; dy++)
          def aux (x, d)
            x - n / 2.0 + d / (ss :> double)
          def dir = Vec (aux (x, dx), aux (y, dy), n).Unitise ()
          g += ray_trace (light, Vec (0,0,-4), dir, scene)
      def res = (0.5 + 255 * g / (ss*ss)) :> int
      s.WriteByte (res :> byte)
      
match (Nemerle.Collections.List.FromArray (System.Environment.GetCommandLineArgs()))
  | [_, level, n] => main (int.Parse (level), int.Parse (n))
  | _ => main (9, 512)