function [ts_out,f, coef, phase] = harmonics(t, data, hnum )

%does the harmonic anaylsis of the data set, retaining hnum harmonics

% OUTPUTS

% f is the frequency of each harmonic
% coef = coefficient of each harmonic from 0,1, ... hnum 
% phase = phase of each harmonic to reconstruct the data, add the phase to the cosine of the 
%               root of the power PLUS the phase  (a sin in the data = cos(t-pi/2)) 
%  !!!!!!!!!!!!!!!!!!!!! t is the time from start of time series  !!!!!

% INPUTS
% data = a vector of the data on an evenly spaced domain

SPD=(length(t)-1)/(max(t)-min(t)); %samples per day 
N=length(data); %number of points in domain

%do the fft
Yfft=fft(data);
f=SPD*(0:1:N/2)/N; %definition of frequencies, should be multiplied by samples per day
NF=length(f);

coe(1:NF)=Yfft(1:NF)/(N/2);
coef(1:NF)=abs(Yfft(1:NF)/(N/2));
coef(1)=mean(data)
phase(1:NF)=angle(coe); %defined such that the original data is cos(t + phase) --- i.e. the sin would have phase -pi/2
phase(1)=0;

%recompose the original data

ts_out=coef(1)*ones(1,length(t));
for n=2:hnum+1;
    ts_out=ts_out+(coef(n)*cos(2*pi*f(n)*(t-min(t))+phase(n))); %add in each component in the time domain
end

figure; clf; LW=2;
plot(t,data,'r', 'LineWidth', LW)
hold on; plot(t, ts_out, 'g', 'LineWidth', LW)
legend('original', 'harmonic out')
set(gca, 'FontSize', 17)



