comment
Nicer 3D plots
endcomment

function f3d (ffunction,n=10)
## f3d("f") or f3d("f",n;...)
## Draw a function defined on [-1,1]^2 with x-, y- and z-axis.
## 2n is the spacing.
	x=-1:1/n:1; y=x;
	{X,Y}=field(x,y);
	if isfunction(ffunction); Z=ffunction(X,Y,args());
	else; Z=expreval(ffunction,X,y=Y);
	endif;
	view(6,3,0.4,0.2);
	z1=totalmin(Z); z2=totalmax(Z); h=(z2+z1)/2; Z=Z-h;
	i=n+1:2*n+1; j=1:n+1;
	hold off; solid(X[i,j],Y[i,j],Z[i,j]); hold on;
	linewidth(3); 
	wire([0,0],[0,0],[-1.5,1.5]);
	wire([-1.5,0],[0,0],[-h,-h]);
	wire([0,0],[0,1.5],[-h,-h]);
	linewidth(1);
	solid(X[i,i],Y[i,i],Z[i,i]);
	solid(X[j,j],Y[j,j],Z[j,j]);
	linewidth(3);
	wire([0,1.5],[0,0],[-h,-h]);
	wire([0,0],[-1.5,0],[-h,-h]);
	linewidth(1);
	solid(X[j,i],Y[j,i],Z[j,i]); hold off;
	return ""
endfunction

function fcd (ffunction,n=10,xmin=-1,xmax=1,ymin=-1,ymax=1,nc=20)
## plots a function of two variables or expression of x and y
## with density and contour.
	x=linspace(xmin,xmax,n); y=linspace(ymin,ymax,n);
	{X,Y}=field(x,y);
	if isfunction(ffunction); Z=ffunction(X,Y,args());
	else; Z=expreval(ffunction,X,y=Y);
	endif;
	z1=totalmin(Z); z2=totalmax(Z);
	density(Z,1);
	h=holding(1); contour(Z,linspace(z1,z2,nc)); holding(h);
	setplot(xmin,xmax,ymin,ymax);
	return ""
endfunction

function fcontour (ffunction,n=10,xmin=-1,xmax=1,ymin=-1,ymax=1,nc=20)
## Draw contour lines of a function or expression in x and y.
## 2n is the spacing.
	x=linspace(xmin,xmax,n); y=linspace(ymin,ymax,n);
	{X,Y}=field(x,y);
	if isfunction(ffunction); Z=ffunction(X,Y,args());
	else; Z=expreval(ffunction,X,y=Y);
	endif;
	z1=totalmin(Z); z2=totalmax(Z);
	contour(Z,linspace(z1,z2,nc));
	setplot(xmin,xmax,ymin,ymax);
	return ""
endfunction

function f3dpolar (ffunction,n=10)
## f3dpolar("f") or f3dpolar("f",n;...)
## Draw a function defined on [-1,1]^2 with x-, y- and z-axis.
## 2n is the spacing.
	r=0:1/n:1; p=linspace(0,2*pi,8*n);
	{P,R}=field(p,r); X=R*cos(P); Y=R*sin(P);
	if isfunction(ffunction); Z=ffunction(X,Y,args());
	else; Z=expreval(ffunction,X,y=Y);
	endif;
	view(6,3,0.4,0.2);
	z1=totalmin(Z); z2=totalmax(Z); h=(z2+z1)/2; Z=Z-h;
	i=1:n+1;
	j=2*n:4*n+1;
	solid(X[i,j],Y[i,j],Z[i,j]); hold on;
	linewidth(3);
	wire([0,0],[0,0],[-1.5,1.5]);
	wire([-1.5,0],[0,0],[-h,-h]);
	wire([0,0],[0,1.5],[-h,-h]);
	linewidth(1);
	j=1:2*n+1;
	solid(X[i,j],Y[i,j],Z[i,j]);
	j=4*n:6*n+1;
	solid(X[i,j],Y[i,j],Z[i,j]);
	linewidth(3);
	wire([0,1.5],[0,0],[-h,-h]);
	wire([0,0],[-1.5,0],[-h,-h]);
	linewidth(1);
	j=6*n:8*n+1;
	solid(X[i,j],Y[i,j],Z[i,j]); hold off;
	return ""
endfunction

function upperwindow (title="")
## select the upper window for plot
	window(150,860,75,500);
	title(title);
	return title
endfunction

function lowerwindow (title="")
## select the upper window for plot
	window(150,860,575,1010);
	ctext(title,512,510);
	return title
endfunction

.. eof