function [knSat,Coh,phase]=phase_cospectra(time,data,nu,varunits,timedt);
%% Magdalena Carranza (magdac@ucar.edu)
% Code to compute Co-spectral Analyses between two concomitant timeseries;
% Plots Coherence and Phase, Quadrature and Co-spectrum
 % Co-spectrum: "defines the degree of co-oscillation for those frequency 
 % constituents of the two time series that fluctuate in-phase"
 % Quad-Spectrum: defines the degree of co-oscillation for frequency 
 % constituents of the two series that co-oscillate but are out of phase by +-90deg
 % Coherency Spectrum (Coherency): squared coherency, represents the fraction of the variance in x1 
 % ascribable to x2 through a linear relationship between x1 and x2.
 
% Reference:
% Data Analysis Methods in Physical Oceanography / Edition 2, by William J.
% Emery and Richard E. Thomson,
% pages 480-490 

%%  Input variables should be as follows:
% time=matlab time vector (same time interval for the two geophysical variables)
% data=[data1 data2]; % geophysical variables, two vectors in columns
% nu=number of dof
% varunits=string specigying var units, e.g '(m/s)^2'
% timedt=string specifying timestep, e.g 'day'

fts=14; %fontsize for plots
%% ------------------- Satellite Spectra ----------------------------------

% nu=20; %degrees of freedom
ns=nu/2; % number of segments
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
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
    
    % For rotary spectra
    am1(:,j)=fft(detrend(data(m:p,1),'constant').*wind.*sqrt(8/3)); % U component
    SE1(:,j)=2*(abs(am1(2:(N/2)+1,j)).^2)/dk/N^2;
    varcheck1(j,1)=var(detrend(data(m:p,1),'constant'));
    varcheck1(j,2)=sum(SE1(:,j)*dk);
    
    am2(:,j)=fft(detrend(data(m:p,2),'constant').*wind.*sqrt(8/3)); % V component
    SE2(:,j)=2*(abs(am2(2:(N/2)+1,j)).^2)/dk/N^2;
    varcheck2(j,1)=var(detrend(data(m:p,2),'constant'));
    varcheck2(j,2)=sum(SE2(:,j)*dk);
    
    m=m+floor(nn);
    p=p+floor(nn);
end


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

meanSE1=mean(SE1,2); % Mean Spectral Estimate
meanSE2=mean(SE2,2); % Mean Spectral Estimate

%% -------------- COHERENCE ANALYSIS --------------------------------------
% -------------------------------------------------------------------------

% ----- Co and Quadrature spectral estimates at nu d.o.f-----------
[N,c]=size(am1);
for i=1:N
    for j=1:c
        proS(i,j)=am1(i,j)*conj(am2(i,j)); % multiply Fourier coefficients
    end
end
CrossSE=2*mean(proS(2:(N/2)+1,:),2)/dk/(N^2); % nu d.o.f. Cross Spectral estimate

% %% Plot Co-spectrum and Quadrature-spectrum
% figure(1)
% % subplot(2,1,1,'Parent',gcf,'FontSize',fts,'CLim',[1 2]);
% semilogx(knSat,real(CrossSE),'color','k','linewidth',1.5); % Co-spectrum
% grid on
% % title('Cross Spectral Analysis','FontSize',fts)
% ylabel(['Co-spectrum [(m/s)^2/cp',timedt,']'],'FontSize',fts)
% xlabel(['Frequency [cp',timedt,']'],'FontSize',fts)
% % xlim([0.03 10])
% xlim([fi fN])
% 
% figure(2)
% % subplot(2,1,2,'Parent',gcf,'FontSize',fts,'CLim',[1 2]);
% semilogx(knSat,imag(CrossSE),'k','linewidth',1.5); % Quadrature spectrum
% hold on
% % xlim([0.03 10])
% grid on
% ylabel(['Quadrature-spectrum [(m/s)^2/cp',timedt,']'],'FontSize',fts)
% xlabel('Frequency [cpday]','FontSize',fts)
% xlim([fi fN])
% 
% % maximize;
% % set(gcf,'PaperPositionMode','auto')
% % print('CrossSpectralAnalysis_CCMPvsBuoy.eps','-depsc2')

%%  Coherence and Phase estimates
%
% Coherence
for j=1:length(CrossSE)
    Coh(j)=abs(CrossSE(j))/sqrt(meanSE1(j)*meanSE2(j));
end

% Phase
phase=atan2(imag(CrossSE),real(CrossSE));

% Confidence
edof=ns; % number of segments (Biltof and Pardyjak, 2009)
deltaCoh=sqrt(1-alpha^(1/(edof-1))); % Thompson (1979)
deltaCohF=2*finv(0.95,10,10)/(edof-2+2*finv(0.95,10,10)); %3.81 for 15 edof according to Biltof

% %% Plot Coherence
% figure(3)
% % subplot(2,1,1,'Parent',gcf,'FontSize',18,'CLim',[1 2]);
% semilogx(knSat,Coh,'color',[0.2 0.2 0.2],'linewidth',2);
% % semilogx(knSat,CohB,'color',[0.1 0.7 0.1],'linewidth',1.5); %green
% hold on
% % semilogx(knSat,CohS,'color','k','linewidth',1.5)
% % hold on
% % plot(knSat,deltaCoh*ones(length(knSat),1),'--','color','k');
% % hold on
% plot(knSat,deltaCohF*ones(length(knSat),1),'--','color',[0.2 0.2 0.2]);
% hold on
% 
% % xlim([0.03 10])
% xlim([fi fN])
% 
% % title('Coherence and Phase Estimates','FontSize',18)
% ylabel('Coherence','FontSize',18)
% xlabel(['Frequency [cp',timedt,']'],'FontSize',18)
% set(gca,'ylim',[0 1.01],'ytick',[0:0.2:1])
% hold on
% grid on

% Plot Phase
figure(4)
% subplot(2,1,2,'Parent',gcf,'FontSize',18,'CLim',[1 2]);
semilogx(knSat,rad2deg(phase),'color',[0.2 0.2 0.2],'linewidth',2);
% hold on
% semilogx(knSat,rad2deg(phaseB),'k','linewidth',1.5);
% xlim([0.03 10])
xlim([fi fN])

ylabel('Phase [degrees]','FontSize',18)
xlabel(['Frequency [cp',timedt,']'],'FontSize',18)
set(gca,'ylim',[-180 180],'ytick',[-180:30:180])
grid on

% maximize;
% set(gcf,'PaperPositionMode','auto')
% print('CoherenceAnalysis_CCMPvsSTRATUS.eps','-depsc2')

%% ---------- Confidence interval for Coherence --------------------------
% Normalized error for Squared Coherence:

e=sqrt(2).*(1-(Coh.^2))./(Coh.*sqrt(ns));

% Normalized error for Coherence (eq. 8.10 B&P)
e=e./2;

% Biased of squared coherence:

b=((1-Coh.^2).^2)./ns;

% Biased of Coherence:
uno=ones(1,length(b));
b=1/ns; % systematic error for nule coherence
b=sqrt(b);
b=uno.*b;

% Confidence limits for Coherence at 68% confidence level:
% (approach valid if random error is small (<0.2 according to B&P))
% Taking as a treshhold e<0.5 (minimum(e) in our data is 0.22!):

for j=1:length(e)
    if e(j)<=0.5
        llim(j)=(Coh(j))/(1+e(j));
        ulim(j)=(Coh(j))/(1-e(j));
    else
        llim(j)=NaN;
        ulim(j)=NaN;
    end
end


% figure(18)
% % Plot Coherence with confidence limits and biased error for zero coherence:
% subplot(211)
% semilogx(kn,(CohS),'color',[0.1 0.7 0.1],'linewidth',1.5)
% hold on
% title('Coherence and Phase','fontsize',12)
% ylabel('Coherence','fontsize',12)
% xlabel('Frequency [cps]','fontsize',12)
% hold on
% grid on
% semilogx(kn,llim,'k--')
% hold on
% semilogx(kn,ulim,'k--')
% % hold on
% % semilogx(kn,b,'k')
% %legend('Coherence','Upper limit 67% CL','Lower limit 67% CL','Bias Error')
%
%
% % Plot phase when Coherence is significantly different from zero
% % (i.e when lower limit is above the biased error value)
% %
%
% for j=1:length(llim)
%     if llim(j)>b
%     cphase(j)=phase(j); % significant phase
%     else
%     cphase(j)=NaN;
%     end
% end
%
% subplot(212)
%
% semilogx(kn,rad2deg(cphase),'r*','markersize',10)
% hold on
% semilogx(kn,rad2deg(phase),'r--','linewidth',0.2)
% legend('significant phase')
% ylabel('Phase [degrees]','fontsize',12)
% xlabel('Frequency [cps]','fontsize',12)
% set(gca,'ylim',[-180 180],'ytick',[-180:30:180])
% grid on
