```--- a/share/dynamics/complex_dynamics.lisp
+++ b/share/dynamics/complex_dynamics.lisp
@@ -20,207 +20,124 @@

(in-package :maxima)

-;; Function hsv2rgb transforms a color specification triplet
-;; (hue,saturation,value), into RGB values. Hue should be an integer (angle
-;; in degrees), while saturation and value should be floating-point numbers
-;; between 0 and 1.
-;; The RGB values returned are six hexadecimal digits: two digits for R,
-;; followed by two digits for G and two digits for B.
-
-(defun hsv2rgb(hue sat val)
-  (let (v rgb h i f u p q r)
-  (setq v (* val 255))
-  (if (= sat 0)
-      (setq rgb (+ v (* 256 (+ v (* 256 v)))))     ; achromatic case
-    (progn
-      (setq h (/ (mod hue 360) 60))                ; 0<= h <6
-      (setq i (truncate h))                        ; largest integer <= h
-      (setq f (- h i))                             ; fractional part of h
-      (setq u (round v))                           ; u, p, q and r are all
-      (setq p (round (* v (- 1 sat))))             ; integers between 0
-      (setq q (round (* v (- 1 (* sat f)))))       ; and 255
-      (setq r (round (* v (- 1 (* sat (- 1 f))))))
-      (cond
-       ((eql i 0) (setq rgb (+ p (* 256 (+ r (* 256 u))))))     ; red
-       ((eql i 1) (setq rgb (+ p (* 256 (+ u (* 256 q))))))     ; yellow
-       ((eql i 2) (setq rgb (+ r (* 256 (+ u (* 256 p))))))     ; green
-       ((eql i 3) (setq rgb (+ u (* 256 (+ q (* 256 p))))))     ; cyan
-       ((eql i 4) (setq rgb (+ u (* 256 (+ p (* 256 r))))))     ; blue
-       ((eql i 5) (setq rgb (+ q (* 256 (+ p (* 256 u)))))))))  ; magenta
-  (format nil "~6,'0x" (round rgb))))
-
-;; Function \$mandelbrot creates an XPM graphic file with the Mandelbrot set
+;; Function \$mandelbrot displays the Mandelbrot set
;; on the region: xmin < x <xmax, ymin < y <ymax.
-;; The input argument m must be a positive integer, which defines the number
-;; of iterations to be made, before considering a point as part of the set.
-
-(defun \$mandelbrot(&rest options)
-(let (a b c d e colorbytes colorcodes colors line num
-      dx dy x xmax xmin y ymax ymin
-      (xcenter 0.0) (ycenter 0.0) (radius 2.0)
-      (m 12) (nx 400) (ny 400) (file "mandelbrot.xpm")
-      (huerange 360) (hue0 -60) (sat 0.76) (val 0.96) (color "000000")
-      (\$numer t) (optname nil)
-      (letters "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ"))
-
-  (dolist (option options)
-    (setf optname (second option))
-    (cond
-     ((eql optname '\$size)
-      (if (> (length option) 3)
-	  (progn (setq nx (third option)) (setq ny (fourth option)))
-	  (progn (setq nx (third option)) (setq ny nx))))
-     ((eql optname '\$levels) (setq m (third option)))
-     ((eql optname '\$huerange) (setq huerange (third option)))
-     ((eql optname '\$hue) (setq hue0 (third option)))
-     ((eql optname '\$saturation) (setq sat (third option)))
-     ((eql optname '\$value) (setq val (third option)))
-     ((eql optname '\$color)
-      (setq color (hsv2rgb (third option) (fourth option) (fifth option))))
-     ((eql optname '\$center)
-      (setq xcenter (coerce-float (third option)))
-      (setq ycenter (coerce-float (fourth option))))
-     ((eql optname '\$filename)
-      (setq file (format nil "~a.xpm"
-                     (print-invert-case (stripdollar (third option))))))))
-
-  (setq xmin (- xcenter radius))
-  (setq xmax (+ xcenter radius))
-  (setq ymin (- ycenter radius))
-  (setq ymax (+ ycenter radius))
-  (setq dx (/ (float (- xmax xmin)) nx))   ; x displacement per pixel
-  (setq dy (/ (float (- ymax ymin)) ny))   ; y displacement per pixel
-
-  ;; setup color tables
-  (if (> m 52) (setq colorbytes 2) (setq colorbytes 1))
-  (setq colors (list color))
-  (setq colorcodes
-	(cond ((eql colorbytes 1) '(".")) ((eql colorbytes 2) '(".."))))
-  (dotimes (i m)
-    (push (hsv2rgb (+ hue0 (* huerange (- 1 (/ i (- m 1))))) sat val) colors)
-    (push (cond
-	   ((eql colorbytes 1) (string (char letters i)))
-	   ((eql colorbytes 2) (format nil "~a~a" (char letters (floor i 50))
-				       (char letters (mod i 50)))))
-	  colorcodes))
-
-  (with-open-file (st file :direction :output :if-exists :supersede)
-  ;; print XPM file header
-    (format st "/* XPM */~%static char *mandelbrot[] = {~%")
-    (format st "\"~a ~a ~a ~a\",~%" nx ny (+ m 1) colorbytes)
-
-  ;; print colors table
-    (dotimes (i (+ m 1))
-      (format st "\"~a c #~a\",~%" (nth i colorcodes) (nth i colors)))
-
-  ;; iteration through all grid points
-    (dotimes (i ny)
-      (setq y (+ ymin (* (- ny (+ i 1)) dy)))
-      (setq line "")
-      (dotimes (j nx)
-	(setq x (+ xmin (* j dx)))
-	(setq a 0) (setq b 0)
-	(setq num m)
-	(dotimes (l m)
-	  (setq c (* a a))
-	  (setq d (* b b))
-	  (setq e (* 2 a b))
-	  (when (> (+ c d) 4) (progn (setq num l) (return)))
-	  (setq a (+ x (- c d)))
-	  (setq b (+ y e)))
-	(setq line (concatenate 'string line (nth num colorcodes))))
-      (if (= i (- ny 1))
-	  (format st "\"~a\"~%" line) (format st "\"~a\",~%" line)))
-    (format st "};~%"))
-  (format t "File ~a was created.~%" file)))
-
-;; Function \$julia creates an XPM graphic file with the Julia set for the
+
+(defun \$mandelbrot (&rest extra-options)
+  (let (file-name output-file (options (copy-tree *plot-options*))
+        a b c d e num  dx dy x xmax xmin y ymax ymin m nx ny)
+    (setf (getf options :type) "plot2d")
+    (unless (getf options :x) (setf (getf options :x) '(-2 2)))
+    (unless (getf options :y) (setf (getf options :y) '(-2 2)))
+    (unless (getf options :iterations) (setf (getf options :iterations) 9))
+    (unless (getf options :xlabel) (setf (getf options :xlabel) "x"))
+    (unless (getf options :ylabel) (setf (getf options :ylabel) "y"))
+    (unless (member :colorbox options) (setf (getf options :colorbox) t))
+    (setf (getf options :palette)
+          '(((mlist) \$gradient \$magenta \$violet \$blue \$cyan \$green \$yellow
+            \$orange \$red \$brown \$black)))
+
+    ;; Parses the options given in the command line
+    (setq options (plot-options-parser extra-options options))
+    (setf file-name (plot-temp-file "maxout.gnuplot"))
+    (unless (getf options :yxratio) (setf (getf options :samexy) t))
+
+    (setq xmin (car (getf options :x)))
+    (setq xmax (cadr (getf options :x)))
+    (setq ymin (car (getf options :y)))
+    (setq ymax (cadr (getf options :y)))
+    (setq m (getf options :iterations))
+    (setq nx (car (getf options :grid)))
+    (setq ny (cadr (getf options :grid)))
+    (setq dx (/ (rationalize (- xmax xmin)) nx)) ; x incr. per pixel
+    (setq dy (/ (rationalize (- ymax ymin)) ny)) ; y incr. per pixel
+
+    (with-open-file
+        (st file-name :direction :output :if-exists :supersede)
+      (setq output-file (gnuplot-print-header st options))
+      (format st "set palette ~a~%"
+              (gnuplot-palette (rest (first (getf options :palette)))))
+      (format st "unset key~%")
+      (format st "plot '-' with image~%")
+      ;; iterates through all grid points
+      (dotimes (i ny)
+        (setq y (+ ymin (/ dy 2) (* i dy)))
+        (dotimes (j nx)
+          (setq x (+ xmin (/ dx 2) (* j dx)))
+          (setq a 0) (setq b 0)
+          (setq num m)
+          (dotimes (l m)
+            (setq c (* a a))
+            (setq d (* b b))
+            (setq e (* 2 a b))
+            (when (> (+ c d) 4) (progn (setq num l) (return)))
+            (setq a (+ (float x) (- c d)))
+            (setq b (+ (float y) e)))
+          (format st "~f ~f ~d~%" x y num)))
+      (format st "e~%"))
+    (gnuplot-process options file-name)
+    (format nil "~a~&~a" file-name output-file)))
+
+;; Function \$julia(x,y) displays the Julia set for the
;; point (x,y) of the complex plane, on the region: xmin < x <xmax,
;; ymin < y <ymax.
-;; The input argument m must be a positive integer, which defines the number
-;; of iterations to be made, before considering a point as part of the set.
-
-(defun \$julia(x y &rest options)
-(let (a b b0 c d e colorbytes colorcodes colors line num
-      dx dy xmax xmin ymax ymin
-      (xcenter 0.0) (ycenter 0.0) (radius 2.0)
-      (m 12) (nx 400) (ny 400) (file "julia.xpm")
-      (huerange 360) (hue0 -60) (sat 0.76) (val 0.96) (color "000000")
-      (\$numer t) (optname nil)
-      (letters "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ"))
-
-  (dolist (option options)
-    (setf optname (second option))
-    (cond
-     ((eql optname '\$size)
-      (if (> (length option) 3)
-	  (progn (setq nx (third option)) (setq ny (fourth option)))
-	  (progn (setq nx (third option)) (setq ny nx))))
-     ((eql optname '\$levels) (setq m (third option)))
-     ((eql optname '\$huerange) (setq huerange (third option)))
-     ((eql optname '\$hue) (setq hue0 (third option)))
-     ((eql optname '\$saturation) (setq sat (third option)))
-     ((eql optname '\$value) (setq val (third option)))
-     ((eql optname '\$color)
-      (setq color (hsv2rgb (third option) (fourth option) (fifth option))))
-     ((eql optname '\$center)
-      (setq xcenter (coerce-float (third option)))
-      (setq ycenter (coerce-float (fourth option))))
-     ((eql optname '\$filename)
-      (setq file (format nil "~a.xpm"
-                     (print-invert-case (stripdollar (third option))))))))
-
-  (setq xmin (- xcenter radius))
-  (setq xmax (+ xcenter radius))
-  (setq ymin (- ycenter radius))
-  (setq ymax (+ ycenter radius))
-  (setq dx (/ (float (- xmax xmin)) nx))   ; x displacement per pixel
-  (setq dy (/ (float (- ymax ymin)) ny))   ; y displacement per pixel
-
-  ;; setup color tables
-  (if (> m 52) (setq colorbytes 2) (setq colorbytes 1))
-  (setq colors (list color))
-  (setq colorcodes
-	(cond ((eql colorbytes 1) '(".")) ((eql colorbytes 2) '(".."))))
-  (dotimes (i m)
-    (push (hsv2rgb (+ hue0 (* huerange (- 1 (/ i (- m 1))))) sat val) colors)
-    (push (cond
-	   ((eql colorbytes 1) (string (char letters i)))
-	   ((eql colorbytes 2) (format nil "~a~a" (char letters (floor i 50))
-				       (char letters (mod i 50)))))
-	  colorcodes))
-
-  (with-open-file (st file :direction :output :if-exists :supersede)
-  ;; print XPM file header
-    (format st "/* XPM */~%static char *mandelbrot[] = {~%")
-    (format st "\"~a ~a ~a ~a\",~%" nx ny (+ m 1) colorbytes)
-
-  ;; print colors table
-    (dotimes (i (+ m 1))
-      (format st "\"~a c #~a\",~%" (nth i colorcodes) (nth i colors)))
-
-  ;; iteration through all grid points
-    (dotimes (i ny)
-      (setq b0 (+ ymin (* (- ny (+ i 1)) dy)))
-      (setq line "")
-      (dotimes (j nx)
-	(setq a (+ xmin (* j dx)))
-	(setq b b0)
-	(setq num m)
-	(dotimes (l m)
-	  (setq c (* a a))
-	  (setq d (* b b))
-	  (setq e (* 2 a b))
-	  (when (> (+ c d) 4) (progn (setq num l) (return)))
-	  (setq a (+ x (- c d)))
-	  (setq b (+ y e)))
-	(setq line (concatenate 'string line (nth num colorcodes))))
-      (if (= i (- ny 1))
-	  (format st "\"~a\"~%" line) (format st "\"~a\",~%" line)))
-    (format st "};~%"))
-  (format t "File ~a was created.~%" file)))
+
+(defun \$julia (x y &rest extra-options)
+  (let (file-name output-file (options (copy-tree *plot-options*))
+        num dx dy xmax xmin ymax ymin a b c d e a0 b0 m nx ny)
+    (setf (getf options :type) "plot2d")
+    (unless (getf options :x) (setf (getf options :x) '(-2 2)))
+    (unless (getf options :y) (setf (getf options :y) '(-2 2)))
+    (unless (getf options :iterations) (setf (getf options :iterations) 9))
+    (unless (getf options :xlabel) (setf (getf options :xlabel) "x"))
+    (unless (getf options :ylabel) (setf (getf options :ylabel) "y"))
+    (unless (member :colorbox options) (setf (getf options :colorbox) t))
+    (setf (getf options :palette)
+          '(((mlist) \$gradient \$magenta \$violet \$blue \$cyan \$green \$yellow
+            \$orange \$red \$brown \$black)))
+
+    ;; Parses the options given in the command line
+    (setq options (plot-options-parser extra-options options))
+    (setf file-name (plot-temp-file "maxout.gnuplot"))
+    (unless (getf options :yxratio) (setf (getf options :samexy) t))
+
+    (setq xmin (car (getf options :x)))
+    (setq xmax (cadr (getf options :x)))
+    (setq ymin (car (getf options :y)))
+    (setq ymax (cadr (getf options :y)))
+    (setq m (getf options :iterations))
+    (setq nx (car (getf options :grid)))
+    (setq ny (cadr (getf options :grid)))
+    (setq dx (/ (rationalize (- xmax xmin)) nx)) ; x incr. per pixel
+    (setq dy (/ (rationalize (- ymax ymin)) ny)) ; y incr. per pixel
+
+    (with-open-file
+        (st file-name :direction :output :if-exists :supersede)
+      (setq output-file (gnuplot-print-header st options))
+      (format st "set palette ~a~%"
+              (gnuplot-palette (rest (first (getf options :palette)))))
+      (format st "unset key~%")
+      (format st "plot '-' with image~%")
+
+      ;; iterates through all grid points
+      (dotimes (i ny)
+        (setq b0 (+ ymin (/ dy 2) (* i dy)))
+        (dotimes (j nx)
+          (setq a0 (+ xmin (/ dx 2) (* j dx)))
+          (setq a (float a0))
+          (setq b (float b0))
+          (setq num m)
+          (dotimes (l m)
+            (setq c (* a a))
+            (setq d (* b b))
+            (setq e (* 2 a b))
+            (when (> (+ c d) 4) (progn (setq num l) (return)))
+            (setq a (+ x (- c d)))
+            (setq b (+ y e)))
+          (format st "~f ~f ~d~%" a0 b0 num)))
+      (format st "e~%"))
+    (gnuplot-process options file-name)
+    (format nil "~a~&~a" file-name output-file)))

;; Function \$rk implements the 4th order Runge-Kutta numerical method.
;;  exprs:   an expression or maxima list with n expressions
@@ -229,8 +146,8 @@
;;  domain:  maxima list with four elements (name of the independent
;;           variable, its initial and final values and its increment)
(defun \$rk (exprs vars initial domain
-             &aux d u fun k1 k2 k3 k4 r1 r2 r3 traj r
-             (it (mapcar #'coerce-float (cddr domain))))
+            &aux d u fun k1 k2 k3 k4 r1 r2 r3 traj r
+              (it (mapcar #'coerce-float (cddr domain))))
(unless (\$listp exprs) (setq exprs `((mlist) ,exprs)))
(unless (\$listp initial) (setq initial `((mlist) ,initial)))
(unless (\$listp vars) (setq vars `((mlist) ,vars)))
@@ -238,8 +155,8 @@
(unless (symbolp var)
(merror (intl:gettext "rk: variable name expected; found: ~M") var)))
-      (merror (intl:gettext "rk: variable name expected; found: ~M")
+    (merror (intl:gettext "rk: variable name expected; found: ~M")
(setq vars (concatenate 'list '((mlist)) (list (cadr domain)) (cdr vars)))
(setq r (concatenate 'list `(,(car it)) (mapcar #'coerce-float (cdr initial))))
(setq fun (mapcar #'(lambda (x) (coerce-float-fun x vars)) (cdr exprs)))
@@ -263,9 +180,9 @@
(mapcar #'(lambda (x) (* 1/3 x)) k3)
(mapcar #'(lambda (x) (* 1/6 x)) k4)))
(setq r
-        (concatenate 'list
-           `(,(+ (car it) (* m (caddr it))))
-            (map 'list #'+ (cdr r) (mapcar #'(lambda (x) (* (caddr it) x)) u))))
+            (concatenate 'list
+                         `(,(+ (car it) (* m (caddr it))))
+                         (map 'list #'+ (cdr r) (mapcar #'(lambda (x) (* (caddr it) x)) u))))
(push (cons '(mlist) r) traj)))
(when (< (car r) (cadr it))
(let ((s (- (cadr it) (car r))))
@@ -289,5 +206,6 @@
(concatenate 'list