function [ranom, fep, sts]=f_Tanom(Ty, Tr, vargin)

%* function f_Tanom(fn, Ty, Tr, vargin)
%
% compute anomaly and best fit linear trend
% assumes 999.9 is missing data
%
% Ty        years for temperature record
% Tr        temperature record to plot
% vargin    start year for mean, 1951 default
% vargin    end year for mean, 1980 default
%
% return values
% anom      temperature anomaly to plot
% fep       best fit endpoints for graphing
% p(1)      slope of trend line 
% rsq       coefficient of determination ("r-squared")


goodies=find(Tr~=999.9);

%if length(vargin)==2
    mstart=find(Ty==vargin(1));
    mend=find(Ty==vargin(2));
%else
%    mstart=find(Ty==1951);
%    mend=find(Ty==1980);
%end

Tm=mean(Tr(mstart:mend));
anom=Tr(goodies)-Tm;

% linear trend
[p,S] = polyfit(Ty(goodies),anom,1);
f=p(1)*Ty(goodies)+p(2);

fep=[Ty(min(goodies)) f(1); Ty(max(goodies)) f(length(f))];

%* metrics
anomb=mean(anom);
sst=sum((anom-anomb).^2);   % total sum of squares
ssg=sum((f-anomb).^2);      % regression sum of squares
ssr=sum((f-anom).^2);    % residual sum of squares

% coefficient of determination (goodness of fit of linear model)
rsq=1-ssg/sst;


%* sgnificance of trend (for t-test)
%  f    is the model 
%  anom is the sample
sse=sqrt(ssr)/(length(anom)-2)/ sqrt(sum((Ty(goodies)-mean(Ty(goodies))).^2))   ;
T=abs(p(1))/sse;


sts=[p(1); rsq; T];

% pad anomaly with NaNs for missing data
ranom=NaN(size(Tr));
ranom(goodies)=anom;