function [index,X,Y,Z,S,data_interp]=SectionValues(md,data,infile,resolution)
%SECTIONVALUES - compute the value of a field on a section
%
%   This routine gets the value of a given field of the model on points
%   given by filname (Argus type file)
%
%   Usage:
%      [elements,x,y,z,s,data]=SectionValues(md,data,filename,resolution)

%read infile:
contempt=expread(infile,1);
nods=contempt.nods;
x=contempt.x;
y=contempt.y;

%get the specified resolution
if isnumeric(resolution(1))
	res_h=resolution(1);
else
	error('SectionValues error message: wrong resolution type. Resolution must be an array [horizontal_resolution vertical_resolution]')
end
if strcmpi(md.type,'3d')
	if (length(resolution)==2 & isnumeric(resolution(2)))
		res_v=resolution(2);
	else
		error('SectionValues error message: wrong resolution type. Resolution must be an array [horizontal_resolution vertical_resolution]')
	end
end

%initialization
X=[]; %X-coordinate
Y=[]; %Y-coordinate
S=0;  %curvilinear coordinate

for i=1:nods-1

	x_start=x(i);
	x_end=x(i+1);
	y_start=y(i);
	y_end=y(i+1);
	s_start=S(end);

	length_segment=sqrt((x_end-x_start)^2+(y_end-y_start)^2);
	portion=ceil(length_segment/res_h);

	x_segment=zeros(portion,1);
	y_segment=zeros(portion,1);
	s_segment=zeros(portion,1);

	for j=1:portion
		x_segment(j)=x_start+(j-1)*(x_end-x_start)/portion;
		y_segment(j)=y_start+(j-1)*(y_end-y_start)/portion;
		s_segment(j)=s_start+j*length_segment/portion;
	end

	%plug into X and Y
	X=[X;x_segment];
	Y=[Y;y_segment];
	S=[S;s_segment];
end
X(end+1)=x(nods);
Y(end+1)=y(nods);

%Number of grids:
numberofgrids=size(X,1);

%Compute Z
Z=zeros(numberofgrids,1);

%New mesh and Data interpolation
if strcmpi(md.type,'2d')

	%Interpolation of data on specified points
	data_interp=griddata_mesh_to_mesh(md.elements,md.x,md.y,data,X,Y);

	%Compute index
	index=[1:1:(numberofgrids-1);2:1:numberofgrids]';

else

	%vertically extrude mesh

	%Get bed and surface for each 2d point, offset to make sure that it is inside the glacier system
	offset=10^-10;
	bed=griddata_mesh_to_mesh_3d(md.elements,md.x,md.y,md.z,md.bed,X,Y,Z,'node')+offset;
	surface=griddata_mesh_to_mesh_3d(md.elements,md.x,md.y,md.z,md.surface,X,Y,Z,'node')-offset;

	%Some useful parameters
	layers=ceil(res_v/md.numlayers);
	gridsperlayer=numberofgrids;
	gridstot=gridsperlayer*layers;
	elementsperlayer=gridsperlayer-1;
	elementstot=(gridsperlayer-1)*(layers-1);

	%initialization
	X3=zeros(gridsperlayer*layers,1); Y3=zeros(gridsperlayer*layers,1); Z3=zeros(gridsperlayer*layers,1); S3=zeros(gridsperlayer*layers,1); index3=zeros(elementstot,4);

	%Get new coordinates in 3d
	for i=1:layers
		X3(i:layers:end)=X;
		Y3(i:layers:end)=Y;
		Z3(i:layers:end)=bed+(i-1)*(surface-bed)/(layers-1);
		S3(i:layers:end)=S;

		if i<layers %Build index3 with quads
			index3((i-1)*elementsperlayer+1:i*elementsperlayer,:)=[i:layers:gridstot-layers; i+1:layers:gridstot-layers; i+layers+1:layers:gridstot; i+layers:layers:gridstot]';
		end
	end

	%Interpolation of data on specified points
	data_interp=griddata_mesh_to_mesh_3d(md.elements,md.x,md.y,md.z,data,X3,Y3,Z3);

	%build outputs
	X=X3; Y=Y3; Z=Z3;  S=S3; index=index3;
end
