;;;
;;; Mandelbrot set and Julia set
;;; 
;;; 2018.2.5, Tomio Arisaka

;; The Mandelbrot set is made by iteration of the complex plane
;; 
;; c = x + y * %i
;; z[0] = 0
;; z[n+1] = z[n]^2 + c
;; abs(z) is the radius of a circle on the complex plane used as a boundary
;; to determine when iteration can stop. (abs(z) > 2.0)
;; z[1]   = c
;; z[1]^2 = c^2
;;        = (x + y*%i)^2
;;        = x^2 + (-y^2) + 2*x*y*%i
;;        = zx + zy*%i
;; zy[n+1] = 2*zx[n]*zy[n] + y
;; zx[n+1] = zx[n]^2 - zy[n]^2 + x
;; abs(z) = (zx^2 + zy^2)^(1/2)
;;        = pr^(1/2)
;;
;; calculate the number of iteration of the complex plane
;;
;;   x           : the real part of the complex number
;;   y           : the imaginary part of the complex number
;;   iteration   : the maximum iteration of the complex plane
;; 
(defun mandelbrot-e1 (x y iteration)
  (let ((n 0) (zx 0.0) (zy 0.0) (pr 0.0) (px 0.0) (py 0.0))
    (loop while (and (<= pr 4.0) (< n iteration))
          do (progn
               (setq zy (+ (* zx zy 2) y)
                     zx (+ (- px py) x)
                     px (* zx zx)
                     py (* zy zy)
                     pr (+ px py))
               (incf n)))
    n))

;; enter the calculated Mandelbrot set into the matrix 
;;
;;   mat         : the value of the matrix variable of Maxima
;;                 in which the calculated Mandelbrot set is set
;;   iteration   : the maximum iteration of the complex plane
;;   x1,x2       : the range of the real part of the complex plane  (x1 < x2)
;;   y1,y2       : the range of the imaginary part of the complex plane  (y1 < y2)
;; 
(defun mandelbrot-set1 (mat iteration x1 x2 y1 y2)
  (let ((i 0) (j 0) (dx 0.0) (dy 0.0)
        (fx1 (float x1)) (fx2 (float x2)) (fy1 (float y1)) (fy2 (float y2)))
    (cond
     (($matrixp mat)
      (setq dx (/ (- fx2 fx1) (1- (length (cdadr mat)))))
      (setq dy (/ (- fy2 fy1) (1- (length (cdr mat)))))
      (setq j 1)
      (loop for y downfrom fy2 downto fy1 by dy
            do (progn
                 (setq i 1)
                 (loop for x from fx1 to fx2 by dx
                       do (progn
                            (setf (nth i (nth j mat)) (mandelbrot-e1 x y iteration))
                            (incf i)))
                 (incf j)))
      t)
     (t nil))))


;; The Julia set is made by iteration of the complex plane
;; 
;; c is a complex number. z is a variable of the complex plane.
;; c = x0 + y0 * %i
;; z = x + y * %i
;; f(z) = z^2 + c
;; abs(f(z)) is the radius of a circle on the complex plane used as a boundary
;; to determine when iteration can stop. (abs(f(z)) > 2.0)
;; f(z) = z^2 + c
;;      = (x + y*%i)^2 + (x0 + y0*%i)
;;      = x^2 + (-y^2) + x0 + (2*x*y + y0)*%i
;;      = zx + zy*%i
;; zy = 2*x*y + y0
;; zx = x^2 - y^2 + x0
;; abs(f(z)) = (zx^2 + zy^2)^(1/2)
;;           = pr^(1/2)
;;
;; calculate the number of iteration of the complex plane
;;
;;   x0          : the real part of a complex number c
;;   y0          : the imaginary part of a complex number c
;;   x           : the real part of the complex number
;;   y           : the imaginary part of the complex number
;;   iteration   : the maximum iteration of the complex plane
;; 
(defun julia-e1 (x0 y0 x y iteration)
  (let* ((zx (float x)) (px (* zx zx)) (zy (float y)) (py (* zy zy)) (pr 0.0) (n 0))
    (loop while (and (<= pr 4.0) (< n iteration))
          do (progn
               (setq zy (+ (* zx zy 2) y0)
                     zx (+ (- px py) x0)
                     px (* zx zx)
                     py (* zy zy)
                     pr (+ px py))
               (incf n)))
    n))

;; enter the calculated Julia set into the matrix 
;;
;;   mat         : the value of the matrix variable of Maxima
;;                 in which the calculated Julia set is set
;;   iteration   : the maximum iteration of the complex plane
;;   x0          : the real part of a complex number c
;;   y0          : the imaginary part of a complex number c
;;   x1,x2       : the range of the real part of the complex plane  (x1 < x2)
;;   y1,y2       : the range of the imaginary part of the complex plane  (y1 < y2)
;; 
(defun julia-set1 (mat iteration x0 y0 x1 x2 y1 y2)
  (let ((i 0) (j 0) (dx 0.0) (dy 0.0) (fx0 (float x0)) (fy0 (float y0))
        (fx1 (float x1)) (fx2 (float x2)) (fy1 (float y1)) (fy2 (float y2)))
    (cond
     (($matrixp mat)
      (setq dx (/ (- fx2 fx1) (1- (length (cdadr mat)))))
      (setq dy (/ (- fy2 fy1) (1- (length (cdr mat)))))
      (setq j 1)
      (loop for y downfrom fy2 downto fy1 by dy
            do (progn
                 (setq i 1)
                 (loop for x from fx1 to fx2 by dx
                       do (progn
                            (setf (nth i (nth j mat)) (julia-e1 fx0 fy0 x y iteration))
                            (incf i)))
                 (incf j)))
      t)
     (t nil))))