function md=planetmixedmesh(md,radius,angleresol)
%PLANETMIXEDMESH: build 2d shell mesh
%
% Usage: md=planetmixedshell(md,radius,angleresol)
%

conv=pi/180;

r=radius;
theta=(0:angleresol:360)';
theta=theta*conv;
phi=(0:angleresol:180)';
phi=phi(2:end-1)*conv;

nphi=length(phi);
nthetha=length(theta);

x=zeros(nphi*nthetha,1);
y=zeros(nphi*nthetha,1);
z=zeros(nphi*nthetha,1);

for i=1:nphi,
	phii=phi(i);
	x((i-1)*nthetha+1:i*nthetha)=r.*cos(theta).*sin(phii);
	y((i-1)*nthetha+1:i*nthetha)=r.*sin(theta).*sin(phii);
	z((i-1)*nthetha+1:i*nthetha)=r.*cos(phii);
end

quads=zeros((nthetha-1)*(nphi-1),4);

%build elements
for i=1:nphi-1,
	for j=1:nthetha-1,
		count=(i-1)*(nthetha-1)+j;

		A=(i-1)*nthetha+j;
		B=(i-1)*nthetha+j+1;
		C=(i)*nthetha+j+1;
		D=(i)*nthetha+j;

		quads(count,:)=[A B C D];
	end
end

%now, add polar elements:
%first north: phi = 0;
x=[r.*cos(0).*sin(0);x];
y=[r.*sin(0).*sin(0);y];
z=[r.*cos(0);z];

%add 1 to quads: 
quads=quads+1;

%add tria elements:
trias=zeros(nthetha-1,4);

for i=1:nthetha-1,
	A=1;
	B=i+1;
	C=i+2;
	trias(i,:)=[A B C NaN];
end

quads=[trias;quads];

% now add south pole: 
phii=180*conv;
x=[x;r.*cos(phii).*sin(phii)];
y=[y;r.*sin(phii).*sin(phii)];
z=[z;r.*cos(phii)];
nods=length(x);

%add tria elements:
trias=zeros(nthetha-1,4);

start=nods-nthetha;
for i=1:nthetha-1,
	A=start+i-1;
	B=start+i;
	C=nods;
	trias(i,:)=[A B C NaN];
end

quads=[quads;trias];

md.mesh.elements=quads;
md.mesh.x=x;
md.mesh.y=y;
md.mesh.z=z;
md.mesh.r=sqrt(x.^2+y.^2+z.^2);
md.mesh.theta=acos(z./md.mesh.r);
md.mesh.phi=atan2(y,x);

md.mesh.numberofvertices=length(md.mesh.x);
md.mesh.numberofelements=size(md.mesh.elements,1);

md.mesh.dimension=3;
