function [ towerObs, towerNames ] = getTowerObservations( obsDir, period, info )
%getTowerObservations Get subset of tower obseration over the specified range of 8-day periods
%
%   [towerObs, towerNames] = getTowerObservations( obsDir, period )
%        This function returns a cell array of observations, with each cell
%        containing a vector, of size N, of CO2 concentrations. Tower's
%        names are also returned. 
%  INPUT :
%     obsDir : the directory where the observations are kept
%     period : a 2-element vector that contains the starting and ending
%              periods (each period is 8-days) of the duration of interest.
%              Only observations over this duration are returned.
%     info   : a struct with miscellaneous info like # of cells in US, # of
%             observations per period etc.
%
%  OUTPUT :
%      towerObs : a cell array, nTowers long. Each contains a vector
%                 containing [CO2] measurements at 3-hr intervals. Each
%                 8-day period as 8 x 8 = 64 observations per tower.
%      towerNames : Cell-array, with the names of the towers. Both the
%                 output cell arrays are nTowers long. 
%
%  Jaideep Ray
%  3/8/2012
% -------------------------------------------------------------------------

%% These are the tower names that we know of; also year of observation.
% Note that if synthetic obs, then the transport matrices are 2008 and
% Gurney fluxes used to generate synthetic observations are 2002.

towerNames = {'AMT', 'ARM', 'BAO', 'BRW', 'CDL', 'CEN', 'CHI', 'CVA', ...
              'EGB', 'ETL', 'FIR', 'FRD', 'GAL', 'HDP', 'HFM', 'KEW', ...
              'LEF', 'LJA', 'LLB', 'MAP', 'MEA', 'MET', 'MOM', 'NGB', ...
              'NWR', 'OZA', 'ROL', 'SBL', 'SCT', 'SNP', 'SPL', 'WBI', ...
              'WGC', 'WKT', 'YAH'};
year = 2008 ;

%% Loop over all towers and read in their observations

% 3-hrly indices corresponding to the period. Each day has 8 observations
% and each period has 8 days
startIndex = info.nObsPerPeriod * (period(1)-1) + 1;
endIndex   = info.nObsPerPeriod * period(2) ;

fprintf('Loading Observations for ');
for i = 1 : length(towerNames)
    fprintf('%s  ', towerNames{i});
    tmp = [obsDir towerNames{i} '_' num2str(year) '_coarse_obs.dat'];
    obs = load(tmp) ;
    towerObs{i} = obs(startIndex : endIndex) ;
end
fprintf('\n') ;

end

