#calcplot package version 1.1 by Tim Murdoch. Originally programmed:
#July/August, 1994. Revised: June, 1995
#Modified by Richard Pulskamp for Release 5 November 1998
#The package consists of thirty procedures. Two procedures are for
#plotting surfaces of rotation, five procedures are for plotting
#non-rectangular domains in the plane (e.g., domains for double
#integrals in cartesian and polar coordinates), four procedures are for
#plotting graphs of functions over non-rectangular domains in the
#plane, and nineteen procedures are for plotting non-rectangluar
#regions (e.g., domains for triple integrals in cartesain, cylindrical,
#and spherical coordinates) in space.
#
#Revisions for v. 1.1: Fixed plot options to include defaults unless
#specifically over-ridden; added regionplot2d and regionplot3d
#routines and help files.

#This file is free. You may do anything you like with it with one exception: 
#if you change anything in this file, then you MUST rename it!

#Last modified: 6/16/95

#This file, can be saved as calcplot.m and read into an active Maple session 
#using the "read" command. See the calcplot manual for details.

#Define a procedure that takes either an expression, a procedural function of a
#single variable, or a parametric curve in the form of a two-element list, 
#along with a range of the form x=a..b and an axis of rotation in the form x=c
#as input.


approxint2d := proc(f,xexpress, yexpress ,n ,m)
           local a,b,c,d, dx,i,j, dy;
        a := op(1, op(2, xexpress));
        b := op(2, op(2, xexpress));
        c := op(1, op(2, yexpress));
        d := op(2, op(2, yexpress));
                         dx := evalf((b-a)/n);  dy := evalf((d-c)/m);
          sum(sum(subs(x=a+i*dx,y=c+j*dy,f),i=1..n),j=1..m)*dx*dy;
end: 

blockapp := proc( f, xexpress, yexpress,m,n)
local x1, x2, y1,y2, i ,j , a, b ,c, d, blocks:
a := op(1, op(2, xexpress));
 b := op(2, op(2, xexpress));
 c := op(1, op(2, yexpress));
  d := op(2, op(2, yexpress));
blocks := {} :
for i from 1 to m do
x1 := a + (i-1)*(b-a)/m:
x2 := a + i*(b-a)/m:
for j from 1 to n do
y1 := c + (j-1)*(d-c)/n:
y2 := c + j*(d-c)/n:
blocks := blocks union {[x1, s*y1+ (1-s)*y2, t*subs(x=x2,y=y2, f)], [x2,
s*y1+(1-s)*y2, t*subs(x=x2, y=y2, f)], [s*x1+ (1-s)*x2, y1,
t*subs(x=x2,y=y2, f)], [ s*x1+ (1-s)*x2, y2, t*subs(x=x2,y=y2, f)],
[s*x1+ (1-s)*x2, t*y1+(1-t)*y2, subs(x=x2,y=y2, f)]} :
od:
od:
plot3d( blocks, s=0..1, t=0..1, grid=[2,2]) ;  end:


rotyplot:=proc()

local a, b, c, d, f, i, xbound, rotaxis, f1, opts, opt_seq, p1, profile, tloc;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

f:=args[1];
xbound:=args[2];
rotaxis:=lhs(args[3])=evalf(rhs(args[3]));
a:=op(1, rhs(xbound));
b:=op(2, rhs(xbound));

### WARNING: semantics of type `string` have changed
if not type (rhs(xbound), range) then
ERROR(`range expression for x is incorrect`) fi;
### WARNING: semantics of type `string` have changed
if not type (rhs(rotaxis), numeric) then
ERROR(`axis expression is not correct`) fi;
if not (lhs(xbound) = 'x') or not (lhs(rotaxis) = 'x') then
ERROR(`the variable name must be x`) fi;
if not type (evalf(a), numeric) or not type(evalf(b), numeric) then 
ERROR(`range limits are not real numbers` ) fi;

c:=rhs(rotaxis);
d:=lhs(xbound);

opts:=style=PATCH, scaling=CONSTRAINED, axes=BOXED, grid=[25,25],labels=[` `,'x','y'];

if not type (f, procedure) or  not type(f, list) then
f1:=unapply(f,d); fi;

if type(f,procedure) then
f1:=f; fi;

if not type(f, list) then
profile:=[(d-c)*sin(tloc),(d-c)*cos(tloc)+c,f1(d)]; fi;

if type(f,list) then f1:=unapply(f,d);
profile:=[(f1(d)[1]-c)*sin(tloc),(f1(d)[1]-c)*cos(tloc)+c,f1(d)[2]]; fi;

if nargs = 3 then
p1:=plot3d(profile,d=a..b,tloc=0..2*Pi,opts); fi;

if nargs > 3 then
opt_seq:=opts,seq(args[i],`i`=4..nargs);
p1:=plot3d(profile,d=a..b,tloc=0..2*Pi),opt_seq; fi;

plots[display3d](p1);

end:

#A y-axis rotation plotter. Same input as xrotplot, except that the axis of 
#rotation is given in the form y=c.

rotxplot:=proc()

local a, b, c, d, f1, f, i, xbound, rotaxis, opts, opt_seq, p1, profile, tloc;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

f:=args[1];
xbound:=args[2];
rotaxis:=lhs(args[3])=evalf(rhs(args[3]));
a:=op(1, rhs(xbound));
b:=op(2, rhs(xbound));


### WARNING: semantics of type `string` have changed
if not type (rhs(xbound), range) then
ERROR(`range expression for x is incorrect`) fi;
## WARNING: semantics of type `string` have changed
if not type (rhs(rotaxis), numeric) then
ERROR(`axis expression is not correct`) fi;
if not (lhs(xbound) = 'x') or not (lhs(rotaxis) = 'y') then
ERROR(`the variable name must be x`) fi;
if not type (evalf(a), numeric) or not type(evalf(b), numeric) then 
ERROR(`range limits are not real numbers` ) fi;

c:=rhs(rotaxis);
d:=lhs(xbound);

opts:=style=PATCH, scaling=CONSTRAINED, axes=BOXED, grid=[25,25],labels=[` `,'x','y'];

if not type (f, procedure) or not type (f, list) then
f1:=unapply(f,d); fi;

if type (f, procedure) then
f1:=f; fi;

if not type(f,list) then
profile:=[(f1(d)-c)*sin(tloc),d,(f1(d)-c)*cos(tloc)+c]; fi;

if type(f,list) then f1:=unapply(f,d);
profile:=[(f1(d)[2]-c)*sin(tloc),f1(d)[1],(f1(d)[2]-c)*cos(tloc)+c]; fi;

if nargs = 3 then
p1:=plot3d(profile,d=a..b,tloc=0..2*Pi,opts); fi;

if nargs > 3 then
opt_seq:=opts,seq(args[i],`i`=4..nargs);
p1:=plot3d(profile,d=a..b,tloc=0..2*Pi),opt_seq; fi;

plots[display3d](p1);

end:

#Procedures for double integration region plotters that return a plot 
#structure. These procedures are analogous to the drdtplot and dtdrplot 
#procedures in mvcal. The procedures take two variable ranges as arguments.

#Type I region.

dydxplot:=proc()

local a, b, c, d, g, h, i, j, f1, f2, opts, opt_seq, p1, pset, v, ybound, xbound;

if nargs < 2 then
ERROR(`there must be at least two arguments`) fi;

ybound:=args[1];
xbound:=args[2];
a:=op(1, rhs(xbound));
b:=op(2, rhs(xbound));

### WARNING: semantics of type `string` have changed
if not type (rhs(xbound), range) then
ERROR(`range expression for x is incorrect`) fi;
### WARNING: semantics of type `string` have changed
if not type (rhs(ybound), range) then
ERROR(`range expression for y is not correct`) fi;
if not (lhs(xbound) = 'x') then
ERROR(`the second variable name must be x`) fi;
if not (lhs(ybound) = 'y')  then
ERROR(`the first variable name must be y`) fi;
if not type (evalf(a), numeric) or not type(evalf(b), numeric) then 
ERROR(`range limits for x are not real numbers` ) fi;

c:=op(1, rhs(ybound));
d:=op(2, rhs(ybound));
g:=lhs(xbound);
h:=lhs(ybound);

opts:=scaling=CONSTRAINED, axes=NORMAL;

f1:=unapply(c,g);
f2:=unapply(d,g);
p1:=plot({f1,f2},a..b);
pset:={[v,f1(v),v=a..b],[v,f2(v),v=a..b]};

for j from 1 to 5 do
pset:=pset union
{[a+(j-1)/4*(b-a),v*f2(a+(j-1)/4*(b-a))+(1-v)*f1(a+(j-1)/4*(b-a)),v=0..1]}; od;

if nargs = 2 then
p1:=plot(pset,opts); fi;

if nargs > 2 then
opt_seq:=opts,seq(args[i],`i`=3..nargs);
p1:=plot(pset),opt_seq; fi;

plots[display](p1);

end:

#Type II region

dxdyplot:=proc()

local a, b, c, d, g, h, i, j, f1, f2, opts, opt_seq, p1, pset, v, ybound, xbound;

if nargs < 2 then
ERROR(`there must be at least two arguments`) fi;

xbound:=args[1];
ybound:=args[2];
a:=op(1, rhs(ybound));
b:=op(2, rhs(ybound));

### WARNING: semantics of type `string` have changed
if not type (rhs(xbound), range) then
ERROR(`range expression for x is incorrect`) fi;
### WARNING: semantics of type `string` have changed
if not type (rhs(ybound), range) then
ERROR(`range expression for y is not correct`) fi;
if not (lhs(xbound) = 'x') then
ERROR(`the second variable name must be x`) fi;
if not (lhs(ybound) = 'y')  then
ERROR(`the first variable name must be y`) fi;
if not type(evalf(a), numeric) or not type(evalf(b), numeric) then 
ERROR(`range limits for x are not real numbers` ) fi;

c:=op(1, rhs(xbound));
d:=op(2, rhs(xbound));
g:=lhs(xbound);
h:=lhs(ybound);

opts:=scaling=CONSTRAINED, axes=NORMAL;

f1:=unapply(c,h); 
f2:=unapply(d,h); 
pset:={[f1(h),h,h=a..b],[f2(h),h,h=a..b]};  

for j from 1 to 5 do
pset:=pset union
{[v*f2(a+(j-1)/4*(b-a))+(1-v)*f1(a+(j-1)/4*(b-a)),
a+(j-1)/4*(b-a),v=0..1]}; od;

if nargs = 2 then
p1:=plot(pset,opts); fi;

if nargs >2 then
opt_seq:=opts,seq(args[i],`i`=3..nargs);
p1:=plot(pset),opt_seq; fi;

plots[display](p1);

end:

#Procedure to plot regions in polar coordinates. These procedures are
#entirely similar to the cartesian coordinate plotters above.

#Radius is a function of angle.

drdtplot:=proc()

local a, b, c, d, dr, dt, g, g1, h, h1, i, opts, opt_seq, p1, pset, rbound,
tbound, tloc;

if nargs < 2 then
ERROR(`there must be at least two arguments`) fi;

rbound:=args[1];
tbound:=args[2];
a:=op(1, rhs(tbound));
b:=op(2, rhs(tbound));

### WARNING: semantics of type `string` have changed
if not type (rhs(rbound), range) then
ERROR(`range expression for r is incorrect`) fi;
### WARNING: semantics of type `string` have changed
if not type (rhs(tbound), range) then
ERROR(`range expression for theta is not correct`) fi;
if not (lhs(rbound) = 'r') then
ERROR(`the first variable name must be r`) fi;
if not (lhs(tbound) = 'theta')  then
ERROR(`the first variable name must be theta`) fi;
if not type(evalf(a), numeric) or not type(evalf(b), numeric) then 
ERROR(`range limits for theta are not real numbers` ) fi;

c:=lhs(tbound);
d:=lhs(rbound);
g:=op(1, rhs(rbound));
h:=op(2, rhs(rbound));

g1:=unapply(g,c);
h1:=unapply(h,c);

opts:=scaling=CONSTRAINED, axes=NORMAL;

pset:={[g1(tloc),tloc,tloc=a..b],[h1(tloc),tloc,tloc=a..b]};

for i from 0 to 5 do
dt:= (5-i)/5*a + i/5*b ;
dr:=evalf(h1(dt) - g1(dt));

if (dr < - 0.000001) then 
ERROR(`the the outer radius is less than the inner radius within theta range`)  fi;

if (dr > 10^(-3) ) then 
pset:=pset union 
{[tloc*h1(dt)+(1-tloc)*g1(dt),dt, tloc=0..1]} fi;
od;

if nargs =2 then
p1:=plot(pset,coords=polar, opts); fi;

if nargs > 2 then
opt_seq:=opts,seq(args[i], `i`=3..nargs);
p1:=plot(pset, coords=polar), opt_seq; fi;

plots[display](p1);

end:

#Angle as a function of radius.

dtdrplot:=proc()

local a, b, c, d, dr, dt, g, g1, h, h1, i, opts, opt_seq, p1, pset, rbound,
tbound, tloc;

if nargs < 2 then
ERROR(`there must be at least two arguments`) fi;

tbound:=args[1];
rbound:=args[2];
a:=op(1, rhs(rbound));
b:=op(2, rhs(rbound));

### WARNING: semantics of type `string` have changed
if not type (rhs(rbound), range) then
ERROR(`range expression for r is incorrect`) fi;
### WARNING: semantics of type `string` have changed
if not type (rhs(tbound), range) then
ERROR(`range expression for theta is not correct`) fi;
if not (lhs(rbound) = 'r') then
ERROR(`the first variable name must be r`) fi;
if not (lhs(tbound) = 'theta')  then
ERROR(`the first variable name must be theta`) fi;
if not type(evalf(a), numeric) or not type(evalf(b), numeric) then 
ERROR(`range limits for r are not real numbers` ) fi;

c:=lhs(rbound);
d:=lhs(tbound);
g:=op(1, rhs(tbound));
h:=op(2, rhs(tbound));

g1:=unapply(g,c);
h1:=unapply(h,c);

opts:=scaling=CONSTRAINED, axes=NORMAL;

pset:={[tloc,g1(tloc),tloc=a..b],[tloc,h1(tloc),tloc=a..b]};

for i from 0 to 5 do
dr:= (5-i)/5*a + i/5*b ;
dt:=evalf(h1(dr) - g1(dr));

if (evalf(dr) < - 0.000001) then 
ERROR(`outer radius is less than inner radius within theta range`)  fi;

if (evalf(dr) > 10^(-3) ) then 
pset:=pset union 
{[dr,tloc*h1(dr)+(1-tloc)*g1(dr), tloc=0..1]} fi;
od;

if nargs = 2 then
p1:=plot(pset,coords=polar, opts); fi;

if nargs > 2 then
opt_seq:=opts,seq(args[i],`i`=3..nargs); 
plot(pset,coords=polar),opt_seq; fi;

plots[display](p1);

end:

#regionplot2d plots the image of a region of the type h(v)<=u<=k(v),a<=v<=b
#under a coordinate transform [x(u,v),y(u,v)].

regionplot2d:=proc()
local a, b, h, h1, k, k1, i, j, face1, face2, face3, face4, leftedge, opts, p1, pset, rightedge, t, ubound, u1, v1, vbound, xtrans, ytrans;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

ubound:=args[1];
vbound:=args[2];
a:=op(1, rhs(vbound));
b:=op(2, rhs(vbound));

### WARNING: semantics of type `string` have changed
if not type (rhs(ubound), range) then
ERROR(`range expression for first argument is not correct`) fi;
### WARNING: semantics of type `string` have changed
if not type (rhs(ubound), range) then
ERROR(`range expression for second argument is not correct`) fi;
if not type(args[3],list) then
ERROR(`the third argument must be a list of transformation functions`) fi;
if not type (evalf(a), numeric) or not type(evalf(b), numeric) then 
ERROR(`range limits for the second argument are not real numbers` ) fi;

h:=op(1, rhs(ubound));
k:=op(2, rhs(ubound));
u1:=lhs(ubound);
v1:=lhs(vbound);

xtrans:=unapply(args[3][1],u1,v1);
ytrans:=unapply(args[3][2],u1,v1);

h1:=unapply(h,v1);
k1:=unapply(k,v1);
leftedge:=t*k1(a)+(1-t)*h1(a);
rightedge:=t*k1(b)+(1-t)*h1(b);
face1:=[xtrans(leftedge,a),ytrans(leftedge,a),t=0..1];
face2:=[xtrans(rightedge,b),ytrans(rightedge,b),t=0..1];
face3:=[xtrans(h1(v1),v1),ytrans(h1(v1),v1),v1=a..b];
face4:=[xtrans(k1(v1),v1),ytrans(k1(v1),v1),v1=a..b];
pset:={face1,face2,face3,face4};

for j from 1 to 6 do
pset:=pset union
{[xtrans(t*k1(a+(j-1)/5*(b-a))+(1-t)*h1(a+(j-1)/5*(b-a)),a+(j-1)/5*(b-a)), ytrans(t*k1(a+(j-1)/5*(b-a))+(1-t)*h1(a+(j-1)/5*(b-a)),a+(j-1)/5*(b-a)), t=0..1]}; od;

if nargs = 3 then
opts:=axes=NORMAL,scaling=CONSTRAINED;
else
opts:=axes=NORMAL,scaling=CONSTRAINED,seq(args[i],`i`=4..nargs);
fi;

plot(pset,opts);

end:

#Procedures to plot graphs of functions of the form z=h(x,y) over
#domains of type I : f(x)<=y<=g(x), a<=x<=b.  II: r(y)<=x<=s(y),
#c<=y<=d Special case h=0 yields plots of regions of the above type in
#the xy-plane. Note that these types of plots can be obtained directly
#with the 3d plot command; however, orientation of the plots is a bit
#of a problem.

#This one handles type I regions.

yxgraphplot:=proc()

local a, b, c, d, g, h, i, k, f1, g1, h1, opts, opt_seq, p1, tloc, sloc,
surface, xbound, ybound;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

h:=args[1];
ybound:=args[2];
xbound:=args[3];
a:=op(1,rhs(xbound));
b:=op(2,rhs(xbound));

### WARNING: semantics of type `string` have changed
if not type (rhs(xbound), range) then
ERROR(`range expression for x is incorrect`) fi;
### WARNING: semantics of type `string` have changed
if not type (rhs(ybound), range) then
ERROR(`range expression for y is not correct`) fi;
if not (lhs(xbound) = 'x') then
ERROR(`the second variable name must be x`) fi;
if not (lhs(ybound) = 'y')  then
ERROR(`the first variable name must be y`) fi;
if not type (evalf(a), numeric) or not type(evalf(b), numeric) then 
ERROR(`range limits for x are not real numbers` ) fi;

c:=op(1,rhs(ybound));
d:=op(2,rhs(ybound));
g:=lhs(xbound);
k:=lhs(ybound);

tloc:=(1-u)*a+u*b;
f1:=unapply(c,g);
g1:=unapply(d,g);
sloc:=(1-v)*f1(tloc)+v*g1(tloc);

opts:=style=PATCH,axes=BOXED,scaling=CONSTRAINED, grid=[20,20], labels=['x','y','z'];

if not type (h, procedure) then 
h1:=unapply(h,(g,k)); fi;

if type(h,procedure) then h1:=h fi;

surface:=[tloc,sloc,h1(tloc,sloc)];

if nargs = 3 then
p1:=plot3d(surface,u=0..1,v=0..1,opts); fi;

if nargs > 3 then
opt_seq:=opts,seq(args[i],`i`=4..nargs);
p1:=plot3d(surface, u=0..1,v=0..1,opt_seq); fi;

plots[display3d](p1);

end:

#This one handles type II regions.

xygraphplot:=proc()

local a, b, c, d, g, h, i, k, f1, g1, h1, opts, opt_seq, p1, sloc, tloc,
surface, xbound, ybound;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

h:=args[1];
xbound:=args[2];
ybound:=args[3];
a:=op(1,rhs(ybound));
b:=op(2,rhs(ybound));

### WARNING: semantics of type `string` have changed
if not type (rhs(xbound), range) then
ERROR(`range expression for x is incorrect`) fi;
### WARNING: semantics of type `string` have changed
if not type (rhs(ybound), range) then
ERROR(`range expression for y is not correct`) fi;
if not (lhs(xbound) = 'x') then
ERROR(`the second variable name must be x`) fi;
if not (lhs(ybound) = 'y')  then
ERROR(`the first variable name must be y`) fi;
if not type (evalf(a), numeric) or not type(evalf(b), numeric) then 
ERROR(`range limits for x are not real numbers` ) fi;

c:=op(1,rhs(xbound));
d:=op(2,rhs(xbound));
g:=lhs(xbound);
k:=lhs(ybound);

tloc:=(1-u)*a+u*b;
f1:=unapply(c,k); 
g1:=unapply(d,k);
sloc:=(1-v)*f1(tloc)+v*g1(tloc);

opts:=style=PATCH,axes=BOXED,scaling=CONSTRAINED, grid=[20,20], labels=['x','y','z'];

if not type (h, procedure) then 
h1:=unapply(h,(g,k)); fi;

if type(h,procedure) then h1:=h fi;

surface:=[sloc,tloc,h1(sloc,tloc)];

if nargs = 3 then
p1:=plot3d(surface,u=0..1,v=0..1,opts); fi;

if nargs > 3 then
opt_seq:=opts,seq(args[i],`i`=4..nargs);
p1:=plot3d(surface,u=0..1,v=0..1,opt_seq); fi;

plots[display3d](p1);

end:

#Graphs over polar regions. Same concept as above, except the domain
#is expressed in polar coordinates. The first argument (a function)
#should be expressed in polar coordinates.

rtgraphplot:=proc()

local a, b, c, d, dr, dt, g, h, i, k, f1, g1, h1, opts, opt_seq, p1, rbound,
tbound, tloc, sloc, surface;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

h:=args[1];
rbound:=args[2];
tbound:=args[3];
a:=op(1,rhs(tbound));
b:=op(2,rhs(tbound));

### WARNING: semantics of type `string` have changed
if not type (rhs(tbound), range) then
ERROR(`range expression for theta is incorrect`) fi;
### WARNING: semantics of type `string` have changed
if not type (rhs(rbound), range) then
ERROR(`range expression for r is not correct`) fi;
if not (lhs(tbound) = 'theta') then
ERROR(`the second variable name must be theta`) fi;
if not (lhs(rbound) = 'r')  then
ERROR(`the first variable name must be r`) fi;
if not type (evalf(a), numeric) or not type(evalf(b), numeric) then 
ERROR(`range limits for theta are not real numbers` ) fi;

c:=op(1,rhs(rbound));
d:=op(2,rhs(rbound));
g:=lhs(rbound);
k:=lhs(tbound);

tloc:=(1-u)*a+u*b;
f1:=unapply(c,k);
g1:=unapply(d,k);
sloc:=(1-v)*f1(tloc)+v*g1(tloc);

opts:=style=PATCH,axes=BOXED,scaling=CONSTRAINED, grid=[20,20], labels=['x','y','z'];

if not type (h, procedure) then 
h1:=unapply(h,(g,k)); fi;

if type(h,procedure) then h1:=h fi;

surface:=[sloc,tloc,h1(sloc,tloc)];

for i from 0 to 10 do
dt:= (10-i)/10*a + i/10*b ;
dr:=evalf(g1(dt) - f1(dt));

if (dr < - 0.000001) then 
ERROR(`outer radius is less than inner radius within theta range`)  fi; od;

if nargs = 3 then
p1:=plots[cylinderplot](surface,u=0..1,v=0..1,opts); fi;

if nargs >3 then
opt_seq:=opts,seq(args[i],`i`=4..nargs);
p1:=plots[cylinderplot](surface, u=0..1,v=0..1,opt_seq); fi;

plots[display3d](p1);

end:

#Reverse the order: angle as a function of radius.

trgraphplot:=proc()

local a, b, c, d, g, k, f1, g1, h1, i, h, opts, opt_seq, p1, rbound, tbound,
tloc, sloc, surface;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

h:=args[1];
tbound:=args[2];
rbound:=args[3];
a:=op(1,rhs(rbound));
b:=op(2,rhs(rbound));

### WARNING: semantics of type `string` have changed
if not type (rhs(tbound), range) then
ERROR(`range expression for theta is incorrect`) fi;
### WARNING: semantics of type `string` have changed
if not type (rhs(rbound), range) then
ERROR(`range expression for r is not correct`) fi;
if not (lhs(tbound) = 'theta') then
ERROR(`the second variable name must be theta`) fi;
if not (lhs(rbound) = 'r')  then
ERROR(`the first variable name must be r`) fi;
if not type (evalf(a), numeric) or not type(evalf(b), numeric) then 
ERROR(`range limits for theta are not real numbers` ) fi;

c:=op(1,rhs(tbound));
d:=op(2,rhs(tbound));
g:=lhs(rbound);
k:=lhs(tbound);

tloc:=(1-u)*a+u*b;
f1:=unapply(c,g);
g1:=unapply(d,g);
sloc:=(1-v)*f1(tloc)+v*g1(tloc);

opts:=style=PATCH,axes=BOXED,scaling=CONSTRAINED, grid=[20,20], labels=['x','y','z'];

if not type (h, procedure) then 
h1:=unapply(h,(g,k)); fi;

if type(h,procedure) then h1:=h fi;

surface:=[tloc,sloc,h1(tloc,sloc)];

if nargs = 3 then
p1:=plots[cylinderplot](surface,u=0..1,v=0..1,opts); fi;

if nargs >3 then
opt_seq:=opts,seq(args[i],`i`=4..nargs);
p1:=plots[cylinderplot](surface, u=0..1,v=0..1,opt_seq); fi;

plots[display3d](p1);

end:

#Six triple integral region plots for domains expressed in Cartesian
#coordinate systems.

#1. dxdydz

dxdydzplot:= proc()

local surface, f, f1, g, g1, h, h1, k, k1, p1, a, b, u, v, uloc, uloc1, uloc2,
vloc, wloc, xloc, yloc, zloc, i, x1, y1, z1, opts, opt_seq, xbound, ybound, zbound;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

xbound:=args[1];
ybound:=args[2];
zbound:=args[3];
a:=op(1, rhs(zbound));
b:=op(2, rhs(zbound));

### WARNING: semantics of type `string` have changed
if not type (rhs(zbound), range) or not (op(1, xbound) = 'x') 
or not (op(1, ybound) = 'y') or not (op(1, zbound) = 'z') 
or not type (evalf(a),numeric) or not type(evalf(b), numeric) then 
ERROR(`input expression is not of dxdydz type` ) fi;

f:=op(1, rhs(xbound));
g:=op(2, rhs(xbound));
h:=op(1, rhs(ybound));
k:=op(2, rhs(ybound));
x1:=lhs(xbound);
y1:=lhs(ybound);
z1:=lhs(zbound);

f1:=unapply(f,(y1,z1));
g1:=unapply(g,(y1,z1));
h1:=unapply(h,z1);
k1:=unapply(k,z1);

uloc:=u*b+(1-u)*a;
uloc1:=u*h1(a)+(1-u)*k1(a);
uloc2:=u*h1(b)+(1-u)*k1(b);
vloc:=v*k1(uloc)+(1-v)*h1(uloc);
wloc:=v*f1(h1(uloc),uloc)+(1-v)*g1(h1(uloc),uloc);
xloc:=v*f1(k1(uloc),uloc)+(1-v)*g1(k1(uloc),uloc);
yloc:=v*f1(uloc1,a)+(1-v)*g1(uloc1,a);
zloc:=v*f1(uloc2,b)+(1-v)*g1(uloc2,b);

surface:={[f1(vloc,uloc),vloc,uloc],[g1(vloc,uloc),vloc,uloc],
[wloc,h1(uloc),uloc],[xloc,k1(uloc),uloc],[yloc,uloc1,a],[zloc,uloc2,b]};

opts:=axes=BOXED, grid =[20,20],scaling=CONSTRAINED, style=PATCH;

if nargs = 3 then
p1:=plot3d(surface, u=0..1, v=0..1, opts);
fi;

if nargs > 3 then
opt_seq:=opts,seq(args[i],`i`=4..nargs);
p1:=plot3d(surface, u=0..1, v=0..1,opt_seq); fi;

plots[display3d](p1);

end :

#2. dydxdz

dydxdzplot:= proc()

local surface, f, f1, g, g1, h, h1, k, k1, p1, a, b, u, v,uloc, uloc1,
uloc2, vloc, wloc, xloc, yloc, zloc, i, x1, y1, z1, opts, opt_seq, xbound,
ybound, zbound;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

ybound:=args[1];
xbound:=args[2];
zbound:=args[3];
a:=op(1, rhs(zbound));
b:=op(2, rhs(zbound));

### WARNING: semantics of type `string` have changed
if not type (rhs(zbound), range) or not (op(1, xbound) = 'x') 
or not (op(1, ybound) = 'y') or not (op(1, zbound) = 'z') 
or not type (evalf(a),numeric) or not type(evalf(b), numeric) then 
ERROR(`input expression is not of dydxdz type` ) fi;

f:=op(1, rhs(ybound));
g:=op(2, rhs(ybound));
h:=op(1, rhs(xbound));
k:=op(2, rhs(xbound));
x1:=lhs(xbound);
y1:=lhs(ybound);
z1:=lhs(zbound);

f1:=unapply(f,(x1,z1));
g1:=unapply(g,(x1,z1));
h1:=unapply(h,z1);
k1:=unapply(k,z1);

uloc:=u*b+(1-u)*a;
uloc1:=u*h1(a)+(1-u)*k1(a);
uloc2:=u*h1(b)+(1-u)*k1(b);
vloc:=v*k1(uloc)+(1-v)*h1(uloc);
wloc:=v*f1(h1(uloc),uloc)+(1-v)*g1(h1(uloc),uloc);
xloc:=v*f1(k1(uloc),uloc)+(1-v)*g1(k1(uloc),uloc);
yloc:=v*f1(uloc1,a)+(1-v)*g1(uloc1,a);
zloc:=v*f1(uloc2,b)+(1-v)*g1(uloc2,b);

surface:={[vloc,f1(vloc,uloc),uloc],[vloc,g1(vloc,uloc),uloc],
[h1(uloc),wloc,uloc],[k1(uloc),xloc,uloc],[uloc1,yloc,a],[uloc2,zloc,b]};

opts:=axes=BOXED, grid =[20,20], scaling=CONSTRAINED, style=PATCH;

if nargs = 3 then
p1:=plot3d(surface, u=0..1, v=0..1,opts);
fi;

if nargs > 3 then
opt_seq:=opts,seq(args[i],`i`=4..nargs);
p1:=plot3d(surface, u=0..1, v=0..1,opt_seq); fi;

plots[display3d](p1);

end :

#3. dxdzdy

dxdzdyplot:= proc()

local surface, f, f1, g, g1, h, h1, k, k1, p1, a, b, u, v, uloc,
uloc1, uloc2, vloc, wloc, xloc, yloc, zloc, i, x1, y1, z1, opts, opt_seq,
xbound, ybound, zbound;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

xbound:=args[1];
zbound:=args[2];
ybound:=args[3];
a:=op(1, rhs(ybound));
b:=op(2, rhs(ybound));

### WARNING: semantics of type `string` have changed
if not type (rhs(zbound), range) or not (op(1, xbound) = 'x') 
or not (op(1, ybound) = 'y') or not (op(1, zbound) = 'z') 
or not type (evalf(a),numeric) or not type(evalf(b), numeric) then 
ERROR(`input expression is not of dxdzdy type` ) fi;

f:=op(1, rhs(xbound));
g:=op(2, rhs(xbound));
h:=op(1, rhs(zbound));
k:=op(2, rhs(zbound));
x1:=lhs(xbound);
y1:=lhs(ybound);
z1:=lhs(zbound);

f1:=unapply(f,(z1,y1));
g1:=unapply(g,(z1,y1));
h1:=unapply(h,y1);
k1:=unapply(k,y1);

uloc:=u*b+(1-u)*a;
uloc1:=u*h1(a)+(1-u)*k1(a);
uloc2:=u*h1(b)+(1-u)*k1(b);
vloc:=v*k1(uloc)+(1-v)*h1(uloc);
wloc:=v*f1(h1(uloc),uloc)+(1-v)*g1(h1(uloc),uloc);
xloc:=v*f1(k1(uloc),uloc)+(1-v)*g1(k1(uloc),uloc);
yloc:=v*f1(uloc1,a)+(1-v)*g1(uloc1,a);
zloc:=v*f1(uloc2,b)+(1-v)*g1(uloc2,b);

surface:={[f1(vloc,uloc),uloc,vloc],[g1(vloc,uloc),uloc,vloc],
[wloc,uloc,h1(uloc)],[xloc,uloc,k1(uloc)],[yloc,a,uloc1],[zloc,b,uloc2]};

opts:=axes=BOXED, grid =[20,20], scaling=CONSTRAINED, style=PATCH;

if nargs = 3 then
p1:=plot3d(surface, u=0..1, v=0..1, opts);
fi;

if nargs > 3 then
opt_seq:=opts,seq(args[i],`i`=4..nargs);
p1:=plot3d(surface, u=0..1, v=0..1,opt_seq); fi;

plots[display3d](p1);

end :

#4. dzdxdy

dzdxdyplot:= proc()

local surface, f, f1, g, g1, h, h1, k, k1, p1, a, b, u, v, uloc,
uloc1, uloc2, vloc, wloc, xloc, yloc, zloc, i, x1, y1, z1, opts, opt_seq,
xbound, ybound, zbound;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

zbound:=args[1];
xbound:=args[2];
ybound:=args[3];
a:=op(1, rhs(ybound));
b:=op(2, rhs(ybound));

### WARNING: semantics of type `string` have changed
if not type (rhs(zbound), range) or not (op(1, xbound) = 'x') 
or not (op(1, ybound) = 'y') or not (op(1, zbound) = 'z') 
or not type (evalf(a),numeric) or not type(evalf(b), numeric) then 
ERROR(`input expression is not of dzdxdy type` ) fi;

f:=op(1, rhs(zbound));
g:=op(2, rhs(zbound));
h:=op(1, rhs(xbound));
k:=op(2, rhs(xbound));
x1:=lhs(xbound);
y1:=lhs(ybound);
z1:=lhs(zbound);

f1:=unapply(f,(x1,y1));
g1:=unapply(g,(x1,y1));
h1:=unapply(h,y1);
k1:=unapply(k,y1);

uloc:=u*b+(1-u)*a;
uloc1:=u*h1(a)+(1-u)*k1(a);
uloc2:=u*h1(b)+(1-u)*k1(b);
vloc:=v*k1(uloc)+(1-v)*h1(uloc);
wloc:=v*f1(h1(uloc),uloc)+(1-v)*g1(h1(uloc),uloc);
xloc:=v*f1(k1(uloc),uloc)+(1-v)*g1(k1(uloc),uloc);
yloc:=v*f1(uloc1,a)+(1-v)*g1(uloc1,a);
zloc:=v*f1(uloc2,b)+(1-v)*g1(uloc2,b);

surface:={[vloc,uloc,f1(vloc,uloc)],[vloc,uloc,g1(vloc,uloc)],
[h1(uloc),uloc,wloc], [k1(uloc),uloc,xloc],[uloc1,a,yloc],[uloc2,b,zloc]};

opts:=axes=BOXED, grid =[20,20], scaling=CONSTRAINED, style=PATCH;

if nargs = 3 then
p1:=plot3d(surface, u=0..1, v=0..1, opts);
fi;

if nargs > 3 then
opt_seq:=opts,seq(args[i],`i`=4..nargs);
p1:=plot3d(surface, u=0..1, v=0..1,opt_seq); fi;

plots[display3d](p1);

end :

#5. dydzdx

dydzdxplot:= proc()

local surface, f, f1, g, g1, h, h1, k, k1, p1, a, b, u, v, uloc,
uloc1, uloc2, vloc, wloc, xloc, yloc, zloc, i, x1, y1, z1, opts, opt_seq,
xbound, ybound, zbound;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

ybound:=args[1];
zbound:=args[2];
xbound:=args[3];
a:=op(1, rhs(xbound));
b:=op(2, rhs(xbound));

### WARNING: semantics of type `string` have changed
if not type (rhs(zbound), range) or not (op(1, xbound) = 'x') 
or not (op(1, ybound) = 'y') or not (op(1, zbound) = 'z') 
or not type (evalf(a),numeric) or not type(evalf(b), numeric) then 
ERROR(`input expression is not of dydzdx type` ) fi;

f:=op(1, rhs(ybound));
g:=op(2, rhs(ybound));
h:=op(1, rhs(zbound));
k:=op(2, rhs(zbound));
x1:=lhs(xbound);
y1:=lhs(ybound);
z1:=lhs(zbound);

f1:=unapply(f,(z1,x1));
g1:=unapply(g,(z1,x1));
h1:=unapply(h,x1);
k1:=unapply(k,x1);

uloc:=u*b+(1-u)*a;
uloc1:=u*h1(a)+(1-u)*k1(a);
uloc2:=u*h1(b)+(1-u)*k1(b);
vloc:=v*k1(uloc)+(1-v)*h1(uloc);
wloc:=v*f1(h1(uloc),uloc)+(1-v)*g1(h1(uloc),uloc);
xloc:=v*f1(k1(uloc),uloc)+(1-v)*g1(k1(uloc),uloc);
yloc:=v*f1(uloc1,a)+(1-v)*g1(uloc1,a);
zloc:=v*f1(uloc2,b)+(1-v)*g1(uloc2,b);

surface:={[uloc,f1(vloc,uloc),vloc],[uloc,g1(vloc,uloc),vloc],
[uloc,wloc,h1(uloc)], [uloc,xloc,k1(uloc)],[a,yloc,uloc1],[b,zloc,uloc2]};

opts:=axes=BOXED, grid =[20,20], scaling=CONSTRAINED, style=PATCH;

if nargs = 3 then
p1:=plot3d(surface, u=0..1, v=0..1, opts);
fi;

if nargs > 3 then
opt_seq:=opts,seq(args[i],`i`=4..nargs);
p1:=plot3d(surface, u=0..1, v=0..1,opt_seq); fi;

plots[display3d](p1);

end :

#6. dzdydx

dzdydxplot:= proc()

local surface, f, f1, g, g1, h, h1, k, k1, p1, a, b, u, v, uloc,
uloc1, uloc2, vloc, wloc, xloc, yloc, zloc, i, x1, y1, z1, opts, opt_seq,
xbound, ybound, zbound;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

zbound:=args[1];
ybound:=args[2];
xbound:=args[3];
a:=op(1, rhs(xbound));
b:=op(2, rhs(xbound));

### WARNING: semantics of type `string` have changed
if not type (rhs(zbound), range) or not (op(1, xbound) = 'x') 
or not (op(1, ybound) = 'y') or not (op(1, zbound) = 'z') 
or not type (evalf(a),numeric) or not type(evalf(b), numeric) then 
ERROR(`input expression is not of dzdydx type` ) fi;

f:=op(1, rhs(zbound));
g:=op(2, rhs(zbound));
h:=op(1, rhs(ybound));
k:=op(2, rhs(ybound));
x1:=lhs(xbound);
y1:=lhs(ybound);
z1:=lhs(zbound);

f1:=unapply(f,(y1,x1));
g1:=unapply(g,(y1,x1));
h1:=unapply(h,x1);
k1:=unapply(k,x1);

uloc:=u*b+(1-u)*a;
uloc1:=u*h1(a)+(1-u)*k1(a);
uloc2:=u*h1(b)+(1-u)*k1(b);
vloc:=v*k1(uloc)+(1-v)*h1(uloc);
wloc:=v*f1(h1(uloc),uloc)+(1-v)*g1(h1(uloc),uloc);
xloc:=v*f1(k1(uloc),uloc)+(1-v)*g1(k1(uloc),uloc);
yloc:=v*f1(uloc1,a)+(1-v)*g1(uloc1,a);
zloc:=v*f1(uloc2,b)+(1-v)*g1(uloc2,b);

surface:={[uloc,vloc,f1(vloc,uloc)],[uloc,vloc,g1(vloc,uloc)],
[uloc,h1(uloc),wloc], [uloc,k1(uloc),xloc],[a,uloc1,yloc],[b,uloc2,zloc]};

opts:=axes=BOXED, grid =[20,20], scaling=CONSTRAINED, style=PATCH;

if nargs = 3 then
p1:=plot3d(surface, u=0..1, v=0..1, opts);
fi;

if nargs > 3 then
opt_seq:=opts,seq(args[i],`i`=4..nargs);
p1:=plot3d(surface, u=0..1, v=0..1,opt_seq); fi;

plots[display3d](p1);

end :

#Six triple-integral region plots for domains expressed in cylindrical
#coordinates, r, theta, z.

#1. drd(theta)dz

drdtdzplot:= proc()

local surface, f, f1, g, g1, h, h1, k, k1, a, b, p1, r1, t1, z1, u, v,
uloc, uloc1, uloc2, vloc, wloc, xloc, yloc, zloc, i, opts, opt_seq, rbound,
tbound, zbound;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

rbound:=args[1];
tbound:=args[2];
zbound:=args[3];
a:=op(1, rhs(zbound));
b:=op(2, rhs(zbound));

### WARNING: semantics of type `string` have changed
if not type (rhs(zbound), range) or not (op(1, rbound) = 'r') 
or not (op(1, tbound) = 'theta') or not (op(1, zbound) = 'z') 
or not type (evalf(a),numeric) or not type(evalf(b), numeric) then 
ERROR(`input expression is not of drdtdz type` ) fi;

f:=op(1, rhs(rbound));
g:=op(2, rhs(rbound));
h:=op(1, rhs(tbound));
k:=op(2, rhs(tbound));
r1:=lhs(rbound);
t1:=lhs(tbound);
z1:=lhs(zbound);

f1:=unapply(f,(t1,z1));
g1:=unapply(g,(t1,z1));
h1:=unapply(h,z1);
k1:=unapply(k,z1);

uloc:=u*b+(1-u)*a;
uloc1:=u*h1(a)+(1-u)*k1(a);
uloc2:=u*h1(b)+(1-u)*k1(b);
vloc:=v*k1(uloc)+(1-v)*h1(uloc);
wloc:=v*f1(h1(uloc),uloc)+(1-v)*g1(h1(uloc),uloc);
xloc:=v*f1(k1(uloc),uloc)+(1-v)*g1(k1(uloc),uloc);
yloc:=v*f1(uloc1,a)+(1-v)*g1(uloc1,a);
zloc:=v*f1(uloc2,b)+(1-v)*g1(uloc2,b);

surface:={[f1(vloc,uloc),vloc,uloc],[g1(vloc,uloc),vloc,uloc],
[wloc,h1(uloc),uloc], [xloc,k1(uloc),uloc],[yloc,uloc1,a],[zloc,uloc2,b]};

opts:=axes=BOXED, grid =[20,20], scaling=CONSTRAINED, style=PATCH;

if nargs = 3 then
p1:=plots[cylinderplot](surface, u=0..1, v=0..1, opts); fi;

if nargs > 3 then
opt_seq:=seq(args[i],`i`=4..nargs);
p1:=plots[cylinderplot](surface, u=0..1, v=0..1),opt_seq; fi;

plots[display3d](p1);

end :

#2. dzd(theta)tdr

dzdtdrplot:= proc()

local surface, f, f1, g, g1, h, h1, k, k1, a, b, p1, r1, t1, z1, u, v,
uloc, uloc1, uloc2, vloc, wloc, xloc, yloc, zloc, i, opts, opt_seq, rbound,
tbound, zbound;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

zbound:=args[1];
tbound:=args[2];
rbound:=args[3];
a:=op(1, rhs(rbound));
b:=op(2, rhs(rbound));

### WARNING: semantics of type `string` have changed
if not type (rhs(rbound), range) or not (op(1, rbound) = 'r') 
or not (op(1, tbound) = 'theta') or not (op(1, zbound) = 'z') 
or not type (evalf(a), numeric) or not type(evalf(b), numeric) then 
ERROR(`input expression is not of dzdtdr type` ) fi;

f:=op(1, rhs(zbound));
g:=op(2, rhs(zbound));
h:=op(1, rhs(tbound));
k:=op(2, rhs(tbound));
r1:=lhs(rbound);
t1:=lhs(tbound);
z1:=lhs(zbound);

f1:=unapply(f,(t1,r1));
g1:=unapply(g,(t1,r1));
h1:=unapply(h,r1);
k1:=unapply(k,r1);

uloc:=u*b+(1-u)*a;
uloc1:=u*h1(a)+(1-u)*k1(a);
uloc2:=u*h1(b)+(1-u)*k1(b);
vloc:=v*k1(uloc)+(1-v)*h1(uloc);
wloc:=v*f1(h1(uloc),uloc)+(1-v)*g1(h1(uloc),uloc);
xloc:=v*f1(k1(uloc),uloc)+(1-v)*g1(k1(uloc),uloc);
yloc:=v*f1(uloc1,a)+(1-v)*g1(uloc1,a);
zloc:=v*f1(uloc2,b)+(1-v)*g1(uloc2,b);

surface:={[uloc,vloc,f1(vloc,uloc)],[uloc,vloc,g1(vloc,uloc)],
[uloc,h1(uloc),wloc], [uloc,k1(uloc),xloc],[a,uloc1,yloc],[b,uloc2,zloc]};

opts:=axes=BOXED, grid =[20,20], scaling=CONSTRAINED, style=PATCH;

if nargs = 3 then
p1:=plots[cylinderplot](surface, u=0..1, v=0..1, opts); fi;

if nargs > 3 then
opt_seq:=opts,seq(args[i],`i`=4..nargs);
p1:=plots[cylinderplot](surface, u=0..1, v=0..1,opt_seq); fi;

plots[display3d](p1);

end :

#3. drdzd(theta)

drdzdtplot:= proc()

local surface, f, f1, g, g1, h, h1, k, k1, a, b, p1, r1, t1, z1, u, v,
uloc, uloc1, uloc2, vloc, wloc, xloc, yloc, zloc, i, opts, opt_seq, rbound,
tbound, zbound;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

rbound:=args[1];
zbound:=args[2];
tbound:=args[3];
a:=op(1, rhs(tbound));
b:=op(2, rhs(tbound));

### WARNING: semantics of type `string` have changed
if not type (rhs(rbound), range) or not (op(1, rbound) = 'r') 
or not (op(1, tbound) = 'theta') or not (op(1, zbound) = 'z') 
or not type (evalf(a), numeric) or not type(evalf(b), numeric) then 
ERROR(`input expression is not of drdzdt type` ) fi;

f:=op(1, rhs(rbound));
g:=op(2, rhs(rbound));
h:=op(1, rhs(zbound));
k:=op(2, rhs(zbound));
r1:=lhs(rbound);
t1:=lhs(tbound);
z1:=lhs(zbound);

f1:=unapply(f,(z1,t1));
g1:=unapply(g,(z1,t1));
h1:=unapply(h,t1);
k1:=unapply(k,t1);

uloc:=u*b+(1-u)*a;
uloc1:=u*h1(a)+(1-u)*k1(a);
uloc2:=u*h1(b)+(1-u)*k1(b);
vloc:=v*k1(uloc)+(1-v)*h1(uloc);
wloc:=v*f1(h1(uloc),uloc)+(1-v)*g1(h1(uloc),uloc);
xloc:=v*f1(k1(uloc),uloc)+(1-v)*g1(k1(uloc),uloc);
yloc:=v*f1(uloc1,a)+(1-v)*g1(uloc1,a);
zloc:=v*f1(uloc2,b)+(1-v)*g1(uloc2,b);

surface:={[f1(vloc,uloc),uloc,vloc],[g1(vloc,uloc),uloc,vloc],
[wloc,uloc,h1(uloc)], [xloc,uloc,k1(uloc)],[yloc,a,uloc1],[zloc,b,uloc2]};

opts:=axes=BOXED, grid =[20,20], scaling=CONSTRAINED, style=PATCH;

if nargs = 3 then
p1:=plots[cylinderplot](surface, u=0..1, v=0..1, opts); fi;

if nargs > 3 then
opt_seq:=opts,seq(args[i],`i`=4..nargs);
p1:=plots[cylinderplot](surface, u=0..1, v=0..1,opt_seq); fi;

plots[display3d](p1);

end :

#4. dzdrd(theta)

dzdrdtplot:= proc()

local surface, f, f1, g, g1, h, h1, k, k1, a, b, p1, r1, t1, z1, u, v,
uloc, uloc1, uloc2, vloc, wloc, xloc, yloc, zloc, i, opts, opt_seq, rbound,
tbound, zbound;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

zbound:=args[1];
rbound:=args[2];
tbound:=args[3];
a:=op(1, rhs(tbound));
b:=op(2, rhs(tbound));

### WARNING: semantics of type `string` have changed
if not type (rhs(zbound), range) or not (op(1, rbound) = 'r') 
or not (op(1, tbound) = 'theta') or not (op(1, zbound) = 'z') 
or not type (evalf(a),numeric) or not type(evalf(b), numeric) then 
ERROR(`input expression is not of dzdrdt type` ) fi;

f:=op(1, rhs(zbound));
g:=op(2, rhs(zbound));
h:=op(1, rhs(rbound));
k:=op(2, rhs(rbound));
r1:=lhs(rbound);
t1:=lhs(tbound);
z1:=lhs(zbound);

f1:=unapply(f,(r1,t1));
g1:=unapply(g,(r1,t1));
h1:=unapply(h,t1);
k1:=unapply(k,t1);

uloc:=u*b+(1-u)*a;
uloc1:=u*h1(a)+(1-u)*k1(a);
uloc2:=u*h1(b)+(1-u)*k1(b);
vloc:=v*k1(uloc)+(1-v)*h1(uloc);
wloc:=v*f1(h1(uloc),uloc)+(1-v)*g1(h1(uloc),uloc);
xloc:=v*f1(k1(uloc),uloc)+(1-v)*g1(k1(uloc),uloc);
yloc:=v*f1(uloc1,a)+(1-v)*g1(uloc1,a);
zloc:=v*f1(uloc2,b)+(1-v)*g1(uloc2,b);

surface:={[vloc,uloc,f1(vloc,uloc)],[vloc,uloc,g1(vloc,uloc)],
[h1(uloc),uloc,wloc], [k1(uloc),uloc,xloc],[uloc1,a,yloc],[uloc2,b,zloc]};

opts:=axes=BOXED, grid =[20,20], scaling=CONSTRAINED, style=PATCH;

if nargs = 3 then
p1:=plots[cylinderplot](surface, u=0..1, v=0..1, opts); fi;

if nargs > 3 then
opt_seq:=opts,seq(args[i],`i`=4..nargs);
p1:=plots[cylinderplot](surface, u=0..1, v=0..1,opt_seq); fi;

plots[display3d](p1);

end :

#5. d(theta)dzdr

dtdzdrplot:= proc()

local surface, f, f1, g, g1, h, h1, k, k1, a, b, p1, r1, t1, z1, u, v,
uloc, uloc1, uloc2, vloc, wloc, xloc, yloc, zloc, i, opts, opt_seq, rbound,
tbound, zbound;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

tbound:=args[1];
zbound:=args[2];
rbound:=args[3];
a:=op(1, rhs(rbound));
b:=op(2, rhs(rbound));

### WARNING: semantics of type `string` have changed
if not type (rhs(zbound), range) or not (op(1, rbound) = 'r') 
or not (op(1, tbound) = 'theta') or not (op(1, zbound) = 'z') 
or not type (evalf(a),numeric) or not type(evalf(b), numeric) then 
ERROR(`input expression is not of dtdzdr type` ) fi;

f:=op(1, rhs(tbound));
g:=op(2, rhs(tbound));
h:=op(1, rhs(zbound));
k:=op(2, rhs(zbound));
r1:=lhs(rbound);
t1:=lhs(tbound);
z1:=lhs(zbound);

f1:=unapply(f,(z1,r1));
g1:=unapply(g,(z1,r1));
h1:=unapply(h,r1);
k1:=unapply(k,r1);

uloc:=u*b+(1-u)*a;
uloc1:=u*h1(a)+(1-u)*k1(a);
uloc2:=u*h1(b)+(1-u)*k1(b);
vloc:=v*k1(uloc)+(1-v)*h1(uloc);
wloc:=v*f1(h1(uloc),uloc)+(1-v)*g1(h1(uloc),uloc);
xloc:=v*f1(k1(uloc),uloc)+(1-v)*g1(k1(uloc),uloc);
yloc:=v*f1(uloc1,a)+(1-v)*g1(uloc1,a);
zloc:=v*f1(uloc2,b)+(1-v)*g1(uloc2,b);

surface:={[uloc,f1(vloc,uloc),vloc],[uloc,g1(vloc,uloc),vloc],
[uloc,wloc,h1(uloc)], [uloc,xloc,k1(uloc)],[a,yloc,uloc1],[b,zloc,uloc2]};

opts:=axes=BOXED, grid =[20,20], scaling=CONSTRAINED, style=PATCH;

if nargs = 3 then
p1:=plots[cylinderplot](surface, u=0..1, v=0..1, opts); fi;

if nargs > 3 then
opt_seq:=opts,seq(args[i],`i`=4..nargs);
p1:=plots[cylinderplot](surface, u=0..1, v=0..1,opt_seq); fi;

plots[display3d](p1);

end :


#6. d(theta)drdz

dtdrdzplot:= proc()

local surface, f, f1, g, g1, h, h1, k, k1, a, b, p1, r1, t1, z1, u, v,
uloc, uloc1, uloc2, vloc, wloc, xloc, yloc, zloc, i, opts, opt_seq, rbound,
tbound, zbound;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

tbound:=args[1];
rbound:=args[2];
zbound:=args[3];
a:=op(1, rhs(zbound));
b:=op(2, rhs(zbound));

### WARNING: semantics of type `string` have changed
if not type (rhs(zbound), range) or not (op(1, rbound) = 'r') 
or not (op(1, tbound) = 'theta') or not (op(1, zbound) = 'z') 
or not type (evalf(a),numeric) or not type(evalf(b), numeric) then 
ERROR(`input expression is not of dtdrdz type` ) fi;

f:=op(1, rhs(tbound));
g:=op(2, rhs(tbound));
h:=op(1, rhs(rbound));
k:=op(2, rhs(rbound));
r1:=lhs(rbound);
t1:=lhs(tbound);
z1:=lhs(zbound);

f1:=unapply(f,(r1,z1));
g1:=unapply(g,(r1,z1));
h1:=unapply(h,z1);
k1:=unapply(k,z1);

uloc:=u*b+(1-u)*a;
uloc1:=u*h1(a)+(1-u)*k1(a);
uloc2:=u*h1(b)+(1-u)*k1(b);
vloc:=v*k1(uloc)+(1-v)*h1(uloc);
wloc:=v*f1(h1(uloc),uloc)+(1-v)*g1(h1(uloc),uloc);
xloc:=v*f1(k1(uloc),uloc)+(1-v)*g1(k1(uloc),uloc);
yloc:=v*f1(uloc1,a)+(1-v)*g1(uloc1,a);
zloc:=v*f1(uloc2,b)+(1-v)*g1(uloc2,b);

surface:={[vloc,f1(vloc,uloc),uloc],[vloc,g1(vloc,uloc),uloc],
[h1(uloc),wloc,uloc], [k1(uloc),xloc,uloc],[uloc1,yloc,a],[uloc2,zloc,b]};

opts:=axes=BOXED, grid =[20,20], scaling=CONSTRAINED, style=PATCH;

if nargs = 3 then
p1:=plots[cylinderplot](surface, u=0..1, v=0..1, opts); fi;

if nargs > 3 then
opt_seq:=opts,seq(args[i],`i`=4..nargs);
p1:=plots[cylinderplot](surface, u=0..1, v=0..1,opt_seq); fi;

plots[display3d](p1);

end :

#Six spherical coordinate region plotters. The spherical coordinates
#are rho (radius), theta (longitudinal angle; i.e., angle made by
#(x,y,0) and x-axis), and phi (latitudinal angle; i.e., angle made by
#(x,y,z) and positive z-axis).

#1. d(rho)d(theta)d(phi)

dpdtdphiplot:= proc()

local surface, f, f1, g, g1, h, h1, k, k1, a, b, p1, q1, t1, phi1, u,
v, uloc, uloc1, uloc2, vloc, wloc, xloc, yloc, zloc, i, opts, opt_seq,
pbound, tbound, phibound;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

pbound:=args[1];
tbound:=args[2];
phibound:=args[3];
a:=op(1, rhs(phibound));
b:=op(2, rhs(phibound));

### WARNING: semantics of type `string` have changed
if not type (rhs(phibound), range) or not (op(1, pbound) = 'rho') 
or not (op(1, tbound) = 'theta') or not (op(1, phibound) = 'phi') 
or not type (evalf(a),numeric) or not type(evalf(b), numeric) then 
ERROR(`input expression is not of dpdtdphi type` ) fi;

f:=op(1, rhs(pbound));
g:=op(2, rhs(pbound));
h:=op(1, rhs(tbound));
k:=op(2, rhs(tbound));
p1:=lhs(pbound);
t1:=lhs(tbound);
phi1:=lhs(phibound);

f1:=unapply(f,(t1,phi1));
g1:=unapply(g,(t1,phi1));
h1:=unapply(h,phi1);
k1:=unapply(k,phi1);

uloc:=u*b+(1-u)*a;
uloc1:=u*h1(a)+(1-u)*k1(a);
uloc2:=u*h1(b)+(1-u)*k1(b);
vloc:=v*k1(uloc)+(1-v)*h1(uloc);
wloc:=v*f1(h1(uloc),uloc)+(1-v)*g1(h1(uloc),uloc);
xloc:=v*f1(k1(uloc),uloc)+(1-v)*g1(k1(uloc),uloc);
yloc:=v*f1(uloc1,a)+(1-v)*g1(uloc1,a);
zloc:=v*f1(uloc2,b)+(1-v)*g1(uloc2,b);

surface:={[f1(vloc,uloc),vloc,uloc],[g1(vloc,uloc),vloc,uloc],
[wloc,h1(uloc),uloc], [xloc,k1(uloc),uloc],[yloc,uloc1,a],[zloc,uloc2,b]};

opts:=axes=BOXED, grid =[20,20], scaling=CONSTRAINED, style=PATCH;

if nargs =3 then
q1:=plots[sphereplot](surface, u=0..1, v=0..1, opts); fi;

if nargs > 3 then
opt_seq:=opts,seq(args[i],`i`=4..nargs);
q1:=plots[sphereplot](surface, u=0..1, v=0..1,opt_seq); fi;

plots[display3d](q1);

end:

#2. d(rho)d(phi)d(theta)

dpdphidtplot:= proc()

local surface, f, f1, g, g1, h, h1, k, k1, a, b, p1, q1, t1, phi1, u,
v, uloc, uloc1, uloc2, vloc, wloc, xloc, yloc, zloc, i, opts, opt_seq,
pbound, tbound, phibound;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

pbound:=args[1];
phibound:=args[2];
tbound:=args[3];
a:=op(1, rhs(tbound));
b:=op(2, rhs(tbound));

### WARNING: semantics of type `string` have changed
if not type (rhs(phibound), range) or not (op(1, pbound) = 'rho') 
or not (op(1, tbound) = 'theta') or not (op(1, phibound) = 'phi') 
or not type (evalf(a),numeric) or not type(evalf(b), numeric) then 
ERROR(`input expression is not of dpdphidt type` ) fi;

f:=op(1, rhs(pbound));
g:=op(2, rhs(pbound));
h:=op(1, rhs(phibound));
k:=op(2, rhs(phibound));
p1:=lhs(pbound);
t1:=lhs(tbound);
phi1:=lhs(phibound);

f1:=unapply(f,(phi1,t1));
g1:=unapply(g,(phi1,t1));
h1:=unapply(h,t1);
k1:=unapply(k,t1);

uloc:=u*b+(1-u)*a;
uloc1:=u*h1(a)+(1-u)*k1(a);
uloc2:=u*h1(b)+(1-u)*k1(b);
vloc:=v*k1(uloc)+(1-v)*h1(uloc);
wloc:=v*f1(h1(uloc),uloc)+(1-v)*g1(h1(uloc),uloc);
xloc:=v*f1(k1(uloc),uloc)+(1-v)*g1(k1(uloc),uloc);
yloc:=v*f1(uloc1,a)+(1-v)*g1(uloc1,a);
zloc:=v*f1(uloc2,b)+(1-v)*g1(uloc2,b);

surface:={[f1(vloc,uloc),uloc,vloc],[g1(vloc,uloc),uloc,vloc],
[wloc,uloc,h1(uloc)], [xloc,uloc,k1(uloc)],[yloc,a,uloc1],[zloc,b,uloc2]};

opts:=axes=BOXED, grid =[20,20], scaling=CONSTRAINED, style=PATCH;

if nargs =3 then
q1:=plots[sphereplot](surface, u=0..1, v=0..1, opts); fi;

if nargs > 3 then
opt_seq:=opts,seq(args[i],`i`=4..nargs);
q1:=plots[sphereplot](surface, u=0..1, v=0..1,opt_seq); fi;

plots[display3d](q1);

end:

#3. d(theta)d(phi)d(rho)

dtdphidpplot:= proc()

local surface, f, f1, g, g1, h, h1, k, k1, a, b, p1, q1, t1, phi1, u,
v, uloc, uloc1, uloc2, vloc, wloc, xloc, yloc, zloc, i, opts, opt_seq,
pbound, tbound, phibound;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

tbound:=args[1];
phibound:=args[2];
pbound:=args[3];
a:=op(1, rhs(pbound));
b:=op(2, rhs(pbound));

### WARNING: semantics of type `string` have changed
if not type (rhs(phibound), range) or not (op(1, pbound) = 'rho') 
or not (op(1, tbound) = 'theta') or not (op(1, phibound) = 'phi') 
or not type (evalf(a),numeric) or not type(evalf(b), numeric) then 
ERROR(`input expression is not of dtdphidp type` ) fi;

f:=op(1, rhs(tbound));
g:=op(2, rhs(tbound));
h:=op(1, rhs(phibound));
k:=op(2, rhs(phibound));
p1:=lhs(pbound);
t1:=lhs(tbound);
phi1:=lhs(phibound);

f1:=unapply(f,(phi1,p1));
g1:=unapply(g,(phi1,p1));
h1:=unapply(h,p1);
k1:=unapply(k,p1);

uloc:=u*b+(1-u)*a;
uloc1:=u*h1(a)+(1-u)*k1(a);
uloc2:=u*h1(b)+(1-u)*k1(b);
vloc:=v*k1(uloc)+(1-v)*h1(uloc);
wloc:=v*f1(h1(uloc),uloc)+(1-v)*g1(h1(uloc),uloc);
xloc:=v*f1(k1(uloc),uloc)+(1-v)*g1(k1(uloc),uloc);
yloc:=v*f1(uloc1,a)+(1-v)*g1(uloc1,a);
zloc:=v*f1(uloc2,b)+(1-v)*g1(uloc2,b);

surface:={[uloc,f1(vloc,uloc),vloc],[uloc,g1(vloc,uloc),vloc],
[uloc,wloc,h1(uloc)], [uloc,xloc,k1(uloc)],[a,yloc,uloc1],[b,zloc,uloc2]};

opts:=axes=BOXED, grid =[20,20], scaling=CONSTRAINED, style=PATCH;

if nargs =3 then
q1:=plots[sphereplot](surface, u=0..1, v=0..1, opts); fi;

if nargs > 3 then
opt_seq:=opts,seq(args[i],`i`=4..nargs);
q1:=plots[sphereplot](surface, u=0..1, v=0..1,opt_seq); fi;

plots[display3d](q1);

end:

#4. d(phi)d(theta)d(rho)

dphidtdpplot:= proc()

local surface, f, f1, g, g1, h, h1, k, k1, a, b, p1, q1, t1, phi1, u,
v, uloc, uloc1, uloc2, vloc, wloc, xloc, yloc, zloc, i, opts, opt_seq,
pbound, tbound, phibound;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

phibound:=args[1];
tbound:=args[2];
pbound:=args[3];
a:=op(1, rhs(pbound));
b:=op(2, rhs(pbound));

### WARNING: semantics of type `string` have changed
if not type (rhs(phibound), range) or not (op(1, pbound) = 'rho') 
or not (op(1, tbound) = 'theta') or not (op(1, phibound) = 'phi') 
or not type (evalf(a),numeric) or not type(evalf(b), numeric) then 
ERROR(`input expression is not of dphidtdp type` ) fi;

f:=op(1, rhs(phibound));
g:=op(2, rhs(phibound));
h:=op(1, rhs(tbound));
k:=op(2, rhs(tbound));
p1:=lhs(pbound);
t1:=lhs(tbound);
phi1:=lhs(phibound);

f1:=unapply(f,(t1,p1));
g1:=unapply(g,(t1,p1));
h1:=unapply(h,p1);
k1:=unapply(k,p1);

uloc:=u*b+(1-u)*a;
uloc1:=u*h1(a)+(1-u)*k1(a);
uloc2:=u*h1(b)+(1-u)*k1(b);
vloc:=v*k1(uloc)+(1-v)*h1(uloc);
wloc:=v*f1(h1(uloc),uloc)+(1-v)*g1(h1(uloc),uloc);
xloc:=v*f1(k1(uloc),uloc)+(1-v)*g1(k1(uloc),uloc);
yloc:=v*f1(uloc1,a)+(1-v)*g1(uloc1,a);
zloc:=v*f1(uloc2,b)+(1-v)*g1(uloc2,b);

surface:={[uloc,vloc,f1(vloc,uloc)],[uloc,vloc,g1(vloc,uloc)],
[uloc,h1(uloc),wloc], [uloc,k1(uloc),xloc],[a,uloc1,yloc],[b,uloc2,zloc]};

opts:=axes=BOXED, grid =[20,20], scaling=CONSTRAINED, style=PATCH;

if nargs =3 then
q1:=plots[sphereplot](surface, u=0..1, v=0..1, opts); fi;

if nargs > 3 then
opt_seq:=opts,seq(args[i],`i`=4..nargs);
q1:=plots[sphereplot](surface, u=0..1, v=0..1,opt_seq); fi;

plots[display3d](q1);

end:

#5. d(theta)d(rho)d(phi)

dtdpdphiplot:= proc()

local surface, f, f1, g, g1, h, h1, k, k1, a, b, p1, q1, t1, phi1, u,
v, uloc, uloc1, uloc2, vloc, wloc, xloc, yloc, zloc, i, opts, opt_seq,
pbound, tbound, phibound;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

tbound:=args[1];
pbound:=args[2];
phibound:=args[3];
a:=op(1, rhs(pbound));
b:=op(2, rhs(pbound));

### WARNING: semantics of type `string` have changed
if not type (rhs(phibound), range) or not (op(1, pbound) = 'rho') 
or not (op(1, tbound) = 'theta') or not (op(1, phibound) = 'phi') 
or not type (evalf(a),numeric) or not type(evalf(b), numeric) then 
ERROR(`input expression is not of dtdpdphi type` ) fi;

f:=op(1, rhs(tbound));
g:=op(2, rhs(tbound));
h:=op(1, rhs(phibound));
k:=op(2, rhs(phibound));
p1:=lhs(pbound);
t1:=lhs(tbound);
phi1:=lhs(phibound);

f1:=unapply(f,(p1,phi1));
g1:=unapply(g,(p1,phi1));
h1:=unapply(h,p1);
k1:=unapply(k,p1);

uloc:=u*b+(1-u)*a;
uloc1:=u*h1(a)+(1-u)*k1(a);
uloc2:=u*h1(b)+(1-u)*k1(b);
vloc:=v*k1(uloc)+(1-v)*h1(uloc);
wloc:=v*f1(h1(uloc),uloc)+(1-v)*g1(h1(uloc),uloc);
xloc:=v*f1(k1(uloc),uloc)+(1-v)*g1(k1(uloc),uloc);
yloc:=v*f1(uloc1,a)+(1-v)*g1(uloc1,a);
zloc:=v*f1(uloc2,b)+(1-v)*g1(uloc2,b);

surface:={[vloc,f1(vloc,uloc),uloc],[vloc,g1(vloc,uloc),uloc],
[h1(uloc),wloc,uloc], [k1(uloc),xloc,uloc],[uloc1,yloc,a],[uloc2,zloc,b]};

opts:=axes=BOXED, grid =[20,20], scaling=CONSTRAINED, style=PATCH;

if nargs =3 then
q1:=plots[sphereplot](surface, u=0..1, v=0..1, opts); fi;

if nargs > 3 then
opt_seq:=opts,seq(args[i],`i`=4..nargs);
q1:=plots[sphereplot](surface, u=0..1, v=0..1,opt_seq); fi;

plots[display3d](q1);

end:

#6. d(phi)d(rho)d(theta)

dphidpdtplot:= proc()

local surface, f, f1, g, g1, h, h1, k, k1, a, b, p1, q1, t1, phi1, u,
v, uloc, uloc1, uloc2, vloc, wloc, xloc, yloc, zloc, i, opts, opt_seq,
pbound, tbound, phibound;

if nargs < 3 then
ERROR(`there must be at least three arguments`) fi;

phibound:=args[1];
pbound:=args[2];
tbound:=args[3];
a:=op(1, rhs(tbound));
b:=op(2, rhs(tbound));

### WARNING: semantics of type `string` have changed
if not type (rhs(phibound), range) or not (op(1, pbound) = 'rho') 
or not (op(1, tbound) = 'theta') or not (op(1, phibound) = 'phi') 
or not type (evalf(a),numeric) or not type(evalf(b), numeric) then 
ERROR(`input expression is not of dphidpdt type` ) fi;

f:=op(1, rhs(phibound));
g:=op(2, rhs(phibound));
h:=op(1, rhs(pbound));
k:=op(2, rhs(pbound));
p1:=lhs(pbound);
t1:=lhs(tbound);
phi1:=lhs(phibound);

f1:=unapply(f,(p1,t1));
g1:=unapply(g,(p1,t1));
h1:=unapply(h,t1);
k1:=unapply(k,t1);

uloc:=u*b+(1-u)*a;
uloc1:=u*h1(a)+(1-u)*k1(a);
uloc2:=u*h1(b)+(1-u)*k1(b);
vloc:=v*k1(uloc)+(1-v)*h1(uloc);
wloc:=v*f1(h1(uloc),uloc)+(1-v)*g1(h1(uloc),uloc);
xloc:=v*f1(k1(uloc),uloc)+(1-v)*g1(k1(uloc),uloc);
yloc:=v*f1(uloc1,a)+(1-v)*g1(uloc1,a);
zloc:=v*f1(uloc2,b)+(1-v)*g1(uloc2,b);

surface:={[vloc,uloc,f1(vloc,uloc)],[vloc,uloc,g1(vloc,uloc)],
[h1(uloc),uloc,wloc], [k1(uloc),uloc,xloc],[uloc1,a,yloc],[uloc2,b,zloc]};

opts:=axes=BOXED, grid =[20,20], scaling=CONSTRAINED, style=PATCH;

if nargs =3 then
q1:=plots[sphereplot](surface, u=0..1, v=0..1, opts); fi;

if nargs > 3 then
opt_seq:=opts,seq(args[i],`i`=4..nargs);
q1:=plots[sphereplot](surface, u=0..1, v=0..1,opt_seq); fi;

plots[display3d](q1);

end:

#regionplot3d plots the image of a three-dimensional region of the type
#h(u,v)<=w<=k(u,v),p(u)<=v<=q(u),a<=u<=b under a coordinate transformation
#[x(u,v,w),y(u,v,w),z(u,v,w)].

regionplot3d:=proc()
local backface, frontface, leftface, rightface, ulface, face1, face2, face3, face4, face5, face6, j, h, k, p, q,a, b, t, s, opts, transfun, ubound, vbound, wbound, xtrans, ytrans, ztrans, u1, v1, w1,p1,p2;

if nargs<4 then ERROR(`three arguments are required`) fi;
if not type(args[4],list) and nargs(args[4])<>3 then ERROR(`fourth argument must be a list of three coordinate transformations`);
### WARNING: semantics of type `string` have changed
elif not type(rhs(args[3]), range) or not type(evalf(op(1,rhs(args[3]))),numeric) or not type(evalf(op(2,rhs(args[3]))), numeric) then 
ERROR(`input expression is not of 3d region type` )
fi;

wbound:=args[1];
vbound:=args[2];
ubound:=args[3];
transfun:=args[4];
a:=op(1,rhs(ubound));
b:=op(2,rhs(ubound));
u1:=lhs(ubound);
v1:=lhs(vbound);
w1:=lhs(wbound);

p:=unapply(op(1,rhs(vbound)),u1);
q:=unapply(op(2,rhs(vbound)),u1);
h:=unapply(op(1,rhs(wbound)),u1,v1);
k:=unapply(op(2,rhs(wbound)),u1,v1);
xtrans:=unapply(transfun[1],u1,v1,w1);
ytrans:=unapply(transfun[2],u1,v1,w1);
ztrans:=unapply(transfun[3],u1,v1,w1);

ulface:=t*q(u1)+(1-t)*p(u1);
frontface:=t*k(u1,p(u1))+(1-t)*h(u1,p(u1));
backface:=t*k(u1,q(u1))+(1-t)*h(u1,q(u1));
leftface:=s*k(a,t*q(a)+(1-t)*p(a))+(1-s)*h(a,t*q(a)+(1-t)*p(a));
rightface:=s*k(b,t*q(b)+(1-t)*p(b))+(1-s)*h(b,t*q(b)+(1-t)*p(b));

face1:=[xtrans(u1,ulface,h(u1,ulface)),ytrans(u1,ulface,h(u1,ulface)),ztrans(u1,ulface,h(u1,ulface))];
face2:=[xtrans(u1,ulface,k(u1,ulface)),ytrans(u1,ulface,k(u1,ulface)),ztrans(u1,ulface,k(u1,ulface))];
face3:=[xtrans(u1,p(u1),frontface),ytrans(u1,p(u1),frontface),ztrans(u1,p(u1),frontface)];
face4:=[xtrans(u1,q(u1),backface),ytrans(u1,q(u1),backface),ztrans(u1,q(u1),backface)];
face5:=[xtrans(a,t*q(a)+(1-t)*p(a),leftface),ytrans(a,t*q(a)+(1-t)*p(a),leftface),ztrans(a,t*q(a)+(1-t)*p(a),leftface)];
face6:=[xtrans(b,t*q(b)+(1-t)*p(b),rightface),ytrans(b,t*q(b)+(1-t)*p(b),rightface),ztrans(b,t*q(b)+(1-t)*p(b),rightface)];

if nargs=4 then
opts:=scaling=constrained,axes=boxed,style=patch,labels=['x','y','z'];
else
opts:=scaling=constrained,axes=boxed,style=patch,labels=['x','y','z'],args[5..nargs];
fi;

p1:=plot3d({face1,face2,face3,face4},u1=a..b,t=0..1,opts):
p2:=plot3d({face5,face6},s=0..1,t=0..1,opts):
plots[display]({p1,p2},opts);

end:

#Begin the on-line help files.

`help/text/calcplot` := TEXT(`HELP FOR: Introduction to the calcplot package.`,
`   `,
`CALLING SEQUENCE:`,
`   `,
`	<function>(args)`,
`   `,
`SYNOPSIS:   `,
`   `,
`-To use a calcplot function, invoke the function with a call <function>(args)\
,`,
`where <function> is the name of the function, and args are the appropriate `,
`arguments for the function.`,
`   `,
`-The functions available are:`,
`	rotxplot rotyplot`,
`	dxdyplot dydxplot`,
`	drdtplot dtdrplot`,
`	xygraphplot yxgraphplot`,
`	rtgraphplot trgraphplot`,
`	dxdydzplot dxdzdyplot dydxdzplot dydzdxplot dzdxdyplot dzdydxplot`,
`	drdtdzplot drdzdtplot dtdrdzplot dtdzdrplot dzdrdtplot dzdtdrplot`,
`	dpdtdphiplot dpdphidtplot dtdpdphiplot dtdphidpplot dphidpdtplot`,
`	dtdphidpplot regionplot2d regionplot3d`,
`   `,
`-For more information on a particular function, invoke help for that function`,` with the command ?<function>`,`  `):

`help/text/rotxplot`:=TEXT(`HELP FOR: rotxplot and rotyplot`,
`   `,
`CALLING SEQUENCE: `,
`	rotxplot(f,x=a..b,y=c, opts)`,
`	rotyplot(f,x=a..b,x=c, opts)`,
`   `,
`PARAMETERS:`,
`	f is either an expression in the variable x,`,
`	a procedure, or a two-element list.`,`   `,
`	a, b, c must be symbolic constants, such as Pi and E, integers,`,
`	and real floating point numbers; i.e., evalf(a), evalf(b) must`,
`	be of type numeric.`,
`	`,
`	opts is a sequence of plot3d options.`,
`   `,
`SYNOPSIS:   `,`   `,
`-The rotxplot, rotyplot procedures plot the surface of rotation obtained`,
` by rotating the curve y=f(x), for x from a to b, about the axis given by`,
` the third argument of the function. The procedures return plot3d structures`,
` accept plot3d options, and may be passed to the display3d procedure from`,
` the plots package.`,
`	`,
`EXAMPLES:   `,
`   `,
`rotxplot(sin(x),x=0..2*Pi,y=0);`,
`	`,
`rotyplot(x^2,x=0..1,x=1);`,
`	`,
`rotxplot(x->sin(x),x=0..2*Pi,y=-1, style=PATCHCONTOUR,`,
`orientation=[-25,71], title=``A whale of a surface``);`,
`	`,
`SEE ALSO: plot3d[options], plots[display3d]`):

`help/text/rotyplot`:=`help/text/rotxplot`:

`help/text/dxdyplot`:=TEXT(`HELP FOR: dxdyplot and dydxplot`,
`   `,
`CALLING SEQUENCE: `,`	`,
`	dxdyplot(x=g(y)..k(y), y=a..b, opts)`,
`	dydxplot(y=h(x)..m(x), x=a..b, opts)`,
`   `,
`PARAMETERS:`,
`	The variable names ``x' and ``y' must be used,`,
`   `,
`	g(y), k(y) are expressions in the variable y,`,
`	h(x), m(x) are expressions in the variable x,`,
`   `,
`	a, b must be symbolic constants, such as Pi and E, integers,`,
`	or real floating point numbers; i.e., evalf(a), evalf(b) must`,
`	be of type numeric.`,
`	`,
`	opts is a sequence of plot options.`,
`	`,
`SYNOPSIS:   `,`   `,
`-The dxdyplot, dydxplot procedures plot regions in the plane with the given `,
` arguments as boundary curves. The procedures draw three line segments that `,
` pass through the interior of the region. The syntax is the same as if int`,
` procedure were invoked for a double integral over the region. The procedures`,
` return plot structures, accept plot options, and may be passed to the display`,
` procedure from the plots package.`,
`	`,
`EXAMPLES:   `,
`   `,
`dxdyplot(x=0..sqrt(y),y=0..1);`,
`	`,
`dydxplot(y=x^2..1,x=0..1);`,
`	`,
`dydxplot(y=0..2*x+1,x=0..1,view=[-1..2,0..4],title=``A trapezoidal region``);`,
`	`,
`SEE ALSO: plot[options], plots[display], xygraphplot, yxgraphplot`):

`help/text/dydxplot`:=`help/text/dxdyplot`:

`help/text/drdtplot`:=TEXT(`HELP FOR: drdtplot and dtdrplot`,
`   `,
`CALLING SEQUENCE: `,
`	`,
`	drdtplot(r=f(theta)..g(theta),theta=a..b, opts)`,
`	dtdrplot(theta=h(r)..k(r),r=a..b,opts)`,
`   `,
`PARAMETERS:`,
`	`,
`	The variable names ``r' and ``theta' must be used,`,
`   `,
`	f(theta), g(theta) are expressions in the variable theta,`,
`	h(r), k(r) are expressions in the variable r,`,
`   `,
`	a, b must be symbolic constants, such as Pi and E, integers,`,
`	or real floating point numbers; i.e., evalf(a), evalf(b) must`,
`	be of type numeric.`,
`	`,
`	opts is a sequence of plot options.`,
`   `,
`SYNOPSIS:  `,
`   `,
`-The drdtplot, dtdrplot procedures plot regions in the plane which are `,
` described in polar coordinates. The procedures draw four line segments that `,
` pass through the interior of the region. The syntax is the same as if int`,
` procedure were invoked for a double integral over the region in polar `,
` coordinates. The procedures return plot structures, accept plot options,`, 
` and may be passed to the display procedure from the plots package.`,
`   `,
`EXAMPLES:   `,
`   `,
`drdtplot(r=1..1+cos(theta),theta=-Pi/2..Pi/2);`,
`	`,
`dtdrplot(theta=0..r,r=Pi/2..Pi);`,
`	`,
`plots[display]({drdtplot(r=sin(3*theta)..sin(theta),theta=Pi/4..Pi/3), `,
`drdtplot(r=0..sin(theta),theta=Pi/3..Pi/2)}, scaling=CONSTRAINED,`,
`title=``A region between r=sin(3*theta) and r=sin(theta)``);`,
`	`,
`SEE ALSO: plot[options], plots[dislay], rtgraphplot, trgraphplot`):

`help/text/dtdrplot`:=`help/text/drdtplot`:

`help/text/regionplot2d`:=TEXT(`HELP FOR: regionplot2d`,
`   `,
`CALLING SEQUENCE: `,
`	`,
`	regionplot2d(v=h(u)..k(u), u=a..b, [x(u,v),y(u,v)], opts)`,
`	`,
`PARAMETERS: `,
`	`,
`	g(u), k(u) are expressions in the variable v,`,
`	[x(u,v),y(u,v)] is a list of expressions in the variables u and v,`,
`	`,
`	a, b must be symbolic constants, such as Pi and E, integers,`,
`	or real floating point numbers; i.e., evalf(a), evalf(b) must`,
`	be of type numeric.`,
`	`,
`	opts is a sequence of plot options.`,
`	`,
`SYNOPSIS:   `,
`	`,
`-The regionplot2d procedure plots the image of a region in the uv-plane of the`,
` form h(u)<= v <= k(u), a<= u <= b under a coordinate transformation given by`,
` the list [x(u,v),y(u,v)]. The procedure also draws four segments passing`,
` through the interior of the region. The variable order is the same as if the`,
` int procedure were invoked for a double integral over the region. The procedure`,
` returns a plot structure, accepts plot options, and may be passed to the`,
` display procedure from the plots package.`,
`	`,
`EXAMPLES:   `,
`	`,
`regionplot2d(r=1..3,t=0..2*Pi,[r*cos(t),r*sin(t)],title=``Polar Coordinates``);`,
`	`,
`regionplot2d(u=v^2..v,v=0..1,[u+v,v-u]);`,
`	`,
`regionplot2d(s=0..1,t=0..1,[s-s*t,s*t],color=blue,title=``A Triangle``);`,
`	`,
`SEE ALSO: plot[options], plots[display]`):

`help/text/xygraphplot` := TEXT(`HELP FOR: xygraphplot and yxgraphplot`,
`   `,
`CALLING SEQUENCE:`,
`   `,
`	xygraphplot(f, x=g(y)..k(y), y=a..b, opts)`,
`	yxgraphplot(f, y=h(x)..m(x), x=a..b, opts)`,
`   `,
`PARAMETERS:`,
`	The variable names ``x' and ``y' must be used,`,
`   `,
`	f must be an expression in x and y or a procedure accepting two`,
`	arguments.`,
`   `,
`	g(y), k(y) are  expressions in the variable y,`,
`	h(x), m(x) are  expressions in the variable x,`,
`   `,
`	a, b must be symbolic constants, such as Pi and E, integers,`,
`	and real floating point numbers; i.e., evalf(a), evalf(b) must`,
`	be of type numeric.`,
`	`,
`	opts is a sequence of plot3d options.`,
`   `,
`SYNOPSIS:   `,
`   `,
`-The xygraphplot, yxgraphplot procedures plot graphs z=f(x,y) regions`,
` in the plane with the given arguments as boundary curves. The syntax`,
` is similar to the dxdyplot, dydxplot procedures. The xygraphplot,`,
` yxgraphplot procedures return plot3d structures and so may be passed to`,
` the display3d procedure from the plots package.`,
`   `,
`EXAMPLES:   `,
`	`,
`xygraphplot(x*y,x=-sqrt(1-y^2)..sqrt(1-y^2),y=0..1);`,
`   `,
`yxgraphplot(x*y,y=0..sqrt(1-x^2),x=-1..1);`,
`   `,
`yxgraphplot((x,y)->sin(Pi*(x^2+y^2)),y=x^2..x^(1/4),x=0..1,`
,`style=PATCHCONTOUR, title=``A slice of sine``);`,
`   `,
`SEE ALSO: plot3d[options], plots[display3d], dxdyplot, dydxplot`):

`help/text/yxgraphplot`:=`help/text/xygraphplot`:

`help/text/rtgraphplot`:=TEXT(`HELP FOR: rtgraphplot and trgraphplot`,
`   `,
`CALLING SEQUENCE:`,`   `,
`	rtgraphplot(f, r=g(theta)..j(theta), theta=a..b, opts)`,
`	trgraphplot(f, theta=h(r)..k(r), r=a..b, opts)`,
`   `,
`PARAMETERS:`,
`	The variable names ``r' and ``theta' must be used,`,
`   `,
`	f is an expression in ``r' and ``theta' or a procedure`,
`	that takes two arguments. The arguments will be treated as (r,theta)`,
`	in that order.`,
`   `,
`	g(theta), j(theta) are expressions in the variable theta,`,
`	h(r), k(r) are expressions in the variable r,`,
`   `,
`	a, b must be symbolic constants, such as Pi and E, integers,`,
`	or real floating point numbers; i.e., evalf(a), evalf(b) must`,
`	be of type numeric.`,
`	`,
`	opts is a sequence of plot3d options.`,
`   `,
`SYNOPSIS:   `,
`   `,
`-The rtgraphplot, trgraphplot procedures plot graphs z = f(r,theta) over `,
` regions in the plane which are described in polar coordinates. `,
` The syntax is similar to the drdtplot, dtdrplot procedures. The rtgraphplot, `,
` trgraphplot procedures return plot3d structures, accept plot3d options, and`,
` may be passed to the display3d procedure from the plots package.`,
`   `,
`EXAMPLES:   `,
`   `,
`rtgraphplot(r^2*sin(theta)*cos(theta),r=1..1+cos(theta),theta=-Pi/2..Pi/2);`,
`   `,
`trgraphplot(theta*log(r),theta=0..r,r=Pi/2..Pi);`,
`   `,
`p1:=rtgraphplot(r^2,r=1..1+cos(theta), theta=-Pi/2..Pi/2):`,
`p2:=rtgraphplot(-1,r=1..1+cos(theta), theta=-Pi/2..Pi/2):`,
`plots[display3d]({p1,p2},title=``A graph of r^2 over its domain``);`,
`   `,
`SEE ALSO: plot3d[options], plots[display3d], drdtplot, dtdrplot`):

`help/text/trgraphplot`:=`help/text/rtgraphplot`:

`help/text/dxdydzplot`:=TEXT(`HELP FOR: dxdydzplot, dxdzdyplot, dydxdzplot, dydzdxplot, dzdxdyplot,`,
`	dzdydxplot`,
`   `,
`CALLING SEQUENCE:`,
`   `,
`	dxdydzplot(x=f(y,z)..g(y,z), y=h(z)..k(z), z=a..b, opts)`,
`	dxdzdyplot(x=f(y,z)..g(y,z), z=h(y)..k(y), y=a..b, opts)`,
`	dydxdzplot(y=f(x,z)..g(x,z), x=h(z)..k(z), z=a..b, opts)`,
`	dydzdxplot(y=f(x,z)..g(x,z), z=h(x)..k(x), x=a..b, opts)`,
`	dzdxdyplot(z=f(x,y)..g(x,y), x=h(y)..k(y), y=a..b, opts)`,
`	dzdydxplot(z=f(x,y)..g(x,y), y=h(x)..k(x), x=a..b, opts)`,
`   `,
`PARAMETERS:`,
`	The variable names ``x', ``y', and ``z' must be used,`,
`   `,
`	f, g are two-variable expressions,`,
`	h, k are one-variable expressions,`,
`   `,
`	a, b must be symbolic constants, such as Pi and E, integers,`,
`	or real floating point numbers; i.e., evalf(a), evalf(b) must`,
`	be of type numeric.`,
`	`,
`	opts is a list of plot3d options.`,
`   `,
`SYNOPSIS:   `,
`   `,
`-All six procedures plot regions in space with the given `,
` arguments providing ranges for the variables. The syntax is the same as if the`,
` int procedure were invoked for a triple integral over the region in cartesian`,
` coordinates. All the procedures return plot3d structures, accept plot3d`,
` options, and may be passed to the display3d procedure from the plots package.`,
`   `,
`EXAMPLES:   `,
`   `,
`dzdydxplot(z=0..1-x-y,y=0..1-x,x=0..1);`,
`   `,
`dzdxdyplot(z=0..1-x-y,x=0..1-y,y=0..1);`,
`   `,
`dzdydxplot(z=x^2+y^2..2,y=0..x,x=0..1, style=PATCH, orientation=[30,70],`,
`axes=NORMAL, scaling=CONSTRAINED, title=``An error in the text revealed``);`,
`   `,
`SEE ALSO: plot3d[options], plots[display3d], drdtdzplot, dpdtdphiplot`):

`help/text/dxdzdyplot`:=`help/text/dxdydzplot`:

`help/text/dydxdzplot`:=`help/text/dxdydzplot`:

`help/text/dydzdxplot`:=`help/text/dxdydzplot`:

`help/text/dzdxdyplot`:=`help/text/dxdydzplot`:

`help/text/dzdydxplot`:=`help/text/dxdydzplot`:

`help/text/drdtdzplot`:=TEXT(`HELP FOR: drdtdzplot, drdzdtplot, dtdrdzplot, dtdzdrplot, dzdrdtplot,`,
`	dzdtdrplot`,
`   `,
`CALLING SEQUENCE:`,
`   `,
`  drdtdzplot(r=f(theta,z)..g(theta,z), theta=h(z)..k(z), z=a..b, opts)`,
`  drdzdtplot(r=f(theta,z)..g(theta,z), z=h(theta)..k(theta), theta=a..b,opts)`,
`  dtdrdzplot(theta=f(r,z)..g(r,z), r=h(z)..k(z), z=a..b, opts)`,
`  dtdzdrplot(theta=f(r,z)..g(r,z), z=h(r)..k(r), r=a..b, opts)`,
`  dzdrdtplot(z=f(r,theta)..g(r,theta), r=h(theta)..k(theta), theta=a..b, opts)`,
`  dzdtdrplot(z=f(r,theta)..g(r,theta), theta=h(r)..k(r), r=a..b, opts)`,
`   `,
`PARAMETERS:`,
`	The variable names ``r', ``theta', and ``z' must be used,`,
`   `,
`	f, g are two-variable expressions,`,
`	h, k are one-variable expressions,`,
`   `,
`	a, b must be symbolic constants, such as Pi and E, integers,`,
`	or real floating point numbers; i.e., evalf(a), evalf(b) must`,
`	be of type numeric.`,
`	`,
`	opts is a sequence of plot3d options.`,
`   `,
`SYNOPSIS:   `,
`   `,
`-All six procedures plot regions in space with the given `,
` arguments providing ranges for the variables in cylindrical coordinates r, `,
` theta, z. The syntax is the same as if the int procedure were invoked for a`,
` triple integral over the region in cylindrical coordinates. All the procedures`,
` return plot3d structures, accept plot3d options, and may be passed to the`,
` display3d procedure from the plots package.`,
`   `,
`EXAMPLES:   `,
`   `,
`drdtdzplot(r = 2..3, theta = Pi/4..Pi/2,z = 0..1);`,
`   `,
`dzdtdrplot(z=-(1-r^2*cos(theta)^2)..(1-r^2*cos(theta)^2),theta=0..2*Pi,r=0..\
1);`,
`   `,
`dzdrdtplot(z=0..sqrt(4-r^2),r=0..1+cos(theta),theta=0..2*Pi);`,
`   `,
`p1:=dzdrdtplot(z=0..sqrt(1-r^2),r=0..1,theta=-Pi/2..Pi/2):`,
`p2:=dzdrdtplot(z=0..sqrt(1-r^2),r=0..1+cos(theta),theta=Pi/2..3*Pi/2):`,
`plots[display3d]({p1,p2});`,
`   `,
`SEE ALSO: plot3d[options], plots[display3d], dxdydzplot, dpdtdphiplot`):

`help/text/drdzdtplot`:=`help/text/drdtdzplot`:

`help/text/dtdrdzplot`:=`help/text/drdtdzplot`:

`help/text/dtdzdrplot`:=`help/text/drdtdzplot`:

`help/text/dzdrdtplot`:=`help/text/drdtdzplot`:

`help/text/dzdtdrplot`:=`help/text/drdtdzplot`:

`help/text/dpdtdphiplot`:=TEXT(`HELP FOR: dpdtdphiplot, dpdphidtplot, dtdpdphiplot, dtdphidpplot, dphidpdtplo\
t,`,
`	dphidtdpplot`,
`   `,
`CALLING SEQUENCE:`,
`   `, 
` dpdtdphiplot(rho=f(theta,phi)..g(theta,phi), theta=h(phi)..k(phi), z=a..b, opts)`,
`   `,
` dpdphidtplot(rho=f(theta,phi)..g(theta,phi), phi=h(theta)..k(theta), the\
ta=a..b, opts)`,
`   `,
` dtdpdphiplot(theta=f(p,phi)..g(p,phi), rho=h(phi)..k(phi), phi=a..b, opts)`,
`   `,
` dtdphidpplot(theta=f(rho,phi)..g(rho,phi), phi=h(rho)..k(rho), rho=a..b, opts)`,
`   `,
` dphidpdtplot(phi=f(rho,theta)..g(rho,theta), rho=h(theta)..k(theta), theta=a..b, opts)`,
`   `,
` dphidtdpplot(phi=f(rho,theta)..g(rho,theta), theta=h(rho)..k(rho), rho=a..b, opts)`,
`   `,
`PARAMETERS:`,
`	The variable names ``rho', ``theta', and ``phi' must be used,`,
`   `,
`	f, g are two-variable expressions,`,
`	h, k are one-variable expressions,`,
`   `,
`	a, b must be symbolic constants, such as Pi and E, integers,`,
`	or real floating point numbers; i.e., evalf(a), evalf(b) must`,
`	be of type numeric.`,
`	`,
`	opts is a sequence of plot3d options.`,
`   `,
`SYNOPSIS:   `,
`   `,
`-All six procedures plot regions in space with the given `,
` arguments providing ranges for the variables in spherical coordinates rho, `,
` theta, phi. The syntax is the same as if the int procedure were invoked for a `,
` triple integral over the region in spherical coordinates. All the procedures`
,` return plot3d structures, accept plot3d options, and may be passed to the`,
` display3d procedure from the plots package.`,
`   `,
`EXAMPLES:   `,
`   `,
`dphidpdtplot(phi=0..Pi/3,rho=0..1,theta=0..2*Pi);`,
`   `,
`P1:=dpdtdphiplot(rho=sec(phi)..2/(cos(phi)+sin(phi)),theta=0..2*Pi,`,
`phi=0..Pi/4):`,
`P2:=dpdtdphiplot(rho=0..cot(phi)*csc(phi),theta=0..2*Pi,phi=Pi/4..Pi/2):`,
`plots[display3d]({P1,P2},style=wireframe);`,
`   `,
`SEE ALSO: plot3d[options], plots[display3d], dxdydzplot, drdtdzplot`):

`help/text/dpdphidtplot`:=`help/text/dpdtdphiplot`:

`help/text/dtdpdphiplot`:=`help/text/dpdtdphiplot`:

`help/text/dtdphidpplot`:=`help/text/dpdtdphiplot`:

`help/text/dphidpdtplot`:=`help/text/dpdtdphiplot`:

`help/text/dphidtdpplot`:=`help/text/dpdtdphiplot`:

`help/text/regionplot3d`:=TEXT(`HELP FOR: regionplot3d`,
`	`,
`CALLING SEQUENCE:`,
`   `, 
` regionplot3d(w=h(u,v)..k(u,v),v=p(u)..q(u),u=a..b,[x(u,v,w),y(u,v,w),z(u,v,w)], opts)`,
`PARAMETERS:`,
`	`,
`	p(u), q(u) are expressions in the variable u,`,
`	h(u,v), k(u,v) are expressions in the variables u,v,`,
`	[x(u,v,w),y(u,v,w),z(u,v,w)] is a list of expressions`,
`	in the variables u,v,w,`,
`	`,
`	a, b must be symbolic constants, such as Pi and E, integers,`,
`	or real floating point numbers; i.e., evalf(a), evalf(b) must`,
`	be of type numeric.`,
`	`,
`	opts is a sequence of plot3d options.`,
`	`,
`SYNOPSIS:   `,
`	`,
`-The regionplot3d procedure plots the image of a three-dimensional region in`,
` uvw-space of the type h(u,v)<= w <=k(u,v), p(u)<= v <=q(u), a<= u <= b,`,
` under a coordinate transformation into xyz-space given by the list `,
` [x(u,v,w),y(u,v,w),z(u,v,w)]. The variable order is the same as if the int`,
` procedure were invoked for a triple integral over the region. The procedure`,
` returns a plot3d structure, accepts plot3d options, and may be passed to the`,
` display3d procedure from the plots package.`,
`	`,
`EXAMPLES: `,
`	`,
`regionplot3d(r=1..sqrt(10-z^2),t=0..2*Pi,z=-3..3,[r*cos(t),r*sin(t),z],title=``Cored Apple``);`,
`	`,
`regionplot3d(r=1..2,t=Pi/2..Pi,p=Pi/2..Pi,[r*cos(t)*sin(p),r*sin(t)*sin(p),`,
`r*cos(p)],title=``Spherical Coordinates``);`,
`	`,
`regionplot3d(w=0..1,v=0..1,u=0..1,[u*v-u*v*w,u*v*w,1-u],style=patch,axes=normal);`,
`	`,
`SEE ALSO: plot3d[options], plots[display3d]`):

#save `calcplot.m`;
#quit