function [ figNo ] = plotTowerCO2Conc(info, towerNames, towerObs, ...
                                      towerPreds, periods, figNo, ...
                                      visualize, saveStuff, oDir, ...
                                      fs, lfs, sc, lc)
%plotTowerCO2Conc Posterior predictive check at a few towers
%   figNo = plotTowerCO2Conc(info, towerObs, towerPreds, periods, ...
%                             visualize, saveStuff) plots the CO2
% concentrations, as observed and as recosntructed at a few towers. 
%
% INPUTS :
%   info     : struct with info
%   towerNames : tower's names
%   towerObs : cell array, with period-averaged [CO2] at 35 towers
%   towerPreds : predictions at towers, to be compared with towerObs
%   periods  : [start, end] of the periods of interest
%   figNo    : figNo to use
%   visualize : Should we plot?
%   saveStuff : save to file?
%   oDir      : output directory
%   fs, lfs   : font size, legend font size
%   sc, lc   : symbol and line color (cell array)
%   
%
% OUTPUTS :
%   figNo : Start downstream with this figNo.
%
%  Jaideep Ray, 6/27/2012
% -------------------------------------------------------------------------

% ---- How many towers to plot on one plot
towerIndices = floor(1 : info.nTowers/3 : info.nTowers) ;
figNo = figNo + 1 ; 
buf = info.buffer * info.nObsPerPeriod ;

% ---- Calculate the discrepancy between obs and posterior predictive. 
% First, an array to store it.
disc = zeros( (periods(2)-periods(1)+1)*info.nObsPerPeriod, ... 
              length(towerIndices) );
          
% ---- Loop over 3 towers. Plot obsered [CO2]; also, store the observation 
% in discrepancy array          
figure(figNo) ;
for i = 1 : length(towerIndices)
    allObs = towerObs{towerIndices(i)} ; 
    obsInPeriodOfInterest = allObs(buf+1 : end-buf) ;
    disc(:, i) = obsInPeriodOfInterest ;  % Save the observations
    plot(obsInPeriodOfInterest, sc{i}) ; hold on ;
    legendString{i} = [towerNames{towerIndices(i)} '; obs'] ;
end
    
% ---- Now, plot the prediction and calculate discrepancy, per tower, per
% period
for i = 1 : length(towerIndices)
    allPred = towerPreds{towerIndices(i)} ;
    predInPeriodOfInterest = allPred(buf+1 : end-buf) ;
    disc(:, i) = (disc(:, i) - predInPeriodOfInterest) ; % Calc discrepancy
    plot(predInPeriodOfInterest, lc{i}) ; hold on ;
    legendString{length(towerIndices)+i} = [towerNames{towerIndices(i)} '; pred'] ;
end

% ---- Increase fonts in graph and put in legend
legend(legendString, 'Location', 'NorthWest') ; legend boxoff ;
str = sprintf('Anthropogenic CO_2 concentrations at 3 towers (ppm) Periods %d - %d\n', ...
               periods(1), periods(2));
title(str, 'FontSize', lfs);
set(gca, 'FontSize', fs);

% ---- Should I save the plot to disc
if (saveStuff == true)
    str = sprintf('TowerPredsP%02d-P%02d.jpg', periods(1), periods(2));
    print('-djpeg', [oDir '/' str]);
end

% ---- Calculate the rms error for each of the towers
fprintf('\n\nRoot mean square posterior predictive performance at towers\n');
for i = 1 : length(towerIndices)
    allObs = towerObs{towerIndices(i)} ; 
    obsInPeriodOfInterest = allObs(buf+1 : end-buf) ;
    rms_obs = sqrt(mean(obsInPeriodOfInterest.^2)) ;
    rms_error = sqrt(mean(disc(:, i).^2) );
    fprintf('Tower %d : Error = %d,  Obs = %d, Rel. err = %f\n', ...
        i, rms_error, rms_obs, rms_error / rms_obs);
end

end

