function index = get_point(geometry,x,y,tol)
%This function finds the index of a point at given coordinates with tolerance tol
%index = get_point(geometry,x,y,tol)

 %set default tolerance=================================================================================================
 if nargin < 4 
    tol = eps;
 end

 % plot geometry =======================================================================================================
 %geomplot(geometry,'pointlabels','on','edgelabels','on');

 %get number of objects=================================================================================================
 gd = geominfo(geometry,'out',{'gd'});
 no = geominfo(geometry,'out',{'no'},'od',0:gd);
 
 %get vertex coordinates================================================================================================
 xy=flgeomvtx(geometry);
 
 %get indice that lie in the region of interest=========================================================================

 % getting elements with x-coordinate > x-tol 
 element_x_g = ceil(find(xy(1,:) >=  x - tol));

 % getting elements with x-coordinate < x+tol 
 element_x_s = ceil(find(xy(1,:) <=  x + tol));

 % creating a matrix for the elements greater x
 coordinate_matrix_x_g = zeros(1,no(1));

 % creating a matrix for the elements smaller x
 coordinate_matrix_x_s = zeros(1,no(1));

 % filling the matrix with a one when element fullfills the condition > x
 coordinate_matrix_x_g(element_x_g) = 1;

 % filling the matrix with a one when element fullfills the condition < x
 coordinate_matrix_x_s(element_x_s) = 1;

 % Merge the x-coordinate_s and x-coordinate_g matrix
 coordinate_matrix_x = coordinate_matrix_x_g & coordinate_matrix_x_s;

 % getting elements with y-coordinate > y-tol 
 element_y_g = ceil(find(xy(2,:) >=  y - tol));
 
 % getting elements with y-coordinate < y + tol 
 element_y_s = ceil(find(xy(2,:) <=  y + tol));

 % creating a matrix for the elements greater y
 coordinate_matrix_y_g = zeros(1,no(1));

 % creating a matrix for the elements smaller y
 coordinate_matrix_y_s = zeros(1,no(1));

 % filling the matrix with an one when element fullfills the condition > y
 coordinate_matrix_y_g(element_y_g) = 1;

 % filling the matrix with an one when element fullfills the condition < y
 coordinate_matrix_y_s(element_y_s) = 1;

 % Merge the y-coordinate_s and y-coordinate_g matrix
 coordinate_matrix_y = coordinate_matrix_y_g & coordinate_matrix_y_s;


 %% Calculating the Object which fullfills all the conditions x y 
 
 %get indice 
 index = find(coordinate_matrix_x & coordinate_matrix_y == 1);
 if length(index) > 1
  gmsgbox(['More than one point found in region around x = ' num2str(x,3) ', y = ' num2str(y,3) ],'Error','error');
 end
 
 if length(index) < 1
  gmsgbox(['No point found in region around x = ' num2str(x,3) ', y = ' num2str(y,3) ],'Error','error');
 end
  
end