#calcplot package version 1.0 by Tim Murdoch. Programmed: July/August, 1994.
#The package consists of twenty-eight procedures. Two procedures are for
#plotting surfaces of rotation, four 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 eighteen procedures are for plotting non-rectangluar
#regions (e.g., domains for triple integrals in cartesain, cylindrical,
#andspherical coordinates) in space.

#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: 8/16/94

#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.

rotyplot:=proc()

local a, b, c, d, f, i, xbound, rotaxis, f1, 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:=args[3];
a:=op(1, rhs(xbound));
b:=op(2, rhs(xbound));

if not type (xbound, string = range) then
ERROR(`range expression for x is incorrect`) fi;
if not type (rotaxis, string = 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:=op(1, rhs(rotaxis));
d:=lhs(xbound);

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,style=PATCH, scaling=CONSTRAINED,
axes=BOXED, grid=[25,25],labels=[` `,x,y]); fi;

if nargs > 3 then
opt_seq:=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, 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:=args[3];
a:=op(1, rhs(xbound));
b:=op(2, rhs(xbound));

if not type (xbound, string = range) then
ERROR(`range expression for x is incorrect`) fi;
if not type (rotaxis, string = 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:=op(1, rhs(rotaxis));
d:=lhs(xbound);

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,style=PATCH, scaling=CONSTRAINED,
axes=BOXED, grid=[25,25], labels=[` `,x,y]); fi;

if nargs > 3 then
opt_seq:=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, 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));

if not type (xbound, string = range) then
ERROR(`range expression for x is incorrect`) fi;
if not type (ybound, string = 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);

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,axes=NORMAL,scaling=CONSTRAINED); fi;

if nargs > 2 then
opt_seq:=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, 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));

if not type (xbound, string = range) then
ERROR(`range expression for x is incorrect`) fi;
if not type (ybound, string = 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);

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,axes=NORMAL,scaling=CONSTRAINED); fi;

if nargs >2 then
opt_seq:=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, 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));

if not type (rbound, string = range) then
ERROR(`range expression for r is incorrect`) fi;
if not type (tbound, string = 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);

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 positions of the angle functions are not correct`)  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, axes=NORMAL, scaling=CONSTRAINED); fi;

if nargs > 2 then
opt_seq:=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, 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));

if not type (rbound, string = range) then
ERROR(`range expression for r is incorrect`) fi;
if not type (tbound, string = 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);

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(`the radius must always be positive`)  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, axes=NORMAL, scaling=CONSTRAINED); fi;

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

plots[display](p1);

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

if not type (xbound, string = range) then
ERROR(`range expression for x is incorrect`) fi;
if not type (ybound, string = 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);

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,style=PATCH,axes=BOXED,scaling=CONSTRAINED, 
grid=[20,20], labels=[x,y,z]); fi;

if nargs > 3 then
opt_seq:=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, 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));

if not type (xbound, string = range) then
ERROR(`range expression for x is incorrect`) fi;
if not type (ybound, string = 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);

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,style=PATCH,axes=BOXED,scaling=CONSTRAINED,
grid=[20,20], labels=[x,y,z]); fi;

if nargs > 3 then
opt_seq:=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, 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));

if not type (tbound, string = range) then
ERROR(`range expression for theta is incorrect`) fi;
if not type (rbound, string = 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);

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(`the radius must always be positive`)  fi; od;

if nargs = 3 then
p1:=plots[cylinderplot](surface,u=0..1,v=0..1,style=PATCH,axes=BOXED,
scaling=CONSTRAINED, grid=[20,20], labels=[x,y,z]); 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:

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

trgraphplot:=proc()

local a, b, c, d, g, k, f1, g1, h1, i, h, 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));

if not type (tbound, string = range) then
ERROR(`range expression for theta is incorrect`) fi;
if not type (rbound, string = 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);

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,style=PATCH,axes=BOXED,
scaling=CONSTRAINED, grid=[10,40], labels=[x,y,z]); 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:

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

if not type (zbound, string = 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]};

if nargs = 3 then
p1:=plot3d(surface, u=0..1, v=0..1, axes=BOXED, grid =[20,20],
scaling=CONSTRAINED, style=PATCH);
fi;

if nargs > 3 then
opt_seq:=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, 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));

if not type (zbound, string = 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]};

if nargs = 3 then
p1:=plot3d(surface, u=0..1, v=0..1, axes=BOXED, grid =[20,20],
scaling=CONSTRAINED, style=PATCH);
fi;

if nargs > 3 then
opt_seq:=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, 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));

if not type (zbound, string = 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,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),uloc,vloc],[g1(vloc,uloc),uloc,vloc],
[wloc,uloc,h1(uloc)], [xloc,uloc,k1(uloc)],[yloc,a,uloc1],[zloc,b,uloc2]};

if nargs = 3 then
p1:=plot3d(surface, u=0..1, v=0..1, axes=BOXED, grid =[20,20],
scaling=CONSTRAINED, style=PATCH);
fi;

if nargs > 3 then
opt_seq:=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, 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));

if not type (zbound, string = 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]};

if nargs = 3 then
p1:=plot3d(surface, u=0..1, v=0..1, axes=BOXED, grid =[20,20],
scaling=CONSTRAINED, style=PATCH);
fi;

if nargs > 3 then
opt_seq:=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, 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));

if not type (zbound, string = 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]};

if nargs = 3 then
p1:=plot3d(surface, u=0..1, v=0..1, axes=BOXED, grid =[20,20],
scaling=CONSTRAINED, style=PATCH);
fi;

if nargs > 3 then
opt_seq:=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, 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));

if not type (zbound, string = 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]};

if nargs = 3 then
p1:=plot3d(surface, u=0..1, v=0..1, axes=BOXED, grid =[20,20],
scaling=CONSTRAINED, style=PATCH);
fi;

if nargs > 3 then
opt_seq:=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, 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));

if not type (zbound, string = 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]};

if nargs = 3 then
p1:=plots[cylinderplot](surface, u=0..1, v=0..1, axes=BOXED, grid
=[20,20],scaling=CONSTRAINED,style=PATCH); 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, 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));

if not type (rbound, string = 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]};

if nargs = 3 then
p1:=plots[cylinderplot](surface, u=0..1, v=0..1, axes=BOXED, grid
=[20,20],scaling=CONSTRAINED,style=PATCH); 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 :

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

if not type (rbound, string = 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]};

if nargs = 3 then
p1:=plots[cylinderplot](surface, u=0..1, v=0..1, axes=BOXED, grid
=[20,20],scaling=CONSTRAINED,style=PATCH); 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 :

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

if not type (zbound, string = 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]};

if nargs = 3 then
p1:=plots[cylinderplot](surface, u=0..1, v=0..1, axes=BOXED, grid
=[20,20],scaling=CONSTRAINED,style=PATCH); 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 :

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

if not type (zbound, string = 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]};

if nargs = 3 then
p1:=plots[cylinderplot](surface, u=0..1, v=0..1, axes=BOXED, grid
=[20,20],scaling=CONSTRAINED,style=PATCH); 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 :


#6. d(theta)drdz

dtdrdzplot:= proc()

local surface, f, f1, g, g1, h, h1, k, k1, a, b, r1, t1, z1, u, v,
uloc, uloc1, uloc2, vloc, wloc, xloc, yloc, zloc, i, 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));

if not type (zbound, string = 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]};

if nargs = 3 then
plots[cylinderplot](surface, u=0..1, v=0..1, axes=BOXED, grid
=[20,20],scaling=CONSTRAINED,style=PATCH); fi;

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

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

if not type (phibound, string = 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]};

if nargs =3 then
q1:=plots[sphereplot](surface, u=0..1, v=0..1, axes=BOXED, grid
=[15,25],style=PATCH, scaling=CONSTRAINED); fi;

if nargs > 3 then
opt_seq:=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, 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));

if not type (phibound, string = 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]};

if nargs =3 then
q1:=plots[sphereplot](surface, u=0..1, v=0..1, axes=BOXED, grid
=[15,25],style=PATCH, scaling=CONSTRAINED); fi;

if nargs > 3 then
opt_seq:=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, 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));

if not type (phibound, string = 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]};

if nargs =3 then
q1:=plots[sphereplot](surface, u=0..1, v=0..1, axes=BOXED, grid
=[15,25],style=PATCH, scaling=CONSTRAINED); fi;

if nargs > 3 then
opt_seq:=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, 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));

if not type (phibound, string = 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]};

if nargs =3 then
q1:=plots[sphereplot](surface, u=0..1, v=0..1, axes=BOXED, grid
=[15,25],style=PATCH, scaling=CONSTRAINED); fi;

if nargs > 3 then
opt_seq:=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, 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));

if not type (phibound, string = 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]};

if nargs =3 then
q1:=plots[sphereplot](surface, u=0..1, v=0..1, axes=BOXED, grid
=[15,25],style=PATCH, scaling=CONSTRAINED); fi;

if nargs > 3 then
opt_seq:=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, 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));

if not type (phibound, string = 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]};

if nargs =3 then
q1:=plots[sphereplot](surface, u=0..1, v=0..1, axes=BOXED, grid
=[15,25],style=PATCH, scaling=CONSTRAINED); fi;

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

plots[display3d](q1);

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`,
`   `,
`-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.`,
`   `,
`SEE ALSO: plot3d[options], plots[display3d]`,
`	`,
`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``);`):

`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.`,
`   `,
`SEE ALSO: plot[options], plots[display], xygraphplot, yxgraphplot`,
`	`,
`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``);`):

`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.`,
`	`,
`SEE ALSO: plot[options], plots[dislay], rtgraphplot, trgraphplot`,
`   `,
`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)``);`):

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

`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.`,`   `,
`SEE ALSO: plot3d[options], plots[display3d], dxdyplot, dydxplot`,
`   `,
`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``);`):

`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 procedures `,
` draw four line segments that pass through the interior of the region. 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.`,
`   `,
`SEE ALSO: plot3d[options], plots[display3d], drdtplot, dtdrplot`,
`   `,
`EXAMPLES:   `,
`   `,
`rtgraphplot(r^2*sin(theta)*cos(theta),r=1..1+cos(theta),theta=-Pi/2..Pi/2);`,
`   `,
`rgraphplot(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``);`):

`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.`,
`   `,
`SEE ALSO: plot3d[options], plots[display3d], drdtdzplot, dpdtdphiplot`,
`   `,
`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``);`):

`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.`,
`   `,
`SEE ALSO: plot3d[options], plots[display3d], dxdydzplot, dpdtdphiplot`,
`   `,
`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});`):

`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.`,
`   `,
`SEE ALSO: plot3d[options], plots[display3d], dxdydzplot, drdtdzplot`,
`   `,
`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);`):

`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`:

#save calcplot.m
#quit