function [kk,ph,kkl,kku,phl,phu,df]=pcoh(I11,I22,I12,kw,pconf)
%
% [kk,ph,kkl,kku,phl,phu] = pcoh(I11,I22,I12,kw,pconf)
%
% Computes coherence and phase estimates from spectral estimates 
%          obtained by smoothing auto and cross periodograms calculated by
%          CROSSP.M 
%          Confidence interval analysis is as given in Jenkins & Watts,
%          Spectral Analysis and its Applications, Holden-Day, 1968, eq.s 
%          9.2.22-9.2.25.
%          (refer to Koopmans' book for the degrees of freedonm calculation)
% 
% kk,kkl,kku  -are squared coherence and lower, upper confidence bands
% ph,phl,phu  -are phase spectrum and conf. bands in degrees (+ve => x2 leads)
% I11,I22,I12 -are periodograms as calculated by CROSSP.M            
% kw          -total (odd) width of (symmetric) spectral window to use
% pconf       -percent conf limit (90 or 95)
%
% NOTE : - CROSSP.M required for periodogram computation.
%        - Phase straightening can be done with STRAIGHT.M
%
% This routine requires access to periodogram smoothing macros 
% csmooth.m, tsmooth.m, rsmooth.m.
% MIM 25 May 1992
% MIM 23.9.97 pconf added as an optional i/p variable 
%                   = % confidence.....can only be 90 or 95, 
%                   defaults to 90 if not set or if not set to 95
%
% JLW: I only have the triangular wind tsmooth so remove the other options


% (1) smooth periodograms for spectral estimates

[f11,df] = tsmooth(I11,kw);
f22 = tsmooth(I22,kw);
f12 = tsmooth(I12,kw);
df = df(1); % No. deg. freedom in Chis-sq. approx. to dist. of spectra

% (2) Coherence estimation

phsfac = 180./pi;
ph = angle(f12);
kk = (f12.*conj(f12)) ./ (f11.*f22);

% (3) Compute 90 or 95% confidence intervals 
%  - Jenkins & Watts[1968], p380 recipe will be used.....

if nargin > 4
  switch pconf
    case 90
      pcon = 90;
      z1 = 1.645;
    case 95
      pcon = 95;
      z1 = 1.96;
    otherwise
      pcon = 90;
      z1 = 1.645;
  end
else
  pcon = 90.0;
  z1 = 1.645; % pt from gaussian distrib
end;

yunc = z1/sqrt(df); % (J/W equation 9.2.23)
disp([num2str(pcon),' c.i. for arctanh(k12) given by  +/-',num2str(yunc)]);

k12 = sqrt(kk);

% coherence uncertainty via Arctanh(k) as in Jenkins & Watts
ok = find(k12~=1);
y12 = -999*(1+0*k12); % trap for k12 = 1 
y12(ok) = .5*log((1+k12(ok))./(1-k12(ok)));
y12u= y12 + yunc;
y12l= y12 - yunc;
k12u=k12;
k12l=k12;
k12u(ok) = tanh(y12u(ok));
bad=find(k12u<0);
k12u(bad)=0*(bad);
k12l(ok) = tanh(y12l(ok));
bad=find(k12l<0);
k12l(bad)=0*(bad);
kku = k12u.*k12u;
kkl= k12l.*k12l;

cph = cos(ph);
temp= 0*(cph);
ok2=find(cph~=0&k12~=1);
temp(ok2)=(1+0*ok2)./cph(ok2);
sec4 = temp.*temp.*temp.*temp;
tpesq= 0*(temp); 
tpesq(ok2)=(sec4(ok2)/df).*( (1+0*ok2)./(kk(ok2)) -1);
tpe  = 0*(ph);
ok3=find(tpesq>0);
tpe(ok3) = sqrt(tpesq(ok3));

tpe=tpe*z1;  %<--------23.9.97 MIM this is an error fixed
 %  previously we had effective z1=1 giving essentially 1 sigma =
 %  68% confidence interval on phase!!!! this bug has been here and in the 
 %  previous fortran code for 15 bloody years!!!!  Fortunately, coherence has 
 %  had the correct uncertainty calculation all along
 %  Note, z1 is Jenkins & watts n(1-a/2) in their equation 9.2.25

%%%USE DERIVATIVE OF PHASE WRT TAN PHASE TO COMPUTE VARIATION IN PHASE:
tanp = tan(ph);
dpdtp = (1+0*tanp)./(1+tanp.*tanp);
pe = dpdtp .* tpe;
phu = (ph + pe)*180/pi;
phl = (ph - pe)*180/pi;
ph  = ph*180/pi;
mask=find(kk<.2);
phu(mask)=ph(mask);
phl(mask)=ph(mask);