function [data datatype]=processdata(md,data,options);
%PROCESSDATA - process data to be plotted
%
%   datatype = 1 -> elements
%   datatype = 2 -> nodes
%   datatype = 3 -> node quivers
%   datatype = 4 -> patch
%
%   Usage:
%      [data datatype]=processdata(md,data,options);
%
%   See also: PLOTMODEL, PROCESSMESH

%check format
if (iscell(data) | isempty(data) | length(data)==0 | (length(data)==1 & ~isstruct(data) & isnan(data))),
	error('plotmodel error message: data provided is empty');
end

%Process Patch
if isstruct(data) 
	if (isfield(data,'index') & isfield(data,'value')),
		if data.interpolation(1)==P1Enum(),
			data=data.value;
			data=data';
			data=data(:);
			datatype=4;
		elseif data.interpolation(1)==P0Enum(),
			data=data.value;
			datatype=5;
		else
			error(['interpolation ' data.interpolation(1) ' not supported yet']);
		end
	else
		error('structure other than Patch not supported yet');
	end
else
	%initialize datatype
	datatype=0;
end

%get datatype
datasize=size(data);

%non patch processing
if datatype~=4 & datatype~=5,

	%transpose data if necessary
	if (size(data,2) > size(data,1)),
		data=data';
	end
	datasize=size(data);

	%convert to double if necessary
	if ~isnumeric(data);
		disp('processdata info message: data is not numeric (logical?). Converted to double');
		data=double(data);
	end

	%check length
	if datasize(1)~=md.numberofgrids & datasize(1)~=md.numberofelements & datasize(1)~=md.numberofgrids*6 & (md.dim==3 & ~(datasize(1)==md.numberofelements2d | datasize(1)==md.numberofgrids2d))
		error('plotmodel error message: data not supported yet');
	end

	%quiver?
	if datasize(2)>1,
		datatype=3;

		%check number of columns, add zeros if necessary,
		if (md.dim==3)
			if datasize(2)==2,
				data=[data, zeros(datasize(1),1)];
			elseif datasize(2)~=3,
				error('plotmodel error message: data provided should have 2 or 3 columns for quiver plot, and 1 for regular plot');
			end
		%elseif ((md.dim==2) & datasize(2)~=2),
		%	error('plotmodel error message: data provided should have 2 columns for quiver plot, and 1 for regular plot');
		end
	end

	%treat the case datasize(1)=6*grids
	if datasize(1)==6*md.numberofgrids
		%keep the only norm of data
		data1=data(1:6:md.numberofgrids*6,:);
		data2=data(2:6:md.numberofgrids*6,:);
		data=sqrt(data1.^2+data2.^2);
		datasize(1)=md.numberofgrids;
		%---> go to grid data
	end

	%treat the case datasize(1)=grids2d
	if (md.dim==3 & datasize(1)==md.numberofgrids2d),
		data=project3d(md,data,'node');
		datasize(1)=md.numberofgrids;
		%---> go to grid data
	end

	%treat the case datasize(1)=grids2d
	if (md.dim==3 & datasize(1)==md.numberofelements2d),
		data=project3d(md,data,'element');
		datasize(1)=md.numberofelements;
		%---> go to grid data
	end

	%smoothing?
	if exist(options,'smooth')
		data=averaging(md,data,getfieldvalue(options,'smooth'));
		datasize(1)=md.numberofgrids;
		%---> go to grid data
	end
end

%element data
if (datasize(1)==md.numberofelements & datasize(2)==1),

	%Initialize datatype if non patch
	if datatype~=4 & datatype~=5,
		datatype=1;
	end

	%ice sheet only?
	if getfieldvalue(options,'iceshelf',1)==0,
		data(find(md.elementoniceshelf),:)=NaN;
	end
	%ice shelf only?
	if getfieldvalue(options,'icesheet',1)==0,
		data(find(~md.elementoniceshelf),:)=NaN;
	end
	%no water?
	if getfieldvalue(options,'water',1)==0,
		data(find(md.elementonwater),:)=NaN;
	end
	%log?
	if exist(options,'log'),
		bounds=getfieldvalue(options,'caxis',[min(data(:)) max(data(:))]);
		data(find(data<bounds(1)))=bounds(1);
		if any(data<=0),
			error('Log option cannot be applied on negative values. Use caxis option (Rignot''s settings: [1.5 max(data)])');
		end
		pos=find(~isnan(data));
		data(pos)=log(data(pos))/log(getfieldvalue(options,'log'));
	end
end

%grid data
if (datasize(1)==md.numberofgrids & datasize(2)==1),
	datatype=2;
	%ice sheet only?
	if getfieldvalue(options,'iceshelf',1)==0,
		pos=find(md.gridoniceshelf);
		data(pos,:)=NaN;
	end
	%ice shelf only?
	if getfieldvalue(options,'icesheet',1)==0,
		pos=find(md.gridonicesheet);
		data(pos,:)=NaN;
	end
	%no water?
	if getfieldvalue(options,'water',1)==0,
		pos=find(md.gridonwater);
		data(pos,:)=NaN;
	end
	%log?
	if exist(options,'log'),
		%if any(data<=0),
		%	error('Log option cannot be applied on negative values. Use caxis option (Rignot''s settings: [1.5 max(data)])');
		%end
		data=log(data)/log(getfieldvalue(options,'log'));
	end

end

%layer projection? 
if getfieldvalue(options,'layer',0)>=1,
	data=project2d(md,data,getfieldvalue(options,'layer')); %project onto 2d mesh
end

%control arrow density if quiverplot
if datatype==3 & exist(options,'density')
	databak=data;
	data=NaN*ones(datasize);
	density=getfieldvalue(options,'density');
	data(1:density:end,:)=databak(1:density:end,:);
	clear databak
end