accumSeasonTS(datain,fn)

%% accumSeasonTS- aggregate timeseries to solstices and equinoxes
%
% This function will output cumulative seasonal totals of data. It's
% applciation is limited default parameterization of ACCUMARRAY matlab
% function (i.e. sum totals).  It's certainly possible to modify this
% function to allow the user to apply different statistics in the
% ACCUMARRAY function, but not done here.
%
% Additionally, this function uses another function (EQNSOL) to calculate the
% summer/winter solstices and spring/fall equinoxes per year, then
% sum up the data within each season.
%
% Function does not descriminate on partial seasons.  Partial seasons that
% are either because of beginning or ending dates, or because of missing
% data are not censored.  Partial seasons will be output along with all
% complete seasons. It's up to the user to either pre-process the data or
% post process the results.
%
% Expected inputs is a timeseries of data in two columns: a decimal date,
% and the values to be accumulated.
%
% datain = n x 2
%   datain(:,1) = decimal date/time
%   datain(:,2) = values to sum by season
% fn = supplied string specifying a valid summary statistic function.
% Example: @sum, @max, @min, etc.
% 
% dataout = n x 2
%   dataout(:,1) = decimal date/time of the start of the season
%   dataout(:,2) = values summed for the season
%
% Acknoledgements:
%
% Cardillo G. (2007) Equinoxes and Solstices: compute the date and time of
% equinoxes and solstices. 
% http://www.mathworks.com/matlabcentral/fileexchange/17977
%   
% 
% Written by:
%   Jeff Burkey
%   King County- DNRP
%   email: [e-mail]
%   May 30, 2009
% 
% Syntax:
%   [dataout] = accumSeasonTS(datain)
function [dataout] = accumSeasonTS(datain,fn)
    % function will sum up values to daily time steps
    [y m] = datevec(datain(:,1));
    yrs = unique(y);
    nyrs = size(yrs,1);
    yseas = zeros(nyrs,5);
    datain = [datain zeros(size(datain,1),2)];
    
    yseas(1,:) = [yrs(1) eqnsol(yrs(1))'];
    
    datain(datain(:,1) < yseas(1,2),3)=1;
    datain(yseas(1,2)<= datain(:,1) & datain(:,1) < yseas(1,3),3)=2;
    datain(yseas(1,3)<= datain(:,1) & datain(:,1) < yseas(1,4),3)=3;
    datain(yseas(1,4)<= datain(:,1) & datain(:,1) < yseas(1,5),3)=4;
    % depending on when the data started the starting season count
    % could start anywhere between 1 and 4.
    
    datain(datain(:,1) < yseas(1,2),4)=4;
    datain(yseas(1,2)<= datain(:,1) & datain(:,1) < yseas(1,3),4)=1;
    datain(yseas(1,3)<= datain(:,1) & datain(:,1) < yseas(1,4),4)=2;
    datain(yseas(1,4)<= datain(:,1) & datain(:,1) < yseas(1,5),4)=3;
    
    for n=2:nyrs
        yseas(n,:) = [yrs(n) eqnsol(yrs(n))'];
        
        % index the season and assign to a running count. This is
        % necessary because the winter season spans calendar years, and
        % the fall season spans water years.
        datain(yseas(n-1,5)<= datain(:,1) & datain(:,1) < yseas(n,2),3)=((n-1)*4)+1;
        datain(yseas(n,2)<= datain(:,1) & datain(:,1) < yseas(n,3),3)=((n-1)*4)+2;
        datain(yseas(n,3)<= datain(:,1) & datain(:,1) < yseas(n,4),3)=((n-1)*4)+3;
        datain(yseas(n,4)<= datain(:,1) & datain(:,1) < yseas(n,5),3)=((n-1)*4)+4;
        
        % set value to season 1=spring, 2=summer, 3=fall, 4=winter
        datain(yseas(n-1,5)<= datain(:,1) & datain(:,1) < yseas(n,2),4)=4;
        datain(yseas(n,2)<= datain(:,1) & datain(:,1) < yseas(n,3),4)=1;
        datain(yseas(n,3)<= datain(:,1) & datain(:,1) < yseas(n,4),4)=2;
        datain(yseas(n,4)<= datain(:,1) & datain(:,1) < yseas(n,5),4)=3;
        
    end
    % finish up assigning seasons for last winter season.
    datain(datain(:,1) >= yseas(nyrs,5),3)=nyrs*4+1;
    datain(datain(:,1) >= yseas(nyrs,5),4)=4;
    
    [b1, m1, n1] = unique(datain(:,3),'rows');
    dataout = [datain(m1,1) accumarray(n1,datain(:,2),[],fn)];
end