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