function Cp = cpmatrix( hdf_filename, tap_ind )

% This function was adapted from 'load_hdf_file.m' by Dave Morrish of
% University of Western Ontario, which was distributed as part of the 
% NIST Database Software. 

sd_id=hdfsd('start',hdf_filename,'read'); % obtain a file ID for specified file
if sd_id <=0
    error(['Could not open >' hdf_filename '< for reading: check paths and permisions.']);
end
[ndatasets,natt,sata]=hdfsd('fileinfo',sd_id); % obtain file information

% List of datasets needed for assembling Cp matrix:
datasets_list = {'Ts_Multiplier','Roof_Height_Velocity_Ratio','Time_Series'};
%datasets_list = {'Ts_Multiplier','Roof_Height_Velocity_Ratio','Time_Series','Tap_Position_List'};

% Loop through datasets and create a field in the structure 'hdf_struct'
% for each dataset in the list above; ignore other datasets as well as any
% attributes provided for the datasets:

hdf_struct.Time_Series = []; % initialize data structure
for sds_idx=0:ndatasets-1
    sds_id=hdfsd('select',sd_id,sds_idx);
    [dsname,dsndim,dsdims,dstype,dsatts,stat]=hdfsd('getinfo',sds_id);
    % check if current dataset name is in 'datasets_list':
    if ismember(1,strcmp(dsname,datasets_list))
        % extract read dataset and store in 'hdf_struct'
        ds_start=zeros(1,dsndim); ds_stride=[]; ds_edges=dsdims;
        [ds_data,status]=hdfsd('readdata',sds_id,ds_start,ds_stride,ds_edges);
        hdf_struct=setfield(hdf_struct,dsname,ds_data);  
    end
end
hdfml('closeall');

int_scl = hdf_struct.Ts_Multiplier; % factor for converting time series from integer format
ht_scl = (hdf_struct.Roof_Height_Velocity_Ratio)^-2; % factor for scaling pressure coefficients to eave height velocities
ts_scl = int_scl*ht_scl; % total scale factor for time series

Cp = ts_scl*double(hdf_struct.Time_Series(:,tap_ind)'); % Pressure coefficient time series
% tap_num = double(hdf_struct.Tap_Position_List(tap_ind)');