1 ;;Copyright William F. Schelter 1990, All Rights Reserved
8 /* plot of z^
(1/3)...
*/
9 plot3d
(r^
.33*cos
(th/3),[r
,0,1],[th
,0,6*%pi
],['grid
,12,80],['transform_xy
,polar_to_xy
],['plot_format
,geomview
]);
11 /* plot of z^
(1/2)...
*/
12 plot3d
(r^
.5*cos
(th/2),[r
,0,1],[th
,0,6*%pi
],['grid
,12,80],['transform_xy
,polar_to_xy
],['plot_format
,xmaxima
]);
15 plot3d
([cos
(x)*(3+y
*cos
(x/2)),sin
(x)*(3+y
*cos
(x/2)),y
*sin
(x/2)],[x
,-%pi
,%pi
],[y
,-
1,1],['grid
,50,15]);
18 plot3d
([5*cos
(x)*(cos(x/2)*cos
(y)+sin
(x/2)*sin
(2*y
)+3.0) -
10.0,
19 -
5*sin
(x)*(cos(x/2)*cos
(y)+sin
(x/2)*sin
(2*y
)+3.0),
20 5*(-sin(x/2)*cos
(y)+cos
(x/2)*sin
(2*y
))],[x
,-%pi
,%pi
],[y
,-%pi
,%pi
],
23 plot3d
([cos
(y)*(10.0
+6*cos
(x)), sin
(y)*(10.0
+6*cos
(x)),-
6*sin
(x)],
24 [x
,0,2*%pi
],[y
,0,2*%pi
],['grid
,40,40]);
27 (defclass gnuplot-plot
()
28 ((data :initarg
:data
:initform
"")
29 (pipe :initarg
:pipe
:initform nil
)))
31 (defclass xmaxima-plot
()
32 ((data :initarg
:data
:initform
"")
33 (pipe :initarg
:pipe
:initform nil
)))
35 (defclass geomview-plot
()
36 ((data :initarg
:data
:initform
"")
37 (pipe :initarg
:pipe
:initform nil
)))
39 (defgeneric plot-preamble
(plot options
)
40 (:documentation
"Plots the preamble for a plot."))
42 (defgeneric plot2d-command
(plot fun options range
)
43 (:documentation
"Writes the command that creates a plot."))
45 (defgeneric plot3d-command
(plot functions options titles
)
46 (:documentation
"Writes the command that creates a plot."))
48 (defgeneric plot-shipout
(plot options
&optional output-file
)
49 (:documentation
"Sends the plot commands to the graphic program."))
51 (defun ensure-string (x)
54 ((symbolp x
) (print-invert-case (stripdollar x
)))
55 (t (maybe-invert-string-case (string (implode (strgrind x
)))))))
60 ((listp (car l
)) (flatten2 (car l
) (flatten2 (cdr l
) z
)))
61 ((atom (car l
)) (cons (car l
) (flatten2 (cdr l
) z
)))))
67 (if (and ($listp x
) ($listp y
))
68 (cons '(mlist) (loop for w in
(cdr x
) for u in
(cdr y
) collect w collect u
))
69 (merror (intl:gettext
"join: both arguments must be lists."))))
71 (defun coerce-float (x) ($float
(meval* x
)))
73 (defvar *maxima-plotdir
* "")
74 (declare-top (special *maxima-tempdir
* *maxima-prefix
*))
76 ;; *ROT* AND FRIENDS ($ROT, $ROTATE_PTS, $ROTATE_LIST) CAN PROBABLY GO AWAY !!
77 ;; THEY ARE UNDOCUMENTED AND UNUSED !!
78 (defvar *rot
* (make-array 9 :element-type
'flonum
))
81 ;; Global plot options list; this is a property list.. It is not a
82 ;; Maxima variable, to discourage users from changing it directly; it
83 ;; should be changed via set_plot_option
85 (defvar *plot-options
*
87 ,(if (string= *autoconf-windows
* "true")
90 :grid
(30 30) :run_viewer t
:axes t
91 ;; With adaptive plotting, 29 nticks should be enough; adapt_depth
92 ;; controls the number of splittings adaptive-plotting will do.
93 :nticks
29 :adapt_depth
5
94 :color
($blue $red $green $magenta $black $cyan
)
95 :point_type
($bullet $box $triangle $plus $times $asterisk
)
96 :palette
(((mlist) $gradient $green $cyan $blue $violet
)
97 ((mlist) $gradient $magenta $violet $blue $cyan $green $yellow
98 $orange $red $brown $black
))
99 :gnuplot_preamble
"" :gnuplot_term $default
))
101 (defvar $plot_options
103 ((mlist) $plot_format
104 ,(if (string= *autoconf-windows
* "true")
108 ;; $plot_realpart option is false by default but *plot-realpart* is true
109 ;; because coerce-float-fun is used outside of plot package too.
110 (defvar *plot-realpart
* t
)
112 (defun maybe-realpart (x)
115 (if (zerop1 ($imagpart x
))
119 (defvar *missing-data-indicator
* "NaN")
121 (defvar *gnuplot-stream
* nil
)
122 (defvar *gnuplot-command
* "")
124 (defvar $gnuplot_command
"gnuplot")
126 (defun start-gnuplot-process (path)
127 ;; TODO: Forward gnuplot's stderr stream to maxima's stderr output
128 #+clisp
(setq *gnuplot-stream
* (ext:make-pipe-output-stream path
))
129 ;; TODO: Forward gnuplot's stderr stream to maxima's stderr output
130 #+lispworks
(setq *gnuplot-stream
* (system:open-pipe path
))
131 #+cmu
(setq *gnuplot-stream
*
132 (ext:process-input
(ext:run-program path nil
:input
:stream
133 :output
*error-output
* :wait nil
)))
134 #+scl
(setq *gnuplot-stream
*
135 (ext:process-input
(ext:run-program path nil
:input
:stream
136 :output
*error-output
* :wait nil
)))
137 #+sbcl
(setq *gnuplot-stream
*
138 (sb-ext:process-input
(sb-ext:run-program path nil
140 :output
*error-output
* :wait nil
142 #+gcl
(setq *gnuplot-stream
*
143 (open (concatenate 'string
"| " path
) :direction
:output
))
145 (setq *gnuplot-stream
* (ext:run-program path nil
:input
:stream
:output
*error-output
* :error
:output
:wait nil
)))
146 #+ccl
(setf *gnuplot-stream
*
147 (ccl:external-process-input-stream
148 (ccl:run-program path nil
149 :wait nil
:output
*error-output
*
151 #+allegro
(setf *gnuplot-stream
* (excl:run-shell-command
152 path
:input
:stream
:output
*error-output
* :wait nil
))
153 #+abcl
(setq *gnuplot-stream
* (system::process-input
(system::run-program path nil
:wait nil
)))
154 #-
(or clisp cmu sbcl gcl scl lispworks ecl ccl allegro abcl
)
155 (merror (intl:gettext
"plotting: I don't know how to tell this Lisp to run Gnuplot."))
157 (if (null *gnuplot-stream
*)
158 (merror (intl:gettext
"plotting: I tried to execute ~s but *GNUPLOT-STREAM* is still null.~%") path
))
160 ;; set mouse must be the first command send to gnuplot
161 (send-gnuplot-command "set mouse"))
163 (defun check-gnuplot-process ()
164 (if (null *gnuplot-stream
*)
165 (start-gnuplot-process $gnuplot_command
)))
167 (defmfun $gnuplot_close
()
168 (stop-gnuplot-process)
171 (defmfun $gnuplot_start
()
172 (check-gnuplot-process)
175 (defmfun $gnuplot_restart
()
179 (defun stop-gnuplot-process ()
180 (unless (null *gnuplot-stream
*)
182 (close *gnuplot-stream
*)
183 (setq *gnuplot-stream
* nil
))))
185 (defun send-gnuplot-command (command &optional recursive
)
186 (if (null *gnuplot-stream
*)
187 (start-gnuplot-process $gnuplot_command
))
188 (handler-case (unless (null command
)
189 (format *gnuplot-stream
* "~a ~%" command
)
190 (finish-output *gnuplot-stream
*))
192 ;; allow gnuplot to restart if stream-error, or just an error is signaled
193 ;; only try to restart once, to prevent an infinite loop
197 (warn "~a~%Trying new stream.~%" e
)
198 (setq *gnuplot-stream
* nil
)
199 (send-gnuplot-command command t
))))))
201 (defmfun $gnuplot_reset
()
202 (send-gnuplot-command "unset output")
203 (send-gnuplot-command "reset"))
205 (defmfun $gnuplot_replot
(&optional s
)
206 (if (null *gnuplot-stream
*)
207 (merror (intl:gettext
"gnuplot_replot: Gnuplot is not running.")))
209 (send-gnuplot-command "replot"))
211 (send-gnuplot-command s
)
212 (send-gnuplot-command "replot"))
214 (merror (intl:gettext
"gnuplot_replot: argument, if present, must be a string; found: ~M") s
)))
217 ;; allow this to be set in a system init file (sys-init.lsp)
219 (defmfun $get_plot_option
(&optional name n
)
221 ;; Converts the options property list into a Maxima list
222 (do* ((list (copy-tree *plot-options
*) (cddr list
))
223 (key (first list
) (first list
))
224 (value (second list
) (second list
)))
226 (let ((max-key (intern (concatenate 'string
"$" (symbol-name key
)))))
228 (push (cons '(mlist) (cons max-key value
)) options
)
229 (push (list '(mlist) max-key value
) options
))))
230 (setf options
(cons '(mlist) (nreverse options
)))
232 (let ((value (find name
(cdr options
) :key
#'second
)))
238 (defun quote-strings (opt)
241 (format nil
"~s" opt
)
243 (cons (quote-strings (car opt
))
244 (quote-strings (cdr opt
)))))
246 (defun get-plot-option-string (option &optional
(index 1))
247 (let* ((val ($get_plot_option option
2))
248 (val-list (if ($listp val
)
251 (ensure-string (nth (mod (- index
1) (length val-list
)) val-list
))))
253 (defmfun $set_plot_option
(&rest value
)
254 (setq *plot-options
* (plot-options-parser value
*plot-options
*))
257 (defmfun $remove_plot_option
(name)
260 ($adapt_depth
:adapt_depth
) ($axes
:axes
) ($azimuth
:azimuth
)
261 ($box
:box
) ($color
:color
) ($color_bar
:color_bar
)
262 ($color_bar_tics
:color_bar_tics
) ($elevation
:elevation
)
263 ($grid
:grid
) ($grid2d
:grid2d
) ($iterations
:iterations
)
264 ($label
:label
) ($legend
:legend
) ($logx
:logx
) ($logy
:logy
)
265 ($mesh_lines_color
:mesh_lines_color
) ($nticks
:nticks
)
266 ($palette
:palette
) ($plotepsilon
:plotepsilon
)
267 ($plot_format
:plot_format
) ($plot_realpart
:plot_realpart
)
268 ($point_type
:point_type
) ($pdf_file
:pdf_file
)
269 ($png_file
:png_file
) ($ps_file
:ps_file
)
270 ($run_viewer
:run_viewer
) ($same_xy
:samexy
)
271 ($same_xyz
:same_xyz
) ($sample
:sample
) ($style
:style
)
272 ($svg_file
:svg_file
) ($t
:t
) ($title
:title
)
273 ($transform_xy
:transform_xy
) ($x
:x
) ($xbounds
:xbounds
)
274 ($xlabel
:xlabel
) ($xtics
:xtics
) ($xy_scale
:xy_scale
)
275 ($y
:y
) ($ybounds
:ybounds
) ($ylabel
:ylabel
) ($ytics
:ytics
)
276 ($yx_ratio
:yx_ratio
) ($z
:z
) ($zlabel
:zlabel
) ($zmin
:zmin
)
278 ($gnuplot_4_0
:gnuplot_4_0
)
279 ($gnuplot_curve_titles
:gnuplot_curve_titles
)
280 ($gnuplot_curve_styles
:gnuplot_curve_styles
)
281 ($gnuplot_default_term_command
:gnuplot_default_term_command
)
282 ($gnuplot_dumb_term_command
:gnuplot_dumb_term_command
)
283 ($gnuplot_out_file
:gnuplot_out_file
)
284 ($gnuplot_pm3d
:gnuplot_pm3d
)
285 ($gnuplot_strings
:gnuplot_strings
)
286 ($gnuplot_preamble
:gnuplot_preamble
)
287 ($gnuplot_postamble
:gnuplot_postamble
)
288 ($gnuplot_pdf_term_command
:gnuplot_pdf_term_command
)
289 ($gnuplot_png_term_command
:gnuplot_png_term_command
)
290 ($gnuplot_ps_term_command
:gnuplot_ps_term_command
)
291 ($gnuplot_svg_term_command
:gnuplot_svg_term_command
)
292 ($gnuplot_term
:gnuplot_term
))))
294 (defun get-gnuplot-term (term)
295 (let* ((sterm (string-downcase (ensure-string term
)))
296 (pos (search " " sterm
)))
301 (defvar $pstream nil
)
303 (defun print-pt1 (f str
)
306 (format str
"~a " *missing-data-indicator
*)))
308 (defstruct (polygon (:type list
)
309 (:constructor %make-polygon
(pts edges
)))
310 (dummy '($polygon simp
))
315 #-gcl
(:compile-toplevel
:execute
)
317 (defmacro z-pt
(ar i
) `(aref ,ar
(the fixnum
(+ 2 (* ,i
3)))))
318 (defmacro y-pt
(ar i
) `(aref ,ar
(the fixnum
(1+ (* ,i
3)))))
319 (defmacro x-pt
(ar i
) `(aref ,ar
(the fixnum
(* ,i
3))))
320 (defmacro rot
(m i j
) `(aref ,m
(the fixnum
(+ ,i
(the fixnum
(* 3 ,j
))))))
322 (defmacro print-pt
(f)
323 `(print-pt1 ,f $pstream
))
325 (defmacro make-polygon
(a b
)
326 `(list '($polygon
) ,a
,b
)))
328 (defun draw3d (f minx maxx miny maxy nxint nyint
)
329 (let* ((epsx (/ (- maxx minx
) nxint
))
331 (epsy (/ (- maxy miny
) nyint
))
335 (ar (make-array (+ 12 ; 12 for axes
336 (* 3 nx ny
)) :fill-pointer
(* 3 nx ny
)
337 :element-type t
:adjustable t
)))
338 (declare (type flonum x y epsy epsx
)
340 (type (cl:array t
) ar
))
342 initially
(setq y miny
)
348 (setf (z-pt ar l
) (funcall f x y
))
353 (make-polygon ar
(make-grid-vertices nxint nyint
))))
355 ;; The following is 3x2 = 6 rectangles
356 ;; call (make-vertices 3 2)
357 ;; there are 4x3 = 12 points.
358 ;; ordering is x0,y0,z0,x1,y1,z1,....,x11,y11,z11
365 (defun make-grid-vertices (nx ny
)
366 (declare (fixnum nx ny
))
367 (let* ((tem (make-array (+ 15 (* 5 nx ny
)) :fill-pointer
(* 5 nx ny
)
369 :element-type
'(mod #x80000000
)))
374 (declare (fixnum i nxpt m
)
375 (type (cl:array
(mod #x80000000
)) tem
))
376 (loop for k below
(length tem
)
378 (setf (aref tem k
) i
)
379 (setf (aref tem
(incf k
))
381 (setf (aref tem
(incf k
))
383 (setf (aref tem
(incf k
)) i
)
384 (setf (aref tem
(incf k
)) 0) ;place for max
392 (defmfun $rotation1
(phi th
)
393 (let ((sinph (sin phi
))
398 ((mlist simp
) ,(* cosph costh
)
399 ,(* -
1.0 cosph sinth
)
401 ((mlist simp
) ,sinth
,costh
0.0)
402 ((mlist simp
) ,(- (* sinph costh
))
406 ;; pts is a vector of bts [x0,y0,z0,x1,y1,z1,...] and each tuple xi,yi,zi is rotated
407 #-abcl
(defmfun $rotate_pts
(pts rotation-matrix
)
408 (or ($matrixp rotation-matrix
) (merror (intl:gettext
"rotate_pts: second argument must be a matrix.")))
411 (x 0.0) (y 0.0) (z 0.0)
413 (declare (type flonum x y z
))
414 (declare (type (cl:array flonum
) rot
))
415 ($copy_pts rotation-matrix
*rot
* 0)
420 (setq x
(aref pts j
))
421 (setq y
(aref pts
(+ j
1)))
422 (setq z
(aref pts
(+ j
2)))
423 (loop for i below
3 with a of-type flonum
= 0.0
425 (setq a
(* x
(aref rot
(+ (* 3 i
) 0))))
426 (setq a
(+ a
(* y
(aref rot
(+ (* 3 i
) 1)))))
427 (setq a
(+ a
(* z
(aref rot
(+ (* 3 i
) 2)))))
428 (setf (aref pts
(+ j i
)) a
))
431 (defmfun $rotate_list
(x)
432 (cond ((and ($listp x
) (not (mbagp (second x
))))
433 ($list_matrix_entries
(ncmul2 $rot x
)))
434 ((mbagp x
) (cons (car x
) (mapcar '$rotate_list
(cdr x
))))))
436 (defmfun $get_range
(pts k
&aux
(z 0.0) (max most-negative-flonum
) (min most-positive-flonum
))
437 (declare (type flonum z max min
))
438 (declare (type (vector flonum
) pts
))
439 (loop for i from k below
(length pts
) by
3
440 do
(setq z
(aref pts i
))
441 (cond ((< z min
) (setq min z
)))
442 (cond ((> z max
) (setq max z
))))
443 (list min max
(- max min
)))
445 (defmfun $polar_to_xy
(pts &aux
(r 0.0) (th 0.0))
446 (declare (type flonum r th
))
447 (declare (type (cl:array t
) pts
))
448 (assert (typep pts
'(vector t
)))
449 (loop for i below
(length pts
) by
3
450 do
(setq r
(aref pts i
))
451 (setq th
(aref pts
(+ i
1)))
452 (setf (aref pts i
) (* r
(cos th
)))
453 (setf (aref pts
(+ i
1)) (* r
(sin th
)))))
455 ;; Transformation from spherical coordinates to rectangular coordinates,
456 ;; to be used in plot3d. Example of its use:
457 ;; plot3d (expr, [th, 0, %pi], [ph, 0, 2*%pi], [transform_xy, spherical_to_xyz])
458 ;; where expr gives the value of r in terms of the inclination (th)
461 (defmfun $spherical_to_xyz
(pts &aux
(r 0.0) (th 0.0) (ph 0.0))
462 (declare (type flonum r th ph
))
463 (declare (type (cl:array t
) pts
))
464 (assert (typep pts
'(vector t
)))
465 (loop for i below
(length pts
) by
3
466 do
(setq th
(aref pts i
))
467 (setq ph
(aref pts
(+ i
1)))
468 (setq r
(aref pts
(+ i
2)))
469 (setf (aref pts i
) (* r
(sin th
) (cos ph
)))
470 (setf (aref pts
(+ i
1)) (* r
(sin th
) (sin ph
)))
471 (setf (aref pts
(+ i
2)) (* r
(cos th
)))))
474 ;; return a function suitable for the transform function in plot3d.
475 ;; FX, FY, and FZ are functions of three arguments.
476 (defmfun $make_transform
(lvars fx fy fz
)
477 (setq fx
(coerce-float-fun fx lvars
))
478 (setq fy
(coerce-float-fun fy lvars
))
479 (setq fz
(coerce-float-fun fz lvars
))
480 (let ((sym (gensym "transform")))
481 (setf (symbol-function sym
)
482 #'(lambda (pts &aux
(x1 0.0)(x2 0.0)(x3 0.0))
483 (declare (type flonum x1 x2 x3
))
484 (declare (type (cl:array t
) pts
))
485 (loop for i below
(length pts
) by
3
487 (setq x1
(aref pts i
))
488 (setq x2
(aref pts
(+ i
1)))
489 (setq x3
(aref pts
(+ i
2)))
490 (setf (aref pts i
) (funcall fx x1 x2 x3
))
491 (setf (aref pts
(+ 1 i
)) (funcall fy x1 x2 x3
))
492 (setf (aref pts
(+ 2 i
)) (funcall fz x1 x2 x3
)))))))
494 ;; Return value is a Lisp function which evaluates EXPR to a float.
495 ;; COERCE-FLOAT-FUN always returns a function and never returns a symbol,
496 ;; even if EXPR is a symbol.
498 ;; Following cases are recognized:
500 ;; name of a Lisp function
501 ;; name of a Maxima function
502 ;; name of a DEFMSPEC function
503 ;; name of a Maxima macro
504 ;; a string which is the name of a Maxima operator (e.g., "!")
505 ;; name of a simplifying function
506 ;; EXPR is a Maxima lambda expression
507 ;; EXPR is a general Maxima expression
509 ;; %COERCE-FLOAT-FUN is the main internal routine for this.
510 ;; COERCE-FLOAT-FUN is the user interface for creating a function that
511 ;; returns floats. COERCE-BFLOAT-FUN is the same, except bfloats are
513 (defun %coerce-float-fun
(float-fun expr
&optional lvars
)
514 (cond ((and (consp expr
) (functionp expr
))
515 (let ((args (if lvars
(cdr lvars
) (list (gensym)))))
516 (coerce-lisp-function-or-lisp-lambda args expr
:float-fun float-fun
)))
517 ;; expr is a string which names an operator
518 ;; (e.g. "!" "+" or a user-defined operator)
519 ((and (stringp expr
) (getopr0 expr
))
520 (let ((a (if lvars lvars
`((mlist) ,(gensym)))))
521 (%coerce-float-fun float-fun
`(($apply
) ,(getopr0 expr
) ,a
) a
)))
522 ((and (symbolp expr
) (not (member expr lvars
)) (not ($constantp expr
)))
525 (let ((args (if lvars
(cdr lvars
) (list (gensym)))))
526 (coerce-lisp-function-or-lisp-lambda args expr
:float-fun float-fun
)))
528 ;; expr is name of a Maxima function defined by := or
532 ((mexpr (mget expr
'mexpr
))
533 (args (cdr (second mexpr
))))
534 (coerce-maxima-function-or-maxima-lambda args expr
:float-fun float-fun
)))
537 ;; expr is the name of a function defined by defmspec
539 ;; expr is the name of a Maxima macro defined by ::=
541 ;; expr is the name of a simplifying function, and the
542 ;; simplification property is associated with the noun
544 (get ($nounify expr
) 'operators
)
545 ;; expr is the name of a simplifying function, and the
546 ;; simplification property is associated with the verb
548 (get ($verbify expr
) 'operators
))
549 (let ((a (if lvars lvars
`((mlist) ,(gensym)))))
550 (%coerce-float-fun float-fun
`(($apply
) ,expr
,a
) a
)))
552 (merror (intl:gettext
"COERCE-FLOAT-FUN: no such Lisp or Maxima function: ~M") expr
))))
554 ((and (consp expr
) (eq (caar expr
) 'lambda
))
555 (let ((args (cdr (second expr
))))
556 (coerce-maxima-function-or-maxima-lambda args expr
:float-fun float-fun
)))
559 (let* ((vars (or lvars
($sort
($listofvars expr
))))
560 (subscripted-vars ($sublist vars
'((lambda) ((mlist) $x
) ((mnot) (($atom
) $x
)))))
561 gensym-vars save-list-gensym subscripted-vars-save
562 subscripted-vars-mset subscripted-vars-restore
)
564 ;; VARS and SUBSCRIPTED-VARS are Maxima lists. Other lists are
566 (when (cdr subscripted-vars
)
567 (setq gensym-vars
(mapcar #'(lambda (ign) (declare (ignore ign
)) (gensym))
568 (cdr subscripted-vars
)))
569 (mapcar #'(lambda (a b
) (setq vars
(subst b a vars
:test
'equal
)))
570 (cdr subscripted-vars
) gensym-vars
)
572 ;; This stuff about saving and restoring array variables
573 ;; should go into MBINDING, and the lambda expression
574 ;; constructed below should call MBINDING. (At present
575 ;; MBINDING barfs on array variables.)
576 (setq save-list-gensym
(gensym))
577 (setq subscripted-vars-save
578 (mapcar #'(lambda (a) `(push (meval ',a
) ,save-list-gensym
))
579 (cdr subscripted-vars
)))
580 (setq subscripted-vars-mset
581 (mapcar #'(lambda (a b
) `(mset ',a
,b
))
582 (cdr subscripted-vars
) gensym-vars
))
583 (setq subscripted-vars-restore
584 (mapcar #'(lambda (a) `(mset ',a
(pop ,save-list-gensym
)))
585 (reverse (cdr subscripted-vars
)))))
589 (declare (special ,@(cdr vars
) errorsw
))
591 ;; Nothing interpolated here when there are no subscripted
593 ,@(if save-list-gensym
`((declare (special ,save-list-gensym
))))
595 ;; Nothing interpolated here when there are no subscripted
597 ,@(if (cdr subscripted-vars
)
598 `((progn (setq ,save-list-gensym nil
)
599 ,@(append subscripted-vars-save subscripted-vars-mset
))))
601 (let (($ratprint nil
)
602 ;; We don't want to set $numer to T when coercing
603 ;; to a bigfloat. By doing so, things like
604 ;; log(400)^400 get converted to double-floats,
605 ;; which causes a double-float overflow. But the
606 ;; whole point of coercing to bfloat is to use
607 ;; bfloats, not doubles.
609 ;; Perhaps we don't even need to do this for
610 ;; double-floats? It would be nice to remove
611 ;; this. For backward compatibility, we bind
612 ;; numer to T if we're not trying to bfloat.
613 ($numer
,(not (eq float-fun
'$bfloat
)))
617 (declare (special errcatch
))
618 ;; Catch any errors from evaluating the
619 ;; function. We're assuming that if an error
620 ;; is caught, the result is not a number. We
621 ;; also assume that for such errors, it's
622 ;; because the function is not defined there,
623 ;; not because of some other maxima error.
625 ;; GCL 2.6.2 has handler-case but not quite ANSI yet.
630 (,float-fun
(maybe-realpart (meval* ',expr
))))
631 ;; Should we just catch all errors here? It is
632 ;; rather nice to only catch errors we care
633 ;; about and let other errors fall through so
634 ;; that we don't pretend to do something when
635 ;; it is better to let the error through.
636 (arithmetic-error () t
)
637 (maxima-$error
() t
))
641 (,float-fun
(maybe-realpart (meval* ',expr
))))
645 ;; Nothing interpolated here when there are no
646 ;; subscripted variables.
647 ,@(if (cdr subscripted-vars
) `((progn ,@subscripted-vars-restore
)))
652 (defun coerce-float-fun (expr &optional lvars
)
653 (%coerce-float-fun
'$float expr lvars
))
655 (defun coerce-bfloat-fun (expr &optional lvars
)
656 (%coerce-float-fun
'$bfloat expr lvars
))
658 (defun coerce-maxima-function-or-maxima-lambda (args expr
&key
(float-fun '$float
))
659 (let ((gensym-args (loop for x in args collect
(gensym))))
661 `(lambda ,gensym-args
(declare (special ,@gensym-args
))
662 (let* (($ratprint nil
)
667 (declare (special errcatch
))
668 ;; Just always try to convert the result to a float,
669 ;; which handles things like $%pi. See also BUG
670 ;; https://sourceforge.net/p/maxima/bugs/1795/
672 ;; Should we use HANDLER-CASE like we do above in
673 ;; %coerce-float-fun? Seems not necessary for what we want
677 (maybe-realpart (mapply ',expr
(list ,@gensym-args
) t
))))))
680 ;; Same as above, but call APPLY instead of MAPPLY.
682 (defun coerce-lisp-function-or-lisp-lambda (args expr
&key
(float-fun '$float
))
683 (let ((gensym-args (loop for x in args collect
(gensym))))
685 `(lambda ,gensym-args
(declare (special ,@gensym-args
))
686 (let* (($ratprint nil
)
689 (result (maybe-realpart (apply ',expr
(list ,@gensym-args
)))))
690 ;; Always use $float. See comment for
691 ;; coerce-maxima-function-ormaxima-lambda above.
692 (,float-fun result
)))
695 (defmacro zval
(points verts i
) `(aref ,points
(+ 2 (* 3 (aref ,verts
,i
)))))
697 ;;sort the edges array so that drawing the edges will happen from the back towards
698 ;; the front. The if n==4 the edges array coming in looks like
699 ;; v1 v2 v3 v4 0 w1 w2 w3 w4 0 ...
700 ;; where vi,wi are indices pointint into the points array specifying a point
701 ;; in 3 space. After the sorting is done, the 0 is filled in with the vertex
702 ;; which is closer to us (ie highest z component after rotating towards the user)
703 ;; and this is then they are sorted in groups of 5.
704 (defun sort-ngons (points edges n
&aux lis
)
705 (declare (type (cl:array
(flonum)) points
)
706 (type (cl:array
(mod #x80000000
)) edges
)
708 (let ((new (make-array (length edges
) :element-type
(array-element-type edges
)))
714 (leng (length edges
))
716 (declare (type (cl:array
(mod #x80000000
)) new
)
717 (fixnum i leng n1 at
)
719 (declare (type flonum z z1
))
722 (loop for i0 below leng by
(+ n
1)
726 (setq z
(zval points edges i
))
729 do
(if (> (setq z1
(zval points edges i
)) z
)
730 (setq z z1 at
(aref edges i
) ))
733 (setf (aref edges i
) at
)
734 collect
(cons z i0
)))
735 (setq lis
(sort lis
#'alphalessp
:key
#'car
))
739 (loop for j from
(cdr v
)
741 do
(setf (aref new i
) (aref edges j
))
744 (copy-array-portion edges new
0 0 (length edges
))
747 (defun copy-array-portion (ar1 ar2 i1 i2 n1
)
748 (declare (fixnum i1 i2 n1
))
749 (loop while
(>= (setq n1
(- n1
1)) 0)
750 do
(setf (aref ar1 i1
) (aref ar2 i2
))
755 (defmfun $concat_polygons
(pl1 pl2
&aux tem new
)
759 for l
= (+ (length v
) (length w
))
760 do
(setq tem
(make-array l
761 :element-type
(array-element-type v
)
766 (setq new
(make-polygon (first new
) (second new
)) )
768 (copy-array-portion (polygon-pts pl1
) (polygon-pts new
)
769 0 0 (length (polygon-pts pl1
)))
770 (copy-array-portion (polygon-pts pl2
) (polygon-pts new
)
771 (length (polygon-pts pl1
))
772 0 (length (polygon-pts pl2
)))
773 (copy-array-portion (polygon-edges pl1
) (polygon-edges new
)
774 0 0 (length (polygon-edges pl1
)))
775 (loop for i from
(length (polygon-edges pl1
))
776 for j from
0 below
(length (polygon-edges pl2
))
777 with lpts1
= (length (polygon-pts pl1
))
778 with ar2
= (polygon-edges pl2
)
779 with arnew
= (polygon-edges new
)
780 do
(setf (aref arnew i
) (+ lpts1
(aref ar2 j
)))))
782 (defmfun $copy_pts
(lis vec start
)
783 (declare (fixnum start
))
785 (declare (type (cl:array flonum
) tem
))
787 (or (typep lis
'flonum
) (setq lis
(float lis
)))
788 (setf (aref tem start
) lis
)
791 ($copy_pts
(cdr lis
) vec
($copy_pts
(car lis
) vec start
)))
792 ((symbolp lis
) start
)
793 (t (merror (intl:gettext
"copy_pts: unrecognized first argument: ~M") lis
)))))
795 ;; Explicit expressions of two variables, for instance, x and y,
796 ;; where expr is of the form f(x,y) = g(x,y).
797 ;; The result is a series of separated line segments.
798 (defun draw2d-implicit (expr options
)
799 (let ((xmin (first (getf options
:x
)))
800 (ymin (first (getf options
:y
)))
801 (xmax (second (getf options
:x
)))
802 (ymax (second (getf options
:y
)))
803 (gridx (or (first (getf options
:sample
)) 50))
804 (gridy (or (second (getf options
:sample
)) 50))
805 (eps (or (getf options
:plotepsilon
) 0.00000001))
806 vx vy dx dy fun
(result nil
))
807 (setq dx
(/ (- xmax xmin
) gridx
))
808 (setq dy
(/ (- ymax ymin
) gridy
))
809 (setq expr
(m- ($lhs expr
) ($rhs expr
)))
810 (setq vx
(getf options
:xvar
))
811 (setq vy
(getf options
:yvar
))
812 (setq fun
(coerce-float-fun expr
`((mlist) ,vx
,vy
)))
814 ;; Implicit functions algorithm by Jaime Villate. 2021
816 ;; The domain is divided into a grid of rectangles,
817 ;; each one with two triangles, with points labelled as follows:
819 ;; lu ______ ru l=left, r=right, d=down, u=up
821 ;; | / | function fun has the following values at those points:
823 ;; ld |/___| rd fld, frd, flu, fru
827 ((interp (xi yi fi xj yj fj
&aux xp yp fp
)
828 (if (< (* fi fj
) eps
)
830 (setq xp
(/ (- (* fi xj
) (* fj xi
)) (- fi fj
)))
831 (setq yp
(/ (- (* fi yj
) (* fj yi
)) (- fi fj
)))
832 (setq fp
(funcall fun xp yp
))
833 (if (and (< (abs fp
) (abs fi
)) (< (abs fp
) (abs fj
)))
838 (push (first p1
) result
)
839 (push (second p1
) result
)
840 (push (first p2
) result
)
841 (push (second p2
) result
)
842 (push 'moveto result
)
843 (push 'moveto result
)))
846 (xr (+ xmin dx
) (+ xr dx
))
847 (flm (funcall fun xmin ymin
) frm
)
848 (frm (funcall fun
(+ xmin dx
) ymin
) (funcall fun
(+ xr dx
) ymin
)))
852 (yu (+ ymin dy
) (+ yu dy
))
855 (flu (funcall fun xl
(+ ymin dy
)) (funcall fun xl
(+ yu dy
)))
856 (fru (funcall fun xr
(+ ymin dy
)) (funcall fun xr
(+ yu dy
))))
859 (if (not (numberp fld
))
865 (numberp frd
) (numberp flu
) (numberp fru
)
866 (setq p1
(interp xr yd frd xr yu fru
))
867 (setq p2
(interp xl yu flu xr yu fru
)))
868 ;; line between segments rd-ru and lu-ru
870 (when (not (numberp fru
))
875 (numberp frd
) (numberp flu
)
876 (setq p1
(interp xl yd fld xr yd frd
))
877 (setq p2
(interp xl yd fld xl yu flu
)))
878 ;; line between segments ld-rd and ld-lu
880 (when (and next
(< (abs fld
) eps
))
885 (setq p2
(interp xl yu flu xr yu fru
)))
886 ;; line from ld to segment lu-ru
887 (plot-line (list xl yd
) p2
))
891 (setq p2
(interp xr yd frd xr yu fru
)))
892 ;; line from lu to segment rd-ru
893 (plot-line (list xl yd
) p2
))
895 (when (and next
(< (abs fru
) eps
))
900 (setq p2
(interp xl yd fld xr yd frd
)))
901 ;; line from ru to segment ld-rd
902 (plot-line (list xr yu
) p2
))
906 (setq p2
(interp xl yd fld xl yu flu
)))
907 ;; line from ru to segment ld-lu
908 (plot-line (list xr yu
) p2
))
911 (if (setq p1
(interp xl yd fld xr yu fru
))
912 ;; zero in segment ld-ru
915 (when (< (abs flu
) eps
)
916 ;; line from segment ld-ru to lu
917 (plot-line p1
(list xl yu
)))
918 (if (setq p2
(interp xl yd fld xl yu flu
))
919 ;; line between segments ld-ru and ld-lu
923 p2
(interp xl yu flu xr yu fru
))
924 ;; line between segments ld-ru and lu-ru
927 (when (< (abs frd
) eps
)
928 ;; line from segment ld-ru to rd
929 (plot-line p1
(list xr yd
)))
930 (if (setq p2
(interp xl yd fld xr yd frd
))
931 ;; line between segments ld-ru and ld-rd
935 p2
(interp xr yd frd xr yu fru
))
936 ;; line between segments ld-ru and rd-ru
937 (plot-line p1 p2
)))))
942 (setq p1
(interp xl yd fld xl yu flu
))
943 (setq p2
(interp xl yu flu xr yu fru
)))
944 ;; line between segments ld-lu and lu-ru
949 (setq p1
(interp xl yd fld xr yd frd
))
950 (setq p2
(interp xr yd frd xr yu fru
)))
951 ;; line between segments ld-rd and rd-ru
952 (plot-line p1 p2
)))))))))
953 (cons '(mlist) (reverse result
))))
955 ;; parametric ; [parametric,xfun,yfun,[t,tlow,thigh],[nticks ..]]
956 ;; the rest of the parametric list after the list will add to the plot options
958 (defun draw2d-parametric-adaptive (param options
&aux range
)
959 (or (= ($length param
) 4)
960 (merror (intl:gettext
"plot2d: parametric plots must include two expressions and an interval")))
961 (setq range
(nth 4 param
))
962 (or (and ($listp range
) (symbolp (second range
)) (eql ($length range
) 3))
963 (merror (intl:gettext
"plot2d: wrong interval for parametric plot: ~M") range
))
964 (setq range
(check-range range
))
965 (let* ((nticks (getf options
:nticks
))
966 (trange (cddr range
))
967 (tvar (second range
))
968 (xrange (or (getf options
:x
) (getf options
:xbounds
)))
969 (yrange (or (getf options
:y
) (getf options
:ybounds
)))
970 (tmin (coerce-float (first trange
)))
971 (tmax (coerce-float (second trange
)))
972 (xmin (coerce-float (first xrange
)))
973 (xmax (coerce-float (second xrange
)))
974 (ymin (coerce-float (first yrange
)))
975 (ymax (coerce-float (second yrange
)))
977 (declare (type flonum ymin ymax xmin xmax tmin tmax
))
978 (setq f1
(coerce-float-fun (third param
) `((mlist), tvar
)))
979 (setq f2
(coerce-float-fun (fourth param
) `((mlist), tvar
)))
981 (let ((n-clipped 0) (n-non-numeric 0)
982 (t-step (/ (- tmax tmin
) (coerce-float nticks
) 2))
983 t-samples x-samples y-samples result
)
984 ;; Divide the range into 2*NTICKS regions that we then
985 ;; adaptively plot over.
986 (dotimes (k (1+ (* 2 nticks
)))
987 (let ((tpar (+ tmin
(* k t-step
))))
988 (push tpar t-samples
)
989 (push (funcall f1 tpar
) x-samples
)
990 (push (funcall f2 tpar
) y-samples
)))
991 (setf t-samples
(nreverse t-samples
))
992 (setf x-samples
(nreverse x-samples
))
993 (setf y-samples
(nreverse y-samples
))
995 ;; Adaptively plot over each region
996 (do ((t-start t-samples
(cddr t-start
))
997 (t-mid (cdr t-samples
) (cddr t-mid
))
998 (t-end (cddr t-samples
) (cddr t-end
))
999 (x-start x-samples
(cddr x-start
))
1000 (x-mid (cdr x-samples
) (cddr x-mid
))
1001 (x-end (cddr x-samples
) (cddr x-end
))
1002 (y-start y-samples
(cddr y-start
))
1003 (y-mid (cdr y-samples
) (cddr y-mid
))
1004 (y-end (cddr y-samples
) (cddr y-end
)))
1009 (cddr (adaptive-parametric-plot
1011 (car t-start
) (car t-mid
) (car t-end
)
1012 (car x-start
) (car x-mid
) (car x-end
)
1013 (car y-start
) (car y-mid
) (car y-end
)
1014 (getf options
:adapt_depth
)
1016 (adaptive-parametric-plot
1018 (car t-start
) (car t-mid
) (car t-end
)
1019 (car x-start
) (car x-mid
) (car x-end
)
1020 (car y-start
) (car y-mid
) (car y-end
)
1021 (getf options
:adapt_depth
)
1023 ;; Fix up out-of-range values and clobber non-numeric values.
1024 (do ((x result
(cddr x
))
1025 (y (cdr result
) (cddr y
)))
1027 (if (and (numberp (car x
)) (numberp (car y
)))
1028 (unless (and (<= ymin
(car y
) ymax
)
1029 (<= xmin
(car x
) xmax
))
1031 (setf (car x
) 'moveto
)
1032 (setf (car y
) 'moveto
))
1034 (incf n-non-numeric
)
1035 (setf (car x
) 'moveto
)
1036 (setf (car y
) 'moveto
))))
1037 ;; Filter out any MOVETO's which do not precede a number.
1038 ;; Code elsewhere in this file expects MOVETO's to
1039 ;; come in pairs, so leave two MOVETO's before a number.
1040 (let ((n (length result
)))
1045 (eq (nth i result
) 'moveto
)
1046 (eq (nth (1+ i
) result
) 'moveto
)
1049 (eq (nth (+ i
2) result
) 'moveto
)))
1050 (setf (nth i result
) nil
)
1051 (setf (nth (1+ i
) result
) nil
))))
1053 (let ((result-sans-nil (delete nil result
)))
1054 (if (null result-sans-nil
)
1056 ((= n-non-numeric
0)
1057 (mtell (intl:gettext
"plot2d: all values were clipped.~%")))
1059 (mtell (intl:gettext
1060 "plot2d: expression evaluates to non-numeric value everywhere in plotting range.~%")))
1062 (mtell (intl:gettext
1063 "plot2d: all values are non-numeric, or clipped.~%"))))
1065 (if (> n-non-numeric
0)
1066 (mtell (intl:gettext
1067 "plot2d: expression evaluates to non-numeric value somewhere in plotting range.~%")))
1069 (mtell (intl:gettext
"plot2d: some values were clipped.~%")))))
1070 (cons '(mlist) result-sans-nil
)))))
1072 ;; draw2d-discrete. Accepts [discrete,[x1,x2,...],[y1,y2,...]]
1073 ;; or [discrete,[[x1,y1]...] and returns [x1,y1,...] or nil, if
1074 ;; non of the points have real values.
1075 ;; Currently any options given are being ignored, because there
1076 ;; are no options specific to the generation of the points.
1077 (defun draw2d-discrete (f)
1078 (let ((x (third f
)) (y (fourth f
)) data gaps
)
1080 (($listp x
) ; x is a list
1082 (($listp
(cadr x
)) ; x1 is a list
1084 ((= (length (cadr x
)) 3) ; x1 is a 2D point
1085 (setq data
(parse-points-xy x
)))
1086 (t ; x1 is not a 2D point
1087 (merror (intl:gettext
"draw2d-discrete: Expecting a point with 2 coordinates; found ~M~%") (cadr x
)))))
1088 (t ; x1 is not a list
1090 (($listp y
) ; y is a list
1092 ((symbolp (coerce-float (cadr y
))); y is an option
1093 (setq data
(parse-points-y x
)))
1094 (t ; y is not an option
1096 (($listp
(cadr y
)) ; y1 is a list
1097 (merror (intl:gettext
"draw2d-discrete: Expecting a y coordinate; found ~M~%") (cadr y
)))
1100 ((= (length x
) (length y
)) ; case [x][y]
1101 (setq data
(parse-points-x-y x y
)))
1103 (merror (intl:gettext
"draw2d-discrete: The number of x and y coordinates do not match.~%")))))))))
1104 (t ; y is not a list
1105 (setq data
(parse-points-y x
)))))))
1106 (t ; x is not a list
1107 (merror (intl:gettext
"draw2d-discrete: Expecting a list of x coordinates or points; found ~M~%") x
)))
1109 ;; checks for non-real values
1111 ((some #'realp data
)
1112 (setq gaps
(count-if #'(lambda (x) (eq x
'moveto
)) data
))
1114 ;; some points have non-real values
1115 (mtell (intl:gettext
"Warning: excluding ~M points with non-numerical values.~%") (/ gaps
2))))
1117 ;; none of the points have real values
1118 (mtell (intl:gettext
"Warning: none of the points have numerical values.~%"))
1122 ;; Two lists [x1...xn] and [y1...yn] are joined as
1123 ;; [x1 y1...xn yn], converting all expressions to real numbers.
1124 ;; If either xi or yi are not real, both are replaced by 'moveto
1125 (defun parse-points-x-y (x y
)
1126 (do ((a (rest x
) (cdr a
))
1127 (b (rest y
) (cdr b
))
1129 ((null b
) (cons '(mlist) (reverse c
)))
1130 (setq af
(coerce-float (car a
)))
1131 (setq bf
(coerce-float (car b
)))
1133 ((or (not (realp af
)) (not (realp bf
)))
1134 (setq c
(cons 'moveto
(cons 'moveto c
))))
1136 (setq c
(cons bf
(cons af c
)))))))
1138 ;; One list [y1...yn] becomes the list [1 y1...n yn],
1139 ;; converting all expressions to real numbers.
1140 ;; If yi is not real, both i and yi are replaced by 'moveto
1141 (defun parse-points-y (y)
1143 (b (rest y
) (cdr b
))
1145 ((null b
) (cons '(mlist) (reverse c
)))
1146 (setq bf
(coerce-float (car b
)))
1149 (setq c
(cons 'moveto
(cons 'moveto c
))))
1151 (setq c
(cons bf
(cons a c
)))))))
1153 ;; List [[x1,y1]...[xn,yn]] is transformed into
1154 ;; [x1 y1...xn yn], converting all expressions to real numbers.
1155 ;; If either xi or yi are not real, both are replaced by 'moveto
1156 (defun parse-points-xy (xy)
1157 (do ((ab (rest xy
) (cdr ab
))
1159 ((null ab
) (cons '(mlist) (reverse c
)))
1160 (setq af
(coerce-float (cadar ab
)))
1161 (setq bf
(coerce-float (caddar ab
)))
1163 ((or (not (realp af
)) (not (realp bf
)))
1164 (setq c
(cons 'moveto
(cons 'moveto c
))))
1166 (setq c
(cons bf
(cons af c
)))))))
1168 ;;; Adaptive plotting, based on the adaptive plotting code from
1169 ;;; YACAS. See http://yacas.sourceforge.net/Algo.html#c3s1 for a
1170 ;;; description of the algorithm. More precise details can be found
1171 ;;; in the file yacas/scripts/plots.rep/plot2d.ys.
1174 ;; Determine if we have a slow oscillation of the function.
1175 ;; Basically, for each 3 consecutive function values, we check to see
1176 ;; if the function is monotonic or not. There are 3 such sets, and
1177 ;; the function is considered slowly oscillating if at most 2 of them
1178 ;; are not monotonic.
1179 (defun slow-oscillation-p (f0 f1 f2 f3 f4
)
1180 (flet ((sign-change (x y z
)
1181 (cond ((not (and (numberp x
) (numberp y
) (numberp z
)))
1182 ;; Something is not a number. Assume the
1183 ;; oscillation is not slow.
1185 ((or (and (> y x
) (> y z
))
1186 (and (< y x
) (< y z
)))
1190 (<= (+ (sign-change f0 f1 f2
)
1191 (sign-change f1 f2 f3
)
1192 (sign-change f2 f3 f4
))
1195 ;; Determine if the function values are smooth enough. This means
1196 ;; that integrals of the functions on the left part and the right part
1197 ;; of the range are approximately the same.
1200 (defun smooth-enough-p (f-a f-a1 f-b f-b1 f-c eps
)
1201 (cond ((every #'numberp
(list f-a f-a1 f-b f-b1 f-c
))
1202 (let ((quad (/ (+ f-a
1208 (quad-b (/ (+ (* 5 f-b
)
1212 ;; According to the Yacas source code, quad is the Simpson
1213 ;; quadrature for the (fb,fb1) subinterval (using points b,b1,c),
1214 ;; subtracted from the 4-point Newton-Cotes quadrature for the
1215 ;; (fb,fb1) subinterval (using points a, a1, b, b1.
1217 ;; quad-b is the Simpson quadrature for the (fb,f1) subinterval.
1219 ;; This used to test for diff <= 0. But in some
1220 ;; situations, like plot2d(0.99,[x,0,5]), roundoff prevents
1221 ;; this from happening. So we do diff < delta instead, for
1222 ;; some value of delta.
1224 ;; XXX: What is the right value for delta? Does this break
1225 ;; other things? Simple tests thus far show that
1226 ;; 100*flonum-epsilon is ok.
1227 (let ((diff (- (abs quad
)
1228 (* eps
(- quad-b
(min f-a f-a1 f-b f-b1 f-c
)))))
1229 (delta (* 150 flonum-epsilon
)))
1232 ;; Something is not a number, so assume it's not smooth enough.
1236 (defun adaptive-plot (fcn a b c f-a f-b f-c depth eps
)
1237 ;; Step 1: Split the interval [a, c] into 5 points
1238 (let* ((a1 (/ (+ a b
) 2))
1240 (f-a1 (funcall fcn a1
))
1241 (f-b1 (funcall fcn b1
))
1243 (cond ((or (not (plusp depth
))
1244 (and (slow-oscillation-p f-a f-a1 f-b f-b1 f-c
)
1245 (smooth-enough-p f-a f-a1 f-b f-b1 f-c eps
)))
1246 ;; Everything is nice and smooth so we're done. Don't
1253 ;; We are not plotting the real part of the function and the
1254 ;; function is undefined at all points - assume it has complex value
1255 ;; on [a,b]. Maybe we should refine it a couple of times just to make sure?
1256 ((and (null *plot-realpart
*)
1257 (null f-a
) (null f-a1
) (null f-b
) (null f-b1
) (null f-c
))
1264 ;; Need to refine. Split the interval in half, and try to plot each half.
1265 (let ((left (adaptive-plot fcn a a1 b f-a f-a1 f-b
(1- depth
) (* 2 eps
)))
1266 (right (adaptive-plot fcn b b1 c f-b f-b1 f-c
(1- depth
) (* 2 eps
))))
1267 (append left
(cddr right
)))))))
1269 (defun adaptive-parametric-plot (x-fcn y-fcn a b c x-a x-b x-c y-a y-b y-c depth eps
)
1270 ;; Step 1: Split the interval [a, c] into 5 points
1271 (let* ((a1 (/ (+ a b
) 2))
1273 (x-a1 (funcall x-fcn a1
))
1274 (x-b1 (funcall x-fcn b1
))
1275 (y-a1 (funcall y-fcn a1
))
1276 (y-b1 (funcall y-fcn b1
)))
1277 (cond ((or (not (plusp depth
))
1278 ;; Should we have a different algorithm to determine
1279 ;; slow oscillation and smooth-enough for parametric
1281 (and (slow-oscillation-p y-a y-a1 y-b y-b1 y-c
)
1282 (slow-oscillation-p x-a x-a1 x-b x-b1 x-c
)
1283 (smooth-enough-p y-a y-a1 y-b y-b1 y-c eps
)
1284 (smooth-enough-p x-a x-a1 x-b x-b1 x-c eps
)))
1285 ;; Everything is nice and smooth so we're done. Don't
1292 ;; We are not plotting the real part of the function and the
1293 ;; function is undefined at all points - assume it has complex value
1294 ;; on [a,b]. Maybe we should refine it a couple of times just to make sure?
1295 ((and (null *plot-realpart
*)
1296 (null y-a
) (null y-a1
) (null y-b
) (null y-b1
) (null y-c
)
1297 (null x-a
) (null x-a1
) (null x-b
) (null x-b1
) (null x-c
))
1304 ;; Need to refine. Split the interval in half, and try to plot each half.
1305 (let ((left (adaptive-parametric-plot x-fcn y-fcn
1309 (1- depth
) (* 2 eps
)))
1310 (right (adaptive-parametric-plot x-fcn y-fcn
1314 (1- depth
) (* 2 eps
))))
1315 ;; (cddr right) to skip over the point that is duplicated
1316 ;; between the right end-point of the left region and the
1317 ;; left end-point of the right
1318 (append left
(cddr right
)))))))
1320 (defun draw2d (fcn range plot-options
)
1321 (if (and ($listp fcn
) (equal '$parametric
(cadr fcn
)))
1323 (draw2d-parametric-adaptive fcn plot-options
)))
1324 (if (and ($listp fcn
) (equal '$discrete
(cadr fcn
)))
1325 (return-from draw2d
(draw2d-discrete fcn
)))
1326 (when (and (listp fcn
) (member 'mequal
(flatten fcn
)))
1327 (return-from draw2d
(draw2d-implicit fcn plot-options
)))
1328 (let* ((nticks (getf plot-options
:nticks
))
1329 (yrange (getf plot-options
:ybounds
))
1330 (depth (getf plot-options
:adapt_depth
)))
1332 (setq fcn
(coerce-float-fun fcn
`((mlist), (second range
))))
1334 (let* ((x-start (coerce-float (third range
)))
1335 (xend (coerce-float (fourth range
)))
1336 (x-step (/ (- xend x-start
) (coerce-float nticks
) 2))
1337 (ymin (coerce-float (first yrange
)))
1338 (ymax (coerce-float (second yrange
)))
1339 (n-clipped 0) (n-non-numeric 0)
1340 ;; What is a good EPS value for adaptive plotting?
1342 x-samples y-samples result
1344 (declare (type flonum ymin ymax
))
1345 ;; Divide the region into NTICKS regions. Each region has a
1346 ;; start, mid and endpoint. Then adaptively plot over each of
1347 ;; these regions. So it's probably a good idea not to make
1348 ;; NTICKS too big. Since adaptive plotting splits the sections
1349 ;; in half, it's also probably not a good idea to have NTICKS be
1351 (when (getf plot-options
:logx
)
1352 (setf x-start
(log x-start
))
1353 (setf xend
(log xend
))
1354 (setf x-step
(/ (- xend x-start
) (coerce-float nticks
) 2)))
1357 (let ((y (if (getf plot-options
:logx
)
1358 (funcall fcn
(exp x
))
1360 (if (and (getf plot-options
:logy
)
1362 (if (> y
0) (log y
) 'und
)
1365 (dotimes (k (1+ (* 2 nticks
)))
1366 (let ((x (+ x-start
(* k x-step
))))
1368 (push (fun x
) y-samples
)))
1369 (setf x-samples
(nreverse x-samples
))
1370 (setf y-samples
(nreverse y-samples
))
1372 ;; For each region, adaptively plot it.
1373 (do ((x-start x-samples
(cddr x-start
))
1374 (x-mid (cdr x-samples
) (cddr x-mid
))
1375 (x-end (cddr x-samples
) (cddr x-end
))
1376 (y-start y-samples
(cddr y-start
))
1377 (y-mid (cdr y-samples
) (cddr y-mid
))
1378 (y-end (cddr y-samples
) (cddr y-end
)))
1380 ;; The region is x-start to x-end, with mid-point x-mid.
1382 ;; The cddr is to remove the one extra sample (x and y value)
1383 ;; that adaptive plot returns. But on the first iteration,
1384 ;; result is empty, so we don't want the cddr because we want
1385 ;; all the samples returned from adaptive-plot. On subsequent
1386 ;; iterations, it's a duplicate of the last point of the
1387 ;; previous interval.
1392 (adaptive-plot #'fun
(car x-start
) (car x-mid
) (car x-end
)
1393 (car y-start
) (car y-mid
) (car y-end
)
1395 (adaptive-plot #'fun
(car x-start
) (car x-mid
) (car x-end
)
1396 (car y-start
) (car y-mid
) (car y-end
)
1399 ;; Fix up out-of-range values
1400 ;; and clobber non-numeric values.
1402 (do ((x result
(cddr x
))
1403 (y (cdr result
) (cddr y
)))
1405 (if (numberp (car y
))
1406 (unless (<= ymin
(car y
) ymax
)
1408 (setf (car x
) 'moveto
)
1409 (setf (car y
) 'moveto
))
1411 (incf n-non-numeric
)
1412 (setf (car x
) 'moveto
)
1413 (setf (car y
) 'moveto
)))
1414 (when (and (getf plot-options
:logx
)
1416 (setf (car x
) (exp (car x
))))
1418 (when (and (getf plot-options
:logy
)
1420 (setf (car y
) (exp (car y
)))))
1422 ;; Filter out any MOVETO's which do not precede a number.
1423 ;; Code elsewhere in this file expects MOVETO's to
1424 ;; come in pairs, so leave two MOVETO's before a number.
1425 (let ((n (length result
)))
1430 (eq (nth i result
) 'moveto
)
1431 (eq (nth (1+ i
) result
) 'moveto
)
1434 (eq (nth (+ i
2) result
) 'moveto
)))
1435 (setf (nth i result
) nil
)
1436 (setf (nth (1+ i
) result
) nil
))))
1438 (let ((result-sans-nil (delete nil result
)))
1439 (if (null result-sans-nil
)
1441 ((= n-non-numeric
0)
1442 (mtell (intl:gettext
"plot2d: all values were clipped.~%")))
1444 (mtell (intl:gettext
"plot2d: expression evaluates to non-numeric value everywhere in plotting range.~%")))
1446 (mtell (intl:gettext
"plot2d: all values are non-numeric, or clipped.~%"))))
1448 (if (> n-non-numeric
0)
1449 (mtell (intl:gettext
"plot2d: expression evaluates to non-numeric value somewhere in plotting range.~%")))
1451 (mtell (intl:gettext
"plot2d: some values were clipped.~%")))))
1452 (cons '(mlist) result-sans-nil
))))))
1454 (defun get-range (lis)
1455 (let ((ymin most-positive-flonum
)
1456 (ymax most-negative-flonum
))
1457 (declare (type flonum ymin ymax
))
1458 (do ((l lis
(cddr l
)))
1460 (or (floatp (car l
)) (setf (car l
) (float (car l
))))
1461 (cond ((< (car l
) ymin
)
1462 (setq ymin
(car l
))))
1463 (cond ((< ymax
(car l
))
1464 (setq ymax
(car l
)))))
1465 (list '(mlist) ymin ymax
)))
1467 #+sbcl
(defvar $gnuplot_view_args
"-persist ~a")
1468 #-sbcl
(defvar $gnuplot_view_args
"-persist ~s")
1470 #+(or sbcl openmcl
) (defvar $gnuplot_file_args
"~a")
1471 #-
(or sbcl openmcl
) (defvar $gnuplot_file_args
"~s")
1473 (defvar $mgnuplot_command
"mgnuplot")
1474 (defvar $geomview_command
"geomview")
1476 (defvar $xmaxima_plot_command
"xmaxima")
1478 (defun plot-temp-file (file &optional
(preserve-file nil
))
1480 (if *maxima-tempdir
*
1481 (format nil
"~a/~a" *maxima-tempdir
* file
)
1483 (declare (special *temp-files-list
*))
1484 (unless preserve-file
1485 (setf (gethash filename
*temp-files-list
*) t
))
1486 (format nil
"~a" filename
)
1489 ;; If no file path is given, uses temporary directory path
1490 (defun plot-file-path (file &optional
(preserve-file nil
))
1491 (if (pathname-directory file
)
1493 (plot-temp-file file preserve-file
)))
1495 (defun gnuplot-process (plot-options &optional file out-file
)
1496 (let ((gnuplot-term (getf plot-options
:gnuplot_term
))
1497 (run-viewer (getf plot-options
:run_viewer
))
1498 #-
(or (and sbcl win32
) (and sbcl win64
) (and ccl windows
))
1500 (string-downcase (getf plot-options
:gnuplot_preamble
))))
1502 ;; creates the output file, when there is one to be created
1503 (when (and out-file
(not (eq gnuplot-term
'$default
)))
1504 #+(or (and sbcl win32
) (and sbcl win64
) (and ccl windows
))
1505 ($system $gnuplot_command
(format nil $gnuplot_file_args file
))
1506 #-
(or (and sbcl win32
) (and sbcl win64
) (and ccl windows
))
1507 ($system
(format nil
"~a ~a" $gnuplot_command
1508 (format nil $gnuplot_file_args file
))))
1510 ;; displays contents of the output file, when gnuplot-term is dumb,
1511 ;; or runs gnuplot when gnuplot-term is default
1515 ;; the options given to gnuplot will be different when the user
1516 ;; redirects the output by using "set output" in the preamble
1517 #+(or (and sbcl win32
) (and sbcl win64
) (and ccl windows
))
1518 ($system $gnuplot_command
"-persist" (format nil $gnuplot_file_args file
))
1519 #-
(or (and sbcl win32
) (and sbcl win64
) (and ccl windows
))
1521 (format nil
"~a ~a" $gnuplot_command
1522 (format nil
(if (search "set out" gnuplot-preamble
)
1523 $gnuplot_file_args $gnuplot_view_args
)
1527 ($printfile
(car out-file
))
1528 (merror (intl:gettext
"plotting: option 'gnuplot_out_file' not defined."))))))))
1530 ;; plot-options-parser puts the plot options given into a property list.
1531 ;; maxopts: a list (not a Maxima list!) with plot options.
1532 ;; options: a property list, or an empty list.
1534 ;; (plot-options-parser (list #$[x,-2,2]$ #$[nticks,30]$) '(:nticks 4))
1536 ;; (:XLABEL "x" :XMAX 2.0 :XMIN -2.0 :NTICKS 30)
1538 (defun plot-options-parser (maxopts options
&aux name
)
1539 (dolist (opt maxopts
)
1540 (unless (or ($listp opt
) (symbolp opt
))
1543 "plot-options-parser: option \"~M\" should be a list or a symbol")
1547 (unless ($symbolp
(setq name
(second opt
)))
1550 "plot-options-parser: Expecting option name as a symbol, found: \"~M\"")
1554 (setf (getf options
:adapt_depth
)
1555 (check-option (cdr opt
) #'naturalp
"a natural number" 1)))
1556 ($axes
(setf (getf options
:axes
)
1557 (check-option-b (cdr opt
) #'axesoptionp
"x, y, solid" 1)))
1558 ($azimuth
(if (caddr opt
)
1559 (setf (caddr opt
) (parse-azimuth (caddr opt
))))
1560 (setf (getf options
:azimuth
)
1561 (check-option (cdr opt
) #'realp
"a real number" 1)))
1562 ($box
(setf (getf options
:box
)
1563 (check-option-boole (cdr opt
))))
1566 (setf (cddr opt
) (mapcar #'coerce-float
(cddr opt
))))
1567 (setf (getf options
:color_bar_tics
)
1568 (check-option-b (cdr opt
) #'realp
"a real number" 3)))
1569 ($color
(setf (getf options
:color
)
1570 (check-option (cdr opt
) #'plotcolorp
"a color")))
1571 ($color_bar
(setf (getf options
:color_bar
)
1572 (check-option-boole (cdr opt
))))
1573 ($elevation
(if (caddr opt
)
1574 (setf (caddr opt
) (parse-elevation (caddr opt
))))
1575 (setf (getf options
:elevation
)
1576 (check-option (cdr opt
) #'realp
"a real number" 1)))
1577 ($grid
(setf (getf options
:grid
)
1578 (check-option (cdr opt
) #'naturalp
"a natural number" 2)))
1579 ($grid2d
(setf (getf options
:grid2d
)
1580 (check-option-boole (cdr opt
))))
1582 (setf (getf options
:iterations
)
1583 (check-option (cdr opt
) #'naturalp
"a natural number" 1)))
1584 ($label
(setf (getf options
:label
)
1585 (check-option-label (cdr opt
))))
1586 ($legend
(setf (getf options
:legend
)
1587 (check-option-b (cdr opt
) #'stringp
"a string")))
1588 ($logx
(setf (getf options
:logx
)
1589 (check-option-boole (cdr opt
))))
1590 ($logy
(setf (getf options
:logy
)
1591 (check-option-boole (cdr opt
))))
1593 (setf (getf options
:mesh_lines_color
)
1594 (check-option-b (cdr opt
) #'plotcolorp
"a color" 1)))
1595 ($nticks
(setf (getf options
:nticks
)
1596 (check-option (cdr opt
) #'naturalp
"a natural number" 1)))
1597 ($palette
(setf (getf options
:palette
)
1598 (check-option-palette (cdr opt
))))
1599 ($plotepsilon
(setf (getf options
:plotepsilon
)
1600 (check-option (cdr opt
) #'realp
"a real number" 1)))
1601 ($plot_format
(setf (getf options
:plot_format
)
1602 (check-option-format (cdr opt
))))
1603 ($plot_realpart
(setf (getf options
:plot_realpart
)
1604 (check-option-boole (cdr opt
))))
1605 ($point_type
(setf (getf options
:point_type
)
1606 (check-option (cdr opt
) #'pointtypep
"a point type")))
1607 ($pdf_file
(setf (getf options
:pdf_file
)
1608 (check-option (cdr opt
) #'stringp
"a string" 1)))
1609 ($png_file
(setf (getf options
:png_file
)
1610 (check-option (cdr opt
) #'stringp
"a string" 1)))
1611 ($ps_file
(setf (getf options
:ps_file
)
1612 (check-option (cdr opt
) #'stringp
"a string" 1)))
1613 ($run_viewer
(setf (getf options
:run_viewer
)
1614 (check-option-boole (cdr opt
))))
1615 ($sample
(setf (getf options
:sample
)
1616 (check-option (cdr opt
) #'naturalp
"a natural number" 2)))
1617 ($same_xy
(setf (getf options
:same_xy
)
1618 (check-option-boole (cdr opt
))))
1619 ($same_xyz
(setf (getf options
:same_xyz
)
1620 (check-option-boole (cdr opt
))))
1621 ($style
(setf (getf options
:style
)
1622 (check-option-style (cdr opt
))))
1623 ($svg_file
(setf (getf options
:svg_file
)
1624 (check-option (cdr opt
) #'stringp
"a string" 1)))
1625 ($t
(setf (getf options
:t
) (cddr (check-range opt
))))
1626 ($title
(setf (getf options
:title
)
1627 (check-option (cdr opt
) #'stringp
"a string" 1)))
1628 ($transform_xy
(setf (getf options
:transform_xy
)
1629 (check-option-b (cdr opt
) #'functionp
"a function make_transform" 1)))
1630 ($x
(setf (getf options
:x
) (cddr (check-range opt
))))
1631 ($xbounds
(setf (getf options
:xbounds
) (cddr (check-range opt
))))
1632 ($xlabel
(setf (getf options
:xlabel
)
1633 (check-option (cdr opt
) #'string
"a string" 1)))
1636 (setf (cddr opt
) (mapcar #'coerce-float
(cddr opt
))))
1637 (setf (getf options
:xtics
)
1638 (check-option-b (cdr opt
) #'realp
"a real number" 3)))
1641 (setf (cddr opt
) (mapcar #'coerce-float
(cddr opt
))))
1642 (setf (getf options
:xy_scale
)
1643 (check-option (cdr opt
) #'realpositivep
1644 "a positive real number" 2)))
1645 ($y
(setf (getf options
:y
) (cddr (check-range opt
))))
1646 ($ybounds
(setf (getf options
:ybounds
) (cddr (check-range opt
))))
1647 ($ylabel
(setf (getf options
:ylabel
)
1648 (check-option (cdr opt
) #'string
"a string" 1)))
1651 (setf (cddr opt
) (mapcar #'coerce-float
(cddr opt
))))
1652 (setf (getf options
:ytics
)
1653 (check-option-b (cdr opt
) #'realp
"a real number" 3)))
1656 (setf (caddr opt
) (coerce-float (caddr opt
))))
1657 (setf (getf options
:yx_ratio
)
1658 (check-option (cdr opt
) #'realp
"a real number" 1)))
1659 ($z
(setf (getf options
:z
) (cddr (check-range opt
))))
1660 ($zlabel
(setf (getf options
:zlabel
)
1661 (check-option (cdr opt
) #'string
"a string" 1)))
1664 (setf (caddr opt
) (coerce-float (caddr opt
))))
1665 (setf (getf options
:zmin
)
1666 (check-option-b (cdr opt
) #'realp
"a real number" 1)))
1669 (setf (cddr opt
) (mapcar #'coerce-float
(cddr opt
))))
1670 (setf (getf options
:ztics
)
1671 (check-option-b (cdr opt
) #'realp
"a real number" 3)))
1672 ($gnuplot_4_0
(setf (getf options
:gnuplot_4_0
)
1673 (check-option-boole (cdr opt
))))
1674 ($gnuplot_curve_titles
1675 (setf (getf options
:gnuplot_curve_titles
)
1676 (check-option (cdr opt
) #'stringp
"a string")))
1677 ($gnuplot_curve_styles
1678 (setf (getf options
:gnuplot_curve_styles
)
1679 (check-option (cdr opt
) #'stringp
"a string")))
1680 ($gnuplot_default_term_command
1681 (setf (getf options
:gnuplot_default_term_command
)
1682 (check-option (cdr opt
) #'stringp
"a string" 1)))
1683 ($gnuplot_dumb_term_command
1684 (setf (getf options
:gnuplot_dumb_term_command
)
1685 (check-option (cdr opt
) #'stringp
"a string" 1)))
1687 (setf (getf options
:gnuplot_out_file
)
1688 (check-option (cdr opt
) #'stringp
"a string" 1)))
1690 (setf (getf options
:gnuplot_pm3d
)
1691 (check-option-boole (cdr opt
))))
1693 (setf (getf options
:gnuplot_strings
)
1694 (check-option-boole (cdr opt
))))
1696 (setf (getf options
:gnuplot_preamble
)
1697 (check-option (cdr opt
) #'stringp
"a string" 1)))
1699 (setf (getf options
:gnuplot_postamble
)
1700 (check-option (cdr opt
) #'stringp
"a string" 1)))
1701 ($gnuplot_pdf_term_command
1702 (setf (getf options
:gnuplot_pdf_term_command
)
1703 (check-option (cdr opt
) #'stringp
"a string" 1)))
1704 ($gnuplot_png_term_command
1705 (setf (getf options
:gnuplot_png_term_command
)
1706 (check-option (cdr opt
) #'stringp
"a string" 1)))
1707 ($gnuplot_ps_term_command
1708 (setf (getf options
:gnuplot_ps_term_command
)
1709 (check-option (cdr opt
) #'stringp
"a string" 1)))
1710 ($gnuplot_svg_term_command
1711 (setf (getf options
:gnuplot_svg_term_command
)
1712 (check-option (cdr opt
) #'stringp
"a string" 1)))
1713 ;; gnuplot_term is a tricky one: when it is just default, dumb or
1714 ;; ps, we want it to be a symbol, but when it is more complicated,
1715 ;; i.e. "ps; size 16cm, 12cm", it must be a string and not a symbol
1717 (let ((s (caddr opt
)))
1719 (cond ((string= s
"default") (setq s
'$default
))
1720 ((string= s
"dumb") (setq s
'$dumb
))
1721 ((string= s
"ps") (setq s
'$ps
))))
1723 (setf (getf options
:gnuplot_term
) s
)
1725 (intl:gettext
"Wrong argument for plot option \"gnuplot_term\". Expecting a string or a symbol but found \"~M\".") s
))))
1728 (intl:gettext
"plot-options-parser: unknown plot option: ~M") opt
))))
1731 ($axes
(setf (getf options
:axes
) t
))
1732 ($box
(setf (getf options
:box
) t
))
1733 ($color_bar
(setf (getf options
:color_bar
) t
))
1734 ($color_bar_tics
(remf options
:color_bar_tics
))
1735 ($grid2d
(setf (getf options
:grid2d
) t
))
1736 ($legend
(remf options
:legend
))
1737 ($logx
(setf (getf options
:logx
) t
))
1738 ($logy
(setf (getf options
:logy
) t
))
1739 ($mesh_lines_color
(remf options
:mesh_lines_color
))
1740 ($plot_realpart
(setf (getf options
:plot_realpart
) t
))
1741 ($run_viewer
(setf (getf options
:run_viewer
) t
))
1742 ($same_xy
(setf (getf options
:same_xy
) t
))
1743 ($same_xyz
(setf (getf options
:same_xyz
) t
))
1744 ($transform_xy
(remf options
:transform_xy
))
1745 ($xtics
(remf options
:xtics
))
1746 ($ytics
(remf options
:ytics
))
1747 ($zmin
(remf options
:zmin
))
1748 ($gnuplot_4_0
(setf (getf options
:gnuplot_4_0
) t
))
1749 ($gnuplot_pm3d
(setf (getf options
:gnuplot_pm3d
) t
))
1750 ($gnuplot_strings
(setf (getf options
:gnuplot_strings
) t
))
1751 ($noaxes
(setf (getf options
:axes
) nil
))
1752 ($nobox
(setf (getf options
:box
) nil
))
1753 ($nocolor_bar
(setf (getf options
:color_bar
) nil
))
1754 ($nocolor_bat_tics
(setf (getf options
:color_bat_tics
) nil
))
1755 ($nogrid2d
(setf (getf options
:grid2d
) nil
))
1756 ($nolegend
(setf (getf options
:legend
) nil
))
1757 ($nologx
(setf (getf options
:logx
) nil
))
1758 ($nology
(setf (getf options
:logy
) nil
))
1759 ($nomesh_lines_color
(setf (getf options
:mesh_lines_color
) nil
))
1760 ($noplot_realpart
(setf (getf options
:plot_realpart
) nil
))
1761 ($norun_viewer
(setf (getf options
:run_viewer
) nil
))
1762 ($nosame_xy
(setf (getf options
:same_xy
) nil
))
1763 ($nosame_xyz
(setf (getf options
:same_xyz
) nil
))
1764 ($notransform_xy
(setf (getf options
:transform_xy
) nil
))
1765 ($noxtics
(setf (getf options
:xtics
) nil
))
1766 ($noytics
(setf (getf options
:ytics
) nil
))
1767 ($nozmin
(setf (getf options
:zmin
) nil
))
1768 ($nognuplot_4_0
(setf (getf options
:gnuplot_4_0
) nil
))
1769 ($nognuplot_pm3d
(setf (getf options
:gnuplot_pm3d
) nil
))
1770 ($nognuplot_strings
(setf (getf options
:gnuplot_strings
) nil
))
1772 (merror (intl:gettext
"Unknown plot option \"~M\".") opt
))))))
1774 ;; plots that create a file work better in gnuplot than gnuplot_pipes
1775 (when (and (eq (getf options
:plot_format
) '$gnuplot_pipes
)
1776 (or (eq (getf options
:gnuplot_term
) '$dumb
)
1777 (getf options
:pdf_file
) (getf options
:png_file
)
1778 (getf options
:ps_file
) (getf options
:svg_file
)))
1779 (setf (getf options
:plot_format
) '$gnuplot
))
1782 ;; natural numbers predicate
1783 (defun naturalp (n) (or (and (integerp n
) (> n
0)) nil
))
1785 ;; positive real numbers predicate
1786 (defun realpositivep (x) (or (and (realp x
) (> x
0)) nil
))
1788 ;; posible values for the axes option
1789 (defun axesoptionp (o) (if (member o
'($x $y $solid
)) t nil
))
1791 ;; the 13 possibilities for the point types
1792 (defun pointtypep (p)
1793 (if (member p
'($bullet $circle $plus $times $asterisk $box $square
1794 $triangle $delta $wedge $nabla $diamond $lozenge
)) t nil
))
1796 ;; Colors can only one of the named colors or a six-digit hexadecimal
1797 ;; number with a # suffix.
1798 (defun plotcolorp (color)
1799 (cond ((and (stringp color
)
1800 (string= (subseq color
0 1) "#")
1801 (= (length color
) 7)
1802 (ignore-errors (parse-integer (subseq color
1 6) :radix
16)))
1804 ((member color
'($red $green $blue $magenta $cyan $yellow
1805 $orange $violet $brown $gray $black $white
))
1809 ;; tries to convert az into a floating-point number between 0 and 360
1810 (defun parse-azimuth (az) (mod ($float
(meval* az
)) 360))
1812 ;; tries to convert el into a floating-poitn number between -180 and 180
1813 (defun parse-elevation (el) (- (mod (+ 180 ($float
(meval* el
))) 360) 180))
1815 ;; The following functions check the value of an option returning an atom
1816 ;; when there is only one argument or a list when there are several arguments
1819 ;; Checks for one or more items of the same type, using the test given
1820 (defun check-option (option test type
&optional count
)
1822 (unless (= (1- (length option
)) count
)
1825 "Wrong number of arguments for plot option \"~M\". Expecting ~M but found ~M.")
1826 (car option
) count
(1- (length option
)))))
1827 (dolist (item (cdr option
))
1828 (when (not (funcall test item
))
1830 (intl:gettext
"Wrong argument for plot option \"~M\". Expecting ~M but found \"~M\".") (car option
) type item
)))
1831 (if (= (length option
) 2)
1835 ;; Accepts one or more items of the same type or false.
1836 ;; When given, n is the maximum number of items.
1837 (defun check-option-b (option test type
&optional count
)
1838 (let ((n (- (length option
) 1)))
1840 (unless (< n
(1+ count
))
1843 "Plot option ~M must have ~M arguments, not ~M.")
1844 (car option
) count
(1- (length option
)))))
1849 "Option ~M should be given arguments, or called by its name (no lists)")
1852 (if (or (funcall test
(cadr option
)) (null (cadr option
))
1853 (eq (cadr option
) t
))
1857 "Value of option ~M. should be ~M or false, not \"~M\".")
1858 (car option
) type
(cadr option
))))
1861 (unless (funcall test
(nth (+ i
1) option
))
1864 "Value of option ~M should be ~M, not \"~M\".")
1865 (car option
) type
(nth (+ i
1) option
))))
1868 ;; Boolean options can be [option], [option,true] or [option,false]
1869 (defun check-option-boole (option)
1870 (if (= 1 (length option
))
1872 (if (and (= 2 (length option
))
1873 (or (eq (cadr option
) t
) (null (cadr option
))))
1875 (merror (intl:gettext
"plot option ~M must be either true or false.")
1878 ;; label can be either [label, string, real, real] or
1879 ;; [label, [string_1, real, real],...,[string_n, real, real]]
1880 (defun check-option-label (option &aux opt
)
1881 (if (not ($listp
(cadr option
)))
1882 (setq opt
(list (cons '(mlist) (cdr option
))))
1883 (setq opt
(cdr option
)))
1885 (when (not (and ($listp item
) (= 4 (length item
)) (stringp (second item
))
1886 (realp (setf (third item
) (coerce-float (third item
))))
1887 (realp (setf (fourth item
) (coerce-float (fourth item
))))))
1890 "Wrong argument ~M for option ~M. Must be either [label,\"text\",x,y] or [label, [\"text 1\",x1,y1],...,[\"text n\",xn,yn]]")
1891 item
(car option
))))
1894 ;; one of the possible formats
1895 (defun check-option-format (option &aux formats
)
1896 (setq formats
'($geomview $gnuplot $gnuplot_pipes $mgnuplot $xmaxima
))
1897 (unless (member (cadr option
) formats
)
1900 "Wrong argument ~M for option ~M. Must one of the following symbols: geomview, gnuplot, mgnuplot, xmaxima (or gnuplot_pipes in Unix)")
1901 (cadr option
) (car option
)))
1904 ; palette most be one or more Maxima lists starting with the name of one
1905 ;; of the 5 kinds: hue, saturation, value, gray or gradient.
1906 (defun check-option-palette (option)
1907 (if (and (= (length option
) 2) (null (cadr option
)))
1910 (dolist (item (cdr option
))
1911 (when (not (and ($listp item
)
1913 '($hue $saturation $value $gray $gradient
))))
1916 "Wrong argument ~M for option ~M. Not a valid palette.")
1917 item
(car option
))))
1920 ;; style can be one or several of the names of the styles or one or several
1921 ;; Maxima lists starting with the name of one of the styles.
1922 (defun check-option-style (option)
1923 (if (and (= (length option
) 2) (null (cadr option
)))
1927 (dolist (item (cdr option
))
1929 (setq name
(second item
))
1931 (when (not (member name
1932 '($lines $points $linespoints $dots $impulses
)))
1935 "Wrong argument ~M for option ~M. Not a valid style")
1937 (setq parsed
(cons item parsed
)))
1938 (reverse parsed
)))))
1940 ;; Transform can be false or the name of a function for the transformation.
1941 (defun check-option-transform (option)
1942 (if (and (= (length option
) 2)
1943 (or (atom (cadr option
)) (null (cadr option
))))
1947 "Wrong argument ~M for option ~M. Should be either false or the name of function for the transformation") option
(car option
))))
1952 plot2d
(sec(x), [x
, -
2, 2], [y
, -
20, 20]);
1954 plot2d
(exp(3*s
), [s
, -
2, 2], logy
);
1956 plot2d
([parametric
, cos
(t), sin
(t), [t
, -%pi
, %pi
]], same_xy
);
1958 xy
:[[10,.6], [20,.9], [30,1.1], [40,1.3], [50,1.4]]$
1959 plot2d
( [ [discrete
, xy
], 2*%pi
*sqrt
(l/980) ], [l
, 0, 50],
1960 [style
, points
, lines
], [color
, red
, blue
], [point_type
, box
],
1961 [legend
, "experiment", "theory"],
1962 [xlabel
, "pendulum's length (cm)"], [ylabel
, "period (s)"]);
1964 plot2d
( x^
2-
1, [x
, -
3, 3], [y
, -
2, 10], nobox
, [color
, red
],
1965 [ylabel
, "x^2-1"], [plot_format
, xmaxima
]);
1967 plot2d
( x^
2+y^
2 = 1, [x
, -
2, 2], [y
, -
2 ,2]);
1970 (fun &optional xrange
&rest extra-options
1972 ($display2d nil
) (*plot-realpart
* *plot-realpart
*)
1973 (options (copy-tree *plot-options
*)) yrange output-file plot
)
1974 ;; fun must be a maxima list with several objects: expressions (simple
1975 ;; functions), maxima lists (parametric or discrete cases).
1976 ;; A single parametric or discrete plot is placed inside a maxima list.
1977 (setf (getf options
:type
) "plot2d")
1978 (when (and (consp fun
)
1979 (or (eq (second fun
) '$parametric
) (eq (second fun
) '$discrete
)))
1980 (setq fun
`((mlist) ,fun
)))
1981 ;; If by now fun is not a maxima list, it is then a single expression
1982 (unless ($listp fun
) (setq fun
`((mlist) ,fun
)))
1983 ;; 2- Get names for the two axis and values for xmin and xmax if needed.
1984 ;; If any of the objects in the fun list is a simple function,
1985 ;; the xrange option is mandatory and will provide the name of
1986 ;; the horizontal axis and the values of xmin and xmax.
1987 (let ((xrange-required nil
) (bounds-required nil
) (yrange-required nil
)
1988 small huge fpfun vars1 vars2 prange
)
1989 #-clisp
(setq small
(- (/ most-positive-flonum
1024)))
1990 #+clisp
(setq small
(- (/ most-positive-double-float
1024.0)))
1991 #-clisp
(setq huge
(/ most-positive-flonum
1024))
1992 #+clisp
(setq huge
(/ most-positive-double-float
1024.0))
1993 (setf (getf options
:ybounds
) (list small huge
))
1994 (dolist (f (rest fun
))
1999 (unless bounds-required
2000 (setq bounds-required t
)
2001 ;; Default X and Y bound large so parametric plots don't get
2002 ;; prematurely clipped. Don't use most-positive-flonum
2003 ;; because draw2d will overflow.
2004 (setf (getf options
:xbounds
) (list small huge
)))
2005 (setq prange
(check-range ($fourth f
)))
2006 ;; The two expressions can only depend on the parameter given
2007 (setq fpfun
(coerce-float-fun ($second f
) ($rest prange -
2)))
2008 (setq vars1
($listofvars
(mfuncall fpfun
($first prange
))))
2009 (setq fpfun
(coerce-float-fun ($third f
) ($rest prange -
2)))
2010 (setq vars2
($listofvars
(mfuncall fpfun
($first prange
))))
2011 (setq vars1
($listofvars
`((mlist) ,vars1
,vars2
)))
2012 (setq vars1
(delete ($first prange
) vars1
))
2013 (when (> ($length vars1
) 0)
2016 "plot2d: parametric expressions ~M and ~M should depend only on ~M")
2017 ($second f
) ($third f
) ($first prange
))))
2022 "plot2d: a keyword 'parametric' or 'discrete' missing in ~M")
2024 ;; The expression represents a function, explicit or implicit
2026 (unless xrange-required
2027 (setq xrange-required t
)
2028 (setq xrange
(check-range xrange
))
2029 (setq xrange-required t
)
2030 (unless (getf options
:xlabel
)
2031 (setf (getf options
:xlabel
) (ensure-string (second xrange
))))
2032 (setf (getf options
:xvar
) (cadr xrange
))
2033 (setf (getf options
:x
) (cddr xrange
)))
2034 (if (and (listp f
) (member 'mequal
(flatten f
)))
2036 ;; Implicit function
2039 (coerce-float-fun (m- ($lhs f
) ($rhs f
)) ($rest xrange -
2)))
2040 (setq vars1
($listofvars
(mfuncall fpfun
($first xrange
))))
2043 (= ($length vars1
) 2)
2044 (not (member ($first xrange
) vars1
)))
2047 "plot2d: ~M is not one of the variables in ~M")
2049 (setq vars1
(delete ($first xrange
) vars1
))
2050 (if (< ($length vars1
) 2)
2054 (or (= ($length vars1
) 0)
2055 (eq ($first yrange
) ($first vars1
)))
2058 "plot2d: ~M should only depend on ~M and ~M")
2059 f
($first xrange
) ($first vars1
)))
2061 (setq yrange-required t
)
2062 (if (null extra-options
)
2065 "plot2d: Missing interval for variable 2."))
2067 (setq yrange
(pop extra-options
))
2068 (setq vars1
(delete ($first yrange
) vars1
))
2069 (unless (= ($length vars1
) 0)
2072 "plot2d: ~M should only depend on ~M and ~M")
2073 f
($first xrange
) ($first yrange
)))
2074 (setq yrange
(check-range yrange
))
2075 (setf (getf options
:yvar
) ($first yrange
))
2076 (setf (getf options
:y
) (cddr yrange
)))))))
2079 "plot2d: ~M should only depend on 2 variables")
2082 ;; Explicit function
2083 (setq fpfun
(coerce-float-fun f
($rest xrange -
2)))
2084 (setq vars1
($listofvars
(mfuncall fpfun
($first xrange
))))
2085 (setq vars1
(delete ($first xrange
) vars1
))
2086 (when (> ($length vars1
) 0)
2089 "plot2d: expression ~M should depend only on ~M")
2090 f
($first xrange
))))))))
2091 (when (not xrange-required
)
2092 ;; Make the default ranges on X nd Y large so parametric plots
2093 ;; don't get prematurely clipped. Don't use most-positive-flonum
2094 ;; because draw2d will overflow.
2095 (setf (getf options
:xbounds
) (list small huge
))
2097 ;; second argument was really a plot option, not an xrange
2098 (setq extra-options
(cons xrange extra-options
)))))
2099 ;; If no global options xlabel or ylabel have been given, choose
2100 ;; a default value for them: the expressions given, converted
2101 ;; to Maxima strings, if their length is less than 50 characters,
2102 ;; or the default "x" and "y" otherwise.
2103 (when (= (length fun
) 2)
2104 (let ((v (second fun
)) xlabel ylabel
)
2106 (setq xlabel
"x") (setq ylabel
($sconcat v
)))
2107 ((eq (second v
) '$parametric
)
2108 (setq xlabel
($sconcat
(third v
)))
2109 (setq ylabel
($sconcat
(fourth v
))))
2110 ((eq (second v
) '$discrete
)
2111 (setq xlabel
"x") (setq ylabel
"y"))
2113 (setq xlabel
"x") (setq ylabel
($sconcat v
))))
2114 (unless (getf options
:xlabel
)
2115 (if (< (length xlabel
) 50) (setf (getf options
:xlabel
) xlabel
)))
2116 (unless (getf options
:ylabel
)
2117 (if (< (length ylabel
) 50) (setf (getf options
:ylabel
) ylabel
)))))
2118 ;; Parse the given options into the options list
2119 (setq options
(plot-options-parser extra-options options
))
2120 (when (getf options
:y
) (setf (getf options
:ybounds
) (getf options
:y
)))
2121 ;; Remove axes labels when no box is used in gnuplot
2122 (when (and (member :box options
) (not (getf options
:box
))
2123 (not (eq (getf options
:plot_format
) '$xmaxima
)))
2124 (remf options
:xlabel
)
2125 (remf options
:ylabel
))
2126 ;; check options given
2127 (let ((xmin (first (getf options
:x
))) (xmax (second (getf options
:x
))))
2129 (and (getf options
:logx
) xmin xmax
)
2132 (let ((revised-xmin (/ xmax
1000)))
2135 "plot2d: lower bound must be positive when using 'logx'.~%plot2d: assuming lower bound = ~M instead of ~M")
2137 (setf (getf options
:x
) (list revised-xmin xmax
))
2138 (setq xrange
`((mlist) ,(second xrange
) ,revised-xmin
,xmax
))))
2141 "plot2d: upper bound must be positive when using 'logx'; found: ~M")
2143 (let ((ymin (first (getf options
:y
)))
2144 (ymax (second (getf options
:y
))))
2145 (when (and (getf options
:logy
) ymin ymax
)
2148 (let ((revised-ymin (/ ymax
1000)))
2151 "plot2d: lower bound must be positive when using 'logy'.~%plot2d: assuming lower bound = ~M instead of ~M")
2153 (setf (getf options
:y
) (list revised-ymin ymax
))))
2156 "plot2d: upper bound must be positive when using 'logy'; found: ~M")
2158 (setq *plot-realpart
* (getf options
:plot_realpart
))
2159 ;; Creates the object that will be passed to the external graphic program
2160 (case (getf options
:plot_format
)
2162 (setq plot
(make-instance 'xmaxima-plot
)))
2164 (setq plot
(make-instance 'gnuplot-plot
)))
2166 (setq plot
(make-instance 'gnuplot-plot
))
2167 (setf (slot-value plot
'pipe
) T
))
2169 (merror (intl:gettext
"plot2d: plot format ~M not supported")
2170 (getf options
:plot_format
))))
2171 ;; Parse plot object and pass it to the graphic program
2172 (setq output-file
(plot-preamble plot options
))
2173 (plot2d-command plot fun options xrange
)
2174 (plot-shipout plot options output-file
))
2177 (and (symbolp x
) (char= (char (symbol-value x
) 0) #\$
)))
2180 (defmfun $tcl_output
(lis i
&optional
(skip 2))
2181 (when (not (typep i
'fixnum
))
2183 (intl:gettext
"tcl_ouput: second argument must be an integer; found ~M")
2185 (when (not ($listp lis
))
2187 (intl:gettext
"tcl_output: first argument must be a list; found ~M") lis
))
2188 (format *standard-output
* "~% {")
2189 (cond (($listp
(second lis
))
2192 (format *standard-output
* "~,10g " (nth i v
))))
2194 (setq lis
(nthcdr i lis
))
2195 (loop with v
= lis while v
2197 (format *standard-output
* "~,10g " (car v
))
2198 (setq v
(nthcdr skip v
)))))
2199 (format *standard-output
* "~% }"))
2200 (defun tcl-output-list ( st lis
)
2206 when
(eql 0 (mod n
5))
2209 (format st
"~,10g " v
))
2211 (t (tcl-output-list st
(car lis
))
2212 (tcl-output-list st
(cdr lis
)))))
2214 (defun check-range (range &aux tem a b
)
2215 (or (and ($listp range
)
2216 (setq tem
(cdr range
))
2217 (or (symbolp (car tem
)) ($subvarp
(car tem
)))
2218 (numberp (setq a
($float
(meval* (second tem
)))))
2219 (numberp (setq b
($float
(meval* (third tem
)))))
2223 (intl:gettext
"plotting: range must be of the form [variable, min, max]; found: ~M")
2226 (intl:gettext
"plotting: no range given; must supply range of the form [variable, min, max]"))))
2227 `((mlist) ,(car tem
) ,(float a
) ,(float b
)))
2229 (defmfun $zero_fun
(x y
) x y
0.0)
2231 (defun output-points (pl &optional m
)
2232 "If m is supplied print blank line every m lines"
2234 (declare (fixnum j
))
2235 (loop for i below
(length (polygon-pts pl
))
2236 with ar
= (polygon-pts pl
)
2237 do
(print-pt (aref ar i
))
2239 (print-pt (aref ar i
))
2241 (print-pt (aref ar i
))
2245 (cond ((eql j
(the fixnum m
))
2250 (defun output-points-tcl (dest pl m
)
2251 (format dest
" {matrix_mesh ~%")
2252 ;; x y z are done separately:
2253 (loop for off from
0 to
2
2254 with ar
= (polygon-pts pl
)
2255 with i of-type fixnum
= 0
2259 while
(< i
(length ar
))
2260 do
(format dest
"~% {")
2262 do
(print-pt (aref ar i
))
2264 (format dest
"}~%"))
2265 (format dest
"}~%"))
2266 (format dest
"}~%"))
2268 (defun show-open-plot (ans file
)
2269 (cond ($show_openplot
2270 (with-open-file (st1 (plot-temp-file (format nil
"maxout~d.xmaxima" (getpid))) :direction
:output
:if-exists
:supersede
)
2272 ($system
(concatenate 'string
*maxima-prefix
*
2273 (if (string= *autoconf-windows
* "true") "\\bin\\" "/bin/")
2274 $xmaxima_plot_command
)
2275 #-
(or (and sbcl win32
) (and sbcl win64
) (and ccl windows
))
2276 (format nil
" ~s &" file
)
2277 #+(or (and sbcl win32
) (and sbcl win64
) (and ccl windows
))
2279 (t (princ ans
) "")))
2282 ;; contour_plot -- set some parameters for Gnuplot and punt to plot3d
2284 ;; We go to some trouble here to avoid clobbering the Gnuplot preamble
2285 ;; specified by the user, either as a global option (via set_plot_option)
2286 ;; or specified in arguments to contour_plot. Just append or prepend
2287 ;; the parameters for contour plotting to the user-specified preamble.
2288 ;; Assume that arguments take precedence over global options.
2290 ;; contour_plot knows how to set parameters only for Gnuplot.
2291 ;; If the plot_format is not a Gnuplot format, complain.
2295 ;; contour_plot (x^2 + y^2, [x, -4, 4], [y, -4, 4]);
2296 ;; contour_plot (sin(y) * cos(x)^2, [x, -4, 4], [y, -4, 4]);
2297 ;; F(x, y) := x^3 + y^2;
2298 ;; contour_plot (F, [u, -4, 4], [v, -4, 4]);
2299 ;; contour_plot (F, [u, -4, 4], [v, -4, 4], [gnuplot_preamble, "set size ratio -1"]);
2300 ;; set_plot_option ([gnuplot_preamble, "set cntrparam levels 12"]);
2301 ;; contour_plot (F, [u, -4, 4], [v, -4, 4]);
2302 ;; set_plot_option ([plot_format, xmaxima]);
2303 ;; contour_plot (F, [u, -4, 4], [v, -4, 4]); => error: must be gnuplot format
2304 ;; contour_plot (F, [u, -4, 4], [v, -4, 4], [plot_format, gnuplot]);
2306 (defmfun $contour_plot
(expr &rest optional-args
)
2308 ((plot-format-in-options (getf *plot-options
* :plot_format
))
2309 (plot-format-in-arguments
2310 (caddar (member '$plot_format optional-args
:key
#'cadr
)))
2311 (preamble-in-plot-options (getf *plot-options
* :gnuplot_preamble
))
2312 (preamble-in-arguments
2313 (caddar (member '$gnuplot_preamble optional-args
:key
#'cadr
)))
2314 (contour-preamble "set contour; unset surface; set view map")
2315 (gnuplot-formats '($gnuplot $gnuplot_pipes
))
2317 ;; Ensure that plot_format is some gnuplot format.
2318 ;; Argument takes precedence over global option.
2321 (and plot-format-in-arguments
2322 (not (member plot-format-in-arguments gnuplot-formats
:test
#'eq
)))
2323 (and (not plot-format-in-arguments
)
2324 (not (member plot-format-in-options gnuplot-formats
:test
#'eq
))))
2326 (mtell(intl:gettext
"contour_plot: plot_format = ~a not recognized; must be a gnuplot format.~%")
2327 (ensure-string (or plot-format-in-arguments plot-format-in-options
)))
2328 (return-from $contour_plot
)))
2330 ;; Prepend contour preamble to preamble in arguments (if given)
2331 ;; and pass concatenated preamble as an argument to plot3d.
2332 ;; Otherwise if there is a global option preamble,
2333 ;; append contour preamble to global option preamble.
2334 ;; Otherwise just set global option preamble to the contour preamble.
2336 ;; All this complication is to avoid clobbering the preamble
2337 ;; if one was specified somehow (either global option or argument).
2339 (if preamble-in-arguments
2342 (remove-if #'(lambda (e) (and ($listp e
) (eq ($first e
) '$gnuplot_preamble
))) optional-args
))
2344 ($sconcat contour-preamble
(format nil
"~%")
2345 preamble-in-arguments
)))
2346 (if preamble-in-plot-options
2348 ($sconcat preamble-in-plot-options
(format nil
"~%")
2350 (setq preamble contour-preamble
)))
2354 (list '((mlist) $palette nil
))
2355 (list `((mlist) $gnuplot_preamble
,preamble
))))))
2360 plot3d
(2^
(-u^
2 + v^
2), [u
, -
3, 3], [v
, -
2, 2], [palette
, false
]);
2362 plot3d
( log
( x^
2*y^
2 ), [x
, -
2, 2], [y
, -
2, 2], [grid
, 29, 29]);
2364 expr_1
: cos
(y)*(10.0
+6*cos
(x))$
2365 expr_2
: sin
(y)*(10.0
+6*cos
(x))$
2367 plot3d
([expr_1
, expr_2
, expr_3
], [x
, 0, 2*%pi
], [y
, 0, 2*%pi
],
2368 ['grid
, 40, 40], [z
,-
8,8]);
2370 plot3d
(cos (-x^
2 + y^
3/4), [x
, -
4, 4], [y
, -
4, 4],
2371 [mesh_lines_color
, false
], [elevation
, 0], [azimuth
, 0], [grid
, 150, 150]);
2373 spherical
: make_transform
([th
, phi
,r
], r
*sin
(phi)*cos
(th),
2374 r
*sin
(phi)*sin
(th), r
*cos
(phi))$
2375 plot3d
( 5, [th
, 0, 2*%pi
], [phi
, 0, %pi
], [transform_xy
, spherical
],
2376 [palette
, [value
, 0.65, 0.7, 0.1, 0.9]], [plot_format
,xmaxima
]);
2378 V
: 1 / sqrt
( (x+1)^
2+y^
2 ) -
1 / sqrt
( (x-1)^
2+y^
2 )$
2379 plot3d
( V
, [x
, -
2, 2], [y
, -
2, 2], [z
, -
4, 4]);
2382 (fun &rest extra-options
2384 lvars xrange yrange titles output-file functions exprn domain tem
2385 (options (copy-tree *plot-options
*)) (*plot-realpart
* *plot-realpart
*)
2386 (usage (intl:gettext
2388 To plot a single function f of 2 variables v1 and v2:
2389 plot3d (f, [v1, min, max], [v2, min, max], options)
2390 A parametric representation of a surface with parameters v1 and v2:
2391 plot3d ([f1, f2, f3], [v1, min, max], [v2, min, max], options)
2392 Several functions depending on the two variables v1 and v2:
2393 plot3d ([f1, f2, ..., fn, [v1, min, max], [v2, min, max]], options)")))
2394 (setf (getf options
:type
) "plot3d")
2395 ;; Ensure that fun is a list of expressions and maxima lists, followed
2396 ;; by a domain definition
2398 (if (= 1 (length (check-list-plot3d fun
)))
2399 ;; fun consisted of a single parametric expression
2400 (setq fun
`(,fun
,(pop extra-options
) ,(pop extra-options
)))
2401 ;; fun was a maxima list with several independent surfaces
2403 ;; fun consisted of a single expression
2404 (setq fun
`(,fun
,(pop extra-options
) ,(pop extra-options
))))
2405 ;; go through all the independent surfaces creating the functions stack
2407 (setq exprn
(pop fun
))
2410 (setq domain
(check-list-plot3d exprn
))
2411 (case (length domain
)
2413 ;; exprn is a parametric representation of a surface
2414 (let (vars1 vars2 vars3
)
2415 ;; list fun should have two valid ranges after exprn
2416 (setq xrange
(check-range (pop fun
)))
2417 (setq yrange
(check-range (pop fun
)))
2418 ;; list of the two variables for the parametric equations
2419 (setq lvars
`((mlist),(second xrange
) ,(second yrange
)))
2420 ;; make sure that the 3 parametric equations depend only
2421 ;; on the two variables in lvars
2423 ($listofvars
(mfuncall
2424 (coerce-float-fun (second exprn
) lvars
)
2425 (second lvars
) (third lvars
))))
2427 ($listofvars
(mfuncall
2428 (coerce-float-fun (third exprn
) lvars
)
2429 (second lvars
) (third lvars
))))
2431 ($listofvars
(mfuncall
2432 (coerce-float-fun (fourth exprn
) lvars
)
2433 (second lvars
) (third lvars
))))
2434 (setq lvars
($listofvars
`((mlist) ,vars1
,vars2
,vars3
)))
2435 (if (<= ($length lvars
) 2)
2436 ;; we do have a valid parametric set. Push it into
2437 ;; the functions stack, along with their domain
2439 (push `(,exprn
,xrange
,yrange
) functions
)
2440 ;; add a title to the titles stack
2441 (push "Parametric function" titles
)
2442 ;; unknown variables in the parametric equations
2443 ;; ----- GNUPLOT 4.0 WORK-AROUND -----
2444 (when (and ($constantp
(fourth exprn
))
2445 (getf options
:gnuplot_4_0
))
2446 (setf (getf options
:const_expr
)
2447 ($float
(meval (fourth exprn
))))))
2449 (intl:gettext
"plot3d: there must be at most two variables; found: ~M")
2452 ;; expr is a simple function with its own domain. Push the
2453 ;; function and its domain into the functions stack
2454 (setq xrange
(second domain
))
2455 (setq yrange
(third domain
))
2456 (push `(,(second exprn
) ,xrange
,yrange
) functions
)
2457 ;; push a title for this plot into the titles stack
2458 (if (< (length (ensure-string (second exprn
))) 36)
2459 (push (ensure-string (second exprn
)) titles
)
2460 (push "Function" titles
)))
2462 ;; syntax error. exprn does not have the expected form
2465 "plot3d: argument must be a list of three expressions; found: ~M")
2468 ;; exprn is a simple function, defined in the global domain.
2469 (if (and (getf options
:xvar
) (getf options
:yvar
))
2470 ;; the global domain has already been defined; use it.
2472 (setq xrange
`((mlist) ,(getf options
:xvar
)
2473 ,(first (getf options
:x
))
2474 ,(second (getf options
:x
))))
2475 (setq yrange
`((mlist) ,(getf options
:yvar
)
2476 ,(first (getf options
:y
))
2477 ,(second (getf options
:y
)))))
2478 ;; the global domain should be defined by the last two lists
2479 ;; in fun. Extract it and check whether it is valid.
2483 (check-list-plot3d (append `((mlist) ,exprn
) (last fun
2))))
2484 (setq fun
(butlast fun
2))
2485 (if (= 3 (length domain
))
2486 ;; domain is valid. Make it the global one.
2488 (setq xrange
(second domain
))
2489 (setq yrange
(third domain
))
2490 (setf (getf options
:xvar
) (second xrange
))
2491 (setf (getf options
:x
) (cddr xrange
))
2492 (setf (getf options
:yvar
) (second yrange
))
2493 (setf (getf options
:y
) (cddr yrange
)))
2495 ;; ----- GNUPLOT 4.0 WORK-AROUND -----
2496 (when (and ($constantp exprn
) (getf options
:$gnuplot_4_0
))
2497 (setf (getf options
:const_expr
) ($float
(meval exprn
))))
2498 ;; push the function and its domain into the functions stack
2499 (push `(,exprn
,xrange
,yrange
) functions
)
2500 ;; push a title for this plot into the titles stack
2501 (if (< (length (ensure-string exprn
)) 36)
2502 (push (ensure-string exprn
) titles
)
2503 (push "Function" titles
))))
2504 (when (= 0 (length fun
)) (return)))
2505 ;; recover the original ordering for the functions and titles stacks
2506 (setq functions
(reverse functions
))
2507 (setq titles
(reverse titles
))
2508 ;; parse the options given to plot3d
2509 (setq options
(plot-options-parser extra-options options
))
2510 (setq tem
(getf options
:transform_xy
))
2511 (when (and (member :gnuplot_pm3d options
) (null (getf options
:gnuplot_pm3d
)))
2512 (setf (getf options
:palette
) nil
))
2513 (setq *plot-realpart
* (getf options
:plot_realpart
))
2514 ;; set up the labels for the axes, unless no box is being shown
2515 (unless (and (member :box options
) (not (getf options
:box
)))
2516 (if (and (getf options
:xvar
) (getf options
:yvar
) (null tem
))
2518 ;; Don't set xlabel (ylabel) if the user specified one.
2519 (unless (getf options
:xlabel
)
2520 (setf (getf options
:xlabel
) (ensure-string (getf options
:xvar
))))
2521 (unless (getf options
:ylabel
)
2522 (setf (getf options
:ylabel
) (ensure-string (getf options
:yvar
)))))
2524 (setf (getf options
:xlabel
) "x")
2525 (setf (getf options
:ylabel
) "y")))
2526 (unless (getf options
:zlabel
) (setf (getf options
:zlabel
) "z")))
2527 ;; x and y should not be bound, when an xy transformation function is used
2528 (when tem
(remf options
:x
) (remf options
:y
))
2529 ;; Set up the plot command
2530 (let (plot (legend (getf options
:legend
)))
2531 ;; titles will be a list. Titles given in the legend option prevail
2532 ;; over titles generated by plot3d. No legend if option [legend,false]
2533 (unless (listp legend
) (setq legend
(list legend
)))
2534 (when (member :legend options
)
2535 (if (first legend
) (setq titles legend
)) (setq titles nil
))
2536 (case (getf options
:plot_format
)
2538 (setq plot
(make-instance 'xmaxima-plot
)))
2540 (setq plot
(make-instance 'gnuplot-plot
)))
2542 (setq plot
(make-instance 'gnuplot-plot
))
2543 (setf (slot-value plot
'pipe
) T
))
2545 (setq plot
(make-instance 'geomview-plot
)))
2547 (merror (intl:gettext
"plot3d: plot format ~M not supported")
2548 (getf options
:plot_format
))))
2549 ;; Parse plot object and pass it to the graphic program
2550 (setq output-file
(plot-preamble plot options
))
2551 (plot3d-command plot functions options titles
)
2552 (plot-shipout plot options output-file
)))
2554 ;; Given a Maxima list with 3 elements, checks whether it represents a function
2555 ;; defined in a 2-dimensional domain or a parametric representation of a
2556 ;; 3-dimensional surface, depending on two parameters.
2557 ;; The return value will be a Maxima list if the test is succesfull or nil
2559 ;; In the case of a function and a 2D domain, it returns the domain, validated.
2560 ;; When it is a parametric representation it returns an empty Maxima list.
2562 (defun check-list-plot3d (lis)
2563 (let (xrange yrange
)
2564 ;; Makes sure list has the form ((mlist) $atom item1 item2)
2565 (unless (and ($listp lis
) (= 3 ($length lis
)) (not ($listp
(second lis
))))
2566 (return-from check-list-plot3d nil
))
2567 ;; we might have a function with domain or a parametric representation
2568 (if ($listp
(third lis
))
2569 ;; lis is probably a function with a valid domain
2570 (if ($listp
(fourth lis
))
2571 ;; we do have a function and a domain. Return the domain
2573 (setq xrange
(check-range (third lis
)))
2574 (setq yrange
(check-range (fourth lis
)))
2575 (return-from check-list-plot3d
`((mlist) ,xrange
,yrange
)))
2576 ;; wrong syntax: [expr1, list, expr2]
2577 (return-from check-list-plot3d nil
))
2578 ;; lis is probably a parametric representation
2579 (if ($listp
(fourth lis
))
2580 ;; wrong syntax: [expr1, expr2, list]
2581 (return-from check-list-plot3d nil
)
2582 ;; we do have a parametric representation. Return an empty list
2583 (return-from check-list-plot3d
'((mlist)))))))
