%Find the phase of the fundamental, expressed as days past the begining of
%the year
function [day] = fund_phase_3d(data);


[NT, NN,N]=size(data);
display(['time dimension is ' num2str(NT)]) 


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

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


