function y=mk_tide(t,coef,n)
% y=mk_tide(t,coef,n) computes a simple tidal time series for the 
% time period t for the 5 major tidal constituents whose amplitude 
% is given in coef [O1,K1,N2,M2,S2]. guassian noise is added to the 
% series with std = n. the mulitaper spectrum of y is plotted.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1/16/99
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

f=[.0387306544;.0417807462;.0789992488;.0805114007;.08333333333];
t=t(:)';coef=coef(:);
b=(2.*pi.*f*t).';
y=cos(b)*coef;
if (nargin>2)
   noise=n.*randn(1,length(t));
   y=y(:)'+noise;
else
   y=y(:)';
end

ny=length(y);
[Px,Pxc,fx]=pmtm(y,2,ny,1);
loglog(fx,Px,'-b',fx,Pxc(:,1),'-r',fx,Pxc(:,2),'-r')
xlabel('frequency (cph)')
ylabel('power density (units^2/cph)')
title(['multitaper estimate of power density spectrum with nw=2  ',date])
% hold on
% [Px,Fx]=psd(y,ny,1,boxcar(ny),0);
% loglog(Fx,Px,'-k')
pause
