function md2=modelzoom(md,areaname)
%MODELZOOM - exctract a smaller model out of a large one. 
%
%   It is used to be able to plot large models, one area at a time. 
%   Be careful: usage is: md2=modelzoom(md,areaname) where md2 is the returned 
%   model, md the large input model, and areaname the name of the area wherefrom 
%   the smaller model should be extracted. 
%   Because this is intended for mainly plotting purposes, the private member segments will 
%   not be extracted. Ie: this model cannot be used for solving.
%
%   Usage:
%      md2=modelzoom(md,areaname)
%
%   Example:
%      md2=modelzoom(md,'pineislandglacier');
%
%   See also: MODELEXTRACT

%some checks on list of arguments
if ((nargin~=2) | (nargout~=1)),
	modelzoomusage();
	error('modelzoom error message');
end

[xlim0 ylim0]=antmodelzoom(areaname);

%Create a list of elements locations from the input model
xel=(md.x(md.elements(:,1))+md.x(md.elements(:,2))+md.x(md.elements(:,3)))/3;
yel=(md.y(md.elements(:,1))+md.y(md.elements(:,2))+md.y(md.elements(:,3)))/3;

%Figure out which elements are within the area: 
elementlist=find( (xel>xlim0(1)) & (xel<xlim0(2)) & (yel>ylim0(1)) & (yel<ylim0(2)) );

%Create a list of elements locations from the input model, which are on the segmentonneumann_diag: 
if(~isempty(md.segmentonneumann_diag)),
	xelneumann=(md.x(md.elements(md.segmentonneumann_diag(:,3),1))+md.x(md.elements(md.segmentonneumann_diag(:,3),2))+md.x(md.elements(md.segmentonneumann_diag(:,3),3)))/3;
	yelneumann=(md.y(md.elements(md.segmentonneumann_diag(:,3),1))+md.y(md.elements(md.segmentonneumann_diag(:,3),2))+md.y(md.elements(md.segmentonneumann_diag(:,3),3)))/3;
	%Figure out which elements are within the area, and on segmentonneumann_diag: 
	elementlistneumann=find( (xelneumann>xlim0(1)) & (xelneumann<xlim0(2)) & (yelneumann>ylim0(1)) & (yelneumann<ylim0(2)) );
end


%Figure out the corresponding grids in the area: 
flags=zeros(md.numberofgrids,1);
flags(md.elements(elementlist,1))=1;
flags(md.elements(elementlist,2))=1;
flags(md.elements(elementlist,3))=1;
gridlist=find(flags);

%Reverse indexing
gridpos=zeros(md.numberofgrids,1);
gridpos(gridlist)=1:1:length(gridlist);
elempos=zeros(md.numberofelements,1);
elempos(elementlist)=1:1:length(elementlist);

%Create output model: 
md2=model;
md2.numberofgrids=length(gridlist);
md2.numberofelements=length(elementlist);

structfields=fields(md);
for i=1:length(structfields),
	field=structfields(i);field=field{1};
	fieldval=getfield(md,field);

	if isnumeric(fieldval),
		if isempty(fieldval),
			md2=setfield(md2,field,fieldval);
		end
		if length(fieldval)==md.numberofgrids,
			fieldval=fieldval(gridlist);
			md2=setfield(md2,field,fieldval);
		end
		if length(fieldval)==md.numberofelements,
			if strcmp(field,'elements'),
				fieldval=fieldval(elementlist,:);
				md2=setfield(md2,field,fieldval);
				%Now, we need to reset the element grid ids to the 
				%new mesh grid ids.
				md2.elements(:)=gridpos(md2.elements(:));
			else
				fieldval=fieldval(elementlist,:);
				md2=setfield(md2,field,fieldval);
			end
		end
		if isscalar(fieldval),
			if (~(strcmp(field,'numberofelements') | strcmp(field,'numberofgrids')))
				md2=setfield(md2,field,fieldval);
			end
		end
	end
	if ischar(fieldval),
		md2=setfield(md2,field,fieldval);
	end
end

%Deal with segmentonneumann_diag: 
if(~isempty(md.segmentonneumann_diag)),
	md2.segmentonneumann_diag=md.segmentonneumann_diag(elementlistneumann,:);
	md2.segmentonneumann_diag(:,1)=gridpos(md2.segmentonneumann_diag(:,1));
	md2.segmentonneumann_diag(:,2)=gridpos(md2.segmentonneumann_diag(:,2));
	md2.segmentonneumann_diag(:,3)=elempos(md2.segmentonneumann_diag(:,3));
end

%Add a note: 
md2=addnote(md2,['This model is a modelzoom on area: ' areaname]);
end %end of function


function modelzoomusage(),
disp('md2=modelzoom(md,areaname');
end %end of function

