function [knSat,meanSESat,dk,varcheckSat,h,SESat]=spectra(time,data,nu,varunits,timedt,color);
%% Magdalena Carranza (magdac@ucar.edu)
% Code to compute and plot Spectral Density (in units of var^2/cpday)
% Reference:
% Data Analysis Methods in Physical Oceanography / Edition 2, by William J. Emery and Richard E. Thomson

% Input variables should be as follows:
% time=matlab time vector of obs; or dt in fraction of day
% data=geophysical variable vector
% nu=number of degrees of freedom (dof)
% varunits=string specifying var units for ploting, e.g '(m/s)^2'
% timedt=string specifying timestep/units, e.g 'day'
% color= for plots
% NOTE: data cannot have NaNs
% NOTE: if working with profile data in columns, nu should be twice the
% numer of profiles.. such that ns=nu/2 is the number of segments (i.e.
% profiles)
%%
fts=14; % fontsize for plots

%% ------------------- Spectral Analysis ----------------------------------

% Arrange data in column vector
data=data(:);
time=time(:);

% nu=20; % degrees of freedom
ns=nu/2; % number of segments
% nA=Adata(:,1)); nQ=length(Qdata(:,1)); nW=length(Wdata(:,1)); % use length of shortest timeseries
nn=length(data)/ns; % min([nA nQ nW])/ns;

m=1;
p=floor(nn);
N=p; % number of data points in the time series of each individual spectra

if length(time)>1 % input is time vector
    dt=nanmean(diff(time));  % compute mean sampling time interval
elseif length(time)==1
    dt=time; % time = dt
end

dk=1/(dt*N);   % sampling wavenumber
n=1:(N/2);
knSat=n./((N-1)*dt);

% Nyquist frequency
fN=1/(2*dt);

% Lowest resolved frequency
fi=1/(N*dt);

% Hanning window
wind=hann(N);

for j=1:ns % compute spectra for each segment
    
    % Compute spectral coefficients
    amSat(:,j)=fft(detrend(data(m:p),'constant').*wind.*sqrt(8/3)); % 
    SESat(:,j)=2.*(abs(amSat(2:(N/2)+1,j)).^2)/dk/N^2;
    
    % Check that the sum of the amplitudes is equal to the total variance
    varcheckSat(j,1)=var(detrend(data(m:p),'constant'));
    varcheckSat(j,2)=sum(SESat(:,j)*dk);
    
    m=m+floor(nn);
    p=p+floor(nn);
end

% Compute Spectral Estimates (averaging over segments)
meanSESat=mean(SESat,2); % Mean Spectral Estimate

% Compute Spectral error bar for nu d.0.f.
alpha=0.05; % for 95% confidence
% Serr_low=nu/chi2inv(alpha/2,nu);
% Serr_high=nu/chi2inv(1-alpha/2,nu);

% Plot Spectra in a loglog scale
h=loglog(knSat,meanSESat,'color',color,'linewidth',1.5);
xlim([fi fN])
% ylim([nanmin(nanmin(SESat)) nanmax(nanmax(SESat))])

ylabel(['Spectral Density [(',varunits,')^2/cp',timedt,']'],'fontsize',fts)
xlabel(['Frequency [cp',timedt,']'],'fontsize',fts)
set(gca,'XGrid','on','fontsize',fts)
hold on

% % Compute Spectral error bar for nu d.0.f.
% alpha=0.05; % for 95% confidence
% err_low=nu/chi2inv(alpha/2,nu);
% err_high=nu/chi2inv(1-alpha/2,nu);
% 
% % Plot Confidence bar
% ylocation=nanmax(meanSESat)-0.8*nanmax(meanSESat); % set's where to place the error bar in the plot
% 
% a=fi+1/10; b=fN-0.8*fN;
% rx = a + (b-a).*rand(1,1); %random number to move error bar x-location
% xlocation=fN-0.8*fN+rx; %sets where to place the error bar in the xaxis
% xlocation=3*10^(-2);
% 
% 
% semilogy([xlocation xlocation],[err_low err_high]*ylocation,'color',color,'linewidth',2)

hold on
% maximize;
% title(['All satellite @ STRATUS VIII Location, ',num2str(nu),' d.o.f '],'fontsize',14)
% % legend with mean values and stds:
% legend(['STRATUS     m=',num2str(nanmean(WSPD10is),'%6.2f'),', std=',num2str(nanstd(WSPD10is),'%6.2f')]);
% legend1=legend('All (QSCAT + WindSat + ASCAT) ','u10 (COARE v3.0)','wspd(3.3m)','Location','southwest');
