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