function [kn,cw,ccw,CR,fN,fi,h1,h2]=rotaryspectra(time,blat,data,nu,varunits,timedt,color);

% Input variables should be as follows:
% time=matlab time vector (same time interval for U and V wind components)
% blat: mean latitude, to compute inertial frequency
% data=[U V]; % U and V wind components in columns
% nu=number of dof
% varunits=string specigying var units, e.g '(m/s)^2'
% timedt=string specifying timestep, e.g 'day'
% color: color for line plots, e.g. color='b'; or rgb triplet color=[0.3 0.3 0.3]

% Inertial frequency @ Buoy location
omega=2*pi/1; % in cpday!!
finertial=2*omega*sin(deg2rad(abs(blat))); % in cpday is 4.22
Tinertial=2*pi/finertial;

%% ------------------- Satellite Spectra ----------------------------------

% Make sure data is in column vector

% 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
dt=nanmean(diff(time));  % mean sampling time interval in fractions of day?
dk=1/(dt*N);   % sampling wavenumber
n=1:(N/2);
kn=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
    
    % For rotary spectra
    amU(:,j)=fft(detrend(data(m:p,1),'constant').*wind.*sqrt(8/3)); % U component
    SEU(:,j)=2*(abs(amU(2:(N/2)+1,j)).^2)/dk/N^2;
    varcheckU(j,1)=var(detrend(data(m:p),'constant'));
    varcheckU(j,2)=sum(SEU(:,j)*dk);
    
    amV(:,j)=fft(detrend(data(m:p,2),'constant').*wind.*sqrt(8/3)); % V component
    SEV(:,j)=2*(abs(amV(2:(N/2)+1,j)).^2)/dk/N^2;
    varcheckV(j,1)=var(detrend(data(m:p),'constant'));
    varcheckV(j,2)=sum(SEV(:,j)*dk);
    
    m=m+floor(nn);
    p=p+floor(nn);
end

meanSEU=mean(SEU,2); % Mean Spectral Estimate
meanSEV=mean(SEV,2); % Mean Spectral Estimate

% 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);


%% --------------- Rotary spectra -----------------------------------------

% ----------- Cross Spectra ----------------------------------------------
% ----- Co and Quadrature spectral estimates ------------------------------
[N,c]=size(amU);
for i=1:N
    for j=1:c
        % autospectra of the scalar components
        pu(i,j)=amU(i,j)*conj(amU(i,j)); % multiply Fourier coefficients
        pv(i,j)=amV(i,j)*conj(amV(i,j)); % multiply Fourier coefficients
        proUV(i,j)=amU(i,j)*conj(amV(i,j));
        
        % cross spectra
        puv(i,j)= real(amU(i,j)).*real(amV(i,j)) + imag(amU(i,j)).*imag(amV(i,j));
        quv(i,j)= -real(amU(i,j)).*imag(amV(i,j)) + real(amV(i,j)).*imag(amU(i,j));
        
    end
end

AutoSEU=2*mean(pu(2:(N/2)+1,:),2)/dk/(N^2); % nu d.o.f. Cross Spectral estimate
AutoSEV=2*mean(pv(2:(N/2)+1,:),2)/dk/(N^2);


CrossSESatUV=2*mean(proUV(2:(N/2)+1,:),2)/dk/(N^2);

CoSEUV=real(CrossSESatUV);
QSEUV=imag(CrossSESatUV);

quv=2*mean(quv(2:(N/2)+1,:),2)/dk/(N^2);
puv=2*mean(puv(2:(N/2)+1,:),2)/dk/(N^2);


% % autospectra of the scalar components
% pu = fu.*conj(fu); pv = fv.*conj(fv);
%
% % cross spectra
% puv =  real(fu).*real(fv) + imag(fu).*imag(fv);
%
% % quadrature spectra
% quv = -real(fu).*imag(fv) + real(fv).*imag(fu);
%
% rotatory components
% cw   = (pu+pv-2*quv)./8;
% ccw  = (pu+pv+2*quv)./8;

% rotatory components
cw   = (AutoSEU+AutoSEV-2*quv)./8;
ccw  = (AutoSEU+AutoSEV+2*quv)./8;


% % Rotary Coefficient
CR=(cw-ccw)./(cw+ccw); % =1 for pure rotary motion, =0 for unidirectional motion
% CR=(cwCCMP-ccwCCMP)./(cwCCMP+ccwCCMP); % =1 for pure rotary motion, =0 for unidirectional motion

%%    Rotary spectra

figure(1)
subplot(1,2,1,'FontSize',20,'CLim',[1 2]); % axes('Parent',gcf,
h1=loglog(kn,cw,'color',color,'linewidth',1.05);
% loglog(kn,cw,'color',[0.1 0.7 0.1],'linewidth',1.05); %green
%loglog(kn,ccw+cw,'k')
xlim([fi fN])
ylim([nanmin([nanmin(nanmin(cw)) nanmin(nanmin(ccw))]) nanmax([nanmax(nanmax(cw)) nanmax(nanmax(ccw))])])

set(gca,'XGrid','on','XTick',[0.1 1 2 4 10])
ylabel(['Rotary Spectra [(',varunits,')^2/cp',timedt,'day]'],'fontsize',20)
xlabel(['Frequency [cp',timedt,']'],'fontsize',20)

% Plot Confidence bar
hold on
ylocation=nanmax(cw)-0.8*nanmax(cw);
xlocation=fN-0.2*fN;
semilogy([xlocation xlocation],[err_low err_high]*ylocation,'color',color,'linewidth',2)
hold on
semilogy(finertial*ones(1,100),[0.00001:100],'--','color','k','linewidth',2)
title('Clockwise ','fontsize',26); %,num2str(nu),' d.o.f']


subplot(1,2,2,'FontSize',20,'CLim',[1 2]); %axes('Parent',gcf
h2=loglog(kn,ccw,'color',color,'linewidth',1.05);
% loglog(kn,ccw,'color',[0.1 0.7 0.1],'linewidth',1.05);
xlim([fi fN])
ylim([nanmin([nanmin(nanmin(cw)) nanmin(nanmin(ccw))]) nanmax([nanmax(nanmax(cw)) nanmax(nanmax(ccw))])])
set(gca,'XGrid','on','XTick',[0.1 1 2 4 10])
%ylabel('Rotary Spectra [(m/s)^2/cpday]','fontsize',20)
xlabel('Frequency [cpday]','fontsize',20)
%legend('ALL cw','ALL ccw','QSCAT cw', 'QSCAT ccw', 'ASCAT cw', 'ASCAT ccw',...
%    'WSAT cw', 'WSAT ccw','buoy cw', 'buoy ccw');
% 
% Plot Confidence bar
hold on
ylocation=nanmax(ccw)-0.8*nanmax(ccw);
xlocation=fN-0.2*fN;
semilogy([xlocation xlocation],[err_low err_high]*ylocation,'color',color,'linewidth',2)

% To plot inertial frequency
hold on;
semilogy(finertial*ones(1,100),[0.00001:100],'--','color','k','linewidth',2)
title('Counterclockwise','fontsize',26)% ,num2str(nu),' d.o.f'

% maximize;
set(gcf,'PaperPositionMode','auto')
% print('RotarySpectra_CW&CCW_ALL&CCMP&STRATUS.eps','-depsc2')

% %%
% 
% figure(2) % Rotary coefficient
% axes('Parent',gcf,'FontSize',20,'CLim',[1 2]);
% semilogx(kn,CR,'color',[0.4 0.4 0.4],'linewidth',2);
% %semilogx(kn,CR,'color',[0.2 0.8 0.2],'linewidth',1.05); %color for all sat
% 
% xlim([fi fN])
% ylim([-1 1])
% set(gca,'XGrid','on')
% ylabel('Rotary Coefficient','fontsize',20)
% xlabel(['Frequency [cp',timedt,']'],'fontsize',20)
% % legend('CCMP','Buoy');
% %legend('All Sat','Buoy');
% % title('STRATUS 2007-2008')
% hold on
% semilogx(kn,zeros(1,length(kn)),'k');
% 
% % maximize;
% set(gcf,'PaperPositionMode','auto')
% % print('RotaryCoeff_CCMP+STRATUS_2007-2008.eps','-depsc2')