function plot_segmentonneumann(md,options_structure,width,i,datai);
%PLOT_SEGMENTONNEUMANN - plot segment on neumann BC
%
%   Usage:
%      plot_segmentonneumann(md,options_structure,width,i);
%
%   See also: PLOTMODEL

%plot mesh boundaries
subplot(width,width,i); 

%units
if ~isnan(options_structure.unitmultiplier),
	md.x=md.x*options_structure.unitmultiplier;
	md.y=md.y*options_structure.unitmultiplier;
	md.z=md.z*options_structure.unitmultiplier;
end

if strcmp(datai,'segmentonneumann_diag') segmentonneumann=md.segmentonneumann_diag;end
if strcmp(datai,'segmentonneumann_prog') segmentonneumann=md.segmentonneumann_prog;end

%we are dealing with the loads. 
length_icefront=sqrt( (md.x(segmentonneumann(:,1))-md.x(segmentonneumann(:,2))).^2 + (md.y(segmentonneumann(:,1))-md.y(segmentonneumann(:,2))).^2 );
normal_icefront(:,1)=cos(atan2( (md.x(segmentonneumann(:,1))-md.x(segmentonneumann(:,2))) , (md.y(segmentonneumann(:,2))-md.y(segmentonneumann(:,1))) ) );
normal_icefront(:,2)=sin(atan2( (md.x(segmentonneumann(:,1))-md.x(segmentonneumann(:,2))) , (md.y(segmentonneumann(:,2))-md.y(segmentonneumann(:,1))) ) );

if strcmpi(md.type,'2d'),
	for n=1:size(segmentonneumann,1),

		plot(md.x(segmentonneumann(n,1:2)),md.y(segmentonneumann(n,1:2)),'b-','linewidth',2);
		hold on
		plot(md.x(segmentonneumann(n,1:2)),md.y(segmentonneumann(n,1:2)),'r.','markersize',12);

		xstart=mean(md.x(segmentonneumann(n,1:2)));
		ystart=mean(md.y(segmentonneumann(n,1:2)));

		xend=xstart+length_icefront(n)*normal_icefront(n,1);
		yend=ystart+length_icefront(n)*normal_icefront(n,2);

		plot([xstart xend],[ystart yend],'r-')
		legend('element edges on ice front','grids on ice front','normal vectors')
	end
else
	for n=1:size(segmentonneumann,1),
		hold on

		%Plot the two triangles of each quad (element edges on ice front)
		A=segmentonneumann(n,1);  B=segmentonneumann(n,2); C=segmentonneumann(n,3); D=segmentonneumann(n,4); E=md.numberofgrids+1;
		coordinates= [md.x md.y md.z; mean(md.x(segmentonneumann(n,1:4))) mean(md.y(segmentonneumann(n,1:4))) mean(md.z(segmentonneumann(n,1:4)))];
		p1=patch( 'Faces', [A B E A], 'Vertices',coordinates,'FaceVertexCData', [1 1 1],'FaceColor','flat','EdgeColor','blue','linewidth',2);
		patch( 'Faces', [A D E A], 'Vertices',coordinates,'FaceVertexCData', [1 1 1],'FaceColor','flat','EdgeColor','blue','linewidth',2);
		patch( 'Faces', [D C E D], 'Vertices',coordinates,'FaceVertexCData', [1 1 1],'FaceColor','flat','EdgeColor','blue','linewidth',2);
		patch( 'Faces', [C B E C], 'Vertices',coordinates,'FaceVertexCData', [1 1 1],'FaceColor','flat','EdgeColor','blue','linewidth',2);

		%Plot the nodes of each quad (grids on ice front)
		p2=plot3(md.x(segmentonneumann(n,1:4)),md.y(segmentonneumann(n,1:4)),md.z(segmentonneumann(n,1:4)),'r.','markersize',12);
		plot3(coordinates(end,1),coordinates(end,2),coordinates(end,3),'g.','markersize',12);

		%Retrieve coordinates of the four nodes of the Quad
		xyz1=[md.x(segmentonneumann(n,1));md.y(segmentonneumann(n,1));md.z(segmentonneumann(n,1));];
		xyz2=[md.x(segmentonneumann(n,2));md.y(segmentonneumann(n,2));md.z(segmentonneumann(n,2));];
		xyz3=[md.x(segmentonneumann(n,3));md.y(segmentonneumann(n,3));md.z(segmentonneumann(n,3));];
		xyz4=[md.x(segmentonneumann(n,4));md.y(segmentonneumann(n,4));md.z(segmentonneumann(n,4));];
		xyz5=coordinates(end,:)';
		
		%Build the norms of each triangle
		V1=cross(xyz5-xyz1,xyz5-xyz2); normal1=1/norm(V1)*V1';
		V2=cross(xyz5-xyz2,xyz5-xyz3); normal2=1/norm(V2)*V2';
		V3=cross(xyz5-xyz3,xyz5-xyz4); normal3=1/norm(V3)*V3';
		V4=cross(xyz5-xyz4,xyz5-xyz1); normal4=1/norm(V4)*V4';

		xstart1=mean([xyz1(1) xyz2(1) xyz5(1)]); ystart1=mean([xyz1(2) xyz2(2) xyz5(2)]); zstart1=mean([xyz1(3) xyz2(3) xyz5(3)]);
		xstart2=mean([xyz2(1) xyz3(1) xyz5(1)]); ystart2=mean([xyz2(2) xyz3(2) xyz5(2)]); zstart2=mean([xyz2(3) xyz3(3) xyz5(3)]);
		xstart3=mean([xyz3(1) xyz4(1) xyz5(1)]); ystart3=mean([xyz3(2) xyz4(2) xyz5(2)]); zstart3=mean([xyz3(3) xyz4(3) xyz5(3)]);
		xstart4=mean([xyz4(1) xyz1(1) xyz5(1)]); ystart4=mean([xyz4(2) xyz1(2) xyz5(2)]); zstart4=mean([xyz4(3) xyz1(3) xyz5(3)]);

		xend1=xstart1+length_icefront(n)*normal1(1);	yend1=ystart1+length_icefront(n)*normal1(2);	zend1=zstart1+length_icefront(n)*normal1(3);
		xend2=xstart2+length_icefront(n)*normal2(1);	yend2=ystart2+length_icefront(n)*normal2(2);	zend2=zstart2+length_icefront(n)*normal2(3);
		xend3=xstart3+length_icefront(n)*normal3(1);	yend3=ystart3+length_icefront(n)*normal3(2);	zend3=zstart3+length_icefront(n)*normal3(3);
		xend4=xstart4+length_icefront(n)*normal4(1);	yend4=ystart4+length_icefront(n)*normal4(2);	zend4=zstart4+length_icefront(n)*normal4(3);
		
		%plot the normals of each triangle (normal vectors)
		p3=plot3([xstart1 xend1],[ystart1 yend1],[zstart1 zend1],'r-');
		plot3([xstart2 xend2],[ystart2 yend2],[zstart2 zend2],'r-');
		plot3([xstart3 xend3],[ystart3 yend3],[zstart3 zend3],'r-');
		plot3([xstart4 xend4],[ystart4 yend4],[zstart4 zend4],'r-');

	end
	legend([p1,p2,p3],'element edges on ice front','grids on ice front','normal vectors');
end

%apply options
if isnan(options_structure.title)
	options_structure.title='Neumann boundary conditions';
end 
if isnan(options_structure.colorbar)
	options_structure.colorbar=0;
end
applyoptions(md,[],options_structure);
