hdrload_mean(files,varargin);

function [H,x1,Ybar] = hdrload_mean(files,varargin);
%hdrload_mean
%
%[H,x,Ybar] = hdrload_mean(files) loads and averages the data matrices in
%cellstring array = {file1; file2; ...}. Files contain tab deliminited
%data and header readable by hdrload. Ybar is the averaged data
%corresponding to x1. x1 is the first column of file 1. Data from other
%files is interpolated on x1. H is a cellstring array that contains the 
%headers of the individual files.
%
%[H,x,Ybar] = hdrload_mean(files,col,range) only opperates on the data whose
%value in column col falls within the values in vector range. Example:
%[H,dbar] = hdrload_mean(files,1,[0,3.4]) loads the data in files and
%averages data that has column 1 value between 0 and 3.4. If the files do
%not have identical data points in column col, it picks file(1) to
%determine the ordinate and interpolates the other data sets to pick out
%the values to averaged with the file(1) data set.
%
%Files are assumed to have same number of data columns, but not
%necessarrily the same number of rows.
%
% requires hdrload.m by Jeff Daniels
%
% Tom Allison 7/22/2008


numfiles = length(files);

%load first file
[h1,d1] = hdrload(files{1});

[r1,c1] = size(d1);

if isempty(varargin)
    col = 1;
    range = [d1(1,1),d1(r1,1)];
elseif length(varargin) == 2
    col = varargin{1};
    range = varargin{2}
else
    disp('Incorrect number of input arguments');
end

range = sort(range);
%construct relevant d matrix

r = find(d1(:,col) >= range(1) & d1(:,col) <= range(2)); %find rows in range

if (col == 1) %pick out starting Y matrix
   Y = d1(r,2:c1);
elseif (col == c1)
   Y = d1(r,1:(c1-1));
else    
   Y = d1(r,[1:(col-1),(col+1):c1]);
end

x1 = d1(r,col); % data points to interpolate to
H = {h1};

if numfiles == 1
    Ybar = Y;
    return
end

for j = 2:numfiles
    [hnow,dnow] = hdrload(files{j});
    xnow = dnow(:,col); % pick out ordinate to interpolate on
    if (col == 1) %pick out Y matrix to be interpolated
        Ynow = dnow(:,2:c1);
    elseif (col == c1)
        Ynow = dnow(:,1:(c1-1));
    else    
        Ynow = dnow(:,[1:(col-1),(col+1):c1]);
    end
    Y = Y + interp1(xnow,Ynow,x1);
    H{j} = hnow;
end

Ybar = Y/numfiles;