function [prof_point]=get_profpoint(xc,yc,prof_lon,prof_lat)
%-----------------------------------------------------------
% funtion [prof_point]=get_profpoint(xc,yc,prof_lon,prof_lat)
% 
% input: xc (if llc, need to be in worldmap view, size [4*nx 4*nx]
%        yc (if llc, need to be in worldmap view, size [4*nx 4*nx]
%        [prof_lon,prof_lat]: pair of [lon,lat] of profiles
% output: nearest neighbor based on delaunay triangulation
% 
% 10.jan.2014: function to get prof_point given a grid xc,yc
% note: syntax to get worldmap view for llc from compact:
% xc=readbin('XC.data',[nx nx*13]);xc=patchface3D(nx,nx*13,1,xc,2);
%-----------------------------------------------------------

%get delaunay trigulation, save to file if needed
oldTri=1;newTri=0;
clear aa;aa=which('DelaunayTri');if(length(aa)>0);oldTri=0;newTri=1;end;
%if(length(dirGrid)>0);
%  fDelaunay=[dirGrid 'TRI_Stri' num2str(newTri) '.mat'];%load(fDelaunay);%Stri,TRI,kvalid
%else;
%  fDelaunay='not_saving_fDelaunay';fprintf('%s\n',fDelaunay);
%end;
%if(exist(fDelaunay)==0);
  if(oldTri);
    kvalid=1:prod(size(xc));
    TRI=delaunay(xc(kvalid),yc(kvalid));nxy=length(kvalid);
    Stri=sparse(TRI(:,[1 1 2 2 3 3]),TRI(:,[2 3 1 3 1 2]),1,nxy,nxy);
  elseif(newTri);
    kvalid=find(~isnan(xc));
    TRI=DelaunayTri(xc(kvalid),yc(kvalid));Stri=[];
  end;
%  if(length(dirGrid)>0);
%    save(fDelaunay,'TRI','Stri','kvalid');
%    fprintf('saving to file %s\n',fDelaunay);
%  end;
%else;
%  load(fDelaunay);
%end;

%get prof_point:
%make sure prof_[lon,lat] are in vector forms
sz=size(prof_lon);if(sz(2)>sz(1));prof_lon=prof_lon';end;
sz=size(prof_lat);if(sz(2)>sz(1));prof_lat=prof_lat';end;
if(oldTri);
  prof_point=dsearch(xc(kvalid),yc(kvalid),TRI,prof_lon,prof_lat,Stri);
elseif(newTri);
  prof_point=TRI.nearestNeighbor(prof_lon,prof_lat);
  prof_point=kvalid(prof_point);
end;

return
