> restart; with(plots):

Checking set-up of iterated integrals

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.

 >

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);

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]

 > 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]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]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]

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(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));

 >

 >

 >

 >