tests/rtest_plot.mac

   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  */
  12
  13 /* Examples copied from plot2d documentation */
  14
  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]]);
  28
  29 "Plots of common functions." $
  30
  31 plot2d (sin(x), [x, -5, 5])$
  32 plot2d (sec(x), [x, -2, 2], [y, -20, 20], [nticks, 200])$
  33
  34 "Plotting functions by name." $
  35
  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) $
  42
  43 "We can plot a circle using a parametric plot ..." $
  44
  45 plot2d ([parametric, cos(t), sin(t), [t,-%pi,%pi]], [nticks,80], [x, -4/3, 4/3])$
  46
  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" $
  49
  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])$
  51
  52 "Combination of an ordinary plot ..." $
  53
  54 plot2d ([x^3+2, [parametric, cos(t), sin(t), [t, -5, 5]]], [x, -3, 3],  [nticks, 80])$
  55
  56 "Example of a logarithmic plot:" $
  57
  58 plot2d (exp(3*s), [s, -2, 2], [logy])$
  59
  60 "To show some examples of discrete plots, ..." $
  61
  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])$
  66
  67 "We will now show the plot with only points, ..." $
  68
  69 plot2d([discrete, xy], [style, points])$
  70
  71 "The plot of the data points can be shown together with ..." $
  72
  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)"])$
  76
  77 "To save a plot ..." $
  78
  79 plot2d (sin(x), [x, 0, 2*%pi], [ps_file, concat (maxima_tempdir, "/sin.eps")])$
  80
  81 "The next example uses the y option ..." $
  82
  83 plot2d ([gamma(x), 1/gamma(x)], [x, -4.5, 5], [y, -10, 10], [gnuplot_preamble, "set key bottom"])$
  84
  85 "We can also use a `gnuplot_preamble' to produce fancy x-axis labels." $
  86
  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]);
  92
  93 "... fancy x-axis labels, and produces PostScript output ..." $
  94
  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")]);
 102
 103 /* Examples copied from plot3d documentation */
 104
 105 "Displays a plot of one or three expressions as functions of two variables." $
 106
 107 plot3d (2^(-u^2 + v^2), [u, -3, 3], [v, -2, 2]);
 108
 109 "The same graph can be plotted using xmaxima ..." $
 110
 111 plot3d (2^(-u^2 + v^2), [u, -3, 3], [v, -2, 2], [plot_format, xmaxima]);
 112
 113 "An example of the third pattern of arguments is" $
 114
 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]);
 117
 118 "This example shows a plot of the real part of `z^1/3' ..." $
 119
 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]);
 122
 123 "Other examples are the Klein bottle:" $
 124
 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]);
 129
 130 "and a torus:" $
 131
 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]);
 136
 137 "Sometimes it is necessary to define a function to plot the expression." $
 138
 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)$
 144
 145 "Here is a three-dimensional plot using the gnuplot pm3d terminal." $
 146
 147 plot3d (atan (-x^2 + y^3/4), [x, -4, 4], [y, -4, 4], [grid, 50, 50], [gnuplot_pm3d, true])$
 148
 149 "And a three-dimensional plot without a mesh and with contours ..." $
 150
 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])$
 155
 156 "Finally, a plot where the z-axis is represented by color only." $
 157
 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])$
 160
 161 /* Examples copied from plot_options documentation */
 162
 163 "Option: `plot_realpart'" $
 164
 165 plot2d (log(x), [x, -5, 5], [plot_realpart, false]);
 166 plot2d (log(x), [x, -5, 5], [plot_realpart, true]);
 167
 168 /* further *plot-realpart* examples from mailing list 2010-08-13 "plot3d: function not defined everywhere in the plotting range" */
 169
 170 plot2d (log(x), [x, -3, 3]);
 171 plot3d (log(x), [x, -3, 3], [y, -3, 3]);
 172
 173 /* (other options listed under plot_options do not have plot2d or plot3d examples) */
 174
 175 /* Examples adapted from demo/demo.dem
 176  * (others are duplicates or not different enough)
 177  */
 178
 179 "two dimensional parametric plot" $
 180
 181 plot2d ([parametric, t*sin(t), t*cos(t), [t, 0, 80]], [nticks, 1000]);
 182
 183 "three dimensional cartesian plot of a bessel function" $
 184
 185 plot3d (bessel_j (0, sqrt(x^2 + y^2)), [x, -12, 12], [y, -12, 12]);
 186
 187 "three dimensional polar plot of the same bessel function" $
 188
 189 plot3d (bessel_j (0, r), [th, 0, 2*%pi], [r, 0, 12], ['transform_xy, polar_to_xy]);
 190
 191 "three dimensional plot of x*exp(-x^2-y^2)" $
 192
 193 plot3d (x*exp(- x^2 - y^2), [x, -2, 2], [y, -2, 2]);
 194
 195 /* Example adapted from demo/plots.mac
 196  * (others are duplicates or not different enough)
 197  */
 198
 199 "Real part of z^1/6" $
 200
 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]);
 202
 203 /* Examples related to bug reports or other specific topics */
 204
 205 /* SF bug [ 1699445 ] plot2d in very narrow range
 206  * triggers Gnuplot bug on Windows (line outside bounding box)
 207  */
 208
 209 "[ 1699445 ] plot2d in very narrow range" $
 210
 211 plot2d (sin(x), [x, 1.57079628, 1.570796326794897]);
 212
 213 "[ 2234113 ] plot2d odd roots of X plots only positive values" $
 214
 215 plot2d (x^(1/3), [x, -5, 5]);
 216
 217 plot2d (u^(1/5), [u, -5, 5]);
 218
 219 "r1.55 src/plot.lisp: ensure that log plots are adequately sampled" $
 220
 221 plot2d (x, [x, 1e-5, 100], [logy]);
 222
 223 plot2d (x, [x, 1e-5, 100], [logx]);
 224
 225 "r1.62 src/plot.lisp: expand the cases recognized by COERCE-FLOAT-FUN" $
 226
 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));
 235
 236 /* BUG: this fails, but should not */
 237 errcatch(
 238   plot2d ([f, g, h, i, sin, k], [u, 0, 1]));
 239
 240 /* "j" in list of functions tickles a bug -- quote marks not sanitized
 241  * work around it by explicit legend spec
 242  */
 243 plot2d (["+", "j"], [u, 0, 1], [legend, "+", "j"]);
 244
 245 "(not a plot example, maybe it should be moved elsewhere)" $
 246
 247 map (lambda ([f], quad_qags (f, u, 0, 1)), [f, g, h, i, "+", "j", sin, k]);
 248 :lisp (unintern '$f)
 249 :lisp (unintern '$h)
 250
 251 "r1.82 src/plot.lisp: floatify non-float numbers" $
 252
 253 plot2d ([discrete, [1, 2, 3, 4], [1, 1/2, 1/3, 1/4]]);
 254
 255 plot2d ([discrete, [1, 2, 3, 4], [1/%e, 1/%pi, 1/%phi, 1/%gamma]]);
 256
 257 plot2d ([discrete, [1, 2, 3, 4], [1b0, 0.5b0, 0.33b0, 0.25b0]]);
 258
 259 "SF bug [ 2672976 ] wxMaxima 0.8.1: set logscale x gives error";
 260
 261 plot2d (sin(x), [x, 0, 1], [logx, true]);
 262
 263 plot2d (sin(u), [u, 0, 1], [logx, true]);
 264
 265 errcatch(plot2d (sin(u), [u, 0, 1], [logx, true], [y, 0, 1], [logy, true]));
 266
 267 "verify that nested numerical integral is handled correctly";
 268
 269 /* BUG: this fails, but should not */
 270 errcatch(
 271   plot2d (w^2 * 'quad_qags (1/(s - w), s, 1, 5) [1], [w, -5, -1], [adapt_depth, 1]));
 272
 273 "another nested numerical integral";
 274
 275 /* BUG: this fails, but should not */
 276 errcatch(kill (W, R, A),
 277  W(t) := 95*sqrt(t)*sin(t/6)^2,
 278  R(t) := 275*sin(t/3)^2,
 279  A(t) := 1200 + quad_qags (W(x) - R(x), x, 0, t) [1],
 280  plot2d (A(t), [t, 0, 18], [adapt_depth, 1]));
 281
 282 "plotting realpart and imagpart";
 283
 284 plot3d (realpart (asin (x + y*%i)), [x, -2, 2], [y, -2, 2], [grid, 20, 20]);
 285
 286 plot3d (imagpart (asin (x + y*%i)), [x, -2, 2], [y, -2, 2], [grid, 20, 20]);
 287
 288 plot3d (realpart (exp (x + y*%i)), [x, -2, 2], [y, -2, 2], [grid, 20, 20]);
 289
 290 plot3d (imagpart (exp (x + y*%i)), [x, -2, 2], [y, -2, 2], [grid, 20, 20]);
 291
 292 "messages about non-numeric values and clipped values";
 293
 294 plot2d (sqrt (x), [x, 0, 1]); /* no message */
 295
 296 plot2d (sqrt (x), [x, -1, 1]); /* some non-numeric */
 297
 298 plot2d (sqrt (x), [x, 0, 1], [y, 0, 1/2]); /* some clipped */
 299
 300 plot2d (sqrt (x), [x, -1, 1], [y, 0, 1/2]); /* some non-numeric, some clipped */
 301
 302 errcatch(plot2d (sqrt (x), [x, -2, -1])); /* all non-numeric */
 303
 304 errcatch(plot2d (sqrt (x), [x, 0, 1], [y, -2, -1])); /* all clipped */
 305
 306 errcatch(plot2d (sqrt (x), [x, -1, 1], [y, -2, -1])); /* all non-numeric or clipped */
 307
 308 "coping with overflow in intermediate results";
 309
 310 /* from the mailing list "plot numerical question" 2009-03-11 */
 311
 312 (r4(s) := block ([s : bfloat(s)], s!/(((s/4)!)^4 * 4^s)),
 313  plot2d (r4, [s, 200, 300], [adapt_depth, 1], [nticks, 5]));
 314
 315 plot2d ('r4(s), [s, 200, 300], [adapt_depth, 1], [nticks, 5]);
 316
 317 /* from mailing list 2009-02-18
 318  * "Re: [Maxima] I want to tell maxima (-1)^0.33333333333333=-1, what should i do?"
 319  */
 320
 321 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
 323 plot2d (foo29 (u), [u, -1, 0]);
 324
 325 plot2d (foo29, [u, -1, 0]);
 326
 327 compile (foo29);
 328
 329 plot2d (foo29, [u, -1, 0]);
 330
 331 /* bug report "Wrong result given by coerce-float-fun" ID: 2880115 */
 332
 333 (kill(f),
 334  f(k) := integrate (exp(%i*k*x)*sin(x)/x, x, minf, inf),
 335  plot2d (f, [x, -3, 3], [adapt_depth, 1], [nticks, 5]));
 336
 337 (translate(f),
 338  plot2d (f, [x, -3, 3], [adapt_depth, 1], [nticks, 5]));
 339
 340 /* bug report "xlabel and ylabel don't change plot3d axis labels" - ID: 3020589 */
 341
 342 (Bxt(x, r) := x^2 + r^2,
 343  [R, L] : [1, 1],
 344  plot3d(Bxt(x,r),[x,0,L],[r,0,R],[xlabel,"x [m]"],[ylabel,"r [m]"],[zlabel,"Bx [T]"],[legend,"Axial field"]));
 345
 346 /* same but now using default axis labels */
 347 plot3d(Bxt(x,r),[x,0,L],[r,0,R],[legend,"Axial field"]);
 348
 349 /* from the mailing list 2011-05-26:
 350  * "Wrong usage" error message when trying to plot with plot3d or contour_plot
 351  * should expect this example to provoke an advisory message and make a plot
 352  */
 353
 354 block (local (fn, Fn, pw, fbn, fbnxcy, loww),
 355  fn(x):=1/sqrt(2*%pi)*exp((-x^2)/2),
 356  Fn(x):=integrate(fn(t),t,-inf,x),
 357  pw(r1,r2,s1,s2):=1-Fn((r2-r1)/sqrt(s1^2+s2^2)),
 358  fbn(x,y,r):=1/(2*%pi*sqrt(1-r^2))*exp((-(x^2-2*r*x*y+y^2))/(2*(1-r^2))),
 359  fbnxcy(x,y,r) := fbn(x,y,r) / fn(y),
 360  loww(rs2, ro, r, s1, s2) := quad_qagi(fbnxcy(rs1, rs2, r)*pw(rs1, ro, s1, s2), rs1, minf, inf),
 361  plot3d( 'loww(rs2, 0, r, 1, 1)[1] , [rs2, -1.5, 2.5], [r, 0.1, 0.9] ));
 362
 363 /* simpler version of the preceding one */
 364
 365 block (local (foo),
 366  foo(x, y) := if numberp(x) and numberp(y) then [x^2 - y^2] else funmake (foo, [x, y]),
 367  plot3d (foo (a, b)[1], [a, -1, 1], [b, -1, 1]));
 368
 369 /* helix -- see mailing list 2012-01-22 "3D curve parametric plot" */
 370
 371 plot3d([sin(t), cos(t), t], [t,-5,5], [y,-5,5], [grid,100,2], [gnuplot_pm3d,false])$
 372
 373 /* plot2d/plot3d with subscripted variable */
 374
 375 /* BUG? this fails, but should not */
 376 errcatch(kill (x, a), plot2d (a[x]^3, [a[x], -1, 1]));
 377
 378 plot3d (a[x]^2 - x[a]^3, [a[x], -1, 1], [x[a], -1, 1]);
 379
 380 /* from wxMaxima forum 2013-02-16: "plot2d: expression evaluates to non-numeric value everywhere in plotting range"
 381  * https://sourceforge.net/p/wxmaxima/discussion/435775/thread/d89ea980/
 382  */
 383
 384 F(omega) := abs((%i * omega + 2E4)^2);
 385 plot2d(F(omega), [omega, 0.01, 2E5]);
 386
 387 "FINIS" $
 388