1
;;; -*- Mode: Lisp -*- ;;;;
2
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
4
;; Author: Andrej Vodopivec <andrejv@users.sourceforge.net> ;;
8
;; implicit_plot(expr, xrange, yrange, [options]); ;;
9
;; Plots the curve `expr' in the region given by `xrange' and `yrange'. ;;
10
;; `expr' is a curve in plane defined in implicit form or a list of such ;;
11
;; curves. If `expr' is not an equality, then `expr=0' is assumed. Works ;;
12
;; by tracking sign changes, so it will fail if expr is something like ;;
14
;; Optional argument `options' can be anything that is recognized by ;;
15
;; `plot2d'. Options can also be set using `set_plot_option'. ;;
16
;; Works only with gnuplot! ;;
19
;; implicit_plot(y^2=x^3-2*x+1, [x,-4,4], [y,-4,4], ;;
20
;; [gnuplot_preamble, "set zeroaxis"])$ ;;
21
;; implicit_plot([x^2-y^2/9=1,x^2/4+y^2/9=1], [x,-2.5,2.5], [x,-3.5,3.5]); ;;
22
;; implicit_plot(x^2+2*y^3=15, [x,-10, 10], [y,-5,5])$ ;;
23
;; implicit_plot(x^2*y^2=(y+1)^2*(4-y^2), [x,-10, 10], [y,-3,3]); ;;
24
;; implicit_plot(x^3+y^3 = 3*x*y^2-x-1, [x,-4,4], [y,-4,4]); ;;
25
;; implicit_plot(x^2*sin(x+y)+y^2*cos(x-y)=1, [x,-10,10], [y,-10,10]); ;;
27
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
30
(macsyma-module implicit_plot)
32
(defmvar $ip_grid `((mlist simp) 75 75)
33
"Grid for the first sampling.")
34
(defmvar $ip_grid_in `((mlist simp) 10 10)
35
"Grid for the second sampling.")
37
(defmvar $ip_epsilon 0.0000001
38
"Epsilon for implicit plot routine.")
40
(defun contains-zeros (i j sample)
41
(not (and (> (* (aref sample i j) (aref sample (1+ i) j )) 0)
42
(> (* (aref sample i j) (aref sample i (1+ j) )) 0)
43
(> (* (aref sample i j) (aref sample (1+ i) (1+ j) )) 0) )))
45
(defun sample-data (expr xmin xmax ymin ymax sample grid)
46
(let* ((xdelta (/ (- xmax xmin) ($first grid)))
47
(ydelta (/ (- ymax ymin) ($second grid))))
48
(do ((x-val xmin (+ x-val xdelta))
51
(do ((y-val ymin (+ y-val ydelta))
53
((> j ($second grid)))
54
(let ((fun-val (funcall expr x-val y-val)))
55
(if (or (eq fun-val t) (>= fun-val $ip_epsilon))
56
(setf (aref sample i j) 1)
57
(setf (aref sample i j) -1)))))))
59
(defun print-mid (file point xmin xdelta ymin ydelta)
60
(let ((x (+ xmin (/ (* xdelta (+ (car point) (caddr point))) 2)))
61
(y (+ ymin (/ (* ydelta (+ (cadr point) (cadddr point))) 2))))
62
(format file "~f ~f~%" x y)))
64
(defun print-square (file xmin xmax ymin ymax sample grid)
65
(let* ((xdelta (/ (- xmax xmin) ($first grid)))
66
(ydelta (/ (- ymax ymin) ($second grid))))
68
((>= (+ i 2) ($first grid)))
70
((>= (+ j 2) ($second grid)))
71
(if (contains-zeros i j sample)
73
(if (< (* (aref sample i j) (aref sample (1+ i) j)) 0)
74
(setq points (cons `(,i ,j ,(1+ i) ,j) points)))
75
(if (< (* (aref sample (1+ i) j) (aref sample (1+ i) (1+ j))) 0)
76
(setq points (cons `(,(1+ i) ,j ,(1+ i) ,(1+ j)) points)))
77
(if (< (* (aref sample i (1+ j)) (aref sample (1+ i) (1+ j))) 0)
78
(setq points (cons `(,i ,(1+ j) ,(1+ i) ,(1+ j)) points)))
79
(if (< (* (aref sample i j) (aref sample i (1+ j))) 0)
80
(setq points (cons `(,i ,j ,i ,(1+ j)) points)))
81
(print-mid file (car points) xmin xdelta ymin ydelta)
82
(print-mid file (cadr points) xmin xdelta ymin ydelta)
83
(format file "~%")))))
86
(defun $implicit_plot (expr xrange yrange &rest options)
87
(let* (($numer t) ($plot_options $plot_options)
90
(xmin ($second xrange))
91
(xmax ($third xrange))
92
(xdelta (/ (- xmax xmin) ($first $ip_grid)))
93
(ymin ($second yrange))
94
(ymax ($third yrange))
95
(ydelta (/ (- ymax ymin) ($second $ip_grid)))
96
(sample (make-array `(,(1+ ($first $ip_grid))
97
,(1+ ($second $ip_grid)))))
98
(ssample (make-array `(,(1+ ($first $ip_grid_in))
99
,(1+ ($second $ip_grid_in)))))
100
file-name gnuplot-out-file gnuplot-term)
102
(dolist (v options) ($set_plot_option v))
103
(setq xrange (check-range xrange))
104
(setq yrange (check-range yrange))
106
(if (not ($listp expr))
107
(setq expr `((mlist simp) ,expr)))
109
(setf gnuplot-term ($get_plot_option '$gnuplot_term 2))
111
(if ($get_plot_option '$gnuplot_out_file 2)
112
(setf gnuplot-out-file (get-plot-option-string '$gnuplot_out_file)))
114
(if (and (eq gnuplot-term '$default)
116
(setf file-name gnuplot-out-file)
117
(setf file-name (plot-temp-file "maxout.gnuplot")))
121
(file file-name :direction :output :if-exists :supersede)
122
(gnuplot-print-header file)
123
(format file "set xrange [~d:~d]~%" (caddr xrange) (cadddr xrange))
124
(format file "set yrange [~d:~d]~%" (caddr yrange) (cadddr yrange))
125
(format file "set style data lines~%")
127
(dolist (v (cdr expr))
132
(setf string (coerce (mstring v) 'string))
133
(setf string (coerce (mstring v) 'string)))
134
(if (< (length string) 80)
136
(format nil "fun~a" i))))
139
(let ((title (get-plot-option-string '$gnuplot_curve_titles i)))
140
(if (equal title "default")
141
(setf title (format nil "title '~a'" plot-name)))
142
(format file " '-' ~a ~a" title
143
(get-plot-option-string '$gnuplot_curve_styles i))))
145
(dolist (e (cdr expr))
146
(setq e (coerce-float-fun ($float (m- ($rhs e) ($lhs e)))
150
(sample-data e xmin xmax ymin ymax sample $ip_grid)
152
((= (1+ i) ($first $ip_grid)))
154
((= (1+ j) ($second $ip_grid)))
155
(if (contains-zeros i j sample)
156
(let* ((xxmin (+ xmin (* i xdelta)))
157
(xxmax (+ xxmin xdelta xdelta))
158
(yymin (+ ymin (* j ydelta)))
159
(yymax (+ yymin ydelta ydelta)))
160
(sample-data e xxmin xxmax yymin yymax
162
(print-square file xxmin xxmax yymin yymax
163
ssample $ip_grid_in) )) ))
164
(format file "e~%") ))
167
(gnuplot-process file-name) ))
added
removed
share/contrib/implicit_plot.lisp
