(require 'asdf)

(push #P"/usr/local/asdf-install/site-systems/" asdf:*central-registry*)

(asdf:operate 'asdf:load-op 'zpng :verbose nil)

(defun box-muller (&optional (random-state *random-state*))
  (let ((u1 (random 1.0 random-state))
	(u2 (random 1.0 random-state)))
    (let ((rr (sqrt (* -2.0 (log u1))))
	  (aa (* 2.0 pi u2)))
      (values (* rr (cos aa))
	      (* rr (sin aa))))))

(defstruct point
  location
  velocity)

(defun make-random-vector ( dims &optional (scale 1.0d0)
			                   (random-state *random-state*))
  (flet ((make-random-list ()
	   (subseq (loop :for ii :from 0 :below dims :by 2
		      :appending (multiple-value-bind (aa bb)
				     (box-muller random-state)
				   (list (* aa scale)
					 (* bb scale))))
		   0
		   dims)))
    (make-array (list dims) :element-type 'double-float
		:initial-contents (make-random-list))))

(defun make-random-point ( &optional (random-state *random-state*) )
  (make-point :location (make-random-vector 5 2.0 random-state)
	      :velocity (make-random-vector 5 0.25 random-state)))

#+nope
(defun make-random-point ( &optional (random-state *random-state*) )
  (make-point :location (make-array 5
				    :initial-contents (list 1.0
							    0.2
							    0.0
							    0.0
							    0.0))
	      :velocity (make-random-vector 5 0.25 random-state)))

(defun point-amplitude ( point )
  (* (aref (point-location point) 0) 0.5))
(defun point-x ( point )
  (aref (point-location point) 1))
(defun point-y ( point )
  (aref (point-location point) 2))
(defun point-h ( point )
  (aref (point-location point) 3))
(defun point-k ( point )
  (aref (point-location point) 4))

(defun make-rgb-value (vv)
  (min (max (round (* vv 255)) 0) 255))

(defun write-image ( pathname buffer )
  (let ((png (make-instance 'zpng:png
			    :color-type (case (array-dimension buffer 2)
					  (1 :grayscale)
					  (2 :grayscale-alpha)
					  (3 :truecolor)
					  (4 :truecolor-alpha))
			    :width (array-dimension buffer 1)
			    :height (array-dimension buffer 0))))
    (loop :for jj :from 0 :below (zpng:height png)
       :do (loop :for ii :from 0 :below (zpng:width png)
	      :do (loop :for kk :from 0 :below (array-dimension buffer 2)
		     :do (setf (aref (zpng:data-array png) jj ii kk)
			       (make-rgb-value (aref buffer jj ii kk))))))
    (zpng:write-png png pathname)))

(defun hsv-to-rgb (hh ss vv)
  (let ((angle (* (mod hh (* 2.0 pi)) 180 (/ pi))))
    (let ((ff (- (/ angle 60) (floor (/ angle 60)))))
      (let ((pp (* vv (- 1 ss)))
	    (qq (* vv (- 1 (* ff ss))))
	    (tt (* vv (- 1 (* (- 1 ff) ss)))))
	(case (mod (floor (/ angle 60)) 6)
	  (0 (values vv tt pp 1.0))
	  (1 (values qq vv pp 1.0))
	  (2 (values pp vv tt 1.0))
	  (3 (values pp qq vv 1.0))
	  (4 (values tt pp vv 1.0))
	  (5 (values vv pp qq 1.0)))))))

(defun normalize (nx ny nz)
  (let ((scale (sqrt (+ (* nx nx) (* ny ny) (* nz nz)))))
    (if (< scale 0.00001)
	(values 0.0 0.0 0.0)
	(values (/ nx scale) (/ ny scale) (/ nz scale)))))

(defun make-normal ( fx fy)
  (normalize fx fy -7.0))

(defun calculate-color ( points lights xx yy )
  (labels ((calc ( points &optional (f 0.0d0) (fx 0.0d0) (fy 0.0d0) )
	     (cond
	       ((null points) (values f fx fy))
	       (t (let ((pp (first points)))
		    (let ((aa (point-amplitude pp))
			  (ax (* 2.0 pi (+ (* xx (point-x pp)) (point-h pp))))
			  (ay (* 2.0 pi (+ (* yy (point-y pp)) (point-k pp)))))
		      (let ((cx (cos ax))
			    (cy (cos ay))
			    (sx (sin ax))
			    (sy (sin ay)))
			(calc (rest points)
			      (+ f (* aa cx cy))
			      (+ fx (* aa 2.0 pi (point-x pp) sx cy))
			      (+ fy (* aa 2.0 pi (point-y pp) cx sy))))))))))
    (multiple-value-bind (f fx fy)
	(calc points)
      (flet ((specular (nx ny nz)
	       (loop :for ll :in lights
		  :summing (let ((dot (+ (* nx (first ll))
					 (* ny (second ll))
					 (* nz (third ll)))))
			     (if (minusp dot)
				 (expt (- dot) (fourth ll))
				 0.0)))))
	(multiple-value-bind (nx ny nz)
	    (make-normal fx fy)
	  (multiple-value-bind (dr dg db)
	      (hsv-to-rgb (* (abs f) 3.0) 0.9 (abs nz))
	    (let ((ss (* (specular nx ny nz) 0.3)))
	      (values (+ dr ss)
		      (+ dg ss)
		      (+ db ss)))))))))

(defun make-image ( pathname points lights &optional (width 512) (height 512) )
  (let ((buffer (make-array (list height width 3))))
    (loop :for jj :from 0 :below height
       :do (loop :for ii :from 0 :below width
	      :do (multiple-value-bind (rr gg bb)
		      (calculate-color points
				       lights
				       (/ (- ii (* width 0.5)) width 0.75)
				       (/ (- jj (* height 0.5)) width 0.75))
		    (setf (aref buffer jj ii 0) rr
			  (aref buffer jj ii 1) gg
			  (aref buffer jj ii 2) bb))))
    (write-image pathname buffer)))

(defvar +max-lead-acceleration+ 8.0)
(defvar +max-lead-velocity+ 12.0)
(defvar +max-acceleration+ 15.0)
(defvar +max-velocity+ 25.0)

(defun advance-points ( points elapsed )
  (labels ((subtract-with-copy (aa bb)
	     (let ((ret (copy-seq aa)))
	       (loop :for ii :from 0 :below (length aa)
		  :do (decf (aref ret ii) (aref bb ii)))
	       ret))
	   (dot (aa bb)
	     (loop :for ii :from 0 :below (length aa)
		:summing (* (aref aa ii) (aref bb ii))))
	   (mag (aa)
	     (sqrt (dot aa aa)))
	   (scale (aa ss)
	     (loop :for ii :from 0 :below (length aa)
		:do (setf (aref aa ii) (* (aref aa ii) ss)))
	     aa)
	   (increment (aa bb ss)
	     (loop :for ii :from 0 :below (length aa)
		:do (incf (aref aa ii) (* (aref bb ii) ss)))
	     aa)
	   (advance (pp target &optional (max-acc +max-acceleration+)
			          (max-vel +max-velocity+))
	     (let ((aa (subtract-with-copy (point-location pp) target)))
	       (scale aa (* (- (mag aa)) max-acc))
	       (increment (point-velocity pp) aa elapsed)
	       (when (> (mag (point-velocity pp)) max-vel)
		 (scale (point-velocity pp)
			(/ max-vel (mag (point-velocity pp)))))
	       (increment (point-location pp) (point-velocity pp) elapsed))))
    ;; the first point is the leader's goal
    (advance (second points)
	     (point-location (first points))
	     +max-lead-acceleration+ +max-lead-velocity+)
    ;; when the leader gets too close to the goal, give him a new goal
    (when (< (mag (subtract-with-copy (point-location (first points))
				      (point-location (second points))))
	     3.0)
      (setf (first points) (make-random-point))
      (format t "~S~%" (point-location (first points))))

    ;; everyone else follows the leader
    (mapc #'(lambda (pp)
	      (advance pp (point-location (second points))))
	  (rest (rest points)))
    points))
	       

(ignore-errors
  (with-open-file (in #P"random-seed.lisp"
		      :direction :input)
    (let ((rs (read in nil nil)))
      (when rs
	(setf *random-state* rs)))))

(loop :for ii :from 1 :to 3600
   :with elapsed = (/ 1800.0)
   :for points = (let ((orig (loop :for ii :from 1 :to 14
				:collecting (make-random-point))))
		   (loop :for ii :from 1 :to 900
		      :do (advance-points orig elapsed))
		   orig)
        :then (advance-points points elapsed)
   :with lights = (list (list 0.3 0.5 (sqrt 0.66) 35.0)
			(list -0.3 (- (sqrt 0.66)) 0.5 45.0)
			(list (sqrt 0.66) 0.3 0.5 15.0))
   :do (make-image (format nil "foo~4,'0D.png" ii)
		   (rest (rest points))
		   lights
		   640 360))

(with-open-file (out #P"random-seed.lisp"
		     :direction :output
		     :if-exists :supersede
		     :if-does-not-exist :create)
  (write *random-state* :stream out))