function Xout = avgoasis2(Xin, Lon, TimeLim)
% AVGOASIS2 - Create a daily average of data in oasis files, high precision version
%
% Use as: avgoasis2(X, Lon, TimeLim)
%         avgoasis2(X, Lon);
%
% Inputs: X   = matrix of data (from READOASIS), 1st column is julian day,
%             last column is year.
%         Lon = Longitude of readings (scalar) used to calculate the offset
%             from GMT. 
%         TimeLims = a scalar (fractional day). The Local noon is calcualted by
%             using NOON_. TimeLim specifies the time on either side of noon to
%             average between. Default = 2/24 (2 our offset)
%
% Outputs: Xo = daily average of column wise data.
%
% Requires: GMTOFFSET, JUL2DATE, Oceans toolbox

% Brian Schlining
% 17 May 1999


if nargin < 2 | isempty(Lon) | isnan(Lon)
   error('(AVGOASIS) usage: avgoasis2(X,Lon) OR avgoasis2(X,Lon,TimeLim)');
end

%=========================================
% Error Checking on the TimeLims argument
%=========================================
if nargin < 3                             % Default
   TimeLim = 2/24;   
end

Offset = gmtoffset(Lon);                    % Calculate the offset in days from GMT

%===================================
% Get the date number of each point
%===================================
d        = jul2date(Xin(:,1),Xin(:,end));   % Convert Julian day to matlab's format

GoodDate = find(~isnan(d));                 % Use only finite numbers
d        = d(GoodDate);       
Xin      = Xin(GoodDate,:);

[d, DateOrder] = sortrows(d);               % Make sure the data is sorted by time
Xin      = Xin(DateOrder,:);

[Year, Month, Day, Hour, Minute, Second] = datevec(d); % GMT

LAN       = noon_(Year,Month,Day,Lon*-1,0)/24;   % Get local noon in days GMT
StartTime = (LAN - abs(TimeLim))* 24;       % in hours
EndTime   = (LAN + abs(TimeLim))* 24;       % in hours

%================
% Setup for loop
%================
n = 1;                                     % Counter for the start of each 15 minute data set
m = 1;                                     % The position we're writing the data to
[NumOfPoints NumOfFields] = size(Xin);
Xout = ones(size(Xin))*NaN;                % Preallocate memory


fprintf(1,'(AVGOASIS2)')
OK = 1;      
while OK
   %======================================================
   % Create a daily average for data between 10am and 2pm
   %======================================================
   
   % Find all data points between 10:00 and 14:00 for each day
   j = find(Year == Year(n) & Month == Month(n) & Day == Day(n) & ...
      Hour >= StartTime(n) & Hour <= EndTime(n));
   
   % 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
   
   if ~rem(m,200)
      fprintf(1,'.')
   end
   
   if n > NumOfPoints
      fprintf(1,'\n')
      % Remove all the trailing NaN's before returning OUT
      NotNaN = find(~isnan(Xout(:,1)));
      Xout = Xout(NotNaN,:);
      return
   end
end