diff options
author | dos-reis <gdr@axiomatics.org> | 2007-08-14 05:14:52 +0000 |
---|---|---|
committer | dos-reis <gdr@axiomatics.org> | 2007-08-14 05:14:52 +0000 |
commit | ab8cc85adde879fb963c94d15675783f2cf4b183 (patch) | |
tree | c202482327f474583b750b2c45dedfc4e4312b1d /src/algebra/plot3d.spad.pamphlet | |
download | open-axiom-ab8cc85adde879fb963c94d15675783f2cf4b183.tar.gz |
Initial population.
Diffstat (limited to 'src/algebra/plot3d.spad.pamphlet')
-rw-r--r-- | src/algebra/plot3d.spad.pamphlet | 541 |
1 files changed, 541 insertions, 0 deletions
diff --git a/src/algebra/plot3d.spad.pamphlet b/src/algebra/plot3d.spad.pamphlet new file mode 100644 index 00000000..804a1207 --- /dev/null +++ b/src/algebra/plot3d.spad.pamphlet @@ -0,0 +1,541 @@ +\documentclass{article} +\usepackage{axiom} +\begin{document} +\title{\$SPAD/src/algebra plot3d.spad} +\author{Clifton J. Williamson, Michael Monagan, Jon Steinbach} +\maketitle +\begin{abstract} +\end{abstract} +\eject +\tableofcontents +\eject +\section{domain PLOT3D Plot3D} +<<domain PLOT3D Plot3D>>= +)abbrev domain PLOT3D Plot3D +++ Author: Clifton J. Williamson based on code by Michael Monagan +++ Date Created: Jan 1989 +++ Date Last Updated: 22 November 1990 (Jon Steinbach) +++ Basic Operations: pointPlot, plot, zoom, refine, tRange, tValues, +++ minPoints3D, setMinPoints3D, maxPoints3D, setMaxPoints3D, +++ screenResolution3D, setScreenResolution3D, adaptive3D?, setAdaptive3D, +++ numFunEvals3D, debug3D +++ Related Constructors: +++ Also See: +++ AMS Classifications: +++ Keywords: plot, parametric +++ References: +++ Description: Plot3D supports parametric plots defined over a real +++ number system. A real number system is a model for the real +++ numbers and as such may be an approximation. For example, +++ floating point numbers and infinite continued fractions are +++ real number systems. The facilities at this point are limited +++ to 3-dimensional parametric plots. +Plot3D(): Exports == Implementation where + B ==> Boolean + F ==> DoubleFloat + I ==> Integer + L ==> List + N ==> NonNegativeInteger + OUT ==> OutputForm + P ==> Point F + S ==> String + R ==> Segment F + O ==> OutputPackage + C ==> Record(source: F -> P,ranges: L R, knots: L F, points: L P) + + Exports ==> PlottableSpaceCurveCategory with + + pointPlot: (F -> P,R) -> % + ++ pointPlot(f,g,h,a..b) plots {/emx = f(t), y = g(t), z = h(t)} as + ++ t ranges over {/em[a,b]}. + pointPlot: (F -> P,R,R,R,R) -> % + ++ pointPlot(f,x,y,z,w) \undocumented + plot: (F -> F,F -> F,F -> F,F -> F,R) -> % + ++ plot(f,g,h,a..b) plots {/emx = f(t), y = g(t), z = h(t)} as + ++ t ranges over {/em[a,b]}. + plot: (F -> F,F -> F,F -> F,F -> F,R,R,R,R) -> % + ++ plot(f1,f2,f3,f4,x,y,z,w) \undocumented + + plot: (%,R) -> % -- change the range + ++ plot(x,r) \undocumented + zoom: (%,R,R,R) -> % + ++ zoom(x,r,s,t) \undocumented + refine: (%,R) -> % + ++ refine(x,r) \undocumented + refine: % -> % + ++ refine(x) \undocumented + + tRange: % -> R + ++ tRange(p) returns the range of the parameter in a parametric plot p. + tValues: % -> L L F + ++ tValues(p) returns a list of lists of the values of the parameter for + ++ which a point is computed, one list for each curve in the plot p. + + minPoints3D: () -> I + ++ minPoints3D() returns the minimum number of points in a plot. + setMinPoints3D: I -> I + ++ setMinPoints3D(i) sets the minimum number of points in a plot to i. + maxPoints3D: () -> I + ++ maxPoints3D() returns the maximum number of points in a plot. + setMaxPoints3D: I -> I + ++ setMaxPoints3D(i) sets the maximum number of points in a plot to i. + screenResolution3D: () -> I + ++ screenResolution3D() returns the screen resolution for a 3d graph. + setScreenResolution3D: I -> I + ++ setScreenResolution3D(i) sets the screen resolution for a 3d graph to i. + adaptive3D?: () -> B + ++ adaptive3D?() determines whether plotting be done adaptively. + setAdaptive3D: B -> B + ++ setAdaptive3D(true) turns adaptive plotting on; + ++ setAdaptive3D(false) turns adaptive plotting off. + numFunEvals3D: () -> I + ++ numFunEvals3D() returns the number of points computed. + debug3D: B -> B + ++ debug3D(true) turns debug mode on; + ++ debug3D(false) turns debug mode off. + + Implementation ==> add + import PointPackage(F) + +--% local functions + + fourth : L R -> R + checkRange : R -> R + -- checks that left-hand endpoint is less than right-hand endpoint + intersect : (R,R) -> R + -- intersection of two intervals + union : (R,R) -> R + -- union of two intervals + join : (L C,I) -> R + parametricRange: % -> R +-- setColor : (P,F) -> F + select : (L P,P -> F,(F,F) -> F) -> F +-- normalizeColor : (P,F,F) -> F + rangeRefine : (C,R) -> C + adaptivePlot : (C,R,R,R,R,I,I) -> C + basicPlot : (F -> P,R) -> C + basicRefine : (C,R) -> C + point : (F,F,F,F) -> P + +--% representation + + Rep := Record( display: L R, _ + bounds: L R, _ + screenres: I, _ + axisLabels: L S, _ + functions: L C ) + +--% global constants + + ADAPTIVE : B := true + MINPOINTS : I := 49 + MAXPOINTS : I := 1000 + NUMFUNEVALS : I := 0 + SCREENRES : I := 500 + ANGLEBOUND : F := cos inv (4::F) + DEBUG : B := false + + point(xx,yy,zz,col) == point(l : L F := [xx,yy,zz,col]) + + fourth list == first rest rest rest list + + checkRange r == (lo r > hi r => error "ranges cannot be negative"; r) + intersect(s,t) == checkRange (max(lo s,lo t) .. min(hi s,hi t)) + union(s:R,t:R) == min(lo s,lo t) .. max(hi s,hi t) + join(l,i) == + rr := first l + u : R := + i = 0 => first(rr.ranges) + i = 1 => second(rr.ranges) + i = 2 => third(rr.ranges) + fourth(rr.ranges) + for r in rest l repeat + i = 0 => union(u,first(r.ranges)) + i = 1 => union(u,second(r.ranges)) + i = 2 => union(u,third(r.ranges)) + union(u,fourth(r.ranges)) + u + parametricRange r == first(r.bounds) + + minPoints3D() == MINPOINTS + setMinPoints3D n == + if n < 3 then error "three points minimum required" + if MAXPOINTS < n then MAXPOINTS := n + MINPOINTS := n + maxPoints3D() == MAXPOINTS + setMaxPoints3D n == + if n < 3 then error "three points minimum required" + if MINPOINTS > n then MINPOINTS := n + MAXPOINTS := n + screenResolution3D() == SCREENRES + setScreenResolution3D n == + if n < 2 then error "buy a new terminal" + SCREENRES := n + adaptive3D?() == ADAPTIVE + setAdaptive3D b == ADAPTIVE := b + + numFunEvals3D() == NUMFUNEVALS + debug3D b == DEBUG := b + +-- setColor(p,c) == p.colNum := c + + xRange plot == second plot.bounds + yRange plot == third plot.bounds + zRange plot == fourth plot.bounds + tRange plot == first plot.bounds + + tValues plot == + outList : L L F := nil() + for curve in plot.functions repeat + outList := concat(curve.knots,outList) + outList + + select(l,f,g) == + m := f first l + if (EQL(m, _$NaNvalue$Lisp)$Lisp) then m := 0 +-- for p in rest l repeat m := g(m,fp) + for p in rest l repeat + fp : F := f p + if (EQL(fp, _$NaNvalue$Lisp)$Lisp) then fp := 0 + m := g(m,fp) + m + +-- normalizeColor(p,lo,diff) == +-- p.colNum := (p.colNum - lo)/diff + + rangeRefine(curve,nRange) == + checkRange nRange; l := lo nRange; h := hi nRange + t := curve.knots; p := curve.points; f := curve.source + while not null t and first t < l repeat + (t := rest t; p := rest p) + c : L F := nil(); q : L P := nil() + while not null t and first t <= h repeat + c := concat(first t,c); q := concat(first p,q) + t := rest t; p := rest p + if null c then return basicPlot(f,nRange) + if first c < h then + c := concat(h,c); q := concat(f h,q) + NUMFUNEVALS := NUMFUNEVALS + 1 + t := c := reverse_! c; p := q := reverse_! q + s := (h-l)/(MINPOINTS::F-1) + if (first t) ^= l then + t := c := concat(l,c); p := q := concat(f l,p) + NUMFUNEVALS := NUMFUNEVALS + 1 + while not null rest t repeat + n := wholePart((second(t) - first(t))/s) + d := (second(t) - first(t))/((n+1)::F) + for i in 1..n repeat + t.rest := concat(first(t) + d,rest t); t1 := second t + p.rest := concat(f t1,rest p) + NUMFUNEVALS := NUMFUNEVALS + 1 + t := rest t; p := rest p + t := rest t + p := rest p + xRange := select(q,xCoord,min) .. select(q,xCoord,max) + yRange := select(q,yCoord,min) .. select(q,yCoord,max) + zRange := select(q,zCoord,min) .. select(q,zCoord,max) +-- colorLo := select(q,color,min); colorHi := select(q,color,max) +-- (diff := colorHi - colorLo) = 0 => +-- error "all points are the same color" +-- map(normalizeColor(#1,colorLo,diff),q)$ListPackage1(P) + [f,[nRange,xRange,yRange,zRange],c,q] + + + adaptivePlot(curve,tRg,xRg,yRg,zRg,pixelfraction,resolution) == + xDiff := hi xRg - lo xRg + yDiff := hi yRg - lo yRg + zDiff := hi zRg - lo zRg +-- xDiff = 0 or yDiff = 0 or zDiff = 0 => curve--!! delete this? + if xDiff = 0::F then xDiff := 1::F + if yDiff = 0::F then yDiff := 1::F + if zDiff = 0::F then zDiff := 1::F + l := lo tRg; h := hi tRg + (tDiff := h-l) = 0 => curve + t := curve.knots + #t < 3 => curve + p := curve.points; f := curve.source + minLength:F := 4::F/resolution::F + maxLength := 1/4::F + tLimit := tDiff/(pixelfraction*resolution)::F + while not null t and first t < l repeat (t := rest t; p := rest p) + #t < 3 => curve + headert := t; headerp := p + st := t; sp := p + todot : L L F := nil() + todop : L L P := nil() + while not null rest rest st repeat + todot := concat_!(todot, st) + todop := concat_!(todop, sp) + st := rest st; sp := rest sp + st := headert; sp := headerp + todo1 := todot; todo2 := todop + n : I := 0 + + while not null todo1 repeat + st := first(todo1) + t0 := first(st); t1 := second(st); t2 := third(st) + if t2 > h then leave + t2 - t0 < tLimit => + todo1 := rest todo1 + todo2 := rest todo2; + if not null todo1 then (t := first(todo1); p := first(todo2)) + sp := first(todo2) + x0 := xCoord first(sp); y0 := yCoord first(sp); z0 := zCoord first(sp) + x1 := xCoord second(sp); y1 := yCoord second(sp); z1 := zCoord second(sp) + x2 := xCoord third(sp); y2 := yCoord third(sp); z2 := zCoord third(sp) + a1 := (x1-x0)/xDiff; b1 := (y1-y0)/yDiff; c1 := (z1-z0)/zDiff + a2 := (x2-x1)/xDiff; b2 := (y2-y1)/yDiff; c2 := (z2-z1)/zDiff + s1 := sqrt(a1**2+b1**2+c1**2); s2 := sqrt(a2**2+b2**2+c2**2) + dp := a1*a2+b1*b2+c1*c2 + s1 < maxLength and s2 < maxLength and _ + (s1 = 0 or s2 = 0 or + s1 < minLength and s2 < minLength or _ + dp/s1/s2 > ANGLEBOUND) => + todo1 := rest todo1 + todo2 := rest todo2 + if not null todo1 then (t := first(todo1); p := first(todo2)) + if n = MAXPOINTS then leave else n := n + 1 + --if DEBUG then + --r : L F := [minLength,maxLength,s1,s2,dp/s1/s2,ANGLEBOUND] + --output(r::E)$O + st := rest t + if not null rest rest st then + tm := (t0+t1)/2::F + tj := tm + t.rest := concat(tj,rest t) + p.rest := concat(f tj, rest p) + todo1 := concat_!(todo1, t) + todo2 := concat_!(todo2, p) + t := rest t; p := rest p + todo1 := concat_!(todo1, t) + todo2 := concat_!(todo2, p) + t := rest t; p := rest p + todo1 := rest todo1; todo2 := rest todo2 + + tm := (t1+t2)/2::F + tj := tm + t.rest := concat(tj, rest t) + p.rest := concat(f tj, rest p) + todo1 := concat_!(todo1, t) + todo2 := concat_!(todo2, p) + t := rest t; p := rest p + todo1 := concat_!(todo1, t) + todo2 := concat_!(todo2, p) + todo1 := rest todo1; todo2 := rest todo2 + if not null todo1 then (t := first(todo1); p := first(todo2)) + else + tm := (t0+t1)/2::F + tj := tm + t.rest := concat(tj,rest t) + p.rest := concat(f tj, rest p) + todo1 := concat_!(todo1, t) + todo2 := concat_!(todo2, p) + t := rest t; p := rest p + todo1 := concat_!(todo1, t) + todo2 := concat_!(todo2, p) + t := rest t; p := rest p + + tm := (t1+t2)/2::F + tj := tm + t.rest := concat(tj, rest t) + p.rest := concat(f tj, rest p) + todo1 := concat_!(todo1, t) + todo2 := concat_!(todo2, p) + todo1 := rest todo1; todo2 := rest todo2 + if not null todo1 then (t := first(todo1); p := first(todo2)) + if n > 0 then + NUMFUNEVALS := NUMFUNEVALS + n + t := curve.knots; p := curve.points + xRg := select(p,xCoord,min) .. select(p,xCoord,max) + yRg := select(p,yCoord,min) .. select(p,yCoord,max) + zRg := select(p,zCoord,min) .. select(p,zCoord,max) + [curve.source,[tRg,xRg,yRg,zRg],t,p] + else curve + + basicPlot(f,tRange) == + checkRange tRange; l := lo tRange; h := hi tRange + t : L F := list l; p : L P := list f l + s := (h-l)/(MINPOINTS-1)::F + for i in 2..MINPOINTS-1 repeat + l := l+s; t := concat(l,t) + p := concat(f l,p) + t := reverse_! concat(h,t) + p := reverse_! concat(f h,p) + xRange : R := select(p,xCoord,min) .. select(p,xCoord,max) + yRange : R := select(p,yCoord,min) .. select(p,yCoord,max) + zRange : R := select(p,zCoord,min) .. select(p,zCoord,max) + [f,[tRange,xRange,yRange,zRange],t,p] + + zoom(p,xRange,yRange,zRange) == + [[xRange,yRange,zRange],p.bounds, + p.screenres,p.axisLabels,p.functions] + + basicRefine(curve,nRange) == + tRange:R := first curve.ranges + -- curve := copy$C curve -- Yet another @#$%^&* compiler bug + curve: C := [curve.source,curve.ranges,curve.knots,curve.points] + t := curve.knots := copy curve.knots + p := curve.points := copy curve.points + l := lo nRange; h := hi nRange + f := curve.source + while not null rest t and first(t) < h repeat + second(t) < l => (t := rest t; p := rest p) + -- insert new point between t.0 and t.1 + tm:F := (first(t) + second(t))/2::F + -- if DEBUG then output$O (tm::E) + pm := f tm + NUMFUNEVALS := NUMFUNEVALS + 1 + t.rest := concat(tm,rest t); t := rest rest t + p.rest := concat(pm,rest p); p := rest rest p + t := curve.knots; p := curve.points + xRange := select(p,xCoord,min) .. select(p,xCoord,max) + yRange := select(p,yCoord,min) .. select(p,yCoord,max) + zRange := select(p,zCoord,min) .. select(p,zCoord,max) + [curve.source,[tRange,xRange,yRange,zRange],t,p] + + refine p == refine(p,parametricRange p) + refine(p,nRange) == + NUMFUNEVALS := 0 + tRange := parametricRange p + nRange := intersect(tRange,nRange) + curves: L C := [basicRefine(c,nRange) for c in p.functions] + xRange := join(curves,1); yRange := join(curves,2) + zRange := join(curves,3) + scrres := p.screenres + if adaptive3D? then + tlimit := 8 + curves := [adaptivePlot(c,nRange,xRange,yRange,zRange, _ + tlimit,scrres := 2*scrres) for c in curves] + xRange := join(curves,1); yRange := join(curves,2) + zRange := join(curves,3) + [p.display,[tRange,xRange,yRange,zRange], _ + scrres,p.axisLabels,curves] + + plot(p:%,tRange:R) == + -- re plot p on a new range making use of the points already + -- computed if possible + NUMFUNEVALS := 0 + curves: L C := [rangeRefine(c,tRange) for c in p.functions] + xRange := join(curves,1); yRange := join(curves,2) + zRange := join(curves,3) + if adaptive3D? then + tlimit := 8 + curves := [adaptivePlot(c,tRange,xRange,yRange,zRange,tlimit, _ + p.screenres) for c in curves] + xRange := join(curves,1); yRange := join(curves,2) + zRange := join(curves,3) +-- print(NUMFUNEVALS::OUT) + [[xRange,yRange,zRange],[tRange,xRange,yRange,zRange], + p.screenres,p.axisLabels,curves] + + pointPlot(f:F -> P,tRange:R) == + p := basicPlot(f,tRange) + r := p.ranges + NUMFUNEVALS := MINPOINTS + if adaptive3D? then + p := adaptivePlot(p,first r,second r,third r,fourth r,8,SCREENRES) +-- print(NUMFUNEVALS::OUT) +-- print(p::OUT) + [ rest r, r, SCREENRES, nil(), [ p ] ] + + pointPlot(f:F -> P,tRange:R,xRange:R,yRange:R,zRange:R) == + p := pointPlot(f,tRange) + p.display:= [checkRange xRange,checkRange yRange,checkRange zRange] + p + + myTrap: (F-> F, F) -> F + myTrap(ff:F-> F, f:F):F == + s := trapNumericErrors(ff(f))$Lisp :: Union(F, "failed") + if (s) case "failed" then + r:F := _$NaNvalue$Lisp + else + r:F := s + r + + plot(f1:F -> F,f2:F -> F,f3:F -> F,col:F -> F,tRange:R) == + p := basicPlot(point(myTrap(f1,#1),myTrap(f2,#1),myTrap(f3,#1),col(#1)),tRange) + r := p.ranges + NUMFUNEVALS := MINPOINTS + if adaptive3D? then + p := adaptivePlot(p,first r,second r,third r,fourth r,8,SCREENRES) +-- print(NUMFUNEVALS::OUT) + [ rest r, r, SCREENRES, nil(), [ p ] ] + + plot(f1:F -> F,f2:F -> F,f3:F -> F,col:F -> F,_ + tRange:R,xRange:R,yRange:R,zRange:R) == + p := plot(f1,f2,f3,col,tRange) + p.display:= [checkRange xRange,checkRange yRange,checkRange zRange] + p + +--% terminal output + + coerce r == + spaces := " " :: OUT + xSymbol := "x = " :: OUT; ySymbol := "y = " :: OUT + zSymbol := "z = " :: OUT; tSymbol := "t = " :: OUT + tRange := (parametricRange r) :: OUT + f : L OUT := nil() + for curve in r.functions repeat + xRange := coerce curve.ranges.1 + yRange := coerce curve.ranges.2 + zRange := coerce curve.ranges.3 + l : L OUT := [xSymbol,xRange,spaces,ySymbol,yRange,_ + spaces,zSymbol,zRange] + l := concat_!([tSymbol,tRange,spaces],l) + h : OUT := hconcat l + l := [p::OUT for p in curve.points] + f := concat(vconcat concat(h,l),f) + prefix("PLOT" :: OUT,reverse_! f) + +----% graphics output + + listBranches plot == + outList : L L P := nil() + for curve in plot.functions repeat + outList := concat(curve.points,outList) + outList + +@ +\section{License} +<<license>>= +--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. +--All rights reserved. +-- +--Redistribution and use in source and binary forms, with or without +--modification, are permitted provided that the following conditions are +--met: +-- +-- - Redistributions of source code must retain the above copyright +-- notice, this list of conditions and the following disclaimer. +-- +-- - Redistributions in binary form must reproduce the above copyright +-- notice, this list of conditions and the following disclaimer in +-- the documentation and/or other materials provided with the +-- distribution. +-- +-- - Neither the name of The Numerical ALgorithms Group Ltd. nor the +-- names of its contributors may be used to endorse or promote products +-- derived from this software without specific prior written permission. +-- +--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER +--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +@ +<<*>>= +<<license>> + +<<domain PLOT3D Plot3D>> +@ +\eject +\begin{thebibliography}{99} +\bibitem{1} nothing +\end{thebibliography} +\end{document} |