\documentclass{article}
%\usepackage{doublespace,jpo,psfig,amssymb}
\usepackage{jpo,psfig,amsmath,moreverb}
\oddsidemargin -0in %left margin on odd-numbered pages is 1in + 0in
\topmargin -.5in %top margin is 1in -0.5in
\textwidth 6.375in %width of text
\textheight 9in % size of page

\begin{document}

\title{Classical Tidal Harmonic Analysis with Error Analysis in MATLAB using T\_TIDE}
 
\author{Rich Pawlowicz\thanks{Dept. Of Earth and Ocean Science, 
University of British Columbia
6270 University Blvd., Vancouver, B.C., Canada, V6T 1Z4. email: rich@ocgy.ubc.ca}
 \\ 
{\small Department of Earth and Ocean Sciences\/}, \\ 
{\small University of British Columbia\/}, \\
{\small Vancouver, Canada, V6T 1Z4} \\
\\Bob Beardsley, Steve Lentz,\\
{\small Woods Hole Oceanographic Institution\/}, \\
{\small Woods Hole, MA, 02543, USA}\\
\\
}
\date{Submitted to: Computers and Geosciences}
 
\maketitle 
\newpage
\begin{abstract} 


A standard part of any oceanic pressure gauge or current meter analysis is the
separation of tidal from non-tidal components of the signal. The tidal signal can
either be discarded, or its characteristics described in some fashion useful for
further analysis. Although tidal
signals can be removed by standard high or band-pass filtering techniques, their relatively
deterministic character and large amplitude make special techniques more
effective. In classical harmonic analysis the tidal signal is modelled as the
sum of a finite set of sinusoids at specific frequencies related to
astronomical parameters. 
A set of programs has been written in MATLAB to a) perform 
classical harmonic analysis for periods of about one year or shorter, b) account 
for (some) unresolved consituents using nodal corrections, and c) compute 
confidence intervals for the analyzed components.


\end{abstract}

%%\section*{Introduction}

As the earth rotates on its axis spatially varying gravitational forces from the moon and the 
sun act on the ocean, generating a forced elevation and current response primarily (but not solely) at
diurnal and semi-diurnal frequencies. Body forces act directly on deep oceanic
waters. Tidal effects in coastal regions are not directly forced by these astronomical
forces. Instead they arise as a side-effect of deep oceanic variability, 
propagating through shallower coastal waters as a wave.  
In a typical oceanic time series tidal 
variability is often the largest signal.
Power spectra for such time series are often characterized by a broad hump with a low-frequency
maximum and a decline at higher frequencies. Superimposed are a number of sharp tidal peaks near
diurnal and semi-diurnal frequencies, and sometimes a 
broader peak associated with Coriolis or inertial effects. 
Dynamical analysis requires the separation of the tidal signal from sub or
super-tidal variations, or in some cases that we separate tidal effects from inertial effects at
a nearby frequency. The tidal information information is either discarded kept for
further analysis.

Standard high/low/bandpass filtering techniques \cite{jackson:86} can be used but
in general these are inefficient because fairly narrow filters with a great deal of rejection
are needed. Also, although these are useful in analyzing non-tidal variability they provide
no compression of the tidal information. Specialized techniques have been devised to take
advantage of the ``deterministic'' nature of tidal processes. In classical harmonic analysis,
the tidal forcing is modelled as a set of spectral lines, i.e., the sum of a finite set of 
sinusoids at specific frequencies. These frequencies are
related directly to astronomical parameters \cite{godin:72}. They can be described as arising
from interactions of 6 basic periods: the lunar day (24.8 hour), the lunar month (27 days), the 
tropical year (1 year), lunar perigee (8.85 years), the precession
of the lunar nodes (18.6 years), and perihelion (21000 years) and are specified by 
various combinations of sums and differences 
of integer multiples of the 6 frequencies. The set of 6 signed integers required to describe a particular
frequency are then called {\it Doodson numbers}. 
Many of the more important frequencies have names
such as ``$M_2$'', ``$K_1$'', etc. From astronomical considerations an ``equilibrium'' response
can be predicted; this is the phase and amplitude that would be observed if the response of the
earth was fast enough that the surface deformation was effectively in equilibrium with the forcing at
all times. The real ocean is very definitely not in equilibrium with the tidal
forcing.  However, as tidal
amplitudes are small compared with total ocean depth the dynamics are very nearly linear,
implying that the forced response contains only those frequencies present in the
forcing. A least squares fit can be
used to determine the relative phase and amplitude of each frequency in the response. This
phase/amplitude data thus provides a compression of the data in the complete time series,
which can then be compared with similar data at other locations to understand the characteristics
of tidal dynamics, or can be used to synthesize time series of tidal effects at other times
for predictive purposes. 

There are several drawbacks to classical harmonic analysis. The first is that, ignoring
the modulation of perihelion which is effectively constant over historical time, an $\approx 18.6$ 
year time
series is required to resolve all of the listed frequencies (that is, the number of wavelengths of
each constituent in the record is at least 1 different from all other constituents). In practise
record lengths are generally one year or shorter.
In order to handle this issue
an assumption is made that the relative phase/amplitudes of sinusoids with similar frequencies
(i.e. those whose first three Doodson numbers are identical)
are identical to predicted in the equilibrium response. In such a cluster large equilibrium
peaks are surrounded by small subsidiary peaks in frequency-space which provide ``nodal
modulations'' (or more correctly ``satellite modulations'') to 
the main peak. The appearance of the total signal will be a sinusoid whose phase and
amplitude varies slowly with time. These changes are slow enough to be considered effectively
constant for record lengths of up to one year.
At much shorter record lengths another problem arises. The frequency resolution further
degrades until even dissimilar constituents are unresolveable. The best solution is to
apply inference. This technique for finding the absolute phase/amplitude requires that the
relative differences in phase/amplitude between the two unresolved constituents is known
from other nearby data.
If this is not the case, it
is thought best to discard the smaller constituents and fit only to the largest in a given
frequency interval.

Another drawback of classical analysis is that it provides no easy way
to determine whether the resulting phase/amplitude of a given sinusoid is meaningful in a
deterministic way (i.e. it is truly a tidal line), or whether it results from fitting to a
component of the broad-spectrum
variability. In general it is likely to include both and some kind of
confidence interval would be useful.
In order to address this issue the ``response'' method was invented \cite{response}.
Although this provides better results than classical harmonic analysis it has
not found widespread use. 

Further problems with classical harmonic analysis arise in coastal regions where the tidal
response is in the form of a wave propagating onshore. In large estuaries
the seasonal change in salinity and flow may change the dynamic response 
but as these changes can vary from year to year this is not really deterministic. 
Instead spectral peaks are broadened so that they are no
longer pure lines. Within smaller estuaries the
tidal height variations are significant compared to water column depth and a variety
of nonlinear effects can occur. For example, flood periods shorten and intensify and ebbs lengthen. As 
long
as these effects are resonably deterministic they may be handled by adding extra ``shallow
water'' constituents. More problematic in these regions are the effects of internal variability.
Tidal interactions with varying topography can produce large internal waves and bores whose
characteristics are highly sensitive to ambient stratification. In such cases the assumption
of ``line'' frequencies becomes questionable and other techniques such as wavelet analysis
have been suggested \cite{jay:99}.

Here we describe T\_TIDE, a package of routines that can be used to
perform classical harmonic analysis with nodal corrections,
inference, and a variety of user-specified options. There are several novel
features. First, although the harmonic analysis algorithm with nodal corrections etc. itself is not
original (other than conversion to complex algebra), being heavily based on \cite{foreman:77} and
\citeN{foreman:78}, it 
is implemented in MATLAB, an analysis package
widely used by oceanographers. This allows for easy use within the framework of a complete
analysis involving plotting of raw data, scatter plots, and so forth. Second, the code is
written directly in matrix terms and thus relatively easy to understand and modify
if required. Finally, in order to differentiate between true deterministic
(line) frequencies and broad-spectrum variability, confidence intervals for the estimated
tidal parameters are computed using one of several user-selectable algorithms. The
package is made up of a number of files each of which contains one or more functions.
User-callable functions generally have a ``{\it t\_}'' prefix to prevent namespace
collisions.

The paper is composed of four parts. In the first the form of the equilibrium
potential is described. In the second part the mathematical basis of the technique
for making phase/amplitude estimates is described. In the third part the generation of confidence
intervals is outlined. Finally an illustrative example is discussed. 

\section{Tidal potential}

The effect of gravitational force vectors from the sun and moon, $\bf F$, can be written as the gradient of
a scalar potential $V$, ${\bf F}= -\nabla V$. The magnitude of this potential at the earth's
surface at any time obviously is dependent on the relative positions of the earth, moon,
and sun. In Doodson's development (described in /cite{godin:72}) the potential is described as a function of
lunar time $\tau$, and other astronomical variables (which are also functions of time):
\begin{eqnarray}
 s & = & \text{mean longitude of moon} \nonumber \\
 h & = & \text{mean longitude of sun} \nonumber \\
 p & = & \text{longitude of perigee} \nonumber \\
 N' & = & \text{negative of the longitude of the ascending node} \nonumber \\
 p' & = & \text{longitude of perihelion }\nonumber 
\end{eqnarray}
where all terms are in units of cycles. These variables can be evaluated for a given Julian date
in subfunction {\it astron} of file {\it t\_vuf.m}.
Their effects are combined with the integer Doodson number set $\{i', j', k', l', m', n' \}$
into the astronomical argument $V_a = i' \tau + j' s + k' h + l' p + m' N' + n' p'$. Sets with a common $i'$ are
called a species (thus the slow, diurnal, and semidiurnal species for $i'=0,1,2$ respectively), 
and sets with common $i'j'k'$ are called a subgroup. The constituent 
frequency $\sigma$ is defined as $\sigma = 2\pi d V_a / dt $.
The tidal potential is then be written in the form
\begin{equation}
V = \sum_{i' = 0}^3 \left[ G_{i'} (\theta) \sum_{j' k' l' m' n'}  A'_{i'j'k'l'm'n'} \cos(2 \pi V_a)
 + G'_{i'} (\theta) \sum_{j' k' l' m' n'}  B'_{i'j'k'l'm'n'} \sin(2 \pi V_a) \right] \label{equilib}
\end{equation}
For a given Doodson number set either $A'$ or $B'$ is nonzero, but not both. These constants are 
tabulated and stored in data structures that can be loaded using {\it t\_getconsts}. The geodetic 
functions $G_{i'}$ and $G'_{i'}$ vary with species $i'$ and latitude $\theta$, and also depend on
such constants as the radius of the earth and the masses and separations of the earth, moon, and
sun. The equilibrium amplitude for a constituent is defined as $G A'/g$ or $G' B'/g$ and can
be generated for a particular latitude using {\it t\_equilib}.


\section{Phase/Amplitude Estimates}

The algorithm used here is generally as described by \citeN{godin:72}, \citeN{foreman:77}, 
and \citeN{foreman:78}. Unlike
these authors we use complex algebra directly rather than deal with sine and cosine
fits separately. This has the advantage of unifying the treatment for scalar (e.g.
pressure) and vector (e.g. horizontal currents) time series which are represented
as $u+iv$.

Consider a time series of observations $y(t), t=t_1, t_2, ... t_M$ arranged in a vector, 
where the observation
times are regularly spaced at an interval $\Delta t$ (default 1 hour) and $M$ is an odd number
(an endpoint is discarded if required).
The time axis is defined such that the origin (or central time) is at $t_{(M+1)/2}$.
Some missing observations can be handled
by using a ``missing data'' marker in the input vector (by MATLAB convention this is NaN, 
the IEEE arithmethic representation for Not-a-Number).
This regular interval restriction does not arise from the least-squares
fit itself but rather from the automated constituent-selection algorithm and is also a
requirement when spectra are estimated in one of the confidence interval algorithms.
This time series may be be  composed of either real or complex numbers. The time series is
passed to the analysis program {\it t\_tide} along with a variety of (mostly optional) parameters.
 
The tidal response is be modelled as
\begin{equation}
x(t) = b_0 + b_1 t + \sum_{k=1...N} \!\! a_k e^{i \sigma_k t } + a_{-k} e^{-i \sigma_k t} \label{complex}
\end{equation}
where $N$ constituents (each with unique Doodson number sets) are being used. 
Each constituent has a frequency $\sigma_k$ which is known from the development of the potential,
and a complex amplitude $a_k$ which is not known, although if $y(t)$ is a real time series
$a_k = \overline{a_{-k}}$. A possible offset and (optional) linear drift is handled by the first two terms.
The traditional approach uses real sinusoids:
\begin{equation}
x(t) = b_0 + b_1 t + \sum_{k=1...N} \!\! A_k \cos (\sigma_k t) + B_k \sin (\sigma_k t) \label{real}
\end{equation}
and is related to (\ref{complex}) by $A_k = a_k + a_{-k}$ and $B_k = i (a_k - a_{-k})$.
The real representation is more convenient for the linear error analysis described later. 

Constituents can be chosen from a list of 45 astronomical and 101 shallow-water
consituents. Data structures containing information about these constituents is
loaded using {\it t\_getconsts}. There are several alternatives for
selecting constituents. For general use an automated selection
algorithm (following \citeN{foreman:77}) is in place, which works as follows. A basis 
of all astronomical and 24 of the most
important shallow-water constituents are gathered together. 
All constituents are listed in order of
predefined importance. Less important constituents whose frequencies are less than
a Rayleigh resolution limit $\alpha (N \Delta t)^{-1}$ (with default $\alpha = 1$)
apart from more important constituents in frequency are 
discarded. Additional shallow-water constituents can be specified if required.
If the relative phase/amplitude of two constituents that are otherwise unresolveable is
known from other sources than an inference procedure can be carried out.
Alternatively, constituent lists can be explicitly specified.

The least squares fit are the coefficients minimizing
\begin{equation}
E = \sum_m | x(t_m) - y(t_m) |^2 = || T a - y ||^2
\end{equation}
where $y = [ y(t_1), y(t_2), ... y(t_M)]'$, $a=[b_0,b_1,a_1,a_{-1},a_2,a_{-2},...a_{-N}]'$, and
T is an $M\times 2N+2$ matrix of sinusoids evaluated at observation times. The solution is found using the 
Matlab ``\verb+\+'' matrix operator.

Once the fit has been performed various corrections are applied. These are generated in {\it t\_vuf}.
First, the phase of the constituent response is usually reported as ``Greenwich phase'' $g_k$, that is, phase 
referenced to the phase of
the equilibrium response at $0^{\circ}$ longitude (the Greenwich meridian). Alternatively this
can be interpreted
as reporting the phase of the response at the time when the equilibrium forcing is at its largest 
positive value at $0^{\circ}$ longitude.
It is simplest to find the fitted phase at the central time of the record ($t=0$); the equilibrium phase $v_k$ 
is then just $V_a$ for the given constituent computed at the julian date corresponding to this
central time, with possible 
adjustments of $1/4, 1/2,$ or $3/4$ cycle depending on whether $A$ or $B$ is nonzero, and their signs.

Second, if a latitude is specified then nodal or satellite corrections are computed as follows. Consider
a main peak of index $k$ with satellites with indices $kj$. The effect of the different satellites will 
be to slowly modulate the phase/amplitude of the main peak over various periods, all more than 8 years. 
Our fitted
response $\hat{a}_k$ over some period can then be written as a modification of the ``true'' response $a_k$,  
$\hat{a}_k = f_k \exp(i u_k ) a_k $ where $f_k$ and $u_k$ are called the nodal correction amplitude and
phase respectively. Then
\begin{eqnarray}
f_k e^{i u_k } a_k e^{i \sigma_k t} &=& a_k e^{i \sigma_k t } + \sum_j a_{kj} e^{i \sigma_{kj} t} \nonumber \\
f_k e^{i u_k }  &=&  1 + \sum_j \frac{a_{kj}}{a_k} e^{i (\sigma_{kj} - \sigma_k) t }  \nonumber \\
f_k e^{i u_k }  &\approx&  1 + \sum_j \frac{a_{kj}}{a_k}   
\end{eqnarray}
The final approximation will hold as long as $(\sigma_{kj} - \sigma_k) t$  remains 
``small'' (i.e. $N\delta t << 8$ years).
In general the true phases and amplitudes of the satellites are not known. However since their
frequencies are very similar to that of the main peak it is standard to assume the ratio of true amplitudes
is the same as the ratio of amplitudes in the equilibrium response,
and the difference in true phases will be equal to the difference in equilibrium phases. The nodal
corrections are thus computed from the equilibrium response (\ref{equilib}). 
A latitudinal 
dependence arises from the geodetic functions. $G'_1$ is zero at the
equator and a crude limiting is used to prevent some diurnal corrections from getting overly large.
The validity of using the latitude-dependent equilibrium response to predict an aspect of the dynamic behavior
in one part of an ocean basin in such a simple way is not clear.


The products of the analysis are now a pair of complex values (possibly corrected for nodal
modulation)  for
each constituent $k$, $\{ a_k, a_{-k} \}$.  These are now converted into standard parameters:
\begin{eqnarray}
U_{k} & = &  |a_k| + | a_{-k}| \label{ell1} \\
u_{k} & = &  |a_k| - | a_{-k}| \\
\theta_k & = & \frac{ang( a_k ) + ang( a_{-k} )}{2} \\
g_k      & = & v_k  - \frac{ang( a_k ) - ang( a_{-k}) }{2}  \label{ell2}
\end{eqnarray}
For horizontal currents these parameters describe the features of an ellipse traced out by
the tip of the velocity vector:
the length of the semi-major and semi-minor axis of the ellipse ($U_{k}, u_{k}$ 
respectively), the inclination
of the northern semi-major axis counter-clockwise from due east $\theta_k$, and the Greenwich phase $g_k$.
For scalar time series the parameter $U_k$ is merely called the amplitude, and 
$u_k, \theta_k \equiv 0$.
Note that the restriction of the definition of inclination to use the northern axis 
means that analyses of constituents
whose ellipses are aligned in an east/west direction may have inclinations that fluctuate between
close to $0^{\circ}$ and nearly $180^{\circ}$ due to noise. These apparently large jumps are an
artifact of the restriction but do not represent similarly large changes in the appearance of the
ellipse. 

Once ellipse parameters are found these can be used for further analysis. They can also be used to
generate predictions at other times using {\it t\_predic}. Nodal corrections in {\it t\_predic}
are computed at the time series mid-point so that it is an exact inverse of {\it t\_tide}.

\section{Confidence Intervals}

One drawback of classical harmonic analysis is that the degree to which a given constituent
represents true tidal energy as opposed to the energy of a broad-band nontidal process is
not determined. This is useful information for two reasons: first, it allows one to make better
estimates of the tidal behavior, and second it can allow one to quantitatively compare
different analyses. There are two steps to producing confidence intervals. First, we must
form an estimate of the characteristics of non-tidal or residual ``noise'' affecting 
the $a_k$ (or $A_k, B_k$). 
Second, we must convert these estimates into confidence intervals for the standard parameters
through a nonlinear mapping. We discuss the case of real time series first.

\subsection{Residual noise (real)}

After the harmonic analysis for an N-point real time 
series $x(t)$ is performed we examine the structure of the residual
series. In the simplest case the noise is gaussian and uncorrelated in time.
If this is the case then the total power $P_T = \sigma_x^2 = P/\Delta t$, where
$P$ is a spectral density. The amplitude of the fit to sine and cosine terms ($A$ and $B$ respectively) 
will be contaminated by errors arising from
unresolved noise components within a frequency interval of $\Delta f = (N \Delta t)^{-1}$
around the line. Thus $\sigma^2_{A} = \sigma^2_{B} = P\Delta f =  \sigma_x^2 /N$. It is unlikely that
a geophysical series will be spectrally flat, and a more sophisticated approach 
used in {\it t\_tide} is to find
a local value of $P$ by making a spectral estimate from the residual time
series (i.e. after the removal of all fitted constituents) and averaging the power over
frequency bins in a window around the frequency of any constituent, neglecting bins
in which fitted constituents reside. Here we chose a sequence of windows of width
0.4 cpd centered on $1,2,3,...$ cpd.

\subsection{Conversion to standard parameters (real)}

A conversion from errors in the cos/sine amplitudes to errors
in standard parameters (amplitude and phase) can be done  using
a linearized analysis. Consider a constituent $k$. Let
Let $\xi = F(A_k,B_k) $ be a nonlinear function of these
parameters, either the amplitude or the Greenwich phase. 
Then if $\{A_k,B_k\}$ are independent random variables, we
can find a linearized estimate of the standard error of $\xi$ in terms of
the standard errors of the sinusoid amplitudes:
\begin{equation}
\sigma_\xi^2 = \left( \frac{\partial F}{\partial A_k} \right)^2 \sigma_A^2 +
\left( \frac{\partial F}{\partial B_k} \right)^2 \sigma_B^2 
\end{equation}
where the partial derivatives can be derived exactly (but tediously) 
from (\ref{ell1}-\ref{ell2}). 

Alternatively the nonlinear mapping can be handled directly using a ``parametric bootstrap'' \cite{bootstrap}.
In this case we assume that the noise contaminating $A,B$ is gaussian 
with a known variance. 
We then simulate a number of realizations or replicated of the analysis by taking our estimate of the
sinusoid amplitudes and adding noise to them. All of these realizations are then converted
nonlinearly to standard parameters using (\ref{ell1}-\ref{ell2}) and an estimate of the 
standard error computed from this replicate dataset directly.

Once a standard error is found, 95\% confidence intervals can be estimated
using standard techniques. Alternatively a signal-to-noise power ratio (SNR) can be computed.
Simulations performed in {\it t\_synth} (and described in the text file {\it t\_errors})
show that 
the linear procedure appears to be adequate for real time series (e.g. 
tidal height), as long as the SNR $<10$, and is probably 
not bad for SNR as low as 2 or 3. The nonlinear procedure gives similar results to the
linearized procedure at high SNR, and is more accurate at low SNR.

\subsection{Residual noise (complex)}

A complex residual time series $u+iv$ can be modelled as bivariate white noise and variances
$\sigma^2_u, \sigma^2_v$ and covariance $\sigma_{uv}$ computed. 
If we assume further
that the noise in both components is not correlated 
($\sigma_{uv} \approx 0$) then a coloured bivariate noise model
can be used and variances assigned to real and imaginary parts of constituent amplitude
separately on the basis of local spectral
densities as described above. If it is suspected that ($\sigma_{uv} \ne 0$) then it is recommended that
the time series be rotated into a coordinate system in which this is true (e.g. into
along/across channel axes).


\subsection{Conversion to standard parameters (complex)}

The linearized
analysis now involves functions of 4 variables, since both $A_k = A_r + iA_i$ and $B_k = B_r+iB_i$ have
real and imaginary parts, and analytic expressions for partial derivatives 
\begin{equation}
\sigma_\xi^2 = \left( \frac{\partial F}{\partial A_r} \right)^2 \sigma_{A_r}^2 +
\left( \frac{\partial F}{\partial A_i} \right)^2 \sigma_{A_i}^2 +
\left( \frac{\partial F}{\partial B_r} \right)^2 \sigma_{B_r}^2 +
\left( \frac{\partial F}{\partial B_i} \right)^2 \sigma_{B_i}^2 
\end{equation}
become large.
Some analytical simplification is possible by
assuming that the all 4 are independent.

The bootstrap approach is now applied to the complex coefficients $a_k$. One 
complication that arises is that unless the noise is circular ($\sigma^2_u = \sigma^2_v, \sigma_{uv}=$)
the errors in $a_k$ and $a_{-k}$ are correlated with each other. The bootstrap process now
requires generating correlated noise replicates.




\section{Example}

The analysis of an example dataset provided in \citeN{foreman:77} is shown in Figure \ref{f1}
and table \ref{t1}. This example is included in datafile {\it t\_example.mat}. 
The example dataset consists of 66 days of hourly 
elevations with a 3 day gap. The automated
constituent selection algorithm selected 35 constituents. One shallower water constituent
(M10) was added and two other constituents (P1 and K2) analyzed via inference.
The coloured bootstrap analysis was used to determine significance and confidence
intervals. The time series can be loaded and analyzed using the demonstration script
{\it t\_demo}.
Table \ref{t1} gives the output of the program. In the first column the name of
the constituent is given. Significant constituents (those with SNR in the last column
$>1$) are marked with a ``*''. Frequencies are listed in cph and phases in degrees.
11 constituents
were judged to be
significant (only 6 would be so at an SNR cutoff of 2). Figure \ref{f1}a
shows the results of the 
analysis and the residual time series after removal of these significant constituents.
In Figure \ref{f1}b are shown the
analyzed phases. Note that the significant constituents generally have reasonably
small phase errors. One would expect insignificant constituents have large
errors but some of them are still reasonably small. 

\section{Summary}

Separation of tidal and nontidal energy in oceanic time series is an important task in
any analysis. Here we discuss the theoretical foundation and some implementation details
of a MATLAB package for classical harmonic analysis. The package can also compute
confidence intervals for the tidal parameters using one of three different sets of
assumptions about the structure of residual noise. An example is provided to show typical
results. The code is available at \verb+http://www.ocgy.ubc.ca/~rich+.

%\section*{acknowledgements}



\bibliographystyle{/users/rich/doc/local/jpo}
\bibliography{/users/rich/doc/bib/oceanog,/users/rich/doc/bib/math,/users/rich/doc/bib/tide}


%\psfig{file=tidepic.eps,width=5in}
\begin{table}
\begin{center}
\small
\begin{boxedverbatim}
file name: PAPEROUT.txt
date: 17-Aug-2001
nobs = 1584,  ngood = 1510,  record length (days) = 66.00
start time: 06-Jul-1975 01:00:00
rayleigh criterion = 1.0
Greenwich phase computed with nodal corrections applied to amplitude 
and phase relative to center time

x0= 1.98, x trend= 0

var(x)= 0.82196   var(xp)= 0.21224   var(xres)= 0.60972
percent var predicted= 25.8 %

     tidal amplitude and phase with 95% CI estimates

tide   freq       amp     amp_err    pha    pha_err     snr
 MM    0.00151    0.2121    0.503   263.34   161.41     0.18
 MSF   0.00282    0.1561    0.526   133.80   188.82    0.088
 ALP1  0.03440    0.0152    0.044   334.95   150.82     0.12
 2Q1   0.03571    0.0246    0.044    82.69   106.21     0.31
 Q1    0.03722    0.0158    0.045    65.74   160.30     0.12
*O1    0.03873    0.0764    0.055    74.23    43.35      1.9
 NO1   0.04027    0.0290    0.035   238.14    74.68     0.69
*P1    0.04155    0.0465    0.045    71.88    70.96      1.1
*K1    0.04178    0.1405    0.059    64.81    23.49      5.7
 J1    0.04329    0.0253    0.050     7.32   129.76     0.25
 OO1   0.04483    0.0531    0.059   235.75    72.96     0.81
 UPS1  0.04634    0.0298    0.055    91.73   137.06     0.29
 EPS2  0.07618    0.0211    0.030   184.59   104.65     0.51
*MU2   0.07769    0.0419    0.034    83.23    48.82      1.5
*N2    0.07900    0.0838    0.035    44.52    25.54      5.9
*M2    0.08051    0.4904    0.035    77.70     4.51  1.9e+02
 L2    0.08202    0.0213    0.037    35.22   113.22     0.33
*S2    0.08333    0.2197    0.038   126.72     9.14       34
*K2    0.08356    0.0598    0.043   149.12    46.60        2
 ETA2  0.08507    0.0071    0.033   246.05   207.25    0.048
*MO3   0.11924    0.0148    0.014   234.97    67.38      1.1
 M3    0.12077    0.0123    0.014   261.57    62.11     0.81
 MK3   0.12229    0.0049    0.012   331.60   144.92     0.18
 SK3   0.12511    0.0023    0.010   237.69   219.86    0.054
 MN4   0.15951    0.0092    0.011   256.47    69.76     0.68
*M4    0.16102    0.0126    0.011   291.78    65.09      1.4
 SN4   0.16233    0.0083    0.011   270.85    91.22     0.54
 MS4   0.16384    0.0010    0.008   339.35   248.82    0.015
 S4    0.16667    0.0047    0.010   299.56   142.32     0.23
 2MK5  0.20280    0.0013    0.005   310.10   181.03    0.067
 2SK5  0.20845    0.0045    0.006   104.00    99.71     0.64
 2MN6  0.24002    0.0035    0.007   271.24   133.30     0.22
 M6    0.24153    0.0017    0.006   158.88   197.43    0.093
 2MS6  0.24436    0.0056    0.008   306.10    90.03     0.54
 2SM6  0.24718    0.0023    0.007   298.92   175.13     0.11
*3MK7  0.28331    0.0086    0.006   212.25    44.21      2.1
 M8    0.32205    0.0030    0.004    42.43    75.29     0.55
 M10   0.40256    0.0009    0.003   198.23   209.99    0.089
\end{boxedverbatim}
\normalsize
\caption{Analysis for Tuktoyuktuk datset}
\label{t1}
\end{center}
\end{table}

\begin{figure}
\begin{center}
\psfig{file=tidepic.eps,width=5in}
\caption{Tuktoyuktuk analysis example. (a) Raw time series. 
(b) Upper curve is residual time series after removal of tidal signal.
Lower curve is a synthesized tidal series using significant constituents.
(c) Amplitude of all analyzed components with 95\% significance level. Note
frequency dependence. Significant constituents are marked with a solid circle.
(d) The phase of significant constituents with a 95\% confidence interval. 
}
\label{f1}
\end{center}
\end{figure}



\end{document}



