function [ nfield, lat2D, long2D ] = insertField( field, lats, longs, ttc, tbc )
%insertField Function to insert a field into a larger domain, discretized
%by a mesh of identical resolution
%   insertField( field, lats, lons, ttc, tbc ) takes a discretized field,
%   field, and inserts it into a larger domain discretized with a mesh of
%   identical resolution as field. It then returns the new field, in the
%   larger domain, along with matrices of the lat and long of the new
%   field.
%
% INPUTS :
%   field : a 2D array, with a cell-centered field
%   lats  : identical 2D array, with the lats of the cell-centers
%   longs : ditto, but longitudes
%   ttc   : [lat, long] to the top corner of the target domain
%   tbc   : [lat, long] of the bottom right corner of the target domain
%
% OUTPUTS:
%    nfield : 2D array, in the domain defined by (ttc, tbc). It'll have the
%    same spatial resolution as field, by the array will be square and
%    dyadic. The field is cell-centered
%
%    lat2D, long2D: identical 2D arrays, with the lat, long of the
%                   cell-centered mesh covering the larger domain
%
% Jaideep Ray, 5/21/2012
% -------------------------------------------------------------------------

% ---- make sure that the domain is bigger than the supplied field
assert( min(lats(:))  >= tbc(1) ) ;
assert( max(lats(:))  <= ttc(1) ) ;
assert( min(longs(:)) >= ttc(2) ) ;
assert( max(longs(:)) <= tbc(2) ) ;

% ---- Size and resolution of the incoming mesh
[nlat, nlong] = size(field) ;
dlat  = ( max(lats(:)) - min(lats(:)) + 1 ) / nlat ;
dlong = ( max(longs(:)) - min(longs(:)) + 1 ) / nlong ;

% ---- Is the incoming target domain discretizable using a dyadic grid at
% this resolution?
ntlongs = ( tbc(2) - ttc(2) + 1) / dlong ; dyadlongs = ceil(log(ntlongs)/log(2)) ;
ntlats  = ( ttc(1) - tbc(1) + 1) / dlat  ; dyadlats  = ceil(log(ntlats) /log(2)) ;
if ( (2^dyadlongs ~= ntlongs) || (2^dyadlats ~= ntlats) )
    error('insertField() : A dyadic grid cannot be put on the target domain') ;
end

% ---- OK, the incoming domain can have a dyadic grid. make the grid
tlats  = tbc(1) : dlat  : ttc(1) ;
tlongs = ttc(2) : dlong : tbc(2) ;
[long2D, lat2D] = meshgrid(tlongs, tlats) ;

% ---- now make the field
nfield = zeros(ntlongs, ntlats) ;
row = lats(:) - tbc(1) + 1 ;
col = longs(:) - ttc(2) + 1 ;
ind = sub2ind( [ntlongs, ntlats], row, col) ;
nfield(ind) = field(:) ;
end

