% Here is the Heapsort algorithm we will need later
	vardef heapsort(suffix d,I) =
		save l,ir,exitin,exitout,rra,k,i,j;
		l:=floor(d/2)+1;
		ir:=d;
		boolean exitin,exitout;
		exitout:=false;
		forever: exitif exitout;
			if(l>1):
				l:=l-1;
 				rra:=d[l];
	 			k:=I[l];
			else:
				rra:=d[ir];
				k:=I[ir]; 
				d[ir]:=d[1];
				I[ir]:=I[1];
				ir:=ir-1;
				if(ir<2):
					d[1]:=rra;
					I[1]:=k;
					exitout:=true;
				fi;
			fi;
			if not exitout:
				i:=l;
				j:=l+1;
				exitin:=false;
				forever: exitif exitin;
					if ((j<ir) and (d[j]<d[j+1])): j:=j+1; fi;
						if (rra<d[j]):
							d[i]:=d[j];
							I[i]:=I[j];
							i:=j;
							j:=2*j;
					else: exitin:=true;
					fi;
					if (j>=ir): exitin:=true; fi;
				endfor;
				d[i]:=rra;
				I[i]:=k;
			fi;
		endfor;
	enddef;

	% theta and phi are for the 3D to 2D projection
	theta:=-37.5;
	phi:=20;

	% Function to surface plot
	vardef f(expr x,y) = sin(pi*x)*cos(pi*y) enddef;

	% x & y range (and number of grid points)
	ax:=0; bx:=2; Nx:=40;
	ay:=0; by:=2; Ny:=40;

	% Approximate width of graphic and x,y,z scale factors;
	w:=3.5in;
	rx:=1; ry:=1; rz:=2/(1+sqrt(5));

	% x & y spacing
	dx:=(bx-ax)/(Nx-1);
	dy:=(by-ay)/(Ny-1);

	% Evaluate f(x,y) at grid points and compute range of f(x,y)
	pair Z[][];
	numeric F[][];
	az:=infinity; bz:=-infinity;
	for nx=0 upto Nx-1:
		for ny=0 upto Ny-1:
			Z[nx][ny]:=(ax+nx*dx,ay+ny*dy);
			F[nx][ny]:=f(xpart Z[nx][ny],ypart Z[nx][ny]);
			if (F[nx][ny] < az):
				az:=F[nx][ny];
			fi;
			if (F[nx][ny] > bz):
				bz:=F[nx][ny];
			fi;
		endfor;
	endfor;

	% Scale factor for x,y,z and define viewpoint
	ux:=rx*w/(bx-ax);
	uy:=ry*w/(by-ay);
	uz:=rz*w/(bz-az);
	color viewpoint;
	viewpoint:=(sind(theta)*cosd(phi),-cosd(theta)*cosd(phi),sind(phi));

	% T is the orthographic projection function
	% P is the same as T but scaled for drawing
	def T primary c = ((redpart c)*cosd(theta)+(greenpart c)*sind(theta),((greenpart c)*cosd(theta)-(redpart c)*sind(theta))*sind(phi)+(bluepart c)*cosd(phi)) enddef;
	def P primary c = T(ux*(redpart c),uy*(greenpart c),uz*(bluepart c)) enddef;

	% Calculate the distance d from each point to the viewpoint plane
	% Redefine Z to be the projected (x,y,f(x,y)) points
	numeric d[][];
	for ny=0 upto Ny-1:
		for nx=0 upto Nx-1:
			d[nx][ny]:=((xpart Z[nx][ny])*(redpart viewpoint))+((ypart Z[nx][ny])*(greenpart viewpoint))+((F[nx][ny])*(bluepart viewpoint));
			Z[nx][ny]:=P(xpart Z[nx][ny],ypart Z[nx][ny],F[nx][ny]);
		endfor;
	endfor;

	% Define the distance from each "element" to the viewpoint plane
	% by using the closest corner of the element
	numeric D,D[],I,I[];
	D:=(Nx-1)*(Ny-1); I:=D;
	for ny=0 upto Ny-2:
		for nx=0 upto Nx-2:
			m:=ny*(Nx-1)+nx+1;
			D[m]:=max(d[nx][ny],d[nx+1][ny],d[nx+1][ny+1],d[nx][ny+1]);
			I[m]:=m;
		endfor;
	endfor;

	% Sort the distances so that the furthest elements are filled/drawn first
	heapsort(D,I);

	% Here's a hardcoded axes set just for this example
	for x=0 step 0.5 until 2:
		draw P(x,0,-1)--P(x,2,-1)--P(x,2,1) dashed evenly scaled 0.5;
	endfor;
	for y=0 step 0.5 until 2:
		draw P(0,y,-1)--P(2,y,-1)--P(2,y,1) dashed evenly scaled 0.5;
	endfor;
	for z=-1 step 0.5 until 1:
		draw P(2,0,z)--P(2,2,z)--P(0,2,z) dashed evenly scaled 0.5;
	endfor;
	draw P(0,2,-1)--P(0,2,1);
	draw P(0,0,-1)--P(0,2,-1);
	draw P(0,0,-1)--P(2,0,-1);

	% Fill and draw the surface grid
	path rgJet,bJet;
	rgJet:=(0,0)--(0,0)--(0,1)--(1,1)--(1,0)--(1/2,0);
	bJet:=(1/2,0)--(1,0)--(1,0)--(0,0)--(0,0)--(0,0);
	JetLen:=length rgJet;
	path elempath;
	color c;
	for m=1 upto I:
		ny:=floor((I[m]-1)/(Nx-1));
		nx:=I[m]-1-ny*(Nx-1);
		elempath:=Z[nx][ny]--Z[nx+1][ny]--Z[nx+1][ny+1]--Z[nx][ny+1]--cycle;
		Fm:=(F[nx][ny]+F[nx+1][ny]+F[nx+1][ny+1]+F[nx][ny+1])/4;
		n:=(Fm-az)/(bz-az)*JetLen;
		c:=(xpart point n of rgJet,ypart point n of rgJet,xpart point n of bJet);
		fill elempath withcolor c;
		draw elempath;
	endfor;