function OUT = avgout(IN,flag,TimeLim)
% AVGOUT    - Process Eqpac data from READOUT
%
% AVGOUT processes data from structures created by READOUT. The data in most 
% fields is returned as an averaged of the 15 seconds of data taken every 15 minutes.
% Alternatively it can be returned as an single daily average of data taken
% between 10am and 2pm (local).
%
% The filename must contain the longitude in the 11-13 positions. 
%
% Use as: out = avgout(in)
%         out = avgout(in, Lon)
%         out = avgout(in, Lon, TimeLims)
% Inputs: in = Data structure with at least the following fields:
%                 Filename (must be the first field)
%                 TIME_TAG
%         Lon = Optional argument. Degrees longitude of sampiling location.
%               If specified, a daily average between 
%               10 and 2pm (Local to the specified Longitude) is returned.
%               If not specified, is empty, or is NaN then an average of data occuring
%               during each minute is returned.
%         TimeLims = Optional argument. Two element vector specifying the min
%               and max time to avg between (default = [10 14] , ie. 10am and 2pm)
%
% NOTE: The time field in the input structure X must be named TIME_TAG and the
% structure must also have the first field named Filename.

% Brian Schlining
% 13 Jan 1998

if nargin < 2 | isempty(flag) | isnan(flag)
   flag = 0;
else
   Lon  = flag;
   flag = 1;
end

%=========================================
% Error Checking on the TimeLims argument
%=========================================
if nargin < 3                             % Default
   StartTime = 10;
   EndTime   = 14;
else
   if TimeLim(1) > TimeLim(2)
      warning('  AVGOUT: TimeLims must be [smaller larger]')
      warning('  AVGOUT: TimeLims have been sorted')
      TimeLim = sort(TimeLim);
   end
   if TimeLim(2) > 24 | TimeLim(1) < 0
      error('  AVGOUT: Invalid Time limit')
   end
   StartTime = TimeLim(1);
   EndTime   = TimeLim(2);
end


%=========================================================================
% Need the date, down to the Minute, to find the ~15 1 sec readings taken
% every 15 minutes
%=========================================================================
OUT.Filename = IN.Filename;            % hang onto the filename
switch flag
case 0
   Offset = 0;
otherwise
   flag = 1;
   Offset = gmtoffset(Lon);            % Calculate the offset in days from GMT
end

Fields       = fieldnames(IN);         % Get all the field names
NumOfFields  = length(Fields);         % Count the number of fields
[TIME_TAG, DateOrder] = unique(IN.TIME_TAG); % Make sure the data is sorted
NumOfPoints  = length(TIME_TAG);    % Count the number of readings
[Year, Month, Day, Hour, Minute, Second] = datevec(TIME_TAG + Offset);


n = 1;                                 % Counter for the start of each 15 minute data set
m = 1;                                 % The position we're writing the data to

fprintf(1,'(AVGOUT)')
OK = 1;                                
while OK
   
   % Preallocate LOTS of memory on the first run through the infinite while loop. If
   % this isn't done the program is REALLY SLOW. Here we're allocating a structure 
   % the same size as IN and filling it with NaN's
   if m == 1
      
      
      for i = 1:NumOfFields
         if ~strcmp(Fields(i),'Filename')
            % Preallocate memory
            OUT = setfield(OUT,{1},Fields{i},{1:NumOfPoints,1},NaN);         
            Xin(:,i) = getfield(IN,Fields{i},{DateOrder});
         end
      end
      % Sort the rows based on date
      %Xin = Xin(DateOrder,:);
   end
   
   switch flag
   case 0   % average the data from every 15 minutes
      
      % Find all data points that occur during any given minute
      j = find(Year == Year(n) & Month == Month(n) & Day == Day(n) & ...
         Hour == Hour(n) & Minute == Minute(n));
      
      % Get the data that occurs during each minute, take its average and then 
      % write it out to the new structure, OUT.
      if length(j) > 1
         Xout(m,:)  = mean(Xin(j,:));
      else
         Xout(m,:) = Xin(j,:);
      end
      
      n = n + length(j);   % position for the next set of replicates
      m = m + 1;           % position index for OUT
      
   case 1   % Create a daily average for data between 10am and 2pm
      
      % Find all data points between 10:00 and 14:00 for each day
      % At some point we'll need to flag this if not enough data
      % points occur during this period
      j = find(Year == Year(n) & Month == Month(n) & Day == Day(n) & ...
         Hour >= StartTime & Hour <= EndTime);
      
      % Move pointer to the next day
      k = find(Year == Year(n) & Month == Month(n) & Day == Day(n));
      
      if ~isempty(j)
         % Get the data that occurs during each minute, take its average and then 
         % write it out to the new structure, OUT.
         if length(j) == 1
            Xout(m,:) = Xin(j,:);
         else
            Xout(m,:)  = nanmean(Xin(j,:));
         end
         n = n + length(k);
         m = m + 1;
      elseif isempty(j) & ~isempty(k)
         n = n + length(k)';
      else
         n = n + 1;
      end
      
      
   end
   
   
   if ~rem(m,200)
      fprintf(1,'.')
   end
   
   
   % Exit the loop when the end of the data is reached. 
   if n > NumOfPoints
      fprintf(1,'\n')
      % Remove all the trailing NaN's before returning OUT
      timeCol = strmatch('TIME_TAG',Fields);
      notNaN = find(~isnan(Xout(:,timeCol)));
      
      for i = 1:NumOfFields  
         if ~strcmp(Fields(i),'Filename')
            OUT = setfield(OUT,Fields{i},Xout(notNaN, i));
         end
         
      end
      if flag == 1
         OUT.TIME_TAG = OUT.TIME_TAG;
      end
      
      
      return      % Exit 
      
   end
   
   
end