Merge branch 'bug-4370-eigevalues-xref-dgeev'
[maxima.git] / tests / rtest_plot.mac
blobe36a6108c80eadde254e3c0a336fe45089f6f319
1 /* execute this script 
2  *
3  * (1) wait for user between examples:
4  *  demo ("tests/rtest_plot.mac");
5  *
6  * (2) or noninteractively:
7  *  batch ("tests/rtest_plot.mac");
8  *
9  * Maxima waits for the user to close the Openmath window
10  * before proceeding, but otherwise it just blazes through.
11  */
13 /* Examples copied from plot2d documentation */
15 "Regression for f0 issue"$
16 plot2d(x^2,[x,-1,1],[yx_ratio,1]);
17 plot2d(x^2,[x,-1,1],[xy_scale,1,1]);
18 plot3d(x^2,[x,-1,1],[y,-1,1],[zmin,1]);
19 plot3d(x^2,[x,-1,1],[y,-1,1],[xtics,1]);
20 plot3d(x^2,[x,-1,1],[y,-1,1],[ytics,1]);
21 plot3d(x^2,[x,-1,1],[y,-1,1],[ztics,1]);
22 plot3d(x^2,[x,-1,1],[y,-1,1],[xtics,1,1]);
23 plot3d(x^2,[x,-1,1],[y,-1,1],[ytics,1,1]);
24 plot3d(x^2,[x,-1,1],[y,-1,1],[ztics,1,1]);
25 plot3d(x^2,[x,-1,1],[y,-1,1],[color_bar_tics,1]);
26 plot3d(x^2,[x,-1,1],[y,-1,1],[label,"tako",1,1]);
27 plot3d(x^2,[x,-1,1],[y,-1,1],[label,["tako",1,1],["surume",1,1],["ika",1,1]]);
29 "Plots of common functions." $
31 plot2d (sin(x), [x, -5, 5])$
32 plot2d (sec(x), [x, -2, 2], [y, -20, 20], [nticks, 200])$
34 "Plotting functions by name." $
36 F(x) := x^2 $
37 :lisp (defun |$g| (x) (m* x x x))
38 H(x) := if x < 0 then x^4 - 1 else 1 - x^5 $
39 plot2d (F, [u, -1, 1])$
40 plot2d ([F, G, H], [u, -1, 1], [y, -1.5, 1.5])$
41 kill(F,H,G) $
43 "We can plot a circle using a parametric plot ..." $
45 plot2d ([parametric, cos(t), sin(t), [t,-%pi,%pi]], [nticks,80], [x, -4/3, 4/3])$
47 "If we repeat that plot with only 8 points ..." $
48 "As parametric plot uses adaptive method, nticks=4 and adapt_depth=1 is for 8 point drawing" $
50 plot2d ([parametric, cos(t), sin(t), [t, -%pi*2, %pi*2]], [nticks, 4], [adapt_depth,1], [x, -2, 2], [y, -1.5, 1.5])$
52 "Combination of an ordinary plot ..." $
54 plot2d ([x^3+2, [parametric, cos(t), sin(t), [t, -5, 5]]], [x, -3, 3],  [nticks, 80])$
56 "Example of a logarithmic plot:" $
58 plot2d (exp(3*s), [s, -2, 2], [logy])$
60 "To show some examples of discrete plots, ..." $
62 xx:[10, 20, 30, 40, 50]$
63 yy:[.6, .9, 1.1, 1.3, 1.4]$
64 xy:[[10,.6], [20,.9], [30,1.1], [40,1.3], [50,1.4]]$
65 plot2d([discrete,xx,yy])$
67 "We will now show the plot with only points, ..." $
69 plot2d([discrete, xy], [style, points])$
71 "The plot of the data points can be shown together with ..." $
73 plot2d([[discrete,xy], 2*%pi*sqrt(l/980)], [l,0,50],
74  [style, [points,5,2,6], [lines,1,1]], [legend,"experiment","theory"],
75  [xlabel,"pendulum's length (cm)"], [ylabel,"period (s)"])$
77 "To save a plot ..." $
79 plot2d (sin(x), [x, 0, 2*%pi], [ps_file, concat (maxima_tempdir, "/sin.eps")])$
81 "The next example uses the y option ..." $
83 plot2d ([gamma(x), 1/gamma(x)], [x, -4.5, 5], [y, -10, 10], [gnuplot_preamble, "set key bottom"])$
85 "We can also use a `gnuplot_preamble' to produce fancy x-axis labels." $
87 my_preamble: "set xtics ('-2pi' -6.283, \
88 '-3pi/2' -4.712, '-pi' -3.1415, '-pi/2' -1.5708, '0' 0, \
89 'pi/2' 1.5708, 'pi' 3.1415,'3pi/2' 4.712, '2pi' 6.283)"$
90 plot2d([cos(x), sin(x), tan(x), cot(x)], [x, -2*%pi, 2.1*%pi], [y, -2, 2],
91  [axes, x], [gnuplot_preamble, my_preamble]);
93 "... fancy x-axis labels, and produces PostScript output ..." $
95 my_preamble: "set xtics ('-2{/Symbol p}' \
96 -6.283, '-3{/Symbol p}/2' -4.712, '-{/Symbol p}' -3.1415, \
97 '-{/Symbol p}/2' -1.5708, '0' 0,'{/Symbol p}/2' 1.5708, \
98 '{/Symbol p}' 3.1415,'3{/Symbol p}/2' 4.712, '2{/Symbol p}' \
99 6.283)"$
100 plot2d ([cos(x), sin(x), tan(x)], [x, -2*%pi, 2*%pi],
101  [y, -2, 2], [gnuplot_preamble, my_preamble], [ps_file, concat (maxima_tempdir, "/trig.eps")]);
103 /* Examples copied from plot3d documentation */
105 "Displays a plot of one or three expressions as functions of two variables." $
107 plot3d (2^(-u^2 + v^2), [u, -3, 3], [v, -2, 2]);
109 "The same graph can be plotted using xmaxima ..." $
111 plot3d (2^(-u^2 + v^2), [u, -3, 3], [v, -2, 2], [plot_format, xmaxima]);
113 "An example of the third pattern of arguments is" $
115 plot3d ([cos(x)*(3 + y*cos(x/2)), sin(x)*(3 + y*cos(x/2)),
116  y*sin(x/2)], [x, -%pi, %pi], [y, -1, 1], ['grid, 50, 15]);
118 "This example shows a plot of the real part of `z^1/3' ..." $
120 plot3d (r^.33*cos(th/3), [r, 0, 1], [th, 0, 6*%pi], ['grid, 12, 80],
121  ['transform_xy, polar_to_xy], [box, false],[legend,false]);
123 "Other examples are the Klein bottle:" $
125 expr_1: 5*cos(x)*(cos(x/2)*cos(y) + sin(x/2)*sin(2*y)+ 3.0) - 10.0$
126 expr_2: -5*sin(x)*(cos(x/2)*cos(y) + sin(x/2)*sin(2*y)+ 3.0)$
127 expr_3: 5*(-sin(x/2)*cos(y) + cos(x/2)*sin(2*y))$
128 plot3d ([expr_1, expr_2, expr_3], [x, -%pi, %pi], [y, -%pi, %pi], ['grid, 40, 40]);
130 "and a torus:" $
132 expr_1: cos(y)*(10.0+6*cos(x))$
133 expr_2: sin(y)*(10.0+6*cos(x))$
134 expr_3: -6*sin(x)$
135 plot3d ([expr_1, expr_2, expr_3], [x, 0, 2*%pi], [y, 0, 2*%pi], ['grid, 40, 40]);
137 "Sometimes it is necessary to define a function to plot the expression." $
139 kill(f) $
140 M: matrix([1, 2, 3, 4], [1, 2, 3, 2], [1, 2, 3, 4], [1, 2, 3, 3])$
141 f(x, y) := float (M [?round(x), ?round(y)])$
142 plot3d (f, [x, 1, 4], [y, 1, 4], ['grid, 4, 4])$
143 kill(M,f)$
145 "Here is a three-dimensional plot using the gnuplot pm3d terminal." $
147 plot3d (atan (-x^2 + y^3/4), [x, -4, 4], [y, -4, 4], [grid, 50, 50], [gnuplot_pm3d, true])$
149 "And a three-dimensional plot without a mesh and with contours ..." $
151 my_preamble: "set pm3d at s;unset surface;set contour;\
152 set cntrparam levels 20;unset key"$
153 plot3d(atan(-x^2 + y^3/4), [x, -4, 4], [y, -4, 4], [grid, 50, 50],
154  [gnuplot_pm3d, true], [gnuplot_preamble, my_preamble])$
156 "Finally, a plot where the z-axis is represented by color only." $
158 plot3d (cos (-x^2 + y^3/4), [x, -4, 4], [y, -4, 4],
159  [gnuplot_preamble, "set view map; unset surface"], [gnuplot_pm3d, true], [grid, 150, 150])$
161 /* Examples copied from plot_options documentation */
163 "Option: `plot_realpart'" $
165 plot2d (log(x), [x, -5, 5], [plot_realpart, false]);
166 plot2d (log(x), [x, -5, 5], [plot_realpart, true]);
168 /* further *plot-realpart* examples from mailing list 2010-08-13 "plot3d: function not defined everywhere in the plotting range" */
170 plot2d (log(x), [x, -3, 3]);
171 plot3d (log(x), [x, -3, 3], [y, -3, 3]);
173 /* (other options listed under plot_options do not have plot2d or plot3d examples) */
175 /* Examples adapted from demo/demo.dem
176  * (others are duplicates or not different enough)
177  */
179 "two dimensional parametric plot" $
181 plot2d ([parametric, t*sin(t), t*cos(t), [t, 0, 80]], [nticks, 1000]);
183 "three dimensional cartesian plot of a bessel function" $
185 plot3d (bessel_j (0, sqrt(x^2 + y^2)), [x, -12, 12], [y, -12, 12]);
187 "three dimensional polar plot of the same bessel function" $
189 plot3d (bessel_j (0, r), [th, 0, 2*%pi], [r, 0, 12], ['transform_xy, polar_to_xy]);
191 "three dimensional plot of x*exp(-x^2-y^2)" $
193 plot3d (x*exp(- x^2 - y^2), [x, -2, 2], [y, -2, 2]);
195 /* Example adapted from demo/plots.mac
196  * (others are duplicates or not different enough)
197  */
199 "Real part of z^1/6" $
201 plot3d (r^(1/6.0)*cos(th/6), [r, 0, 1], [th, 0, 2*6*%pi], ['grid, 12, 121], ['transform_xy, polar_to_xy]);
203 /* Examples related to bug reports or other specific topics */
205 /* SF bug [ 1699445 ] plot2d in very narrow range
206  * triggers Gnuplot bug on Windows (line outside bounding box)
207  */
209 "[ 1699445 ] plot2d in very narrow range" $
211 plot2d (sin(x), [x, 1.57079628, 1.570796326794897]);
213 "[ 2234113 ] plot2d odd roots of X plots only positive values" $
215 plot2d (x^(1/3), [x, -5, 5]);
217 plot2d (u^(1/5), [u, -5, 5]);
219 "r1.55 src/plot.lisp: ensure that log plots are adequately sampled" $
221 plot2d (x, [x, 1e-5, 100], [logy]);
223 plot2d (x, [x, 1e-5, 100], [logx]);
225 "r1.62 src/plot.lisp: expand the cases recognized by COERCE-FLOAT-FUN" $
227 :lisp (defun $f (x) (+ (cl:sin x) 0.1))
228 :lisp (defmspec $h (x) (+ (cl:sin (cadr x)) 0.3))
229 (g(x) := sin(x) + 0.2,
230  i(x) ::= sin(x) + 0.4,
231  prefix ("j"),
232  "j"(x) := sin(x) + 0.5,
233  matchdeclare (x, floatnump),
234  tellsimpafter (k(x), sin(x) + 0.6));
236 plot2d ([f, g, h, i, sin, k], [u, 0, 1]);
238 /* "j" in list of functions tickles a bug -- quote marks not sanitized
239  * work around it by explicit legend spec
240  */
241  plot2d (["+", "j"], [u, 0, 1], [legend, "+", "j"]);
243 /* Bug reported to mailing list by G. Königsmann, 2023-10-04
244    Legend question for plot() */
245 plot2d(sin(t),[t,0,8],[legend,"Position"]);
247 "(not a plot example, maybe it should be moved elsewhere)" $
249 map (lambda ([f], quad_qags (f, u, 0, 1)), [f, g, h, i, "+", "j", sin, k]);
250 :lisp (unintern '$f)
251 :lisp (unintern '$h)
253 "r1.82 src/plot.lisp: floatify non-float numbers" $
255 plot2d ([discrete, [1, 2, 3, 4], [1, 1/2, 1/3, 1/4]]);
257 plot2d ([discrete, [1, 2, 3, 4], [1/%e, 1/%pi, 1/%phi, 1/%gamma]]);
259 plot2d ([discrete, [1, 2, 3, 4], [1b0, 0.5b0, 0.33b0, 0.25b0]]);
261 "SF bug [ 2672976 ] wxMaxima 0.8.1: set logscale x gives error";
263 plot2d (sin(x), [x, 0, 1], [logx, true]);
265 plot2d (sin(u), [u, 0, 1], [logx, true]);
267 plot2d (sin(u), [u, 0, 1], [logx, true], [y, 0, 1], [logy, true]);
269 "verify that nested numerical integral is handled correctly";
271 plot2d (w^2 * 'quad_qags (1/(s - w), s, 1, 5) [1], [w, -5, -1], [adapt_depth, 1]);
273 "another nested numerical integral";
275 (kill (W, R, A),
276  W(t) := 95*sqrt(t)*sin(t/6)^2,
277  R(t) := 275*sin(t/3)^2,
278  A(t) := 1200 + quad_qags (W(x) - R(x), x, 0, t) [1],
279  plot2d (A(t), [t, 0, 18], [adapt_depth, 1]));
281 "plotting realpart and imagpart";
283 plot3d (realpart (asin (x + y*%i)), [x, -2, 2], [y, -2, 2], [grid, 20, 20]);
285 plot3d (imagpart (asin (x + y*%i)), [x, -2, 2], [y, -2, 2], [grid, 20, 20]);
287 plot3d (realpart (exp (x + y*%i)), [x, -2, 2], [y, -2, 2], [grid, 20, 20]);
289 plot3d (imagpart (exp (x + y*%i)), [x, -2, 2], [y, -2, 2], [grid, 20, 20]);
291 "messages about non-numeric values and clipped values";
293 plot2d (sqrt (x), [x, 0, 1]); /* no message */
295 plot2d (sqrt (x), [x, -1, 1]); /* some non-numeric */
297 plot2d (sqrt (x), [x, 0, 1], [y, 0, 1/2]); /* some clipped */
299 plot2d (sqrt (x), [x, -1, 1], [y, 0, 1/2]); /* some non-numeric, some clipped */
301 errcatch(plot2d (sqrt (x), [x, -2, -1])); /* all non-numeric */
303 errcatch(plot2d (sqrt (x), [x, 0, 1], [y, -2, -1])); /* all clipped */
305 errcatch(plot2d (sqrt (x), [x, -1, 1], [y, -2, -1])); /* all non-numeric or clipped */
307 "coping with overflow in intermediate results";
309 /* from the mailing list "plot numerical question" 2009-03-11 */
311 (r4(s) := block ([s : bfloat(s)], s!/(((s/4)!)^4 * 4^s)),
312  plot2d (r4, [s, 200, 300], [adapt_depth, 1], [nticks, 5]));
314 plot2d ('r4(s), [s, 200, 300], [adapt_depth, 1], [nticks, 5]);
316 /* from mailing list 2009-02-18
317  * "Re: [Maxima] I want to tell maxima (-1)^0.33333333333333=-1, what should i do?"
318  */
320 foo29(x):=(sqrt(-16*x^4-16*x^3+20*x^2+12*x+23)/(6*sqrt(3))+(16*x^3-12*x^2-6*x-25)/54)^(1/3)$
322 plot2d (foo29 (u), [u, -1, 0]);
324 plot2d (foo29, [u, -1, 0]);
326 compile (foo29);
328 plot2d (foo29, [u, -1, 0]);
330 /* bug report "Wrong result given by coerce-float-fun" ID: 2880115 */
332 (kill(f),
333  f(k) := integrate (exp(%i*k*x)*sin(x)/x, x, minf, inf),
334  plot2d (f, [x, -3, 3], [adapt_depth, 1], [nticks, 5]));
336 (translate(f),
337  plot2d (f, [x, -3, 3], [adapt_depth, 1], [nticks, 5]));
339 /* bug report "xlabel and ylabel don't change plot3d axis labels" - ID: 3020589 */
341 (Bxt(x, r) := x^2 + r^2,
342  [R, L] : [1, 1],
343  plot3d(Bxt(x,r),[x,0,L],[r,0,R],[xlabel,"x [m]"],[ylabel,"r [m]"],[zlabel,"Bx [T]"],[legend,"Axial field"]));
345 /* same but now using default axis labels */
346 plot3d(Bxt(x,r),[x,0,L],[r,0,R],[legend,"Axial field"]);
348 /* from the mailing list 2011-05-26:
349  * "Wrong usage" error message when trying to plot with plot3d or contour_plot
350  * should expect this example to provoke an advisory message and make a plot
351  */
353 block (local (fn, Fn, pw, fbn, fbnxcy, loww),
354  fn(x):=1/sqrt(2*%pi)*exp((-x^2)/2),
355  Fn(x):=integrate(fn(t),t,-inf,x),
356  pw(r1,r2,s1,s2):=1-Fn((r2-r1)/sqrt(s1^2+s2^2)),
357  fbn(x,y,r):=1/(2*%pi*sqrt(1-r^2))*exp((-(x^2-2*r*x*y+y^2))/(2*(1-r^2))),
358  fbnxcy(x,y,r) := fbn(x,y,r) / fn(y),
359  loww(rs2, ro, r, s1, s2) := quad_qagi(fbnxcy(rs1, rs2, r)*pw(rs1, ro, s1, s2), rs1, minf, inf),
360  plot3d( 'loww(rs2, 0, r, 1, 1)[1] , [rs2, -1.5, 2.5], [r, 0.1, 0.9] ));
362 /* simpler version of the preceding one */
364 block (local (foo),
365  foo(x, y) := if numberp(x) and numberp(y) then [x^2 - y^2] else funmake (foo, [x, y]),
366  plot3d (foo (a, b)[1], [a, -1, 1], [b, -1, 1]));
368 /* helix -- see mailing list 2012-01-22 "3D curve parametric plot" */
370 plot3d([sin(t), cos(t), t], [t,-5,5], [y,-5,5], [grid,100,2], [gnuplot_pm3d,false])$
372 /* plot2d/plot3d with subscripted variable */
374 (kill (x, a), plot2d (a[x]^3, [a[x], -1, 1]));
376 plot3d (a[x]^2 - x[a]^3, [a[x], -1, 1], [x[a], -1, 1]);
378 /* from wxMaxima forum 2013-02-16: "plot2d: expression evaluates to non-numeric value everywhere in plotting range"
379  * https://sourceforge.net/p/wxmaxima/discussion/435775/thread/d89ea980/
380  */
382 F(omega) := abs((%i * omega + 2E4)^2);
383 plot2d(F(omega), [omega, 0.01, 2E5]);
385 /* SF bug 3776: "plot3d(\"*\",...) gives internal error" */
387 (kill(mul, x, y), mul(a,b):=a*b)$
388 /* next two should be the same */
389 plot3d(mul,[x,0,1],[y,0,1]);
390 plot3d("*",[x,0,1],[y,0,1]);
392 /* SF bug #3807: "plot2d heuristic to detect unbound variables excludes some valid cases" */
394 plot2d (quad_qags (sin(u), u, 0, x)[1], [x, 0, 1]);
395 plot2d (find_root (sin(u) - x, u, 0, %pi/2), [x, 0, 1]);
396 plot2d(sum(w, w, 1, abs(x)), [x, 1, 10]);
398 /* Bug #3881: plot2d not ignoring errors within functions */
400 plot2d (gamma (x), [x, -1, 1]);
401 plot2d (lambda ([x], gamma (x)), [x, -1, 1]);
403 /* SF bug #1055: "parameter in parametric plot must be named t" */
405 kill (foo); 
406 plot2d ([parametric, sin(foo), cos(foo), [foo, 0, %pi]]);
408 /* SF bug #3958: "plot2d with multiple discrete plots fails" */
409 /* SF bug #3959: "plot2d + Gnuplot 4 with `plot title noenhanced`" */
411 plot2d ([[discrete, [1, 2, 3, 4], [1, 1/2, 1/3, 1/4]], [discrete, [1, 2, 3, 4], [1/%e, 1/%pi, 1/%phi, 1/%gamma]]], [legend, "foo_1", "foo_2"], [gnuplot_strings, false]);
413 plot2d ([[discrete, [1, 2, 3, 4], [1, 1/2, 1/3, 1/4]], [discrete, [1, 2, 3, 4], [1/%e, 1/%pi, 1/%phi, 1/%gamma]]], [legend, "baz_1", "baz_2"], [gnuplot_strings, true]);
415 /* gnuplot_svg_background */
417 /* SF bug #4294: plot_format = gnuplot_pipes (plot2d default) leaks settings from one plot to the next.
418  * Set plot_format = gnuplot in order to ensure each plot is entirely separate,
419  * otherwise the outputs for nognuplot_svg_background are incorrect.
420  */
422 save_plot_format: get_plot_option (plot_format);
423 set_plot_option ([plot_format, gnuplot]);
425 get_plot_option (gnuplot_svg_background);
427 plot2d ([discrete, [1, 2, 3], [4, 5, 6]], [svg_file, "./tmp-svg-bg-white.svg"]);
428 plot2d ([discrete, [1, 2, 3], [4, 5, 6]], [svg_file, "./tmp-svg-bg-no-bg-1.svg"], nognuplot_svg_background);
429 plot2d ([discrete, [1, 2, 3], [4, 5, 6]], [svg_file, "./tmp-svg-bg-gray.svg"], [gnuplot_svg_background, "gray"]);
430 plot2d ([discrete, [1, 2, 3], [4, 5, 6]], [svg_file, "./tmp-svg-bg-ffeedd-1.svg"], [gnuplot_svg_background, "#FFEEDD"]);
432 get_plot_option (gnuplot_svg_background);
433 set_plot_option ([gnuplot_svg_background, "yellow"]);
435 plot2d ([discrete, [1, 2, 3], [4, 5, 6]], [svg_file, "./tmp-svg-bg-yellow.svg"]);
436 plot2d ([discrete, [1, 2, 3], [4, 5, 6]], [svg_file, "./tmp-svg-bg-no-bg-2.svg"], nognuplot_svg_background);
437 plot2d ([discrete, [1, 2, 3], [4, 5, 6]], [svg_file, "./tmp-svg-bg-ffeedd-2.svg"], [gnuplot_svg_background, "#FFEEDD"]);
439 get_plot_option (gnuplot_svg_background);
440 set_plot_option ([gnuplot_svg_background, false]);
442 plot2d ([discrete, [1, 2, 3], [4, 5, 6]], [svg_file, "./tmp-svg-bg-no-bg-3.svg"]);
444 get_plot_option (gnuplot_svg_background);
446 set_plot_option (save_plot_format);
448 /* Allow an option to set the pointinterval gnuplot setting with the 6th argument of the style list. */
450 plot2d(x^3,[x,-2,2],[style, [linespoints, 1, 1.2, 1, 1, 30]]);
452 "FINIS" $