# 3D Plot Region Package Version 1.0

with(plots):

# plot3d_region(x= a..b,y = f1(x) .. f2(x),z=g1(x,y) .. g2(x,y)) or
# plot3d_region(x= a..b,y = f1(x) .. f2(x), num_bars)


# P1 := plot2d_region(x=.1..1,y=x^3..x^2,8):

#display(P1);

plot3d_region := proc(u_range,v_range,w_range)
	local i,j,var1,a,b,ran1,ran2,var2,var3,ran3,v_lower,v_upper,w_lower,w_upper,i,num_bars,t_lower,t_upper,v_lower_fcn,v_upper_fcn,w_lower_fcn,w_upper_fcn,gr1,gr2,gr3,nu,nv,nw,u_val,v_val,w_val:
	num_bars := 6:
	if nargs = 4 then num_bars := args[4] fi:

	var1 := op(1,u_range):
	var2 := op(1,v_range):
	var3 := op(1,w_range):

	ran1 := op(2,u_range):
	ran2 := op(2,v_range):
	ran3 := op(2,w_range):

	a := evalf(op(1,ran1)):
	b := evalf(op(2,ran1)):

	v_lower := op(1,ran2):
	v_upper := op(2,ran2):

	v_lower_fcn := unapply(v_lower,var1):
	v_upper_fcn := unapply(v_upper,var1):


	w_lower := op(1,ran3):
	w_upper := op(2,ran3):

	w_lower_fcn := unapply(w_lower,var1,var2):
	w_upper_fcn := unapply(w_upper,var1,var2):

	nu := num_bars-1:
	nv := num_bars-1:
	nw := num_bars-1:

	u_val := array(1..nu+1):
	v_val := array(1..nu+1,1..nv+1):

	for i from 1 to nu+1 do
	    for j from 1 to nv+1 do
		u_val[i]:=evalf(((i-1)/nu)*a+(1 - (i-1)/nu)*b):
		v_val[i,j]:=evalf(((j-1)/nv)*v_lower_fcn(u_val[i])+
			  (1 - (j-1)/nv)*v_upper_fcn(u_val[i])):
	    od:
	od:

	gr1 := seq(seq([[u_val[i],
		         v_val[i,j],
		         w_lower_fcn(u_val[i],v_val[i,j])],
		        [u_val[i],
		         v_val[i,j+1],
		         w_lower_fcn(u_val[i],v_val[i,j+1])],
		        [u_val[i+1],
		         v_val[i+1,j+1],
		         w_lower_fcn(u_val[i+1],v_val[i+1,j+1])],
		        [u_val[i+1],
		         v_val[i+1,j],
		         w_lower_fcn(u_val[i+1],v_val[i+1,j])]
		        ],
		     i=1..nu),
                j=1..nv),
                seq(seq([[u_val[i],
		         v_val[i,j],
		         w_upper_fcn(u_val[i],v_val[i,j])],
		        [u_val[i],
		         v_val[i,j+1],
		         w_upper_fcn(u_val[i],v_val[i,j+1])],
		        [u_val[i+1],
		         v_val[i+1,j+1],
		         w_upper_fcn(u_val[i+1],v_val[i+1,j+1])],
		        [u_val[i+1],
		         v_val[i+1,j],
		         w_upper_fcn(u_val[i+1],v_val[i+1,j])]
		        ],
		     i=1..nu),
               j=1..nv):

	u_val := array(1..nu+1):
	#v_val(i,1) uses v_lower.
	#v_val(i,2) uses v_upper.
	v_val := array(1..nu+1,1..2):
	# w_val(i,1,1) uses v_lower and w_lower.
	# w_val(i,2,1) uses v_upper and w_lower.
	# w_val(i,1,2) uses v_lower and w_upper.
	# w_val(i,2,2) uses v_upper and w_upper.
	w_val := array(1..nu+1,1..2,1..2):

	for i from 1 to nu+1 do
		u_val[i]:=evalf(((i-1)/nu)*a+(1 - (i-1)/nu)*b):
		v_val[i,1]:= evalf(v_lower_fcn(u_val[i])):
		v_val[i,2]:= evalf(v_upper_fcn(u_val[i])):
		w_val[i,1,1]:=evalf(w_lower_fcn(u_val[i],v_val[i,1])):
		w_val[i,1,2]:=evalf(w_upper_fcn(u_val[i],v_val[i,1])):
		w_val[i,2,1]:=evalf(w_lower_fcn(u_val[i],v_val[i,2])):
		w_val[i,2,2]:=evalf(w_upper_fcn(u_val[i],v_val[i,2])):
	od:

	gr2 := seq(seq([[u_val[i],
		         v_val[i,1],
		         ((j-1)/nw)*w_val[i,1,1]+
			  (1-(j-1)/nw)*w_val[i,1,2]],
			[u_val[i],
		         v_val[i,1],
		         (j/nw)*w_val[i,1,1]+
			  (1-j/nw)*w_val[i,1,2]],
			[u_val[i+1],
		         v_val[i+1,1],
		         (j/nw)*w_val[i+1,1,1]+
			  (1-j/nw)*w_val[i+1,1,2]],
			[u_val[i+1],
		         v_val[i+1,1],
		         ((j-1)/nw)*w_val[i+1,1,1]+
			  (1-(j-1)/nw)*w_val[i+1,1,2]]
			],i=1..nu),
               j=1..nw),
	       seq(seq([[u_val[i],
		         v_val[i,2],
		         ((j-1)/nw)*w_val[i,2,1]+
			  (1-(j-1)/nw)*w_val[i,2,2]],
			[u_val[i],
		         v_val[i,2],
		         (j/nw)*w_val[i,2,1]+
			  (1-j/nw)*w_val[i,2,2]],
			[u_val[i+1],
		         v_val[i+1,2],
		         (j/nw)*w_val[i+1,2,1]+
			  (1-j/nw)*w_val[i+1,2,2]],
			[u_val[i+1],
		         v_val[i+1,2],
		         ((j-1)/nw)*w_val[i+1,2,1]+
			  (1-(j-1)/nw)*w_val[i+1,2,2]]
			],i=1..nu),
               j=1..nw):



	#v_val(1,i) uses a.
	#v_val(2,i) uses b.
	v_val := array(1..2,1..nv+1):
	# w_val(1,i,1) uses a and w_lower.
	# w_val(2,i,1) uses b and w_lower.
	# w_val(1,i,2) uses a and w_upper.
	# w_val(2,i,2) uses b and w_upper.
	w_val := array(1..2,1..nv+1,1..2):

	for i from 1 to nv+1 do
		v_val[1,i]:= evalf(((i-1)/nv)*v_lower_fcn(a)+
			     (1- (i-1)/nv)*v_upper_fcn(a)):
		v_val[2,i]:= evalf(((i-1)/nv)*v_lower_fcn(b)+
			     (1- (i-1)/nv)*v_upper_fcn(b)):
		w_val[1,i,1]:=evalf(w_lower_fcn(a,v_val[1,i])):
		w_val[2,i,1]:=evalf(w_lower_fcn(b,v_val[2,i])):
		w_val[1,i,2]:=evalf(w_upper_fcn(a,v_val[1,i])):
		w_val[2,i,2]:=evalf(w_upper_fcn(b,v_val[2,i])):

	od:

	gr3 := seq(seq([[a,
			 v_val[1,i],
			 ((j-1)/nw)*w_val[1,i,1]+
			  (1-(j-1)/nw)*w_val[1,i,2]],
			[a,
			 v_val[1,i+1],
			 ((j-1)/nw)*w_val[1,i+1,1]+
			  (1-(j-1)/nw)*w_val[1,i+1,2]],
			[a,
			 v_val[1,i+1],
			 (j/nw)*w_val[1,i+1,1]+
			  (1-j/nw)*w_val[1,i+1,2]],
			[a,
			 v_val[1,i],
			 (j/nw)*w_val[1,i,1]+
			  (1-j/nw)*w_val[1,i,2]]
			],
			i=1..nv),
		j=1..nw),
                seq(seq([[b,
			 v_val[2,i],
			 ((j-1)/nw)*w_val[2,i,1]+
			  (1-(j-1)/nw)*w_val[2,i,2]],
			[b,
			 v_val[2,i+1],
			 ((j-1)/nw)*w_val[2,i+1,1]+
			  (1-(j-1)/nw)*w_val[2,i+1,2]],
			[b,
			 v_val[2,i+1],
			 (j/nw)*w_val[2,i+1,1]+
			  (1-j/nw)*w_val[2,i+1,2]],
			[b,
			 v_val[2,i],
			 (j/nw)*w_val[2,i,1]+
			  (1-j/nw)*w_val[2,i,2]]
			],
			i=1..nv),
		j=1..nw):
	
[polygonplot3d([gr1],axes=boxed,shading=zhue,style=patch,labels=[var1,var2,var3]),polygonplot3d([gr2],axes=boxed,shading=zhue,style=patch,labels=[var1,var2,var3]),polygonplot3d([gr3],axes=boxed,shading=zhue,style=patch,labels=[var1,var2,var3])]:
	#{seq(polygonplot3d([gr.i],axes=boxed,shading=zhue,style=patch,labels=[var1,var2,var3]),i=1..3)}:

end:

plot3d_top_bot := proc(u_range,v_range,w_range)
	local i,j,var1,a,b,ran1,ran2,var2,var3,ran3,v_lower,v_upper,w_lower,w_upper,i,num_bars,t_lower,t_upper,v_lower_fcn,v_upper_fcn,w_lower_fcn,w_upper_fcn,gr1,gr2,gr3,nu,nv,nw,u_val,v_val,w_val:
	num_bars := 6:
	if nargs = 4 then num_bars := args[4] fi:

	var1 := op(1,u_range):
	var2 := op(1,v_range):
	var3 := op(1,w_range):

	ran1 := op(2,u_range):
	ran2 := op(2,v_range):
	ran3 := op(2,w_range):

	a := evalf(op(1,ran1)):
	b := evalf(op(2,ran1)):

	v_lower := op(1,ran2):
	v_upper := op(2,ran2):

	v_lower_fcn := unapply(v_lower,var1):
	v_upper_fcn := unapply(v_upper,var1):


	w_lower := op(1,ran3):
	w_upper := op(2,ran3):

	w_lower_fcn := unapply(w_lower,var1,var2):
	w_upper_fcn := unapply(w_upper,var1,var2):

	nu := num_bars-1:
	nv := num_bars-1:
	nw := num_bars-1:

	u_val := array(1..nu+1):
	v_val := array(1..nu+1,1..nv+1):

	for i from 1 to nu+1 do
	    for j from 1 to nv+1 do
		u_val[i]:=evalf(((i-1)/nu)*a+(1 - (i-1)/nu)*b):
		v_val[i,j]:=evalf(((j-1)/nv)*v_lower_fcn(u_val[i])+
			  (1 - (j-1)/nv)*v_upper_fcn(u_val[i])):
	    od:
	od:

	gr1 := seq(seq([[u_val[i],
		         v_val[i,j],
		         w_lower_fcn(u_val[i],v_val[i,j])],
		        [u_val[i],
		         v_val[i,j+1],
		         w_lower_fcn(u_val[i],v_val[i,j+1])],
		        [u_val[i+1],
		         v_val[i+1,j+1],
		         w_lower_fcn(u_val[i+1],v_val[i+1,j+1])],
		        [u_val[i+1],
		         v_val[i+1,j],
		         w_lower_fcn(u_val[i+1],v_val[i+1,j])]
		        ],
		     i=1..nu),
                j=1..nv),
                seq(seq([[u_val[i],
		         v_val[i,j],
		         w_upper_fcn(u_val[i],v_val[i,j])],
		        [u_val[i],
		         v_val[i,j+1],
		         w_upper_fcn(u_val[i],v_val[i,j+1])],
		        [u_val[i+1],
		         v_val[i+1,j+1],
		         w_upper_fcn(u_val[i+1],v_val[i+1,j+1])],
		        [u_val[i+1],
		         v_val[i+1,j],
		         w_upper_fcn(u_val[i+1],v_val[i+1,j])]
		        ],
		     i=1..nu),
               j=1..nv):

	u_val := array(1..nu+1):
	#v_val(i,1) uses v_lower.
	#v_val(i,2) uses v_upper.
	v_val := array(1..nu+1,1..2):
	# w_val(i,1,1) uses v_lower and w_lower.
	# w_val(i,2,1) uses v_upper and w_lower.
	# w_val(i,1,2) uses v_lower and w_upper.
	# w_val(i,2,2) uses v_upper and w_upper.
	w_val := array(1..nu+1,1..2,1..2):

	for i from 1 to nu+1 do
		u_val[i]:=evalf(((i-1)/nu)*a+(1 - (i-1)/nu)*b):
		v_val[i,1]:= evalf(v_lower_fcn(u_val[i])):
		v_val[i,2]:= evalf(v_upper_fcn(u_val[i])):
		w_val[i,1,1]:=evalf(w_lower_fcn(u_val[i],v_val[i,1])):
		w_val[i,1,2]:=evalf(w_upper_fcn(u_val[i],v_val[i,1])):
		w_val[i,2,1]:=evalf(w_lower_fcn(u_val[i],v_val[i,2])):
		w_val[i,2,2]:=evalf(w_upper_fcn(u_val[i],v_val[i,2])):
	od:

	polygonplot3d([gr1],axes=boxed,shading=zhue,style=patch,labels=[var1,var2,var3]):

end:

 help_plot3d_region := proc();
	printf(`\n\t\t\t\tRegion Plot Package Version 1.0\n`);
	printf(`\nThis set of functions aids in picturing elementary 3 dimensional regions of integration.\n`);
	printf(`\n`);
	printf(`Basic examples include:\n`);
	printf(`\t\tP1 := plot3d_region(x = -1 .. 1, y = x .. x^2+5*x,\n`);
	printf(`\t\tz = 0 .. x^2+y^2);\n`);
	printf(`and\n`);
	printf(`\t\tP1 := plot3d_region(x = -1 .. 1, y = x .. x^2+5*x,\n`);
	printf(`\t\tz = 0 .. x^2+y^2, 8);\n`);
	printf(`to specify an 8 by 8 grid on each side of the region instead of the default 5 by 5.\n`);
	printf(`\n`);
	printf(`A list of three plot3d structures is returned.\n \n`);
	printf(`All six sides may be viewed by, for example:\n`);
	printf(`\t\tdisplay(P1,orientation=[45,80]);\n \n`);
	printf(`Just the top and bottom may be seen by:\n`);
	printf(`\t\tdisplay(P1[1],orientation=[45,80]);\n \n`);
	printf(`Other pairs of sides may be viewed by:\n`);
	printf(`\t\tdisplay(P1[2],orientation=[45,80]);\n`);
	printf(`or\n`);
	printf(`\t\tdisplay(P1[3],orientation=[45,80]);\n \n`);
	printf(`Another routine is:\n`);
	printf(`\t\thelp_plot3d_region();\n`);
	printf(`to display this message.\n`);
	printf(`\n`);


end:

help_plot3d_region();