function [x y z elements is2d isplanet]=processmesh(md,data,options)
%PROCESSMESH - process mesh to be plotted
%
%   Usage:
%      [x y z elements is2d]=processmesh(md,data,options)
%
%   See also: PLOTMODEL, PROCESSDATA

%some checks
if md.mesh.numberofvertices==0,
	error('plot error message: mesh is empty')
end
if md.mesh.numberofvertices==md.mesh.numberofelements
	error(['plot error message: the number of elements is the same as the number of nodes...']);
end

if (isempty(data) | ~isstruct(data)),
	%first load x,y, etc ... to speed up plot

	if ~strcmpi(getfieldvalue(options,'coord','xy'),'latlon'),
		x=md.mesh.x;
		if isprop(md.mesh,'x2d'), x2d=md.mesh.x2d; end
		y=md.mesh.y;
		if isprop(md.mesh,'y2d'), y2d=md.mesh.y2d; end
	else
		x=md.mesh.long;
		y=md.mesh.lat;
	end

	z_field=getfieldvalue(options,'z',md.mesh.z);
	if ischar(z_field),
		z=md.(z_field);
	elseif isnumeric(z_field),
		z=z_field;
	else
		z=md.mesh.z;
	end

	if isprop(md.mesh,'elements2d'), elements2d=md.mesh.elements2d; end
	elements=md.mesh.elements;

	%is it a 2d plot?
	if (md.mesh.dimension==2);
		is2d=1;
	elseif (md.mesh.dimension==3),
		if getfieldvalue(options,'layer',0)>=1,
			is2d=1;
		else
			is2d=0;
		end
	else
		error(['dim = ' num2str(md.mesh.dimension) ' not supported yet']);
	end

	%layer projection? 
	if getfieldvalue(options,'layer',0)>=1,
		if strcmpi(getfieldvalue(options,'coord','xy'),'latlon'),
			error('processmesh error message: cannot work with 3D meshes for now');
		end
		%we modify the mesh temporarily to a 2d mesh from which the 3d mesh was extruded. 
		x=x2d;
		y=y2d;
		z=zeros(size(x2d));
		elements=elements2d;
	end
else
	%Process Patch
	if (md.mesh.dimension==2),
		elements=transpose(reshape(1:3*md.mesh.numberofelements,3,md.mesh.numberofelements));
		x=transpose(reshape(md.mesh.x(data.index)',1,3*md.mesh.numberofelements));
		y=transpose(reshape(md.mesh.y(data.index)',1,3*md.mesh.numberofelements));
		z=zeros(3*md.mesh.numberofelements,1);
		is2d=1;
	else
		elements=transpose(reshape(1:6*md.mesh.numberofelements,6,md.mesh.numberofelements));
		x=transpose(reshape(md.mesh.x(data.index)',1,6*md.mesh.numberofelements));
		y=transpose(reshape(md.mesh.y(data.index)',1,6*md.mesh.numberofelements));
		z=transpose(reshape(md.mesh.z(data.index)',1,6*md.mesh.numberofelements));
		is2d=0;
	end
end

%units
if exist(options,'unit'),
	unit=getfieldvalue(options,'unit');
	x=x*unit;
	y=y*unit;
	z=z*unit;
end

if isa(md,'planet'),
	isplanet=1;
else
	isplanet=0;
end
