function [ercx,ersx]=errcoef(xres,fu,dt)
% [ercx,ersx]=errcoef(xres,fu,dt)computes errors in the cos and sin 
% coefficients of the least square fit to data time series x, based 
% on multitaper spectral analysis of the residual series xres. 
%
%  INPUT:  xres=residual time series for x     
%          fu=tidal frequencies used in fit (cycles/hr)
%          dt=sample time (decimal hours)
%          
%  OUTPUT: ercx=error in cos coefficients
%          ersx=error in sin coefficients

% mean values of the noise spectrum are computed for the following 
% 8 frequency bands defined by their center frequency and band width:
% M0 +.1 cpd; M1 +-.2cpd; M2 +-.2cpd; M3 +-.2cpd; M4 +-.2cpd; M5 +-.2cpd;
% M6 +-.21cpd; M7 (.26-.29cpd); and M8 (.30-.50cpd). noise variance within
% each band split equally to cos and sin coefficients.  assumes errors in
% coefficients have normal distributions, so that 95% confidence limits
% are given by +- 1.96 std.  also assumes fu(1)>0.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1/16/99; 1/20/99
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

xres=xres(:).'; % makes xres a row vector
fu=fu(:).'; % makes fu a row vector
mu=length(fu);
[Px,fx]=pmtm(xres,4,length(xres),1/dt);
Px=2.*Px; % double power so all in positive frequencies
df=fx(3)-fx(2);
fband=[.00010 .00417;.03192 .04859;.07218 .08884;.11243 .12910;.15269 .16936;...
       .19295 .20961;.23320 .25100;.260 .500];
Pxave=zeros(8,1);
for k=1:8
   jband=find(fx>=fband(k,1) & fx<=fband(k,2));
   Pxave(k)=mean(Px(jband));
end
erx=1.96.*sqrt(Pxave.*df); % assumes normal distrubution
ercx=zeros(mu,1);
for k1=1:8;
   k=find(fu>=fband(k1,1) & fu<=fband(k1,2));
   ercx(k)=erx(k1);
end
ersx=ercx;

