;;;; mset-attractors.lisp -- finding the points of attraction in the mandelbrot set
;;;; 
;;;; Copyright (c) 2005
;;;;
;;;; By David Steuber <david@david-steuber.com>
;;;;
;;;; Permission is hereby granted, free of charge, to any person obtaining
;;;; a copy of this software and associated documentation files (the
;;;; "Software"), to deal in the Software without restriction, including
;;;; without limitation the rights to use, copy, modify, merge, publish,
;;;; distribute, sublicense, and/or sell copies of the Software, and to
;;;; permit persons to whom the Software is furnished to do so, subject to
;;;; the following conditions:
;;;; 
;;;; The above copyright notice and this permission notice shall be
;;;; included in all copies or substantial portions of the Software.
;;;; 
;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
;;;; MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
;;;; LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
;;;; OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
;;;; WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
;;;; 
;;;; ***********************************************************************

(defun make-axis-array (n)
  (let ((aa (make-array n :element-type 'double-float))
        (m (1- n)))
    (dotimes (i n)
      (setf (aref aa i)
            (+ -2.0d0 (* 4.0d0 (/ i m)))))
    aa))

(defun make-attractor-map (resolution max-iterations)
  (declare (type fixnum resolution max-iterations)
           (optimize speed))
  (let* ((r3 (* 3 resolution))
         (a3 (make-axis-array r3))
         (aa (make-axis-array resolution))
         (am (make-array (list resolution resolution) :element-type 'fixnum :initial-element 0))
         (interval (- (aref aa 1) (aref aa 0)))
         (min-val (aref aa 0)))
    (declare (type double-float interval min-val)
             (type fixnum r3))
    (labels ((index-from-float (float)
               (declare (type double-float float))
               (when (<= -2.0d0 float 2.0d0)
                 (round (- float min-val) interval)))
             (escape-fn (cr ci)
               (declare (type double-float cr ci))
               (loop
                  for iterations of-type fixnum below max-iterations
                  for z-imag of-type double-float         = 0.0D0 then (+ (* 2.0d0 z-real z-imag) ci)       
                  for z-real of-type double-float         = 0.0D0 then (+ (- z-real-squared z-imag-squared) cr)
                  for z-real-squared of-type double-float = 0.0D0 then (* z-real z-real)
                  for z-imag-squared of-type double-float = 0.0D0 then (* z-imag z-imag)
                  with points = ()
                  while (< (+ z-real-squared z-imag-squared) 4.0D0) do
                    (when (< 0 iterations)
                      (let ((x (index-from-float z-real))
                            (y (index-from-float z-imag)))
                        (and x y (incf (aref am x y))
                             (push (cons x y) points))))
                  finally (when (< (length points) (1- max-iterations))
                            (mapc (lambda (p)
                                    (decf (aref am (car p) (cdr p))))
                                  points)))))
      (dotimes (y r3)
        (dotimes (x r3)
          (escape-fn (aref a3 x) (aref a3 y))))
      am)))

(defun write-pgm-file-from-map (map filespec)
  ;; http://astronomy.swin.edu.au/~pbourke/dataformats/ppm/
  (let* ((width (array-dimension map 1))
         (height (array-dimension map 0))
         (max-val (loop for row below height
                     maximize (loop for col below width
                                 maximize (aref map row col))))
         (b (expt (coerce max-val 'double-float) #.(/ 1.0d0 255.0d0))))
    (with-open-file (pgm filespec :direction :output :if-exists :supersede)
      (format pgm "P2 ~A ~A ~A~%" width height 255)
      (dotimes (y height)
        (dotimes (x width)
          (format pgm "~A~%" (let ((n (coerce (aref map x y) 'double-float)))
                               (if (zerop n)
                                   0
                                   (round (log n b))))))))))