(declaim (inline vec *v +v -v dot magnitude unitise ray-sphere intersect
		 ray-trace))

(defconstant infinity most-positive-double-float)
(defconstant delta (sqrt double-float-epsilon))

(defstruct (vec (:conc-name nil)
		(:constructor vec (x y z))
		(:type (vector double-float)))
  x y z)

(defvar zero (vec 0d0 0d0 0d0))

(defun *v (s a) (vec (* s (x a)) (* s (y a)) (* s (z a))))
(defun +v (a b) (vec (+ (x a) (x b)) (+ (y a) (y b)) (+ (z a) (z b))))
(defun -v (a b) (vec (- (x a) (x b)) (- (y a) (y b)) (- (z a) (z b))))
(defun dot (a b) (+ (* (x a) (x b)) (* (y a) (y b)) (* (z a) (z b))))
(defun magnitude (a) (sqrt (the (double-float 0d0) (dot a a))))
(defun unitise (a) (*v (/ 1d0 (magnitude a)) a))

(defstruct ray (orig zero) (dir zero))
(defstruct sphere (center zero) (radius 0d0 :type double-float))
(defstruct (group (:include sphere)) (children nil :type list))

(defun ray-sphere (orig dir center radius)
  (let* ((v (-v center orig))
	 (b (dot v dir))
	 (disc (+ (- (* b b) (dot v v)) (* radius radius))))
    (if (< disc 0d0)
        infinity
        (let* ((disc (sqrt disc))
               (t2 (+ b disc))
               (t1 (- b disc)))
          (cond ((< t2 0d0) infinity)
                ((> t1 0d0) t1)
                (t t2))))))

(defun intersect (orig dir scene)
  (labels ((aux (lam normal scene)
             (declare (double-float lam))
             (let* ((center (sphere-center scene))
                    (lamt (ray-sphere orig
                                      dir
                                      center
                                      (sphere-radius scene))))
               (if (>= lamt lam)
                   (values lam normal)
                   (etypecase scene
		     (group
		      (dolist (kid (group-children scene))
			(setf (values lam normal)
			      (aux lam normal kid)))
		      (values lam normal))
		     (sphere
		      (values lamt (unitise
				    (-v (+v orig (*v lamt dir)) center)))))))))
    (aux infinity zero scene)))

(defun ray-trace (neg-light orig dir scene)
  (multiple-value-bind (lam normal)
      (intersect orig dir scene)
    (if (= lam infinity)
        0d0
        (let ((g (dot normal neg-light)))
          (if (< g 0d0)
              0d0
              (let ((p (+v (+v orig (*v lam dir))
                           (*v delta normal))))
		(multiple-value-bind (slam snormal)
				     (intersect p neg-light scene)
				     (if (< slam infinity)
					 0d0
				       g))))))))

(defun create (n c r)
  (let ((obj (make-sphere :center c :radius r)))
    (if (= n 1)
        obj
        (let ((rt (* 3d0 (/ r (sqrt 12d0)))))
          (flet ((kid (x z)
                   (create (1- n)
                           (+v c (vec x rt z))
                           (/ r 2d0)))
                 (bound (c scenes)
                   (loop for s in scenes
                         maximize (+ (magnitude (-v c (sphere-center s)))
                                     (sphere-radius s)))))
            (let ((objs (list obj
                              (kid (- rt) (- rt))
                              (kid (- rt) rt)
                              (kid rt (- rt))
                              (kid rt rt)))
                  (center (vec (x c) (+ (y c) r) (z c))))
              (make-group :center center
                          :radius (max (bound center objs) r)
                          :children objs)))))))

(defun main (level file-name n)
  (declare (fixnum n) (optimize speed (safety 0)))
  (let* ((scene (create level (vec 0d0 -1d0 4d0) 1d0))
         (ss 4d0)
         (neg-light (unitise (vec 1d0 3d0 -2d0)))
         (-n/2 (/ n -2d0))
         (1-n/2 (1- (/ n 2d0))))
    (with-open-file (s file-name
                       :element-type '(unsigned-byte 8)
                       :if-exists :supersede
                       :if-does-not-exist :create
                       :direction :output)
      (loop for c across (format nil "P5~%~A ~A~%255~%" n n) do
            (write-byte (char-code c) s))
      (loop for y from 1-n/2 downto -n/2 do
            (loop for x from -n/2 to 1-n/2
                  for g of-type double-float = 0d0 do
                  (loop for dx from x below (1+ x) by (/ ss) do
                        (loop for dy from y below (1+ y) by (/ ss)
                              for dir = (unitise (vec dx dy (float n 0d0))) do
                              (incf g (ray-trace neg-light zero dir scene))))
                  (write-byte (round (* 255d0 (/ g (* ss ss)))) s))))))