// 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


[Record] \
class Ray
  public orig : Vec
  public dir : Vec
  
[Record] \
class Hit
  public mutable lambda : double
  public mutable dir : Vec
  
abstract class Scene
  public abstract intersect (r : Ray, h : Hit) : void
  
  [Record] \
  internal sealed class Sphere : Scene
    public center : Vec
    public radius : double
    public override intersect (r : Ray, h : Hit) : void
      def l = M.ray_sphere (r, center, radius)
      when (l < h.lambda)
        h.lambda = l
        h.dir = center
        h.dir = (r.orig + l * r.dir - center).Unitise ()
  
  [Record] \
  internal sealed class Group : Scene
    public center : Vec
    public radius : double
    public scenes : list [Scene]
    
    public override intersect (r : Ray, h : Hit) : void
        def l = M.ray_sphere (r, center, radius)
        when (l < h.lambda)
          foreach (s in scenes)
            s.intersect (r, h)
  
module M
  delta = 1.49011611938476562e-08
  inf : double = 1.0 / 0.0

  internal ray_sphere (r : Ray, center : Vec, radius : double) : double
    def v = center - r.orig
    def b = v ** r.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

  ray_trace (light:Vec, orig:Vec, dir:Vec, scene:Scene) : double
    def h = Hit (inf, Vec())
    scene.intersect (Ray (orig, dir), h)
    if (h.lambda == inf) 0.0
    else
      h.dir = (orig + h.lambda * dir - h.dir).Unitise ()
      def g = h.dir ** light
      if (g >= 0) 0.0
      else
        def p = orig + h.lambda * dir + delta * h.dir
        h.lambda = inf
        scene.intersect (Ray (p, -1 * light), h)
        if (h.lambda < inf) 0.0
        else -g

  internal main (level : int, n : int) : void
    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 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] => M.main (int.Parse (level), int.Parse (n))
  | _ => M.main (9, 512)