% program for processing Lake Temperatures
% The original data should be contained in a variable 'data' that is 
% typically generated automatically when a .csv file is imported with 
% Matlab Import Wizard
% The 'data' varible has 4 columns, where the first 3 are year, month, and
% day, and the 4th column is temperature (may contain NaNs).


jd0=juliandate(2002,1,1);                      % <------ set starting date      
jd=juliandate(data(:,1),data(:,2),data(:,3)) - jd0+1;
temp=data(:,4);   %temp = original temperature data
%%

% eliminate all points below a certain value
d=max(size(temp));
for i=1:d
    if temp(i)<0                              % <-------- min T value
        temp(i)=NaN;
    end
end

%Remove NaNs
inarr(:,1)=jd;
inarr(:,2)=temp;
Torig=sergeiRemoveNaNs(inarr,2);     % Torig (jd,temp)= original data without NaNs
%%

% Merge multiple years into one
Tnew=Torig;
Tnew(:,1)=mod(Torig(:,1),365);
Tnew(:,2)=Torig(:,2);
Tnewsorted=sortrows(Tnew,1);
Torig=Tnewsorted;
Tbackup=Tnewsorted;

%%

%Set averaging window size (days)
ws=30;                                    % < ------- averaging window

%Iterate point eliminations 
for elim=1:9                              % <------ # of elimination rounds

% Compute running mean and variance: run(jd,rmean,npoints,rstd)
run = sergeiRunning (Torig,2,ws);
run=sortrows(run,1);%

%Remove outliers: Tfil(jd,Torig)
crit=mod(elim-1,2);                  % <-----  0-> 2*sigma test; 1 -> CHAUVENET'S CRITERION
Tfil = sergeiRemoveOutliers (Torig,run,2,ws,crit);

Torig=sortrows(Tfil,1);
end %end of elimination rounds
%%

% plot the results 
plot(Tbackup(:,1),Tbackup(:,2),'x')   %original series
axis([0 365 0 40]);
hold on

run=sortrows(run,1);

plot(Tfil(:,1),Tfil(:,2),'o','MarkerEdgeColor','g')         % filtered series

plot(run(:,1),run(:,2),'r','LineWidth',2)     % running mean
plot(run(:,1),run(:,2)+run(:,4),'r')   % mean plus-minus 1 sigma
plot(run(:,1),run(:,2)-run(:,4),'r')
plot(run(:,1),run(:,2)+2*run(:,5),'b')   % mean plus-minus 2 sigma-mu
plot(run(:,1),run(:,2)-2*run(:,5),'b')
hold off
