function [ emissions, lat2D, long2D ] = getEmissionsOnBaseGrid( edir, ...
                                                   periods, latlon)
%getEmissionsOnBaseGrid() Get emissions on the grid where the bases are
%described.
%   getEmissionsOnBaseGrid( edir, periods, bases, latlong) reads in Gurney
%   fluxes from file, and puts them on the dyadic grid, in the large domain
%   where the bases are defined.
%
% INPUTS :
%   edir : directory where the emissions are stored
%   periods : a vector, 2 long, specifying the range of periods of interest
%   latlon : a matrix, ncells X 2 in size, containing the cell-center
%             lat/long of the gridcells on the dyadic grid where the bases
%             are described
% 
% OUTPUTS:
%   emissions : ncells x nperiods matrix of emissions in each period
%   lat2D, long2D : 2D matrices containing the lat/long of the cell-centers
%                   where emissions are described. dyadic matrices.
%
% Jaideep Ray, 5/21/2012
% -------------------------------------------------------------------------

% ---- Get emission field, for each of the period; estimate total emission
[emissions_tmp, lats2D_tmp, longs2D_tmp] = getEmissionsFromFile(edir, periods) ;

% ---- Put the emission field on the same dyadic mesh as the nightlights
ttc = [max(latlon(:, 1)), min(latlon(:, 2))] ;
tbc = [min(latlon(:, 1)), max(latlon(:, 2))] ;
gridsize = sqrt( size(latlon, 1) ) ;
emissions = zeros(gridsize, gridsize, size(emissions_tmp, 2)) ;

% ---- Loop over all time-slices (periods and fill the 3D emissions up
for iperiod = 1 : size(emissions_tmp, 2)
    tmp = emissions_tmp(:, iperiod) ;
    tmp = reshape(tmp, size(lats2D_tmp, 1), size(lats2D_tmp, 2)) ;
    [emissions(:, :, iperiod), lat2D, long2D] = insertField(tmp, ...
                                       lats2D_tmp, longs2D_tmp, ttc, tbc) ;
end
emissions = reshape(emissions, gridsize*gridsize, size(emissions_tmp, 2));
end

