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

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

struct Vec {
  double x, y, z;
  Vec* set(ref Vec b) { return set(b.x, b.y, b.z); }
  Vec* set(double a, double b, double c) { x=a; y=b; z=c; return this; }
  Vec* add(ref Vec b) { x+=b.x; y+=b.y; z+=b.z; return this; }
  Vec* sub(ref Vec b) { x-=b.x; y-=b.y; z-=b.z; return this; }
  Vec* mul(double a) { x*=a; y*=a; z*=a; return this; }
  double dot(ref Vec b) { return x*b.x + y*b.y + z*b.z; }
  Vec* unitise() { return mul(1f / sqrt(dot(*this))); }
}

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;
    v.sub(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;
  }

  override void intersect(ref Hit hit, ref Ray ray) {
    auto lambda = ray_sphere(ray);
    if (lambda >= hit.dist) return;
    hit.dist = lambda;
    hit.pos.set(ray.dir).mul(lambda).add(ray.orig).sub(center).unitise();
  }
}

final class Group : Scene {
  Sphere bound;
  Scene[] child;
  mixin StdConstructor;
  override 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;
}

/* modifies ray for temporary use */
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.;
  /* hit is not needed here already, use as a temporary */
  ray.orig.add(*ray.dir.mul(hit.dist)).add(*hit.pos.mul(delta));
  ray.dir.set(light).mul(-1);
  return (intersect(ray, 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, *Vec(dx, 1, dz).mul(rn).add(c), r/2f);
  return new Group(new Sphere(c, 3*r), child);
}

extern(C) uint clock();

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