function  [out] = sergeiRunning (inarr,dim,ws)
% This function scans a two-column array (julian day, temperature) along the dimension dim (1 or 2) and
% computes the running mean and running variance. The averaging window ws should be specified 
% in terms of time (e.g., number of days).   
% For each averaging interval, the function records the number of data 
% points that fall within that interval.
% The "julian day" column in the input array should contain values between
% 0 and 365.
% The output is an array with 365 rows and columns that correspond to:
% 1) julian day 2) running mean 3) number of points found within the
% averaging window 4) running variance 5) "error on the mean"  

d=max(size(inarr));
% make the data cyclic by extending the sereies for another year  
extended=cat(1,inarr,inarr);
jd=extended(:,1);
jd(d+1:2*d,1)=jd(1:d)+365;
T=extended(:,2);

%cycle over the days of the year
i=1;
for m=1:365
    %find the closest matching days of the year in the data array
    while (jd(i)>m && i>1)
        i=i-1;
    end
    while jd(i)<m
        i=i+1;
    end
           
    %calculate running mean by summing up data within the averaging window
    k=0;sum=0;count=0;
    while jd(i+k)<=jd(i)+ws
        sum=sum+T(i+k);
        count=count+1;
        k=k+1;
    end
    time(m)=jd(i);
    rmean(m)=sum/count;
    npoints(m)=count;
end 

% now use the running mean to calculate the running variance
i=1;
for m=1:365
    while (jd(i)>m && i>1)
        i=i-1;
    end
    while jd(i)<m
        i=i+1;
    end   
      
    k=0;sum2=0;count=0;
    while jd(i+k)<=jd(i)+ws
        sum2=sum2+(T(i+k)-rmean(m))^2;
        k=k+1;
    end
    
    if npoints(m)<=1
        rstd(m)=NaN;
    else
        rstd(m)=sqrt(sum2/(npoints(m)-1));
    end

end   

% center the running values to the middle of the averaging window
time(1:365)=mod(time(1:365)+ws/2,365);

% form the output array
out(:,1)=time(1:365);
out (:,2) = rmean(1:365);
out(:,3)=npoints(1:365);
out(:,4)=rstd(1:365);
out(:,5)=rstd(1:365)/sqrt(npoints(1:365));