function [ towerPreds, emissionsPreds, emissionsTrue2D, emissionsRecons2D, ...
           figNo ] = estimateEmissionReconsError( prob, sol, emissionsLower48, ...
                               latlonLower48, topcorner, botcorner, periods, ...
                               figNo, visualize, saveStuff, oDir, fs)
%estimateEmissionReconsError Estimate emission reconstruction error and
%plot and save it
%
%  [ towerPreds, emissionsPreds, emissionsTrue2D, emissionsRecons2D, ...
%    figNo ] = estimateEmissionReconsError( prob, sol, emissionsLower48, ...
%                               latlonLower48, emissionsPreds, topcorner, ...
%                               botcorner, periods, figNo, visualize, saveStuff)
%
%  This function estimates the error in the reconstruced emissions over
%  Lower48. It can plot the error in emissions aggregated over US and save
%  the data to disk if needed.
%
% INPUTS :
%   prob : struct, with information in it
%   sol  : long vector with all wavelet weights over all periods, inclusing
%          buffer periods
%   emissionsLower48     : true emissions in lower 48
%   latlongLower48       : grid-box centers in the lower 48 [lat, long]
%   topcorner, botcorner : [lat, lon] of topleft and bottom right corners
%   periods              : 2 element vector with the starting and ending 
%                          periods of interest
%   figNo                : fig number to start with
%   visualize            : boolean, plot stuff?
%   saveStuff            : boolean, save files (data and figures)
%   oDir                 : output directory
%   fs                   : font size
%
%  OUTPUTS :
%    towerPres         : predicted [CO2] at the towers using reconstructed
%                        emissions. Cell array, nTowers long
%    emissionPreds     : Reconstructed emissions, in one long array
%    emissionsTrue2D   : 3D matrices, emissions in a 2D field with 3rd dimension 
%                        the time axis. True emissions
%    emissionsRecons2D : Ditto, but reconstructed emissions
%    figNo             : Figure number for the downstream plots to use
%
%  Jaideep Ray, 6/27/2012
% -------------------------------------------------------------------------

% ---- Reconstruct emissions and predict tower concentrations. posterior
% predictive check
[towerPreds, emissionsPreds] = getReconstruction(sol, prob) ;

% ---- Convert the emissions into 2D fields, over time i.e. 3D arrays
[emissionsTrue2D, emissionsRecons2D] = convertTo2DArray(emissionsLower48, ...
                               latlonLower48, emissionsPreds, topcorner, ...
                               botcorner, prob) ;

% ---- Aggregate emissions over Lower48 and check the relative error, true
% versus reconstructed. Print them to screen, for the periods of interest

periodList = periods(1) : periods(2) ;
frac = zeros(length(periodList), 1) ;
info = prob.info ;
for i = info.buffer + 1 : (info.nPeriods - info.buffer)
    % What fraction of emission is reconstructed?
    f = (emissionsRecons2D(:, :, i) - emissionsTrue2D(:, :, i)) ; % flux difference
    
    % Flux difference over US / total emissions
    frac(i -info.buffer) = sum(f(:)) / sum( sum(emissionsTrue2D(:, :, i)) ) ;
end

% --- Print out total emissions
fprintf('  Period     US total emissions est. err [per cent]\n');
for i = 1 : length(periodList)
    fprintf('   %02d      %7.2e\n', periodList(i), 100*frac(i)) ;
end

% ---- Plot the reconstruction errors as a function of time. Always do
% this, regardless of whether visualize is true or false
figNo = figNo + 1; figure(figNo);
plot(periodList, 100*frac, '-k'); 
title('US total emissions estimation error (%)', 'FontSize', fs) ; 
xlabel('Period #', 'FontSize', fs) ; 
ylabel(' 100*(Reconstruction - True)/True ', 'FontSize', fs) ;
set(gca, 'FontSize', fs);   % Increase fonts

%  ---- Save reconstruction
if (saveStuff == true)
    str = sprintf('ReconsError_P%02d-P%02d.jpg', periods(1), periods(2));
    print('-djpeg', [oDir '/' str]);
    
    tmp = [periodList', frac];
    str = sprintf('ReconsError_P%02d-P%02d_data.dat', periods(1), periods(2));
    save([oDir '/' str], 'tmp', '-ascii');
end        
end

