function [Mvalue Mx My]=CreateDataMatrix(x_m,y_m,track_coord,track_values),
%CREATEDATAMATRIX - Create a map with average values of map 
%
%   This routine creates a map with average values of tracks.
%   x_m1 and y_m1 are two vectors containing the coordinates of the matrix
%   trac_coord is an exp file containing the coordinates of the tracks (x and y)
%   trav_values is a vector with the values along the track coordinates
%
%   Usage:
%      [Mvalue Mx My]=CreateDataMatrix(x_m,y_m,track_coord,track_values),
%
%   Example:
%      [Mvalue Mx My]=CreateDataMatrix(x_m,y_m,'trackcoord.exp',thickness_track)
%
%   See also: CREATEDATABOUNDARIES, TRACKSTOMATRIX

%Read the points of the tracks
stru=expread(track_coord,1);
nods=stru.nods;
xtracks=stru.x';
ytracks=stru.y';

%First check that the parameters are ok:
if (size(track_values,1)~=nods)  || (size(xtracks,2)~=nods) || (size(ytracks,2)~=nods),
	error('CreateDataMatrix error message : track coordinates and track values must have the same size');
end

%Compute number of rows and columns
numrow=size(y_m,1)-1;
numcol=size(x_m,1)-1;

%Remove useless points of the track
points=find(track_values==0);
track_values(points)=[];
xtracks(points)=[];
ytracks(points)=[];
points=find(isnan(track_values));
track_values(points)=[];
xtracks(points)=[];
ytracks(points)=[];
points=find(track_values<0);
track_values(points)=[];
xtracks(points)=[];
ytracks(points)=[];

points=find(xtracks<x_m(1) | xtracks>x_m(end) | ytracks<y_m(1) | ytracks>y_m(end));
track_values(points)=[];
xtracks(points)=[];
ytracks(points)=[];

%initialize some matrices
numpoints=zeros(numrow,numcol);
value=zeros(numrow,numcol);
coordx=zeros(numrow,numcol);
coordy=zeros(numrow,numcol);

%Loop over the points of the track
nel=size(track_values,1);
fprintf('%s','      track processing progress:   0.00 %');
for i=1:nel;
	if mod(i,1000)==0,
		fprintf('\b\b\b\b\b\b\b')
		fprintf('%5.2f%s',i/nel*100,' %');
	end

	x=xtracks(i);
	y=ytracks(i);

	%get indices for the matrix
	indexx=max(find(x_m<x));
	indexy=max(find(y_m<y));

	%get weighing coefficient
	val=track_values(i);

	%update numoverlap and weights
	numpoints(indexy,indexx)=numpoints(indexy,indexx)+1;
	value(indexy,indexx)=value(indexy,indexx)+val;
	coordx(indexy,indexx)=coordx(indexy,indexx)+x;
	coordy(indexy,indexx)=coordy(indexy,indexx)+y;

end
if nel>1000,
	fprintf('\b\b\b\b\b\b\b\b')
	fprintf('%4.2f%s\n',100,' %');
end

%Change the values of numoverlap to 1 if 0 since we are going to devide by this matrix
numpoints(find(~numpoints))=1;

%Create the center of mass for coordiantes and values.
Mvalue=value./numpoints;
Mx=coordx./numpoints;
My=coordy./numpoints;