1)abbrev domain PLOT3D Plot3D 2++ Author: Clifton J. Williamson based on code by Michael Monagan 3++ Date Created: Jan 1989 4++ Basic Operations: pointPlot, plot, zoom, refine, tRange, tValues, 5++ minPoints3D, setMinPoints3D, maxPoints3D, setMaxPoints3D, 6++ screenResolution3D, setScreenResolution3D, adaptive3D?, setAdaptive3D, 7++ numFunEvals3D, debug3D 8++ Related Constructors: 9++ Also See: 10++ AMS Classifications: 11++ Keywords: plot, parametric 12++ References: 13++ Description: Plot3D supports parametric plots defined over a real 14++ number system. A real number system is a model for the real 15++ numbers and as such may be an approximation. For example, 16++ floating point numbers and infinite continued fractions are 17++ real number systems. The facilities at this point are limited 18++ to 3-dimensional parametric plots. 19Plot3D() : Exports == Implementation where 20 B ==> Boolean 21 F ==> DoubleFloat 22 I ==> Integer 23 L ==> List 24 N ==> NonNegativeInteger 25 OUT ==> OutputForm 26 P ==> Point F 27 S ==> String 28 R ==> Segment F 29 O ==> OutputPackage 30 C ==> Record(source : F -> P, ranges : L R, knots : L F, points : L P) 31 32 Exports ==> PlottableSpaceCurveCategory with 33 34 pointPlot : (F -> P, R) -> % 35 ++ pointPlot(f, g, h, a..b) plots x = f(t), y = g(t), z = h(t) as 36 ++ t ranges over [a, b]. 37 pointPlot : (F -> P, R, R, R, R) -> % 38 ++ pointPlot(f, x, y, z, w) \undocumented 39 plot : (F -> F, F -> F, F -> F, F -> F, R) -> % 40 ++ plot(f, g, h, a..b) plots x = f(t), y = g(t), z = h(t) as 41 ++ t ranges over [a, b]. 42 plot : (F -> F, F -> F, F -> F, F -> F, R, R, R, R) -> % 43 ++ plot(f1, f2, f3, f4, x, y, z, w) \undocumented 44 45 plot : (%, R) -> % -- change the range 46 ++ plot(x, r) \undocumented 47 zoom : (%, R, R, R) -> % 48 ++ zoom(x, r, s, t) \undocumented 49 refine : (%, R) -> % 50 ++ refine(x, r) \undocumented 51 refine : % -> % 52 ++ refine(x) \undocumented 53 54 tRange : % -> R 55 ++ tRange(p) returns the range of the parameter in a parametric plot p. 56 tValues : % -> L L F 57 ++ tValues(p) returns a list of lists of the values of the parameter for 58 ++ which a point is computed, one list for each curve in the plot p. 59 60 minPoints3D : () -> I 61 ++ minPoints3D() returns the minimum number of points in a plot. 62 setMinPoints3D : I -> I 63 ++ setMinPoints3D(i) sets the minimum number of points in a plot to i. 64 maxPoints3D : () -> I 65 ++ maxPoints3D() returns the maximum number of points in a plot. 66 setMaxPoints3D : I -> I 67 ++ setMaxPoints3D(i) sets the maximum number of points in a plot to i. 68 screenResolution3D : () -> I 69 ++ screenResolution3D() returns the screen resolution for a 3d graph. 70 setScreenResolution3D : I -> I 71 ++ setScreenResolution3D(i) sets the screen resolution for a 3d graph to i. 72 adaptive3D? : () -> B 73 ++ adaptive3D?() determines whether plotting be done adaptively. 74 setAdaptive3D : B -> B 75 ++ setAdaptive3D(true) turns adaptive plotting on; 76 ++ setAdaptive3D(false) turns adaptive plotting off. 77 numFunEvals3D : () -> I 78 ++ numFunEvals3D() returns the number of points computed. 79 debug3D : B -> B 80 ++ debug3D(true) turns debug mode on; 81 ++ debug3D(false) turns debug mode off. 82 83 Implementation ==> add 84 import from PointPackage(F) 85 86--% local functions 87 88 fourth : L R -> R 89 checkRange : R -> R 90 -- checks that left-hand endpoint is less than right-hand endpoint 91 intersect : (R, R) -> R 92 -- intersection of two intervals 93 union : (R, R) -> R 94 -- union of two intervals 95 join : (L C, I) -> R 96 parametricRange : % -> R 97-- setColor : (P, F) -> F 98 select : (L P, P -> F, (F, F) -> F) -> F 99-- normalizeColor : (P, F, F) -> F 100 rangeRefine : (C, R) -> C 101 adaptivePlot : (C, R, R, R, R, I, I) -> C 102 basicPlot : (F -> P, R) -> C 103 basicRefine : (C, R) -> C 104 point : (F, F, F, F) -> P 105 106--% representation 107 108 Rep := Record( display : L R, _ 109 bounds : L R, _ 110 screenres : I, _ 111 axisLabels : L S, _ 112 functions : L C ) 113 114--% global constants 115 116 ADAPTIVE : B := true 117 MINPOINTS : I := 49 118 MAXPOINTS : I := 1000 119 NUMFUNEVALS : I := 0 120 SCREENRES : I := 500 121 ANGLEBOUND : F := cos inv (4::F) 122 DEBUG : B := false 123 124 point(xx, yy, zz, col) == point(l : L F := [xx, yy, zz, col]) 125 126 fourth list == first rest rest rest list 127 128 checkRange r == (low(r) > high(r) => error "ranges cannot be negative"; r) 129 intersect(s, t) == checkRange(max(low(s), low(t))..min(high(s), high(t))) 130 union(s : R, t : R) == min(low(s), low(t))..max(high(s), high(t)) 131 join(l, i) == 132 rr := first l 133 u : R := 134 i = 0 => first(rr.ranges) 135 i = 1 => second(rr.ranges) 136 i = 2 => third(rr.ranges) 137 fourth(rr.ranges) 138 for r in rest l repeat 139 i = 0 => u := union(u, first(r.ranges)) 140 i = 1 => u := union(u, second(r.ranges)) 141 i = 2 => u := union(u, third(r.ranges)) 142 u := union(u, fourth(r.ranges)) 143 u 144 parametricRange r == first(r.bounds) 145 146 minPoints3D() == MINPOINTS 147 setMinPoints3D n == 148 if n < 3 then error "three points minimum required" 149 if MAXPOINTS < n then MAXPOINTS := n 150 MINPOINTS := n 151 maxPoints3D() == MAXPOINTS 152 setMaxPoints3D n == 153 if n < 3 then error "three points minimum required" 154 if MINPOINTS > n then MINPOINTS := n 155 MAXPOINTS := n 156 screenResolution3D() == SCREENRES 157 setScreenResolution3D n == 158 if n < 2 then error "buy a new terminal" 159 SCREENRES := n 160 adaptive3D?() == ADAPTIVE 161 setAdaptive3D b == ADAPTIVE := b 162 163 numFunEvals3D() == NUMFUNEVALS 164 debug3D b == DEBUG := b 165 166-- setColor(p, c) == p.colNum := c 167 168 xRange plot == second plot.bounds 169 yRange plot == third plot.bounds 170 zRange plot == fourth plot.bounds 171 tRange plot == first plot.bounds 172 173 tValues plot == 174 outList : L L F := [] 175 for curve in plot.functions repeat 176 outList := concat(curve.knots, outList) 177 outList 178 179 select(l, f, g) == 180 m := f first l 181 -- if (EQL(m, _$NaNvalue$Lisp)$Lisp) then m := 0 182-- for p in rest l repeat m := g(m, fp) 183 for p in rest l repeat 184 fp : F := f p 185 -- if (EQL(fp, _$NaNvalue$Lisp)$Lisp) then fp := 0 186 m := g(m, fp) 187 m 188 189-- normalizeColor(p, lo, diff) == 190-- p.colNum := (p.colNum - lo)/diff 191 192 rangeRefine(curve, nRange) == 193 checkRange nRange; l := low(nRange); h := high(nRange) 194 t := curve.knots; p := curve.points; f := curve.source 195 while not(empty?(t)) and first t < l repeat 196 (t := rest t; p := rest p) 197 c : L F := []; q : L P := [] 198 while not(empty?(t)) and first t <= h repeat 199 c := concat(first t, c); q := concat(first p, q) 200 t := rest t; p := rest p 201 if empty?(c) then return basicPlot(f, nRange) 202 if first c < h then 203 c := concat(h, c); q := concat(f h, q) 204 NUMFUNEVALS := NUMFUNEVALS + 1 205 t := c := reverse! c; p := q := reverse! q 206 s := (h-l)/(MINPOINTS::F-1) 207 if (first t) ~= l then 208 t := c := concat(l, c); p := q := concat(f l, p) 209 NUMFUNEVALS := NUMFUNEVALS + 1 210 while not(empty?(rest(t))) repeat 211 n := wholePart((second(t) - first(t))/s) 212 d := (second(t) - first(t))/((n+1)::F) 213 for i in 1..n repeat 214 t.rest := concat(first(t) + d, rest t); t1 := second t 215 p.rest := concat(f t1, rest p) 216 NUMFUNEVALS := NUMFUNEVALS + 1 217 t := rest t; p := rest p 218 t := rest t 219 p := rest p 220 xRange := select(q, xCoord, min) .. select(q, xCoord, max) 221 yRange := select(q, yCoord, min) .. select(q, yCoord, max) 222 zRange := select(q, zCoord, min) .. select(q, zCoord, max) 223-- colorLo := select(q, color, min); colorHi := select(q, color, max) 224-- (diff := colorHi - colorLo) = 0 => 225-- error "all points are the same color" 226-- map(normalizeColor(#1, colorLo, diff), q)$ListPackage1(P) 227 [f, [nRange, xRange, yRange, zRange], c, q] 228 229 230 adaptivePlot(curve, tRg, xRg, yRg, zRg, pixelfraction, resolution) == 231 xDiff := high(xRg) - low(xRg) 232 yDiff := high(yRg) - low(yRg) 233 zDiff := high(zRg) - low(zRg) 234-- xDiff = 0 or yDiff = 0 or zDiff = 0 => curve--!! delete this? 235 if xDiff = 0::F then xDiff := 1::F 236 if yDiff = 0::F then yDiff := 1::F 237 if zDiff = 0::F then zDiff := 1::F 238 l := low(tRg); h := high(tRg) 239 (tDiff := h-l) = 0 => curve 240 t := curve.knots 241 #t < 3 => curve 242 p := curve.points; f := curve.source 243 minLength : F := 4::F/resolution::F 244 maxLength := 1/4::F 245 tLimit := tDiff/(pixelfraction*resolution)::F 246 while not(empty?(t)) and first t < l repeat (t := rest t; p := rest p) 247 #t < 3 => curve 248 headert := t; headerp := p 249 st := headert; sp := headerp 250 n : I := 0 251 252 while not(empty?(rest(rest(st)))) repeat 253 t0 := first(st); t1 := second(st); t2 := third(st) 254 if t2 > h then break 255 t2 - t0 < tLimit => 256 st := rest st 257 sp := rest sp 258 x0 := xCoord first(sp); y0 := yCoord first(sp); z0 := zCoord first(sp) 259 x1 := xCoord second(sp); y1 := yCoord second(sp); z1 := zCoord second(sp) 260 x2 := xCoord third(sp); y2 := yCoord third(sp); z2 := zCoord third(sp) 261 a1 := (x1-x0)/xDiff; b1 := (y1-y0)/yDiff; c1 := (z1-z0)/zDiff 262 a2 := (x2-x1)/xDiff; b2 := (y2-y1)/yDiff; c2 := (z2-z1)/zDiff 263 s1 := sqrt(a1^2+b1^2+c1^2); s2 := sqrt(a2^2+b2^2+c2^2) 264 dp := a1*a2+b1*b2+c1*c2 265 s1 < maxLength and s2 < maxLength and _ 266 (s1 = 0 or s2 = 0 or 267 s1 < minLength and s2 < minLength or _ 268 dp/s1/s2 > ANGLEBOUND) => 269 st := rest st 270 sp := rest sp 271 if n = MAXPOINTS then break else n := n + 1 272 --if DEBUG then 273 --r : L F := [minLength, maxLength, s1, s2, dp/s1/s2, ANGLEBOUND] 274 --output(r::E)$O 275 t := st 276 p := sp 277 tj := (t0+t1)/2::F 278 qsetrest!(t, cons(tj, rest t)) 279 qsetrest!(p, cons(f tj, rest p)) 280 t := rest t; p := rest p 281 t := rest t; p := rest p 282 283 tj := (t1+t2)/2::F 284 qsetrest!(t, cons(tj, rest t)) 285 qsetrest!(p, cons(f tj, rest p)) 286 if n > 0 then 287 NUMFUNEVALS := NUMFUNEVALS + n 288 t := curve.knots; p := curve.points 289 xRg := select(p, xCoord, min) .. select(p, xCoord, max) 290 yRg := select(p, yCoord, min) .. select(p, yCoord, max) 291 zRg := select(p, zCoord, min) .. select(p, zCoord, max) 292 [curve.source, [tRg, xRg, yRg, zRg], t, p] 293 else curve 294 295 basicPlot(f, tRange) == 296 checkRange tRange; l := low(tRange); h := high(tRange) 297 t : L F := list l; p : L P := list f l 298 s := (h-l)/(MINPOINTS-1)::F 299 for i in 2..MINPOINTS-1 repeat 300 l := l+s; t := concat(l, t) 301 p := concat(f l, p) 302 t := reverse! concat(h, t) 303 p := reverse! concat(f h, p) 304 xRange : R := select(p, xCoord, min) .. select(p, xCoord, max) 305 yRange : R := select(p, yCoord, min) .. select(p, yCoord, max) 306 zRange : R := select(p, zCoord, min) .. select(p, zCoord, max) 307 [f, [tRange, xRange, yRange, zRange], t, p] 308 309 zoom(p, xRange, yRange, zRange) == 310 [[xRange, yRange, zRange], p.bounds, 311 p.screenres, p.axisLabels, p.functions] 312 313 basicRefine(curve, nRange) == 314 tRange : R := first curve.ranges 315 -- curve := copy$C curve -- Yet another @#$%^&* compiler bug 316 curve : C := [curve.source, curve.ranges, curve.knots, curve.points] 317 t := curve.knots := copy curve.knots 318 p := curve.points := copy curve.points 319 l := low(nRange); h := high(nRange) 320 f := curve.source 321 while not(empty?(rest(t))) and first(t) < h repeat 322 second(t) < l => (t := rest t; p := rest p) 323 -- insert new point between t.0 and t.1 324 tm : F := (first(t) + second(t))/2::F 325 -- if DEBUG then output$O (tm::E) 326 pm := f tm 327 NUMFUNEVALS := NUMFUNEVALS + 1 328 t.rest := concat(tm, rest t); t := rest rest t 329 p.rest := concat(pm, rest p); p := rest rest p 330 t := curve.knots; p := curve.points 331 xRange := select(p, xCoord, min) .. select(p, xCoord, max) 332 yRange := select(p, yCoord, min) .. select(p, yCoord, max) 333 zRange := select(p, zCoord, min) .. select(p, zCoord, max) 334 [curve.source, [tRange, xRange, yRange, zRange], t, p] 335 336 refine p == refine(p, parametricRange p) 337 refine(p, nRange) == 338 NUMFUNEVALS := 0 339 tRange := parametricRange p 340 nRange := intersect(tRange, nRange) 341 curves : L C := [basicRefine(c, nRange) for c in p.functions] 342 xRange := join(curves, 1); yRange := join(curves, 2) 343 zRange := join(curves, 3) 344 scrres := p.screenres 345 if adaptive3D? then 346 tlimit := 8 347 curves := [adaptivePlot(c, nRange, xRange, yRange, zRange, _ 348 tlimit, scrres := 2*scrres) for c in curves] 349 xRange := join(curves, 1); yRange := join(curves, 2) 350 zRange := join(curves, 3) 351 [p.display, [tRange, xRange, yRange, zRange], _ 352 scrres, p.axisLabels, curves] 353 354 plot(p : %, tRange : R) == 355 -- re plot p on a new range making use of the points already 356 -- computed if possible 357 NUMFUNEVALS := 0 358 curves : L C := [rangeRefine(c, tRange) for c in p.functions] 359 xRange := join(curves, 1); yRange := join(curves, 2) 360 zRange := join(curves, 3) 361 if adaptive3D? then 362 tlimit := 8 363 curves := [adaptivePlot(c, tRange, xRange, yRange, zRange, tlimit, _ 364 p.screenres) for c in curves] 365 xRange := join(curves, 1); yRange := join(curves, 2) 366 zRange := join(curves, 3) 367-- print(NUMFUNEVALS::OUT) 368 [[xRange, yRange, zRange], [tRange, xRange, yRange, zRange], 369 p.screenres, p.axisLabels, curves] 370 371 pointPlot(f : F -> P, tRange : R) == 372 p := basicPlot(f, tRange) 373 r := p.ranges 374 NUMFUNEVALS := MINPOINTS 375 if adaptive3D? then 376 p := adaptivePlot(p, first r, second r, third r, fourth r, 8, SCREENRES) 377-- print(NUMFUNEVALS::OUT) 378-- print(p::OUT) 379 [rest(r), r, SCREENRES, [], [p]] 380 381 pointPlot(f : F -> P, tRange : R, xRange : R, yRange : R, zRange : R) == 382 p := pointPlot(f, tRange) 383 p.display:= [checkRange xRange, checkRange yRange, checkRange zRange] 384 p 385 386 myTrap : (F-> F, F) -> F 387 myTrap(ff : F-> F, f : F) : F == 388 s := trapNumericErrors(ff(f))$Lisp :: Union(F, "failed") 389 if (s) case "failed" then 390 r : F := 0 391 else 392 r : F := s 393 r 394 395 plot(f1 : F -> F, f2 : F -> F, f3 : F -> F, col : F -> F, tRange : R) == 396 p := basicPlot((z1 : F) : P +-> point(myTrap(f1, z1), myTrap(f2, z1), myTrap(f3, z1), col(z1)), tRange) 397 r := p.ranges 398 NUMFUNEVALS := MINPOINTS 399 if adaptive3D? then 400 p := adaptivePlot(p, first r, second r, third r, fourth r, 8, SCREENRES) 401-- print(NUMFUNEVALS::OUT) 402 [rest(r), r, SCREENRES, [], [p]] 403 404 plot(f1 : F -> F, f2 : F -> F, f3 : F -> F, col : F -> F, _ 405 tRange : R, xRange : R, yRange : R, zRange : R) == 406 p := plot(f1, f2, f3, col, tRange) 407 p.display:= [checkRange xRange, checkRange yRange, checkRange zRange] 408 p 409 410--% terminal output 411 412 coerce r == 413 spaces := message(" ") 414 xSymbol := message("x = ") 415 ySymbol := message("y = ") 416 zSymbol := message("z = ") 417 tSymbol := message("t = ") 418 tRange := (parametricRange r) :: OUT 419 f : L OUT := [] 420 for curve in r.functions repeat 421 xRange := coerce curve.ranges.1 422 yRange := coerce curve.ranges.2 423 zRange := coerce curve.ranges.3 424 l : L OUT := [xSymbol, xRange, spaces, ySymbol, yRange, _ 425 spaces, zSymbol, zRange] 426 l := concat!([tSymbol, tRange, spaces], l) 427 h : OUT := hconcat l 428 l := [p::OUT for p in curve.points] 429 f := concat(vconcat concat(h, l), f) 430 prefix(message("PLOT" ), reverse! f) 431 432----% graphics output 433 434 listBranches plot == 435 outList : L L P := [] 436 for curve in plot.functions repeat 437 outList := concat(curve.points, outList) 438 outList 439 440--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. 441--All rights reserved. 442-- 443--Redistribution and use in source and binary forms, with or without 444--modification, are permitted provided that the following conditions are 445--met: 446-- 447-- - Redistributions of source code must retain the above copyright 448-- notice, this list of conditions and the following disclaimer. 449-- 450-- - Redistributions in binary form must reproduce the above copyright 451-- notice, this list of conditions and the following disclaimer in 452-- the documentation and/or other materials provided with the 453-- distribution. 454-- 455-- - Neither the name of The Numerical ALgorithms Group Ltd. nor the 456-- names of its contributors may be used to endorse or promote products 457-- derived from this software without specific prior written permission. 458-- 459--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS 460--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED 461--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A 462--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER 463--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 464--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 465--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 466--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 467--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 468--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 469--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 470