module tracetest;
import std.math, std.string: atoi;
import std.c.stdio: printf;
import std.gc: malloc;

const double delta; static this() { delta=sqrt(double.epsilon); }

struct Vec {
  double x, y, z;
  Vec opAdd(ref Vec b) { return Vec(x+b.x, y+b.y, z+b.z); }
  Vec opSub(ref Vec b) { return Vec(x-b.x, y-b.y, z-b.z); }
  Vec opMul(double a) { return Vec(x*a, y*a, z*a); }
  double dot(ref Vec b) { return x*b.x + y*b.y + z*b.z; }
  Vec unitise() { return opMul(1f / sqrt(dot(*this))); }
  static Vec opCall(double a, double b, double c) { Vec res=void; res.x=a; res.y=b; res.z=c; return res; }
}

struct Hit {
  double dist;
  Vec pos;
}
struct Ray { Vec orig, dir; }

class Scene {
  abstract void intersect(ref Hit hit, ref Ray ray);
}

template StdConstructor() {
  this(typeof(this.tupleof) p) {
    foreach (index, bogus; this.tupleof) this.tupleof[index]=p[index];
  }
}

final class Sphere : Scene {
  Vec center;
  double radius;
  mixin StdConstructor;
  double ray_sphere(ref Ray ray) {
    auto v = center - ray.orig;
    double b = v.dot(ray.dir), disc = b*b - v.dot(v) + radius * radius;
    if (disc < 0) return double.infinity;
    double d = sqrt(disc), t2 = b + d;
    if (t2 < 0) return double.infinity;
    auto t1 = b - d;
    return t1 > 0 ? t1 : t2;
  }

  void intersect(ref Hit hit, ref Ray ray) {
    auto lambda = ray_sphere(ray);
    if (lambda >= hit.dist) return;
    hit = Hit(lambda, (ray.orig + lambda*ray.dir - center).unitise);
  }
}

final class Group : Scene {
  Sphere bound;
  Scene[] child;
  mixin StdConstructor;
  void intersect(ref Hit hit, ref Ray ray) {
    double l = bound.ray_sphere(ray);
    if (l < hit.dist)
      for (int i=0; i<child.length; ++i)
        child[i].intersect(hit, ray);
  }
}

Hit intersect(ref Ray ray, Scene s) {
  auto res = Hit(double.infinity, Vec(0, 0, 0));
  s.intersect(res, ray); return res;
}

double ray_trace(ref Vec light, ref Ray ray, Scene s) {
  auto hit = intersect(ray, s);
  if (hit.dist == double.infinity) return 0;
  auto g = hit.pos.dot(light);
  if (g >= 0) return 0.;
  auto p = ray.orig + ray.dir*hit.dist + hit.pos*delta;
  return (intersect(Ray(p, -1. * light), s).dist < double.infinity ? 0 : -g);
}

Scene create(int level, ref Vec c, double r) {
  auto s = new Sphere(c, r);
  if (level == 1) return s;
  Scene[] child; child~=s; //child.length=1; child[0]=s;
  auto rn = 3f*r/sqrt(12f);
  for (auto dz=-1; dz<=1; dz+=2)
    for (auto dx=-1; dx<=1; dx+=2)
      child~=create(level-1, c + rn*Vec(dx, 1, dz), r/2f);
  return new Group(new Sphere(c, 3*r), child);
}

int main(/*string[] args*/) {
  int level = 6, n = 512, ss = 4;
  //if (args.length == 2) level = atoi(args[1]);
  auto light = Vec(-1, -3, 2).unitise;
  auto s=create(level, Vec(0, -1, 0), 1);
  printf("P5\n%i %i\n255\n", n, n);
  for (int y=n-1; y>=0; --y)
    for (int x=0; x<n; ++x) {
      double g=0;
      for (int dx=0; dx<ss; ++dx)
        for (int dy=0; dy<ss; ++dy) {
          auto dir=Vec(x+dx*1./ss-n/2., y+dy*1./ss-n/2., n).unitise;
          g += ray_trace(light, Ray(Vec(0, 0, -4), dir), s);
        }
      if (g!=g) g=0f;
      printf("%c", cast(char)(.5f + 255f * g / (ss*ss)));
    }
  return 0;
}