function plot_rifpenetration(md,options,nlines,ncols,index);
%PLOT_RIFTPENETRATION - plot rift penetration
%
%   Usage:
%      plot_rifpenetration(md,options,width,i);
%
%   See also: PLOTMODEL

%process data and model
[x y z elements is2d isplanet]=processmesh(md,[],options);

subplot(nlines,ncols,index); 
hold on

%plot mesh boundaries
for i=1:size(md.mesh.segments,1),
	plot(x(md.mesh.segments(i,1:2)),y(md.mesh.segments(i,1:2)),'k-');
end

isp1=0;
isp2=0;

if isstruct(md.rifts.riftstruct),
	%plot mesh boundaries
	for i=1:size(md.mesh.segments,1),
		h1=plot(x(md.mesh.segments(i,1:2)),y(md.mesh.segments(i,1:2)),'b-');
	end
	for i=1:size(md.rifts.riftstruct,1),
		penaltypairs=md.rifts.riftstruct(i).penaltypairs;

		segments=md.rifts.riftstruct(i).segments;
		for j=1:size(segments,1),
			plot(x(segments(j,1:2)),y(segments(j,1:2)),'b-');
		end

		normal=zeros(2,1);
		for j=1:size(penaltypairs,1),
			normal(1)=penaltypairs(j,5);
			normal(2)=penaltypairs(j,6);

			vx1=md.initialization.vx(penaltypairs(j,1)); 
			vx2=md.initialization.vx(penaltypairs(j,2));
			vy1=md.initialization.vy(penaltypairs(j,1)); 
			vy2=md.initialization.vy(penaltypairs(j,2));
			penetration=(vx2-vx1)*normal(1)+(vy2-vy1)*normal(2);
			%if penetration is negative, plot in black, positive, plot in red;: ie: if rift is closing, black, if rift is opening, red.
			if(penetration>0),
				p2=plot(x(penaltypairs(j,1)) ,y(penaltypairs(j,1)),'ro-','LineWidth',1);
				set(p2,'MarkerSize',3);
				isp2=1;
			else
				p1=plot(x(penaltypairs(j,1)) ,y(penaltypairs(j,1)),'ko-','LineWidth',1);
				set(p1,'MarkerSize',3);
				isp1=1;
			end
		end

		%point out the tips
		h2=plot(x(md.rifts.riftstruct(i).tips(1)),y(md.rifts.riftstruct(i).tips(1)),'g*');
		plot(x(md.rifts.riftstruct(i).tips(2)),y(md.rifts.riftstruct(i).tips(2)),'g*');
	end
	if strcmpi(getfieldvalue(options,'legend','on'),'on'),
		if isp1 & isp2
			l=legend([h1,h2,p1,p2],'mesh boundaries','crack tips','faults','rifts');
		elseif isp1
			l=legend([h1,h2,p1],'mesh boundaries','crack tips','faults');
		elseif isp2
			l=legend([h1,h2,p2],'mesh boundaries','crack tips','rifts');
		else
			l=legend([h1,h2],'mesh boundaries','crack tips');
		end
		set(l,'Location',getfieldvalue(options,'legend_location','NorthEast'));
	end
else
	error('plot error message: no rifts available!');
end
hold off

%apply options
options=addfielddefault(options,'title','Rift/Fault location');
options=addfielddefault(options,'colorbar',0);
applyoptions(md,[],options);