%Find the phase of the fundamental, expressed as days past the begining of
%the year
function [day] = fund_phase(data);
[NN, NT]=size(data);
if NT == 1;
    data=data';
    [NN, NT]=size(data);
end

for nn=1:NN;
H=fft(data(nn,:));
coef(nn,1)=abs(H(2)/(NT/2));
phase(nn)=angle(H(2));
end

%convert the phase to days past the begining of the year
offset=1; %day offset, 15 if the data is monthly average,  1 if it's for the first of the month
daysperyear=365.25;
for nn=1:NN;
if phase(nn)<=0;
    day(nn)=-phase(nn)*daysperyear/(2*pi) +offset;
else phase(nn) > 0;
    day(nn) = (2*pi - phase(nn))*daysperyear/(2*pi) +offset;
end
end
%day=phase;
%  if phase<=0 & phase >-pi;
%      day=-phase*daysperyear/(2*pi) +offset;
%  elseif phase > 0;
%      day = (2*pi - phase)*daysperyear/(2*pi) +offset;
%  else
%      day=-(phase+2*pi)*daysperyear/(2*pi) +offset;
%  end


