(asdf:operate 'asdf:load-op '#:png)
(use-package :png)

(defun make-pixel (red green blue alpha)
  (let ((pixel 0))
    (setf pixel (dpb red   (byte 8  0) pixel))
    (setf pixel (dpb green (byte 8  8) pixel))
    (setf pixel (dpb blue  (byte 8 16) pixel))
    (setf pixel (dpb alpha (byte 8 24) pixel))))

#|
3D mandelbrot uncovers the 4D version:
the imaginary plane arises from the square root of negative 
unit from the real axis
generalizing to a third axis: the imaginary volume arises from
the square root of the negative unit from the imaginary axis.
ie:  plane is  x + root(-1)y       is  x + iy
and volume is  x + iy + root(-i)z  
             = x + iy + i.root(i).z
             = x + iy + jz       
                 (where j is i.root(i) is root(-i))
squaring a point in the volume:
  (x + iy + jz)^2
= x^2 - y^2 + j^2.z^2 + 2ixy + 2jxz + 2ijyz
= x^2 - y^2 - iz^2    + 2ixy + 2jxz + 2i(i.root(i))z
= x^2 - y^2 - iz^2    + 2ixy + 2jxz - 2kz (where k is root(i))

this reveals a fourth axis: root(i).  j and k are 
orthogonal to each other and to the other axes.

so attempting to generalise from a fractal embedded in a plane
to a volume uncovers something 4D.  i'm sure this must have been
done before. so i use kt (for time) and jz (for third spatial 
dimension)
 
squaring and adding in 2D:
(x+iy)^2 +x0+iy0
x^2 + 2ixy - y^2 +x0+iy0
X=x^2-y^2+x0
Y=2xy+y0

squaring and adding in 4D:   ij=i.irooti
(x+iy+jz+kt)^2 + x0+iy0+jz0+kt0
= x^2-y^2-iz+it
 +2ixy+2jxz+2kxt- 2iyz+2jyt- 2zt   + x0+iy0+jz0+kt0
= x^2-y^2+x0     (1)
 +2ixy-2iyz+iy0  (i=root(-1))
 +2jxz+2jyt+jz0  (j=root(-i))
 +2kxt+kt0       (k=root(i ))

now, i want to create a movie of  the mandelblob by starting z=0
and varying t, checking magnitude of four-vector:

(defun mb (x y max disq)
  (do ((i 0 (+ i 1))
       (xn 0 (+ (- (expt xn 2) (expt yn 2)) x))
       (yn 0 (+ (* 2 xn yn) y)))
      ((or (> (+ (expt xn 2) (expt yn 2)) disq)
	   (>= i max))
       i)))
|#
(defun mb4 (x y z tm max disq)
  (do ((i 0 (+ i 1))
       (xn 0 (+ (- (expt xn 2) (expt yn 2)) x))
       (yn 0 (+ (- (* 2 xn yn) (* 2 yn zn)) y))
       (zn 0 (+ (* 2 xn zn) (* 2 yn tn) z))
       (tn 0 (+ (* 2 xn tn) tm)))
      ((or (> (+ (expt xn 2) (expt yn 2) (expt zn 2) (expt tn 2)) disq)
	   (>= i max))
       i)))
#| <----- z -->   z=100  scale=0.01
   +------|---+   x,y ~= 60,30 ~= 0.1,0.2
   |          |
   |          |   x,-y ~= 60,-30
   - (x,y).   -   0.1=(px-z/2)*scale-xo
   |          |   0.2=(z/2-py))*scale-yo
   |          |
   |          |
   +------|---+
|#
(defun mod-a (a m) (aref a (mod m (array-dimension a 0))))
(defun mbpix (xo yo zo to xz yz scale max disq imap)
  (let ((image (make-array (list yz xz))))
    (dotimes (xi xz)
      (dotimes (yi yz)
	(let* ((x (- (* (- xi (/ xz 2)) scale) xo))
	       (y (- (* (- (/ yz 2) yi) scale) yo))
               (m (mb4 x y zo to max disq)) )
          ;;;(format t "~S ~S ~S~%" x y m)
	  (setf (aref image yi xi) (mod-a imap m)))))
    image))
(defun div (a b) (values (floor a b)))
(defun rat (x y z) (div (* x y) z))
(defun mk-imap (z)
  (do ((imap (make-array z))
       (i 0 (+ i 1)))
      ((>= i z) imap)
    (let ((x (rat 255 i z)))
      (setf (aref imap i) (make-pixel 0 0 x 0)))))
(defvar *png* "/home/jack/mb.png")
(defvar *imap* (mk-imap 16))
(format t "~S~%" *imap*)
(defvar *zo* 0)
(defvar *to* 0.03)
(dotimes (n 99)
 (let* ((fm (format nil "/home/jack/out/4Dmb~2,'0D" n))
       (fn (format nil "~A0.png" fm))
       (ti (* *to* n)))
  (png:encode-file (mbpix 0.6 0 *zo* ti 150 100 0.02 30 900 *imap*) fn)
  (do ((nn 1 (+ nn 1)))
      ((> nn 9))
    (run-program "/bin/ln" 
		 (list "-s" fn (format nil "~A~A.png" fm nn))))))