function [nameu,fu,tidecon,xout]=tide(xin,dt,stime,lat,ray);
% [nameu,fu,tidecon,xout]=tide(xin,dt,stime,lat,ray);
% computes the tidal analysis of the complex time series xin=x+i*y.
%
%  INPUT:  xin=time series                            required
%          dt=sample time (decimal hours)             optional
%          stime=start time  [year,month,day,hour]    optional
%          lat=latitude (decimal degrees)             optional
%          ray=rayleigh criteria (default ray=1)      optional
%
%  OUTPUT: nameu=list of constituents used
%          fu=frequency of tidal constituents (cycles/hr)
%          tidecon=[fmaj,emaj,fmin,emin,finc,einc,pha,epha] for vector xin
%                 =[fmaj,emaj,pha,epha] for scalar (real) xin
%            fmaj,fmin - constituent major and minor axes (same units as xin)       
%            emaj,emin - 95% confidence intervals for fmaj,fmin
%            finc - ellipse orientations (degrees)
%            einc - 95% confidence intervals for finc
%            pha - constituent phases (degrees relative to Greenwich)
%            epha - 95% confidence intervals for pha
%          xout=tidal prediction

% CALLS:  constituents.m - loads 146 constituents, frequency and
%                          comparison constituent
%         find_const.m - finds match of constituent name to list of
%                        constituents
%         vuf.m - computes correction for phase relative to Greenwich
%                 and nodal amplitude and phase corrections. only called 
%                 if stime and lat are given. vuf.m additionally calls 
%                 opnvuf.m, astron.m, julian.m, and find_const.m.
%         errcoef - computes errors in sin,cos coefficients.
%         errell - computes errors in ellipse parameters.
%
% DATA:   constituents.dat - contains list of constituent names, 
%                            frequencies and comparison constituents 
%                            (total presently 146).
%         following three called by opnvuf.m
%         doodson.dat - Doodson numbers (6) for 45 major constituents
%         satellite.dat - satellite constituent information for major
%                         constituents
%         shallow.dat - shallow water constituent information 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1/23/99
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (nargin<5)
 ray=1; % rayleigh criterion set to 1 if not specified
 disp(['   rayleigh criterion set to 1.0   '])
end

% load tidal constituents (presently 146)
% variables: name, freq, kmpr (raleigh comparison constituents)
[name,freq,kmpr]=constituents;

xin=xin(:); % makes xin a column vector
nobs=length(xin);
nobsu=nobs-mod(nobs-1,2);% makes series odd to give a center point
if (nargin==1); dt=1; end; % sample time default is 1 hour

% find useable constituents (nameu) and frequencies (fu) based
% on rayleigh criterion

% coefficients without comparison constituent are not used
% in the present configuration

kbl=['    '];
  for k=1:length(freq);
   if (strcmp(kmpr(k,:),kbl))
    df(k)=0;
   else
    j1=find_const(kmpr(k,:),name);
    df(k)=abs(freq(j1)-freq(k));
   end
  end
ju=find(df.*dt.*nobsu>=ray);
mu=length(ju);  
nameu=name(ju,:);nameu=nameu(2:mu,:);
fu=freq(ju);fu=fu(2:mu);
mu=mu-1; % removes z0 (freq=0) constituent from arrays
disp(['   number of constituents used: ',int2str(mu)])

% check if neighboring chosen constituents violate the rayleigh criteria.
jck=find(diff(fu).*dt.*nobsu<ray);
if (length(jck)>0)
   disp(['  Warning! Following constituent pairs violate Rayleigh criterion']);
   for ick=1:length(jck);
   disp(['     ',nameu(jck(ick),:),'  ',nameu(jck(ick)+1,:)]);
   end;
end

% least squares fit
% solve for least squares fit coefficients for cos and sin, mean, trend, 
% and residual time series for the real (x) and imaginary (y) components of
% the data time series xin. fu are the tidal frequencies used in the fit.
% cx are cosine coefficients, sx sine coefficients, x0 mean, dx0 trend, 
% xp model prediction, and xres the residual time series for x.  similar 
% calculations done for y is xin is complex, otherwise cy and sy set to zero.
% (note that xin may have gaps or missing data denoted by NaN's; these are
% removed before the fit.

x=real(xin); y=imag(xin);
jp=1:nobs;
t=dt.*(jp-ceil(nobs./2));
one=ones(length(t),1);
tr=2.*t'./(dt.*nobs); % normalized by record length
cs=cos(2.*pi.*fu*t).';
sn=sin(2.*pi.*fu*t).';
tc=[one,tr,cs,sn];
jp=1:nobsu;
j=jp(find(~isnan(x(jp)+y(jp)))); % removes NaN's from analysis 
ngood=length(j); % number of good data used in fit
disp(['   number of good data used: ',int2str(ngood)])
coef=tc(j,:)\x(j);
x0=coef(1); dx0=coef(2);
cx=coef(3:mu+2);
sx=coef(mu+3:2.*mu+2);
xp=tc*coef;
xres=x-xp;
varx=cov(x(j));varxp=cov(xp(j));varxr=cov(xres(j));
disp(['   percent of x var removed by fit: ',num2str(100.*(1-varxr./varx)),' %'])

% if xin is complex, do imaginary time series.
svec=sum(abs(y(j)));
if (svec>0)
  coef=tc(j,:)\y(j);
  y0=coef(1);dy0=coef(2);
  cy=coef(3:mu+2);
  sy=coef(mu+3:2.*mu+2);
  yp=tc*coef;
  yres=y-yp;
  vary=cov(y(j)); varyp=cov(yp(j)); varyr=cov(yres(j));
  xout=xp+i.*yp;
  disp(['   percent of y var removed by fit: ',num2str(100.*(1-varyr./vary))])
 else
  cy=zeros(mu,1);sy=cy;
  xout=xp;
end

% estimate uncertainities in cos, sin coefficients for x and y.
jg=find(~isnan(xres));
jt=[min(jg):max(jg)];
xg=interp1(jg,xres(jg),jt); % to fill small gaps within series
reclen=dt.*length(jt)./24;
[ercx,ersx]=errcoef(xg,fu,dt);
if (svec>0);
   yg=interp1(jg,yres(jg),jt);
   [ercy,ersy]=errcoef(yg,fu,dt);   
  else
   ercy=zeros(size(ercx));
   ersy=ercy;
end

% estimate uncertainities in ellipse parameters.
[emaj,emin,einc,epha]=errell(cx,sx,cy,sy,ercx,ersx,ercy,ersy);

% convert to polar coordinates.
t1=cx+sy;
t2=cy-sx;
ap=sqrt(t1.*t1+t2.*t2)./2.0;
epsp=atan2(t2,t1)./(2.0.*pi);
t1=cx-sy;
t2=cy+sx;
am=sqrt(t1.*t1+t2.*t2)./2.0;
epsm=atan2(t2,t1)./(2.0.*pi);

% apply nodal corrections and calculate phase relative to Greenwich.
if (nargin>=3) 
  % for some reason foreman code limits latitude to be greater 5 degrees.
  if (abs(lat)<5); lat=sign(lat).*5; end
  disp(['nodal corrections applied and phase is relative to Greenwich']);
  [v,u,f]=vuf(stime,nobsu,nameu,lat);
  vu=v+u; % total phase correction
  nodcor=['nodal corrections applied to amplitude and phase relative to center time'];  
  disp(['   ',nodcor]);
else
  f=1;vu=0;
  nodcor=['nodal corrections not applied'];  
  disp(['   ',nodcor]);
end

aap=ap./f;
aam=am./f;
fmaj=aap+aam; % major axis
fmin=aap-aam; % minor axis
gp=mod((-epsp+vu).*360.0,360); % positive phase in degrees 
gm=mod((epsm+vu).*360.0,360);  % negative phase in degrees 
finc=mod(360.*(epsp+epsm),360)./2.0; % ellipse inclination in degrees
pha=mod((gp+finc),360); % phase in degrees
snr=(fmaj.^2)./(emaj.^2); % signal to noise ratio

% output results

filen=input('   input file name (hit return to skip): ','s');
sav=input('   enter y to save results to file, skip to display on screen','s');
if strcmpi(sav,'y'), fid=fopen([filen,'.txt'],'w'); else fid=1; end
fprintf(fid,'\n%s\n',['file name: ',filen]);
fprintf(fid,'%s\n',['date: ',date]);
fprintf(fid,'%s\n',['nobs = ',int2str(nobs),'   ngood = ',int2str(ngood),...
      '   record length (days) = ', num2str(reclen)]);
if (nargin>2); fprintf(fid,'%s\n',['start time: ',int2str(stime)]); end
fprintf(fid,'%s\n',['rayleigh criterion = ',num2str(ray)]);
fprintf(fid,'%s\n',nodcor);
fprintf(fid,'\n%s\n',['     coefficients from least squares fit of x']);
fprintf(fid,'\n%s\n',['tide    freq       cosx      ercx      sinx      ersx']);
for k=1:length(fu);
  fprintf(fid,'%s  %8.5f  %8.3f  %8.3f  %8.3f  %8.3f\n',nameu(k,:),fu(k),cx(k),ercx(k),sx(k),ersx(k));
end
fprintf(fid,'\n%s\n',['x0= ',num2str(x0),'      x trend= ',num2str(dx0)]);
fprintf(fid,'%s\n',['var(x)= ',num2str(varx),'   var(xp)= ',num2str(varxp),'   var(xres)= ',num2str(varxr)]);
fprintf(fid,'%s\n',['percent var predicted= ',num2str(100.*varxp./varx),' %']);

if (svec==0)
  tidecon=[fmaj,emaj,pha,epha];
  fprintf(fid,'\n%s\n',['     tidal amplitude and phase with 95% error estimates']);
  fprintf(fid,'\n%s\n',['tide   freq       amp    amp_err    pha    pha_err     snr']);
  for k=1:length(fu);
  fprintf(fid,'%s %8.5f %8.3f %8.3f %8.2f %8.2f %8.2f\n',nameu(k,:),fu(k),fmaj(k),emaj(k),pha(k),epha(k),snr(k));
  end
end

if (svec>0)
  tidecon=[fmaj,emja,fmin,emin,finc,einc,pha,epha];
  fprintf(fid,'\n%s\n',['     coefficients from least squares fit of y']);
  fprintf(fid,'\n%s\n',['tide    freq       cosy      ercy      siny      ersy']);
  for k=1:length(fu);
  fprintf(fid,'%s  %8.5f  %8.3f  %8.3f  %8.3f  %8.3f\n',nameu(k,:),fu(k),cy(k),ercy(k),sy(k),ersy(k));
  end
  fprintf(fid,'\n%s\n',['y0= ',num2str(y0),'      y trend= ',num2str(dy0)]);
  fprintf(fid,'%s\n',['var(y)= ',num2str(vary),'    var(yp)= ',num2str(varyp),'  var(yres)= ',num2str(varyr)]);
  fprintf(fid,'%s\n',['percent var predicted= ',num2str(100.*varyp./vary),' %']);
  fprintf(fid,'\n%s\n',['ellipse parameters with 95% error estimates']);
  fprintf(fid,'\n%s\n',['tide   freq      major  emaj    minor  emin     inc    einc     pha    epha      snr']);
  for k=1:length(fu);
  fprintf(fid,'%s %8.5f %7.2f %6.2f %7.2f %6.2f %8.2f %6.2f %8.2f %6.2f %8.2f\n',...
      nameu(k,:),fu(k),fmaj(k),emaj(k),fmin(k),emin(k),finc(k),einc(k),pha(k),epha(k),snr(k));
  end
  fprintf(fid,'\n%s\n',['total var= ',num2str(varx+vary),'   pred var= ',num2str(varxp+varyp)]);
  fprintf(fid,'%s\n',['percent total var predicted= ',num2str(100.*(varxp+varyp)./(varx+vary)),' %']);
end

if strcmpi(sav,'y'), st=fclose(fid); end


