function [ nf, nflatlon, ind ] = getFieldInLower48( field, flatlon, tlatlon )
%getFieldInLower48 Return a vector or field on a set of points at a subset
%of points
%   getFieldInLower48( field, flatlon, tlatlon ) returns a set of values
%   field, described at a set of points, flatlon, at a subset of points
%   tlatlon.
%
% INPUTS :
%   field : a field described as a vector ncells long. if we have a series
%           of vectors, M, long, described on the same set of points then field is
%           a ncells X M matrix. This could be used to represent a bunch of fields
%           or a time-varying fields
%   flatlon : ncells x 2 matrix, giving the lat/long of the cells where
%             field is described
%   tlatlon : A vector q x 2 in size, containing the lat-lon of a subset of
%             points contained in flatlon. These points must be exact since
%             no interpolation will be done. If the supplied target
%             lat-long are not a subset, then the intersection of flatlon
%             and tlatlon will be supplied.
% OUTPUTS 
%   nf    : Q x M vector containing the value of field at the points that
%           lie at the intersection of flatlon and tlatlon
%   nflatlon : the Q x 2 vector containing the lat/lon of the cells whose
%              values we found
%   ind    : row numbers in field that give nf i.e nf = field(ind, :)
%
% Jaideep Ray, 5/23/2012
% =========================================================================

% ---- which rows in flatlon are in tlatlon
[~, ind, i2] = intersect(flatlon, tlatlon, 'rows') ;

% ---- Make sure that the rows have not been sorted i.e. flatlon(ind) ==
% tlatlon(i2). 
[~, i] = setdiff(flatlon(ind), tlatlon(i2)) ;
assert( numel(i) == 0 );

% ---- Also make sure we found all the points
if (numel(i2) ~= size(tlatlon, 1))
    warning('getFieldInLower48() : All points in target set were not found') ;
end

% ---- Extract the right set of rows from field and flatlon
nf = field(ind, :) ;
nflatlon = flatlon(ind, :) ;
end

