function [ Hmatrices, Hlatlon ] = getHMatrices( hDir, range, towerNames, ...
                                                info)
%getHMatrices Load H matrices for Lower48 for the 8-day periods specified
%by range.
%   Hmatrices = getHMatrices( Hdir, range, buffer, towernames ) returns H matrices for
%   35 towers. The H matrices are for the emissions in the Lower48 states.
%   The H matrices are for only the 8-day periods in 2008 specified by
%   range(1) : range(2), both inclusive.
%
%   The H matrices are for 2008. They are the sensitivity of tower
%   observations to emissions, averaged over 8-day periods, over a region
%   1deg X 1deg. There are 360/8 = 45 periods, plus a 8-day period in 2007,
%   for lag purposes. Thus the H matrix for each tower has 46 x 816 = 37536
%   cols, 816 being the # of 1deg x 1deg cells in Lower48. The H matrix has
%   8 readings per day, for 360 days in 2008.
%
% INPUTS :
%   hDir : which directory are the H matrices in?
%
%   range: Which range of 8-day periods are you interested in? Period # 1
%          is Jan 1 - Jan 8, 2008. Period # 0, also valid, is Dec 24 - Dec 31st,
%          2007. This determines which rows of the year-long H matrix is
%          returned to you. It also partially determines the columns.
%
%   towerNames : names of the towers you are interested in.
%
%   info : a struct with miscellaneous info like # of cells in US, # of
%   observations per period etc.
%
% OUTPUTS :
%    Hmatrices : a cell array, nTowers long, with each cell containing  a
%    matrix [8 * 8 * nperiods] X [ (nperiods + buffer)*816 ] H matrix for a
%    tower.
%
%    Hlatlon : a 816 x 2 matrix containing the latitude and longitude of
%    the locations where H matrices expect their emissions.
%   
%   Jaideep Ray
%
%  3/15/2012
% -------------------------------------------------------------------------

%% Note that if synthetic obs, then the transport matrices are 2008 and
% Gurney fluxes used to generate synthetic observations are 2002.
year = 2008 ;

% Loop over all towers and read in their H matrices
fprintf('Loading H matrices for ');
for i = 1 : length(towerNames)
    fprintf('%s  ', towerNames{i});
    tmp = [hDir towerNames{i} '_' num2str(year) '_coarse_Lower48.mat'];
    load(tmp) ; % Contains a matrix called Hlower48 ;
    Hmatrices{i} = extractH(Hlower48, range, info) ;
end
fprintf('\n') 

load([hDir '/USlatlon_Lower48.mat']) ; % Load a var called USlatlon
Hlatlon = USlatlon ;
end

%% Private function to extract the right part of the H matrix
%
function H = extractH(Hfull, range, info) 
% extractH() Extracts the right range of 8-day periods from given year-long
% H matrix
%
%   info : a struct with miscellaneous info like # of cells in US, # of
%   observations per period etc.
%
nper = size(Hfull, 1) / info.nObsPerPeriod ;
ncells = size(Hfull, 2) / (nper+1) ;

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

startIndexCol = ncells * range(1)  + 1 ;
endIndexCol   = ncells * (range(2) + 1) ;

H = Hfull(startIndexRow:endIndexRow, startIndexCol:endIndexCol) ;
end
