function md=meshyams(md,varargin);
%MESHYAMS - Build model of Antarctica by refining according to observed velocity error estimator
%
%   Usage:
%      md=meshyams(md,varargin);
%      where varargin is a lit of paired arguments. 
%      arguments can be: 'domainoutline': Argus file containing the outline of the domain to be meshed
%      arguments can be: 'velocities': matlab file containing the velocities [m/yr]
%      optional arguments: 'groundeddomain': Argus file containing the outline of the grounded ice
%                          this option is used to minimize the metric on water (no refinement)
%      optional arguments: 'resolution': initial mesh resolution [m]
%      optional arguments: 'nsteps': number of steps of mesh adaptation
%      optional arguments: 'epsilon': average interpolation error wished [m/yr]
%      optional arguments: 'hmin': minimum edge length
%      optional arguments: 'hmanx': maximum edge
%      optional arguments: 'riftoutline': if rifts are present, specifies rift outline file.
%      
%
%   Examples:
%      md=meshyams(md,'domainoutline','Domain.exp','velocities','vel.mat');
%      md=meshyams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp');
%      md=meshyams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp','nsteps',6,'epsilon',2,'hmin',500,'hmax',30000);

%recover options
options=pairoptions(varargin{:});
options=deleteduplicates(options,1);

%recover some fields
disp('MeshYams Options:')
domainoutline=getfieldvalue(options,'domainoutline');
disp(sprintf('   %-15s: ''%s''','DomainOutline',domainoutline));
riftoutline=getfieldvalue(options,'riftoutline','N/A');
disp(sprintf('   %-15s: ''%s''','riftoutline',riftoutline));
groundeddomain=getfieldvalue(options,'groundeddomain','N/A');
disp(sprintf('   %-15s: ''%s''','GroundedDomain',groundeddomain));
velocities=getfieldvalue(options,'velocities');
disp(sprintf('   %-15s: ''%s''','Velocities',velocities));
resolution=getfieldvalue(options,'resolution',5000);
disp(sprintf('   %-15s: %f','Resolution',resolution));
nsteps=getfieldvalue(options,'nsteps',6);
disp(sprintf('   %-15s: %i','nsteps',nsteps));
gradation=getfieldvalue(options,'gradation',2*ones(nsteps,1));
disp(sprintf('   %-15s: %g','gradation',gradation(1)));
epsilon=getfieldvalue(options,'epsilon',3);
disp(sprintf('   %-15s: %f','epsilon',epsilon));
hmin=getfieldvalue(options,'hmin',500);
disp(sprintf('   %-15s: %f','hmin',hmin));
hmax=getfieldvalue(options,'hmax',150*10^3);
disp(sprintf('   %-15s: %f\n','hmax',hmax));

%mesh with initial resolution
disp('Initial mesh generation...');
if strcmpi(riftoutline,'N/A');
	md=mesh(md,domainoutline,resolution);
else
	md=mesh(md,domainoutline,riftoutline,resolution);
	md=meshprocessrifts(md,domainoutline);
end
disp(['Initial mesh, number of elements: ' num2str(md.numberofelements)]);

%load velocities 
disp('loading velocities...');
Names=VelFindVarNames(velocities);
Vel=load(velocities);

%start mesh adaptation
for i=1:nsteps,
	disp(['Iteration #' num2str(i) '/' num2str(nsteps)]);

	%interpolate velocities onto mesh
	disp('   interpolating velocities...');
	if strcmpi(Names.interp,'grid'),
		vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.x,md.y,0);
		vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.x,md.y,0);
	else
		vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.x,md.y,0);
		vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.x,md.y,0);
	end
	field=sqrt(vx_obs.^2+vy_obs.^2);

	%set gridonwater field
	if ~strcmp(groundeddomain,'N/A'),
		gridground=ContourToMesh(md.elements,md.x,md.y,groundeddomain,'node',2);
		md.gridonwater=ones(md.numberofgrids,1);
		md.gridonwater(find(gridground))=0;
	else
		md.gridonwater=zeros(md.numberofgrids,1);
	end

	%adapt according to velocities
	disp('   adapting...');
	md=YamsCall(md,field,hmin,hmax,gradation(i),epsilon);

	%if we have rifts, we just messed them up, we need to recreate the segments that constitute those 
	%rifts, because the segments are used in YamsCall to freeze the rifts elements during refinement.
	if md.numrifts, 
		md.nodeconnectivity=NodeConnectivity(md.elements,md.numberofgrids);
		md.elementconnectivity=ElementConnectivity(md.elements,md.nodeconnectivity);
		md.segments=findsegments(md);
		md=meshyamsrecreateriftsegments(md);
	end

end
	
disp(['Final mesh, number of elements: ' num2str(md.numberofelements)]);

%Now, build the connectivity tables for this mesh.
md.nodeconnectivity=NodeConnectivity(md.elements,md.numberofgrids);
md.elementconnectivity=ElementConnectivity(md.elements,md.nodeconnectivity);

%recreate segments
md.segments=findsegments(md);
md.gridonboundary=zeros(md.numberofgrids,1); md.gridonboundary(md.segments(:,1:2))=1;

%Fill in rest of fields:
md.z=zeros(md.numberofgrids,1);
md.gridonbed=ones(md.numberofgrids,1);
md.gridonsurface=ones(md.numberofgrids,1);
md.elementonbed=ones(md.numberofelements,1);
md.elementonsurface=ones(md.numberofelements,1);
if ~strcmp(groundeddomain,'N/A'),
	gridground=ContourToMesh(md.elements,md.x,md.y,groundeddomain,'node',2);
	md.gridonwater=ones(md.numberofgrids,1);
	md.gridonwater(find(gridground))=0;
else
	md.gridonwater=zeros(md.numberofgrids,1);
end
if strcmpi(Names.interp,'grid'),
	md.vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.x,md.y,0);
	md.vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.x,md.y,0);
else
	md.vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.x,md.y,0);
	md.vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.x,md.y,0);
end
md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);

%deal with rifts 
if md.numrifts,
	%first, recreate rift segments
	md=meshyamsrecreateriftsegments(md);

	%using the segments, recreate the penaltypairs
	for j=1:md.numrifts,
		rift=md.rifts(j);

		%build normals and lengths of segments:
		lengths=sqrt((md.x(rift.segments(:,1))-md.x(rift.segments(:,2))).^2 + (md.y(rift.segments(:,1))-md.y(rift.segments(:,2))).^2 );
		normalsx=cos(atan2((md.x(rift.segments(:,1))-md.x(rift.segments(:,2))) , (md.y(rift.segments(:,2))-md.y(rift.segments(:,1)))));
		normalsy=sin(atan2((md.x(rift.segments(:,1))-md.x(rift.segments(:,2))) , (md.y(rift.segments(:,2))-md.y(rift.segments(:,1)))));

		%ok, build penaltypairs: 
		numpenaltypairs=length(rift.segments)/2-1;
		rift.penaltypairs=zeros(numpenaltypairs,7);

		for i=1:numpenaltypairs,
			rift.penaltypairs(i,1)=rift.segments(i,2);
			rift.penaltypairs(i,2)=rift.segments(end-i,2);
			rift.penaltypairs(i,3)=rift.segments(i,3);
			rift.penaltypairs(i,4)=rift.segments(end-i,3);
			rift.penaltypairs(i,5)=normalsx(i)+normalsx(i+1);
			rift.penaltypairs(i,6)=normalsy(i)+normalsy(i+1);
			rift.penaltypairs(i,7)=(lengths(i)+lengths(i+1))/2;
		end
		%renormalize norms: 
		norms=sqrt(rift.penaltypairs(:,5).^2+rift.penaltypairs(:,6).^2);
		rift.penaltypairs(:,5)=rift.penaltypairs(:,5)./norms;
		rift.penaltypairs(:,6)=rift.penaltypairs(:,6)./norms;

		md.rifts(j)=rift;
	end

end