function [f,Cxy,C,px1,px2,py,px1y,px2y,px1x2] = multcoh(varargin)

%
%  [f,Cxy,C,px1,px2,py,px1y,px2y,px1x2] = multcoh(L,x1,x2,y,eps)
%  
%  Calculated multiple coherence between x1, x2 and y. The assumption is
%  that x1 and x2 are simultaneously related to some variability in y (i.e.,
%  perhaps they are two independent forcings on y). See Priestly (1981)
%  for a discussion on the theory.
%
%   f    - frequencies
%   Cxy  - coherence (not squared) between x1, x2 and y
%   C    - 5% significance level for coherence
%   p(.) - power spectrum of (.)
%   p(..)- cross sepctrum of (. and .)
%   eps  - (optional) fraction (between 0 and 1) which is multiplied by the max
%          of the spectras of x1 and x2 in the calculation of Cxy. Numbers of
%          order 0.02 allow for elimination of errors due to spectra which go to
%          zero at high frequencies.
%  
%  The spectra are averaged with an L-point triangularly weighted window.
%  If L=-1 then the window will be defined as 10% of the length of the
%  series.

% grab inputs, nargs = 4 is eps = 0, nargs = 5 defines eps from input
[cax,args,nargs] = axescheck(varargin{:});
L  = args{1};
x1 = args{2};
x2 = args{3};
y  = args{4};
if nargs == 5
    eps = args{5};
else
    eps = 0;
end

% get spectra for x1, x2, y
[f,px1,py,px1y,cx1y,C] = spectra(L,x1,y);
[f,px2,py,px2y,cx2y,C] = spectra(L,x2,y);
[f,tmp,tmp,px1x2,tmp,tmp,tmp,tmp,tmp,Lh] = spectra(L,x1,x2);

% combine these to get multiple coherence
Cxy = sqrt( (1./(py.*( (px1+eps*max(px1)) .* (px2+eps*max(px2)) -abs(px1x2).^2))).*( (px2+eps*max(px2)).*abs(px1y).^2 + (px1+eps*max(px1)).*abs(px2y).^2 - px1y.*conj(px2y).*conj(px1x2) - px2y.*conj(px1y).*px1x2 ) );

% statistical significance
alpha = 0.05;
degFreed = 2*Lh*length(y)/(2*length(py));
F = finv(1-alpha,2*2,2*(degFreed-2));
C = F./( 0.5*(degFreed-2) + F);
C = sqrt(C);
