intdraw2010.html

> restart;
with(plots):

Checking set-up of iterated integrals

About this worksheet

Author and Date

Shannon Holland and Matthias Kawski

kawski@asu.edu

http://math.asu.edu/~kawski

original version MAPLE V, release 3, about 1995.

The original version is introduced and discussed in detail to the journal article

" 3D-graphics for iterated integrals" , MAPLE-Tech, 4 no.1, (1997) pp. 92-98.
http://math.asu.edu/~kawski/preprints/97mapletech.pdf

Supported in part by the National Science Foundation through EEC 92-21460 : The Foundation Coalition.
All rights reserved.

>

Content, Purpose and Use

This code does NOT use computer algebra -- but it is convenient to be in the same language
in which iterated integrals may be evaluated.

The procedures intdraw and intdraw3d were originally introduced as the MAPLE V release 4
package "asu". Long ago that package structure vecome defunct, and the worksheet now defines
the requisite functions in the first section -- just execute it in toto by hitting ENTER a few times.

The procedures intdraw and intdraw3d serve to visualize two- and three dimensional regions

{ (u,v,w) : u0 < u < u1, v0(u) < v < v1(u), w0(u,v) < w < w1(u,v) } defined by inequalities of

the form encountered in iterated integrals. The regions are sliced or represented by columns,
as indicated by the user, in a way that helps conceptualize the mental images behind iterated
integrals. The procedures recognize Cartesian, polar, cylindrical, and spherical coordinates
from the standard names x, y, z, r, theta, rho, phi, and they provide detailed feedback for
common mistakes (e.g. ui not being constant).

The main purpose is to give immediate feedback to students, with enough detail to help them
identify and correct the wrong limits. Unlike a popular package written at TX A&M package,
these procedures "slice" the regions, and thus give much more precise feeback.
Ultimately, students are expected to "work until the puicture is right", i.e. no grading is

required. To make this work, it is important that regions are defined for which the user
has a mental image (or an outline of the shape should be given).

Note that the order of arguments is the same one as the usual order from left to right (outside to
inside) in an iterated integral. To evaluate the iterated integral using nested calls to int(...) the order
of the limits must be reversed! int(int(int(f(u,v,w),w=w0(u,v)..w1(u,v)),v=v0(u)..v1(u)),u=u0..u1);

Updates and log of modifications.

The original package held up extremely well and was still running in release 7. In 2010, we

removed all references to the original "package" structure, and it is still running like a charm.
Also in March 2010, we added some sectioning, this intro section, and we reformatted the
input with usual indentations. However, the code itself is virtually unchanged from 1995.
We welcome any proficient MAPLE user to modernize notation, and suggest more
efficient coding. Also, we welcome more examples.

For a description of optional paramters such as grid, spacing, and preferences such as
slices see the original help section -- still in the old "package" format.

The examples are largely the original ones chosen by Shannon, although we updated
some choices for finer partitons as today's compouers are a little faster.
Eventually, the default grid sizes in the programs should be increased, too.

Mar 7, 2010: Updated examples.

>

Execute this section to define intdraw/intdraw3d

MakeGrid2D creates centerpoints for columns and is called from intdraw

>    MakeGrid2D:=proc(a,a0,a1,b,b0,b1,CountA,CountB,SpaceA,
                 SpaceB,Coords,g,CountBArg,SpaceBArg,
                 NASpaceing,NBSpaceing)

> local i, j, k, aTable, bTable, b0Table, b1Table,
bMin, bMax, CenterGrid, DeltaA, DeltaB:

>    if evalf(a0)>evalf(a1)
   then ERROR(`lower limit exceeds upper limit`)
   fi:

>    if evalf(a0)=evalf(a1)
   then ERROR(`lower limit equal to upper limit`)
   fi:

>    if NASpacing = true or CountA = 1
   then
   DeltaA:=evalf((a1-a0)/(CountA))
   else
   DeltaA:=evalf((a1-a0)/(CountA-SpaceA))
   fi:

>    for k to (CountA) do
    aTable[k]:=evalf(a0+(k-SpaceA/2)*DeltaA-DeltaA/2):
    od:

>    if type(evalf(subs(a=aTable[1],b0)),constant)<>true
   or
   type(evalf(subs(a=aTable[1],b1)),constant)<>true
   then
   ERROR(`invalid boundaries`)
   fi:

>    for k to (CountA) do
    b0Table[k]:=evalf(subs(a=aTable[k],b0)):
    b1Table[k]:=evalf(subs(a=aTable[k],b1)):
    od:

> bMax:=max(seq(op(op(k,[entries(b1Table)])),
k=1..nops([entries(b1Table)])),

> evalf(subs(a=a0,b1)),evalf(subs(a=a1,b1))):

> bMin:=min(seq(op(op(k,[entries(b0Table)])),
k=1..nops([entries(b0Table)])),

> evalf(subs(a=a0,b0)),evalf(subs(a=a1,b0))):

>    if bMin>bMax
   then ERROR(`lower limit exceeds upper limit`)
   fi:

>    if NBSpacing = true or CountB=1
   then
   DeltaB:=evalf((bMax-bMin)/(CountB))
   else
   DeltaB:=evalf((bMax-bMin)/(CountB-SpaceB))
   fi:

>    for j to (CountB)
    do
    bTable[j]:=evalf(bMin+(j-SpaceB/2)*DeltaB-DeltaB/2)
    od:

>    if op(2,g) = theta or SpaceBArg = true or CountBArg = true
   then
   for k to (CountA) do
       for j to (CountB) do
           if (b0Table[k]<bTable[j])and(bTable[j]<b1Table[k])
              then
              CenterGrid[k,j]:=aTable[k],bTable[j]
              fi
           od
       od:
   else
   for k to (CountA) do
         CenterGrid[k]:=aTable[k]
       od
   fi:

>    if eval(CenterGrid)= CenterGrid
   then ERROR(`invalid boundaries`)
   fi:

> eval(CenterGrid),DeltaA,DeltaB:

> end:

MakeCorners2D takes centerpoints and creates corners of boxes. It is called from intdraw.

>    MakeCorners2D:=proc(CentersTable, a, b, b0, b1, SpaceA, SpaceB,
      DeltaA, DeltaB, xPlace, yPlace, g, CountBArg, SpaceBArg)
local f,i,ctr,u0,u1,w0,w1,CornersTemp:
for i from 1 to nops([entries(CentersTable)])
    do
    ctr:=op(i,[entries(CentersTable)]):
    u0:=evalf(op(1, ctr)-(1-SpaceA)*DeltaA/2):
    u1:=evalf(op(1, ctr)+(1-SpaceA)*DeltaA/2):
    if op(2,g)<>theta and CountBArg <> true
       and SpaceBArg <> true
       then
       w0:=subs(a=((u1+u0)/2),b0):
       w1:=subs(a=((u1+u0)/2),b1):
       if type(w0,constant)<> true or type(w1,constant)<> true
          then
          ERROR(`Second/third boundary error`)
          fi:
       else
       w0:=eval(op(2,ctr)-(1-SpaceB)*DeltaB/2):
       w1:=eval(op(2,ctr)+(1-SpaceB)*DeltaB/2):
       fi:
    f:=[[u0,u1],[w0,w1]]:
    CornersTemp[i]:=[[op(1,op(xPlace,f)),op(2,op(yPlace,f))],
                     [op(2,op(xPlace,f)),op(2,op(yPlace,f))],
                     [op(1,op(xPlace,f)),op(1,op(yPlace,f))],
                     [op(2,op(xPlace,f)),op(1,op(yPlace,f))]]
    od:
eval(CornersTemp)
end:

MakePoly2D takes corners and creates the 6 sides of the boxes. called from intdraw.

>    MakePoly2D:=proc(Corners)
local TempPolys,i:
for i to nops([entries(Corners)]) do
    TempPolys[i]:=[op(1, Corners[i]),op(2,Corners[i]),
                   op(4,Corners[i]),op(3,Corners[i])]:
   od:
eval(TempPolys)
end:

ChangeCoords2D changes coordinates from Cartesian to Coords called by `intdraw`

>    ChangeCoords2D:=proc(Corners,Coords)
local j,k,CornersTemp:
if Coords=cylindrical
   then
   for k to nops([entries(Corners)]) do
       CornersTemp[k]:=[seq(
           [op(1,op(j,Corners[k]))*cos(op(2,op(j,Corners[k]))),
            op(1,op(j,Corners[k]))*sin(op(2,op(j,Corners[k])))],
       j = 1..4)]
       od:
   elif Coords=Cartesian
      then
      CornersTemp:=eval(Corners)
   fi:
eval(CornersTemp)
end:

intdraw : Main Program, draws the sliced 2-D regions.

>    intdraw:=proc()
local a ,a0, a1,b ,b0,b1, CountA,CountB,SpaceA, SpaceB,TempArgs,
   CentersTable,CornersTable,Coords,TempTable,PolygonsTable,
   PlotTitle,i,j,gg,DeltaA,DeltaB,xPlace,yPlace,Axes,Scaling,Style,
   xTick,yTick,CountArgs,SpaceArgs,SpaceBArg, CountBArg,TitleArgs,
   ScalingArgs,kk,Options,ggg,NASpacing, NBSpacing:

> ScalingArgs:=false: CountArgs:=false:SpaceArgs:=false:TitleArgs:=false:

> SpaceBArg:=false: CountBArg:Options:=[]:NASpacing:=false:

> NBSpacing :=false:

>    if nargs < 2
   then
   ERROR(`Boundaries must be entered as first 2 arguments`)
   fi:

>    if type(args[1],`=`)=false or
   type(args[2],`=`)=false or
   type(op(2,args[1]),range)=false or
   type(op(2,args[2]),range)=false or
   type(op(2,args[2]),range)=false
   then
   ERROR(`Boundaries must be written as var=lower limit..upper limit`)
   fi:
a:=op(1,args[1]):
a0:=evalf(op(1, op(2,args[1]))):
a1:=evalf(op(2, op(2,args[1]))):
b:=op(1,args[2]):
b0:=evalf(op(1, op(2,args[2]))):
b1:=evalf(op(2, op(2,args[2]))):
if nargs>2
   then
   for kk from 3 to nargs do
       if op(1,args[kk]) = grid
          then
          CountArgs:=true:
          CountA:=round(evalf(op(1, op(2,args[kk])))):
          if nops(op(2,args[kk]))>1
             then
             CountBArg := true:
             CountB:=round(evalf(op(2, op(2,args[kk])))):
             else CountB:=15
             fi:
          elif op(1,args[kk]) = spaces
             then
             SpaceArgs:=true:
             SpaceA:=evalf(op(1, op(2,args[kk]))):
             if nops(op(2,args[kk]))>1
                then
                SpaceBArg:=true:
                SpaceB:=evalf(op(2, op(2,args[kk]))):
                else SpaceB:=.3
                fi:
             elif op(1,args[kk])=title
                  then
                  TitleArgs:=true:
                  PlotTitle :=op(2,args[kk])
             elif op(1,args[kk])=scaling
                 then
                 ScalingArgs:=true:
                 Scaling:=op(2,args[kk])
                 else
                 Options:=[op(Options),args[kk]]
           fi
       od
   fi:
if TitleArgs=false
   then PlotTitle:=cat(a,` `,b)
   fi:
if ScalingArgs=false
   then Scaling :=CONSTRAINED
   fi:
if CountArgs=false
   then
   CountA:=15:
   CountB:=15
   fi:
if SpaceArgs=false
   then
   SpaceA:=.3:
   SpaceB:=0 fi:
if SpaceA<0 or SpaceA>1 or SpaceB<0 or SpaceB>1
   then
   ERROR(`spaces arguments must be between or equal to 0 and 1`) fi:
ggg:=[a,b,c]:
if member(x,ggg)=true and member(y,ggg)=true
   then
   Coords:=Cartesian:
   if x <> 'x'or y <>'y'
      then
      ERROR(`x, and y are key variables and must be`,
            `cleared before running program, (ex x:='x')`)
      fi:
   member(x,ggg,'xPlace'):
   member(y,ggg,'yPlace')
   elif member(r,ggg)=true and member(theta,ggg)=true
      then
      Coords:=cylindrical:
      member(r,ggg,'xPlace'):
      member(theta,ggg,'yPlace'):
      if r <> 'r' or theta <>'theta'
         then
         ERROR(`r, and theta are key variables and must be`,
               `cleared before running program, (ex r:='r')`)
         fi:
      if yPlace = 1 and evalf(a1-a0)=evalf(2*Pi)
         then NASpacing :=true
         fi:
      if yPlace = 2 and evalf(b1-b0)=evalf(2*Pi)
         then NBSpacing :=true
         fi:
      else ERROR(`Cartesian coords accepts only x, and y variables`,
                 `polar coords accepts only r, and theta variables`)
      fi:
TempTable:=MakeGrid2D(a,a0,a1,b,b0,b1,CountA,CountB,SpaceA,
                      SpaceB,Coords,ggg, CountBArg, SpaceBArg):
CentersTable:=TempTable[1]:DeltaA:=TempTable[2]:
DeltaB:=TempTable[3]:
CornersTable:=MakeCorners2D(CentersTable,a,b,b0,b1,SpaceA,
                            SpaceB,DeltaA,DeltaB,xPlace,yPlace,
                            ggg, CountBArg, SpaceBArg):
CornersTable:=ChangeCoords2D(CornersTable,Coords):
PolygonsTable:=MakePoly2D(CornersTable,CentersTable):
polygonplot({(seq(PolygonsTable[i],
                  i=1 .. nops([entries(CentersTable)])))},
                  title=PlotTitle,scaling = Scaling, op(Options))
end:

MakeGrid creates centerpoints for columns. called from `intdraw3d`

>    MakeGrid:=proc(a, a0, a1, b, b0, b1, c, c0, c1, CountA,
               CountB, CountC, SpaceA, SpaceB, SpaceC, Coords,
               g, NASpacing, NBSpacing, NCSpacing,CSArg)
local i, j, k, aTable, bTable, b0Table, b1Table, cTable,
      c0Table, c1Table, bMin, bMax, cMax, cMin, CenterGrid,
      DeltaA, DeltaB, DeltaC:
if evalf(a0)>evalf(a1)
   then
   ERROR(`lower limit exceeds upper`)
   fi:
if evalf(a0)=evalf(a1)
   then
   ERROR(`lower equal to upper`)
   fi:
if NASpacing = true
   then
   DeltaA:=evalf((a1-a0)/(CountA))
   else
   DeltaA:=evalf((a1-a0)/(CountA-SpaceA))
   fi:
for k to (CountA) do
   aTable[k]:=evalf(a0+(k-SpaceA/2)*DeltaA-DeltaA/2):
   od:
if type(evalf(subs(a=aTable[1],b0)),constant)<>true or
   type(evalf(subs(a=aTable[1],b1)),constant)<>true
   then
   ERROR(`invalid first or second boundary`)
   fi:
for k to (CountA) do
   b0Table[k]:=evalf(subs(a=aTable[k],b0)):
   b1Table[k]:=evalf(subs(a=aTable[k],b1)):
   od:
bMax:=max(seq(op(op(k,[entries(b1Table)])),
              k=1..nops([entries(b1Table)])),
          evalf(subs(a=a0,b1)),
          evalf(subs(a=a1,b1))):
bMin:=min(seq(op(op(k,[entries(b0Table)])),
              k=1..nops([entries(b0Table)])),
          evalf(subs(a=a0,b0)),
          evalf(subs(a=a1,b0))):
if bMin>bMax
   then
   ERROR(`lower limit exceeds upper`)
   fi:
if NBSpacing = true
   then
   DeltaB:=evalf((bMax-bMin)/(CountB))
   else
   DeltaB:=evalf((bMax-bMin)/(CountB-SpaceB))
   fi:
for j to (CountB) do
    bTable[j]:=evalf(bMin+(j-SpaceB/2)*DeltaB-DeltaB/2)
    od:
if op(3,g)<>theta and op(3,g)<>phi and CSArg = false
   then
   DeltaC:=0:
   for k to (CountA) do
       for j to (CountB) do
          if (b0Table[k]<bTable[j])and(bTable[j]<b1Table[k])
              then
              CenterGrid[k,j]:=aTable[k],bTable[j]
              fi
          od
       od:
   else
   if type(evalf(subs(a=aTable[1], b=bTable[1],c0)),constant)<>true or
      type(evalf(subs(a=aTable[1], b=bTable[1],c1)) ,constant) <> true
      then
      ERROR(`invalid boundaries`)
      fi:
   for k to (CountA) do
      for j to (CountB) do
         if (b0Table[k]<bTable[j])and(bTable[j]<b1Table[k])
            then
            c0Table[j,k]:=evalf(subs(a=aTable[k], b=bTable[j],c0)):
            c1Table[j,k]:=evalf(subs(a=aTable[k], b=bTable[j],c1)):
            fi
         od
      od:
   cMax:=max(seq(op(op(k,[entries(c1Table)])),
                 k=1..nops([entries(c1Table)])),
             evalf(subs(a=a0,subs(b=b0,c1))),
             evalf(subs(a=a0,subs(b=b1,c1))),
             evalf(subs(a=a1,subs(b=b0,c1))),
             evalf(subs(a=a1,subs(b=b1,c1)))):
   cMin:=min(seq(op(op(k,[entries(c0Table)])),
                 k=1..nops([entries(c0Table)])),
             evalf(subs(a=a0,subs(b=b0,c0))),
             evalf(subs(a=a1,subs(b=b1,c0))),
             evalf(subs(a=a1,subs(b=b0,c0))),
             evalf(subs(a=a1,subs(b=b1,c0)))):
   if cMin>cMax
      then
      ERROR(`lower limit exceeds upper`)
      fi:
   if NCSpacing = true
      then
      DeltaC:=(cMax-cMin)/(CountC)
      else
      DeltaC:=evalf((cMax-cMin)/(CountC-SpaceC))
      fi:
for k to CountC do
   cTable[k]:=evalf(cMin+(k-SpaceC/2)*DeltaC-DeltaC/2):
   od:
for i to (CountA) do
   for j to (CountB) do
      for k to (CountC) do
         if (b0Table[i]<bTable[j]) and
            (bTable[j]<b1Table[i]) and
            (evalf(subs(a=aTable[i],b=bTable[j],c0))<cTable[k]) and
            (cTable[k]<evalf(subs(a=aTable[i],b=bTable[j],c1)))
            then
            CenterGrid[i,j,k]:=aTable[i],bTable[j],cTable[k]
            fi
         od
      od
   od
fi:
if eval(CenterGrid)= CenterGrid
   then
   ERROR(`invalid boundaries`)
   fi:
eval(CenterGrid),DeltaA,DeltaB,DeltaC:
end:

MakeCorners takes centerpoints and creates the corners of the boxes called from `intdraw3d`

>    MakeCorners:=proc(CentersTable, a, b, c, c0, c1,
                  SpaceA, SpaceB, SpaceC,
                  DeltaA, DeltaB, DeltaC,
                  xPlace, yPlace, zPlace, g, CSArg)
local f, i, ctr, u0, u1, v0, v1, w0, w1,CornersTemp:
for i from 1 to nops([entries(CentersTable)]) do
    ctr:=op(i,[entries(CentersTable)]):
    u0:=evalf(op(1, ctr)-(1-SpaceA)*DeltaA/2):
    u1:=evalf(op(1, ctr)+(1-SpaceA)*DeltaA/2):
    v0:=evalf(op(2, ctr)-(1-SpaceB)*DeltaB/2):
    v1:=evalf(op(2, ctr)+(1-SpaceB)*DeltaB/2):
    if op(3,g)<>theta and op(3,g)<>phi and CSArg = false
       then
       w0:=subs(a=((u1+u0)/2),b=((v1+v0)/2),c0):
       w1:=subs(a=((u1+u0)/2),b=((v1+v0)/2),c1):
       if type(w0,constant)<> true or type(w1,constant)<> true
          then
          ERROR(`Second/third boundry error`)
          fi:
          else
       w0:=eval(op(3,ctr)-(1-SpaceC)*DeltaC/2):
       w1:=eval(op(3,ctr)+(1-SpaceC)*DeltaC/2):
       fi:
     f:=[[u0,u1],[v0,v1],[w0,w1]]:
     CornersTemp[i]:=[
        [op(1,op(xPlace,f)),op(1,op(yPlace,f)),op(1,op(zPlace,f))],
        [op(1,op(xPlace,f)),op(1,op(yPlace,f)),op(2,op(zPlace,f))],
        [op(1,op(xPlace,f)),op(2,op(yPlace,f)),op(1,op(zPlace,f))],
        [op(1,op(xPlace,f)),op(2,op(yPlace,f)),op(2,op(zPlace,f))],
        [op(2,op(xPlace,f)),op(1,op(yPlace,f)),op(1,op(zPlace,f))],
        [op(2,op(xPlace,f)),op(1,op(yPlace,f)),op(2,op(zPlace,f))],
        [op(2,op(xPlace,f)),op(2,op(yPlace,f)),op(1,op(zPlace,f))],
        [op(2,op(xPlace,f)),op(2,op(yPlace,f)),op(2,op(zPlace,f))]]
     od:

> eval(CornersTemp)

>

> end:

MakePoly takes corners and creates the 6 sides of the boxes. called from `intdraw3d`

>    MakePoly:=proc(Corners)
local i, TempPolys:
for i to nops([entries(Corners)]) do
   TempPolys[i,1]:=[op(1, Corners[i]),op(2,Corners[i]),
                    op(4,Corners[i]),op(3,Corners[i])]:
   TempPolys[i,2]:=[op(5, Corners[i]),op(6,Corners[i]),
                    op(8,Corners[i]),op(7,Corners[i])]:
   TempPolys[i,3]:=[op(2, Corners[i]),op(4,Corners[i]),
                    op(8,Corners[i]),op(6,Corners[i])]:
   TempPolys[i,4]:=[op(4, Corners[i]),op(3,Corners[i]),
                    op(7,Corners[i]),op(8,Corners[i])]:
   TempPolys[i,5]:=[op(3, Corners[i]),op(1,Corners[i]),
                    op(5,Corners[i]),op(7,Corners[i])]:
   TempPolys[i,6]:=[op(1, Corners[i]),op(2,Corners[i]),
                    op(6,Corners[i]),op(5,Corners[i])]:
   od:
eval(TempPolys)
end:

ChangeCoords changes coordinates from Cartesian to Coords called by `intdraw3d`

>    ChangeCoords:=proc(Corners,Coords)
local j, k, CornersTemp:
if Coords=cylindrical
   then
   for k to nops([entries(Corners)]) do
       CornersTemp[k]:=[seq([
           op(1,op(j,Corners[k]))*cos(op(2,op(j,Corners[k]))),
           op(1,op(j,Corners[k]))*sin(op(2,op(j,Corners[k]))),
           op(3,op(j,Corners[k]))],j = 1..8)]
       od:
elif Coords=spherical
   then
   for k to nops([entries(Corners)]) do
       CornersTemp[k]:=[seq([
           op(1,op(j,Corners[k]))*sin(op(2,op(j,Corners[k])))*cos(op(3,op(j,Corners[k]))),
           op(1,op(j,Corners[k]))*sin(op(2,op(j,Corners[k])))*sin(op(3,op(j,Corners[k]))),
           op(1,op(j,Corners[k]))*cos(op(2,op(j,Corners[k])))],j = 1..8)]
       od:
elif Coords=Cartesian
   then
   CornersTemp:=eval(Corners)
   fi:
eval(CornersTemp)
end:

intdraw3d Main Program, draws region of 3-D integral

> intdraw3d:=proc()

>    local a ,a0, a1,b ,b0,b1, c, c0, c1, CountA, CountB, CountC,
      SpaceA, SpaceB, SpaceC, TempArgs, CentersTable,
      CornersTable, Coords, TempTable, PolygonsTable,
      PlotTitle, i, j, gg, DeltaA, DeltaB, DeltaC,
      xPlace, yPlace, zPlace, Scaling, Style,
      xTick, yTick, zTick, aTemp, bTemp, cTemp,
      CountArgs, SpaceArgs, TitleArgs, ScalingArgs,
      kk, Options, ggg, NASpacing, NBSpacing, NCSpacing,
      CSArg, Labels, LabelsArgs, Axes, AxesArgs:
ScalingArgs:= false:
CountArgs := false:
SpaceArgs := false:
AxesArgs   := false:
TitleArgs := false:
Options    := []:
NASpacing := false:
NBSpacing := false:
NCSpacing := false:
CSArg      := false:
LabelsArgs := false:
if nargs < 3
   then
   ERROR(`Boundaries must be entered as first 3 arguments`)
   fi:
if type(args[1],`=`)=false or
   type(args[2],`=`)=false or
   type(args[3],`=`)=false or
   type(op(2,args[1]),range)=false or
   type(op(2,args[2]),range)=false or
   type(op(2,args[2]),range)=false
   then
   ERROR(`Boundaries must be written as var=range`)
   fi:
a:=op(1,args[1]):
a0:=evalf(op(1, op(2,args[1]))):
a1:=evalf(op(2, op(2,args[1]))):
b:=op(1,args[2]):
b0:=evalf(op(1, op(2,args[2]))):
b1:=evalf(op(2, op(2,args[2]))):
c:=op(1,args[3]):
c0:=evalf(op(1, op(2,args[3]))):
c1:=evalf(op(2, op(2,args[3]))):
if nargs>3
   then
   for kk from 4 to nargs do
       if op(1,args[kk]) = grid
          then
          CountArgs:=true:
          CountA:=evalf(op(1, op(2,args[kk]))):
          if nops(op(2,args[kk]))>1
             then
             CountB:=evalf(op(2, op(2,args[kk]))):
             else
             CountB:=5:
             fi:
          if nops(op(2,args[kk]))>2
             then
             CountC:=evalf(op(3,op(2,args[kk]))):
             CSArg := true
             else
             CountC:=5:
             fi:
        elif op(1,args[kk]) = spaces
             then
             SpaceArgs:=true:
             if op(1,op(2,args[kk]))= slices
                then
                SpaceA:=.8:
                SpaceB:=.1:
                SpaceC:=0:
                else
                SpaceA:=evalf(op(1, op(2,args[kk]))):
                if nops(op(2,args[kk]))>1
                   then
                   SpaceB:=evalf(op(2, op(2,args[kk]))):
                   else SpaceB:=.3
                   fi:
                if nops(op(2,args[kk]))>2
                   then
                   SpaceC:=evalf(op(3,op(2,args[kk]))):
                   CSArg:=true
                   else
                   SpaceC:=0
                   fi:
                fi:
             elif op(1,args[kk])=title
                then
                TitleArgs:=true:
                PlotTitle :=op(2,args[kk])
             elif op(1,args[kk])=scaling
                then
                ScalingArgs:=true:
                Scaling:=op(2,args[kk])
             elif op(1,args[kk])= axes
                then
                AxesArgs:=true:
                Axes := op(2,args[kk])
             elif op(1,args[kk])= labels
                then
                LabelsArgs := true:
                Labels :=op(2,args[kk])
                else
                Options:=[op(Options),args[kk]]
                fi
             od
          fi:
       if LabelsArgs=false
          then
          Labels:=[`x`,`y`,`z`]
          fi:
if AxesArgs=false
   then
   Axes:=FRAME
   fi:
if TitleArgs=false
   then
   PlotTitle:=cat(a,` `,b,` `,c)
   fi:
if ScalingArgs=false
   then
   Scaling :=CONSTRAINED
   fi:
if CountArgs=false
   then
   CountA:=5:
   CountB:=5:
   CountC:=5
   fi:
if SpaceArgs=false
   then
   SpaceA:=.3:
   SpaceB:=.3:
   SpaceC:=0
   fi:
if SpaceA<0 or SpaceA>1 or SpaceB<0 or
   SpaceB>1 or SpaceC<0 or SpaceC>1
   then
   ERROR(`spaces arguments must range from 0 and 1`)
   fi:
if CountA<=0 or CountB<=0 or CountC<=0
   then
   ERROR(`grid arguments must be integers that are greater than 0`)
   fi:
ggg:=[a,b,c]:
if member(x,ggg)=true and
   member(y,ggg)=true and
   member(z,ggg)=true
   then
   Coords:=Cartesian:
   if x <> 'x'or y <>'y' or z <> 'z'
      then
      ERROR(`x, y, and z are key variables and must be`,
            `cleared before running program, (ex x:='x')`)
      fi:
   member(x,ggg,'xPlace'):
   member(y,ggg,'yPlace'):
   member(z,ggg,'zPlace'):
elif member(r,ggg)=true and
     member(theta,ggg)=true and
     member(z,ggg)= true
   then
   Coords:=cylindrical:
   if theta <> 'theta' or z <> 'z' or r<>'r'
      then
      ERROR(`z, theta, and r are key variables and must be`,
            `cleared before running program, (ex r:='r')`)
      fi:
   member(r,ggg,'xPlace'):
   member(theta,ggg,'yPlace'):
   member(z,ggg,'zPlace'):
if yPlace = 1 and evalf(a1-a0)=evalf(2*Pi)
   then
   NASpacing :=true fi:
if yPlace = 2 and evalf(b1-b0)=evalf(2*Pi)
   then NBSpacing :=true
   fi:
if yPlace = 3 and evalf(c1-c0)=evalf(2*Pi)
   then NCSpacing :=true
   fi:
elif member(rho,ggg)=true and
     member(phi,ggg)=true and
     member(theta,ggg)=true
   then
   Coords:=spherical:
   member(rho,ggg,'xPlace'):
   member(phi,ggg,'yPlace'):
   member(theta,ggg,'zPlace'):
   if theta <> 'theta' or rho <>'rho' or phi <> 'phi'
       then
       ERROR(`theta, rho, and phi are key variables and must be`,
             `cleared before running program, (ex rho:='rho')`)
       fi:
   if (yPlace = 1 or zPlace = 1) and evalf(a1-a0) = evalf(2*Pi)
      then
      NASpacing :=true
      fi:
   if (yPlace = 2 or zPlace = 2) and evalf(b1-b0) =evalf(2*Pi)
      then
      NBSpacing:=true
      fi:
   if (yPlace = 3 or zPlace = 3) and evalf(c1-c0) = evalf(2*Pi)
      then NCSpacing :=true
      fi:
   else ERROR(`Cartesian coords accepts only x, y, and z variables`,
              `cylindrical coords accepts only r, theta, and z variables`,
              `spherical coords accepts only rho, theta, and phi variables`)
   fi:
TempTable:=MakeGrid(
      a,a0,a1,b,b0,b1,c,c0,c1,CountA,CountB,CountC,
      SpaceA,SpaceB,SpaceC,Coords,ggg,
      NASpacing,NBSpacing, NCSpacing, CSArg):
CentersTable:=TempTable[1]:
DeltaA:=TempTable[2]:
DeltaB:=TempTable[3]:DeltaC:=TempTable[4]:
CornersTable:=MakeCorners(
      CentersTable,a,b,c,c0,c1,
      SpaceA,SpaceB,SpaceC,DeltaA,DeltaB,DeltaC,
      xPlace,yPlace,zPlace,ggg,CSArg):
CornersTable:=ChangeCoords(CornersTable,Coords):
PolygonsTable:=MakePoly(CornersTable,CentersTable):
polygonplot3d({(seq(seq(PolygonsTable[i,j],
                        j=1..6),
                    i=1 .. nops([entries(CentersTable)])))},
               title=PlotTitle,scaling=Scaling,
               axes=Axes,labels=Labels,op(Options))

> end:

>

Old help pages for the MAPLE V release 4 "package" asu

> `help/text/intdraw3d`:=TEXT(

> `FUNCTION: intdraw3d -creates a picture of a parameterized `,

> ` three-dimensional region that is sliced according to`,

> ` the selected parameterization`,

> ` `,

> `CALLING SEQUENCE:`,

> ` intdraw3d(boundaries);`,

> ` intdraw3d(boundaries, intdraw3d options, plot3d options);`,

> ` `,

> `PARAMETERS:`,

> `boundaries - three ranges in the form u=u0..u1,v=v0..v1,w=w0..w1`,

> ` u0 and u1 must be constants, v0 and v1 must be expressions `,

> ` depending on the variable u only, w0 and w1 are expressions`,

> ` that may depend on the variables u and v. `,

> ` The variables u,v,w may represent either rectangular, cylindrical`,

> ` or spherical coordinates, see below. `,

> ` `,

> `INTDRAW3D OPTIONS:`,

> `grid - determines the numbers of boxes to be drawn. The call grid=[m,n]`,

>

> ` is interpreted as grid=[m,n,1], i.e. drawing columns, unless the `,

> ` variable w is either theta or phi, in which case it is interpreted as `,

> ` grid=[m,n,5]. The default is grid=[5,5,5]. `,

> ` `,

> `spaces - determines the spacing between the boxes (columns). `,

> ` Its argument is a list of three constants between 0 and 1. `,

> ` Default is spaces=[.3,.3,0]. The option spaces = slices is `,

> ` equivalent to spaces= [.8,.1,0]. `,

> ` `,

> `SYNOPSIS:`,

> `-these defaults were choosen because they produced the nicest and `,

> ` most informative pictures, they may be changed easily in the plot3d `,

> ` window or in the intdraw3d calling statement (see examples bellow)`,

> `-for Cartesian coordinates, variables used must be x, y, and z.`,

> `-for cylindrical coordinates, variables used must be r, theta and z.`,

> `-for spherical coordinates, variables used must be rho, phi, and theta.`,

> `-scaling (plot3d option) default in intdraw3d is CONSTRAINED`,

> `-labels (plot3d option) default in intdraw3d is [x,y,z]`,

> `-axes (plot3d option) default in intdraw3d is FRAME`,

> ` `,

> `#if you run all examples at once your computer will be working for a while`,

> `#Cartesian coordinates examples (wedges)`,

> `intdraw3d(x= -1..0 , y= -sqrt(1-x^2) ..sqrt(1-x^2) , z= -x..0 );`,

> `intdraw3d(z= -1..0 , y= -sqrt(1-z^2) ..sqrt(1-z^2) , x= -z..0 , spaces=[.5,0],`,

> `axes=BOX);`,

> `intdraw3d(z= -1..1,x= -sqrt(1-z^2)..sqrt(1-z^2),y= -sqrt(1-x^2-z^2)..`,

> `sqrt(1-z^2-x^2),spaces=[0,0], style = PATCHNOGRID);`,

> ` `,

> `#cylindrical coordinates examples`,

> `intdraw3d(theta = 0..2*Pi, r=0..1/sqrt(2), z=r..sqrt(1-r^2), spaces=[0,0]);`,

> `intdraw3d(z= 0..Pi, r=0..Pi, theta=0..Pi, grid=[3,2], orientation = [-70,60], `,

> `spaces=[.4,.4]);`,

> `intdraw3d(r=3..5, theta = 0..2*Pi, z=-sqrt(1-(4-r)^2).. sqrt(1-(4-r)^2), `,

> `grid = [5,10 ], axes = BOX );`,

> ` `,

> `I1:= intdraw3d(r=0..3, theta = 0..3*Pi/4, z=-1..1, spaces=[.4,.4],grid = 4,`,

> `color = green, style = patchnogrid, light = [0,50 , .9,.9,.9], `,

> `light = [50,0,.6,.6,.6], ambientlight= [.4,.4,.4]):`,

> `I2:= intdraw3d(r=3..3*sqrt(2), z=-1..1, theta = arccos(3/r)..arcsin(3/r),`,

> `grid=[5,3,8],spaces=[.4,.4],color = blue, style = patchnogrid):`,

> `I3:= intdraw3d(r=3..3*sqrt(2),z=-1..1, theta=Pi-arcsin(3/r)..3/4*Pi,`,

> `spaces=[.4,.4],grid = [5,3,6], color = red, style = patchnogrid):`,

> `display3d({I1,I2,I3}, orientation = [-70,50]);`,

> ` `,

> `#spherical coordinates examples`,

> `intdraw3d(rho=0..1, phi = 0..Pi, theta=0..2*Pi,grid=[1],spaces=[.5,.5]);`,

> `intdraw3d(rho=0..1,theta= -Pi/2..Pi/2,phi=arctan(1/cos(theta))..Pi/2,`,

> `grid=[4,6,12]);`,

> `# changing the spacing`,

> `intdraw3d(rho=0..1,theta= -Pi/2..Pi/2,phi=arctan(1/cos(theta))..Pi/2,`,

> `grid=[4,6,12],spaces=[.9,0]);`,

> `intdraw3d(rho=0..1, theta = 0..2*Pi, phi=0..Pi,grid=[1,8,8], spaces=[0,0]);`,

> ` `,

> `#movie examples (may take a while)`,

> `with(plots):display([seq(intdraw3d(x=-1..1,y=-sqrt(1-x^2)..sqrt(1-x^2), `,

> `z=-sqrt(1-x^2-y^2)..sqrt(1-x^2-y^2),grid=[k,k],spaces =[0,0] ),k=1..10)],`,

> `insequence = true);`,

> ` `,

> `#this one takes 45 seconds on a Pentium60`,

> `int1:=display3d([seq(intdraw3d(x=0..1 , y=-sqrt(1-x^2) ..sqrt(1-x^2) , `,

> `z=0..x , grid = [5,5], spaces=[k/10,0]), k=0..9)], insequence = true):`,

> `int2:=display3d([seq(intdraw3d(x=0..1 , y=-sqrt(1-x^2) ..sqrt(1-x^2) ,`,

> ` z=0..x , grid = [5,5], spaces=[1,k/10]), k=0..9)], insequence = true):`,

> `int3:=display3d([seq(intdraw3d(x=0..1 , y=-sqrt(1-x^2) ..sqrt(1-x^2) ,`,

> ` z=0..x , grid = [5,5], spaces=[1-k/10,1]), k=0..9)], insequence = true):`,

> `int4:=display3d([seq(intdraw3d(x=0..1 , y=-sqrt(1-x^2) ..sqrt(1-x^2) , `,

> `z=0..x , grid = [5,5], spaces=[0,1-k/10]), k=0.10)], insequence = true):`,

> `display3d([int1, int2, int3, int4],insequence = true, orientation = [-134,62]);` ):

> `help/text/intdraw`:=TEXT(

> `FUNCTION: asu[intdraw] -creates a picture of the region of a`,

> ` two-dimentional integral`,

> ` `,

> `CALLING SEQUENCE:`,

> ` intdraw(boundaries);`,

> ` intdraw(boundaries, intdraw options, plot options);`,

> ` `,

> `PARAMETERS:`,

> `boundaries -two 2D integral boundaries written as `,

> ` var = Lower Limit..Upper Limit.`,

> ` Outer integral must be the first `,

> ` argument, then inner integral.`,

> ` `,

> `INTDRAW OPTIONS:`,

> `grid -used to vary number of columns; for polar `,

> ` coordinates, grid's second argument is used `,

> ` to vary curvature of columns as well. Numbers `,

> ` entered into grid may not always correlated to number `,

> ` of columns shown on graph. Default is grid=[15,15].`,

> ` `,

> `spaces -used to vary spacing between columns or boxes. `,

> ` Spaces arguments must range from 0 to 1. Default is `,

> ` spaces=[.3,0]. `,

> `SYNOPSIS:`,

> `-for Cartesian coordinates, variables used must be x, and y.`,

> `-for polar coordinates, variables used must be r, and theta.`,

> `-scaling (plot option) default in intdraw is CONSTRAINED`,

> `-labels (plot option) default in intdraw is [x,y]`,

> `-these defaults were choosen because they produced the nicest and `,

> ` most informative pictures, they may be changed easily in the plot `,

> ` window or in the intdraw calling statement (see examples bellow)`,

> ` `,

> `#examples`,

> `intdraw(x=-1..1, y=-sqrt(1-x^2)..sqrt(1-x^2),spaces =[.5], color = blue, `,

> `axes = BOX);`,

> `intdraw(x=-1..1, y=-sqrt(1-x^2)..sqrt(1-x^2),spaces =[.3,.3],color = yellow);`,

> `intdraw(r= 0..1, theta=0..2*Pi,grid = [3,25]); `,

> ` `,

> `region:=plot([3,t,t=0..3]),plot([t,3,t=-3..3]),plot([t,-t,t=-3..-sqrt(2)]),`,

> `plot([2*cos(t),2*sin(t),t=0..3*Pi/4]):`,

> ` `,

> `region:=plot([3,t,t=0..3]),plot([t,3,t=-3..3]),plot([t,-t,t=-3..-sqrt(2)]),`,

> `plot([2*cos(t),2*sin(t),t=0..3*Pi/4]):`,

> ` `,

> `int1:=intdraw(theta= 0..Pi/4,r=2..3/cos(theta), color = blue):`,

> `int2:=intdraw(theta= Pi/4..3*Pi/4,r=2..3/sin(theta), color = green):`,

> `display([int1,int2,region]);`,

> ` `,

> `inta:=intdraw(r=2..3, theta = 0..3/4*Pi, grid = 4,color = green):`,

> `intb:=intdraw(r=3..3*sqrt(2), theta = arccos(3/r)..arcsin(3/r),`,

> `grid = [10,25], color = blue):`,

> `intc:=intdraw(r=3..3*sqrt(2), theta = Pi-arcsin(3/r)..3/4*Pi, grid = 6,`,

> `color = red): `,

> `display([inta,intb,intc,region]);`):

>

Examples in two dimensions

Minimal syntax

> intdraw(x=0..1,y=x^2..x);

Overlaying the outline of the region, more slices, using color

> region0:=plot([x,x^2],x=0..1,thickness=3):
display(region0,
intdraw(x=0..1,y=x^2..x,grid=23,color=blue));

Matt's orignal favorite region, leading to common mistake, featured in the article
" 3D-graphics for iterated integrals" , MAPLE-Tech, 4 no.1, (1997) pp. 92-98.
http://math.asu.edu/~kawski/preprints/97mapletech.pdf

>    region:=display(
        pointplot([[2,0],[3,0],[3,3],[-3,3],[-sqrt(2),sqrt(2)]],
                  style=line),
        plot([t,3,t=-3..3]),
        plot([t,-t,t=-3..-sqrt(2)]),
        plot([2*cos(t),2*sin(t),t=0..3*Pi/4]),
        thickness=3,scaling=constrained):
region;

> int1:=intdraw(theta= 0..Pi/4,r=2..3/cos(theta), color = blue):
int2:=intdraw(theta= Pi/4..3*Pi/4,r=2..3/sin(theta), color = green):
display([int1,int2,region]);

>    inta:=intdraw(r=2..3, theta = 0..3/4*Pi,
              grid = 4,color = green):
intb:=intdraw(r=3..3*sqrt(2), theta = arccos(3/r)..arcsin(3/r),
              grid = [10,25], color = blue):
intc:=intdraw(r=3..3*sqrt(2), theta = arcsin(3/r)..3/4*Pi,
              grid = 6, color = red):
display([inta,intb,intc,region]);

>    inta:=intdraw(r=2..3, theta = 0..3/4*Pi, grid = 4,color = green):
intb:=intdraw(r=3..3*sqrt(2), theta = arccos(3/r)..arcsin(3/r),
              grid = [10,25], color = blue):
intc:=intdraw(r=3..3*sqrt(2), theta = Pi-arcsin(3/r)..3/4*Pi,
              grid = 6, color = red):
display([inta,intb,intc,region]);

>

Examples in three dimensions

It is not always triangles in this triangular pyramid

>    pyra:=pointplot3d([
   [0,0,0],[0,1,0],[1,1,0],[1,0,0],[0,0,0],[0,0,1],
   [0,1,1],[1,1,1],[1,0,1],[0,0,1],[0,0,0],[0,1,0],
   [0,1,1],[1,1,1],[1,1,0],[1,0,0],[1,0,1]],
       style=line,color=black,thickness=1),
   pointplot3d([
   [0,1,1],[1,0,1],[1,1,1],[0,1,1],[1,0,0],[1,1,1],
   [1,0,1],[1,0,0]],
       style=line,color=black,thickness=3):
display([pyra],orientation=[-125,70]);

>    display([pyra,
        intdraw3d(x=0..1,y=1-x..1,z=y..1,
          grid=[7,31,31],spaces=[0.95,0,0],
          shading=ZHUE,orientation=[-125,70],
          axes=frame,tickmarks=[2,2,2])]);

>    display([pyra,
        intdraw3d(x=0..1,y=1-x..1,z=y..1,
          grid=[31,7,31],spaces=[0,0.95,0],
          shading=ZHUE,orientation=[-125,70],
          axes=frame,tickmarks=[2,2,2])]);

>    intdraw3d(rho=0..2,theta=0..2*Pi,phi=0..Pi/6,
          grid=[12,12,12],spaces=[0.1,0.3,0.7],
          shading=ZHUE,orientation=[80,75]);

> intdraw3d(x=-1..0,y=-sqrt(1-x^2)..sqrt(1-x^2),z=-x..0,
grid=[10,20]);

>    intdraw3d(z=-1..1,x=-sqrt(1-z^2)..sqrt(1-z^2),
          y=-sqrt(1-x^2-z^2)..sqrt(1-z^2-x^2),
          spaces=[0.1,0.1], style = PATCHNOGRID,grid=[20,20]);

cylindrical coordinates examples

>    intdraw3d(theta=0..2*Pi,r=0..1,z=sqrt(3)*r..sqrt(4-r^2),
          grid=[24,9,13],spaces=[0,0,0.95],
          axes=none,orientation=[-50,70]);

> intdraw3d(r=3..5,theta=0..2*Pi,z=-sqrt(1-(4-r)^2).. sqrt(1-(4-r)^2),
grid = [7,12 ], axes = none,orientation=[-140,60] );

spherical coordinates examples

> intdraw3d(rho=0..1,phi=0..Pi,theta=0..2*Pi,
grid=[1],spaces=[.5,.5]);

> intdraw3d(rho=0..1,theta=0..2*Pi,phi=0..Pi,
grid=[1,13,13],spaces=[0,0.4,0.4]);

movie examples (may take a while)

>    display([seq(
    intdraw3d(x=-1..1,y=-sqrt(1-x^2)..sqrt(1-x^2),
              z=-sqrt(1-x^2-y^2)..sqrt(1-x^2-y^2),
              grid=[k,k,k],spaces =[0,0],axes=none),
       k=1..13)],insequence = true,
       title=`This is an animation. Play it.`,
       titlefont=[TIMES,BOLD,16]);

>    display3d([
     seq(intdraw3d(rho=0..2,theta=0..2*Pi,phi=0..Pi/6,
              grid=[9,24,5],spaces=[k/10,0],axes=NONE,title=` `),k=0..9),
     seq(intdraw3d(rho=0..2,theta=0..2*Pi,phi=0..Pi/6,
              grid=[9,24,5],spaces=[0,0,1-k/10],axes=NONE,title=` `),k=0..9)
     ],insequence = true,orientation=[60,75],shading=ZHUE,
       title=`This is an animation. Play it.`,
       titlefont=[TIMES,BOLD,16]);

1994: "this one takes 45 seconds on a Pentium60"

>    display([
    seq(intdraw3d(x=0..1,y=-sqrt(1-x^2) ..sqrt(1-x^2),z=0..x,
           grid=[11,11],spaces=[k/10,0],axes=none),k=0..9),
    seq(intdraw3d(x=0..1,y=-sqrt(1-x^2)..sqrt(1-x^2),z=0..x,
           grid=[11,11],spaces=[1,k/10],axes=none),k=0..9),
    seq(intdraw3d(x=0..1,y=-sqrt(1-x^2)..sqrt(1-x^2),z=0..x,
           grid = [11,11], spaces=[1-k/10,1],axes=none),k=0..9),
    seq(intdraw3d(x=0..1,y=-sqrt(1-x^2)..sqrt(1-x^2),z=0..x,
           grid = [11,11], spaces=[0,1-k/10],axes=none),k=0..10)
    ],insequence = true,orientation = [-134,62],
       title=`This is an animation. Play it.`,
       titlefont=[TIMES,BOLD,16]);

>

start working here (after executing the second to top section)

> intdraw(x=2..3,y=1/x..sqrt(x));

>