;;; As in The Tinkertoy Computer

(defparameter *resolution* 500)
(defparameter *iterations* 4000)
(defparameter *warmup* 600)

(defun lyapunov-exponent (rs start &optional (pattern '(0)))
  (let ((k (length pattern))
        (zeroes 0))
    (/ (iter (for i from 0 below (+ *warmup* *iterations*))
             (for r = (elt rs (elt pattern (mod i k))))
             (for x first start then (* r x (- 1 x)))
             (when (>= i *warmup*)
               (for diff = (abs (- r (* 2 r x))))
               (if (> diff 0)
                   (sum (log diff 2) into exponent)
                   (incf zeroes))) 
             (finally (return exponent)))
       (if (= *iterations* zeroes)
           1
           (- *iterations* zeroes)))))

(defun write-image (filename pattern &key (minx 2) (maxx 4) (miny 2) (maxy 4) (start 0.5))
  (with-open-file (s filename :direction :output :if-exists :supersede)
    (format s "P2~%~d ~d~%255~%" *resolution* *resolution*)
    (iter (for j from 0 below *resolution*)
          (for y = (+ maxy (* (- miny maxy) (/ j (1- *resolution*)))))
          (iter (for i from 0 below *resolution*)
                (for x = (+ minx (* (- maxx minx) (/ i (1- *resolution*)))))      
                (for e = (lyapunov-exponent (list x y) start pattern))
                (if (>= e 0)
                    (format s "0 ")
                    (format s "~d " (round (* 255 (expt 2 e))))))
          (terpri s))))

#+nil
(write-image "/tmp/test1.pgm" '(0 1))

#+nil
(write-image "/tmp/test2.pgm" '(1 1 0 1 0))

#+nil
(write-image "/tmp/test3.pgm" '(1 1 1 1 1 1 0 0 0 0 0 0) :minx 2.6 :maxx 3.5 :miny 3.1 :maxy 4)